Cluster-based density-functional approach to quantum transport through molecular and atomic contacts

We present a cluster-based density-functional approach to model charge transport through molecular and atomic contacts. The electronic structure of the contacts is determined in the framework of density functional theory, and the parameters needed to describe transport are extracted from finite clusters. A similar procedure, restricted to nearest-neighbor interactions in the electrodes, has been presented by Damle et al. [Chem. Phys. 281, 171 (2002)]. Here, we show how to systematically improve the description of the electrodes by extracting bulk parameters from sufficiently large metal clusters. In this way we avoid problems arising from the use of nonorthogonal basis functions. For demonstration we apply our method to electron transport through Au contacts with various atomic-chain configurations and to a single-atom contact of Al.


Introduction
Advances in the experimental techniques for manipulating atomic-sized objects have turned the vision of molecular-scale electronic circuits into a realistic goal [1]- [5]. This has intensified the interdisciplinary efforts to study charge transport in nanostructures. Ideally, the circuits would be constructed in a bottom-up fashion with functional units and all the wiring on the molecular scale. To approach this goal, present-day experiments in the area of molecular electronics concentrate on measuring the current-voltage response of single molecules contacted to metallic electrodes. In these studies, also purely metallic atomic contacts serve as important reference systems [6].
In order to support the experiments and to stimulate further technological advance, theoretical modeling of charge transport at the atomic scale is needed. Here, one faces the challenge to describe infinitely extended, low-symmetry quantum systems that may, in addition, be far from equilibrium and involve strong electronic correlations. While a complete theoretical understanding is still lacking, sophisticated ab initio methods have been developed for approximate but parameter-free numerical simulations. In order to study the prototypical metalmolecule-metal systems or metallic atomic contacts, many groups employ density functional theory (DFT) combined with nonequilibrium Green's function (NEGF) techniques [7]- [22]. Some shortcomings related to the use of DFT in this context have been pointed out, and solutions are being sought [23]- [26]. On the other hand, DFT presently appears to be one of the few practicable ab initio electronic structure methods, since studies of quantum transport require dealing with a large number of atoms. Furthermore, it can handle the hybrid metalmolecule-metal contacts, where the central regions, comprising the molecule, are frequently rather insulator-like, whereas the electrodes are metallic. For a more complete discussion, we refer to [27]- [29]. 3 The DFT approaches can mainly be divided into two types. In the first type, atomic-sized contacts are modeled by periodically repeated supercells, and computer codes developed for solid-state calculations are employed [7,15,22]. The use of periodic boundary conditions facilitates the electrode description. However, typically the conductance is determined for an array of parallel junctions and may be affected by artificial interactions between them. The second type is based on finite clusters and originates more from the chemistry community [8,9,14,16,21]. It has the advantage that genuinely single-atom or single-molecule contacts are studied. The drawback is commonly the description of the electrodes, since it is difficult to treat bulk properties based on finite clusters. Furthermore, the coupling between the central region and the electrodes can be complicated by finite-size and surface mismatch effects.
To arrive at an ab initio DFT description it is desirable to treat the whole system consistently by using the same basis set and exchange-correlation functional everywhere. The problem of the cluster-based approaches regarding the electrodes is apparent, for example, from the work of [8,9,21], where the authors resort to a separate tight-binding parametrization obtained from the literature [30]. In such cases, uncertainties mainly arise from the coupling between the central region and the electrodes, which requires knowledge about level alignments and basis functions. Damle et al proposed to resolve this issue by extracting bulk parameters from finite clusters computed within DFT [16,31]. However, their treatment of the electrodes should be seen as a first approximation, since only couplings between nearest-neighbor atoms were considered. Furthermore, they finally use energy-independent self-energies, which is well justified only for electrode materials with a constant density of states (DOS) near the Fermi energy E F .
In this work, we present a cluster-based DFT approach for the atomistic description of quantum transport. We follow the ideas of Damle et al [31], but place special emphasis on the treatment of the electrodes. In particular, we show that extracting bulk parameters from small metal clusters can lead to an unphysical behavior of the overlap of the nonorthogonal basis functions in k-space. The description of the electrodes can be improved systematically by employing metal clusters of increasing size. Our implementation is based on the quantumchemistry package TURBOMOLE, which allows us to study clusters of several hundred atoms. In this way, we obtain an ab initio formulation of quantum transport in atomic-sized contacts, where the whole system is treated on an equal footing. It has the advantage that we can employ well-tested, high-quality quantum-chemical Gaussian basis sets. Since energies in our isolated systems are measured with respect to the vacuum level, we can directly compare values for E F to experimental work functions 8 . Furthermore, our method offers a high degree of geometric flexibility. Thus, we can describe contact configurations with molecules of large transverse extent and with electrodes that are laterally offset or tilted relative to each other.
The theoretical framework of our approach is presented in section 2. Several technical details related to the use of nonorthogonal basis functions and the electrode treatment can be found in appendices A-C. To demonstrate the usefulness of our method, we study in section 3 the transport properties of atomic contacts of Au and Al. The choice of these materials is motivated by the fact that Au exhibits a rather energy independent DOS near E F , whereas Al does not. Furthermore, for these systems, we can compare our results with the literature. We find good agreement, and demonstrate the sufficient robustness of our calculations. Further applications have been presented in [32]- [36]. We summarize our results in section 4.

