Range-dependent adiabatic connections

Range-dependent adiabatic connections

A. M. Teale,^1,a兲 S. Coriani,^1,2and T. Helgaker¹

1Department of Chemistry, Centre for Theoretical and Computational Chemistry, University of Oslo, P.O. Box 1033, Blindern, N-0315 Oslo, Norway

2Dipartimento di Scienze Chimiche, Università degli Studi di Trieste, Via Licio Giorgieri 1, I-34127 Trieste, Italy

共Received 11 May 2010; accepted 17 August 2010; published online 27 October 2010兲

Recently, we have implemented a scheme for the calculation of the adiabatic connection linking the Kohn–Sham system to the physical, interacting system. This scheme uses a generalized Lieb functional, in which the electronic interaction strength is varied in a simple linear fashion, keeping the potential or the density fixed in the process. In the present work, we generalize this scheme further to accommodate arbitrary two-electron operators, allowing the calculation of adiabatic connections following alternative paths as outlined by Yang关J. Chem. Phys. 109, 10107 共1998兲兴.

Specifically, we examine the error-function and Gaussian-attenuated error-function adiabatic connections. It is shown that while the error-function connection displays some promising features, making it amenable to the possible development of new exchange-correlation functionals by modeling the adiabatic connection integrand, the Gaussian-attenuated error-function connection is less promising. We explore the high-density and strong static correlation regimes for two-electron systems. Implications of this work for the utility of range-separated schemes are discussed.

©2010 American Institute of Physics.关doi:10.1063/1.3488100兴

I. INTRODUCTION

The adiabatic-connection 共AC兲 formula for the exchange-correlation energy^1–5 in density-functional theory 共DFT兲 has motivated the construction of orbital-dependent functionals,^6,7 which represent some of the most successful approximations in widespread use. The AC formula arises from a consideration of the link between the Kohn–Sham noninteracting system and the physical, interacting system as a function of the interaction strength. A number of studies have examined the AC using approximate methods,^8–16 and some high-accuracy studies have been carried out for few- electron atomic systems.^17–22 Recently, we presented an implementation of a scheme allowing the calculation of ac- curate AC curves fromab initio densities^22,23 via optimiza- tion of Lieb functionals.²⁴Our implementation considers not only the usual density-fixed AC, of relevance in DFT, but also the potential-fixed AC,²² of relevance in potential- functional theory共PFT兲,²⁵in which the fully interacting sys- tem is related to the noninteracting, bare-nucleus system 共with the potential fixed at the external potential from the nuclei兲. The same connection was considered independently by Gross and Proetto,²⁶ who also discussed the differences between potential-functional theory variants based on the bare-nucleus noninteracting system 共as examined here兲 and potential-functional theories based on the Kohn–Sham non- interacting system as put forward by Yang and co-workers.²⁵ The relationship between the density- and potential-fixed ACs is particularly clear from the point of view of the Lieb formulation of DFT and will be further elucidated in the present work.

Most previous studies of the AC consider only the case in which the electron-electron repulsion is modulated in a simple linear fashion, by introducing a straightforward scal- ing of the two-electron interaction. However, as was pointed out by Yang,²⁷ this choice is not unique. In fact, the elec- tronic interaction may be modified by any function that smoothly connects the noninteracting and physical systems.

These generalized ACs are of particular relevance to theories constructed to combine the Kohn–Sham DFT and wave- function approaches as proposed by Savin.²⁸With an appro- priate modification of the electronic interaction, it is possible to attempt the construction of hybrid theories, in which short-range interactions are treated by DFT and long-range interactions by a suitable choice of wave-function methodol- ogy. Recently, a variety of short-range DFT functionals have been developed within the local-density approximation,^28,29 the generalized gradient approximation共GGA兲,^30–32 and the meta-GGA.³³ Several implementations of these hybrid schemes exist combining short-range DFT functionals with long-range Hartree–Fock 共HF兲,³⁴ configuration-interaction 共CI兲,^35,36 second-order Møller–Plesset,³⁴ coupled-cluster,³¹ multiconfigurational self-consistent field,^37,38 andn-electron valence second-order perturbation-theory³⁹methods. We also note that range separation of only the exchange interaction has been explored in the context of developing new DFT exchange-correlation functionals. Notable examples are the long-range corrected 共LC兲 functionals developed by Hirao and co-workers,^40,41the␻-PBEh and HSE functionals devel- oped by Scuseria and co-workers,^42,43and the CAM-B3LYP functional developed by Handy and co-workers.^44,45 These functionals emphasize either short-range^42,43 or long-range^40,41,44 interactions; a variant emphasizing the middle range has also been reported.⁴⁶ Finally, in a different

a兲Electronic mail: [email protected]. FAX:⫹47 228 55441.

0021-9606/2010/133共16兲/164112/14/$30.00 133, 164112-1 © 2010 American Institute of Physics

context, we note the use of a family of similar interactions by Gill and co-workers^47–49to remove the long-range tails of the Coulomb interaction.

While some studies of nonlinear ACs have been carried out from the point of view of calculating short-range DFT exchange-correlation energies and potentials,^21,50 no direct studies of the generalized, range-dependent AC integrand have been presented, in contrast to the linear case. Given the central role that the AC formulation plays in the theory un- derlying range-separated approaches, we are motivated to consider the generalization of our previously introduced scheme to this task. These generalized ACs are also of a wider interest than range-separated methods. This point was clear in the work of Yang,²⁷who considered the complemen- tary error function for modulation of the electronic interac- tion. By choosing an alternative form for the AC, the shape of the integrand is altered. For the complementary error func- tion, this means that both the noninteracting and physical points are known to be simple constants. This behavior sharply contrasts the corresponding linear case, where the noninteracting point is the orbital-exchange energy func- tional of the Hartree–Fock theory and the interacting point an expectation value of the full CI共FCI兲wave function.

