Correlation potentials for molecular bond dissociation within the self-consistent random phase approximation

0.24 MB·English

by  M. Hellgren

Correlation potentials for molecular bond dissociation within the self-consistent
random phase approximation
Maria Hellgren1, Daniel R. Rohr1,2, and E. K. U. Gross1
1Max-Planck-Institut fu¨r Mikrostrukturphysik, Weinberg 2, 06120 Halle (Saale), Germany
2Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 2-4, 14195 Berlin, Germany
(Dated: January 24, 2012)
Self-consistent correlation potentials for H2 and LiH for various inter-atomic separations are ob-
tained within the random phase approximation (RPA) of density functional theory. The RPA
2 correlation potential shows a peak at thebond midpoint, which is an exact feature of thetruecor-
1 relationpotential,butlacksanotherexactfeature: thestepimportanttopreserveintegerchargeon
0
the atomic fragments in the dissociation limit. An analysis of the RPA energy functional in terms
2
of fractional charge is given which confirms these observations. We find that the RPA misses the

n derivativediscontinuityatoddintegerparticlenumbersbutexplicitlyeliminatesthefractional spin
a error in the exact-exchangefunctional. The latter findingexplains the accurate total energy in the
J dissociation limit.

3
2
I. INTRODUCTION ments. Thiserrorisknownintheliteratureasdelocaliza-
] tion error or fractional charge error.22,23,29 The delocal-
h
izationerroroftheRPAhasbeenstudiedonlyforthedis-
The random phase approximation (RPA) in Kohn-
-p Sham (KS) density functional theory (DFT) has in re- sociationofopen-shellH+2 andHe+2.30 Itwasfoundtobe
m cent years received considerable attention in quantum rather severe leading to too low total energies. Whether
e chemistry1–3 andmaterialscience.4–11 TheRPAincorpo- the RPA suffers from fractionalcharge errorin the cases
h rates exchangeeffects exactly andthe correlationenergy of heteronuclear molecules is presently unknown.
c istreatednon-perturbativelybysummingasubsetofpo- The exactfunctionalensuresneutraldissociationfrag-
s. larization diagrams to infinite order.12,13 Furthermore, mentsbyvirtueoftheKSpotential.31–33 Inthedissocia-
c the RPA can be systematically improved being the first tionlimitthehighestoccupiedorbitalofeachofthefrag-
i
s approximation within the so-called adiabatic connection ments must be aligned (or degenerate). Consequently,
y fluctuation dissipation (ACFD) framework.14–16 the highest occupied molecular orbital (HOMO) of the
h
The RPA is an implicit functional of the density and wholesystem(includingbothfragments)isalinearcom-
p
[ canthereforeincludenon-localcorrelationeffectslikee.g. binationoftheorbitalsofeachfragment. Toobtainequal
the van der Waals interactions. That these are, indeed, weights, i.e. integer charge,at the fragments the KS po-
2 accurately captured by the RPA has been demonstrated tentialexhibits a sharpstep atthe bond midpoint, shift-
v
in many recent works.17–19 For systems described by ing the energy levels of only one of the fragments. This
2
strong Hubbard-like correlations, the RPA is, however, feature is a directconsequenceofthe derivativedisconti-
6
0 stillnotfullyinvestigated. Apopulartestcaseinthisre- nuity in the correlation part of the energy.34,35 Another
6 gardis the dissociation of diatomic molecules with cova- feature of the exact KS correlation potential of dissoci-
. lent bonds.20 All density functional approximationscon- ated molecules is a peak at the bond midpoint.36 The
0
1 structedsofarfailinthiscontextifproperspin-symmetry peak emerges with increasing inter atomic distance, and
1 is enforced. The total energy in the dissociation limit is acts to further localize the electrons.
1 toohighandspuriousfractionalchargesarefoundatthe The RPA correlation potential for atoms was recently
: fragments. obtained37,38 and showed a close resemblance to the ex-
v
i The largeerrorinthe totalenergyhasbeencharacter- act correlation potential. However, all RPA calculations
X izedasstaticcorrelationerrororso-calledfractionalspin onmoleculeshavesofarbeencarriedoutusingpotentials
r error,studiedindetailinthepioneeringworksofCohen, originating from other functionals and hence precluding
a
Mori-S`anchez and Yang.21–23 It has been demonstrated a full assessment of the RPA. The aim of this work is
thattheRPAstronglyimprovesthedissociationlimitfor to provide a more complete analysis of the RPA. To this
homoatomic systems such as the H molecule.24–27 purpose, we have calculated self-consistent RPA poten-
2
Spuriousfractionalcharges,ontheotherhand,appear tials for molecules and investigated the RPA correlation
only in the dissociation limit of heteroatomic molecules potentials for H2 and LiH at different inter-atomic dis-
and are related to an incorrect behavior of the total en- tances. These systems allow us to study both the static
ergy as a function of particle number. The exact en- correlationerror and the delocalization error. Moreover,
ergy functional exhibits a kink or derivative disconti- the LiHis anexamplewhere aself-consistentcalculation
nuity at integer particle numbers along with a straight isessential. WehavealsoanalyzedtheRPAenergyfunc-
line behavior between the integers.28 The smooth and tionalin terms offractionalchargeby studying the RPA
non-linear behavior of approximate functionals leads to functional on an extended domain of spin-compensated
chargedfragmentsinthedissociationlimit,meaningthat densities allowing for non-integer number of particles.
oneelectronistoodelocalized,i.e.,spreadoverbothfrag- We concludethatthe RPApotentials exhibitthepeak
2
at the bond midpoint but lack the step feature, where Here, we have used the notation (r ,t ) = 1 etc. and
1 1
the latter is related to a missing intra-shell derivative introduced Λ(3,2;1) = iG (3,1)G (1,2). The correla-
s s
−
discontinuity. The total energy of LiH is, however, still tion part of the self-energy Σ in the RPA is given by
c
largely improved in the dissociation limit as compared
δE
to,e.g.,theexact-exchange(EXX)functional,suggesting Σ =i c =ivχRPAvG (6)
c s
only a smaller delocalization error. δGs
where
II. RPA CORRELATION ENERGY AND χRPA =χs+χsvχRPA. (7)
POTENTIAL
In the appendix we show the expression for n , i.e., the
c
right hand side of Eq. (5), in terms of KS orbitals and
WithintheACFDframeworktheexactcorrelationen-
KS eigenvalues.
ergy is written as39–41
i 1 ∞ dω III. FRACTIONAL CHARGE AND SPIN
E = dλ Tr v[χ (ω) χ (ω)] (1)
c λ s
2Z Z 2π { − }
0 −∞
TheRPAfunctionalproducesaccuratedissociationen-
whereχsisthenon-interactingKSdensityresponsefunc- ergiesforH2,incontrasttoallcommondensityfunction-
tion and χλ is the scaled interacting density response als which yield a far too high energy due to a spurious
function. The scaling refers to a system with a linearly self-interactionintheHfragments. Itiswell-knownthat
scaled Coulomb interaction λv(r,r′) plus a fictitious po- EXXisself-interactionfreeinthecaseofaspin-polarized
tentialwhichkeepsthedensityfixedasλischanged. The Hatom. InthedissociationlimitofH ,theHatomsare,
2
parameter λ runs between 0 (non-interacting KS case) however, not spin-polarized, but rather described by a
and 1 (fully interacting case). We have used the short mixture of a spin-up and a spin-down H atom, in which
handnotationTrfg = drdr′f(r,r′)g(r′,r)foranytwo- caseEXXdoessufferfromself-interaction. Toobtainthe
pointfunctionsf andgR. WithinTDDFTthefunctionχλ correct dissociation limit a significant correlationcontri-
reads42 butionis thusneeded. Inthefollowingwewillshowthat
theRPAcorrelationfunctionalexactlycancelsthespuri-
χλ =χs+χs λv+fxλc χλ. (2) ous self-interactionin the EXX functional in the dissoci-
(cid:2) (cid:3) ationlimit. Thetotalenergywillstillnotbeexactasthe
The scaled XC kernel fλ is defined as the functional correlationenergycontainsaself-correlationterm,which
xc
derivative of the scaled XC potential vλ with respect to does not vanish in the dissociation limit.
xc
the density n, evaluated at the ground state density.
In the RPA fλ = 0 and thus corresponds to the
xc
simplest approximation within the ACFD formalism. A. RPA in the dissociation limit
Within the RPA the λ-integral in Eq. (1) can be evalu-
ated analytically with the result The RPA energy is an explicit functional of the KS
GreenfunctionG andforspin-compensatedsystemsG
s s
i dω is proportional to the identity matrix in spin-space. In
E = Tr ln[1 vχ ]+vχ . (3)
c −2Z 2π { − s s} the frequency domain a spin-component of Gs reads
occ uocc
DiagrammaticallyEq. (3)is equalto aninfinite summa-
G (r,r′,ω)= G−(r,r′,ω)+ G+(r,r′,ω) (8)
tion of ring-diagrams. s k k
Xk Xk
The RPA correlation potential v can be obtained as
c
the functional derivative of Eq. (3) with respect to the where we have defined
density. IfweletVs signifythetotalKSpotential,Gs the ϕ (r)ϕ (r′)
non-interacting KS Green function and χ = iG G , G±(r,r′,ω)= k k (9)
s − s s k ω ε iη
the functional derivative can be obtained via the chain − k±
rule with ϕ being a KS spin-orbital and ε the correspond-
k k
ing eigenvalue. Consider now a stretched homonuclear
δE δn δE δG
n c = c s. (4) diatomic moleculecomposedofatomAandatomB.For
c
≡ δn δV δG δV
s large but finite inter-atomic separation R the molecular
KS orbitals can approximately be written as symmetric
The result is the well-known linearized Sham-Schlu¨ter
and antisymmetric linear combination of the atomic KS
(LSS) equation43,44
orbitals
1
Z χs(1,2)vc(2)d2=Z Λ(3,2;1)Σc(2,3)d2d3. (5) ϕk(r)= √2[ϕAk(r)±ϕBk(r)] (10)
3
where ϕA(r) is a KS orbital of atom A. This expression where 0 referstothegroundstateofN evennumberof
k | i 0
becomes exact as R . It is easy to show that in the electrons, referestoaspin-compensatedN +2-state
0
→∞ |↑↓i
dissociation limit and ( )aspin-down(up)polarizedstatewithN +1
0
|↓i |↑i
electrons. The definition of the non-interacting spin-up
ERPA[Gs]→ERPA[GAs]+ERPA[GBs ] (11) G↑α and spin-down G↓α Green functions are
G↑(rt,r′t′)= αT ψ↑(rt)ψ↑(r′t′)† α (19)
where α h | { }| i
G↓(rt,r′t′)= αT ψ↓(rt)ψ↓(r′t′)† α (20)
occ uocc α h | { }| i
1
GAs = Xk6=0G−k +kX6=0G+k + 2(cid:2)G−0 +G+0(cid:3). (12) winhtehree KαScasynstbeema.nHyeoref,tTheissttahteestim|0ei,-o|rd↓ie,ri|n↑giooprer|a↑to↓ri
andψβ(rt)thefieldoperatorwiththepropertyofadding
Here, GAs contains only states of the isolated atom A. ψβ(rt)† or removing ψβ(rt) an electron with spin β in
The orbital k = 0 is the special orbital of the highest r at time t. Evaluating the Green functions as ensem-
occupiedstatewhichhastobeconsideredpartiallyoccu- ble expectation values using Eqs. (17-18) we find Green
piedandpartiallyunoccupied. Forahomoatomicsystem functions exactly of the form as in Eq. (15-16). Thus,
the fraction is always 1/2. the ensembles proposed correspond to how an atom in
Inthecaseofacovalentbondedheteronucleardiatomic the RPA is described in the dissociation limit.
molecule a similar analysis can be carried out. The only TheRPAcorrelationenergyisusuallyevaluatedusing
differenceisthatnowonlythe highestoccupiedandlow- the time-ordered KS response function χ . Using the
s
est unoccupied arise from a degeneracy of the isolated ensemble Green functions we find
atoms. Due to the lack of symmetry we also have to al-
χE(r,r′,ω) = iπδ(ω)p(2 p)ϕ (r)2 ϕ (r′)2
low for a more generallinear combination of KS orbitals s − − | 0 | | 0 |
+χ˜p(r,r′,ω) (21)
s
1
ϕL0UMO(r) = √2[√pϕA0(r)− (2−p)ϕB0(r)] (13) where
p
ϕHOMO(r) = 1 [ (2 p)ϕA(r)+√pϕB(r)],(14) χ˜p(r,r′,ω) = 2 fq(r)fq(r′) fq(r)fq(r′) (22)
0 √2 p − 0 0 s Xq ω−εq+2iη − ω+εq−2iη
where p [0,2]. The KS Green functions to be inserted and we have performed a spin summation. The or-
∈
in Eq. (11) become bitals ϕ are now referring to the orbitals of an iso-
k
lated atom. The index q = (k,k′) is a composite in-
occ uocc p 2 p dex of occupied (k) and unoccupied (k′) states. The
GAs = G−k + G+k + 2G−0 + −2 G+0 (15) special transition q = (0,0) is excluded since it has
Xk6=0 Xk6=0 been taken out explicitly (first term in Eq. (21)). The
GBs = occ G−k +uoccG+k + 2−2 pG−0 + 2pG+0. (16) ftuionncstioannsdfwq(er)n=oteϕkt(hra)tϕkw′h(re)naϕre0ciaslloedcceuxpcieitdatiitonshfuounlcd-
Xk6=0 Xk6=0 be multiplied by p/2 and when it is unoccupied by
(2 p)/2. Whenpcalculatingthecorrelationenergythe
In the exact system the energy assumes its minimum at −
pinteracting RPA response function has to first be evalu-
p = 1. This is, however, not guaranteed using an ap-
ated χE = χE +λχEvχE. The final expression of the
proximate functional. In fact, most functionals yield a λ s s λ
RPA correlation energy reads
so-calledfractionalchargeerroronthe atomicfragments
(i.e. the minimum is found at p=1). p(2 p)
6 Ec = − drdr′ ϕ0(r)2v(r,r′)ϕ0(r′)2
− 4 Z | | | |
i 1 ∞ dω
B. Fractional charge + dλ Tr v[χ˜p(ω) χ˜p(ω)] (23)
2Z Z 2π { λ − s }
0 −∞
In this section we derive an expression for the RPA where χ˜p(ω) is the scaled interacting response function
λ
correlation energy in terms of fractionally charged spin- calculated using Eq. (22). This equation will allow us
compensatedsystems. Tothisend,weproposeensembles to investigate the fractional spin error in the following
of the following form sectionand to obtain numericalresults for the fractional
charge error in Section V.
p
γˆ< = (1 p)0 0 + ( + ) (17)
− | ih | 2 |↑ih↑| |↓ih↓|
for p [0,1) C. Fractional spin error
∈
2 p
γˆ> = − ( + )+(p 1) (18)
2 |↑ih↑| |↓ih↓| − |↑↓ih↑↓| Thefractionalspinerrorisdefinedastheenergydiffer-
for p (1,2] ence of a systemwith N +1 electrons with one electron
0
∈
4
spinpolarized(sp)ononehandandthesamesystembut As a reference potential we use the sum of the nuclear
fully spin compensated (sc) on the other hand. Let us potential and the Fermi-Amaldi potential, vFA(r), eval-
splitthe RPAinteractionenergyinto the sumofHartree uated with the Hartree-Fock density
and exchange energy (Hx) and the correlation energy. N 1 n(r′)
The sum of Hx energy in the sp case reads vFA(r)= − dr′. (28)
N Z r r′
| − |
N0/2 N0/2 For closed shell systems the derivative of the total en-
EHx = ij ij + 0i 0i (24) ergyfunctionalwithrespecttotheexpansioncoefficients,
sp h || i h || i
Xi,j Xi bt, is readily evaluated via Eq. (4). For systems with
fractional charge we evaluated the gradient with the fi-
where we introduced the anti- nite difference method. The gradient is then fed to the
symmetrized integral ij kl = Broyden-Fletcher-Goldfarb-Shannooptimization routine
drdr′ϕ∗i(r)ϕ∗j(r′)v(r,r′)(2ϕk(r)ϕl(hr′)||−iϕl(r)ϕk(r′)). to find the minimum totalenergy. The calculationswere
RThe second term accounts for all interactions of the consideredconvergedwhenthe vectornormofthe gradi-
singly occupied atomic orbital ϕ0. The expression for ent was less than 10−3. The algorithmwas implemented
the spin compensated case reads in a local version of the PyQuante module.46
It is a well-known fact that finite basis set OEP cal-
1
EHx =EHx+ drdr′ ϕ (r)2v(r,r′)ϕ (r′)2. (25) culations can become numerically unstable if a too large
sc sp 4Z | 0 | | 0 | potentialbasissetisused.47–49Wecarefullychosethepo-
tential basis sets to be balancedto the respective orbital
The additional term is a nonzero self-interaction term,
basis sets. The proper choice is reflected in the smooth
which is equal to the fractional spin error of EXX. It is
potentials that we obtain. For our calculations we used
this term alone that is responsible for the wrong disso-
the standard orbital basis cc-pVTZ from the Dunning
ciation limit of EXX for H . It also introduces a large
2 family.
additional error in the dissociation limit of any system,
For the potential bases we used even tempered gaus-
on top of the errorinherent in the EXX functional when
sians. Each set is characterized by three numbers: the
describing the fragments.
smallestexponent,α ,thelargestexponent,α ,and
min max
We now turn to the RPA correlation energy. The cor-
the number of basis functions, N. With this set of pa-
relation energy for the sc system is taken from Eq. (23)
rameters we construct the exponent α as
i
with p=1.
(i−1)
α N−1
max
Escc = −41Z drdr′|ϕ0(r)|2v(r,r′)|ϕ0(r′)|2 αi =αmin(cid:20)αmin(cid:21) , (29)
whereirunsfrom1toN. Fortheatomiccalculationswe
i 1 ∞ dω
+ dλ Tr v[χ˜1(ω) χ˜1(ω)] (26) use only s-type functions. For molecular calculations we
2Z0 Z−∞ 2π { λ − s } add a set of p-type functions using the same exponents
and omitting the largest.
We now see that the first term exactly cancels the frac-
For the one dimensional (1D) calculations we use a
tionalspinerrorofEXX.The secondterm inEq. (23)is
basis set of cubic splines, which permits us to solve
identical to the correlation energy obtained in the spin-
Eq. (5) by inverting the static KS response function.
polarized case. Consequently, the RPA functional does
The spline-basisis described indetailin severalprevious
not suffer from fractional spin error. There will, how-
works where it has been used successfully for OEP-type
ever, still be a self-correlation due to the second term of equations.37,50,51
but this is, in general, expected to be smaller. These
conclusions are confirmed by the numerical results ob-
tained previously24–26 and here in Sec. V. V. RESULTS
As a first step we verify our implementation. To this
IV. COMPUTATIONAL DETAILS end, we compare the total energy and the correlation
potential with accurate reference data for He, Be, and
Our goal is to find the local potential that mini- Ne.38 The total energies and ionization potentials (IP)
mizes the RPA energy functional. For the three dimen- arelistedintableI. Thepotentialbasissetsaregivenby
sional (3D) calculations we utilize the direct minimiza- αmin =0.01, αmax =10, and N =5 for He, αmin =0.01,
tion scheme for the optimized effective potential (OEP) αmax =10,andN =9forNe,andαmin =0.001,αmax =
proposed by Yang and Wu.45 The potential is expanded 100, and N =11 for Ne (compare equation (29)).
in a basis plus a reference potential The total energies of our new implementation are
somewhat higher than the reference energies. The dif-
v(r)=v (r)+ b g (r). (27) ferences range from 9 mH to 158 mH. This result is ex-
0 t t
pected since the reference calculations were performed
Xt
5
this work benchmark
Atom IP energy IP energy RPA correlation potential for Be
He -0.907 -2.936 -0.902 -2.945 RPA
0.0
Be -0.349 -14.681 -0.354 -14.754 Benchmark
Ne -0.773 -128.945 -0.796 -129.103 Exact
0.1
(cid:1)
TABLE I: Energies in Hartree.
0.2
(cid:1)
RPA correlation potential for He
0.04
0.3
(cid:1)
0.02
0.00 (cid:1)0.4
0.02 10-2 10-1 100 101
(cid:0)
Distance to nucleus in Bohr
0.04
(cid:0) FIG. 2: The correlation potential of Be for RPA38 (dotted),
0.06
(cid:0) for RPA in this work (dashed), and the exact correlation
(cid:0)0.08 RPA potential52 (solid).
Benchmark
0.10
(cid:0) Exact
0.12 RPA correlation potential for Ne
(cid:0) 10-2 10-1 100 101 0.15
Distance to nucleus in Bohr
FIG. 1: The correlation potential of He for RPA38 (dot- 0.10
ted),forRPAinthiswork(dashed)andtheexactcorrelation
potential52 (solid). 0.05
0.00
with a virtually complete orbital basis. The differences
0.05
in IP are minor (from 5 mH to 23 mH), thus indicating (cid:2)
RPA
an accurate potential. This conclusion is supported by
0.10 Benchmark
(cid:2)
figures 1 - 3, which show our RPA correlation potentials
Exact
(dashed) of He, Be, and Ne comparing to accurate RPA 0.15
correlation potentials (dotted).38 The RPA correlation (cid:2) 10-2 10-1 100 101
Distance to nucleus in Bohr
potentials shown in the figures are calculated as
FIG. 3: The correlation potential of Ne for RPA38 (dotted),
vRPA(r)=vRPA[n ](r) vEXX[n ](r). (30) for RPA in this work (dashed), and the exact correlation
c KS RPA − KS EXX potential52 (solid).
Inallfigureswequalitativelyreproducethereferencepo-
tentials. The largest deviation is found close to the nu-
cleus (the x-axis is on a logarithmic scale). For com- With our new implementation we are able to investi-
parison we also show the exact KS correlation potential gate the H molecule in all three dimensions and with
2
(solid)52 and we see that the RPA correlation potential thefullcoulombinteraction(potentialbasis: α =0.1,
min
closely resembles the exact correlation potential. α = 1.0, and N = 5). In Figure 5 we show the RPA
max
To analyze the RPA molecular correlation poten- correlation potential (see equation (30)) along the bond
tial we first investigate a 1D model system with axis,whereagainthebondmidpointisatzero. Weshow
soft coulomb interactions. The nuclear potentials are the potentials for equilibrium distance (1.4 Bohr; solid
Z/ (x R/2)2+1,whereZ isthenuclearchargeand curve),for5.0(dotted),andfor10.0Bohr(dashed). The
− ±
R ispthe internuclear distance. The electron-electron in- potentials qualitatively resemble the potentials obtained
teraction is set to 1/ (x x )2+1. The first system from the 1D calculations (c.f. 4). Only the peak at the
1 2
−
consists of two H atopms (Z = 1.2). Figure 4 shows the bond midpoint is absent for small (1.4) and large (10.0)
RPA correlation potentials along the bond axis at bond interatomic separation. For large separationthe absence
distances of2,4, and6 Bohr,withthe bond midpoint at is easily explained. The orbitaland potential basis func-
zero. We notice that the minimum ofthe correlationpo- tionsaresimplynotdiffuseenoughtoextendtothebond
tential is shifted away from the atom, but that the shift midpoint. This results in a vanishing correlation poten-
decreases with increasing interatomic distance. At the tial at the bond midpoint. The situation is different for
bond midpoint a positive peak emerges with increasing a bond distance of 1.4 Bohr. With a small atomic sepa-
bond distance. Both features are also observed for the rationthe orbitalandpotentialbasisextendto the bond
exact KS correlation potential.31,33 midpointasisevidentfromthenon-vanishingcorrelation
6
RPA correlation potential of H RPA correlation potential of LiH
0.4 2 0.8
1D v1DRPA 2 Bohr v1DRPA at 3 Bohr
0.3 1D vc1DRPA 4 Bohr 0.6 vc1DRPA at 6 Bohr
c c
0.2 1D v1DRPA 6 Bohr v1DRPA at 8 Bohr
c 0.4 c
0.1
0.2
0.0
0.0
0.1
(cid:3)
(cid:3)0.2 (cid:5)0.2
0.3 0.4
(cid:3) 0 2 4 6 8 10 (cid:5) 10 5 0 5 10
Distance to bond midpoint in Bohr (cid:5) (cid:5)Distance to bond midpoint in Bohr
FIG. 4: The correlation potential along the bond axis with FIG. 6: The correlation potential along the bond with the
the bond midpoint at zero. The system is 1D H2 with a soft bond midpoint at zero. The system is 1D LiH with a soft
coulomb potential. coulomb potential. The Li atom is located at -1, -1.5, and -4
Bohr, respectively. The H atom is located at 1, 1.5, and 4
Bohr, respectively.
RPA correlation potential of H
0.05 2
0.00
However,astep,asispresentintheexactKScorrelation
(cid:4)0.05 potential, is not observed.
0.10 Figure 7 shows the correlation potential (equa-
(cid:4)
0.15 tion(30))ofthethreedimensionalLiHforbonddistances
(cid:4)
3 (solid), 6 (dotted), and 8 Bohr (dashed). The parame-
0.20
(cid:4) vcRPA at 1.4 tersforthepotentialbasisoftheHatomareαmin =0.01,
(cid:4)0.25 vcRPA at 5.0 αmax = 1.0, and N = 2. For the Li atom they read
(cid:4)0.30 vcRPA at 10.0 αmin =0.001,αmax =10,andN =3. The same features
as in the 1D caseare also found in the results for the 3D
0.35
(cid:4) 0 2 4 6 8 10
Distance to bond midpoint in Bohr correlation potentials. In the region of the Li atom (-10
to0)the potentialqualitativelyresemblesthatofthe 1D
FIG. 5: The correlation potential of H2 along the bond axis. system in Figure 6. We see a well with a peak close to
The bond midpoint is at zero. We show theRPA correlation the nucleus. Also in the region of the H atom (0 to 10)
potential for interatomic distance 1.4 (solid), 5.0 (dotted), figures 6 and 7 show the same structure. A well emerges
and 10.0 (dashed).
with increasing bond distance. The difference between
1D and 3D is found only at the bond midpoint. Like in
the case of H a peak emerges only for the 1D system
2
potential. However,the potential basis functions arenot (figure 6). In contrast, the 3D system (figure 7) exhibits
compactenoughtoproduceapeak. Wehaveverifiedthis a peak at zero only for small and intermediate bond dis-
by using more compact potential basis functions. With tances (3 and 5 Bohr). The explanation for the absence
this set of basis functions, however, we find unphysical of the peak at the bond midpoint for 8 Bohr is the same
oscillationsforlargeratomicseparations. Wealsoplaced as in the case of H .
2
orbitalandpotentialfunctionsatvariouspointsalongthe
We further analyze the RPA functional in the context
bond axis. We observed a peak in the potential at each
of fractional charge. The total energy of H and Li (po-
ofthepoints,whichleadsustoconcludethatthesepeaks
tentialbases like inLiH) is plotted (figures 8 and9) as a
are artifacts of the basis rather than genuine features of
functionofthe numberofelectrons,whereweallowfrac-
the functional.
tional values, according to Eq. 23. Both atoms show
Atthispointwerestrainfromshowingthedissociation a smooth behavior, thus showing that the RPA (red)
curve and rather point the reader to a future publica- misses the kink at integer electron numbers. Comparing
tionwithadetaileddiscussionofthedissociationcurves. totheexactcurveswesee,however,alargeimprovement
We would only like to mention that the well-known comparedtotheEXX(blue)functionalregardingtheto-
”bump”24–27 is still present, but somewhat weaker. tal energy. The agreement is particularly well at integer
We nowturnto the LiH molecule. Infigure6 weshow number of electrons. This explains the gooddissociation
the RPA correlation potential (defined in equation (30)) energiesfor homoatomic systems found in the RPA. We
of 1D LiH (Z = 3.6,Z = 1.2). The solid, dotted and can also combine these two figures to analyze the RPA
Li H
dashed curves represent the correlation potentials for 2, energy of LiH in the dissociation limit. In order to do
3, and 8 Bohr interatomic distances, respectively. The this we add the energies for the H atom to the energies
build up of the peak at the bond midpoint is apparent. oftheLiatomsothatthetotalnumberofelectronssums
7
RPA correlation potential of LiH Fractional charge for Li
0.05 7.30
(cid:8)
EXX
0.00
RPA
(cid:6)0.05 e)(cid:8)7.35 exact
(cid:6)0.10 artre
(cid:6)0.15 y (h(cid:8)7.40
g
0.20 er
(cid:6)(cid:6)0.25 vvccRRPPAA aatt 36..00 BBoohhrr otal en(cid:8)7.45
(cid:6)0.30 vcRPA at 9.0 Bohr T
(cid:6)10 (cid:6)D5istance to bond0 midpoint in Boh5r 10 (cid:8)7.50
2.4 2.6 2.8 3.0 3.2 3.4 3.6
FIG.7: Thecorrelation potentialofLiHalongthebondaxis. Number of electrons
Thebondmidpointisatzero. TheLiatomislocatedat-1.5,
-2.5,and-4Bohr,respectively. TheHatom islocated at1.5, FIG. 9: The total energy of Li as a function of number of
2.5, and 4 Bohr, respectively. electrons for EXX (blue),RPA (red), and exact (black).
Fractional charge for Li + H
Fractional charge for H 7.6
0.20 (cid:9)
(cid:7) EXX
EXX
(cid:7)0.25 RPA 7.7 RexPaAct
gy (hartree)(cid:7)(cid:7)000...433005 exact ergy (hartree)(cid:9)(cid:9)7.8
al ener(cid:7)(cid:7)0.45 otal en(cid:9)7.9
ot T
T 0.50
(cid:7)
8.0
0.55 (cid:9)
(cid:7)
0.4 0.6 0.8 1.0 1.2 1.4 1.6
0.4 0.6 0.8 1.0 1.2 1.4 1.6 Number of electrons at the H atom
Number of electrons
FIG. 10: The total energy of dissociated LiH as a function
FIG. 8: The total energy of H as a function of number of
of number of electrons at the H atom for EXX (blue), RPA
electrons for EXX (blue), RPA(red), and exact (black).
(red),and exact (black).
uptofour. Figure10showsthetotalenergyasafunction VI. CONCLUSIONS
of the number of electrons at the H atom. The number
of electrons at the Li atom will then be four minus the
We have obtained self-consistent RPA correlation po-
x value. The exact functional (black) has a minimum
tentialsfordiatomicmoleculesandstudiedtheirbehavior
at 1.0, because it dissociates LiH into a neutral H atom
asafunctionofinteratomicdistances. Atlargedistances
and a neutral Li atom. In Figure 10 we see that EXX
the RPA potential correctly exhibits a peak at the bond
(blue)andRPA(red)donotdissociateLiHintotheneu-
midpoint but misses the step feature.
tral atoms. In both cases there is a surplus of electronic
We have also analyzed the RPA functional for frac-
charge at the H atom. However, in the case of RPA the
tional charges by evaluating the ensemble averaged KS
surplus (0.16 electrons) is much smaller than in the case
Green function using spin-compensated ensembles with
of EXX (0.4 electrons). The large improvement may be
non-integernumberofparticles. Thisprocedurecaneas-
relatedtothepeakthatispresentintheRPAcorrelation
ilybecarriedovertoanyfunctionalconstructedfromthe
potential.
Finally, in table II we collect the total energies in the
dissociationlimit of H2, Li2, and LiH for the exact func- exact EXX RPA
tional,EXX,andRPA.PleasenotethattheLiHenergies H2 -1.000 -0.713 -1.035
of EXX and RPA are taken as the minimum of the re- Li2 -14.948 -14.749 -14.961
spective curves in figure 10 and not the values at 1.0. LiH -7.974 -7.759 -8.007
This is to account for the fractional charge error present
in both functionals. TABLEII:TotalenergiesinHartreeinthedissociation limit.
8
KS Green function. rootoftheeigenvaluescorrespondstothetrueexcitation
The numerical results show that the kink at integer energiesZ andthe excitationfunctions aretransformed
q
number of electrons is missed in the RPA. As a conse- according to
quence,wefindspuriousfractionalchargesonthedissoci-
ated fragments. The chargesare, however,much smaller Fq(r)= fq′(r)Uq′q, (A4)
compared to other functionals. On the other hand, the Xq′
RPAwillmostlikelynotbeabletodescribetheso-called
field counteracting effect in the correlation potential of where U is the matrix diagonalizing V. The interacting
hydrogen chains which has been discussed to have its time-ordered response function in the RPA can then be
origin in the derivative discontinuity.53 written as
Insummary,wehavefoundthattheRPAaccomplishes
1
the following: χRPA(r,r′,ω) = F (r)F (r′)
q q
2Z
The dissociation limit of closed-shell molecules is Xq q
•
well reproduced in the RPA, due to the explicit 1 1
(A5)
elimination of the self-interaction term present in ×(cid:20)ω Z +iη − ω+Z iη(cid:21)
q q
− −
the EXX functional.
and the integral in Eq. (A1) becomes
RPA separates the chargesin bond dissociation by
•
virtue of a peak at the bond midpoint, an exact
feature of the true correlationpotential. Σc(r,r′,ω)= ϕk(r) dr1v(r,r1)Fq(r1)
Z
Xkq
RPA exhibits only a small fractional charge er-
• ror in the cases of closed-shell covalently bonded ϕ (r′) dr v(r′,r )F (r )
k 1 1 q 1
× Z
molecules.
1 1 n n
k k
These findings consolidate the high expectations on the − + . (A6)
×2Z (cid:20)ω ε Z +iη ω ε +Z iη(cid:21)
q k q k q
RPA that are currently prevalent. There is, however, − − − −
still room for improvement as the discontinuity at odd
Next, we evaluate
integer particle number is missing. We believe that im-
provementsontheRPAwithintheACFDframeworkwill dω
n = i Σ (ω)G (ω)G (ω). (A7)
help to overcome this shortcoming. c c s s
− Z 2π
Afterperformingstandardcontourintegrationswegetin
Appendix A: Derivative of the RPA energy totalsixterms,whichaftersymmetryconsiderationscan
be reduced to four. In summary, we find
In this Appendix we return to Eq. (5) and evaluate
nc(r) in terms of KS orbitals, ϕk, and KS eigenvalues, n (r) = 1 (...)ϕ∗(r)ϕ (r)
ε . The self-energy (Eq. (6)) involves an integration c Z s p
k XkspXq q
over the frequency of the following form
dr dr F (r )v(r ,r )ϕ (r )ϕ∗(r )
Σ (ω)=i dω′G (ω′+ω)vχRPA(ω′)v, (A1) × Z 1 2 q 1 1 2 k 2 p 2
c s
Z 2π
dr dr ϕ (r )ϕ∗(r )v(r ,r )F∗(r ). (A8)
× Z 1 2 s 1 k 1 1 2 q 2
wherewehavesuppressedthespace-coordinates. Toper-
form the integration we write the KS Green function in
wherethedotssignifiestheinsertionofthefollowingfour
its Lehmann representation
terms
G (r,r′,ω) = ϕ (r)ϕ (r′)
s k k (1 n )n n
k p s
Xk −(ε +Z −ε )(ε +Z ε ) (A9)
k q p k q s
n 1 n − −
k + − k ,(A2) 2(1 nk)np(1 ns)
×(cid:20)ω ε iη ω ε +iη(cid:21) + − − (A10)
− k− − k (εp Zq εk)(εp εs)
− − −
where nk is the occupation number. The response func- 2nknp(1 ns)
tion in the RPA is given by Eq. (7). This equation can + − (A11)
(ε +Z ε )(ε ε )
s q k p s
easily be rewritten as an eigenvalue problem in terms of − −
n (1 n )(1 n )
the matrix + k − s − p . (A12)
(ε +Z ε )(ε +Z ε )
s q k p q k
where f˜q(r)=V2q√q′ω=qfωqq(2rδ)q,q′w+ithhf˜qf|qv(|rf˜)q′bi,eing a KS ex(cAit3a)- The expression in Eq.−(A8) has, in−this work, been im-
tionfunction and ω a KSexcitationenergy. The square plemented both in the 3D and in the 1D case.
q
9
1 F. Furche,Phys. Rev.B 64, 195120 (2001). 28 J.P.Perdew,R.G.Parr,M.Levy,andJ.L.Balduz,Phys.
2 F. Furche,J. Chem. Phys. 129, 114105 (2008). Rev.Lett. 49, 1691 (1982).
3 H.Eshuis,J.Yarkony,andF.Furche,J.Chem.Phys.129, 29 P. Mori-S´anchez, A. J. Cohen, and W. Yang, Phys. Rev.
234114 (2010). Lett. 100, 146401 (2008).
4 L. Schimka, J. Harl, A. Stroppa, A. Grneis, M. Marsman, 30 P. Mori-S´anchez, A. Cohen, and W. Yang,
F. Mittendorfer, and G. Kresse, Nature Materials 9, 741 ArXiv:0903.4403v1 (2009).
(2010). 31 O.V.GritsenkoandE.J.Baerends,Phys.Rev.A54,1957
5 S. Leb`egue, J. Harl, T. Gould, J. G. A´ngy´an, G. Kresse, (1996).
and J. F. Dobson, Phys. Rev.Lett. 105, 196401 (2010). 32 P. R. T. Schipper, O. V. Gritsenko, and E. J. Baerends,
6 J. Harl and G. Kresse, Phys. Rev. Lett. 103, 056401 Phys. Rev.A 57, 1729 (1998).
(2009). 33 D. Tempel, T. J. Martinez, and N. T. Maitra,
7 X. Ren, A. Tkatchenko, P. Rinke, and M. Scheffler, Phys. J. Chem. Theor. Comp. 5, 770 (2009).
Rev.Lett. 106, 153003 (2011). 34 J. P. Perdew, in Density Functional Methods in Physics,
8 X. Ren, P. Rinke, and M. Scheffler, Phys. Rev. B 80, edited by R. M. Dreizler and J. da Providencia (Plenum,
045402 (2009). New York,1985).
9 R. W. Godby,M. Schlu¨ter, and L. J. Sham, Phys. Rev.B 35 E. Sagvolden and J. P. Perdew, Phys. Rev. A 77, 012517
36, 6497 (1987). (2008).
10 R. W. Godby, M. Schlu¨ter, and L. J. Sham, 36 M. A. Buijse, E. J. Baerends, and J. G. Snijders,
Phys. Rev.Lett. 56, 2415 (1986). Phys. Rev.A 40, 4190 (1989).
11 M.Gru¨ning,A.Marini,andA.Rubio,J.Chem.Phys.124, 37 M. Hellgren and U. von Barth, Phys. Rev. B 76, 075107
154108 (2006). (2007).
12 A.L.FetterandJ.D.Walecka,QuantumTheoryofMany- 38 M.HellgrenandU.vonBarth,J.Chem.Phys.132,044101
Particle Systems (DoverPubl., New York,2003), 2nd ed. (2010).
13 F. Aryasetiawan, T. Miyake, and K. Terakura, 39 D. C. Langreth and J. P. Perdew, Solid State Commun.
Phys. Rev.Lett. 88, 166401 (2002). 17, 1425 (1975).
14 M. Fuchsand X.Gonze, Phys.Rev. B 65, 235109 (2002). 40 O.GunnarssonandB.I.Lundqvist,Phys.Rev.B13,4274
15 W. Kohn, Y. Meir, and D. E. Makarov, Phys. Rev. Lett. (1976).
80, 4153 (1998). 41 M. Marques, C. A. Ullrich, F. Nogueira, A. Rubio,
16 J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. K. Burke, and E. K. U. Gross, eds., Time-Dependent
A´ngy´an,Phys. Rev.Lett. 102, 096404 (2009). Density Functional Theory (Springer, Berlin, Heidelberg,
17 J. F. Dobson, in Topics in Condensed Matter Physics, 2006).
edited by M. P. Das (Nova, New York, 1994), cond- 42 E. K. U. Gross and W. Kohn, Phys. Rev. Lett. 55, 2850
mat/0311371. (1985).
18 H.-V. Nguyen and G. Galli, J. Chem. Phys. 132, 044109 43 L. J. Sham and M. Schlu¨ter, Phys. Rev. Lett. 51, 1888
(2010). (1983).
19 S. Leb`egue, J. Harl, T. Gould, J. G. A´ngy´an, G. Kresse, 44 U. von Barth, N. Dahlen, R. van Leeuwen, and G. Ste-
and J. F. Dobson, Phys. Rev.Lett. 105, 196401 (2010). fanucci, Physical Review B 72, 235109 (2005).
20 N. Helbig, I. V. Tokatly, and A. Rubio, J. Chem. Phys. 45 W. Yangand Q. Wu,Phys. Rev.Lett. 89, 143002 (2002).
131, 224105 (2009). 46 URL pyquante.sourceforge.net.
21 A. J. Cohen, P. Mori-S´anchez, and W. Yang, 47 D. R. Rohr, O. V. Gritsenko, and E. J. Baerends,
J. Chem. Phys.129, 121104 (2008). Chem. Phys. Lett. 432, 336 (2006).
22 P. Mori-S´anchez, A. J. Cohen, and W. Yang, Phys. Rev. 48 D. R.Rohrand A.Savin, J. Mol. Struct.943, 90 (2010).
Lett. 102, 066403 (2009). 49 A. Heßelmann, A. W. G¨otz, F. D. Salla, and A. G¨orling,
23 A.J.Cohen,P.Mori-S´anchez,andW.Yang,Science321, J. Chem. Phys.127, 054102 (2007).
792 (2008). 50 H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, and
24 N. E. Dahlen, R. van Leeuwen, and U. von Barth, F. Martin, Rep.Prog. Phys. 64, 1815 (2001).
Int.J. QuantumChem. 101, 512 (2005). 51 M.HellgrenandU.vonBarth,J.Chem.Phys.131,044110
25 M. Fuchs, Y.-M. Niquet, X. Gonze, and K. Burke, (2009).
J. Chem. Phys.122, 094116 (2005). 52 C.J.UmrigarandX.Gonze,Phys.Rev.A50,3827(1994).
26 N. E. Dahlen, R. van Leeuwen, and U. von Barth, 53 S. J. A. van Gisbergen, P. R. T. Schipper, O. V. Grit-
Phys. Rev.A 73, 012511 (2006). senko, E.J. Baerends, J. G. Snijders, B. Champagne, and
27 A. Heßelmann and A. G¨orling, Phys. Rev. Lett. 106, B. Kirtman, Phys. Rev.Lett. 83, 694 (1999).
093001 (2011).