Electronic structure
Our ab initio calculations are based on the implementation of DFT in TURBOMOLE 5.9 [37]. By ab initio we mean that the simulations require no system-specific parameters. Self-consistent DFT calculations of large systems are generally very time-consuming. TURBOMOLE, however, is specialized in handling such systems, and offers several possibilities to reduce the computational effort. Most essential are the techniques of the 'resolution of the identity in J ' (RI-J ) [38,39] and the 'multipole-accelerated resolution of the identity J ' (MARI-J ) [40], which are both implemented in the ridft module of TURBOMOLE. The approximations help to reduce the effort to compute the Coulomb integrals J , which are particularly expensive to evaluate. With the help of the RI-J approximation, known also under the name 'density fitting', the four-center-two-electron integrals can be expressed as three-center-two-electron ones [41]. Calculations are faster by a factor of 10-100 as compared with standard DFT, but equally accurate. The MARI-J technique adds to RI-J an efficient treatment of Coulomb interactions between distant atoms. The interactions are divided into a near-field and a far-field part, where the near field is treated with RI-J and the far field by a multipole approximation. Compared with RI-J , it can accelerate the calculations by another factor of 2-7 [40]. Besides, if the contact configuration admits, one can exploit point group symmetries, including non-Abelian ones. In this way, the calculations speed up further by a factor given by the order of the point group.
DFT requires the choice of an exchange-correlation functional [42]. We select the generalized-gradient functional BP86 [43,44], which is known to yield good results for large metal clusters [45]- [48]. Since it contains no contribution of Hartree-Fock exchange, it is efficiently evaluated within RI-J or MARI-J . To express the orbital wave functions, Gaussian basis sets of split-valence-polarization (SVP) quality are used [38,39,49], which are the TURBOMOLE standard. Within the closed-shell formalism of DFT, total energies of all our clusters are converged to a precision better than 10 −6 au. In order to obtain ground-state structures, the total energy needs to be minimized with respect to the nuclear coordinates. We perform such geometry optimizations or 'relaxations' until the maximum norm of the Cartesian gradient has fallen below 10 −4 au.

Transport formalism
We compute transport properties of atomic-sized contacts using the Landauer-Büttiker theory and Green's functions expressed in a nonorthogonal basis of atomic-like orbitals [50,51]. The local nature of the basis allows us to partition the basis states |i, α into left (L), central (C) and right (R) ones, according to a division of the contact geometry. In the basis states, α refers to the type of orbital at the position of atom i. For reasons of brevity, we will frequently suppress the orbital index. The Hamiltonian (or Kohn-Sham) matrix H iα, jβ = i, α|H | j, β , and analogously the overlap matrix S iα, jβ = i, α| j, β , can thus be written in block form, Both S and H are real-valued and are hence symmetric. In addition, we assume the C region to be large enough to have S LR = H LR = 0. Within the Landauer-Büttiker theory [52], the linear 5 conductance can be expressed as where f (E, T ) = {exp[(E − µ)/k B T ] + 1} −1 is the Fermi function, and the chemical potential µ is approximately equal to E F . Using the standard NEGF technique, the transmission function is given by with the transmission matrix Here, we define Green's functions and and the hopping-rate matrices where g r X X = (E S X X − H X X ) −1 is the electrode Green's function for lead X = L, R. At low temperatures, the expression for the conductance simplifies to with G 0 = 2e 2 / h the quantum of conductance and τ n the eigenvalues of t † t. The latter are the transmission probabilities of the transmission eigenchannels n. 9 Also other observables, such as the thermopower or the photoconductance, can be studied based on the knowledge of τ (E) [32,36,53,54]. Information on the position of energy levels of a system may help to identify conduction mechanisms. Such information can be extracted from the spectral density [48] Using this, we define the local density of states (LDOS) at atom i and its decomposition into orbitals α via In appendix A, we discuss the approximations involved in this definition. There, we also consider further the issues related to the use of nonorthogonal basis sets to evaluate the singleparticle Green's functions and the electric current. 9 The fact that the τ n are probabilities, i.e. that 0 τ n 1, can be proven by defining r = 1 + i √ L G a CC √ L such that r † r + t † t = 1. Using relations presented in appendix A.3, it is easy to show the positive-semidefiniteness of L , necessary for the determination of both r and t (equation (4)). Figure 1. Quantum transport scheme. The atomic-sized contact (a) is divided into a C region and two semi-infinite L and R electrodes. Using a similar division for the ECC (b), information on the electronic structure of the C region (S CC , H CC ) as well as the CL and CR couplings (S CL , H CL and S CR , H CR ) is extracted. In order to obtain the self-energies r L and r R , the remaining task is to determine the electrode surface Green's functions g r LL and g r RR . This procedure is described further below in the text.