In the present work, we generalize our optimization scheme for the Lieb functionals²² to electronic interactions weighted by the error function and the Gaussian-attenuated error function. We commence, in Sec. II, by introducing the theory of ACs with general two-electron operators and then briefly review our approach to optimization of the Lieb func- tionals and calculation of the AC integrands. Here we focus on details specific to this generalized scheme, referring the reader to our previous paper²² and the work of Wu and Yang²³for details of the optimization scheme. In Sec. III, we present results for the calculation of ACs corresponding to FCI densities for some simple two-electron systems, present- ing potential-fixed as well as density-fixed connections. We also discuss the prospects for approximating these alternative connections by simple forms suitable for a self-consistent implementation. Finally in Sec. IV, we make some conclud- ing remarks and discuss directions for future work.

II. THEORY

A. Lieb’s convex conjugate theory

Consider anN-electron system described by the Hamil- tonian

Hˆ

␭关v兴=Tˆ+Wˆ

␭+

兺

_i ^{v共ri兲, 0ⱕ ␭ ⱕ1, 共1兲}

where v共r兲 is the external potential at r, Tˆ is the kinetic- energy operator

Tˆ= −1

2

兺

_i ^{ⵜi2, 共2兲}

andWˆ

␭is a generalized electron interaction operator depend- ing on a coupling-strength parameter ␭ that varies between

␭= 0共the noninteracting system兲and␭= 1共the fully interact- ing system兲,

Wˆ

␭= 1

2

兺

_i⫽j^{w␭共rij兲, w0共rij兲= 0, w1共rij兲= 1/rij. 共3兲}

We now introduce the ground-state energy E_␭关v兴 as a func- tional of the external potential and the energy F_␭关␳兴 as a functional of the electron density by the following con- strained minimizations^24,51–53 over density matrices␥^ˆ:

E_␭关v兴= inf

␥ˆ→N

TrHˆ

␭关v兴␥^{ˆ= TrHˆ}_␭关v兴␥^ˆ_␭^v, 共4兲

F_␭关␳兴= inf

␥ˆ→␳TrHˆ

␭关0兴␥^ˆ= TrHˆ

␭关0兴␥^ˆ_␭^␳, 共5兲

where we denote the minimizers by␥^ˆ_␭^{v and}␥^ˆ_␭^␳, respectively.

Whereas a minimizer␥^ˆ_␭^␳ always exists in Eq.共5兲, this is not so for the minimization in Eq.共4兲, where␥^ˆ_␭^v only exists for those potentialsvthat support an electronic ground state for a given interaction strength ␭. In the following, we shall always assume that a minimizer exists.

As first discussed by Lieb,²⁴the ground-state energy as a functional of the external potentialE_␭关v兴and the energy as a functional of the density F_␭关␳兴 are conjugate functionals 共mutual Legendre–Fenchel transforms兲,

E_␭关v兴= inf

␳苸X共F_␭关␳兴+共v兩␳兲兲, 共6兲

F_␭关␳兴= sup

v苸X^ⴱ

共E_␭关v兴−共v兩␳兲兲, 共7兲

where the domains X and X^ⴱ are reflexive Banach spaces such that 共v兩␳兲=兰v共r兲␳共r兲dr is finite for all ␳苸X and v 苸X^ⴱ. In general, we obtain from Eqs.共6兲and共7兲the Fenchel inequality

E_␭关v兴ⱕF_␭关␳兴+共v兩␳兲, 共8兲 which holds for allvand␳. In the absence of degeneracies, the conditions for a minimizing density␳^{in Eq.}共6兲and for a maximizing potentialvin Eq.共7兲are equivalent and may be expressed in the following manner:

E_␭关v兴=F_␭关␳兴+共v兩␳兲⇔␦^E␭关v兴

␦v共r兲 =␳共r兲⇔␦^F␭关␳兴

␦␳共r兲 = −v共r兲, 共9兲 where it is assumed that兰␦␳共r兲dr= 0. An external potentialv and a density ␳ that together satisfy Eq. 共9兲 are said to be conjugate. For a given potential v, one or more conjugate densities ␳ may be found provided the potential supports a 共possibly degenerate兲 N-electron ground state. Conversely, an N-electron density␳ has a conjugate potentialv 共unique to within an additive constant兲provided␳^isv-representable.

Substituting Eqs.共4兲and共5兲in Eq.共9兲, we note the relation

␥^ˆ_␭^v=␥^ˆ_␭^␳ 共v and ␳ conjugate at ␭兲, 共10兲 which is valid for conjugate v and ␳ in the absence of de- generacies. In the present work, we consider ACs in which we fix the共nondegenerate兲density at its physical value共the AC of DFT兲and an alternative connection, in which we fix the potential共the AC of PFT兲at the nuclear-attraction poten- tial since this potential corresponds to the ␭= 1 system. In

other words, we consider connections that have conjugatev and␳^at␭= 1.

B. The adiabatic connection

Let us now relate the functionals E_␭关v兴 and F_␭关␳兴 for

␭⬎0 to the corresponding noninteracting quantities E₀关v兴 andF₀关␳兴, respectively,

E_␭关v兴=E₀关v兴+

冕

0

␭

E_␯⬘关v兴d␯^, 共11兲

F_␭关␳兴=F₀关␳兴+

冕

0

␭

F_␯⬘关␳兴d␯^, 共12兲

