First-principles many-body calculations of electronic conduction in thiol- and amine-linked molecules

0.68 MB·English

by  M. Strange

First-principles many-body calculations of electronic conduction in thiol- and
amine-linked molecules
M. Strange,1,2 C. Rostgaard,1 H. Ha¨kkinen,2 and K. S. Thygesen1,∗
1Center for Atomic-scale Materials Design, Department of Physics,
Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark
2Departments of Physics and Chemistry, University of Jyvaskyla, Jyvaskyla, Finland
(Dated: January 28, 2011)
1 Theelectronicconductanceofabenzenemoleculeconnectedtogoldelectrodesviathiol,thiolate,
1 andaminoanchoringgroupsiscalculatedusingnonequilibriumGreenfunctionsincombinationwith
0 thefullyselfconsistent GW approximation. Thecalculated conductanceof benzenedithiolandben-
2 zenediamineisfivetimeslowerthanpredictedbystandarddensityfunctionaltheory(DFT)invery
n good agreement with experiments. In contrast, the widely studied benzenedithiolate structure is
a foundtohaveasignificantlyhigherconductanceduetotheunsaturatedsulfurbonds. Thesefindings
J suggest that more complex gold/thiolate structures where the thiolate anchors are chemically pas-
7 sivated by Au adatoms are responsible for the measured conductance. Analysis of the energy level
2 alignmentobtainedwithDFT,Hartree-FockandGWrevealstheimportanceofself-interactioncor-
rections (exchange) on the molecule and dynamical screening at the metal-molecule interface. The

] maineffectoftheGWself-energyistorenormalizethelevelpositions, however,itsinfluenceonthe
ll shape of molecular resonances also affects the conductance. Non-selfconsistent G0W0 calculations,
a startingfromeitherDFTorHartree-Fock,yieldconductancevalueswithin50%oftheselfconsistent
h
GW results.
-
s
e PACSnumbers: 73.63.-b,73.40.Gk,85.65.+h
m
t. I. INTRODUCTION reported for the gold/benzenediamine junction2,6. This
a
was,however,not confirmed by independent break junc-
m
tion experiments in vacuum16.
The problem of first-principles calculation of elec-
-
d tronic conduction in molecular systems is of longstand- Theuncertaintiesrelatedtothejunctionatomicstruc-
n ing interest. Over the last decade, advances in exper- turerenderstheoreticalbenchmarkcalculationsmoreim-
o
imental techniques have allowed for fundamental stud- portant and more challenging at the same time. Here
c
[ ies of electron transport through few or even individ- progress has been hampered by the inability of conven-
ually contacted molecules. The transport mechanisms tional density functional theory (DFT) methods, which
1 observed in molecular junctions range from ballistic1 for long have been the workhorse and the state of
v over off-resonant tunneling2–6 to the strong correlation the art for quantum transport calculations17,18, to re-
9
Kondo and Coulomb blockade regime7,8 to vibration as- produce the conductance measured for even the sim-
0
3 sisted hopping9. The former two belong to the phase- plest molecular tunnel junctions19,20. As a consequence,
5 coherent transport regime characteristic of relatively the most succesful studies have focused on qualita-
. short molecules (as opposed to molecular wires9) with tive trends in the dependence of conductance on e.g.
1
0 good”chemical”contacttotheelectrodesandisthemain molecularconformation21,molecularlength22,sidegroup
1 focus of the present work. functionalizations23, or have focused on properties inde-
1 pendent of the numerical value of the conductance like
The lack of control over the atomistic details of the
: molecular vibrations24,25. The main shortcoming of the
v metal-molecule interface introduces a strong statistical
DFTapproachhasbeenattributedtoitsbandgapprob-
i element in measurements on single-molecule junctions
X lem,i.e. thefactthatDFTtendstounderestimateenergy
which masks the relation between atomic structure and
r gaps26,andthereforeoverestimatestheconductance. At-
a measured conductance. For example, published conduc-
tempts to overcomethis problemwithin a single-particle
tance values for the gold/benzenedithiol model system
framework have mainly focused on self-interaction cor-
vary by several orders of magnitude, although recent in-
rection schemes27–29.
dependent studies of this system seem agree on a value
around 0.01G (G = 2e2/h)3–6. It is, however, not The well known success of the many-body GW
0 0
clearwhichstructureisresponsibleforthis”typical”con- method30 for quasiparticle (QP) bandstructure calcu-
ductance and recent experimental and theoretical work lations has recently inspired its application to (simpli-
point to complex gold/thiolate structures involving two fied) transport problems34–41. The fact that the GW
molecules bonding the same Au adatom10–15. Scanning approximation succesfully describes systems with highly
tunnelingmicroscopeexperimentsinsolutionhaveshown diverse screening properties ranging from metals42 over
thatthe use ofamino ratherthanthiolanchoringgroups semiconductors44–50 to molecules32,51 is essential for a
leadstomorewelldefinedjunctionproperties2andacon- correct description of metal-molecule interfaces where
ductance around (or just below) 0.01G has also been the electronic character changes from metallic to in-
0
2
sulating over a few angstrom. In particular, screen- anddiscusstheminrelationtothepossibleatomicstruc-
ing by the metal electrons can have large influence on ture of the junctions, we stress that the main purpose of
the QP energies of the adsorbed molecule52–55 – an ef- thisstudyisthebenchmarkingofquantumtransportcal-
fect completely missed by both local and hybrid den- culations for specific, idealized junction with particular
sity functionals53. This has recently motivated the focus on the role of electronic correlationeffects.
use of semi-empirical schemes for correcting the DFT
eigenvalues by a scissors operator prior to transport
calculations23,56,57. While such schemes can be justified II. METHOD
for weakly coupled molecules, they become uncontrolled
in the relevant regime of covalently bonded molecules
We consider a quantum conductor consisting of a
where the screening effects mix with charge transfer and
molecule connected to left (L) and right (R) electrodes.
hybridization54.
We shallassume thatoutside a certainregioncontaining
In this work we combine nonequilibrium Green func- the molecule and part of the electrodes (the ”extended
tionmethodsforelectrontransportwiththefullyselfcon- molecule”), the electron-electron interactions can be de-
sistent GW approximation for exchange and correlation scribed by a mean field potential. The current through
to establish a theoretical benchmark for the electronic the molecule is then given by58,59
structure and conductance of gold/1,4-benzenedithiolate
(BDT), -benzenedithiol (BDT+H) and -benzenediamine i
I = Tr[(Γ Γ )G<+(f Γ f Γ )(Gr Ga)]dE
(BDA) molecular junctions. We find a conductance of 4π Z L− R L L− R R −
0.0042G for BDA and 0.010G for BDT+H in very (1)
0 0
good agreement with experimental data. In compari- where the energy dependence of all quantities has been
son, the conductance obtained from DFT is about five suppressedforsimplicity. InthisequationGistheGreen
times higher while non-selfconsistent G W calculations function matrix of the extended molecule evaluated in
0 0
produce conductances within 50% of the selfconsistent a localized basis, Γ = i[Σr Σa ] is the cou-
L/R L/R − L/R
GW result. We argue that the BDT+H structure can pling strength between the extended molecule and the
be viewed as a simple model of recently proposed RS- left/right electrode, and f are the Fermi Dirac dis-
L/R
Au(I)-SRgold/thiolatestructuresinvolingtwomolecules tribution functions of the two leads. Our implemen-
attached to the same Au adatom. The conductance of a tation applies to the general case of a finite bias volt-
simple BDT molecule between Au(111) surfaces is pre- age, but in this work we focus on the zero bias conduc-
dicted to be on the order of 1G0 by both DFT, Hartree- tancewhichcanbeexpressedintermsofthetransmission
Fock and GW. The origin of the high conductance is function58,59
due to an unsaturated sulfur p orbital with energy just
below the Fermi energy. In the BDT+H and RS-Au(I)- T(E)=Tr[Gr(E)Γ (E)Ga(E)Γ (E)] (2)
L R
SR structures, the sulfurs are fully passivated and the
p-orbitalmovesawayfromE leadingto aneffective de- as G = G T(E ). This formula was originally derived
F 0 F
coupling of the C H moiety from the gold electrodes. fornon-interactingelectrons,butisinfactvalidforinter-
6 4
We find that the main effect of the GW self-energy is to acting electrons in the low-bias limit60. We have verified
shift the molecular levels and can be modelled by a sim- thatthis is indeed fulfilled to highnumericalaccuracyin
ple scissors operator. However, the energy dependence our GW calculations by comparing the conductance ob-
of the GW self-energy may also affect the shape of the tained from Eq. (1) for small finite voltages with T(E )
F
transmissionresonancesandthiscanchangetheconduc- evaluated in equilibrium.
tance by almost a factor of two. The retardedGreenfunction ofthe extendedmolecule
Most implementations of the GW method invoke one is calculated from
orseveraltechnicalapproximationsliketheplasmonpole
approximation, neglect of off-diagonal matrix elements Gr(E)=[(E+iη)S H0+Vxc ∆VH[G]
− −
in the GW self-energy, analytic continuations from the Σ (E) Σ (E) Σ [G](E)]−1 (3)
L R xc
− − −
imaginary to the real frequency axis, neglect of core
states contributions to the self-energy, neglect of self- where η is a positive infinitesimal and
consistency. The range of validity of these approxima-
tionshasbeenexploredforsolidstatesystemsbyanum- Sij = φi φj (4)
h | i
berofauthors42,44–48,however,muchlessisknownabout 1
H = φ 2+v (r)+v (r)+v (r)φ (5)
their applicability to molecular and metal-molecule sys- 0,ij h i|− 2∇ ion H xc | ji
tems. For this reason our implementation of the GW V = φ v (r)φ (6)
xc,ij i xc j
method avoids all of these approximations and as such h | | i
represents an exact treatment of the GW self-energy denotethe overlapmatrix,Kohn-ShamHamiltonianand
within the space of the employed atomic orbital basis exchange-correlation potential, respectively. The matri-
set. ces are evaluated in terms of a basis consisting of nu-
Althoughwecompareourresultstoexperimentaldata merical atomic orbitals61, and are obtained from a DFT
3
supercell calculation performed with the real-space pro- where the factor 2 accounts for spin. Setting P = 0
jector augmented wave method GPAW62. The electrode yields the Hartree-Fockapproximationwhichthus corre-
self-energies Σ are obtained from the Kohn-Sham sponds to complete neglect of screening or equivalently
L/R
Hamiltonian of a bulk DFT calculation using standard complete neglect of correlations. Note that Eqs. (9-10)
techniques63. involve time-ordered quantities defined on the Keldysh
The boundary conditions in the plane normal to the time contour. For completeness we provide the expres-
transport direction enter only via the electrode self- sions for the real time components in Appendix A and
energieswhichareconstructedfromtheelectrodesurface refer the reader to Ref. 35 for more details.
Greenfunction63. The lattershouldrepresentaninfinite
surface but is here approximated by that of a periodic
supercell with 4 4 Au atoms in the surface plane. We
×
have found that this is a very good approximationwhen
thesurfaceGreenfunctionofthe4 4cellisevaluatedat
×
a generalk -point(we use k =(0.1,0.2)in coordinates
⊥ ⊥
of the surface Brillouin zone basis vectors). Using high
symmetrypoints,inparticulartheGammapoint,canbe
problematic64.
Theterm∆V isthedeviationoftheHartreepotential
H
FIG.1: (Coloronline)Schematicofamolecularjunctionwith
from the groundstate DFT Hartree potential contained
the ”molecule” and ”extended molecule” regions indicated
in H ,
0 (small and large box, respectively). The GW self-energy (Σ)
is calculated on the molecule while the polarization (P) and
∆VH,ij =2 vij,kl(̺kl[G]−̺0kl[G0]) (7) screened interaction (W) are evaluated for the extended re-
Xkl giontoensurepropertreatmentofnon-localscreeningeffects.
The electrodes and electrode-molecule coupling is described
In this expression ̺ = i G<(E)dE and ̺0 = attheDFTlevel. AsecondorderFeynmandiagramisshown.
i G<(E)dE are the intera−ctiRng and Kohn-Sham den-
− 0
sityRmatrices, respectively, and As explained above, the matrices P, W, and Σ are
GW
calculatedforthe extendedmolecule. Onthe otherhand
φ∗(r)φ (r)φ (r′)φ∗(r′)
v = i j k l drdr′ (8) it is clear that W, and thus ΣGW, should have contribu-
ij,kl Z Z r r′ tions from polarization diagrams outside this region, see
| − |
Fig. 1. Physically these diagrams describe the poten-
is the bare Coulomb interaction in the atomic orbital tial acting on an electron propagating on the molecule
basis. The factor 2 is due to spin. Ref. 65 describes how due to the polarization that it induces in the lead. For
the all electron Coulomb elements are obtained within this reason the self-energy will not be fully convergedat
the PAW formalism. The last term in Eq. (3) is the the ends of the extended molecule region. To overcome
many-bodyexchange-correlationself-energywhichinthis this problem we only use the part of Σ correspond-
GW
workcanbeeitherthebareexchangepotentialortheGW ing to the molecule and replace the remaining parts of
self-energy. Wenotethatsetting∆VH =0andΣxc =Vxc the xc self-energy of the extended molecule by the DFT
inEq. (3)weobtaintheKohn-ShamGreenfunction,G0. xc-potential. Symbolically,
v v v
xc xc xc
Σ (E)= v Σ (E) v , (11)
A. GW self-energy xc xc GW xc
v v v
 xc xc xc 
The (time-ordered) GW self-energy is given by We stress that although we only include ΣGW on the
molecule, the interactions between electrons on the
Σ (τ)=i G (τ+)W (τ) (9) molecule and electrons in the electrodes (leading e.g. to
GW,ij kl ki,lj
X imagechargerenormalizationofthemolecularlevels)are
kl
included via diagrams of the form shown in Fig. 1. We
where W is the screened interaction and the indices alsonotethattheform(11)impliesthatallmetalatoms,
i,j,k,l run over the atomic basis functions of the ex- boththoseinside andoutsidethe extendedmolecule,are
tended molecule. The screened interaction is calculated consistently described at the same (DFT) level. For the
in the frequency domain as the matrix product W(ω)= size of the extended molecule we have found it sufficient
ǫ−1(ω)v with the dielectric function, ǫ(ω) = 1 vP(ω), to include the goldatomswhicharenearestneighborsto
−
evaluatedin the randomphase approximation. The irre- the sulfur or nitrogen atoms, i.e two gold atoms for the
ducible densityreponsefunctioniscomputedinthe time tip structures and six for the flat structure depicted in
domain, Fig. 2,seeAppendix B. We expectthatthis ratherlocal
screening response is special for covalent metal-molecule
P (τ)= 2iG (τ)G ( τ). (10) bonds.
ij,kl ik lj
− −
4
B. Time/frequency dependence bution to the valence exchange self-energy coming from
the core electrons. As the density matrix is simply the
Thetime/frequencydependenceofG,P,W,andΣ identitymatrixinthesubspaceofatomiccorestates,this
GW
isrepresentedonauniformgridrangingfrom 200to200 valence-core exchange reads
−
eV with a grid spacing of ∆ω = 0.01eV. We have veri-
core
fied that the results are coverged with respect to both Σcore = v , (13)
the size and spacing of this grid. Fast Fourier trans- x,ij −X in,jn
n
formis usedto switchbetweenenergyandtime domains
where i,j represent valence basis functions and n rep-
to avoid convolutions. The calculations are parallelized
resent atomic core states. This contribution is added
both over basis functions and over the time/frequency
to Σ describing the valence-valence interactions. We
grid points. One should always have η ∆ω to ensure GW
≥ limit the inclusion of valence-coreinteractions to the ex-
properrepresentationofpossibleboundstates. However,
changepotential, neglecting it inthe correlation. This is
we havefoundthatthe conductance,andmoregenerally
reasonable, because the polarization bubble, P, involv-
the transmission function at any given energy, can be
linearly extrapolated to the η = 0+ limit. This extrapo- ing core andvalence states will be smalldue to the large
energydifferenceandsmallspatialoverlapofthe valence
lationhas beenperformedfor allthe resultspresentedin
and core states. In generalwe have found that the effect
this work.
of Σcore on molecular energy levels can be up to 1 eV32.
The memory requirements for the GW calculations x,ij
For the benzene-like molecules considered in this work
(defined mainly by the size of the P and W matrices)
the effect is generally less than 0.4 eV for the frontier
are approximately a factor 50 larger than for a corre-
orbitals.
sponding 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
E. Self-consistency
2 hours per selfconsisistency iteration. In comparison a
DFTcalculationforthesamesystemtookaround5hours
on 8 cores. SinceΣGW and∆VH dependonG,theDysonequation
(3) must be solved self-consistently in conjunction with
the self-energies. In practice, this self-consistency prob-
C. Product basis lem is solved by iteration. We have found that a linear
mixing of the Green functions,
ThecalculationofalloftheCoulombmatrixelements,
Gn(E)=(1 α)Gn−1(E)+αGn (E) (14)
vij,kl, is prohibitively costly for larger basis sets. For- in − in out
tunately the matrix is to a large degree dominated by with α = 0.15 generally leads to selfconsistency within
negligible elements. To systematically define the most 20-30 iterations.
significant Coulomb elements, we use the product basis Fully selfconsistentGW calculationsarenotstandard,
technique of Aryasetiawan and Gunnarsson66. In this and in fact only few previous calculations of this type
approach, the pair orbital overlap matrix havebeenreported32,42,51. ConventionalGWbandstruc-
ture calculations typically apply a one-shot technique
S = n n , (12)
ij,kl ij kl
h | i where the self-energyis evaluatedwith a non-interacting
where nij(r) = φ∗i(r)φj(r) is used to screen for the sig- Green function, G0, usually obtained from DFT43. In
nificant elements of v. comparison, the selfconsistent approach has the imme-
The eigenvectors of the overlap matrix Eq. (12) rep- diate advantage of removing the G -dependence, i.e. it
0
resents a set of “optimized pair orbitals” and the eigen- leads to a unique solution.
values their norm. Optimized pair orbitals with insignif- For an approach, like the present, where the chemical
icantnormmustalsoyielda reducedcontributionto the potential is fixed by the external boundary conditions,
Coulomb matrix, and are omitted in the calculation of some kind of selfconsistency (though not necessarily the
v. We have found that the basis for v can be limited to fullGWselfconsistencyemployedhere)isessentialtoen-
optimizedpairorbitalswitha normlargerthan10−5a−3 sure charge neutrality. This is particularly important
0
without sacrificing accuracy. This gives a significant re- for cases where a molecular resonance lies close to the
duction in the number of Coulomb elements that needs chemical potential. A shift in the energy of such a reso-
to be evaluated, and it reduces the matrix size of P(ω) nancecouldleadtoalargechangeinitsoccupation. Ina
and W(ω) correspondingly, see Appendix A. selfconsistent calculation this shift would be counterbal-
anced, mainly by a change in the Hartree potential. On
the other hand, the one-shot approach does not account
D. Valence-core exchange forthiseffectandunphysicallevelalignmentscouldoccur
as a result.
Since both core and all-electron valence states are Finally,thefullyselfconsistentGWapproximationisa
availableinthePAWmethod,wecanevaluatethecontri- conserving approximation in the sense of Baym67. This
5
becomes particularly important in the context of trans- stretched by 1 ˚A. We stress that our calculations do not
port where it ensures that the continuity equation is includeeffectsofentropywhichbecomesrelevantatfinite
satisfied35,67. We mention that the recently introduced temperature, and furthermore assumes that hydrogenin
quasi-selfconsistentGWmethod(nottobeconfusedwith the gas phase is the proper reference for the solvated
the fully selfconsistent GW approximation used here), protonand an electron in the electrode, i.e. the reaction
in which G is chosen such as to mimick the interact- H++e− H (g) is in equilibrium (in electrochemistry
0 2
↔
ing Green function as closely as possible, have shown language we assume the standard hydrogen potential at
that selfconsistency in general improves the band gaps pH=0). For these reasons our calculations are not suffi-
of semiconductors as compared to standard one-shot cient to address the relative stability of thiols vs. thio-
calculations.68 lates under the relevant experimental conditions.
Importantly we note that new experimental evidence
for the chemical structure of the gold-thiolate inter-
III. RESULTS face at the Au(111) surface10–12 or at Au nanocluster
surfaces13,14 has recently emerged, pointing to the ex-
In this section we discuss the results of our self- istence of polymeric SR-Au(I)-SR units where the for-
consistent GW calculations for the electronic struc- mally oxidised Au(I) adatoms are chemically bound to
ture and conductance of the prototypical gold/1,4- thiolates and form a part of a more complex structure
benzenedithiolate (BDT), -benzenedithiol (BDT+H), (see Ref. 15 and references therein). These complexes
and -benzenediamine (BDA) molecular junctions. We arecurrentlytoochallengingtotreatsatisfactorilyatthe
argue that the thiol structure can be considered as a GW level. However, we have found that the electronic
simple model for more complex gold/thiolate structures structureofsuchcomplexesisquitesimilartothatofthe
which have been proposed recently15 but which are cur- dithiol structure, see Sec. IIIG. This is because the hy-
rently too large to be treated satisfactorily at the GW drogenatomsplayarolesimilartothatoftheAuadatom
level. The transportresults are rationalizedby consider- in passivating the sulfur atoms. Therefore the transport
ing the alignment of molecular energy levels in the junc- properties of the dithiol structure should also be similar
tion. Here we show that both DFT and Hartree-Fock to those of the SR-Au(I)-SR units.
providequantitativelyandqualitativelywrongresultsby
predictingagapopeningratherthanreductionwhenthe
moleculeisattachedtoelectrodes. Finally,weinvestigate B. Energy levels of isolated molecules
to what extent the GW results can be reproduced by a
simplescissorsoperatorappliedtotheDFTHamiltonian. A natural requirement for a method intended to de-
scribetheenergylevelsofmoleculesincontactwithelec-
trodes, is that it should be able to describe the energy
A. Junction geometries levels of isolatedmolecules. As we show below,the DFT
approachfailscompletelyinthisrespectunderestimating
The junction geometrieswereoptimized using the real theHOMO-LUMOgapsofthethreemoleculesby5-6eV
space projector augmented wave method GPAW62. We while GW energies lie within 0.5 eV ofthe targetvalues.
used a grid spacing of 0.2 ˚A and the PBE functional for The gas-phase molecular structures have been relaxed
exchangeandcorrelation69. Themoleculeswereattached in a 16 ˚A cubic cell using GPAW grid calculations as
toAu(111)surfaces,modelledbyasevenlayerthick4 4 describedintheprevioussection. Forconsistency,allen-
×
slab, either directly (in the case of BDT) or via tips (in ergy levels have been calculated using the same double-
the case of BDT+H and BDA) as shown in Fig. 2. The zeta (DZ) atomic orbitals basis set. This is the same
surface Brillouin zone was sampled on a 4 4 k-point basis as used for the molecules in the transport calcula-
×
grid,andthestructuresincluding molecule,Autips, and tions presented in the next section. For the DFT and
outermostAu surface layerswere relaxeduntil the resid- Hartree-Fock calculations we have found that the ener-
ual force was below 0.05 eV/˚A. The three structures are giesobtainedwith the DZ basisagreewith accurategrid
shown in the upper panel of Fig. 2 and some key bond calculations to within 0.2 eV. For the GW calculations
lengths are given here70. the energies are within 0.1 eV of those obtained with a
It is generally accepted that the hydrogen atoms dis- DZ+polarization basis set.
sociatefromthe thiolendgroupsforminga gold-thiolate Table I summarizes the results for the HOMO and
structure71. Nevertheless, our total energy calculations LUMO orbital energies obtained from the DFT-PBE
show that the benzenedithiol structure has a slightly eigenvalues, selfconsistent Hartree-Fock, and selfconsis-
lower energy than the benzenedithiolate when inserted tent GW. The energy levels have been identified as the
between two gold tips as shown in Fig. 2(middle). In peaks inthe spectralfunctionTr[ImGr(E)] extrapolated
these calculations the hydrogens were either taken to be toη =0+. Due to lackofaccurateexperimentaldata we
inthe gas-phaseorareadsorbedonthe Autips. Inboth have also performed DFT-PBE total energy calculations
cases the energy gain is less than 0.1 eV at the equi- for the neutral, cation, and anion species to obtain the
librium distance but increases to 0.3 eV for a junction addition/removal energies (last column). This approach
6
1.0 0.1 0.1 1.0 DFT GW HF
DFT 1.0 DFT
0.01
n n GW n
missio GW 0.001 missio 0.01 HF missio
ans0.5 -1 0 HF1 ans0.5 0.001-1 0 1 ans0.5
Tr Tr Tr
0.0 0.0 0.0
-4 -2 0 2 -4 -2 0 2 -2 -1 0 1 2
E-E (eV) E-E (eV) E-E (eV)
F F F
FIG. 2: (Color online) Atomic structure of the BDA (left), BDT+H (middle), and BDT (right) molecular junctions. The
lowerpanelsshowthetransmissionfunctionscalculatedusingDFT-PBE(red),Hartree-Fock(blue),andtheselfconsistentGW
approximation (green). Theinsetsshowazoom ofthetransmission functionsaroundtheFermienergy(set tozero). Thegrey
boxes indicate the conductance windows 0.007G0±50% and 0.01G0±50% which cover the the experimental values reported
in Refs. 2,6 and 3–6, respectively.
has been shown to produce very accurate estimates of Molecule Orbital DFT-PBE HF GW ∆Etot
the experimental ionization and affinity levels of small
BDA HOMO -4.1 -7.2 -6.2 -6.8
molecules32.
(C6H8N2) LUMO -0.9 3.9 2.9 2.3
Relative to this reference,the DFT eigenvalues under-
H-L Gap 3.2 11.1 9.1 9.1
estimate the HOMO-LUMO gap of all three molecules
BDT+H HOMO -5.1 -8.0 -6.9 -7.5
by 5 6 eV, Hartree-Fock overestimates it by 2 3 eV,
while−thegapobtainedwithGWlieswithin0.3eV−.These (C6H6S2) LUMO -1.3 3.3 2.2 1.3
trends are consistent with a recent study of ionization H-L Gap 3.8 11.3 9.1 8.8
potentials of a large number of molecules32. The main BDT HOMO -5.7 -8.6 -7.9 -8.3
reasonforthelargeunderestimationofthegapbyDFTis (C6H4S2) LUMO -5.1 -1.6 -2.3 -2.7
the presenceofself-interactionsinthe PBEfunctional.28
H-L Gap 0.6 7.0 5.6 5.6
On the other hand Hartree-Fock is self-interaction free;
here,byvirtueofKoopmans’theorem,theoverstimation
TABLEI:Calculatedfrontierorbitalenergiesofthemolecules
ofthe gapcanbe seenasaresultofneglectoforbitalre- inthegas-phase. AllenergiesareineVandmeasuredrelative
laxations. The effect of the latter is included in GW via to the vacuum level. DFT-PBE refers to the Kohn-Sham
thescreenedinteractionandthisreducesthegaprelative eigenvalueswhile∆Etot representsaddition/removalenergies
the unscreened Hartree-Fock result. obtained from self-consistent total energy calculations of the
We furthermore note that DFT places the LUMO of neutral, anion and cations at theDFT-PBE level.
BDA and BDT+H below the vacuum level thus incor-
rectly predicting the anions to be stable. For BDT, the
LUMO level is predicted to be negative by all methods anddouble-zetawithpolarizationfortheAuatoms. The
indicating the radical nature of this species. We note results are shown in the lower panels of Fig. 2, and the
that our GW energies for BDT are in good agreement correspondingconductances are summarizedin Table II.
with previous MP2 calculations72.
TheconductanceofBDAandBDT+Hcalculatedwith
the fully selfconsistent GW approximation agree well
withtheexperimentalvaluesreportedinRefs. 2and6for
C. Conductance calculations benzenediamine and Refs. 3–6 for benzenedithiol as in-
dicatedbythegreyboxesinFig. 2(left+middle). Incon-
The transmissionfunctions of the relaxedjunction ge- trast, DFT and Hartree-Fock respectively overestimates
ometries were calculated as described in Sec. II us- and underestimates the experimental conductances by
ing three different approximations for Σ , namely the factors 5-10. Our DFT result for the BDA junction is
xc
PBE xc-potential, the bare exchange potential (leading in good agreement with previous calculations23,56.
to Hartree-Fock theory), and GW. The former choice It is striking that the conductance of the ”classical”
corresponds to the standard DFT-approach. All calcu- BDT junction shown in Fig. 2(right), is predicted by all
lations employ a double-zeta basis set for the molecules three methods to be significantlyhigherthanthe experi-
7
2
Method BDA BDT+H BDT
)
V
DFT-PBE 2.1·10−2 5.4·10−2 2.8·10−1 e HS SH H 2 N NH2
(
HF 4.0·10−4 2.7·10−3 5.7·10−1 ge 1
GG0WW0(PBE) 48..20··1100−−33 11..60··1100−−22 78..53··1100−−11 p chan 0 HF DFT HF DFT
G0W0(HF) 2.2·10−3 9.8·10−3 8.7·10−1 Ga GW GW
-1
TABLE II: Calculated conductance in units of G0 for the
4
three junctions shown in Fig. 2. The last two rows refer to
non-selfconsistent GW calculations based on the DFT-PBE )
V 2
or Hartree-Fock Green function, respectively. e
(
y Vacuum
g 0
r
mental value (we obtain the same high conductances for ne Φ= 4.9
e -2
BDT between tips, i.e. the structure in Fig. 2(middle) O
M
withouthydrogenonsulfur). OurDFTconductanceisin U -4
E
overallgoodagreementwiththelargenumberofprevious L F
d
calculations for the same or similar similar structure73. n -6
a
The high conductance is clearly due to a strong peak in O
M -8
the transmission function close to the Fermi level. The O Au Au
peak moves to higher energies when going from DFT- H
-10
PBE over GW to Hartree-Fock, opposite to the trend
In gas−phase
normally seen for occupied states. The peak comes from -12 In junction In junction
anunsaturatedp-orbitalonthesulphuratomsandisdis-
cussed in more detail in Sec. IIIG. FIG. 3: (Color online) Upper panel shows the change in the
Theseresultssuggestthatthestructuresprobedinex- HOMO-LUMO gap as the molecules are brought from the
periments on benzenedithiol junctions involve a chemi- gas-phaseintothejunction. Lowerpanelshowstheenergyof
cally passivated form of the thiolate linker group. As theHOMOandLUMOlevels(relativetovacuum)inthegas
we show in Sec. IIIG, the high conductance of BDT phaseandinthejunction. Theleftandrightsetoflevelscor-
is due to unsaturated p-states on the sulphurs with en- respond to BDT+H and BDA, respectively. The calculated
workfunction ofthegold junction(with tips)inabsenceofa
ergy close to the Fermi level (the energy of this orbital
molecule is Φ=4.9 eV as indicated.
does not change considerable with DFT, Hartree-Fock
andGW).InthethiolandSR-Au(I)-SRstructures,these
statesformbondstoHandthe Auadatom,respectively,
and are thereby shifted away from the Fermi level. On functionsthroughoutthiswork61)weobtainaworkfunc-
the other hand, the electronic structure and transmis- tion for Au(111) of 5.4 eV in good agreement with the
sion functions of BDT+H and SR-Au(I)-SR structure experimental value of 5.31 eV75. At the position of the
are rather similar indicating that the BDT+H can be molecule,i.e. intheregionbetweenthetwotips,theelec-
viewedasasimplemodelofthemorecomplexSR-Au(I)- trostaticpotentialfromacalculationwherethemolecule
SR structure. hasbeenremoved,convergestoaconstantvalueof4.9eV
It should of course be kept in mind that experiments abovethemetalFermilevel. Thisvaluehasbeenusedas
are performed in solution and at room temperature and reference for the vacuum level in Fig. 3. In the junction
are subject to variations in the detailed atomic struc- wherethelevelsarebroadenedduetohybridizationwith
ture. Thus the measured conductance values should not the metalstates,the levelpositions havebeen defined as
beconsideredashighlyaccuratereferencesfortheoretical the firstmomentof the projecteddensity ofstates ofthe
calculations on idealized junctions. relevant molecular orbital, Im ψn Gr(E)ψn . Here the
h | | i
ψ are obtained by diagonalizing the KS Hamiltonian
n
| i
within the molecular subspace.
D. Energy level alignment The orbital energies obtained from a GW calculation
include the dynamical response of the electron system
In Fig. 3 we show the calculated HOMO and LUMO to the added electron/hole via the correlation part of
energy levels of BDT+H and BDA in the gas-phase and the self-energy. In general the correlations will shift the
inthejunction. Allenergieshavebeenalignedrelativeto occupied levels up and the empty levels down relative
thevacuumlevel. Atthispointwenotethatanaccurate to the bare Hartree-Fock energies. When a molecule is
description of the vacuum level, i.e. the work function, broughtfromthegas-phaseintoajunctiontheelectronic
caningeneralbedifficulttoobtainwithanatomicorbital character of its environment changes from insulating to
basis74. However, by using more diffuse basis functions metallic. The enhanced screening should thus cause the
(anenergyshiftof0.01eVhasbeenusedforallAubasis gap to shrink (neglecting shifts due to hybridization) as
8
comparedtoitsgas-phasevalue. However,ithasrecently in Fig. 3). As an example Fig. 4 shows the calculated
been shown for molecules weaklybonded to a metal sur- transmission functions for the BDA junction in a larger
face, that this effect is completely missing in effective energy range around E .
F
single-particletheoriesbasedona(semi)localdescription For the BDA and BDT+H junctions, G W [PBE]
0 0
of correlations52,53. overestimates the conductance while G W [HF] under-
0 0
Asaresult,inourDFTandHartree-Fockcalculations, estimates the conductance relative to GW. This can be
thechangeinthefrontierorbitalenergiesinducedbythe understood as follows: Since DFT-PBE (Hartree-Fock)
coupling to gold is completely governed by the effect of underestimates (overestimates) the HOMO-LUMO gap,
hybridization which tends to open the gap by 0.7 eV theuseoftheseGreenfunctionstoevaluatetheGWself-
for both molecules, see top panel of Fig. 3. In contrast, energy will lead to an overestimation (underestimation)
the GW gap is reduced by 1 eV due to the enhanced of the screening. As discussed above the screening con-
screening in the contact. Since the hybridization shift is tainedinthecorrelationpartoftheGWself-energytend
ofcoursealsopresentinhe GWcalculation,we conclude to reduce the HOMO-LUMO gap. This reduction of the
that the enhanced screening due to the metal contacts gapis thus overestimatedin the G0W0[PBE]calculation
reduces the HOMO-LUMO gap by 1.7 eV relative to and underestimated in the G0W0[HF] calculation which
the value in the gas-phase. explains the trends in conductance, i.e. more screening
Note that we refrain from using the term ”image smaller gap higher conductance.
→ →
charge effect” to describe the renormalization of molec-
ular orbitals. This term is appropriate for weakly cou-
1
pled molecules where the screening of the added elec-
tron/holetakesplacewithinthe metal. Forintermediate G W (DFT) 0.01
0 0
or strongly coupled molecules, there is no clear distinc-
n GW
tion between metal states and molecular states, and the o
screening is more appropriately described as dynamical si
s 0.005
charge transfer54. mi
From Fig. 3 it follows that the HOMO level of the s0.5
n
molecules in the junction is predicted by DFT-PBE to a G W (HF)
r 0 0
T
lie 0.5 0.7 eV higher than obtained with GW. This -0.4 0 0.4
−
agreeswellwitharecentstudyofbenzenediaminederiva-
tives on gold(111) which showed that DFT places the
HOMO level about 1 eV too high compared to UPS
measurements76. The fact that the DFT-PBE descrip- 0 -4 -3 -2 -1 0 1 2
tion of the energy levels is better for the adsorbed E-E (eV)
F
ratherthan isolatedmolecules may be seenasa resultof
the metallic screening build into the DFT xc-functional
FIG. 4: (Color online) Transmission functions for the
via its origin in the homogeneous electron gas77. It
BDA junction calculated with selfconsistent GW and non-
should also be noted that the error (compared to GW)
selfconsistent G0W0 using either the Hartree-Fock or the
of the DFT eigenvalues is significantly larger for the DFT-PBE Green functions as input.
LUMO than for the HOMO. This is in good agreement
withpreviousplanewavecalculationsformolecule/metal
The G W results for the strongly coupledBDT junc-
interfaces53. 0 0
tionshowlessvariation. Thisisperhapssurprisinggiven
the large difference between the DFT-PBE and Hartree-
Fock results shown in Fig. 2. However, DFT-PBE and
E. G0W0 calculations Hartree-Fock give almost equal density if states at the
Fermi level (the transmission functions at the Fermi en-
To test the role of selfconsistency (in the GW ergy are also rather similar), and therefore the screening
and Hartree self-energies) we have performed non- contributionobtainedwith the two choicesfor G0 is also
selfconsistent G W calculations using both the DFT- very similar.
0 0
PBE Green function and the (selfconsistent) Hartree-
Fock Green function as the initial G . The results for
0
theconductanceofallsystemsaresummarizedinthelast F. Scissors operator
two rows of Table II. We notice is that the conductance
can vary by a factor of three depending on G . We also In this section we investigate to what extent the GW
0
note that the Hartree-Fock starting point comes closer transmission function can be reproduced by a DFT cal-
to the selfconsistent result. This is because the Hartree- culation where the occupied and unoccupied molecular
Fock Green function is an overall better approximation levels have been shifted rigidly to match the GW ener-
to the GW Green function than is the DFT Green func- gies. This is illustrated by applying a scissors operator
tion (see e.g. the comparison of frontier orbital energies (SO) to the Au-BDA-Au junction shown in Fig. 2(left).
9
After shifting the molecular levels we perform a non- levelenergies can explain the main part of the difference
selfconsistent calculation of the transmission function. between the DFT and GW conductance, the different
A similar procedure has previously been successfully shape of the transmission resonances also plays a role.
used to include image charge effects and correct for self- This is clear from the inset which shows the GW trans-
interaction errors in DFT-transport calculations23,56,57. missionfunction (green)andthe DFT transmissionwith
We refer the reader to Ref. 23 for more details on the SO chosen to match the GW HOMO and LUMO levels
SO technique applied here. (the dot in the main figure). The lower conductance ob-
tained with GW is seen to be a consequence of a faster
decay of the HOMO resonance towards the Fermi level.
Thisdifference comesfromthe energydependence ofthe
4.0
10 GW self-energy.
DFT-SO (-0.6, 3.8) 7.2e-03
3.5 on 1 GW
si
s
3.0 mi 0.1
V) ans G. Thiol vs. thiolate structures
of LUMO (e22..05 Tr0.01-3 -2E-E-1F (eV0) 1 e-03 4e-03 2.0e-02 tioInns Ffoirg:. (a6)wtheec“ocmlapsasriecatl”hestDruFcTturteraonfsmbeinsszieonnedfuitnhci--
hift 1.5 4.0 6. olate [structure in Fig. 2(right)] (b) benzenedithiol be-
S tween tips on Au(111) [structure in Fig. 2(middle)], and
1.0
(c) benzenedithiolate in a SR-Au(I)-SR molecular unit
0.5 form [structure 2 of Fig. 1 in Ref. 15]. The conductance
(essentially the transmission function at the Fermi en-
0.03.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0.5 1.0 ergy) obtained for structures (b) and (c) is very similar
(cid:0) (cid:0) (cid:0) (cid:0)Shift(cid:0) of HO(cid:0)MO (e(cid:0)V) but markedly different from (a). The high conductance
ofthestructure(a)isduetothe strongtransmissionres-
FIG. 5: (Color online) DFTconductance of theAu-BDA-Au onance lying just below the Fermi level. This resonance
junction for different shifts of the occupied and unoccupied has also been found in many previous studies and seems
molecular orbital energies. The three iso-contours with val-
to be a characteristic and robust feature of this junc-
ues of 4.2x10-3, 6.4x10-3 and 2.1x10-2 in units of G0 corre-
tion. As shown in this work this transmission resonance
spondtotheGW,experimentalandDFTconductancevalues,
is also present in GW and Hartree-Fock calculations. In
respectively. The dot indicates the conductance obtained by
contrast,for both structure (b) and (c) the transmission
fittingtheDFTHOMOandLUMOlevelpositionstothecor-
responding GW level positions. Inset: Transmission function function is rather flat in an energy window of 1 eV
±
calculated with GW (green) and DFT+SO (red). The SO around the Fermi level.
shiftsare-0.6eVand3.8eVfortheoccupiedandunoccupied To examine the origin of the different transmission
molecular orbital, respectively, which makes the HOMO and functions found for structures (a) and (b,c) we found it
LUMOenergiesoftheDFTcalculationcoincidewiththeGW usefultoconsiderthemolecularorbitalsoftheC H moi-
6 4
energies.
ety of the three molecules. By diagonalizing the (Kohn-
Sham)Hamiltonianofthispartofthemoleculeswefound
In Fig. 5 we show contour plots of the DFT conduc- thatthe orbitaldepictedto theleftinthe upperpanelof
tance for the BDA junction where the shift of the occu- Fig. 7, and which constitute the HOMO-1 of the C6H4
pied and unoccupied molecular orbitalenergies has been moiety,isresponsibleforallthestructureinthetransmis-
varied independently over4 eV. The three values for the sion functions below the Fermi level. Note, that the or-
contour lines shown correspond to the conductance ob- bitalsobtainedinthis wayaredifferentfromthe HOMO
tained with GW (0.0042G ), the experimental conduc- and LUMO levels shown in Fig. 3 which were obtained
0
tance (0.0064G ) and the DFT conductance (0.021G ), bydiagonalizingtheHamiltonianfortheentiremolecules
0 0
respectively. The dot indicates the energy shifts which including the SH and NH2 end groups.
make the HOMO and LUMO levels of the DFT calcu- The different panels of Fig. 7 show the projected den-
lation match the corresponding levels of the GW cal- sityofstates(PDOS)oftheHOMO-1forthethreestruc-
culation; the required SO shifts are -0.6eV and 3.8eV tures together with the PDOS of the sulfur p-orbital
for the occupied and unoccupied molecular orbitals, re- to which the HOMO-1 couples. Within the Newns-
spectively. Shifting the levels by this amount reduces Anderson model78, the sulphur p-orbital is the so-called
the DFT conductance by a factor of 3 from 0.02G to group orbital. Note, that the PDOS of the p-orbital
0
0.0072G . Interestingly, the GW conductance is not re- has been calculated in the absence of coupling to the
0
produced by these shifts; it is even lower by a factor of HOMO-1 as it should in the Newns-Anderson model.
1.5. In fact, to reproduce the GW conductance a shift The transmission function is then essentially the prod-
of about -1.3eV of the DFT HOMO is required (keep- uct of the PDOS of the HOMO-1 and the PDOS of the
ing the LUMO position fixed at the GW position). This grouporbital79. Itisclearthattheoriginofthe(double)
shows that while the renormalization of the molecular transmission peak around -1 eV for structure (a) is due
10
a) b) c) strong peak around -1 eV due to the sulphur p-orbital.
Although the above picture is based on the DFT elec-
tronic structure, the qualitative similarities of the DFT
and GW transmission functions in Fig. 2 implies that
the same picture applies to the GW calculations.
a)
1.5
0.4
n
missio 1 a b -1eV) 0 b)
(
ns S 0.4
a O
Tr0.5 c D
P 0
c)
0.4
0
-4 -3 -2 -1 0 1 2 3
E-E (eV)
F 0
-8 -6 -4 -2 0 2 4
FIG. 6: (Color online) DFT-PBE transmission functions for E-E (eV)
F
threedifferentstructuresoftheAu-BDT-Aujunction: a)the
”classical” benzenedithiolate, b) benzenedithiol, and c) the
FIG. 7: (Color online) Projected density of states for the
SR-Au(I)-SRcomplex (with two benzene molecules replaced
HOMO-1 of the C6H4 moeity and its group orbital (sulfur p
byCH3 unitsforsimplicity). Thetransmission functionofb) orbital). Theseorbitalsessentiallydeterminethetransmission
and c) are rather similar while that of a) shows a peak close
functionaroundtheFermienergyforallthethreestructuers.
to theFermi level.
It should be noticed that for structure (a) (the ”classical”
BDT structure) the sulfur p-orbital has a peak in the PDOS
just below the Fermi level which is responsible for the high
to a resonance formed by the HOMO-1 and the sulfur conductance. In(b)and(c)thesulfurischemicallypassivated
and the PDOS of the p-orbital splits into bonding and anti-
p-orbital. The chemical passivation of sulfur, by either
bondingstates at ±3eV therebylowering the conductance.
hydrogen or the Au adatom, implies that the PDOS of
the p-orbitalsplits 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 particu-
lar its magnitude around the Fermi level is lowered and IV. SUMMARY
as a consequencethe transmissionfunction (being essen-
tially the product of the two curves) is suppressed in a
We have presented a first-principles method for mod-
broad energy window around E . Thus chemical pas-
F ellingquantumtransportinmolecularnanostructuresbe-
sivation of the sulfur is the main reason for the lower
yond the single-particle approximation. The method is
conductance observedin structures (b,c) as comparedto
basedonnon-equilibriumGreenfunctions andappliesto
(a). A secondary effect, giving rise to differences in the
the general case of a finite bias voltage, but in this work
transmission of structure (b) and (c) is the different in-
we focused on the zero bias regime. The conductance of
terface dipoles which shift the electrostatic potential at
benzenedithiolate(BDT),benzenedithiol(BDT+H),and
the C H moiety by different amounts. This shift, how-
6 4 benzenediamine (BDA) wascalculatedusing the selfcon-
ever, has little influence on the conductance due to the
sistentGWapproximation. IncontrasttostandardDFT
flatness of the transmission functions around E .
F and Hartree-Fock methods, the GW approximation was
To verify this scenario, we have applied a scissors op- found to yield consistently accurate values for the en-
erator of 1 eV to the C H moiety of structure (b) in ergy levels of both isolated and contacted molecules due
6 4
order to align the onsite energy of the HOMO-1 to that to its proper treatment of self-interaction and dynami-
ofstructure(c). Theresultingtransmissionfunction(not cal screening. In general, results obtained with GW for
shown)isverysimilartothatof(c)andtheconductance the electronic conductance and energy gaps of contacted
is 3.5 10−2G compared to 3.3 10−2G obtained for molecules lie in between the values obtained with DFT
0 0
· ·
structure (c). On the other hand, the conductance of andHartree-Fock. The latter methods respectivelyover-
BDTcannotbebroughtbelow0.1G bysiftingthelevels estimates and underestimates the screening and none of
0
ofthe C H moiety down(by up to 2 eV) because ofthe themareabletodetectthechangeinthemolecule’selec-
6 4