Implementation of the transport method
2.3.1. Central system. In order to determine the transmission function τ (E), we need a practical scheme to obtain the necessary information on the electronic structure. In figure 1, we present our procedure. The goal is to describe the whole atomic-sized contact (figure 1(a)) consistently, by treating the L, C and R regions with the same basis set and exchange-correlation functional. We obtain the parameters S CC and H CC as well as the couplings to the electrodes S CX and H CX with X = L, R from the extended central system (ECC) ( figure 1(b)), in which we include large parts of the tips of the metallic electrodes. The division of the ECC into the L, C and R regions is performed so that the C region is identical to that in figure 1(a). The atoms in the L and R parts of the ECC (blue-shaded regions in figure 1(a)) correspond to that part of the electrodes, which is assumed to couple to C. The partitioning of the ECC is commonly made somewhere in the middle of the metal tips, and we will also refer to it as 'cut'. The electrodes (L and R regions in figure 1(a)) are modeled as surfaces of semi-infinite crystals, described by the surface Green's functions g r X X . They are constructed from bulk parameters, extracted from large metal clusters. Let us now discuss how this is accomplished.

Electrodes.
We extract bulk parameters describing perfect crystals from large metal clusters. The complete procedure, which aims at determining the surface Green's functions g r X X with X = L, R is summarized in figure 2. In this work, we study exclusively electrode materials with an fcc structure, of which Au and Al are examples.
In a first step (figure 2(a)), we construct spherical metal clusters, henceforth called 'spheres'. They are made up of atoms at positions { R j | R j = 3 n=1 j n a n ∧ | R j | R sphere } with the standard primitive vectors a n of the fcc lattice and the sphere radius R sphere . We will generally use the vector of integer indices j = ( j 1 , j 2 , j 3 ) to characterize the atomic position R j . We do not optimize the geometry of the spheres, but set the lattice constant a 0 to its experimental literature value [56]. Increasing the radius R sphere should make the electronic structure in the center resemble that of a crystal. From the clusters we extract the overlap and Hamiltonian between the central atom at position 0 and the neighboring ones at position j (including j = 0). Figure 2. In order to obtain the electrode Green's functions g r X X for lead X = L, R, we determine bulk parameters from large metal clusters. In a first step (a) we extract overlap and hopping elements, , from the cluster's central atom to all its neighbors. They are (b) symmetrized by imposing the space group of the electrode lattice. After (c) a rotation to adapt them to the orientation of the respective electrode, (d) g r X X is constructed with the help of a decimation procedure.
This yields the matrix elements S . A rotation may still be necessary to arrive at parameters S (X ) j0 , H (X ) j0 , which are adapted to the orientation of the electrode X = L, R (figure 2(c)). Details on the symmetrization procedure and the transformation of bulk parameters under rotations are presented in appendix B. The parameters S (X ) j0 , H (X ) j0 are finally employed to construct the semi-infinite crystals and to obtain the surface Green's function g r X X (figure 2(d)). Due to the finite range of the couplings S CX , H CX , we need to determine g r X X for the first few surface layers only (blue-shaded regions in figure 1(a)). We compute these with the help of the decimation technique of [57], which we have generalized to deal with the nonorthogonal basis sets [58]. The complete procedure is explained in detail in appendix C. The parameters S j0 , H j0 can be computed once for a given metal and can then be used in transport calculations with electrodes of various spatial orientations.
For Au (a 0 = 4.08 Å), we have analyzed spheres ranging between 13 and 429 atoms, whereas for Al (a 0 = 4.05 Å) they vary between 13 and 555 atoms. Since we want to describe bulk, the parameters extracted from the largest clusters will obviously provide the best description. There is, however, an additional criterion, which necessitates the use of large metal clusters for a reliable description of the electrodes. As discussed in appendix B.1, it is based on the positive-semidefiniteness of the bulk overlap matrix. We find a strong violation of this criterion, if the extraction of parameters is performed such that only the couplings of the central atom to its nearest neighbors are considered. As a further demonstration of the quality of our description, we show in appendix B.2 the convergence of the DOS with respect to R sphere .
For the transport calculations, we need a value for the Fermi energy. The biggest Au and Al spheres computed, Au 429 and Al 555 , respectively, are very metallic. They exhibit differences between the highest occupied molecular orbital and the lowest unoccupied molecular orbital of less than 0.06 eV. Therefore, we set E F halfway between these energies. In this way, we obtain E F = −5.0 eV for Au and E F = −4.3 eV for Al. The values will be used in all the results below. Notice that the negative values of E F agree well with experimental work functions of 5.31-5.47 eV for Au and 4.06-4.26 eV for Al [59].

Discussion.
In our approach, we assume the metal tips included in the ECC (figure 1(b)) to be large enough to satisfy basically two criteria. Firstly, all the charge transfer between the L and R electrodes and the C part of the contact should be accounted for. This ensures the proper alignment of the electronic levels in C with E F . Secondly, most of the metal tips, especially the L and R parts of the ECC, should resemble bulk as closely as possible. Under these circumstances, the highest occupied molecular orbital of the ECC should be very similar to E F , and we can evaluate the surface Green's functions by using bulk parameters of an infinitely extended crystal. Owing to the finite size of the cluster, the criteria can be satisfied only approximately. The mismatch between the parameters in the L and R regions of the contact and the ECC will thus lead to spurious scattering at the LC and CR interfaces. In principle this resistance can be eliminated systematically by including more atoms in the metal tips of the ECC. On the other hand, if the resistance in the C region is much larger than the spurious LC and CR interface resistances, they will have little influence on the results.
In order to describe the organic molecule in a single-molecule contact, it might be desirable to use a hybrid functional with a certain amount of Hartree-Fock exchange. However, the exchange interactions greatly increase the computational effort due to their nonlocal character, and may even lead to problems in the self-consistent iterations of DFT for rather metallic systems [60]. To treat the ECC and spherical metal cluster consistently, the exchange may be evaluated using a screened Coulomb potential [60]. A more approximate solution would be to mix two functionals, by including a portion of Hartree-Fock exchange for the ECC, but not for the sphere. In this way, the exchange would be taken into account, where it is most essential.