where the prime denotes differentiation with respect to␯^{. On} the right-hand side of these equations, we insert the expres- sions for the noninteracting energies E₀关v兴 and F₀关␳兴 ob- tained by setting ␭= 0 in Eqs.共4兲 and共5兲. Next, we deter- mine the derivatives E_␯⬘关v兴 and F_␯⬘关␳兴 by differentiation of Eqs. 共4兲 and共5兲 followed by application of the Hellmann–

Feynman theorem, leading to the following AC expressions:

E_␭关v兴=H_s关v兴+

冕

0

␭

W_␯关v兴d␯^, 共13兲

F_␭关␳兴=T_s关␳兴+

冕

0

␭W_␯关␳兴d␯^. 共14兲

We have here introduced the noninteracting bare-nucleus and kinetic-energy functionals

H_s关v兴= inf

␥ˆ→N

TrHˆ

0关v兴␥^ˆ= TrHˆ

0关v兴␥^ˆ0v

, 共15兲

T_s关␳兴= min

␥ˆ→␳TrHˆ

0关0兴␥^{ˆ= TrHˆ}0关0兴␥^ˆ0␳, 共16兲

and the potential- and density-fixed AC integrands as expec- tation values of the differentiated two-electron operatorWˆ

␯⬘^, W_␯关v兴= TrWˆ

␯⬘␥^ˆ_␯^v, 共17兲

W_␯关␳兴= TrWˆ

␯⬘␥^ˆ_␯^␳, 共18兲

with respect to the density matrices ␥^ˆ_␯^{v and}␥^ˆ_␯^␳ optimized at interaction strength␯^{from Eqs.}共4兲and共5兲, respectively. The perturbative expansion of Eqs.共11兲and共12兲in␭leads to the bare-nucleus⁵⁴ and Görling–Levy^55,56 perturbation theories, respectively, as discussed in Ref.22.

Let us now consider the relationship between the potential- and density-fixed connections. From Fenchel’s in- equality Eq.共8兲applied at␭= 0, we obtain

H_s关v兴ⱕT_s关␳兴+共v兩␳兲 共v and ␳ arbitrary兲, 共19兲 where equality occurs when ␳ is conjugate to v at ␭= 0.

Substituting Eqs.共13兲and共14兲into the stationary condition Eq.共9兲and invoking Eq.共19兲, we obtain the inequality

冕

0

␭W_␯关v兴d␯ⱖ

冕

0

␭W_␯关␳兴d␯ 共v and ␳ conjugate at ␭兲,

共20兲 which is valid providedvand␳ are conjugate at␭. Finally, introducing Eq. 共10兲 in Eqs. 共17兲 and 共18兲, we note that W_␯关v兴=W_␯关␳兴 when v and ␳ are conjugate for interaction strength ␯, in the absence of degeneracies. In the present work, the potential-fixed connection hasvequal to the physi- cal external potential due to the nuclei for all interaction strengths␯, making it relevant to PFTs based on this external potential. An alternative PFT was recently discussed by Yang and co-workers,²⁵in which the energy is expressed as a func- tional of the Kohn–Sham potential. In the present context we note that construction of an AC for such a theory would mean using the maximizing potential of Eq.共7兲at all␯^{. This} potential is different at each value of␯ but always conjugate to the physical density. As a consequence this alternative potential-based AC and the density-fixed AC become identi- cal.

From the concavity ofE_␭关v兴andF_␭关␳兴 in␭, it follows that these functions can always be represented in the form of Eqs. 共13兲 and共14兲, where the integrands W_␯关v兴andW_␯关␳兴 are nonincreasing right-continuous functions in␯. Under the assumption of adiabaticity, the two integrands become equal to the derivatives E_␯⬘关v兴 and F_␯⬘关␳兴 in Eqs. 共11兲 and 共12兲, respectively.

C. Coulomb, exchange, and correlation energies It is customary to decompose the total interaction ener- gies in Eqs.共13兲and共14兲in the manner

冕

0

␭

W_␯关v兴d␯^=J␭关v兴+E_x,␭关v兴+E_c,␭关v兴, 共21兲

冕

0

␭

W_␯关␳兴d␯^=J_␭关␳兴+E_x,␭关␳兴+E_c,␭关␳兴, 共22兲

where we have introduced the classical Coulomb functionals 共␳0v is the density associated with␥^ˆ0v兲

J_␭关v兴=1

2

冕冕

^{w␭共r12兲␳0v共r1兲␳0v共r2兲dr1dr2, 共23兲}

J_␭关␳兴=1

2

冕冕

^{w␭共r12兲␳共r1兲␳共r2兲dr1dr2, 共24兲}

the exchange functionals E_x,␭关v兴= TrWˆ

␭␥^ˆ0v−J_␭关v兴, 共25兲

E_x,␭关␳兴= TrWˆ

␭␥^ˆ0␳−J_␭关␳兴, 共26兲

and the correlation functionals

E_c,␭关v兴=

冕

0

␭Wc,␯关v兴d␯^, Wc,␯关v兴= TrWˆ

␯⬘共␥^ˆ_␯^v−␥^ˆ0v兲, 共27兲 E_c,␭关␳兴=

冕

0

␭

W_c,␯关␳兴d␯^, W_c,␯关␳兴= TrWˆ

␯⬘共␥^ˆ_␯^␳−␥^ˆ0␳兲. 共28兲 To show Eqs.共21兲and共22兲, we substitute Eqs.共23兲–共28兲in these equations and use the relation 兰₀^␭Wˆ

␯⬘^d␯^=Wˆ_␭^{. The ex-} change and correlation energies may be combined to give the exchange-correlation energies, which by combination of Eqs.

共25兲and共26兲with Eqs.共27兲and共28兲are given by E_xc,␭关v兴=

冕

