Range-separated density-functional theory with random phase approximation: detailed formalism and illustrative applications
aa r X i v : . [ phy s i c s . c h e m - ph ] J u l Range-separated density-functional theory with random phase approximation:detailed formalism and illustrative applications
Julien Toulouse , ∗ Wuming Zhu , † J´anos G. ´Angy´an , ‡ and Andreas Savin § Laboratoire de Chimie Th´eorique, UPMC Univ Paris 06 and CNRS, 75005 Paris, France CRM2, Institut Jean Barriol, Nancy University and CNRS, 54506 Vandoeuvre-l`es-Nancy, France (Dated: November 8, 2018)Using Green-function many-body theory, we present the details of a formally exact adiabatic-connectionfluctuation-dissipation density-functional theory based on range separation, which was sketched in Toulouse,Gerber, Jansen, Savin and ´Angy´an, Phys. Rev. Lett. , 096404 (2009). Range-separated density-functionaltheory approaches combining short-range density functional approximations with long-range random phase ap-proximations (RPA) are then obtained as well-identified approximations on the long-range Green-function self-energy. Range-separated RPA-type schemes with or without long-range Hartree-Fock exchange response kernelare assessed on rare-gas and alkaline-earth dimers, and compared to range-separated second-order perturbationtheory and range-separated coupled-cluster theory.
I. INTRODUCTION
Range-separated density-functional theory has emerged asa powerful approach for improving the accuracy of standardKohn-Sham (KS) density-functional theory [1, 2] applied withusual local or semi-local density-functional approximations,in particular for electronic systems with strong (static) orweak (van der Waals) correlation e ff ects. Based on a sepa-ration of the electron-electron interaction into long-range andshort-range components, it permits a rigorous combination ofa long-range explicit many-body approximation with a short-range density-functional approximation (see, e.g., Ref. 3 andreferences therein). Several many-body approximations havebeen considered for the long-range part: configuration in-teraction [4, 5], multi-configuration self-consistent-field the-ory [6–8], second-order perturbation theory [9–13], coupled-cluster theory [14–18], multi-reference second-order pertur-bation theory [19], and several variants of the random phaseapproximation (RPA) [20–24].In the context of the recent revived interest in RPA-type ap-proaches to the electron correlation problem in atomic, molec-ular and solid-state systems [25–48], several range-separatedapproaches using long-range RPA-type approximations haveindeed been proposed and show promising results, in particu-lar for describing weak intermolecular interactions. Toulouse et al. [20] have presented a range-separated RPA-type the-ory including the long-range Hartree-Fock exchange responsekernel. Janesko et al. [21–23] have proposed a simpler range-separated RPA scheme with no exchange kernel and in whichthe RPA correlation energy has been rescaled by an empir-ical coe ffi cient. Paier et al. [24] have added the so-calledsecond-order screened exchange to the latter scheme, whichappears to correct the self-interaction error. In all these cases,range separation tends to improve the corresponding full-range RPA-type approach, avoiding the inaccurate description ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] and slow basis-set convergence of short-range correlations inRPA.In Ref. 20, only the main lines of range-separated density-functional theory with long-range RPA were presented. Inthis work, we give now all the missing details of the the-ory. Using Green-function many-body theory, we constructa formally exact adiabatic-connection fluctuation-dissipationdensity-functional theory based on range separation, with-out the need of maintaining the one-particle density constant.Range-separated RPA-type schemes are then obtained as well-identified approximations on the long-range Green-functionself-energy. The range-separated RPA-type methods with orwithout long-range Hartree-Fock exchange response kernelare assessed on rare-gas and alkaline-earth dimers, and com-pared to range-separated second-order perturbation theory andrange-separated coupled-cluster theory. The most tedious de-tails of the theory are given in the appendices.
II. THEORYA. Range-separated density-functional theory
In range-separated density-functional theory (see, e.g.,Ref. 3), the exact ground-state energy of an N -electron systemis expressed as the following minimization over multidetermi-nant wave functions Ψ E = min Ψ n h Ψ | ˆ T + ˆ V ne + ˆ W lr ee | Ψ i + E srH xc [ n Ψ ] o , (1)where ˆ T is the kinetic energy operator, ˆ V ne isthe nuclei-electron interaction operator, ˆ W lr ee = (1 / ! d r d r w lr ee ( r )ˆ n ( r , r ) is a long-range (lr)electron-electron interaction written with w lr ee ( r ) = erf( µ r ) / r and the pair-density operator ˆ n ( r , r ), and E srH xc [ n ] isthe corresponding µ -dependent short-range (sr) Hartree-exchange-correlation (Hxc) density functional that Eq. (1)defines. The parameter µ in the error function controls therange of the separation. The minimizing wave function,denoted by Ψ lr , yields the exact density. Several approxi-mations [3, 7, 14, 18, 49–51] have been proposed for theshort-range exchange-correlation (xc) functional E sr xc [ n ], andan approximate scheme must be used for the long-range wavefunction part of the calculation.In a first step, the minimization in Eq. (1) is restricted tosingle-determinant wave functions Φ , leading to the range-separated hybrid (RSH) approximation [9] E RSH = min Φ n h Φ | ˆ T + ˆ V ne + ˆ W lr ee | Φ i + E srH xc [ n Φ ] o , (2)which does not include long-range correlation. The mini-mizing determinant Φ is given by the self-consistent Euler-Lagrange equation ˆ H | Φ i = E | Φ i , (3)where E is the Lagrange multiplier for the normalization con-straint and ˆ H is the RSH reference Hamiltonianˆ H = ˆ T + ˆ V ne + ˆ V lrH x , HF [ Φ ] + ˆ V srH xc [ n Φ ] , (4)which includes the Hartree-Fock (HF)-type long-rangeHartree-exchange (Hx) potential ˆ V lrH x , HF [ Φ ] and the short-range local Hxc potential ˆ V srH xc [ n ] = R d r v srH xc [ n ]( r )ˆ n ( r )written with v srH xc [ n ]( r ) = δ E srH xc [ n ] /δ n ( r ) and the den-sity operator ˆ n ( r ). As usual, ˆ V lrH x , HF is the sumof a local Hartree part ˆ V lrH = R d r v lrH ( r )ˆ n ( r ) with v lrH ( r ) = R d r w lr ee ( r ) h Φ | ˆ n ( r ) | Φ i , and a non-local ex-change part ˆ V lr x , HF = ! d x d x v lr x ( x , x )ˆ n ( x , x ) writtenwith v lr x ( x , x ) = − w lr ee ( r ) h Φ | ˆ n ( x , x ) | Φ i and the one-particle density-matrix operator ˆ n ( x , x ) expressed withspace-spin coordinates x = ( r , s ) and x = ( r , s ).The RSH scheme does not yield the exact energy and den-sity, even with the exact short-range functional E srH xc [ n ]. Nev-ertheless, the RSH approximation can be used as a referenceto express the exact energy as E = E RSH + E lr c , (5)defining the long-range correlation energy E lr c , for which wewill now give an adiabatic connection formula. We introducethe following energy expression with a formal coupling con-stant λ E λ = min Ψ n h Ψ | ˆ T + ˆ V ne + ˆ V lrH x , HF [ Φ ] + λ ˆ W lr | Ψ i + E srH xc [ n Ψ ] o , (6)where the minimization is done over multideterminant wavefunctions Ψ , ˆ W lr is the long-range Møller-Plesset-type fluctu-ation perturbation operatorˆ W lr = ˆ W lr ee − ˆ V lrH x , HF [ Φ ] , (7)and E srH xc is the previously-defined λ -independent short-rangeHxc functional. The minimizing wave function, denoted by Ψ lr λ , is given by the self-consistent Euler-Lagrange equationˆ H lr λ | Ψ lr λ i = E lr λ | Ψ lr λ i , (8) where E lr λ is the Lagrange multiplier for the normalization con-straint and ˆ H lr λ is the long-range interacting e ff ective Hamilto-nian along the adiabatic connectionˆ H lr λ = ˆ T + ˆ V ne + ˆ V lrH x , HF [ Φ ] + ˆ V srH xc [ n Ψ lr λ ] + λ ˆ W lr , = ˆ H + λ ˆ W lr + (cid:16) ˆ V srH xc [ n Ψ lr λ ] − ˆ V srH xc [ n Φ ] (cid:17) . (9)For λ =
1, Eq. (6) reduces to Eq. (1), and so the physicalenergy E = E λ = and density are recovered. For λ =
0, theminimizing wave function is the RSH determinant Ψ lr λ = = Φ and the Hamiltonian of Eq. (9) reduces to the RSH referenceHamiltonian, ˆ H lr λ = = ˆ H . Note that, because the density at λ = E λ with respectto λ , noting that E λ is stationary with respect to Ψ lr λ , and rein-tegrating between λ = λ = E = E λ = + Z d λ h Ψ lr λ | ˆ W lr | Ψ lr λ i , (10)with E λ = = h Φ | ˆ T + ˆ V ne + ˆ V lrH x , HF [ Φ ] | Φ i + E srH xc [ n Φ ] = E RSH −h Φ | ˆ W lr | Φ i . Thus, the long-range correlation energy is E lr c = Z d λ h h Ψ lr λ | ˆ W lr | Ψ lr λ i − h Φ | ˆ W lr | Φ i i , (11)or, equivalently, E lr c = Z d λ Z d x d x d x ′ d x ′ w lr ( x , x ; x ′ , x ′ ) × P lr c ,λ ( x , x ; x ′ , x ′ ) , (12)where w lr ( x , x ; x ′ , x ′ ) = w lr ee ( r ) δ ( x − x ′ ) δ ( x − x ′ ) − / ( N − h v lrH ( r ) δ ( x − x ′ ) + v lr x ( x , x ′ ) i δ ( x − x ′ ) is the po-tential corresponding to the perturbation operator ˆ W lr and P lr c ,λ ( x , x ; x ′ , x ′ ) is the correlation part of the two-particledensity matrix along the adiabatic connection. B. Long-range many-body perturbation theory
We now derive a formally exact many-body perturbationtheory to calculate the long-range correlation two-particledensity matrix P lr c ,λ . Details are given in Appendix A.The one-particle Green function G lr λ (1 ,
2) along the adia-batic connection of Eq. (9) in terms of space-spin-time co-ordinates 1 = ( x , t ) and 2 = ( x , t ) satisfies the followingDyson equation (cid:16) G lr λ (cid:17) − (1 , = G − (1 , − Σ lr λ (1 , − ∆Σ sr λ (1 , , (13)where G (1 ,
2) is the reference Green function correspond-ing to the RSH Hamiltonian ˆ H , Σ lr λ (1 ,
2) is the self-energycorresponding to the long-range perturbation operator λ ˆ W lr and ∆Σ sr λ (1 ,
2) is the self-energy correction associated withthe short-range potential variation term ˆ V srH xc [ n Ψ lr λ ] − ˆ V srH xc [ n Φ ]due to the variation of the density [52]. The long-range self-energy corresponding to the perturbation operator λ ( ˆ W lr ee − ˆ V lrH x , HF [ Φ ]) is decomposed into Hartree, exchange and cor-relation contributions as Σ lr λ (1 , = Σ lrH xc ,λ [ G lr λ ](1 , − Σ lrH x ,λ [ G ](1 , = λ n Σ lrH x [ G lr λ ](1 , − Σ lrH x [ G ](1 , o +Σ lr c ,λ [ G lr λ ](1 , , (14)where Σ lrH x [ G ](1 ,
2) is the sum of a long-range Hartree self-energy Σ lrH [ G ](1 , = − i Z d d w lr ee (1 , δ (1 , δ (3 , G (4 , + ) = − i δ (1 , Z d w lr ee (1 , G (3 , + ) = δ (1 , Z d r w lr ee ( r ) n ( r ) = δ (1 , v lrH [ n ]( r ) , (15)with the instantaneous electron-electron interaction w lr ee (1 , = δ ( t − t ) w lr ee ( r ) and the density extractedfrom the Green function n ( r ) = − i P s G (3 , + ) (where 3 + stands for x t + with t + = t + η and η is an infinitesimalpositive shift), and a long-range exchange self-energy Σ lr x [ G ](1 , = i Z d d w lr ee (1 , δ (1 , δ (2 , G (4 , + ) = iw lr ee (1 , G (1 , + ) = − δ ( t − t ) w lr ee ( r ) n ( x , x ) = δ ( t − t ) v lr x [ n ]( x , x ) , (16)with the one-particle density matrix extracted from the Greenfunction n ( x , x ) = − iG ( x t , x t + ). The short-range self-energy correction corresponding to the operator ˆ V srH xc [ n Ψ lr λ ] − ˆ V srH xc [ n Φ ] is written as ∆Σ sr λ (1 , = Σ srH xc [ G lr λ ](1 , − Σ srH xc [ G ](1 , , (17)where Σ srH xc [ G ](1 , = δ (1 , v srH xc [ n ]( r ) is the local short-range Hxc self-energy.The long-range four-point polarization propagator χ lr λ (1 ,
2; 1 ′ , ′ ) along the adiabatic connection is given by thesolution of the following Bethe-Salpeter-type equation whichcan be derived from the Dyson equation (13) by consideringvariations with respect to G lr λ [see Appendix A, Eq. (A13)] (cid:16) χ lr λ (cid:17) − (1 ,
2; 1 ′ , ′ ) = (cid:16) χ lrIP ,λ (cid:17) − (1 ,
2; 1 ′ , ′ ) − λ f lrH x (1 ,
2; 1 ′ , ′ ) − f lr c ,λ (1 ,
2; 1 ′ , ′ ) , (18)where χ lrIP ,λ (1 ,
2; 1 ′ , ′ ) = − iG lr λ (1 , ′ ) G lr λ (2 , ′ ) is anindependent-particle (IP) polarization propagator, and λ f lrH x (1 ,
2; 1 ′ , ′ ) = i λ δ Σ lrH x [ G lr λ ](1 , ′ ) /δ G lr λ (2 ′ ,
2) and f lr c ,λ (1 ,
2; 1 ′ , ′ ) = i δ Σ lr c ,λ [ G lr λ ](1 , ′ ) /δ G lr λ (2 ′ ,
2) are long-range Hartree-exchange and correlation kernels. Note that these kernels only stem from the self-energy term Σ lrH xc ,λ [ G lr λ ]in Eq. (13) that corresponds to the two-electron interaction λ ˆ W lr ee , the other self-energy contributions which come fromthe one-electron terms are absorbed in the definition of χ lr λ (1 ,
2; 1 ′ , ′ ). The Hartree kernel is obtained from Eq. (15) f lrH (1 ,
2; 1 ′ , ′ ) = w lr ee (1 , δ (1 , ′ ) δ (2 , ′ ) = w lr ee ( r ) δ ( t − t ) δ (1 , ′ ) δ (2 , ′ ) , (19)while the HF-like exchange kernel is obtained from Eq. (16) f lr x (1 ,
2; 1 ′ , ′ ) = − w lr ee (1 , δ (1 , ′ ) δ (1 ′ , = − w lr ee ( r ) δ ( t − t ) δ (1 , ′ ) δ (1 ′ , . (20)The fluctuation-dissipation theorem is then used to express P lr c ,λ as [see Appendix A, Eq. (A24)] P lr c ,λ ( x , x ; x ′ , x ′ ) = − Z ∞−∞ d ω π i e i ω + h χ lr λ ( x , x ; x ′ , x ′ ; ω ) − χ ( x , x ; x ′ , x ′ ; ω ) i +∆ lr λ ( x , x ; x ′ , x ′ ) , (21)where χ lr λ ( x , x ; x ′ , x ′ ; ω ) is the frequency-dependent Fouriertransform of the one-time-interval polarization propagator χ lr λ ( x , x ; x ′ , x ′ ; τ = t − t ) = χ lr λ ( x t , x t ; x ′ t + , x ′ t + ), χ ( x , x ; x ′ , x ′ ; ω ) is the equivalent quantity for the RSH ref-erence Hamiltonian (at λ = ∆ lr λ ( x , x ; x ′ , x ′ ) is thecontribution coming from the variation of the one-particledensity matrix along the adiabatic connection. The expressionof ∆ lr λ in terms of the Green functions G lr λ and G is straightfor-ward but it is su ffi cient to write it as ∆ lr λ = Γ [ G lr λ ] − Γ [ G ] where Γ is a known functional given in Appendix A [Eq. (A22)].So far, the theory is in principle exact . In the following weconsider two possible approximations. The RPA approxima-tion Σ lr xc ,λ = , (22)corresponds to neglecting long-range exchange-correlation inall one-electron properties. Indeed, with this approximation,one can check that G lr λ = G is a solution of the Dyson equa-tion (13), i.e. the Green function remains unchanged alongthe adiabatic connection. It follows that ∆ lr λ = f lr xc ,λ = χ lrIP ,λ (1 ,
2; 1 ′ , ′ ) = − iG (1 , ′ ) G (2 , ′ ) = χ (1 ,
2; 1 ′ , ′ ). Sim-ilarly, the RPAx approximation Σ lr c ,λ = , (23)corresponds to neglecting long-range correlation only in allone-electron properties. Again, this approximation impliesthat the Green function remains unchanged along the adiabaticconnection, i.e. G lr λ = G and it follows that ∆ lr λ = f lr c ,λ = χ lrIP ,λ = χ . As di ff erent terminologies are used in the quan-tum chemistry and condensed-matter physics literature, let usstress that what we call RPA here corresponds to a responseequation (18) with no exchange-correlation kernel (and it isalso sometimes called linear response time-dependent Hartreetheory or direct RPA), and what we call RPAx correspondsto a response equation with an additional HF-like exchangekernel (and it is also sometimes called linear response time-dependent Hartree-Fock theory or full
RPA).
C. Expressions in an orbital basis
The RPA or RPAx equations in an orbital basis are derivedin details in Appendix B. In the basis of RSH spin orbitals,the long-range RPA or RPAx correlation energy writes E lr c = Z d λ X ia , jb h ib | ˆ w lr ee | a j i ( P lr c ,λ ) ia , jb , (24)where i and j refer to occupied spin orbitals, and a and b tovirtual spin orbitals, h ib | ˆ w lr ee | a j i are the two-electron integralswith long-range interaction, and ( P lr c ,λ ) ia , jb are the matrix el-ements of the correlation two-particle density matrix. Theone-electron terms v lrH and v lr x in the perturbation operator inEq. (12) do not contribute to E lr c because of the occupied-virtual / occupied-virtual structure of the two-particle densitymatrix in RPA or RPAx. Following the technique proposed byFurche [26], P lr c ,λ can be obtained as P lr c ,λ = (cid:16) A lr λ − B lr λ (cid:17) / (cid:16) M lr λ (cid:17) − / (cid:16) A lr λ − B lr λ (cid:17) / − . (25)with M lr λ = (cid:16) A lr λ − B lr λ (cid:17) / (cid:16) A lr λ + B lr λ (cid:17) (cid:16) A lr λ − B lr λ (cid:17) / , and the or-bital rotation Hessians (cid:16) A lr λ (cid:17) ia , jb = ( ǫ a − ǫ i ) δ i j δ ab + λ h h ib | ˆ w lr ee | a j i − ξ h ib | ˆ w lr ee | ja i i , (26a) (cid:16) B lr λ (cid:17) ia , jb = λ h h ab | ˆ w lr ee | i j i − ξ h ab | ˆ w lr ee | ji i i . (26b)where ǫ i are the RSH orbital eigenvalues, and ξ = ξ = E lr c = Z d λ X ia , jb h ib | ˆ w lr ee | a j i ( P lr c ,λ ) ia , jb , (27)where i and j now refer to occupied spatial orbitals, and a and b to virtual spatial orbitals, and P lr c ,λ is the spin-singlet-adapted correlation two-particle density matrix obtained as P lr c ,λ = (cid:20)(cid:16) A lr λ − B lr λ (cid:17) / (cid:16) M lr λ (cid:17) − / (cid:16) A lr λ − B lr λ (cid:17) / − (cid:21) , (28)with M lr λ = (cid:16) A lr λ − B lr λ (cid:17) / (cid:16) A lr λ + B lr λ (cid:17) (cid:16) A lr λ − B lr λ (cid:17) / , andthe singlet orbital rotation Hessians (cid:16) A lr λ (cid:17) ia , jb = ( ǫ a − ǫ i ) δ i j δ ab + λ h h ib | ˆ w lr ee | a j i − ξ h ib | ˆ w lr ee | ja i i , (29a) (cid:16) B lr λ (cid:17) ia , jb = λ h h ab | ˆ w lr ee | i j i − ξ h ab | ˆ w lr ee | ji i i . (29b)Only singlet excitations contribute to Eq. (27), since the two-electron integrals involved vanish for triplet excitations. In Eq. (25), it is assumed that A lr λ + B lr λ and A lr λ − B lr λ are pos-itive definite. In RPA, this is always the case. On the contrary,in RPAx, this is not always the case, i.e. instabilities can beencountered, and Eq. (25) can fail. In spin-restricted closed-shell formalism, one may encounter singlet instabilities in theRPAx theory defined here, for example when dissociating abond, but not triplet instabilities since triplet excitations donot contribute at all. In practice, singlet instabilities are usu-ally not encountered for weakly-interacting closed-shell sys-tems. Note that other variants of RPA-type correlation energyexpressions using a HF exchange response kernel, such as theplasmon formula [38, 53, 54] or the equivalent ring coupled-cluster-doubles theory [38], require contributions from bothsinglet and triplet excitations, and are thus subject to tripletinstabilities (e.g. in a system such as Be ).Similarly to the notation used in Ref. 20, the range-separated method obtained by adding to the RSH energy thelong-range RPAx correlation energy [ ξ = + lrRPAx. For con-sistency, the range-separated method obtained by addingto the RSH energy the long-range RPA correlation energy[ ξ = + lrRPA, although it is equivalent to the method called“LC- ω LDA + dRPA” in Refs. 21–24 in the special case of theshort-range LDA functional. At second order in the electron-electron interaction, the RSH + lrRPAx method reduces to therange-separated method of Ref. 9 based on long-range second-order Møller-Plesset perturbation theory, to which we willrefer as RSH + lrMP2. Since RPA approaches can be seenas simple approximations to coupled-cluster theory [38], theRSH + lrRPA and RSH + lrRPAx methods bear some resem-blance to the range-separated method of Ref. 14 where thelong-range correlation energy is evaluated by coupled-clustertheory (with single, double and perturbative triple excita-tions), to which we will refer as RSH + lrCCSD(T).We note that one can develop long-rang many-body pertur-bation theories starting from other references than the RSHreference. For example, starting from the usual (approximate)Kohn-Sham reference could be appropriate for solid-state sys-tems. For the finite systems considered here, RSH is a goodreference, as confirmed by other authors [23]. III. COMPUTATIONAL DETAILS
All calculations have been performed with a developmentversion of MOLPRO 2008 [55], implementing equations (27)-(29). We first perform a self-consistent RSH calculation withthe short-range PBE xc functional of Ref. 14 (this RSH cal-culation could also be referred to as “lrHF + srPBE”, a nota-tion closer to the one used by other authors [14]) and addthe long-range MP2, RPA, RPAx or CCSD(T) correlation en-ergy calculated with RSH orbitals. For RPA or RPAx, the λ -integration in Eq. (27) is done by a 7-point Gauss-Legendrequadrature [26]. The range separation parameter is taken at µ = . − , in agreement with previous studies [56], with-out trying to adjust it for each system. To show the depen-dence on the orbitals, the full-range RPA calculations havebeen done with PBE [57] and HF orbitals, which will be de-noted by PBE + RPA and HF + RPA, respectively [58]. Thefull-range MP2, RPAx and CCSD(T) calculations have beendone with HF orbitals, and thus, for notation consistency, willbe denoted by HF + MP2, HF + RPAx and HF + CCSD(T), re-spectively. We use large Dunning basis sets [59–65]. Coreelectrons are kept frozen in all the full-range and range-separated MP2, RPA, RPAx and CCSD(T) calculations (i.e.only excitations of valence electrons are considered). The ba-sis set superposition error (BSSE) is removed by the coun-terpoise method. For the alkaline-earth dimers, it has beenchecked than adding di ff use basis functions or core excita-tions do not change significantly the results. Extrapolationsto the complete basis set (CBS) limit have also been consid-ered for some systems. For the full-range methods, the stan-dard three-point exponential formula for the HF (or KS) ref-erence E HF ( n ) = E HF (CBS) + Ae − Bn with the cardinal number n = , ,
5, and two-point formula for the correlation energy E c ( n ) = E c (CBS) + C / n with n = , C dispersioncoe ffi cients, the interaction energy E int is calculated at sevenextra distances R i from 30 to 60 bohr, and the coe ffi cient isestimated by averaging with the following formula C = exp X i = (ln | E int ( R i ) | + R i )) , (30)similarly to what has been done in Ref. 22. IV. APPLICATIONSA. Basis set dependence
The convergence of the equilibrium binding energy ofAr with respect to the basis set size up to the CBSlimit for the full-range methods HF + MP2, PBE + RPA,HF + RPA, HF + CCSD(T) and for the range-separated methodsRSH + lrMP2, RSH + lrRPA, RSH + lrRPAx, RSH + lrCCSD(T)is represented in Fig. 1. Full-range RPA with PBE orbitals hasa very strong dependence on the basis size, as already noted(e.g. Refs. 20, 26). Full-range RPA with HF orbitals has abit weaker basis dependence, similar to full-range HF + MP2,HF + RPAx and HF + CCSD(T). All the range-separated meth-ods have essentially identical, very favorable basis set conver-gence. Since the slow convergence of full-range methods is
60 65 70 75 80 85 90 95 100 aVTZ aVQZ aV5Z CBS P erce n t ag e o f C B S b i nd i n g e n er gy One-particle basisHF+MP2PBE+RPAHF+RPAHF+RPAxHF+CCSD(T)RSH+lrMP2RSH+lrRPARSH+lrRPAxRSH+lrCCSD(T)Ar FIG. 1: (Color online) Basis set dependence of the equilibrium bind-ing energy of Ar for di ff erent full-range and range-separated meth-ods, presented as the percentage of the binding energy recovered withrespect to the CBS limit (aVTZ, aVQZ and aV5Z stand for aug-cc-pVTZ, aug-cc-pVQZ and aug-cc-pV5Z, respectively). related to the explicit description of short-range correlation,it is not surprising that range-separated methods have a fasterconvergence because they leave the description of short-rangecorrelation to the short-range density functional. These resultsare consistent with other studies, e.g. Refs. 22, 24. Note that,with the aug-cc-pV5Z basis set, all the range-separated meth-ods are essentially converged (98-99% of the CBS bindingenergy), therefore we will not use CBS extrapolations in thefollowing. However, one should keep in mind that with thisbasis set the full-range methods are not yet fully converged,with about 90% of the CBS binding energy. B. Rare-gas dimers
In Fig. 2, the interaction energy curves of He , Ne , Ar andKr , obtained with the full-range and range-separated meth-ods are compared. As already known, full-range HF + MP2underestimates the interaction energy for the smallest sys-tems He and Ne , and overestimates it for the largest sys-tems Ar and Kr . Full-range PBE + RPA gives an almost dis-sociative curve for He , and largely underestimates the in-teraction energy for Ne , Ar and Kr . Using HF orbitalsin full-range RPA drastically improves the interaction energycurve for He , and to a least extend for Ne , but gives lessbinding for Ar and Kr . Full-range HF + RPAx significantlyimproves over full-range HF + RPA, but still gives underesti-mated interaction energies. It can be noted that full-rangeHF + RPAx yields interaction energy curves almost identical tothe full-range HF + MP2 curves for He and Ne , and almostidentical to the full-range PBE + RPA curves for Ar and Kr .Full-range HF + CCSD(T) gives systematically quite accurateinteraction energies. Quite similarly to full-range HF + MP2,the range-separated RSH + lrMP2 underestimates the interac-tion energy for He and Ne , and overestimates it for Ar and Kr . RSH + lrRPA tends to improve over both full-rangePBE + RPA and HF + RPA but still leads to significantly under-estimated interaction energies. RSH + lrRPAx improves over -0.04-0.03-0.02-0.01 0 0.01 5 6 7 8 9 10 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateHF+MP2PBE+RPAHF+RPAHF+RPAxHF+CCSD(T)He -0.04-0.03-0.02-0.01 0 0.01 5 6 7 8 9 10 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateRSH+lrMP2RSH+lrRPARSH+lrRPAxRSH+lrCCSD(T)He -0.15-0.1-0.05 0 0.05 5 6 7 8 9 10 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateHF+MP2PBE+RPAHF+RPAHF+RPAxHF+CCSD(T)Ne -0.15-0.1-0.05 0 0.05 5 6 7 8 9 10 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateRSH+lrMP2RSH+lrRPARSH+lrRPAxRSH+lrCCSD(T)Ne -0.6-0.5-0.4-0.3-0.2-0.1 0 0.1 6 7 8 9 10 11 12 13 14 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateHF+MP2PBE+RPAHF+RPAHF+RPAxHF+CCSD(T)Ar -0.6-0.5-0.4-0.3-0.2-0.1 0 0.1 6 7 8 9 10 11 12 13 14 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateRSH+lrMP2RSH+lrRPARSH+lrRPAxRSH+lrCCSD(T)Ar -0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 0.1 6 8 10 12 14 16 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateHF+MP2PBE+RPAHF+RPAHF+RPAxHF+CCSD(T)Kr -0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 0.1 6 8 10 12 14 16 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateRSH+lrMP2RSH+lrRPARSH+lrRPAxRSH+lrCCSD(T)Kr FIG. 2: (Color online) Interaction energy curves of He , Ne , Ar and Kr calculated by di ff erent full-range (left) and range-separated (right)methods. The basis is aug-cc-pV5Z. The accurate curves are from Ref. 66. TABLE I: Hard-core radii σ (bohr), equilibrium distances R e (bohr), equilibrium binding energies D e (mhartree), harmonic vibrational fre-quencies ω e (cm − ) and dispersion coe ffi cients C for ten homonuclear and heteronuclear rare-gas dimers from di ff erent full-range and range-separated methods with aug-cc-pV5Z basis. Mean absolute percentage errors (MA%E) are also given. HF + MP2 PBE + RPA HF + RPA HF + RPAx HF + CCSD(T) RSH + lrMP2 RSH + lrRPA RSH + lrRPAx RSH + lrCCSD(T) Estimated exact a He σ R e D e ω e C σ R e D e ω e C σ R e D e ω e C σ R e D e ω e C σ R e D e ω e C σ R e D e ω e C σ R e D e ω e C σ R e D e ω e C σ R e D e ω e C σ R e D e ω e C
159 116 86 105 132 162 109 134 163 129.6MA%E (%) σ R e D e
23 62 56 39 10 16 36 14 11 0.0 ω e
12 36 33 21 3.4 9.5 23 10 4.6 0.0 C
13 7.0 36 22 4.1 14 9.2 10 29 0.0 a From Ref. 66 both RSH + lrRPA and full-range HF + RPAx; it still systemati-cally underestimates the interaction energy at equilibrium, butappears quite accurate at medium and large distances. On thecontrary, RSH + lrCCSD(T) systematically overestimates the interaction energy at medium and large distances.The hard-core radii, equilibrium distances, equilibriumbinding energies, harmonic vibrational frequencies and dis-persion coe ffi cients C for ten homonuclear and heteronu- -6-4-2 0 2 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateHF+MP2PBE+RPAHF+RPAHF+RPAxHF+CCSD(T)Be -6-4-2 0 2 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateRSH+lrMP2RSH+lrRPARSH+lrRPAxRSH+lrCCSD(T)Be -2-1.5-1-0.5 0 0.5 6 7 8 9 10 11 12 13 14 15 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateHF+MP2PBE+RPAHF+RPAHF+RPAxHF+CCSD(T)Mg -2-1.5-1-0.5 0 0.5 6 7 8 9 10 11 12 13 14 15 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateRSH+lrMP2RSH+lrRPARSH+lrRPAxRSH+lrCCSD(T)Mg -5-4-3-2-1 0 1 6 7 8 9 10 11 12 13 14 15 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateHF+MP2PBE+RPAHF+RPAHF+RPAxHF+CCSD(T)Ca -5-4-3-2-1 0 1 6 7 8 9 10 11 12 13 14 15 I n t er a c t i o n e n er gy ( m h a r t ree ) Internuclear distance (bohr)AccurateRSH+lrMP2RSH+lrRPARSH+lrRPAxRSH+lrCCSD(T)Ca FIG. 3: (Color online) Interaction energy curves of Be , Mg and Ca calculated by full-range (left) and range-separated (right) methods. Thebasis is cc-pV5Z. The accurate curves are from Refs. 67, 68 and 69. clear rare-gas dimers calculated with the full-range and range-separated methods are given in Table I. The trends seen inFig. 2 are confirmed. Full-range RPA (with PBE or HF or-bitals) yields very inaccurate equilibrium properties. Full-range HF + RPAx improves over full-range HF + RPA (withthe exception of C coe ffi cients which turn out to be quitegood in PBE + RPA for these systems) but the errors remainlarge. Range separation largely improves RPA and RPAx.RSH + lrRPAx gives much better equilibrium properties thanRSH + lrRPA, with mean absolute percentage errors smallerby more than a factor of two, while these two methods givesimilar accuracy for C coe ffi cients. Full-range HF + MP2 isreasonably accurate and range separation has a much smallerimpact on it. For these systems, RSH + lrMP2 gives an over- all similar accuracy than RSH + RPAx, although the C coef-ficients tend to be globally more accurate in RSH + lrRPAx.Full-range HF + CCSD(T) gives the best results. Surpris-ingly, range separation tends to deteriorate the accuracy ofCCSD(T), especially for C coe ffi cients. Nevertheless, amongthe range-separated methods, RSH + lrCCSD(T) still gives thebest equilibrium properties. C. Alkaline-earth dimers
In Fig. 3, the interaction energy curves of Be , Mg and Ca , obtained with the full-range and range-separatedmethods are compared. These systems have static corre- TABLE II: Hard-core radii σ (bohr), equilibrium distances R e (bohr), equilibrium binding energies D e (mhartree), harmonic vibrationalfrequencies ω e (cm − ) and dispersion coe ffi cients C for Be , Mg and Ca from di ff erent full-range and range-separated methods with cc-pV5Zbasis. Mean absolute percentage errors (MA%E) are also given. HF + MP2 PBE + RPA HF + RPA HF + RPAx HF + CCSD(T) RSH + lrMP2 RSH + lrRPA RSH + lrRPAx RSH + lrCCSD(T) Estimated exactBe σ a R e a D e a ω e
139 297 34 37 242 199 152 198 315 267 a C
256 164 138 180 195 232 149 213 274 214 d Mg σ b R e b D e b ω e
47 7.9 31 35 48 45 30 42 52 51.1 b C
686 405 364 485 616 571 349 494 671 627 d Ca σ c R e c D e c ω e
56 — 44 47 64 60 50 57 68 63.7 c C d MA%E (%) σ R e D e
32 — 69 61 19 26 63 33 21 0.0 ω e
23 — 53 48 5.3 14 35 18 9.1 0.0 C
15 33 40 21 5.0 7.7 41 16 12 0.0 a From Ref. 67 b From Ref. 68 c From Ref. 69 d From Ref. 70 lation e ff ects, especially Be , and are thus more challeng-ing for the single-reference methods tested here. Full-rangePBE + RPA gives unphysical interaction energy curves, witha large bump for Be , and with essentially no bond forMg and Ca . Full-range HF + RPA yields an almost dis-sociative curve for Be with no bump (which is consistentwith Ref. 43), and physically reasonable curves for Mg andCa . Full-range HF + RPAx moderately improves over full-range HF + RPA. Among the full-range methods, HF + MP2and HF + CCSD(T) clearly give the best interaction energycurves. As for rare-gas dimers, RSH + lrRPA always largelyunderestimates the interaction energy. RSH + lrMP2 andRSH + lrRPAx give much less underestimated interaction en-ergies, with RSH + lrMP2 being a bit more accurate for Mg and Ca . While RSH + lrCCSD(T) largely overestimates theinteraction energy for Be , it is remarkably accurate for Mg and Ca . We note that RSH + lrCCSD(T) could be made moreaccurate for Be by choosing a larger range-separation param-eter µ [71].The hard-core radii, equilibrium distances, equilibriumbinding energies, harmonic vibrational frequencies and dis-persion coe ffi cients C for Be , Mg and Ca are given inTable II. It is confirmed that range separation largely im-proves the equilibrium properties of RPA and RPAx. Again,RSH + lrRPAx is much more accurate than RSH + lrRPA, withmean absolute percentage errors smaller by about a factor oftwo. Range separation also overall brings a significant im-provement in MP2. Among the range-separated methods,RSH + lrCCSD(T) gives the best equilibrium properties. V. CONCLUSIONS
We have expounded the details of a formally exactadiabatic-connection fluctuation-dissipation density-functional theory based on range separation. Range-separateddensity-functional theory with random phase approximationsincluding or not the long-range Hartree-Fock exchange re-sponse kernel (referred to as RSH + lrRPA and RSH + lrRPAx,respectively) are then obtained as well-identified approx-imations on the long-range Green-function self-energy[Eqs. (22) and (23)]. The long-range Green function does notvary along the adiabatic connection at the RSH + lrRPA andRSH + lrRPAx levels, which makes these schemes relativelysimple compared to the exact theory. In practice, RSH + lrRPAand RSH + lrRPAx have been applied in a spin-restrictedclosed-shell formalism, in which both schemes only includespin-singlet orbital excitations, and thus are not subject totriplet instabilities.These range-separated RPA-type schemes have been testedon rare-gas and alkaline-earth dimers, featuring challengingweak (van der Waals) interactions. Both range separation andinclusion of the exact Hartree-Fock response kernel largelyimprove the accuracy of RPA. The RSH + lrRPAx method ap-pears as a reasonably accurate method for weak interactions,but globally less accurate for equilibrium properties thanthe more intensive range-separated coupled-cluster method.Although, for the small systems considered here, range-separated second-order perturbation theory (RSH + lrMP2)turns out to yield results similarly as accurate as those fromRSH + lrRPAx (and in fact more accurate for Mg and Ca ),a recent investigation [72] shows that RSH + lrRPAx corrects0the overestimation of the binding energy in RSH + lrMP2 forlarger weakly-interacting stacked complexes, such as the ben-zene dimer. Acknowledgments
We thank J. F. Dobson and T. Gould for numerous discus-sions on RPA during the FAST (French-Australian Scienceand Technology) workshop at Gri ffi th University, Australia.We also thank G. Jansen for discussions. This work was partlysupported by ANR (French national research agency) via con-tract number 07-BLAN-0272 (Wademecom). Appendix A: Adiabatic-connection fluctuation-dissipationdensity-functional theory
In this appendix, we outline a general, formallyexact adiabatic-connection fluctuation-dissipation density-functional theory, using Green-function many-body theory.For further details on standard Green’s function theory, seee.g. Refs. 73–76.
1. Adiabatic connection
We consider the following adiabatic connection defined bythe λ -dependent energy E λ = min Ψ n h Ψ | ˆ K + λ ˆ W | Ψ i + F [ n Ψ ] o , (A1)where ˆ K is an arbitrary one-particle Hamiltonian, ˆ W is aperturbation operator (generally, the sum of a two-particleoperator ˆ W ee and an one-particle operator) and F [ n ] is a λ -independent density functional. The minimizing multideter-minant wave function Ψ λ satisfies the Euler-Lagrange equa-tion ˆ H λ | Ψ λ i = E λ | Ψ λ i , (A2)where E λ is the Lagrange multiplier for the normalization con-straint, and ˆ H λ is the e ff ective Hamiltonian along the adiabaticconnection ˆ H λ = ˆ K + λ ˆ W + ˆ V λ , (A3)where ˆ V λ = R d r δ F [ n Ψ λ ] /δ n ( r ) ˆ n ( r ) is a self-consistent one-particle potential operator. Note that ˆ H λ = is not necessar-ily the physical Hamiltonian. This adiabatic connection linksthe energy of interest E λ = to the reference energy E λ = = h Φ | ˆ K | Φ i + F [ n Φ ] calculated with the single-determinantwave function Φ = Ψ λ = of the reference Hamiltonian ˆ H = ˆ K + ˆ V . The one-particle density is not kept constant withrespect to λ .An adiabatic connection formula for E λ = is found by tak-ing the derivative of E λ with respect to λ , noting that E λ is stationary with respect to Ψ λ , and reintegrating between λ = λ = E λ = = E λ = + Z d λ h Ψ λ | ˆ W | Ψ λ i . (A4)The correlation energy, defined as E c = E λ = − E λ = − ( dE λ / d λ ) λ = where ( dE λ / d λ ) λ = = h Φ | ˆ W | Φ i is the first-order energy correction, is thus given by E c = Z d λ h h Ψ λ | ˆ W | Ψ λ i − h Φ | ˆ W | Φ i i , (A5)or, equivalently, in the representation of space-spin coordi-nates x = ( r , s ) E c = Z d λ Z d x d x d x ′ d x ′ w ( x , x ; x ′ , x ′ ) × P c ,λ ( x , x ; x ′ , x ′ ) , (A6)where w ( x , x ; x ′ , x ′ ) is the interaction potential correspond-ing to the operator ˆ W and P c ,λ ( x , x ; x ′ , x ′ ) is the correlationpart of the two-particle density matrix along the adiabatic con-nection.This exposition encompasses both standard full-rangemany-body theory and range-separated density-functionaltheory. Indeed, if ˆ K is the Hartree-Fock Hamiltonian (i.e.,ˆ K = ˆ T + ˆ V ne + ˆ V H x , HF ), ˆ W is the standard Møller-Plessetfluctuation perturbation operator (i.e., ˆ W = ˆ W ee − ˆ V H x , HF ) and F [ n ] = K = ˆ T + ˆ V ne + ˆ V lrH x , HF and ˆ W = ˆ W lr ee − ˆ V lrH x , HF and the short-range density functional F [ n ] = E srH xc [ n ], Eq. (A6) yields nowthe long-range correlation energy, defined with respect to theRSH energy [Eq. (5)].
2. One-particle Green function
The one-particle Green function along the adiabatic con-nection is defined as G λ (1 , = − i h Ψ λ | T [ ˆ ψ λ (1) ˆ ψ † λ (2)] | Ψ λ i , (A7)where 1 = ( x , t ) and 2 = ( x , t ) refer to space-spin andtime coordinates, ˆ ψ λ (1) = e i ˆ H λ t ˆ ψ ( x ) e − i ˆ H λ t and ˆ ψ † λ (2) = e i ˆ H λ t ˆ ψ † ( x ) e − i ˆ H λ t are the annihilation and creation operatorsin the Heisenberg picture, and T is the Wick time-orderingoperator.A Dyson-type equation connects the inverse of G λ to theinverse of the Green function associated with the one-electronHamiltonian ˆ K + ˆ V λ , denoted by G V ,λ , G − λ (1 , = G − V ,λ (1 , − Σ λ (1 , , (A8)which can be considered as the definition of the self-energy Σ λ . In turn, the inverse of G V ,λ can be expressed from the in-verse of the Green function G of the reference Hamiltonian1ˆ H = ˆ K + ˆ V as G − V ,λ = G − − [ v λ − v ], where v λ and v arethe one-electron potentials associated with ˆ V λ and ˆ V , respec-tively.For time-independent Hamiltonians, the Green functiononly depends on the time di ff erence τ = t − t , so one de-fines G λ ( x , x ; τ ) = G λ ( x t , x t ), which has a discontinu-ity at τ =
0. The one-particle density matrix n ,λ ( x , x ) = h Ψ λ | ˆ n ( x , x ) | Ψ λ i , with ˆ n ( x , x ) = ˆ ψ † ( x ) ˆ ψ ( x ), can be ob-tained from the limit τ → − n ,λ ( x , x ) = − iG λ ( x , x ; τ = − ) . (A9)
3. Four-point polarization propagator
The four-point polarization propagator along the adiabaticconnection is defined as χ λ (1 ,
2; 1 ′ , ′ ) = i (cid:2) G ,λ (1 ,
2; 1 ′ , ′ ) − G λ (1 , ′ ) G λ (2 , ′ ) (cid:3) , (A10)where G ,λ is the two-particle Green function G ,λ (1 ,
2; 1 ′ , ′ ) = −h Ψ λ | T [ ˆ ψ λ (1) ˆ ψ λ (2) ˆ ψ † λ (2 ′ ) ˆ ψ † λ (1 ′ )] | Ψ λ i , (A11)Alternatively, using the Schwinger derivative technique, χ λ can be expressed as the functional derivative of the one-particle Green function with respect to the two-point potential v λ (see, e.g., Refs. 73, 76) χ λ (1 ,
2; 1 ′ , ′ ) = − i δ G V ,λ (1 , ′ ) δ v λ (2 ′ , . (A12)The four-point polarization propagator satisfies a so-calledBethe-Salpeter equation that directly stems from the Dysonequation of Eq. (A8). Considering variations with respect to iG λ (achieved through variations of v λ ) yields − i δ G − λ (1 , ′ ) δ G λ (2 ′ , = − i δ G − V ,λ (1 , ′ ) δ G λ (2 ′ , + i δ Σ λ (1 , ′ ) δ G λ (2 ′ , . (A13)The term on the left-hand side of Eq. (A13) gives straightfor-wardly − i δ G − λ (1 , ′ ) δ G λ (2 ′ , = iG − λ (1 , ′ ) G − λ (2 , ′ ) = χ − ,λ (1 ,
2; 1 ′ , ′ ) , (A14)where χ IP ,λ (1 ,
2; 1 ′ , ′ ) = − iG λ (1 , ′ ) G λ (2 , ′ ) is a so-calledindependent-particle (IP) polarization propagator [77]. Thefirst term on the right-hand side of Eq. (A13) gives the in-verse of the four-point polarization propagator, according toEq. (A12), − i δ G − V ,λ (1 , ′ ) δ G λ (2 ′ , = i δ v λ (1 , ′ ) δ G λ (2 ′ , = χ − λ (1 ,
2; 1 ′ , ′ ) , (A15)and the second term is the so-called Bethe-Salpeter four-pointkernel i δ Σ λ (1 , ′ ) δ G λ (2 ′ , = f λ (1 ,
2; 1 ′ , ′ ) , (A16) and finally, using Eqs. (A14)-(A16) in Eq. (A13), the Bethe-Salpeter equation for χ λ writes χ − λ (1 ,
2; 1 ′ , ′ ) = χ − ,λ (1 ,
2; 1 ′ , ′ ) − f λ (1 ,
2; 1 ′ , ′ ) . (A17)
4. Fluctuation-dissipation theorem
Similarly to the expression of the one-particle density ma-trix in terms of the one-particle Green function [Eq. (A9)],the two-particle density matrix can be extracted from thepolarization propagator. Defining χ λ ( x , x ; x ′ , x ′ ; τ ) = χ λ ( x t , x t ; x ′ t + , x ′ t + ), i.e. the polarization propagator withtimes t ′ → t + and t ′ → t + which depends only on the time dif-ference τ = t − t , it is easy to check that in the limit τ → − ,after applying the time-ordering operator in Eq. (A11) and us-ing Eq. (A9), one has the following relation i χ λ ( x , x ; x ′ , x ′ ; τ = − ) = h Ψ λ | ˆ n ( x , x ′ )ˆ n ( x , x ′ ) | Ψ λ i− n ,λ ( x , x ′ ) n ,λ ( x , x ′ ) . (A18)The two-particle density matrix n ,λ ( x , x ; x ′ , x ′ ) = h Ψ λ | ˆ ψ † ( x ′ ) ˆ ψ † ( x ′ ) ˆ ψ ( x ) ˆ ψ ( x ) | Ψ λ i can thus be expressed as n ,λ ( x , x ; x ′ , x ′ ) = h Ψ λ | ˆ n ( x , x ′ )ˆ n ( x , x ′ ) | Ψ λ i− δ ( x ′ − x ) n ,λ ( x , x ′ ) = i χ λ ( x , x ; x ′ , x ′ ; τ = − ) + n ,λ ( x , x ′ ) n ,λ ( x , x ′ ) − δ ( x ′ − x ) n ,λ ( x , x ′ ) . (A19)The correlation part of the two-particle density matrix P c ,λ = n ,λ − n ,λ = is thus P c ,λ ( x , x ; x ′ , x ′ ) = i χ λ ( x , x ; x ′ , x ′ ; τ = − ) − i χ ( x , x ; x ′ , x ′ ; τ = − ) +∆ λ ( x , x ; x ′ , x ′ ) , (A20)where χ is the polarization propagator of the non-interactingreference system for λ =
0, and ∆ λ is a term coming from thevariation of the one-particle density matrix along the adiabaticconnection ∆ λ ( x , x ; x ′ , x ′ ) = n ,λ ( x , x ′ ) n ,λ ( x , x ′ ) − δ ( x ′ − x ) n ,λ ( x , x ′ ) − n , ( x , x ′ ) n , ( x , x ′ ) + δ ( x ′ − x ) n , ( x , x ′ ) . (A21)Using Eq. (A9), one can also express this term with the Greenfunction as ∆ λ = Γ [ G λ ] − Γ [ G ] where we define the functional Γ as Γ [ G ] = − G ( x , x ′ ; τ = − ) G ( x , x ′ ; τ = − ) + δ ( x ′ − x ) iG ( x , x ′ ; τ = − ) . (A22)Finally, introducing the Fourier transform of χ λ ( x , x ; x ′ , x ′ ; τ ) in terms of the frequency ω , i χ λ ( x , x ; x ′ , x ′ ; τ = − ) = − Z ∞−∞ d ω π i e i ω + × χ λ ( x , x ; x ′ , x ′ ; ω ) , (A23)2we arrive at the form of the fluctuation-dissipation that we use P c ,λ ( x , x ; x ′ , x ′ ) = − Z ∞−∞ d ω π i e i ω + h χ λ ( x , x ; x ′ , x ′ ; ω ) − χ ( x , x ; x ′ , x ′ ; ω ) i +∆ λ ( x , x ; x ′ , x ′ ) . (A24) Appendix B: Random phase approximation in an orbital basis
In this appendix, we give the working equations in an or-bital basis resulting from the many-body theory outlined inAppendix A, in the special case of a random phase approxi-mation (RPA)-type simplification. For further details, see e.g.Refs. 26, 53, 78, 79.
1. Expressions in a spin-orbital basis
In the RPA and RPAx approximations, the Green functiondoes not vary along the adiabatic connection, i.e. G λ = G ,which implies that the independent-particle polarization prop-agator [Eq. (A14)] is just the non-interacting reference polar-ization propagator, χ IP ,λ (1 ,
2; 1 ′ , ′ ) = − iG (1 , ′ ) G (2 , ′ ) = χ (1 ,
2; 1 ′ , ′ ), and in the fluctuation-dissipation theorem ofEq. (A24) the term coming from the variation of the one-particle density matrix vanishes, ∆ λ = χ ( x , x ; x ′ , x ′ ; ω ) = X ia φ ∗ i ( x ′ ) φ a ( x ) φ ∗ a ( x ′ ) φ i ( x ) ω − ( ǫ a − ǫ i ) + i + − X ia φ ∗ i ( x ′ ) φ a ( x ) φ ∗ a ( x ′ ) φ i ( x ) ω + ( ǫ a − ǫ i ) − i + , (B1)where φ p ( x ) and ǫ p are the spin orbitals and correspondingeigenvalues of the reference system, and i and a run over oc-cupied and virtual spin orbitals, respectively. Hence, χ canbe completely represented in the basis of spin-orbital prod-ucts, φ ∗ p ( x ′ ) φ q ( x ), where p refer to an occupied orbital and q to a virtual orbital, and vice versa, with matrix elements( (cid:5) ( ω )) pq , rs = Z d x d x d x ′ d x ′ φ p ( x ′ ) φ ∗ q ( x ) × χ ( x , x ; x ′ , x ′ ; ω ) φ ∗ r ( x ) φ s ( x ′ ) . (B2)Assuming orthonormality of the spin orbitals, the matrix ele-ments are easily calculated( (cid:5) ( ω )) ia , jb = δ i j δ ab ω − ( ǫ a − ǫ i ) + i + , (B3a)( (cid:5) ( ω )) ai , b j = − δ i j δ ab ω + ( ǫ a − ǫ i ) − i + , (B3b)( (cid:5) ( ω )) ai , jb = ( (cid:5) ( ω )) ia , b j = , (B3c) where both i and j refer to occupied orbitals and both a and b to virtual orbitals. The matrix is thus diagonal, and the inverseof χ has the following 2 × (cid:5) ( ω ) − = − " ∆ ǫ ∆ ǫ ! − ω − ! , (B4)where ∆ ǫ ia , jb = ( ǫ a − ǫ i ) δ i j δ ab , each block matrices being re-indexed with the composite indices ia and jb .In the RPA and RPAx approximations, the Bethe-Salpeterkernel of Eq. (A16) is approximated as the frequency-independent Hartree(-Fock) form [Eqs. (19) and (20)] f λ ( x , x ; x ′ , x ′ ) = λ w ee ( r )[ δ ( x − x ′ ) δ ( x − x ′ ) − ξ δ ( x − x ′ ) δ ( x ′ − x )] , (B5)where w ee ( r ) is a two-particle interaction, and ξ = ξ = F λ ) pq , rs = Z d x d x d x ′ d x ′ φ p ( x ′ ) φ ∗ q ( x ) × f λ ( x , x ; x ′ , x ′ ) φ ∗ r ( x ) φ s ( x ′ ) = λ (cid:2) h qr | ˆ w ee | ps i − ξ h qr | ˆ w ee | sp i (cid:3) , (B6)where h qr | ˆ w ee | ps i are the two-electron integrals. The su-permatrix representation of the interacting polarization prop-agator χ λ is then found from the Bethe-Salpeter equation[Eq. (A17)] written in the spin-orbital basis (cid:5) λ ( ω ) − = (cid:5) ( ω ) − − F λ = − " A λ B λ B ∗ λ A ∗ λ ! − ω − ! , (B7)where A λ and B λ are the so-called orbital rotation Hessians( A λ ) ia , jb = ( ǫ a − ǫ i ) δ i j δ ab + λ (cid:2) h ib | ˆ w ee | a j i − ξ h ib | ˆ w ee | ja i (cid:3) , (B8a)( B λ ) ia , jb = λ (cid:2) h ab | ˆ w ee | i j i − ξ h ab | ˆ w ee | ji i (cid:3) . (B8b)We need to consider the linear response non-Hermitian eigen-value equation A λ B λ B ∗ λ A ∗ λ ! X n ,λ Y n ,λ ! = ω n ,λ − ! X n ,λ Y n ,λ ! , (B9)whose solutions come in pairs: positive excitation ener-gies ω n ,λ with eigenvectors (cid:0) X n ,λ , Y n ,λ (cid:1) , and opposite (de-)excitation energies − ω n ,λ with eigenvectors (cid:16) Y ∗ n ,λ , X ∗ n ,λ (cid:17) .Choosing the normalization of the eigenvectors so that X † n ,λ X m ,λ − Y † n ,λ Y m ,λ = δ nm , the supermatrix (cid:5) λ ( ω ) can beexpressed as the following spectral representation (where thesum is over eigenvectors with positive excitation energies) (cid:5) λ ( ω ) = X n " ω − ω n ,λ + i + X n ,λ Y n ,λ ! (cid:16) X † n ,λ Y † n ,λ (cid:17) (B10) − ω + ω n ,λ − i + Y ∗ n ,λ X ∗ n ,λ ! (cid:16) Y ∗† n ,λ X ∗† n ,λ (cid:17) . (B11)3The fluctuation-dissipation theorem [Eq. (A24)] leads to thesupermatrix representation of the correlation part of the two-particle density matrix P c ,λ (using contour integration in theupper half of the complex plane) P c ,λ = − Z ∞−∞ d ω π i e i ω + [ (cid:5) λ ( ω ) − (cid:5) ( ω )] = X n Y ∗ n ,λ Y ∗† n ,λ Y ∗ n ,λ X ∗† n ,λ X ∗ n ,λ Y ∗† n ,λ X ∗ n ,λ X ∗† n ,λ − ! , (B12)the simple contribution coming from (cid:5) ( ω ) resulting fromits diagonal form [Eqs. (B3)], and the correlation energy[Eq. (A6)] has the following expression in spin-orbital basis E c = Z d λ X pq , rs h ps | ˆ w | qr i ( P c ,λ ) pq , rs = Z d λ X ia , jb X n ( h ib | ˆ w ee | a j i ( Y n ,λ ) ∗ ia ( Y n ,λ ) jb + h i j | ˆ w ee | ab i ( Y n ,λ ) ∗ ia ( X n ,λ ) jb + h ab | ˆ w ee | i j i ( X n ,λ ) ∗ ia ( Y n ,λ ) jb + h a j | ˆ w ee | ib i h ( X n ,λ ) ∗ ia ( X n ,λ ) jb − δ i j δ ab i ) , (B13)where out of the integrals h ps | ˆ w | qr i associated with thegeneral perturbation operator only the integrals of the type h ib | ˆ w ee | a j i associated with the two-electron contribution ofthe perturbation operator survive because of the occupied-virtual / occupied-virtual structure of the two-particle densitymatrix. Using now real spin orbitals, the correlation energycan be simplified to E c = Z d λ X ia , jb h ib | ˆ w ee | a j i ( P c ,λ ) ia , jb , (B14)where( P c ,λ ) ia , jb = X n (cid:0) X n ,λ + Y n ,λ (cid:1) ia (cid:0) X n ,λ + Y n ,λ (cid:1) jb − δ i j δ ab , (B15)or, in matrix form, P c ,λ = X n (cid:0) X n ,λ + Y n ,λ (cid:1) (cid:0) X n ,λ + Y n ,λ (cid:1) T − . (B16)Using the well-known fact that, if A λ + B λ and A λ − B λ arepositive definite, the non-Hermitian eigenvalue equation (B9)with real spin orbitals can be transformed into the followinghalf-size symmetric eigenvalue equation M λ Z n ,λ = ω n ,λ Z n ,λ , (B17)where M λ = ( A λ − B λ ) / ( A λ + B λ ) ( A λ − B λ ) / and witheigenvectors Z n ,λ = √ ω n ,λ ( A λ − B λ ) − / (cid:0) X n ,λ + Y n ,λ (cid:1) , andusing the spectral decomposition M − / λ = P n ω − n ,λ Z n ,λ Z T n ,λ ,the correlation two-particle density matrix P c ,λ can be ex-pressed as P c ,λ = ( A λ − B λ ) / M − / λ ( A λ − B λ ) / − . (B18)
2. Expressions for spin-restricted closed-shell calculations
For spin-restricted closed-shell calculations, the eigenvec-tors ( X n ,λ , Y n ,λ ) can be transformed into spin-singlet excita-tion / diexcitation vectors( x n ,λ ) ia = √ h(cid:0) X n ,λ (cid:1) i ↑ a ↑ + (cid:0) X n ,λ (cid:1) i ↓ a ↓ i , (B19a)( y n ,λ ) ia = √ h(cid:0) Y n ,λ (cid:1) i ↑ a ↑ + (cid:0) Y n ,λ (cid:1) i ↓ a ↓ i , (B19b)and spin-triplet excitation / diexcitation vectors( , x n ,λ ) ia = √ h(cid:0) X n ,λ (cid:1) i ↑ a ↑ − (cid:0) X n ,λ (cid:1) i ↓ a ↓ i , (B20a)( , y n ,λ ) ia = √ h(cid:0) Y n ,λ (cid:1) i ↑ a ↑ − (cid:0) Y n ,λ (cid:1) i ↓ a ↓ i , (B20b)( , − x n ,λ ) ia = (cid:0) X n ,λ (cid:1) i ↑ a ↓ , (B20c)( , − y n ,λ ) ia = (cid:0) Y n ,λ (cid:1) i ↓ a ↑ , (B20d)( , x n ,λ ) ia = (cid:0) X n ,λ (cid:1) i ↓ a ↑ , (B20e)( , y n ,λ ) ia = (cid:0) Y n ,λ (cid:1) i ↑ a ↓ , (B20f)the indices i , a , j , b referring now to spatial orbitals. With thistransformation, the linear response eigenvalue equation (B9)decouples into a singlet eigenvalue equation A λ B λ B ∗ λ A ∗ λ ! x n ,λ y n ,λ ! = ω n ,λ − ! x n ,λ y n ,λ ! , (B21)with the singlet orbital rotation Hessians (cid:16) A λ (cid:17) ia , jb = ( ǫ a − ǫ i ) δ i j δ ab + λ (cid:2) h ib | ˆ w ee | a j i − ξ h ib | ˆ w ee | ja i (cid:3) , (B22a) (cid:16) B λ (cid:17) ia , jb = λ (cid:2) h ab | ˆ w ee | i j i − ξ h ab | ˆ w ee | ji i (cid:3) , (B22b)and three identical triplet eigenvalue equations A λ B λ B ∗ λ A ∗ λ ! x n ,λ y n ,λ ! = ω n ,λ − ! x n ,λ y n ,λ ! , (B23)with the triplet orbital rotation Hessians (cid:16) A λ (cid:17) ia , jb = ( ǫ a − ǫ i ) δ i j δ ab − λξ h ib | ˆ w ee | ja i , (B24a) (cid:16) B λ (cid:17) ia , jb = − λξ h ab | ˆ w ee | ji i . (B24b)4Performing the sums over spins in the correlation energy ex-pression of Eq. (B14), one gets, for real spatial orbitals, E c = Z d λ X ia , jb h ib | ˆ w ee | a j i ( P c ,λ ) ia , jb , (B25)where remains only the contribution from the spin-singlet-adapted correlation two-particle density matrix ( P c ,λ ) ia , jb = P σ = ↑ , ↓ P σ = ↑ , ↓ ( P c ,λ ) i σ a σ , j σ b σ , which can be calculated sim-ilarly as before P c ,λ = X n (cid:16) x n ,λ + y n ,λ (cid:17) (cid:16) x n ,λ + y n ,λ (cid:17) T − = (cid:20)(cid:16) A λ − B λ (cid:17) / M − / λ (cid:16) A λ − B λ (cid:17) / − (cid:21) , (B26) where M λ = (cid:16) A λ − B λ (cid:17) / (cid:16) A λ + B λ (cid:17) (cid:16) A λ − B λ (cid:17) / . [1] P. Hohenberg and W. Kohn, Phys. Rev. , B 864 (1964).[2] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[3] J. Toulouse, F. Colonna, and A. Savin, Phys. Rev. A , 062505(2004).[4] T. Leininger, H. Stoll, H.-J. Werner, and A. Savin, Chem. Phys.Lett. , 151 (1997).[5] R. Pollet, A. Savin, T. Leininger, and H. Stoll, J. Chem. Phys. , 1250 (2002).[6] J. K. Pedersen and H. J. A. Jensen, unpublished.[7] E. Fromager, J. Toulouse, and H. J. A. Jensen, J. Chem. Phys. , 074111 (2007).[8] E. Fromager, F. R´eal, P. Wåhlin, U. Wahlgren, and H. J. A.Jensen, J. Chem. Phys. , 054107 (2009).[9] J. G. ´Angy´an, I. C. Gerber, A. Savin, and J. Toulouse, Phys.Rev. A , 012510 (2005).[10] I. C. Gerber and J. G. ´Angy´an, Chem. Phys. Lett. , 370(2005).[11] I. C. Gerber and J. G. ´Angy´an, J. Chem. Phys. , 044103(2007).[12] E. Goll, T. Leininger, F. R. Manby, A. Mitrushchenkov, H.-J. Werner, and H. Stoll, Phys. Chem. Chem. Phys. , 3353(2008).[13] B. G. Janesko and G. E. Scuseria, Phys. Chem. Chem. Phys. ,9677 (2009).[14] E. Goll, H.-J. Werner, and H. Stoll, Phys. Chem. Chem. Phys. , 3917 (2005).[15] E. Goll, H.-J. Werner, H. Stoll, T. Leininger, P. Gori-Giorgi, andA. Savin, Chem. Phys. , 276 (2006).[16] E. Goll, H. Stoll, C. Thierfelder, and P. Schwerdtfeger, Phys.Rev. A , 032507 (2007).[17] E. Goll, H.-J. Werner, and H. Stoll, Chem. Phys. , 257(2008).[18] E. Goll, M. Ernst, F. Moegle-Hofacker, and H. Stoll, J. Chem.Phys. , 234112 (2009).[19] E. Fromager, R. Cimiraglia, and H. J. A. Jensen, Phys. Rev. A , 024502 (2010).[20] J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. ´Angy´an,Phys. Rev. Lett. , 096404 (2009).[21] B. G. Janesko, T. M. Henderson, and G. E. Scuseria, J. Chem.Phys. , 081105 (2009).[22] B. G. Janesko, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. , 034110 (2009).[23] B. G. Janesko and G. E. Scuseria, J. Chem. Phys. , 154106(2009).[24] J. Paier, B. G. Janesko, T. M. Henderson, G. E. Scuseria,A. Gr¨uneis, and G. Kresse, J. Chem. Phys. , 094103 (2010).[25] Z. Yan, J. P. Perdew, and S. Kurth, Phys. Rev. B. , 16430(2000).[26] F. Furche, Phys. Rev. B , 195120 (2001).[27] F. Aryasetiawan, T. Miyake, and K. Terakura, Phys. Rev. Lett. , 166401 (2002).[28] T. Miyake, F. Aryasetiawan, T. Kotani, M. v. Schilfgaarde,M. Usuda, and K. Terakura, Phys. Rev. B. , 245103 (2002).[29] M. Fuchs and X. Gonze, Phys. Rev. B , 235109 (2002).[30] Y. M. Niquet and X. Gonze, Phys. Rev. B , 245115 (2004).[31] M. Fuchs, Y. M. Niquet, X. Gonze, and K. Burke, J. Chem.Phys. , 094116 (2005).[32] F. Furche and T. V. Voorhis, J. Chem. Phys. , 164106 (2005).[33] N. E. Dahlen, R. van Leeuwen, and U. von Barth, Phys. Rev. A , 012511 (2006).[34] A. Marini, P. Garc´ıa-Gonz´alez, and A. Rubio, Phys. Rev. Lett. , 136404 (2006).[35] H. Jiang and E. Engel, J. Chem. Phys. , 184108 (2007).[36] J. Harl and G. Kresse, Phys. Rev. B. , 045136 (2008).[37] F. Furche, J. Chem. Phys. , 114105 (2008).[38] G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, J. Chem.Phys. , 231101 (2008).[39] X. Ren, P. Rinke, and M. Sche ffl er, Phys. Rev. B , 045402(2009).[40] D. Lu, Y. Li, D. Rocca, and G. Galli, Phys. Rev. Lett. ,206411 (2009).[41] J. Harl and G. Kresse, Phys. Rev. Lett. , 056401 (2009).[42] H.-V. Nguyen and S. de Gironcoli, Phys. Rev. B. , 205114(2009).[43] H.-V. Nguyen and G. Galli, J. Chem. Phys. , 044109 (2010).[44] A. Gr¨uneis, M. Marsman, J. Harl, L. Schimka, and G. Kresse,J. Chem. Phys. , 154115 (2009).[45] M. Hellgren and U. von Barth, J. Chem. Phys. , 044101(2010).[46] J. Harl, L. Schimka, and G. Kresse, Phys. Rev. B , 115126(2010).[47] S. Ismail-Beigi, Phys. Rev. B , 195126 (2010). [48] A. Ruzsinszky, J. P. Perdew, and G. I. Csonka, J. Chem. TheoryComput. , 127 (2010).[49] J. Toulouse, A. Savin, and H.-J. Flad, Int. J. Quantum Chem. , 1047 (2004).[50] J. Toulouse, F. Colonna, and A. Savin, J. Chem. Phys. ,014110 (2005).[51] S. Paziani, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet, Phys.Rev. B , 155111 (2006).[52] The short-range self-energy correction ∆Σ sr λ is wrongly missingin Eq. (11) of Ref. 20. However, in practice, this term vanishesin the RPA or RPAx approximation so that the results of Ref. 20are correct.[53] A. D. McLachlan and M. A. Ball, Rev. Mod. Phys. , 844(1964).[54] A. Szabo and N. S. Ostlund, J. Chem. Phys. , 4351 (1977).[55] H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby, M. Sch¨utz,et al., Molpro, version 2008.2, a package of ab initio programs , 100(2005).[57] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. ,3865 (1996).[58] In the context of density-functional theory RPA is usually de-rived from the Kohn-Sham reference, while in the context ofmany-body perturbation theory (see appendices) RPA is usu-ally derived from the HF reference. Therefore, both PBE + RPAand HF + RPA are theoretically justified.[59] T. H. Dunning, J. Chem. Phys. , 1007 (1989).[60] D. Woon and T. Dunning, J. Chem. Phys. , 1358 (1993).[61] D. Woon and T. Dunning, J. Chem. Phys. , 2975 (1994).[62] D. Feller, J. Comput. Chem. , 1571 (1996). [63] A. Wilson, D. Woon, K. Peterson, and T.H.Dunning, J. Chem.Phys. , 7667 (1999).[64] J. Koput and K. A. Peterson, J. Phys. Chem. A , 9595(2002).[65] K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Guru-moorthi, J. Chase, J. Li, and T. L. Windus, J. Chem. Inf. Model. , 1045 (2007).[66] K. T. Tang and J. P. Toennies, J. Chem. Phys. , 4976 (2003).[67] I. Røeggen and L. Veseth, Int. J. Quantum Chem. , 201(2005).[68] W. J. Balfour and A. E. Douglas, Can. J. Phys. , 901 (1970).[69] O. Allard, A. Pashov, H. Kn¨ockel, and E. Tiemann, Phys. Rev.A , 042503 (2002).[70] S. G. Porsev and A. Derevianko, Phys. Rev. A , 020701(R)(2002).[71] P. Reinhardt, J. Toulouse, J. G. ´Angy´an, and A. Savin, unpub-lished.[72] W. Zhu, J. Toulouse, A. Savin, and J. G. ´Angy´an, J. Chem.Phys. , 244108 (2010).[73] G. Strinati, Rivista del Nuovo Cimento , 1 (1988).[74] E. K. U. Gross, E. Runge, and O. Heinonen, Many-particle the-ory (Verlag Adam Hilger, Bristol, 1991).[75] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. , 601(2002).[76] F. Bruneval, Ph.D. thesis, Ecole Polytechnique (2005).[77] The inverse of a 4-point function χ (1 ,
2; 1 ′ , ′ ) is defined ac-cording to R d ′ d ′ χ (1 ,
2; 1 ′ , ′ ) χ − (2 ′ , ′ ; 4 , = δ (1 , δ (2 , Methods of Molecular Quantum Mechanics. Sec-ond Edition (Academic Press, London, 1992).[79] F. Furche, J. Chem. Phys.114