Metallic atomic contacts
In this section, we explore the conduction properties of metallic atomic contacts of Au and Al. These systems, in particular atomic-sized Au contacts, have been studied in detail both experimentally and theoretically, and can therefore be used to test our method. We start by discussing the transport properties of the Au contacts, consisting of a four-atom chain, a threeatom chain, and a two-atom chain or 'dimer'. Since Al does not form chains of more than two atoms [61], we consider only a single-atom contact. Results for Al dimer contacts were already reported in [33]. For all systems, we analyze the transmission, its decomposition into eigenchannels, and, in order to obtain knowledge about the conduction mechanism, the LDOS for atoms in the narrowest part of the contact. Moreover, we investigate the robustness of our transmission curves with respect to different partitionings of the large ECCs.

Gold contacts
Let us first discuss the electronic structure of Au, where we display the DOS in figure 3. The Fermi energy at E F = −5.0 eV is located in a fairly structureless, flat region somewhat above the d band. Based on the electronic configuration [Xe] 4f 14 5d 10 6s 1 of the atom, one might have expected a strong contribution only from the s orbitals at E F . But, as is visible from figure 3(b), s, p and d contributions are all comparable 10 . This signifies that valence orbitals hybridize strongly in the metal. When an atomic contact of Au is in a dimer or atomic chain configuration, a conductance of around 1G 0 is expected from experimental measurements [62]- [66] as well as from theoretical studies [7], [67]- [70]. The analysis shows that this value of the conductance is due to a single, almost fully transparent transmission eigenchannel. It arises dominantly from the s orbitals of the noble metal Au, since the electronic structure in the narrowest part of the contact resembles more the electronic configuration of the atom [66,67].