0

␭Wxc,␯关v兴d␯^, Wxc,␯关v兴= TrWˆ

␯⬘␥^ˆ_␯^v−J_␯⬘关v兴, 共29兲 E_xc,␭关␳兴=

冕

0

␭

W_xc,␯关␳兴d␯^, W_xc,␯关␳兴= TrWˆ

␯⬘␥^ˆ_␯^␳−J_␯⬘关␳兴.

共30兲 In the following, we shall study the potential- and density- fixed AC connections and their contributions for the helium isoelectronic series and H₂ at different internuclear separa- tions, with different choices ofWˆ

␭.

D. One- and two-electron contributions

The Hamiltonian in Eq.共1兲provides a natural decompo- sition of the total electronic energy into one- and two- electron contributions,

E_␭关v兴= TrHˆ

0关v兴␥^ˆ_␭^{v+ TrWˆ}_␭␥^ˆ_␭^v, 共31兲 F_␭关␳兴= TrHˆ

0关0兴␥^ˆ_␭^␳+ TrWˆ

␭␥^ˆ_␭^␳, 共32兲

which may be further decomposed into uncorrelated and cor- related parts. The uncorrelated energy is obtained by the sub- stitution␥^ˆ_␭^v→␥^ˆ₀^vin Eq. 共31兲and the substitution ␥^ˆ_␭^␳→␥^ˆ₀^␳ in Eq.共32兲,

E_u,␭关v兴= TrHˆ

0关v兴␥^ˆ0v+ TrWˆ

␭␥^ˆ0v=H_s关v兴+J_␭关v兴+E_x,␭关v兴, 共33兲 F_u,␭关␳兴= TrHˆ

0关0兴␥^ˆ0␳+ TrWˆ

␭␥^ˆ0␳=T_s关␳兴+J_␭关␳兴+E_x,␭关␳兴.

共34兲 The uncorrelated one-electron energies are thus simply the noninteracting energies H_s关v兴 and T_s关␳兴 in Eqs. 共15兲 and 共16兲, respectively, whereas the uncorrelated two-electron en- ergies are the Coulomb and exchange energies evaluated from␥^ˆ0v and␥^ˆ0␳, respectively. The correlation energy is next obtained by the substitutions ␥^ˆ_␭^v→␥^ˆ_␭^v−␥^ˆ0v in Eq. 共31兲 and

␥^ˆ_␭^␳→␥^ˆ_␭^␳−␥^ˆ₀^␳ in Eq.共32兲, yielding

E_c,␭关v兴= TrHˆ

0关v兴共␥^ˆ_␭^v−␥^ˆ0

v兲+ TrWˆ

␭共␥^ˆ_␭^v−␥^ˆ0 v兲

=H_c,␭关v兴+E_c,^2el_␭关v兴, 共35兲 E_c,␭关␳兴= TrHˆ

0关0兴共␥^ˆ_␭^␳−␥^ˆ0␳兲+ TrWˆ

␭共␥^ˆ_␭^␳−␥^ˆ0␳兲

=T_c,␭关␳兴+E_c,^2el_␭关␳兴, 共36兲 where we use the conventional notationE_c,␭关␳兴=F_c,␭关␳兴. For the standard connection, w_␭共rij兲=␭/rij and it follows that Wˆ

1⬘^=Wˆ1. Comparing the integrands in Eqs. 共27兲 and 共28兲 with the two-electron parts in Eqs. 共35兲 and 共36兲, we then find that the two-electron correlation energy is equal to the AC integrand at␭= 1 :E_c,1^2el关v兴=Wc,1关v兴andE_c,1^2el关␳兴=Wc,1关␳兴. However, these relations are not valid for all possible con- nectionsw_␭.

Finally, we note that the one-electron contributions共ki- netic and interaction with the external potential兲 to the total energy must be the same for conjugatevand␳^,

H_s关v兴+H_c,␭关v兴=T_s关␳兴+T_c,␭关␳兴+共v兩␳兲

共v and ␳ conjugate at ␭兲. 共37兲 Combining this result with Fenchel’s inequality for the non- interacting system Eq.共19兲, we obtain the inequality

H_c,␭关v兴ⱖT_c,␭关␳兴ⱖ0 共v and ␳ conjugate at ␭兲, 共38兲 where the non-negativity of T_c,␭关␳兴follows from the defini- tion ofT_s,␭关␳兴as the lowest kinetic-energy expectation value consistent with the density␳^.

E. Range-independent and range-dependent connections

Thus far, we have established the potential-fixed inte- grandsW_␭关v兴,Wxc,␭关v兴, andWc,␭关v兴and density-fixed inte- grands W_␭关␳兴, Wxc,␭关␳兴, and Wc,␭关␳兴, whose coupling- constant integration yields the total interaction energy, the exchange-correlation energy, and the correlation energy, re- spectively. Yang²⁷observed that since these integrals are de- termined entirely by the functional values at the end points of the integration共␭2⬎␭1兲,

冕

_␭^␭₁^{2W␭关v兴d␭=E␭2关v兴−E␭1关v兴, 共39兲}

冕

_␭^␭₁^{2W␭关␳兴d␭=F␭2关␳兴−F␭1关␳兴, 共40兲}

we may choose Wˆ

␭ freely in Eq. 共1兲 provided its end-point values共typically 0 and 1兲are unaffected. This idea provides a justification for the proposal of Savin²⁸ to construct a va- riety of hybrid theories that merge wave-function approaches with DFT from the viewpoint of a generalized AC.²⁷ For further discussion see Sec. III E.

While some studies have appeared examining the inte- grated quantities E_xc关␳兴 and E_c关␳兴 for the density-fixed connection,^21,50no explicit study of the integrands involved, varying the path between the noninteracting and interacting

