First-principles many-body calculations of electronic conduction in thiol- and amine-linked molecules
Abstract
The electronic conductance of a benzene molecule connected to gold electrodes via thiol, thiolate, and amino anchoring groups is calculated using nonequilibrium Green functions in combination with the fully selfconsistent GW approximation. The calculated conductance of benzenedithiol and benzenediamine is five times lower than predicted by standard density functional theory (DFT) in very good agreement with experiments. In contrast, the widely studied benzenedithiolate structure is found to have a significantly higher conductance due to the unsaturated sulfur bonds. These findings suggest that more complex gold/thiolate structures where the thiolate anchors are chemically passivated by Au adatoms are responsible for the measured conductance. Analysis of the energy level alignment obtained with DFT, Hartree-Fock and GW reveals the importance of self-interaction corrections (exchange) on the molecule and dynamical screening at the metal-molecule interface. The main effect of the GW self-energy is to renormalize the level positions, however, its influence on the shape of molecular resonances also affects the conductance. Non-selfconsistent GW calculations, starting from either DFT or Hartree-Fock, yield conductance values within 50% of the selfconsistent GW results.
pacs:
73.63.-b,73.40.Gk,85.65.+hI Introduction
The problem of first-principles calculation of electronic conduction in molecular systems is of longstanding interest. Over the last decade, advances in experimental techniques have allowed for fundamental studies of electron transport through few or even individually contacted molecules. The transport mechanisms observed in molecular junctions range from ballisticsmit02 over off-resonant tunnelingvenkataraman ; song ; xiao ; tsutsui ; kiguchi to the strong correlation Kondo and Coulomb blockade regimekubatkin ; park to vibration assisted hoppingchoi . The former two belong to the phase-coherent transport regime characteristic of relatively short molecules (as opposed to molecular wireschoi ) with good ”chemical” contact to the electrodes and is the main focus of the present work.
The lack of control over the atomistic details of the metal-molecule interface introduces a strong statistical element in measurements on single-molecule junctions which masks the relation between atomic structure and measured conductance. For example, published conductance values for the gold/benzenedithiol model system vary by several orders of magnitude, although recent independent studies of this system seem agree on a value around ()song ; xiao ; tsutsui ; kiguchi . It is, however, not clear which structure is responsible for this ”typical” conductance and recent experimental and theoretical work point to complex gold/thiolate structures involving two molecules bonding the same Au adatomref18 ; ulstrup ; voznyy ; ref19 ; ref20 ; ref21 . Scanning tunneling microscope experiments in solution have shown that the use of amino rather than thiol anchoring groups leads to more well defined junction propertiesvenkataraman and a conductance around (or just below) has also been reported for the gold/benzenediamine junctionvenkataraman ; kiguchi . This was, however, not confirmed by independent break junction experiments in vacuumruitenbeek .
The uncertainties related to the junction atomic structure renders theoretical benchmark calculations more important and more challenging at the same time. Here progress has been hampered by the inability of conventional density functional theory (DFT) methods, which for long have been the workhorse and the state of the art for quantum transport calculationsbrandbyge ; xue , to reproduce the conductance measured for even the simplest molecular tunnel junctionsmagnus ; evers . As a consequence, the most succesful studies have focused on qualitative trends in the dependence of conductance on e.g. molecular conformationmischenko , molecular lengthhybertsen , side group functionalizationsmowbray , or have focused on properties independent of the numerical value of the conductance like molecular vibrationsokabayashi ; raman . The main shortcoming of the DFT approach has been attributed to its band gap problem, i.e. the fact that DFT tends to underestimate energy gapsgodby1 , and therefore overestimates the conductance. Attempts to overcome this problem within a single-particle framework have mainly focused on self-interaction correction schemesburke ; sanvito ; baranger .
The well known success of the many-body GW methodhedin for quasiparticle (QP) bandstructure calculations has recently inspired its application to (simplified) transport problemsthygesen_jcp ; gw_prb ; darancet ; thygesen_prl ; leeuwen_trans ; leeuwen_trans2 ; spataru1 ; spataru2 . The fact that the GW approximation succesfully describes systems with highly diverse screening properties ranging from metalsbarth over semiconductorsusuda_hamada ; ku ; schilfgaarde_prb ; Shishkin2006 ; rinke ; orr ; kresse to moleculesrostgaard ; stan_eulett is essential for a correct description of metal-molecule interfaces where the electronic character changes from metallic to insulating over a few angstrom. In particular, screening by the metal electrons can have large influence on the QP energies of the adsorbed moleculeneaton ; juanma ; thygesen_image ; kaasbjerg – an effect completely missed by both local and hybrid density functionalsjuanma . This has recently motivated the use of semi-empirical schemes for correcting the DFT eigenvalues by a scissors operator prior to transport calculationsquek ; mowbray ; quek2 . While such schemes can be justified for weakly coupled molecules, they become uncontrolled in the relevant regime of covalently bonded molecules where the screening effects mix with charge transfer and hybridizationthygesen_image .
In this work we combine nonequilibrium Green function methods for electron transport with the fully selfconsistent GW approximation for exchange and correlation to establish a theoretical benchmark for the electronic structure and conductance of gold/1,4-benzenedithiolate (BDT), -benzenedithiol (BDT+H) and -benzenediamine (BDA) molecular junctions. We find a conductance of for BDA and for BDT+H in very good agreement with experimental data. In comparison, the conductance obtained from DFT is about five times higher while non-selfconsistent GW calculations produce conductances within 50% of the selfconsistent GW result. We argue that the BDT+H structure can be viewed as a simple model of recently proposed RS-Au(I)-SR gold/thiolate structures involing two molecules attached to the same Au adatom. The conductance of a simple BDT molecule between Au(111) surfaces is predicted to be on the order of by both DFT, Hartree-Fock and GW. The origin of the high conductance is due to an unsaturated sulfur orbital with energy just below the Fermi energy. In the BDT+H and RS-Au(I)-SR structures, the sulfurs are fully passivated and the -orbital moves away from leading to an effective decoupling of the CH moiety from the gold electrodes. We find that the main effect of the GW self-energy is to shift the molecular levels and can be modelled by a simple scissors operator. However, the energy dependence of the GW self-energy may also affect the shape of the transmission resonances and this can change the conductance by almost a factor of two.
Most implementations of the GW method invoke one or several technical approximations like the plasmon pole approximation, neglect of off-diagonal matrix elements in the GW self-energy, analytic continuations from the imaginary to the real frequency axis, neglect of core states contributions to the self-energy, neglect of self-consistency. The range of validity of these approximations has been explored for solid state systems by a number of authorsbarth ; usuda_hamada ; ku ; schilfgaarde_prb ; Shishkin2006 ; rinke , however, much less is known about their applicability to molecular and metal-molecule systems. For this reason our implementation of the GW method avoids all of these approximations and as such represents an exact treatment of the GW self-energy within the space of the employed atomic orbital basis set.
Although we compare our results to experimental data and discuss them in relation to the possible atomic structure of the junctions, we stress that the main purpose of this study is the benchmarking of quantum transport calculations for specific, idealized junction with particular focus on the role of electronic correlation effects.
Ii Method
We consider a quantum conductor consisting of a molecule connected to left (L) and right (R) electrodes. We shall assume that outside a certain region containing the molecule and part of the electrodes (the ”extended molecule”), the electron-electron interactions can be described by a mean field potential. The current through the molecule is then given byref30 ; nonorthogonal
(1) |
where the energy dependence of all quantities has been suppressed for simplicity. In this equation is the Green function matrix of the extended molecule evaluated in a localized basis, is the coupling strength between the extended molecule and the left/right electrode, and are the Fermi Dirac distribution functions of the two leads. Our implementation applies to the general case of a finite bias voltage, but in this work we focus on the zero bias conductance which can be expressed in terms of the transmission functionref30 ; nonorthogonal
(2) |
as . This formula was originally derived for non-interacting electrons, but is in fact valid for interacting electrons in the low-bias limitgodby . We have verified that this is indeed fulfilled to high numerical accuracy in our GW calculations by comparing the conductance obtained from Eq. (1) for small finite voltages with evaluated in equilibrium.
The retarded Green function of the extended molecule is calculated from
(3) |
where is a positive infinitesimal and
(4) | |||||
(5) | |||||
(6) |
denote the overlap matrix, Kohn-Sham Hamiltonian and exchange-correlation potential, respectively. The matrices are evaluated in terms of a basis consisting of numerical atomic orbitalsGPAW-LCAO , and are obtained from a DFT supercell calculation performed with the real-space projector augmented wave method GPAWref29 . The electrode self-energies are obtained from the Kohn-Sham Hamiltonian of a bulk DFT calculation using standard techniquesref10 .
The boundary conditions in the plane normal to the transport direction enter only via the electrode self-energies which are constructed from the electrode surface Green functionref10 . The latter should represent an infinite surface but is here approximated by that of a periodic supercell with Au atoms in the surface plane. We have found that this is a very good approximation when the surface Green function of the cell is evaluated at a general -point (we use in coordinates of the surface Brillouin zone basis vectors). Using high symmetry points, in particular the Gamma point, can be problematick-point .
The term is the deviation of the Hartree potential from the groundstate DFT Hartree potential contained in ,
(7) |
In this expression and are the interacting and Kohn-Sham density matrices, respectively, and
(8) |
is the bare Coulomb interaction in the atomic orbital basis. The factor 2 is due to spin. Ref. Walter2008, describes how the all electron Coulomb elements are obtained within the PAW formalism. The last term in Eq. (3) is the many-body exchange-correlation self-energy which in this work can be either the bare exchange potential or the GW self-energy. We note that setting and in Eq. (3) we obtain the Kohn-Sham Green function, .
ii.1 GW self-energy
The (time-ordered) GW self-energy is given by
(9) |
where is the screened interaction and the indices run over the atomic basis functions of the extended molecule. The screened interaction is calculated in the frequency domain as the matrix product with the dielectric function, , evaluated in the random phase approximation. The irreducible density reponse function is computed in the time domain,
(10) |
where the factor 2 accounts for spin. Setting yields the Hartree-Fock approximation which thus corresponds to complete neglect of screening or equivalently complete neglect of correlations. Note that Eqs. (9-10) involve time-ordered quantities defined on the Keldysh time contour. For completeness we provide the expressions for the real time components in Appendix A and refer the reader to Ref. gw_prb, for more details.
As explained above, the matrices , , and are calculated for the extended molecule. On the other hand it is clear that , and thus , should have contributions from polarization diagrams outside this region, see Fig. 1. Physically these diagrams describe the potential acting on an electron propagating on the molecule due to the polarization that it induces in the lead. For this reason the self-energy will not be fully converged at the ends of the extended molecule region. To overcome this problem we only use the part of corresponding to the molecule and replace the remaining parts of the xc self-energy of the extended molecule by the DFT xc-potential. Symbolically,
(11) |
We stress that although we only include on the molecule, the interactions between electrons on the molecule and electrons in the electrodes (leading e.g. to image charge renormalization of the molecular levels) are included via diagrams of the form shown in Fig. 1. We also note that the form (11) implies that all metal atoms, both those inside and outside the extended molecule, are consistently described at the same (DFT) level. For the size of the extended molecule we have found it sufficient to include the gold atoms which are nearest neighbors to the sulfur or nitrogen atoms, i.e two gold atoms for the tip structures and six for the flat structure depicted in Fig. 2, see Appendix B. We expect that this rather local screening response is special for covalent metal-molecule bonds.
ii.2 Time/frequency dependence
The time/frequency dependence of , , , and is represented on a uniform grid ranging from to eV with a grid spacing of eV. We have verified that the results are coverged with respect to both the size and spacing of this grid. Fast Fourier transform is used to switch between energy and time domains to avoid convolutions. The calculations are parallelized both over basis functions and over the time/frequency grid points. One should always have to ensure proper representation of possible bound states. However, we have found that the conductance, and more generally the transmission function at any given energy, can be linearly extrapolated to the limit. This extrapolation has been performed for all the results presented in this work.
The memory requirements for the GW calculations (defined mainly by the size of the and matrices) are approximately a factor 50 larger than for a corresponding DFT calculation. The GW calculations for the benzene junctions considered in the present work were performed in parallel on 100-400 cores and took about 2 hours per selfconsisistency iteration. In comparison a DFT calculation for the same system took around 5 hours on 8 cores.
ii.3 Product basis
The calculation of all of the Coulomb matrix elements, , is prohibitively costly for larger basis sets. Fortunately the matrix is to a large degree dominated by negligible elements. To systematically define the most significant Coulomb elements, we use the product basis technique of Aryasetiawan and Gunnarsson aryasetiawan . In this approach, the pair orbital overlap matrix
(12) |
where is used to screen for the significant elements of .
The eigenvectors of the overlap matrix Eq. (12) represents a set of “optimized pair orbitals” and the eigenvalues their norm. Optimized pair orbitals with insignificant norm must also yield a reduced contribution to the Coulomb matrix, and are omitted in the calculation of . We have found that the basis for can be limited to optimized pair orbitals with a norm larger than without sacrificing accuracy. This gives a significant reduction in the number of Coulomb elements that needs to be evaluated, and it reduces the matrix size of and correspondingly, see Appendix A.
ii.4 Valence-core exchange
Since both core and all-electron valence states are available in the PAW method, we can evaluate the contribution to the valence exchange self-energy coming from the core electrons. As the density matrix is simply the identity matrix in the subspace of atomic core states, this valence-core exchange reads
(13) |
where represent valence basis functions and represent atomic core states. This contribution is added to describing the valence-valence interactions. We limit the inclusion of valence-core interactions to the exchange potential, neglecting it in the correlation. This is reasonable, because the polarization bubble, , involving core and valence states will be small due to the large energy difference and small spatial overlap of the valence and core states. In general we have found that the effect of on molecular energy levels can be up to 1 eVrostgaard . For the benzene-like molecules considered in this work the effect is generally less than 0.4 eV for the frontier orbitals.
ii.5 Self-consistency
Since and depend on , the Dyson equation (3) must be solved self-consistently in conjunction with the self-energies. In practice, this self-consistency problem is solved by iteration. We have found that a linear mixing of the Green functions,
(14) |
with generally leads to selfconsistency within 20-30 iterations.
Fully selfconsistent GW calculations are not standard, and in fact only few previous calculations of this type have been reportedrostgaard ; barth ; stan_eulett . Conventional GW bandstructure calculations typically apply a one-shot technique where the self-energy is evaluated with a non-interacting Green function, , usually obtained from DFThybertsen2 . In comparison, the selfconsistent approach has the immediate advantage of removing the -dependence, i.e. it leads to a unique solution.
For an approach, like the present, where the chemical potential is fixed by the external boundary conditions, some kind of selfconsistency (though not necessarily the full GW selfconsistency employed here) is essential to ensure charge neutrality. This is particularly important for cases where a molecular resonance lies close to the chemical potential. A shift in the energy of such a resonance could lead to a large change in its occupation. In a selfconsistent calculation this shift would be counterbalanced, mainly by a change in the Hartree potential. On the other hand, the one-shot approach does not account for this effect and unphysical level alignments could occur as a result.
Finally, the fully selfconsistent GW approximation is a conserving approximation in the sense of Baymbaym62 . This becomes particularly important in the context of transport where it ensures that the continuity equation is satisfiedbaym62 ; gw_prb . We mention that the recently introduced quasi-selfconsistent GW method (not to be confused with the fully selfconsistent GW approximation used here), in which is chosen such as to mimick the interacting Green function as closely as possible, have shown that selfconsistency in general improves the band gaps of semiconductors as compared to standard one-shot calculations.schilfgaarde
Iii Results
In this section we discuss the results of our selfconsistent GW calculations for the electronic structure and conductance of the prototypical gold/1,4-benzenedithiolate (BDT), -benzenedithiol (BDT+H), and -benzenediamine (BDA) molecular junctions. We argue that the thiol structure can be considered as a simple model for more complex gold/thiolate structures which have been proposed recentlyref21 but which are currently too large to be treated satisfactorily at the GW level. The transport results are rationalized by considering the alignment of molecular energy levels in the junction. Here we show that both DFT and Hartree-Fock provide quantitatively and qualitatively wrong results by predicting a gap opening rather than reduction when the molecule is attached to electrodes. Finally, we investigate to what extent the GW results can be reproduced by a simple scissors operator applied to the DFT Hamiltonian.
iii.1 Junction geometries
The junction geometries were optimized using the real space projector augmented wave method GPAWref29 . We used a grid spacing of 0.2 Å and the PBE functional for exchange and correlationref25 . The molecules were attached to Au(111) surfaces, modelled by a seven layer thick slab, either directly (in the case of BDT) or via tips (in the case of BDT+H and BDA) as shown in Fig. 2. The surface Brillouin zone was sampled on a -point grid, and the structures including molecule, Au tips, and outermost Au surface layers were relaxed until the residual force was below 0.05 eV/Å. The three structures are shown in the upper panel of Fig. 2 and some key bond lengths are given herenote .
It is generally accepted that the hydrogen atoms dissociate from the thiol end groups forming a gold-thiolate structurebdt_s . Nevertheless, our total energy calculations show that the benzenedithiol structure has a slightly lower energy than the benzenedithiolate when inserted between two gold tips as shown in Fig. 2(middle). In these calculations the hydrogens were either taken to be in the gas-phase or are adsorbed on the Au tips. In both cases the energy gain is less than 0.1 eV at the equilibrium distance but increases to 0.3 eV for a junction stretched by 1 Å. We stress that our calculations do not include effects of entropy which becomes relevant at finite temperature, and furthermore assumes that hydrogen in the gas phase is the proper reference for the solvated proton and an electron in the electrode, i.e. the reaction + is in equilibrium (in electrochemistry language we assume the standard hydrogen potential at pH=0). For these reasons our calculations are not sufficient to address the relative stability of thiols vs. thiolates under the relevant experimental conditions.
Importantly we note that new experimental evidence for the chemical structure of the gold-thiolate interface at the Au(111) surfaceref18 ; ulstrup ; voznyy or at Au nanocluster surfacesref19 ; ref20 has recently emerged, pointing to the existence of polymeric SR-Au(I)-SR units where the formally oxidised Au(I) adatoms are chemically bound to thiolates and form a part of a more complex structure (see Ref. ref21, and references therein). These complexes are currently too challenging to treat satisfactorily at the GW level. However, we have found that the electronic structure of such complexes is quite similar to that of the dithiol structure, see Sec. III.7. This is because the hydrogen atoms play a role similar to that of the Au adatom in passivating the sulfur atoms. Therefore the transport properties of the dithiol structure should also be similar to those of the SR-Au(I)-SR units.
iii.2 Energy levels of isolated molecules
A natural requirement for a method intended to describe the energy levels of molecules in contact with electrodes, is that it should be able to describe the energy levels of isolated molecules. As we show below, the DFT approach fails completely in this respect underestimating the HOMO-LUMO gaps of the three molecules by 5-6 eV while GW energies lie within 0.5 eV of the target values.
The gas-phase molecular structures have been relaxed in a 16 Å cubic cell using GPAW grid calculations as described in the previous section. For consistency, all energy levels have been calculated using the same double-zeta (DZ) atomic orbitals basis set. This is the same basis as used for the molecules in the transport calculations presented in the next section. For the DFT and Hartree-Fock calculations we have found that the energies obtained with the DZ basis agree with accurate grid calculations to within 0.2 eV. For the GW calculations the energies are within 0.1 eV of those obtained with a DZ+polarization basis set.
Table 1 summarizes the results for the HOMO and LUMO orbital energies obtained from the DFT-PBE eigenvalues, selfconsistent Hartree-Fock, and selfconsistent GW. The energy levels have been identified as the peaks in the spectral function extrapolated to . Due to lack of accurate experimental data we have also performed DFT-PBE total energy calculations for the neutral, cation, and anion species to obtain the addition/removal energies (last column). This approach has been shown to produce very accurate estimates of the experimental ionization and affinity levels of small moleculesrostgaard .
Relative to this reference, the DFT eigenvalues underestimate the HOMO-LUMO gap of all three molecules by eV, Hartree-Fock overestimates it by eV, while the gap obtained with GW lies within eV. These trends are consistent with a recent study of ionization potentials of a large number of moleculesrostgaard . The main reason for the large underestimation of the gap by DFT is the presence of self-interactions in the PBE functional.sanvito On the other hand Hartree-Fock is self-interaction free; here, by virtue of Koopmans’ theorem, the overstimation of the gap can be seen as a result of neglect of orbital relaxations. The effect of the latter is included in GW via the screened interaction and this reduces the gap relative the unscreened Hartree-Fock result.
We furthermore note that DFT places the LUMO of BDA and BDT+H below the vacuum level thus incorrectly predicting the anions to be stable. For BDT, the LUMO level is predicted to be negative by all methods indicating the radical nature of this species. We note that our GW energies for BDT are in good agreement with previous MP2 calculationsmp2_BDT .
Molecule | Orbital | DFT-PBE | HF | GW | |
---|---|---|---|---|---|
BDA | HOMO | -4.1 | -7.2 | -6.2 | -6.8 |
(CHN) | LUMO | -0.9 | 3.9 | 2.9 | 2.3 |
H-L Gap | 3.2 | 11.1 | 9.1 | 9.1 | |
BDT+H | HOMO | -5.1 | -8.0 | -6.9 | -7.5 |
(CHS) | LUMO | -1.3 | 3.3 | 2.2 | 1.3 |
H-L Gap | 3.8 | 11.3 | 9.1 | 8.8 | |
BDT | HOMO | -5.7 | -8.6 | -7.9 | -8.3 |
(CHS) | LUMO | -5.1 | -1.6 | -2.3 | -2.7 |
H-L Gap | 0.6 | 7.0 | 5.6 | 5.6 |
iii.3 Conductance calculations
The transmission functions of the relaxed junction geometries were calculated as described in Sec. II using three different approximations for , namely the PBE xc-potential, the bare exchange potential (leading to Hartree-Fock theory), and GW. The former choice corresponds to the standard DFT-approach. All calculations employ a double-zeta basis set for the molecules and double-zeta with polarization for the Au atoms. The results are shown in the lower panels of Fig. 2, and the corresponding conductances are summarized in Table 2.
The conductance of BDA and BDT+H calculated with the fully selfconsistent GW approximation agree well with the experimental values reported in Refs. venkataraman, and kiguchi, for benzenediamine and Refs. song, ; xiao, ; tsutsui, ; kiguchi, for benzenedithiol as indicated by the grey boxes in Fig. 2(left+middle). In contrast, DFT and Hartree-Fock respectively overestimates and underestimates the experimental conductances by factors 5-10. Our DFT result for the BDA junction is in good agreement with previous calculationsquek ; mowbray .
It is striking that the conductance of the ”classical” BDT junction shown in Fig. 2(right), is predicted by all three methods to be significantly higher than the experimental value (we obtain the same high conductances for BDT between tips, i.e. the structure in Fig. 2(middle) without hydrogen on sulfur). Our DFT conductance is in overall good agreement with the large number of previous calculations for the same or similar similar structurebdt_calc . The high conductance is clearly due to a strong peak in the transmission function close to the Fermi level. The peak moves to higher energies when going from DFT-PBE over GW to Hartree-Fock, opposite to the trend normally seen for occupied states. The peak comes from an unsaturated -orbital on the sulphur atoms and is discussed in more detail in Sec. III.7.
These results suggest that the structures probed in experiments on benzenedithiol junctions involve a chemically passivated form of the thiolate linker group. As we show in Sec. III.7, the high conductance of BDT is due to unsaturated -states on the sulphurs with energy close to the Fermi level (the energy of this orbital does not change considerable with DFT, Hartree-Fock and GW). In the thiol and SR-Au(I)-SR structures, these states form bonds to H and the Au adatom, respectively, and are thereby shifted away from the Fermi level. On the other hand, the electronic structure and transmission functions of BDT+H and SR-Au(I)-SR structure are rather similar indicating that the BDT+H can be viewed as a simple model of the more complex SR-Au(I)-SR structure.
It should of course be kept in mind that experiments are performed in solution and at room temperature and are subject to variations in the detailed atomic structure. Thus the measured conductance values should not be considered as highly accurate references for theoretical calculations on idealized junctions.
Method | BDA | BDT+H | BDT |
---|---|---|---|
DFT-PBE | |||
HF | |||
GW | |||
GW(PBE) | |||
GW(HF) |
iii.4 Energy level alignment
In Fig. 3 we show the calculated HOMO and LUMO energy levels of BDT+H and BDA in the gas-phase and in the junction. All energies have been aligned relative to the vacuum level. At this point we note that an accurate description of the vacuum level, i.e. the work function, can in general be difficult to obtain with an atomic orbital basislorente . However, by using more diffuse basis functions (an energy shift of 0.01 eV has been used for all Au basis functions throughout this workGPAW-LCAO ) we obtain a work function for Au(111) of 5.4 eV in good agreement with the experimental value of 5.31 eVau_work . At the position of the molecule, i.e. in the region between the two tips, the electrostatic potential from a calculation where the molecule has been removed, converges to a constant value of 4.9 eV above the metal Fermi level. This value has been used as reference for the vacuum level in Fig. 3. In the junction where the levels are broadened due to hybridization with the metal states, the level positions have been defined as the first moment of the projected density of states of the relevant molecular orbital, . Here the are obtained by diagonalizing the KS Hamiltonian within the molecular subspace.
The orbital energies obtained from a GW calculation include the dynamical response of the electron system to the added electron/hole via the correlation part of the self-energy. In general the correlations will shift the occupied levels up and the empty levels down relative to the bare Hartree-Fock energies. When a molecule is brought from the gas-phase into a junction the electronic character of its environment changes from insulating to metallic. The enhanced screening should thus cause the gap to shrink (neglecting shifts due to hybridization) as compared to its gas-phase value. However, it has recently been shown for molecules weakly bonded to a metal surface, that this effect is completely missing in effective single-particle theories based on a (semi)local description of correlations neaton ; juanma .
As a result, in our DFT and Hartree-Fock calculations, the change in the frontier orbital energies induced by the coupling to gold is completely governed by the effect of hybridization which tends to open the gap by eV for both molecules, see top panel of Fig. 3. In contrast, the GW gap is reduced by eV due to the enhanced screening in the contact. Since the hybridization shift is of course also present in he GW calculation, we conclude that the enhanced screening due to the metal contacts reduces the HOMO-LUMO gap by eV relative to the value in the gas-phase.
Note that we refrain from using the term ”image charge effect” to describe the renormalization of molecular orbitals. This term is appropriate for weakly coupled molecules where the screening of the added electron/hole takes place within the metal. For intermediate or strongly coupled molecules, there is no clear distinction between metal states and molecular states, and the screening is more appropriately described as dynamical charge transferthygesen_image .
From Fig. 3 it follows that the HOMO level of the molecules in the junction is predicted by DFT-PBE to lie eV higher than obtained with GW. This agrees well with a recent study of benzenediamine derivatives on gold(111) which showed that DFT places the HOMO level about 1 eV too high compared to UPS measurementsref3 . The fact that the DFT-PBE description of the energy levels is better for the adsorbed rather than isolated molecules may be seen as a result of the metallic screening build into the DFT xc-functional via its origin in the homogeneous electron gasrohlfing . It should also be noted that the error (compared to GW) of the DFT eigenvalues is significantly larger for the LUMO than for the HOMO. This is in good agreement with previous plane wave calculations for molecule/metal interfacesjuanma .
iii.5 GW calculations
To test the role of selfconsistency (in the GW and Hartree self-energies) we have performed non-selfconsistent GW calculations using both the DFT-PBE Green function and the (selfconsistent) Hartree-Fock Green function as the initial . The results for the conductance of all systems are summarized in the last two rows of Table 2. We notice is that the conductance can vary by a factor of three depending on . We also note that the Hartree-Fock starting point comes closer to the selfconsistent result. This is because the Hartree-Fock Green function is an overall better approximation to the GW Green function than is the DFT Green function (see e.g. the comparison of frontier orbital energies in Fig. 3). As an example Fig. 4 shows the calculated transmission functions for the BDA junction in a larger energy range around .
For the BDA and BDT+H junctions, GW[PBE] overestimates the conductance while GW[HF] underestimates the conductance relative to GW. This can be understood as follows: Since DFT-PBE (Hartree-Fock) underestimates (overestimates) the HOMO-LUMO gap, the use of these Green functions to evaluate the GW self-energy will lead to an overestimation (underestimation) of the screening. As discussed above the screening contained in the correlation part of the GW self-energy tend to reduce the HOMO-LUMO gap. This reduction of the gap is thus overestimated in the GW[PBE] calculation and underestimated in the GW[HF] calculation which explains the trends in conductance, i.e. more screening smaller gap higher conductance.
The GW results for the strongly coupled BDT junction show less variation. This is perhaps surprising given the large difference between the DFT-PBE and Hartree-Fock results shown in Fig. 2. However, DFT-PBE and Hartree-Fock give almost equal density if states at the Fermi level (the transmission functions at the Fermi energy are also rather similar), and therefore the screening contribution obtained with the two choices for G is also very similar.
iii.6 Scissors operator
In this section we investigate to what extent the GW transmission function can be reproduced by a DFT calculation where the occupied and unoccupied molecular levels have been shifted rigidly to match the GW energies. This is illustrated by applying a scissors operator (SO) to the Au-BDA-Au junction shown in Fig. 2(left). After shifting the molecular levels we perform a non-selfconsistent calculation of the transmission function. A similar procedure has previously been successfully used to include image charge effects and correct for self-interaction errors in DFT-transport calculationsquek ; mowbray ; quek2 . We refer the reader to Ref. mowbray, for more details on the SO technique applied here.
In Fig. 5 we show contour plots of the DFT conductance for the BDA junction where the shift of the occupied and unoccupied molecular orbital energies has been varied independently over 4 eV. The three values for the contour lines shown correspond to the conductance obtained with GW (0.0042G), the experimental conductance (0.0064G) and the DFT conductance (0.021G), respectively. The dot indicates the energy shifts which make the HOMO and LUMO levels of the DFT calculation match the corresponding levels of the GW calculation; the required SO shifts are -0.6eV and 3.8eV for the occupied and unoccupied molecular orbitals, respectively. Shifting the levels by this amount reduces the DFT conductance by a factor of 3 from 0.02G to 0.0072G. Interestingly, the GW conductance is not reproduced by these shifts; it is even lower by a factor of 1.5. In fact, to reproduce the GW conductance a shift of about -1.3eV of the DFT HOMO is required (keeping the LUMO position fixed at the GW position). This shows that while the renormalization of the molecular level energies can explain the main part of the difference between the DFT and GW conductance, the different shape of the transmission resonances also plays a role. This is clear from the inset which shows the GW transmission function (green) and the DFT transmission with SO chosen to match the GW HOMO and LUMO levels (the dot in the main figure). The lower conductance obtained with GW is seen to be a consequence of a faster decay of the HOMO resonance towards the Fermi level. This difference comes from the energy dependence of the GW self-energy.
iii.7 Thiol vs. thiolate structures
In Fig. 6 we compare the DFT transmission functions for: (a) the “classical” structure of benzenedithiolate [structure in Fig. 2(right)] (b) benzenedithiol between tips on Au(111) [structure in Fig. 2(middle)], and (c) benzenedithiolate in a SR-Au(I)-SR molecular unit form [structure 2 of Fig. 1 in Ref. ref21, ]. The conductance (essentially the transmission function at the Fermi energy) obtained for structures (b) and (c) is very similar but markedly different from (a). The high conductance of the structure (a) is due to the strong transmission resonance lying just below the Fermi level. This resonance has also been found in many previous studies and seems to be a characteristic and robust feature of this junction. As shown in this work this transmission resonance is also present in GW and Hartree-Fock calculations. In contrast, for both structure (b) and (c) the transmission function is rather flat in an energy window of eV around the Fermi level.
To examine the origin of the different transmission functions found for structures (a) and (b,c) we found it useful to consider the molecular orbitals of the CH moiety of the three molecules. By diagonalizing the (Kohn-Sham) Hamiltonian of this part of the molecules we found that the orbital depicted to the left in the upper panel of Fig. 7, and which constitute the HOMO-1 of the CH moiety, is responsible for all the structure in the transmission functions below the Fermi level. Note, that the orbitals obtained in this way are different from the HOMO and LUMO levels shown in Fig. 3 which were obtained by diagonalizing the Hamiltonian for the entire molecules including the SH and NH end groups.
The different panels of Fig. 7 show the projected density of states (PDOS) of the HOMO-1 for the three structures together with the PDOS of the sulfur -orbital to which the HOMO-1 couples. Within the Newns-Anderson modelnewns , the sulphur -orbital is the so-called group orbital. Note, that the PDOS of the -orbital has been calculated in the absence of coupling to the HOMO-1 as it should in the Newns-Anderson model. The transmission function is then essentially the product of the PDOS of the HOMO-1 and the PDOS of the group orbitalwannier_thygesen . It is clear that the origin of the (double) transmission peak around -1 eV for structure (a) is due to a resonance formed by the HOMO-1 and the sulfur -orbital. The chemical passivation of sulfur, by either hydrogen or the Au adatom, implies that the PDOS of the -orbital splits into bonding and anti-bonding states around -3 eV and 3 eV, respectively. This in turn shifts the PDOS of the HOMO-1 down in energy. In particular its magnitude around the Fermi level is lowered and as a consequence the transmission function (being essentially the product of the two curves) is suppressed in a broad energy window around . Thus chemical passivation of the sulfur is the main reason for the lower conductance observed in structures (b,c) as compared to (a). A secondary effect, giving rise to differences in the transmission of structure (b) and (c) is the different interface dipoles which shift the electrostatic potential at the CH moiety by different amounts. This shift, however, has little influence on the conductance due to the flatness of the transmission functions around .
To verify this scenario, we have applied a scissors operator of 1 eV to the CH moiety of structure (b) in order to align the onsite energy of the HOMO-1 to that of structure (c). The resulting transmission function (not shown) is very similar to that of (c) and the conductance is compared to obtained for structure (c). On the other hand, the conductance of BDT cannot be brought below by sifting the levels of the CH moiety down (by up to 2 eV) because of the strong peak around -1 eV due to the sulphur -orbital.
Although the above picture is based on the DFT electronic structure, the qualitative similarities of the DFT and GW transmission functions in Fig. 2 implies that the same picture applies to the GW calculations.
Iv Summary
We have presented a first-principles method for modelling quantum transport in molecular nanostructures beyond the single-particle approximation. The method is based on non-equilibrium Green functions and applies to the general case of a finite bias voltage, but in this work we focused on the zero bias regime. The conductance of benzenedithiolate (BDT), benzenedithiol (BDT+H), and benzenediamine (BDA) was calculated using the selfconsistent GW approximation. In contrast to standard DFT and Hartree-Fock methods, the GW approximation was found to yield consistently accurate values for the energy levels of both isolated and contacted molecules due to its proper treatment of self-interaction and dynamical screening. In general, results obtained with GW for the electronic conductance and energy gaps of contacted molecules lie in between the values obtained with DFT and Hartree-Fock. The latter methods respectively overestimates and underestimates the screening and none of them are able to detect the change in the molecule’s electronic environment when it is coupled to electrodes.
Non-selfconsistent GW calculations were found to yield conductance values within 50% of the GW results depending on the initial G. It was shown that the main difference between the GW and DFT calculations comes from the renormalization of the position of the molecular energy levels. However, changes in the shape of transmission resonances, and thus the conductance, also occur due to the energy-dependence of the GW self-energy.
The GW conductance of benzenediamine was found to be in good agreement with experiments. For benzenedithiolate between Au(111) surfaces we found a conductance well above and thus conclude that this structure cannot be responsible for the measured conductance around 0.01. On the other hand, a conductance close to this value was found for a hydrogenated benzenedithiol junction which, as we demonstrated, represents a reasonable model of more complex gold/thiolate structures where the chemical passivation of sulfur is provided by a gold adatom rather than hydrogen. These gold/thiolate structures are presently too demanding for our selfconsistent GW method, however, our results suggest that such structures are responsible for the 0.01 conductance reported experimentally.
In conclusion, we showed that a consistent and quantitatively accurate description of energy level alignment and charge transport in phase-coherent molecular conductors can be obtained from first-principles with the (computationally feasible) GW approximation. We believe this development is important for increased the synergy between theory and experiments in molecular electronics which is essential for continued progress in the field.
V Acknowledgements
We would like to thank Angel Rubio and Karsten Jacobsen for inspiring discussions. The authors acknowledge support from the Lundbeck Foundation’s Center for Atomic-scale Materials Design (CAMD), the Danish Center for Scientific Computing and the Academy of Finland. Part of the computations were done at the CSC - the Finnish IT Center for Science in Espoo.
Appendix A The GW self-energy
Let denote the rotation matrix that diagonalizes the pair orbital overlap , i.e. . The columns of are truncated to those which have corresponding eigenvalues . We then only calculate the reduced number of Coulomb elements
(15) |
where are the optimized pair orbitals
(16) |
which are mutually orthonormal, i.e. .
The determination of the GW self-energy proceeds by calculating first the lesser and greater components of the polarization matrix in the time domain
(17) | ||||
(18) |
where the factor 2 accounts for spin and we have used . This is then downfolded to the reduced representation
(19) |
The screened interaction can be determined from the lesser and greater polarization matrices, and the static interaction , via the relations
(20) | ||||
(21) | ||||
(22) | ||||
(23) |
where all quantities are matrices in the optimized pair orbital basis and matrix multiplication is implied. We obtain the screened interaction in the original orbital basis from
(24) |
which is an approximation due to the truncation of the columns of . Finally the GW self-energy can be determined from
(25) | ||||
(26) |
The exchange and Hartree potentials are given by
(27) | ||||
(28) |
The latter equals the first term in Eq. (7).
Appendix B Size of the ”extended molecule”
As explained in Sec. II.1 the GW self-energy of the molecule has contributions from screening diagrams extending into the electrodes. This is included by calculating and for the extended molecule which includes the Au atoms closest to the molecule. In Fig. 8 we show the dependence of the transmission function on the size of the extended region. The latter has been varied from zero Au atoms (extended molecule coincides with molecule), to the Au tip atoms, to the four-atom Au pyramids on each side of the molecule. From this we conclude that it is sufficient to include polarization diagrams for the Au tip atoms as was done for the calculations presented in the main text. We mention that, for simplicity, the calculations shown in Fig. 8 have been performed at the GW(PBE) level and using a minimal single-zeta basis set. However, the conclusions regarding the convergence of the GW self-energy with respect to the size of the extended molecule, should be equally valid for the case of selfconsistent GW with the DZ/DZP basis set used for all calculations in the main text.
References
- [1]
- [2] R. H. M Smit, Y. Noat, C. Untiedt, N. D. Lang, M. C. van Hemert, and J. M. van Ruitenbeek, Nature (London) 419, 906 (2002)
- [3] L. Venkataraman, J. E. Klare, I. W. Tam, C. Nuckolls, M. S. Hybertsen, and M. L. Steigerwald Nano Lett. 6, 458 (2006)
- [4] H. Song, Y. Kim, Y. H. Jang, H. Jeong, M. A. Reed, and T. Lee Nature 462, 1039 (2009)
- [5] X. Xiao, B. Xu, and N. J. Tao, Nano Lett. 4, 267 (2004)
- [6] M. Tsutsui, M. Taniguchi, T. Kawai Nano Lett. 9, 2433-2439 (2009).
- [7] M. Kiguchi, H. Nakamura, Y. Takahashi, T. Takahashi, and T. Ohto, J. Phys. Chem. C 114, 22254 (2010)
- [8] S. Kubatkin, A. Danilov, M. Hjort, J. Cornil, J.-L. Bredas, N. Stuhr-Hansen, P. Hedeård, and T. Bjørnholm, Nature 425, 698 (2003)
- [9] W. J. Liang, M. P. Shores, M. Bockrath, J. R. Long, and H. Park, Nature 417, 725 (2002)
- [10] S. H. Choi, B. Kim,C. D. Frisbie Science 320, 1482-1486 (2008).
- [11] A. Cossaro et al. Science 321, 943-946 (2008).
- [12] Y. Wang et. al, J. Phys. Chem. C 113, 19601 (2009)
- [13] O. Voznyy et al. J. Am. Chem. Soc. 131, 12989 (2009)
- [14] P. D Jadzinsky et al. Science 318, 430-433 (2007).
- [15] M. Walter et al. Proc. Natl. Acad. Sci. (USA) 105, 9157-9162 (2008).
- [16] M. Strange, O. Lopez-Acevedo, H. Häkkinen, J. Phys. Chem. Lett. 1, 1528-1532 (2010).
- [17] C. A. Martin, D. Ding, H. S. J. van der Zant, and J. M. van Ruitenbeek, New J. Phys. 10, 065008 (2008)
- [18] M. Brandbyge, J. L. Mozos, P. Ordejon, J. Taylor, and K. Stokbro, Phys. Rev. B 65, 165401 (2002)
- [19] Y. Q. Xue, S. Datta, and M. Ratner Chem. Phys. 281, 151 (2002)
- [20] M. Paulsson, C. Krag, T. Frederiksen, and M. Brandbyge Nano Lett. 9, 117 (2009)
- [21] C. Li, I. Pobelov, T. Wandlowski, A. Bragrets, A. Arnold, and F. Evers J. Am. Chem. Soc. 130, 318 (2008)
- [22] A. Mishchenko et. al Nano Lett. 10, 156 (2010)
- [23] M. S. Hybertsen et al. J. Phys:Cond. mat. 20, 374115 (2008)
- [24] D. J. Mowbray, G. Jones, and K. S. Thygesen J. Chem. Phys. 128, 111103 (2008)
- [25] N. Okabayashi, M. Paulsson, H. Ueba, Y. Konda, and T. Komeda Nano Lett. 10, 2950 (2010)
- [26] P. Y. Heayoung et. al Nano Lett. 10, 2897 (2010)
- [27] R. W. Godby, M. Schlüter, and L. J. Sham, Phys. Rev. Lett. 56, 2415 (1986)
- [28] C. Toher, A. Filippetti, S. Sanvito, and K. Burke Phys. Rev. Lett. 95, 146402 (2005)
- [29] C. Toher and S. Sanvito, Phys. Rev. Lett. 99, 056801 (2007)
- [30] S. H. Ke, H. U. Baranger, and W. Yang J. Chem. Phys. 126, 201102 (2007)
- [31] L. Hedin Phys. Rev. 139, A796 (1965)
- [32] M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
- [33] C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen Phys. Rev. B 81, 085103 (2010)
- [34] A. Stan, N. E. Dahlen, and R. van Leeuwen, Europhys. Lett. 76, 298 (2006)
- [35] K. S. Thygesen and A. Rubio, J. Chem. Phys. 126, 091101 (2007).
- [36] K. S. Thygesen, A. Rubio Phys. Rev. B 77, 115333 (2008).
- [37] P. Darancet, A. Ferretti, D. Mayou, and V. Olevano, Phys. Rev. B 75, 075102 (2007)
- [38] K. S. Thygesen, Phys. Rev. Lett. 100, 166804 (2008)
- [39] P. Myöhanen and A. Stan and G. Stefanucci and R. van Leeuwen, Euro. Phys. Lett. 84, 67001 (2008)
- [40] P. Myohanen, A. Stan, G. Stefanucci, R. van Leeuwen, Phys. Rev. B 80, 115107 (2009)
- [41] C. D. Spataru and M. S. Hybertsen and S. G. Louie and A. J. Millis, Phys. Rev. B 79, 155110 (2009)
- [42] C. D. Spataru, Phys. Rev. B 82, 195111 (2010)
- [43] B. Holm and U. von Barth Phys. Rev. B 57, 2108 (1998)
- [44] M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
- [45] M. Usuda, N. Hamada, T. Kotani, and M. van Schiflgaarde, Phys. Rev. B 66, 125101 (2002)
- [46] W. Ku and A. G. Eguiluz, Phys. Rev. Lett. 89, 126401 (2002)
- [47] M. van Schiflgaarde, T. Kotani, and S. V. Faleev, Phys. Rev. B 74, 245125 (2006)
- [48] M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (2006)
- [49] P. Rinke, A. Qteish, J. Neugebauer, C. Freysoldt, and M. Scheffler New J. Phys. 7 126 (2005)
- [50] G. Onida, L. Reining, and A. Rubio Rev. Mod. Phys. 74, 601 (2002)
- [51] J. Paier, M. Marsman, K. Hummer, G. Kresse, I. C. Gerber, and J. G. Angyan J. Chem. Phys. 124, 154709 (2006)
- [52] A. Stan, N. E. Dahlen, and R. van Leeuwen, Europhys. Lett. 76, 298 (2006)
- [53] J. B. Neaton, M. S. Hybertsen, and S. G. Louie, Phys. Rev. Lett. 97, 216405 (2006)
- [54] J. M. Garcia-Lastra et al. Phys. Rev. B 80, 245427 (2009).
- [55] K. S. Thygesen and A. Rubio Phys. Rev. Lett. 102, 046802 (2009)
- [56] K. Kaasbjerg and K. Flensberg Nano Lett. 8, 3809 (2008)
- [57] Quek, S. Y. et al. Nano Lett. 7, 3477 (2007)
- [58] Quek, S. Y. et al. Nano Lett. 9, 3949 (2009)
- [59] Y. Meir, N. S. Wingreen Phys. Rev. Lett. 65, 2512 (1992).
- [60] K. S. Thygesen, Phys. Rev. B 73, 035309 (2006)
- [61] H. Ness, L. K. Dash, and R. W. Godby, Phys. Rev. B 82, 085426 (2010)
- [62] A. H. Larsen, M. Vanin, J. J. Mortensen, K. S. Thygesen, and K. W. Jacobsen Phys. Rev. B 80, 195112 (2009)
- [63] J. Enkovaara et al. J. Phys.: Condens. Matter 22, 253202 (2010).
- [64] M. Strange. et al. J. Chem. Phys. 128, 114714 (2008).
- [65] K. S. Thygesen and K. W. Jacobsen Phys. Rev. B 72, 033401 (2005)
- [66] M. Walter, H. Häkkinen, L. Lehtovaara, M. Puska, J. Enkovaara, C. Rostgaard, and J. J. Mortensen, J Chem. Phys. 128, 244101 (2009)
- [67] F. Aryasetiawan and O. Gunnarsson, Phys. Rev. B 49, 16214 (1994)
- [68] G. Baym Phys. Rev. 127, 1391 (1962)
- [69] M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev. Lett. 96, 226402 (2006)
- [70] P. Perdew,K. Burke, M. Ernzerhof M. Phys. Rev. Lett. 77, 3865-3868 (1996).
- [71] For BDA the distance between the Au tip atoms was 9.90 Å, and the Au-N bondlength was 2.55 Å. For BDT+H the distance between the Au tip atoms was 9.80 Å and the Au-S bondlength was 2.47 Å. The BDT was adsorbed on the Au(111) surface in a fcc-hollow position slightly tilted towards the bridge site. The distance between the Au(111) surfaces was 10.20 Å.
- [72] D. J. Lavrich et al. J. Phys. Chem. B 102, 3456 (1998); W. Andreoni et al. Int. J. Quant. Chem. 80, 598 (2000); P. Maksymovych et al. J. Am. Chem. Soc. 130, 7518 (2008); K. Stokbro et al. Comput. Mater. Sci. 27, 151-160 (2003); J. Tomfohr and O. F. Sankey, Chem. Phys. 120, 1542-1554 (2004).
- [73] T. Shimazaki and K. Yamashita, Int. J. Quant. Chem. 109, 1834 (2009)
- [74] K. Stokbro et al. Comput. Mater. Sci. 27, 151-160 (2003); J. Tomfohr and O. F. Sankey, Chem. Phys. 120, 1542-1554 (2004); M. Strange. et al. J. Chem. Phys. 128, 114714 (2008); S. V. Faleev, F. Leonard, D. A. Stewart, and M. van Schilfgaarde, Phys. Rev. B 71, 195422 (2005); F. Evers, F. Weigend, and M. Koentopp, Phys. Rev. B 69, 235411 (2004); E. G. Emberly and G. Kirczenow, Phys. Rev. Lett. 91, 188301 (2003); H. Kondo, H. Kino, J. Nara, T. Ozaki, and T. Ohno, Phys. Rev. B 73, 235323 (2006); X. Qian, J. Li, and S. Yip, Phys. Rev. B 82, 195442 (2010)
- [75] S. Garcia-Gil, A. Garcia, N. Lorente, and Pablo Ordejon, Phys. Rev. B 79, 075441 (2009)
- [76] H. B. Michaelson, J. Appl. Phys. 48, 4729 (1977)
- [77] M. Dell’Angela, et al. Nano Lett. 10, 2470-2474 (2010).
- [78] M. Rohlfing, arXiv:1008.3492
- [79] D. M. Newns, Phys. Rev. 178, 1123 (1969).
- [80] K. S. Thygesen and K. W. Jacobsen Chem. Phys. 319, 111-125 (2005)