Determination of contact geometries.
Despite the consensus that the conductance of atomic chains of Au is around 1G 0 , the precise atomic positions play an important role [70,71]. Therefore it is necessary to construct reference geometries that have been studied with a well-established transport method. We choose to compare with results obtained with TRANSIESTA [7]. The ECCs investigated are shown in figure 4. The four-atom Au chain with electrodes oriented in the [100]-direction, called Au100c4, corresponds to a contact geometry examined in [72] (see figure 1(b) therein). The three-atom Au chain, Au111c3, is similar to a configuration in figure 9(d) of [7]. In addition, we study a Au dimer contact, Au111c2, where a two-atom chain is forming the narrowest part. In contrast to Au100c4, for the latter two contacts, the electrodes are along the [111]-direction. In each ECC, the main crystallographic direction is aligned with the z-axis, which is the transport direction.
Let us briefly explain, how we determine these geometries ( figure 4). For Au100c4, we construct two ideal, atomically sharp Au [100] pyramids, with two atoms in between. The pyramids end with the layer consisting of 25 atoms. The distance between the layers containing four atoms is set to 12.68 Å (figure 4(a)), as in [72]. Next, we relax the four-chain atoms without imposing symmetries, keeping all other atoms fixed. After geometry optimization, we find that the configuration agrees well with symmetry D 4h . We add two more Au layers with 16 and 9 atoms on each side, where the ECC now consists of 162 atoms, and perform a final DFT calculation, exploiting the symmetry D 4h . Compared with [72], all bond distances indicated in figure 4(a) agree within 0.01 Å, except for the distance between the central chain atoms, where our distance is shorter by 0.07 Å. For Au111c3, we proceed similarly to Au100c4 ( figure 4(b)  atoms to 9.91 Å [7], and cut the pyramid off at the layers containing 10 atoms. Then we add one atom in the middle, relax the three-chain atoms, add two layers on each side with 12 and 6 atoms, and perform a calculation in symmetry D 3d . Au111c3 consists of 77 atoms in total. Our bond distances agree with those in figure 9(d) of [7] to within 0.02 Å. For Au111c2, we include also the first Au layer in the geometry optimization process. The distance between the fixed layers with 6 atoms is 12.12 Å. Otherwise the steps are the same as for Au111c3. The ECC is computed in symmetry D 3d and consists of 76 atoms. In the parts excluded from the geometry optimization, atoms are all positioned on the fcc lattice, where we set the lattice constant to the experimental value of 4.08 Å, which corresponds to a nearest-neighbor distance of 2.88 Å.

Four-atom gold chain.
Let us now study the conduction properties for the contact Au100c4 ( figure 5(a)). There are different possibilities to partition the ECC into the L, C and R regions. The cuts should be done so that L and R are unconnected (S LR = 0 and H LR = 0, see section 2.2). Hence, the C region must be long enough. In order to describe well the coupling to the electrode surface (figure 1), it is furthermore necessary to have sufficiently many layers in the L and R regions. We observe that at least two of them are needed to obtain reasonable transmission curves. For the two different cuts of the ECC τ (E) is plotted in figure 5(b). In both cases it is found to be almost identical, indicating the reliability of our method. The transmissions at the Fermi energy are τ (E F ) = 0.93 and 0.98 for cuts 1 and 2, respectively. These values correspond well to the result τ (E F ) = 0.99 of [72]. For cut 2, it is visible in figure 5(c) that the transmission at E F is dominated by a single eigenchannel, in good agreement with experimental observations [66] and previous theoretical studies [7,68]. In general, the electronic structure at the narrowest part should have the most decisive influence on the conductance of an atomic contact. Therefore, we plot in figure 5(d) the LDOS of the atom indicated by the arrow in figure 5(a), resolved in its orbital contributions. Compared with the DOS of figure 3, it is dominated by s at E F , where the contributions of all other orbitals than s and p z are suppressed. These two orbitals form the almost fully transparent transmission eigenchannel, which is rotationally symmetric with respect to the z-axis.

Three-atom gold chain.
Exactly the same analysis will now be carried out for the contact Au111c3. In figure 6, the geometry of the ECC, the transmission for different partitionings, the transmission eigenchannels, and the LDOS of the central chain atom are shown. As for Au100c4, we observe that the different cuts yield very similar transmission curves ( figure 6(b)). Furthermore all the basic features in τ (E) are the same as in the TRANSIESTA calculation (see figure 11(d) of [7]). Above the d band, which exhibits a very narrow and high final peak, there is a dip in τ (E) in both cases. The transmission recovers, however, and a flat region with a value of around 1 is visible. At the Fermi energy cuts 1 and 2 yield τ (E F ) = 0.96 and 0.99, respectively. This is in reasonable agreement with τ (E F ) = 0.94 in [7], considering the differences in the electrode geometry, basis set, and exchange-correlation functional. We observe from figure 6(c) for cut 1 that the transmission at E F is dominated by a single eigenchannel, and the LDOS indicates a dominant contribution of s orbitals (figure 6(d)).

Two-atom gold chain.
The transmission and LDOS resolved into eigenchannels and orbital components, respectively, are shown in figure 7 for the dimer contact Au111c2. As for Au100c4 and Au111c3, we observe a single dominant eigenchannel at E F , and τ (E F ) = 0.96. The finding of such a dominant channel for chains of two or more atoms is in good agreement with our analysis of less symmetric contacts, which were based on a combination of a tightbinding model and classical molecular dynamics simulations [70]. However, that τ (E) increases partly even above one in the vicinity of E F signals that the influence of other channels is increased as compared with Au100c4 and Au111c3. Indeed, the LDOS of the atom in the narrowest part of the constriction (figures 7(a) and (c)) shows in particular increased p x and p y contributions. Also, the d states exhibit a less pronounced peak structure than was visible in figures 5(d) and 6(d). This is due to the higher coordination number of the atom and the enhanced coupling to the electrodes.

Aluminium contacts
As is visible from the DOS in figure 8, the electronic structure of Al differs substantially from that of Au. While the latter is a noble metal with an s valence, the Al atom has the electronic configuration [Ne] 3s 2 3p 1 with an open p shell, and the metal is hence considered sp-valent. The strong contribution of s and p states is also observable in the DOS, where d states play only a minor role. As compared with Au, the DOS exhibits a noticeable energy dependence around E F . For Al, we study an ideal fcc [111] pyramid, consisting of 251 atoms (figure 9(a)), henceforth referred to as Al111c1. Ideal means that the atoms are positioned on an fcc lattice with the experimental lattice constant a 0 = 4.05 Å. figure 9, the transmission is displayed for three different partitionings of the ECC Al111c1. Also shown are the transmission eigenchannels and the LDOS of the atom in the narrowest part of the contact for a selected cut. For energies below −6 eV, there are practically no differences visible between the curves for the three different partitionings. Nevertheless, some deviations arise at E F , and we obtain τ (E F ) = 2.36 (cut 1), 1.88 (cut 2) and 2.23 (cut 3). We attribute these variations to spurious scattering at the LC and CR interfaces, similar to Au. Our values for τ (E F ) of around 2 agrees nicely with those reported for single-atom contacts in [12]. Compared with Au, the structure of the transmission eigenchannels has changed in an obvious way. There are three channels at E F , which is in line with experimental observations of Scheer et al [66,73]. Due to the D 3d symmetry of the ECC, degeneracies arise, where in particular τ 2 = τ 3 . As is visible from the LDOS, these additional eigenchannel contributions mainly stem from the p x and p y orbitals, whereas s and p z form the nondegenerate first channel [66,67,74].

Conclusions
In conclusion, we have developed a cluster-based method to study the charge transport properties of molecular and atomic contacts. We treat the electronic structure at the level of DFT, and describe transport in terms of the Landauer formalism expressed with standard Green's function techniques. Special emphasis is placed on the modeling of the electrodes and the construction of the associated bulk parameters from spherical metal clusters. We showed that these clusters need to be sufficiently large to produce a reliable description of bulk properties, where a criterion for the extent of the spherical clusters is set by the overlap of the nonorthogonal basis functions. In our studies we crucially rely on the accurate and efficient quantum-chemical treatment of systems consisting of several hundred atoms, made possible by use of the quantum chemistry package TURBOMOLE. Compared with supercell approaches, our method has the advantage that we genuinely describe single-atom or single-molecule contacts.
As an application of our method, we analyzed Au and Al atomic contacts. Studying a four-, a three-, and a two-atom chain with varied electrode lattice orientations for Au, we found a conductance close to 1G 0 , carried by a single transmission channel. Next, we investigated an ideal Al single-atom contact, and found three transmission channels to contribute significantly to the conductance of around 2G 0 . These results are in good agreement with previous experimental and theoretical investigations. Although we observe some finite-size effects, we have demonstrated both for Au and Al a sufficient robustness of our transmission curves with respect to partitionings of the contact systems. The results illustrate the applicability of our approach to various electrode materials.
Beside the metallic atomic contacts examined here, we have applied the method in the field of molecular electronics. Studies include the dc conduction properties of dithiolatedoligophenylene and diamino-alkane junctions [32], [34]- [36] as well as oxygen adsorbates in Al contacts [33]. In addition, the thermopower [36] and photoconductance [32] of molecular junctions have been investigated in this way. Our studies demonstrate the value of parameterfree modeling for understanding transport at the molecular and atomic scale.

Acknowledgments
We thank R Ahlrichs for providing us with TURBOMOLE and acknowledge stimulating discussions with him and members of his group, in particular N Crawford, F Furche, M Kattannek, P Nava, D Rappoport, C Schrodt, M Sierka and F Weigend. This work was supported by the Helmholtz Gemeinschaft (contract no VH-NG-029), by the EU network BIMORE (grant no MRTN-CT-2006-035859), and by the DFG within the CFN and SPP 1243. MH acknowledges funding by the Karlsruhe House of Young Scientists and FP that of a Young Investigator Group at KIT.

Appendix A. Nonorthogonal, local basis sets
For practical reasons, one often employs nonorthogonal basis sets in quantum-chemical calculations, consisting for example of a finite set of Gaussian functions. The electronic structure is then described in the spirit of the linear combination of atomic orbitals (LCAO) [41,75,76], and this is also how TURBOMOLE is implemented. While it is in principle always possible to transform to an orthogonal basis, it may be more convenient to work directly with the nonorthogonal states.
A concise mathematical description using nonorthogonal basis states can be formulated in terms of tensors. The formalism is presented in a fairly general form in [77], where also aspects of second quantization are addressed. Below, we discuss some of the subtleties related to the use of nonorthogonal basis functions that are important for our method [58]. Since the basis functions are real-valued in our case, the full complexity of the tensor formalism is not needed [78,79]. Furthermore, we use a simplified notation, where all tensor indices appear as subscripts of matrices.

A.1. Current formula
The most important quantity for transport calculations is the electric current. In the NEGF formalism, its determination requires a separation of the contact into subsystems similar to figure 1(a) [52,80,81]. However, due to the overlap of the basis functions in a nonorthogonal basis, the charges of the subsystems are not well defined. Different ways of determining them exist, e.g. the Mulliken or Löwdin population analyses [75]. Despite these additional complications, the Landauer formula (equation (2)) can be derived in a similar fashion as for an orthogonal basis. Recent discussions of the derivation can be found in [51,82].

A.2. Single-particle Green's functions
Consider the single-particle Hamiltonian H describing the entire system. The retarded Green's operator is defined as G r (E) = [(E + i0 + ) − H ] −1 . Now consider the local, nonorthogonal basis |i with the matrix elements of the overlap S i j = i| j and the Hamiltonian H i j = i|H | j . Compared with section 2.2 the index i, used throughout this appendix, is a collective index, denoting both the position at which the basis state is centered and its type. The components of the retarded Green's function, defined by G r = i, j |i G r i j j|, 11 satisfy the equation [50].
The Green's function G r CC is defined as G r i j restricted to the central region C. It can be calculated according to equation (5). Due to the nonorthogonal basis the perturbation that couples C to the lead X = L, R and enters the self-energy (equation (6)) is given by H CX − E S CX and thus includes also an overlap contribution. It is interesting to observe that as E → ∞ the self-energies and the Green's function behave as Equation (A.4) describes the 'renormalization' of the inverse overlap matrix of C due to the coupling to the leads.

A.3. LDOS
We have defined the LDOS at atom i and its decomposition into orbitals α in equations (10) and (11). Let us discuss these definitions further. Consider the energy eigenstates |µ of the entire system, satisfying H |µ = ε µ |µ . In this basis, the spectral density of equation (9) has the components ρ µν (E) = µ|ρ(E)|ν = δ(E − ε µ )δ µν . Clearly, they fulfill the normalization If, instead, we consider the components defined by Performing a Löwdin orthogonalization of the basis the normalized LDOS is given as a diagonal element of the integrand. The function LDOS iα (E) of equation (11) (with iα → i) is an approximation to this, where all indices are restricted to C. Since ρ CC (E) = −Im[G r CC (E)]/π is a positive-semidefinite matrix, it is easy to show that LDOS iα (E) is positive for all E. However, the normalization ∞ −∞ dELDOS iα (E) = 1 is only approximately fulfilled. This could be corrected by multiplying in equation (11) with (S −1 ) −1/2 CC (equation (A.4)) instead of S 1/2 CC . But since the contributions X =L,R S CX (S X X ) −1 S X C constitute only a surface correction, their neglect may be justified for atoms in the middle of C.

B.1. Size requirement for cluster construction
How large do the spherical metal clusters (figure 2) need to be for a convergence of the bulk parameters? Since the matrix elements of the Hamiltonian and the overlap decay similarly with increasing interatomic distance (for exchange-correlation functionals without a contribution of Hartree-Fock exchange), we can concentrate on the overlap. Then, a rather well-defined criterion can be found. The clusters should be so large that the extracted bulk overlap matrix is positive-semidefinite.
We define states | k, α = j e i k· R j | j, α in k-space. Since S iα, jβ = i, α| j, β is a positivesemidefinite matrix [75], the same is true for the overlap in k-space , where we used that S lα,mβ = S (l−m)α,0β . In the expression, N is the number of atoms in the crystal and In order to study the positive-semidefiniteness of S αβ ( k, k ) it is hence sufficient to investigate the behavior of S αβ ( k). To do so for a complex quantum-chemistry basis set, we define the positive-definiteness measure In this equation, S( k) is the smallest eigenvalue of the matrix S αβ ( k), where S αβ ( k) is constructed from the bulk parameters extracted from a cluster with radius R sphere [S αβ ( k) = R j ;| R j | R sphere e −i k· R j S jα,0β ]. In the discrete Fourier transformations (see appendix C.3), we assume periodic boundary conditions with a finite periodicity length along the standard primitive lattice vectors [56,58]. R sphere must be chosen large enough for ξ to be positive or, if ξ remains negative, it must at least be sufficiently small in absolute value.
Let us first illustrate the behavior of ξ at the example of an s-orbital model. Gaussian s-functions are described by with an exponent γ , characterizing the radial decay. Hence, the overlap between two atoms decays with their distance R j = | R j | like a Gaussian function. We consider an infinitely extended chain with atoms at equally spaced positions along the x-axis ( R j = j 1 a 0 e x ). The overlap from a selected atom to its neighbors drops off exponentially ( figure B.1(a)). The Fourier transformation will again result in a Gaussian with purely positive values S ss (k x ) ( figure B.1(b)). If, however, overlap matrix elements are taken into account only up to a certain maximum value | R j | R sphere , as in a finite cluster, a rough sin(k x )/k x -behavior results, where S ss (k x ) becomes negative at certain k-values. Upon an increase of R sphere , S ss (k x ) will evolve into a Gaussian function and ξ will thus approach zero from below. The negative tails of S ss (k x ) are unphysical, and our observation implies that the clusters used to extract bulk parameters (figure 2) need to be of a sufficiently large radius R sphere , in order to obtain a reliable description of a crystal. Obviously, with our overlap-based criterion the required magnitude of R sphere depends on the basis set chosen.
In figure B.2, we plot the behavior of ξ as a function of R sphere for Au and Al. Beside the results for the SVP basis set, we display ξ for Au also for the basis set LANL2DZ used in [16,31]. It is visible that ξ is positive for a single atom (R sphere = 0), but negative for small spheres. With increasing R sphere , ξ approaches 0 from below similarly to the s-orbital model. We find that the elimination of diffuse functions reduces the radius R sphere for ξ to become positive or negligibly small. For practical reasons, it may happen that R sphere cannot be chosen large enough. In such a case, negative eigenvalues of S αβ ( k) can lead to negative eigenvalues of the hopping-rate matrices X (equation (7)), since ρ X X = −Im[g r X X ]/π may no longer be positive-semidefinite (see also the discussion in appendix A.3).

B.2. DOS
The DOS can be used as another measure for the convergence to a solid-state description. With a k-space Hamiltonian in an orthogonal basis set H orth ( k), it is given as The positive-definiteness measure ξ for Au and Al as a function of R sphere . Beside the SVP basis set, the behavior of ξ is shown for LANL2DZ used in [16,31]. S( k) (equation (B.2)) is evaluated at 32 3 k-points. The radius R sphere has been scaled with the respective lattice constants (a 0 = 4.08 Å for Au and a 0 = 4.05 Å for Al).
Here, α runs over all basis functions on a bulk atom, and G orth,r 00 (E) = 1BZ d 3 kG orth,r ( k, E)/V 1BZ with the volume V 1BZ of the first Brillouin zone (1BZ) 12  and by imposing the fcc space group (see also figure 2), and to carry out the Fourier transformation last. For parameters extracted from large enough clusters, we observe the equivalence of the DOS construction with respect to the two different orthogonal Hamiltonians H orth ( k). If ξ remains (slightly) negative due to a too small R sphere , then the construction of the DOS from the Löwdin orthogonalization in real space (procedure (ii)) is of a higher quality than that resulting from the orthogonalization in k-space (procedure (i)).
In figure B.3 we show the DOS as constructed via procedure (ii) with parameters extracted from different Au and Al spheres with 141 to 555 atoms. We observe that the DOS seems well converged with respect to R sphere both for Au and Al for the largest spherical clusters Au 429 and Al 555 .

B.3. Transformation of bulk parameters under rotations
We assume that two coordinate systems are connected by the rotation , where r = r . The transformation properties of the bulk parameters with Y = S, H are determined by those of the basis functions r | j, α = φ α ( r − R j ). 13 The Gaussian basis functions used by TURBOMOLE are characterized by the angular momentum l and the multiplicity ν = 1, . . . , 2l + 1, and α is a collective index for both. The rotated basis functions of angular momentum l can be expressed as [83] with the representation D l µν ( ) of the rotation . Using Y ( r ) = Y ( −1 r ), it can be shown that the bulk parameters of the two coordinate systems are related by where D( ) is the representation of in the employed basis set. By knowledge of the D l µν ( ), D( ) can be constructed by the addition of representations [83]. If there are n l basis functions of angular momentum l in the basis set describing Y jα,0β , we have where ⊕ denotes a direct sum. Let us now give the explicit formulae for the D l ( ). In this work only s, p, and d basis functions are used, and hence we restrict ourselves to l = 0, 1 and 2. Since s functions just depend on the radius, ψ 0 ( r ) = ψ 0 (r ), we have For l = 1, there are three p functions, p 1 = p x = f 1 (r )x, p 2 = p y = f 1 (r )y and p 3 = p z = f 1 (r )z, with a certain radial dependence f 1 (r ) [76]. Exploiting −1 = T , we obtain p i = 3 j=1 p j ji . Thus the 3 × 3 representation of the rotation for the p functions is − y 2 )/2 with some radial dependence f 2 (r ). The transformed d functions are given as d i = 3 j=1 d j D 2 ( ) ji with the 5 × 5 representation