systems, has been carried out. In the present work, we con- sider the following general forms forw_␭共rij兲in Eq.共3兲:

w_␭^s共r_ij兲= ␭ rij

共standard兲, 共41兲

w_␭^e共rij兲=

erf

冉

1 −^␭␭rij

冊

r_ij 共error function兲, 共42兲

w_␭^g共rij兲=

erf

冉

1 −^␭␭rij

冊

rij

− 2

冑

_␲

冉

1 −^␭␭

冊

⫻exp

冉

⁻¹3

冉

^{1 −␭}␭

冊

^2rij2

冊

共Gaussian-attenuated error function兲, 共43兲 whose␭ derivatives are given by

w_␭^s⬘共r_ij兲= 1

r_ij, 共44兲

w_␭^e⬘共rij兲=

2 exp

冉

⁻

冉

1 −^␭␭

冊

^2rij2

冊

冑

_␲共1 −␭兲² , 共45兲 w_␭^g⬘共r_ij兲=

2 exp

冉

⁻

冉

^{1 −␭}␭

冊

^2rij2

冊

冑

_␲共1 −␭兲² +

2

冉

²³

冉

^{1 −␭}␭

冊

^{2rij2− 1}

冊

^exp

冉

^{− 13}

冉

^{1 −␭}␭

冊

^2rij2

冊

冑

_␲共1 −␭兲² . 共46兲 The choice w_␭^s in Eq. 共41兲 represents the standard range- independent AC, depending linearly on␭. As␭increases, the interaction is turned on uniformly for all interelectronic sepa- rationsrij. By contrast, with the error-function connectionw_␭^e in Eq. 共42兲 and Gaussian-attenuated error-function connec- tion w_␭^g in Eq.共43兲, the interaction is turned on in a range- dependent, nonuniform manner by the use of the functions erf共␮rij兲, and exp共−␮²r_ij²/3兲ofrij, where␮=␭/共1 −␭兲varies over the range 0ⱕ␮ⱕ⬁when␭increases from 0 to 1. As a result, with these two connections, long-range interactions are accounted for first and short-range interactions last. To illustrate the difference between the above connections, we have in Fig. 1 plotted the functions in Eqs. 共41兲–共43兲 and their derivatives in Eqs.共44兲–共46兲as functions ofrij, for four different values of␭. The␭dependence of the derivatives is relevant since in the evaluation of the AC integrands in Eqs.

共17兲and共18兲, we calculate the expectation value of the den- sity matrix with these derivatives.

To evaluate the AC integrands corresponding to the dif- ferent choices ofWˆ

␭in Eqs.共41兲–共43兲, we must calculate the expectation values ofWˆ

␭⬘ with a wave function correspond- ing to a fixed potential or a fixed density, determined by the

optimizations of Eqs. 共6兲 and共7兲. The minimization of Eq.

共6兲requires only standard techniques with two-electron inte- grals modified as described, for example, in Ref. 57. The maximization of Eq.共7兲is more difficult but can be achieved quite efficiently by the method in Refs. 22and23. Expecta- tion values of the derivatives in Eqs.共44兲–共46兲with the op- timized wave functions necessary for the calculation of the generalized AC integrands require the evaluation of two- electron integrals of the types 共ab兩exp共−␥^rij2兲兩cd兲 and 共ab兩rij2 exp共−␥^rij2兲兩cd兲, where the exponent ␥ is determined by ␭. Such integrals occur in R12 theories and as such are available in a variety of codes; in the present work, we use the integrals implemented for R12 theories by Samson et al.,⁵⁸ specifically the I₂ and I₄ integrals of that paper. The procedure is then to choose a suitable wave function for ac- curate determination of the Lieb functional, as described in Refs. 22 and 23 using modified two-electron integrals ac- cording to the choice of two-electron interaction from Eqs.

共41兲–共43兲. Once optimized, the expectation values of the de- rivatives required for the calculation of the AC integrands in Eqs. 共44兲–共46兲 are calculated, the linear AC being particu- larly simple since only the standard expectation value of the usual two-electron operator is required, see Eq.共44兲.

III. RESULTS

A. Computational details

The Lieb maximization of Eq. 共7兲 has been performed using the algorithm proposed by Wu and Yang²³共see Ref.18 for an alternative approach兲, which has been implemented recently²² in the DALTON quantum chemistry program⁵⁹ for arbitrary interaction strengths and the generalized ACs dis- cussed in Sec. II E. The reader is referred to Refs.22and23 for details of the implementation; here we note that the key to the approach of Wu and Yang²³to perform the Lieb maxi- mization is the parametrization of the potential in the man- ner,

vb共r兲=v_ext共r兲+v_ref,␭共r兲+

兺

t

btgt共r兲, 共47兲

where the first termv_ext共r兲is the external potential due to the nuclei, the second termv_ref,␭共r兲is a fixed reference potential chosen to ensure the correct asymptotic behavior, and the final term is a linear expansion in Gaussian functions gt共r兲 with coefficientsbt. Inserting this expansion into Eq.共7兲and using the gradient and Hessian with respect to the coeffi- cients bt, the Lieb maximization may be performed using

FIG. 1. Attenuated operators共top row兲and their␭derivatives共bottom row兲 as functions ofr₁₂for␭= 0共pink line兲,␭= 1/4共blue line兲,␭= 1/2共green line兲,␭= 3/4共red line兲, and␭= 1共black line兲.

standard quasi-Newton or Newton techniques. Here we have used the Newton method employing both the gradient and Hessian with a truncated singular-value decomposition cutoff of 10⁻⁶ and a convergence target of less than 10⁻⁶ on the gradient norm; for further details see Refs.22and23. All of the energiesE_␭关v兴for the two-electron systems in the present work are calculated at the FCI level.