B.4. Imposing the fcc-space-group symmetry
In this section, we consider how to impose the fcc space group on the parameters H Owing to surface effects, the translational symmetry is not fulfilled for the parameters H sphere jα,0β , as illustrated in figure B.4. Hence, although the deviations decrease with growing radius of the spheres, the translational symmetry needs to be enforced in order to describe a crystal. To avoid numerical errors, we impose at the same time the point-group symmetry O h , although that symmetry is already present due to the shape of our clusters.  In the first step, Green's functions G r 1,0 , G r 3,0 , G r 5,0 , G r 7,0 , . . . are eliminated, in the second step Green's functions G r 2,0 , G r 6,0 , G r 10,0 , G r 14,0 , . . . and so on. The renormalized couplings τ (n) 1 and τ (n) 2 are effective couplings between superlayers 0 and 2, 2 and 4, and so on for n = 1, between 0 and 4, 4 and 8 and so on for n = 2, and generally between superlayers with an index difference of 2 n .

C.2. Decimation procedure
For the calculation of the lead self-energy for side X = L, R (equation (6)), Green's function g r X X is needed only for the first few atomic layers of the surface, since the coupling elements H CX − E S CX have a finite range. Therefore, the semi-infinite crystal is divided into superlayers, each consisting of an equal number N 3 of atomic layers (figure C.1). The surface superlayer is given the index 0, and the number N 3 is determined by requiring the hopping elements to be finite only between the nearest-neighbor superlayers and between C and superlayer 0. We assume from now on that the two-dimensional Fourier transformation to k-space has been performed for all quantities. For notational simplicity, we will suppress the dependence on κ λ and on the index (X ). Furthermore, the atomic layer indices p, p and the atomic orbital (or basis function) indices α, β are replaced by superlayer indices m, m , thus where W = E + S 0,0 − H 0,0 , τ 1 = E + S 0,1 − H 0,1 , τ 2 = E + S † 1,0 − H 1,0 and S 0,0 = S † 0,0 , H 0,0 = H † 0,0 , S 1,0 = S † 0,1 , H 1,0 = H † 0,1 . Note that due to the imaginary part of E + , τ 1 and τ 2 are not the Hermitian conjugates of each other. In order to solve equations (C.4) and (C.5) for G r 0,0 , we now use the iterative decimation scheme [57], illustrated in figure C.1. In the nth step of the iteration, this procedure eliminates all Green's functions G r 2 n−1 (2 j−1),0 with j ∈ N\{0}. The remaining components are then determined through new effective 'onsite' and 'hopping' blocks, given by where the matrices of the 0th step are W (0) b = W (0) s = W , τ (0) 1 = τ 1 , and τ (0) 2 = τ 2 . In the iterative procedure, the couplings τ (n) 1 and τ (n) 2 describe effective couplings between superlayers with an 26 by choosing both electrodes along the same crystallographic direction and by orienting them in the same way. In such a case, the vectors c (X ) i are identical for X = L, R. Consequently, W , τ 1 , and τ 2 need to be constructed only once, and the decimation can be done simultaneously for both sides. For the convergence criterion, we use i, j τ (n−1) In this way, we have a direct control on how much the inverse surface Green's function (equation (C.6)) is still modified. Here, ε = 10 −6 H turned out to be sufficient to converge the transmission curves.