B. The choice of basis sets

In order to perform the Lieb maximization we therefore must choose both a primary orbital basis set and an auxiliary potential basis set. When the potential basis set is chosen to be very different to that of the orbital basis set, unphysical oscillatory potentials can be obtained. This problem has been widely discussed in the literature in the context of the opti- mized effective potential method^60–71 and more recently in the context of constrained-search procedures⁷² at ␭= 0. To illustrate these effects we examine the exchange-correlation potentials for the helium atom along the density-fixed range- independent AC in Fig. 2. The uncontracted aug-cc-pVXZ basis sets have been employed for both the orbital and po- tential expansions. The potentials plotted in Fig.2represent the combinations X=兵6 , 2其, X=兵4 , 4其, and X=兵2 , 6其, where the first number refers to the orbital-basis cardinal number, X_orb, and the second to the potential-basis cardinal number, X_pot. For each combination, the potentials for interaction strengths␭ from 0.0 to 1.0 are shown in increments of 0.1.

The combination X=兵6 , 2其 may be regarded as unbal- anced in the sense that the potential basis is much smaller than that of the orbital-basis; while the potentials for this combination are smooth, the lack of flexibility in the poten- tial expansion may limit the variational freedom of the cal- culation. The combination X=兵4 , 4其 represents a balanced choice; here a small peak can been seen close to the nucleus, although the potential is predominantly smooth. It can be removed by application of the smoothing norm procedure of Ref.72as was done in our previous work.²²The final com- binationX=兵2 , 6其is unbalanced in the sense that the poten- tial basis set is much larger than the orbital basis set; for this combination the unphysical feature at the nucleus becomes much larger. Further increasing the size of the potential basis set can cause these oscillations to grow further.

In TableIwe explore the impact of different choices for the auxiliary potential basis set on the expectation value W_␭关␳兴, which is central to the calculation of the density- fixed AC in Eq.共14兲. Results are presented for ␭= 0.0, 0.5, and 1.0. Each row represents a choice of orbital basis, with each column corresponding to a different potential basis. For each value of␭and all choices of orbital basis, it is apparent that the expectation value is remarkably stable with respect to variations in the auxiliary basis. Furthermore, these small variations are largest for␭= 0 where the potential is largest 共see Fig.2兲and reduce steadily to zero for␭= 1.0, where the potential is zero.

The expectation value has much more significant varia- tion with respect to the orbital basis cardinal number, as would be expected. For each value of ␭, we present W_␭关␳兴

for orbital basis sets with 2ⱕX_orbⱕ6, along with an estimate of the basis-set-limit value calculated using the two-point extrapolation formula

EXY=X³EX−Y³EY

X³−Y³ 共48兲

in Ref.73withX= 5 andY= 6. These extrapolated results are denoted 关56兴; for a discussion of the application of this for- mula in the context of AC calculations, see Ref. 22. In the final row for each value of␭we have presented the deviation of the X_orb= 4 values from the estimated orbital basis set limit, denoted by ⌬. For all interaction strengths between 0 and 1 the absolute deviation from this limit is less than 1 mhartree. In light of this analysis we choose the uncontracted aug-cc-p共C兲VQZ basis sets for both the orbital and potential

FIG. 2. Exchange-correlation potentials for the helium atom along the density-fixed AC with␭= 0.0– 1.0 in steps of 0.1. Pane共a兲corresponds to the basis set combinationX=兵6 , 2其for the orbital basis and potential basis cardinal numbers, respectively. Pane 共b兲corresponds to the combination X=兵4 , 4其and pane共c兲to the combinationX=兵6 , 2其.

expansions for all of the two-electron systems in this study and quote all energetic values to a precision of 1 mhartree.

This choice of basis represents a good compromise between computational efficiency, accuracy, and adequate representa- tion of the exchange-correlation potential. Finally, we note that for the potential-fixed AC when calculating the expecta- tion value W_␭关v兴, we fix the potential vat the physical ex- ternal potential共due to the nuclei兲for all values of the inter- action strength. As such, only the orbital basis set plays a role in these calculations and similar accuracy is achieved in the uncontracted aug-cc-p共C兲VQZ basis sets.

C. The helium isoelectronic series

The helium isoelectronic series has been extensively studied and poses a significant challenge for approximate exchange-correlation functionals, particularly as Z increases.17,18,74–78

In the present work, we examine the sys- tems with 1ⱕZⱕ10 using the uncontracted aug-cc-pCVQZ basis set,^79–82noting that uncontraction and the use of core- correlating functions are essential to describe the compact densities accurately. The total energy and its components are listed in Table II, for the density-fixed connection共columns 3–9兲and the potential-fixed connection共columns 10–12兲. In Fig. 3, we have plotted the total AC integrandsW_␭关v兴 and W_␭关␳兴 for the three connections w_␭^s, w_␭^e, and w_␭^g in Eqs.

共41兲–共43兲, respectively.

As Z increases in the isoelectronic series, the density becomes more compact and may, to a good approximation, be expressed by a scaling of the density in H⁻:␳Z共r兲

⬇Z³␳H⁻共Zr兲. Consequently, the energy and its components increase in magnitude with increasingZ, in an approximately linear manner—see TableII. The only exception to this be- havior are the correlation energies, which remain approxi- mately constant withZ. These observations are in agreement with well-known scaling relations, such as J_␭关␳Z兴=ZJ_␭关␳H⁻兴 for the classical Coulomb energy.

Concerning the quality of the one-electron basis set, we note that the calculated bare-nucleus energyH_s关v兴in TableII differs from the exact bare-nucleus energyH_s关v兴=Z²by less than 0.001 E_h. Moreover, the virial theorem is satisfied to better than 1% for H⁻and better than 0.1% for the remaining systems.

Comparing the density- and potential-fixed results in Table II, we first note thatH_s关v兴is lower thanT_s关␳兴+共v兩␳兲 by 13% forZ= 1, by 3% forZ= 2, and by 0.1% forZ= 10, in agreement with Eq. 共19兲. Likewise, the positive quantity H_c关v兴 is several times larger than T_c关␳兴 for all Z, in agree- ment with Eq. 共38兲. Finally, comparing the positive quanti- ties 兰W_␭关v兴d␭ⱖ兰W_␭关␳兴d␭ 关see Eq. 共20兲兴, we find that the former is larger than the latter by 27% forZ= 1, by 10% for Z= 2, and by 2% forZ= 10. As expected, the energy changes that occur in the potential-fixed system with increasing␭are

TABLE I. The variation of the expectation valueW␭关␳兴with choice of orbital and auxiliary potential expansion basis sets for the density fixed AC of the helium atom. The uncontracted aug-cc-pVXZ basis sets have been used for both expansions. Each row represents the change in the expectation value with for a given orbital basis cardinal numberX_orbas the cardinal number of the potential basis,X_pot, is changed. For the definition of the quantities关56兴and⌬see text.

Exchange-correlation potentials corresponding to the values marked in bold are shown in Fig.2. All values in atomic units.

Orbital basisX_orb

Potential basisX_pot

2 3 4 5 6

␭=0.0

2 1.017 678 1.017 700 1.017 709 1.017 707 1.017 707

3 1.023 222 1.023 224 1.023 232 1.023 229 1.023 230

4 1.024 084 1.024 084 1.024 084 1.024 086 1.024 085

5 1.024 359 1.024 359 1.024 359 1.024 359 1.024 359

6 1.024 472 1.024 472 1.024 472 1.024 472 1.024 472

关56兴 1.024 628 1.024 628 1.024 628 1.024 628 1.024 628

⌬=W_␭⁴−W_␭^关56兴 ⫺0.000 544 ⫺0.000 544 ⫺0.000 543 ⫺0.000 542 ⫺0.000 543

␭=0.5

2 0.981 488 0.981 504 0.981 511 0.981 509 0.981 510

3 0.982 082 0.982 084 0.982 086 0.982 085 0.982 086

4 0.981 591 0.981 585 0.981 591 0.981 592 0.981 591

5 0.981 393 0.981 393 0.981 393 0.981 393 0.981 393

6 0.981 293 0.981 293 0.981 293 0.981 293 0.981 293

关56兴 0.981 155 0.981 155 0.981 155 0.981 155 0.981 155

⌬=W_␭⁴−W_␭^关56兴 0.000 436 0.000 430 0.000 436 0.000 437 0.000 436

␭=1.0

2 0.949 929 0.949 929 0.949 929 0.949 929 0.949 929

3 0.947 658 0.947 658 0.947 658 0.947 658 0.947 658

4 0.946 580 0.946 580 0.946 580 0.946 580 0.946 580

5 0.946 225 0.946 225 0.946 225 0.946 225 0.946 225

6 0.946 065 0.946 065 0.946 065 0.946 065 0.946 065

关56兴 0.945 844 0.945 844 0.945 844 0.945 844 0.945 844

⌬=W_␭⁴−W_␭^关56兴 0.000 735 0.000 735 0.000 735 0.000 735 0.000 735

larger in the corresponding density-fixed system. Physically, these relations may be understood from the observation that the density of the potential-fixed noninteracting system is more compact than that of the density-fixed system. We also note that the differences between the potential- and density- fixed quantities are largest forZ= 1, which has the most dif- fuse electron density. This behavior may be understood from the observation that in anionic systems, the diffuseness of the electron density arises from electron repulsion, which is ne- glected in the noninteracting limit, generating a too compact density in the potential-fixed connection.

A further comparison of the density- and potential-fixed ACs is given in Fig. 3, where we have plotted W_␭关v兴 and W_␭关␳兴 against␭ for 1ⱕZⱕ10, where w_␭^s is in pane共a兲,w_␭^e in pane 共b兲, andw_␭^g in pane共c兲. The density- and potential- fixed integrands are very similar and we note that for all values of ␭ and all three connections, W_␭关v兴ⱖW_␭关␳兴, as may be rationalized by observing that the density becomes more compact in the potential-fixed AC as the electronic interactions are turned off.

Comparing the AC curves arising from the different choices ofw_␭in Fig.3, we recall that each curve represents the expectation value of Wˆ

␭⬘ ^with ␥^ˆ_␭^{␳ or} ␥^ˆ_␭^v, optimized with the two-electron operatorWˆ

␭. The standard connectionw_␭^s in pane共a兲yields nearly straight lines, with a larger slope in the potential-fixed case 共dashed lines兲 than in the density-fixed case 共full lines兲, representing a situation where the interac- tions are turned on uniformly for all interelectronic separa- tions. Thew_␭^e curves in pane共b兲give the same total interac- tions as those in pane共a兲but have very different shapes since the interactions are now first turned on for large interelec- tronic separations and subsequently for short separations.

The AC curves are therefore no longer linear but contain a peak at that value of ␭ where most of the interactions are recovered. ForZ= 1, the peak is broad and occurs already at

␭⬇0.1, reflecting the large range of interelectronic separa- tions that contribute to the interactions in this diffuse system.

For the most compact system with Z= 10, there is a sharp peak at ␭⬇0.87, indicating that most interactions occur at about 0.2a₀– 0.3a₀. The w_␭^g plots in pane 共c兲 are similar to those in pane共b兲but have shaper peaks, reflecting the higher locality ofw_␭^g⬘, see Fig.1.

For the density-fixed AC, we consider separately also the

TABLE II. Energy components of the helium isoelectronic series in the uncontracted aug-cc-pCVQZ basis共atomic units兲.

Z E_tot共Z兲 T_s关␳兴 共v兩␳兲 兰W␭关␳兴d␭ J关␳兴 E_x关␳兴 E_c关␳兴 T_c关␳兴 H_s关v兴 兰W␭关v兴d␭ H_c关v兴

1 ⫺0.527 0.502 ⫺1.374 0.345 0.773 ⫺0.386 ⫺0.041 0.028 ⫺1.000 0.473 0.155

2 ⫺2.903 2.865 ⫺6.751 0.983 2.048 ⫺1.024 ⫺0.041 0.036 ⫺4.000 1.097 0.150

3 ⫺7.279 7.238 ⫺16.126 1.609 3.302 ⫺1.651 ⫺0.042 0.039 ⫺9.000 1.721 0.152

4 ⫺13.654 13.613 ⫺29.501 2.234 4.553 ⫺2.277 ⫺0.043 0.041 ⫺16.000 2.346 0.153

5 ⫺22.030 21.987 ⫺46.875 2.858 5.804 ⫺2.902 ⫺0.043 0.041 ⫺25.000 2.970 0.153

6 ⫺32.405 32.361 ⫺68.249 3.483 7.054 ⫺3.527 ⫺0.044 0.042 ⫺36.000 3.595 0.154

7 ⫺44.780 44.735 ⫺93.623 4.108 8.304 ⫺4.152 ⫺0.044 0.042 ⫺49.000 4.220 0.154

8 ⫺59.155 59.109 ⫺122.998 4.733 9.554 ⫺4.777 ⫺0.044 0.043 ⫺64.000 4.845 0.154

9 ⫺75.530 75.484 ⫺156.372 5.359 10.805 ⫺5.402 ⫺0.044 0.043 ⫺81.000 5.470 0.154

10 ⫺93.905 93.859 ⫺193.747 5.984 12.055 ⫺6.027 ⫺0.044 0.043 ⫺100.000 6.095 0.154

FIG. 3. AC curves共atomic units兲 W␭关␳兴 共full lines兲 andW␭关v兴 共dashed lines兲for the helium isoelectronic series with 1ⱕZⱕ10 forw_␭^s in pane共a兲, forw_␭^ein pane共b兲, and forw_␭^gin pane共c兲. In all panes, the curves increase with increasingZ.

exchange-correlation and correlation contributions to the full AC curve: W_␭关␳兴=W_xc,␭关␳兴+W_c,␭关␳兴. In Fig. 4, we have plotted W_␭^s关␳兴, W_␭^e关␳兴, and W_␭^g关␳兴 and their exchange- correlation and correlation contributions for the helium iso- electronic series. The W_␭^s关␳兴 curves in pane共a兲are positive and nearly constant since the dominant Coulomb and ex- change energies increase linearly with ␭:W_␭^s关␳兴⬇J1关␳兴 +E_x,1关␳兴=J₁关␳兴/2. The exchange-correlation curves in pane 共b兲 are dominated by the exchange energy W_xc,␭^s 关␳兴

⬇E_x,1关␳兴and are approximate mirror images of the curves in pane共a兲.

Pane共c兲 in Fig.4 shows the correlation-only integrand W_c,␭^s 关␳兴for the helium isoelectronic series, on a much larger scale than that used in panes共a兲and共b兲. The curvature of the H⁻ curve is much more pronounced than for the other spe- cies; as Z increases, the density accumulates close to the nucleus and the curves become more linear. This behavior can be understood from the relation E_c关␳Z兴⬇Z²E_c,1/Z关␳H⁻兴, which follows from the observed scaling of␳Zwith increas- ing Z and a general scaling relation of the correlation

energy.⁸³ For largerZ, the AC effectively explores a smaller

␭ interval of some approximately universal AC curve valid for all Z. Consequently, these curves become more linear with increasing charge as the system approaches the high- density limit. This trend toward linearity and the rate at which it occurs are clear in Fig.4. We note that linearity of the correlation AC curve means that the correlation energy increases quadratically with␭, as expected from the validity of second-order Görling–Levy perturbation theory^55,56 for these systems.

In the range-dependent error-function curves in the sec- ond row of Fig.4, long-range interactions are recovered for small values of ␭, while short-range interactions are recov- ered for large␭values. For the total integrandW_␭^e关␳兴in pane 共d兲, the height of the peak increases and moves to the right with increasing Z, as the density contracts and the interac- tions become more short-ranged. As for the standard connec- tion w_␭^s, the exchange-correlation in pane 共e兲 is an approxi- mate mirror image of the total curve in pane 共d兲: W_xc,␭^e 关␳兴

⬇−W_␭^e关␳兴. The range separation induced by the error func-

FIG. 4. AC curves共atomic units兲for the helium isoelectronic series with 1ⱕZⱕ10 forw_␭^sin panes共a兲–共c兲, forw_␭^ein panes共d兲–共f兲, and forw_␭^gin panes共g兲–共i兲. For each AC, we have plotted the total curveW␭关␳兴to the left, the exchange-correlation curveWxc,␭关␳兴in the middle, and the correlation curveWc,␭关␳兴to the right. In panes共a兲,共d兲, and共g兲, the curves increase with increasingZ; in the other panes, the curves may be distinguished by noting that the same color scheme is used in all panes.