János G Ángyán

Iann C Gerber

Andreas Savin

van der Waals forces in density functional theory: perturbational long-range electron interaction corrections

Keywords: numbers: 31.15.Ew, 31.15.Ar, 31.15.Md, 34.20.-b, 71.15.Mb

Long-range exchange and correlation effects, responsible for the failure of currently used approximate density functionals in describing van der Waals forces, are taken into account explicitly after a separation of the electronelectron interaction in the Hamiltonian into short-and long-range components. We propose a "range-separated hybrid" functional based on a local density approximation for the short-range exchange-correlation energy, combined with a long-range exact exchange energy. Long-range correlation effects are added by a second-order perturbational treatment. The resulting scheme is general and is particularly well-adapted to describe van der Waals complexes, like rare gas dimers.

I. INTRODUCTION

Van der Waals (dispersion) interactions are universal attractive forces due to long-range correlation of electrons between weakly-or non-overlapping electron groups [1]. They play an important role in the cohesive energy of practically all kinds of materials: intermolecular complexes, extended systems, like molecular crystals, liquids or biological macromolecules. Although, in principle, density functional theory (DFT) [2] within the Kohn-Sham (KS) scheme [3] is able to provide the exact ground state energy of an electronic system, present approximate density functionals are inappropriate to describe long-range electron correlation and consequently fail for van der Waals interactions, manifested by their incapacity of reproducing the correct R -6 asymptotic behaviour of the intermolecular potential [4].

Several propositions have been published recently to add the missing long-range correlation contribution or to use asymptotically correct correlation energy expressions in DFT. Most of these methods require a partitioning of the system into interacting parts and are valid only for large separations [5][6][7]. Seamless dispersion energy functionals [8], that are valid for the whole range of possible intermolecular separations have also been proposed [9]. A general problem in schemes that use an additive correction to standard functionals is the double counting of a part of the correlation effects that are already present in the original functional.

Moreover, it is not enough to add missing correlation effects to traditional density functionals. Many of the present approximate functionals, like the local density approximation (LDA) which is well-known for its notorious overbinding tendency, and also many popular generalized gradient approximations (GGA), already predict a more-or-less pronounced bound state for simple van der Waals complexes, like rare gas dimers [10]. As it has been pointed out by Harris twenty years ago [11], this behaviour is related to the erroneous distance dependence of approximate exchange functionals. In effect, in self-interaction corrected calculations the minimum on the potential curve disappears [5]. Therefore, in order to describe correctly both the minimum and the asymptotic region of van der Waals potential energy surfaces it is mandatory to remove the unphysical bonding by appropriately correcting the exchange functional.

Here, we propose a scheme based on a long-range/shortrange decomposition of the electron interaction which meets the above requirements and remedies the description of van der Waals forces in the framework of a first-principles approach, which takes into account simultaneously long-range correlation and exchange effects, avoids double counting and is size-extensive.

Our scheme is based on the hypothesis that for the description of van der Waals (London) dispersion forces one should improve the representation of long-range electron interaction (exchange and correlation) effects. At a first level of approximation, we treat the long-range exchange energy explicitly while maintaining a density functional approximation for the short-range exchange-correlation energy. This step defines a "range-separated hybrid" (RSH) scheme, which is corrected in a second step for the long-range correlation effects by a second-order perturbation theory, leading to size extensive Møller-Plesset (MP2)-like correction. This method will be referred to as RSH+MP2.

The idea of a long-range/short-range decomposition of the electron interaction is not new (see, e.g., Refs. 12-16). In the context of DFT, this approach has been used to construct multi-determinantal extensions of the KS scheme [17][START_REF] Savin | Beyond the Kohn-Sham determinant[END_REF][START_REF] Savin | On degeneracy, near degeneracy and density functional theory[END_REF][START_REF] Leininger | [END_REF][21][22][23][24]. A density functional scheme with correct asymptotic behaviour has been proposed along these lines very recently by Baer and Neuhauser [25] and the correct 1/r asymptotic behaviour of the long-range exact exchange has been also exploited in timedependent DFT calculations of polarizabilities [26,27], constituting the major motivation of the recent development of "Coulomb attenuated" hybrid functionals [28]. Heyd, Scuseria and Ernzerhof applied an inverse range-separation in order to get rid of the convergence problems of the exact exchange in solid-state calculations [29,30]. Their HSE03 functional is a generalization of the PBE0 hybrid functional [31] where the long-range portion of the exact exchange is replaced by the long-range component of the PBE exchange functional [32].

In the context of the calculation of van der Waals energies, the idea of separating the electron interaction operator to short-and long-range components has already been explored by the work of Kohn, Meier and Makarov, who applied the adiabatic connection -fluctuation-dissipation approach for long-range electron interactions [33], leading to an asymptotically correct expression of the dispersion forces. It has also been shown [34] that the artificial minimum of the rare gas dimer potential curves can be removed by an exact treatment of the long-range exchange.

The second order perturbational treatment of the full Coulomb interaction has already been used by several authors for the van der Waals problem [5,35,36], and it was shown that the resulting asymptotic potential has the qualitatively correct 1/R 6 form. As shown very recently, quantitatively reliable asymptotic form of the potential energy curve can be expected from adiabatic connection -fluctuation-dissipation theory calculations [37].

The general theoretical framework is outlined in Section II, describing the RSH scheme and the second-order perturbational treatment of long-range correlation effects. As described in Section III, our approach has been tested on rare gas dimers. These systems are typical van der Waals complexes, where the attractive interactions are exclusively due to London dispersion forces. They constitute a stringent test of the method, since the potential curves have very shallow minima of the order of about 100 µH.

Unless otherwise stated, atomic units is assumed throughout this work.

II. THEORY A. Multi-determinantal extension of the Kohn-Sham scheme

We first recall the principle of the multi-determinantal extension of the KS scheme based on a long-range/short-range decomposition (see, e.g., Ref. 24

and references therein).

The starting point is the decomposition the Coulomb electron-electron interaction w ee (r) = 1/r as w ee (r) = w lr,µ ee (r) + w sr,µ ee (r),

where w lr,µ ee (r) = erf(µr)/r is a long-range interaction and w sr,µ ee (r) is the complement short-range interaction. This decomposition is controlled by a single parameter µ. For µ = 0, the long-range interaction vanishes, w lr,µ=0 ee (r) = 0, and the short-range interaction reduces to the Coulomb interaction, w sr,µ=0 ee (r) = w ee (r). Symmetrically, for µ → ∞, the short-range interaction vanishes, w sr,µ→∞ ee (r) = 0, and the long-range interaction reduces to the Coulomb interaction, w lr,µ→∞ ee (r) = w ee (r). Physically, 1/µ represents the distance at which the separation is made.

The Coulombic universal density functional F [n] = min Ψ→n Ψ| T + Ŵee |Ψ [START_REF] Levy | Proc. Natl. Acad. Sci[END_REF], where T is the kinetic energy operator, Ŵee = (1/2) dr 1 dr 2 w ee (r 12 )n 2 (r 1 , r 2 ) is the Coulomb electron-electron interaction operator expressed with the pair-density operator n2 (r 1 , r 2 ), is then decomposed as

F [n] = F lr,µ [n] + E sr,µ Hxc [n], (2)
where

F
E = min n→N F lr,µ [n] + E sr,µ Hxc [n] + drv ne (r)n(r) = min Ψ→N Ψ| T + Ŵ lr,µ ee |Ψ + drv ne (r)n Ψ (r) + E sr,µ Hxc [n Ψ ] , (3)
where the last search is carried out over all N -electron normalized (multi-determinantal) wave functions Ψ. In Eq. ( 3), n Ψ (r) is the density coming from the wave function Ψ, i.e. n Ψ (r) = Ψ|n(r)|Ψ where n(r) is the density operator.

The minimizing wave function Ψ µ in Eq. ( 3) is given by the corresponding Euler-Lagrange equation

T + Ŵ lr,µ ee + Vne + V sr,µ Hxc [n Ψ µ ] |Ψ µ = E µ |Ψ µ , (4)
where Vne = drv ne (r)n(r), V sr,µ Hxc [n] = drv sr,µ Hxc (r)n(r) with the short-range Hartree-exchange-correlation potential v sr,µ Hxc (r) = δE sr,µ Hxc [n]/δn(r), and E µ is the Lagrange multiplier associated to the constraint of the normalization of the wave function. Eq. ( 4) defines a long-range interacting effective Hamiltonian Ĥµ = T + Ŵ lr,µ ee + Vne + V sr,µ Hxc [n Ψ µ ] that must be solved iteratively for its multi-determinantal ground-state wave function Ψ µ which gives, in principle, the exact physical ground-state density n(r) = n Ψ µ (r) = Ψ µ |n(r)|Ψ µ , independently of µ. Finally, the exact ground-state energy expression is thus

E = Ψ µ | T + Ŵ lr,µ ee + Vne |Ψ µ + E sr,µ Hxc [n Ψ µ ]. (5)
This exact formalism enables to combine a long-range wave function calculation with a short-range density functional. In the special case of µ = 0, the KS scheme is recovered, while the limit µ → ∞ corresponds to the usual wave function formulation of the electronic problem.

A short-range LDA [START_REF] Toulouse | [END_REF] and other beyond-LDA [24,40] approximations have been constructed to successfully describe the functional E sr,µ xc [n]. In previous applications of the method, the long-range part of the calculation has been handled by configuration interaction [21] or multi-configurational self-consistent field (MCSCF) [23] methods. We propose in this work to use instead perturbation theory.

B. Range-separated hybrid

At a first level of approximation, we introduce the RSH scheme by restricting the search in Eq. ( 3) to N -electron normalized one-determinant wave functions Φ

E µ,RSH = min Φ→N Φ| T + Ŵ lr,µ ee |Φ + drv ne (r)n Φ (r) + E sr,µ Hxc [n Φ ] . (6)
The associated minimizing one-determinant wave function Φ µ satisfies the Euler-Lagrange equation

T + Vne + V lr,µ Hx,HF [Φ µ ] + V sr,µ Hxc [n Φ µ ] |Φ µ = E µ 0 |Φ µ , (7)
[Φ] = -(1/2) dr 1 dr 2 w lr,µ ee (r 12 ) Φ|n 1 (r 2 , r 1 )|Φ n1 (r 1 , r 2 )
, where n1 (r 1 , r 2 ) is the first-order density matrix operator. Eq. ( 7) defines the RSH non-interacting effective Hamiltonian

Ĥµ 0 = T + Vne + V lr,µ Hx,HF [Φ µ ] + V sr,µ Hxc [n Φ µ ]
that must be solved iteratively for its one-determinant ground-state wave function Φ µ . Of course, Φ µ does not give the exact physical density:

n Φ µ = n.
The RSH energy expression is finally

E µ,RSH = Φ µ | T + Vne |Φ µ + E lr,µ Hx,HF [Φ µ ] + E sr,µ Hxc [n Φ µ ], (8)
where E lr,µ Hx,HF [Φ] = Φ| Ŵ lr,µ ee |Φ is the HF-like long-range Hartree-exchange energy. Eq. ( 8) defines a single-parameter hybrid scheme combining a long-range HF calculation with a short-range density functional. The case µ = 0 still corresponds to the KS scheme while the method reduces now to a standard HF calculation in the limit µ → ∞.

We note that an equivalent to the RSH scheme has been investigated recently by Pedersen and Jensen [23] as a special case of the combination of a long-range MCSCF calculation with a short-range density functional.

C. Long-range correlation corrections by perturbation theory

We develop now a long-range perturbation theory, using the RSH determinant Φ µ as the reference. To do so, we introduce the following energy expression with a formal coupling constant λ E µ,λ = min

Ψ→N Ψ| T + Vne + V lr,µ Hx,HF [Φ µ ] + λ Ŵlr,µ |Ψ +E sr,µ Hxc [n Ψ ] , (9)
where the search is carried out over all N -electron normalized (multi-determinantal) wave functions Ψ and Ŵlr,µ is the longrange fluctuation potential operator

Ŵlr,µ = Ŵ lr,µ ee -V lr,µ Hx,HF [Φ µ ]. (10)
The minimizing wave function Ψ µ,λ in Eq. ( 9) is given by the Euler-Lagrange equation

T + Vne + V lr,µ Hx,HF [Φ µ ] + λ Ŵlr,µ + V sr,µ Hxc [n Ψ µ,λ ] |Ψ µ,λ = E µ,λ |Ψ µ,λ , (11)
where E µ,λ is the Lagrange multiplier associated to the normalization constraint. For λ = 1, the physical energy is recovered, E = E µ,λ=1 , in principle independently of µ, and Eq. ( 11) reduces to Eq. ( 4):

Ψ µ,λ=1 = Ψ µ , E µ,λ=1 = E µ .
For λ = 0, Eq. ( 11) reduces to the RSH effective Schrödinger equation of Eq. ( 7):

Ψ µ,λ=0 = Φ µ , E µ,λ=0 = E µ 0 . We expand E µ,λ in powers of λ, E µ,λ = ∞ k=0 E µ,(k) λ k ,
and apply the general results of the non-linear Rayleigh-Schrödinger perturbation theory [41][42][43] outlined in the Appendix. It is easy to verify that the sum of zeroth-and firstorder energy contributions gives back the RSH total energy

E µ,(0) + E µ,(1) = E µ,RSH . (12)
The second-order correction can be written as

E µ,(2) = -Φ µ | Ŵlr,µ 1 + Rµ 0 Ĝµ 0 -1 Rµ 0 Ŵlr,µ |Φ µ , (13)
where Rµ 0 is the reduced resolvent

Rµ 0 = I |Φ µ I Φ µ I | E µ 0,I -E µ 0 , (14)
in terms of the excited eigenfunctions Φ µ I and eigenvalues E µ 0,I of the RSH effective Hamiltonian Ĥµ 0 , and Ĝµ 0 is a short-range screening operator Ĝsr,µ

0 = 2 drdr ′ n(r)|Φ µ f sr,µ Hxc [n Φ µ ](r, r ′ ) Φ µ |n(r ′ ), (15)
with the short-range Hartree-exchange-correlation kernel f sr,µ Hxc [n](r, r ′ ) = δ 2 E sr,µ Hxc [n]/δn(r)δn(r ′ ). Let insert the spectral resolution of Eq. ( 14) in Eq. ( 13). Since Ŵlr is a two-electron operator only singly and doubly excited determinants, Φ µ i→a and Φ µ ij→ab where i, j refer to occupied spin-orbitals and a, b to virtual spin-orbitals of Φ µ , can a priori contribute to E µ, (2) . Actually, singly excited determinants gives vanishing matrix elements, Φ µ i→a | Ŵlr |Φ µ = 0, since it can be easily verified that

Φ µ i→a | Ŵ lr ee |Φ µ = Φ µ i→a | V lr,µ
Hx,HF [Φ µ ]|Φ µ , as in standard HF theory. Consequently, the product Rµ 0 Ĝµ 0 in Eq. ( 13) involves vanishing matrix elements, Φ µ |n(r)|Φ µ ij→ab = 0, i.e. the non-linear terms are zero with the present choice of the perturbation operator Ŵlr . The second-order energy correction is thus

E µ,(2) = -Φ µ | Ŵlr,µ Rµ 0 Ŵlr,µ |Φ µ = i<j a<b | Φ µ ij→ab | Ŵ lr,µ ee |Φ µ | 2 E µ 0 -E µ 0,ij→ab = i<j a<b | φ µ i φ µ j | ŵlr,µ ee |φ µ a φ µ b -φ µ i φ µ j | ŵlr,µ ee |φ µ b φ µ a | 2 ε µ i + ε µ j -ε µ a -ε µ b , (16)
where φ µ k is a spin-orbital of Φ µ and ε µ k is its associated eigenvalue, φ µ i φ µ j | ŵlr,µ ee |φ µ a φ µ b are the two-electron integrals associated to the long-range interaction w lr,µ ee (r 12 ), and we recall that the indexes i, j refer to occupied spin-orbitals and a, b to virtual spin-orbitals. Eq. ( 16) is fully analogous to the conventional MP2 energy correction. The total RSH+MP2 energy is E µ,RSH+MP2 = E µ,RSH + E µ, (2) . ¿From a practical point of view, once the RSH orbitals and one-electron eigenvalues are available, any standard MP2 implementation can be used, provided that the long-range electron repulsion integrals corresponding to the RSH orbitals are plugged in. Due to the long-range nature of these integrals one can take advantage of efficient modern algorithms, like the local MP2 [44], multipolar integral approximations, which have particularly favorable convergence properties for long-range part of the split Coulomb interaction [45], or the resolution of identity approach [46]. It means that in appropriate implementations the extra cost of the MP2 corrections can be made negligible for large systems with respect to the resolution of the self-consistent RSH equations, similar to a usual KS calculations with a hybrid functional. Solid state applications for semi-conductors can also be envisaged on Wannier orbital-based MP2 implementations [47].

III. RESULTS AND CONCLUSIONS

The above described RSH+MP2 approach has been applied to rare gas dimers, using a LDA-based short-range exchangecorrelation functional with a range-separation parameter of µ = 0.5. This latter value corresponds to the smallest mean average error of the atomization energies calculated by the RSH scheme for the G2-1 set (a subset of 55 molecules of the G3 set [48,49]) of small molecules [50]. This value is in agreement with the intuitive picture predicting that 1/µ should be close to the physical dimensions of a valence electron pair. The interaction energies were calculated with a modified version of the MOLPRO package [START_REF] Werner | MOLPRO, a package of ab initio programs[END_REF]. The basis set superposition error (BSSE) has been removed by the counterpoise method.

The results are presented as reduced potentials, U * (r * ) = U (r * • d m )/ε m , where the reduced variables U * = U/ε m and r * = r/d m , are defined with respect to the equilibrium distance d m and the well-depth ε m of accurate "experimental" potential curves [52] (cf. Table I). The calculated potentials are characterized by the hard core radius, σ * defined by U * (σ * ) = 0 (experimentally σ * ≈ 0.89), the reduced well depth, U * m , and the equilibrium distance, r * m , (experimentally, by construction, U * m = -1 and r * m = 1). The minimum region is also characterized by the harmonic vibrational frequencies, ω, related to the second derivative of the potential at the minimum.

The long-range behaviour of the potential energy curves can be appreciated from the C 6 coefficients. Experimental C 6 coefficients are usually obtained from optical data (dipole oscillator strength distributions) [53] and characterize the purely dipolar contribution to the long-range interaction energy. Since we had no access to such a decomposition of the interaction energy, we have determined an effective C 6 coefficient by a logarithmic fit of the interaction energies between 30 and 60 Bohrs. This quantity, which includes higher order multipolar effects too, is presented in the form of a reduced variable, C * 6 = C 6 /C fit 6 . Here C fit 6 has been obtained from an analogous fit to the points of the reference potential reported in Table I. For the sake of comparison, the experimental C exp 6 (purely dipolar) values are also reported.

The RSH and RSH+MP2 potential curves, as well as the HF, the standard MP2 and the coupled-cluster CCSD(T) ones, calculated with the aug-cc-pVTZ basis set are represented for the four dimers in Figure 1, and compared to the experimental curves. Note that the reduced representations of the ex- perimental potentials of different rare gas dimers are practically indistinguishable. The calculated RSH potentials are always repulsive, like the HF ones. The RSH+MP2 potentials are slightly too repulsive at short interatomic distances, as reflected by the values of the hard core radii, systematically higher than the experiment (around 0.89). The RSH+MP2 and CCSD(T) curves are almost the same for Ne 2 with a well depth of around U * m = 0.6, while the RSH+MP2 minima of the Ar 2 and Kr 2 systems are even better (U * m > 0.9) than the CCSD(T) ones (U * m ≈ 0.7). The position of the minimum is predicted within 1-4% in the RSH+MP2 approximation. The 6-8% deviation found for the He 2 RSH+MP2 minimum can be explained by an exaggerated repulsion, reflected by the highest σ * found in this case. In comparison with the usual MP2 potential curves, the RSH+MP2 follow similar trends, being systematically more stabilizing and closer to the experimental curve.

The main quantitative features of the RSH+MP2 potentials obtained by the aug-cc-pVTZ and aug-cc-pV5Z basis sets [54][55][56][57] are summarized in Table II and compared to the results of standard MP2 and CCSD(T) supermolecule calculations with the same basis sets. The basis set has a nonnegligible effect on the calculated parameters of the potential curves, which converge systematically towards the experimental values for all the properties. For He 2 the double augmented d-aug-cc-pV5Z basis set results are also included, representing a further improvement of the well depth, but having practically no effect on the equilibrium distance.

The basis set superposition error of the equilibrium distances and of the interaction energies are reported in Table III, as the difference in the parameters of the BSSE-contaminated and BSSE-free reduced potential energy curves. The BSSE corrections on the bond lengths and on the interaction energies are always negative, i.e. the BSSE-contaminated distances are too short and the energies are too low. In some case, like the Ne 2 dimer with aug-cc-pVTZ basis, the binding energy correction may attain 55 or 67% of the well depth at the MP2 and CCSD(T) level of approximation. The corresponding RSH+MP2 BSSE effect is considerably smaller, but it is still 34%. The BSSE effect on the bond lengths are much less spectacular, but still more pronounced in the MP2 and CCSD(T) methods than in the RSH+MP2 approach. As a general trend we can conclude that the RSH+MP2 has usually less than the half of the MP2 or CCSD(T) basis set superposition errors. This is a considerable advantage for an efficient and reliable exploration of potential energy surfaces, especially when the lack of well-defined subsystems make impossible to perform a counterpoise correction.

Effective C * 6 coefficients obtained from the RSH+MP2 approach agree with the experiment within 5% for He 2 and Ne 2 , and are overestimated by 15-20% for Ar 2 and Kr 2 . It means that the asymptotic behaviour of the RSH+MP2 potential curves is reasonable. We recall that the exact C 6 coefficient is given by the Casimir-Polder relation [1]

C 6 = 3 π ∞ 0 dωα 1 (iω)α 2 (iω) (17)
where α 1 (iω) and α 2 (iω) are the exact dynamical polarizabilities of the monomers. It is known that the asymptotic form of the MP2 energy expression corresponds to an uncoupled HF-type, non-interacting approximation of the monomer polarizabilities. This means that MP2 calculations do not reproduce the exact C * 6 coefficients: usually they tend to overestimate them. For instance, in the case of the benzene dimer, this overestimation in the complete basis limit may reach a factor of 2; for less polarizable systems the situation is less critical. An analogous behaviour is expected for RSH+MP2. Note however, that in this case one-electron excitations are obtained from the self-consistent RSH one-electron states, which include, in addition to the long-range exact exchange, short-range exchange-correlation effects too. A more reliable approximation can be developed on the basis of the adiabatic connection -fluctuation-dissipation approach [33,37] which would ensure, in principle, the exact asymptotic limit of the potential energy curves. The development of a rangeseparated version of this method is under progress.

In conclusion, the RSH+MP2 approach provides an efficient DFT-based description of weak intermolecular complexes bound by dispersion forces. Even in its simplest, LDAbased implementation, it represents a huge improvement over KS calculations, which lead to unreliable potential curves in the minimum region with a qualitatively wrong asymptotic behaviour. Range-separated extensions of other density functionals, like the gradient-corrected PBE functional, are in progress. By removing systematic errors of currently used approximate DFT functionals and introducing corrections which grasp the essential physics of van der Waals interactions, the RSH+MP2 approach extends the applicability of density functional calculations to weak intermolecular forces. Further tests should decide whether this method is generally applicable to the important domains of the physisorption, or cohesion in molecular crystals and in layered semi-conductors.

Appendix A: Nonlinear Rayleigh-Schrödinger perturbation theory

Let us consider the following general total energy expression, involving a Hamiltonian Ĥ(0) , a perturbation operator Ŵ and a density functional F [n],

E λ = min Ψ→N Ψ| Ĥ(0) + λ Ŵ |Ψ + F [n Ψ ] , (A1)
where the search is carried out over all N -electron normalized wave functions Ψ, Ψ|Ψ = 1, and n Ψ is the density coming from Ψ, n Ψ (r) = Ψ|n(r)|Ψ , where n(r) is the density operator. In Eq. (A1), λ is a formal coupling constant; we are ultimately interested in the case λ = 1. The minimizing wave function Ψ λ satisfies the Euler-Lagrange equation

Ĥ(0) + λ Ŵ + Ωλ |Ψ λ = E λ |Ψ λ , (A2)
where the eigenvalue E λ comes from the normalization condition and Ωλ is a potential operator coming from the variation of F [n], non linear in λ,

Ωλ = dr δF [n λ ] δn(r) n(r), (A3)
where n λ is the density coming from Ψ λ , n λ (r) = Ψ λ |n(r)|Ψ λ . Starting from the reference λ = 0, we develop a perturbation theory in λ. We introduce the intermediate normalized wave function Ψλ

| Ψλ = |Ψ λ Ψ λ=0 |Ψ λ , (A4)
and expand Ψλ , n λ , Ωλ and E λ in powers of λ:

Ψλ = ∞ k=0 Ψ(k) λ k , n λ = ∞ k=0 n (k) λ k , Ωλ = ∞ k=0 Ω(k) λ k and E λ = ∞ k=0 E (k) λ k . The coefficients n (k) are obtained from the expansion of Ψλ through n λ (r) = Ψλ |n(r)| Ψλ Ψλ | Ψλ , (A5)
and the coefficients Ω(k) are found from the expansion of n λ , after expanding Ωλ around n (0) ,

Ωλ = dr δF [n (0) ] δn(r) n(r) + drdr ′ δ 2 F [n (0) ] δn(r)δn(r ′ ) ∆n λ (r ′ )n(r) + • • • . (A6)
where ∆n λ = n λ -n (0) . The zeroth-order equation is

Ĥ(0) + Ω(0) | Ψ(0) = E (0) | Ψ(0) , (A7)
and of course Ψ(0) = Ψ λ=0 . For the general order k ≥ 1,

Ĥ(0) + Ω(0) -E (0) | Ψ(k) + Ŵ | Ψ(k-1) + k i=1 Ω(i) | Ψ(k-i) = k i=0 E (i) | Ψ(k-i) . (A8)
The corresponding eigenvalue correction of order k is (A14)

E (k) = Ψ(0) | Ŵ | Ψ(k-1) + k i=1 Ψ(0) | Ω(i) | Ψ(k-i) , (
The zeroth-order total energy is simply ) ] -dr δF [n (0) ] δn(r) n (0) (r), (A15)

E (0) = E (0) + F [n (0
The general correction of order k ≥ 1 writes

E (k) = Ψ(0) | Ŵ | Ψ(k-1) + ∆ (k) (A16)
where ∆ (k) is

∆ (k) = k i=1 Ψ(0) | Ω(i) | Ψ(k-i) + D (k) . ( A17
)
At first order, it can be verified that the nonlinearity term of the eigenvalue and the double counting correction cancel each other, i.e. ∆ (1) = 0, and we obtain the conventional first-order energy correction

E (1) = Ψ(0) | Ŵ | Ψ(0) . ( A18
)
At second order, the situation is analogous, i.e. ∆ (2) = 0, and again the conventional form of the energy correction is retrieved

E (2) = Ψ(0) | Ŵ | Ψ(1) . ( A19
)
The nonlinearity effects are "hidden" in the first-order wave function correction, which can be obtained from the selfconsistent equation:

| Ψ(1) = -R0 Ŵ | Ψ(0) -R0 Ω(1) | Ψ(0) (A20)
Since the first-order potential operator is, for real wave functions, Ω (A23) The final expression of the second-order energy correction can be written as the series

E (2) = -Ψ(0) | Ŵ 1 + R0 Ĝ0 -1 R0 Ŵ | Ψ(0) = -Ψ(0) | Ŵ R0 Ŵ | Ψ(0) + Ψ(0) | Ŵ R0 Ĝ0 R0 Ŵ | Ψ(0) -• • • . (A24)
Further details and higher-order expressions will be given in a forthcoming publication.

  [n] = T s [n], and the short-range functional to the usual Hartree-exchange-correlation functional, E sr,µ=0 Hxc [n] = E Hxc [n]. In the limit µ → ∞, the longrange functional reduces to the Coulombic universal functional, F lr,µ→∞ [n] = F [n], and the short-range functional vanishes, E sr,µ→∞ Hxc [n] = 0. The exact ground-state energy of a N -electron system in an external nuclei-electron potential v ne (r), E = min n→N F [n] + drv ne (r)n(r) where the search is over all N -representable densities, can be re-expressed using the long-range/short-range decomposition of F [n]

8 FIG. 1 :

 81 FIG.1: Reduced HF (dotted repulsive), RSH (dashed repulsive), MP2 (long dashed), CCSD(T) (dashed), RSH+MP2 (dotted) and Tang-Toennis reference (full) potential curves for He2, Ne2, Ar2 and Kr2 dimers.

I′ δ 2 F

 2 A9)containing, besides the usual first term, a "non-linearity" term as well. Introducing the reduced resolvent, R0 , are the excited eigenfunctions and eigenvalues of Ĥ(0) , the wave function correction of order k writes| Ψ(k) = -R0 Ŵ | Ψ(k-1) -R0 Ω(k) | Ψ(0) ) -E (i) | Ψ(k-i) . (A11)The total energy can be re-expressed in terms of the eigenvalue E λ and the "double counting correction" D λE λ = E λ + D λ ,(A12)whereD λ = F [n λ ] -dr δF [n λ ] δn(r) n λ (r).(A13)We expand E λ and D λ in powers of λ:E λ = ∞ k=0 E (k) λ k and D λ = ∞ k=0 D (k) λ k . The coefficients D (k) are found from the expansion of n λ , after expanding D λ around n (0) , D λ = F [n (0) ] + dr δF [n (0) ] δn(r) ∆n λ [n (0) ] δn(r)δn(r ′ ) ∆n λ (r ′ )∆n λ (r) + • • • -dr δF [n (0) ] δn(r) n λ (r) -drdr ′ δ 2 F [n (0) ] δn(r)δn(r ′ ) ∆n λ (r ′ )n λ (r) -• • • .

= 2 drdr ′ δ 2 F

 2 [n(0) ] δn(r)δn(r ′ ) Ψ(0) |n(r ′ )| Ψ(1) n(r), (A21) Eq. (A20) can be re-expressed as| Ψ(1) = -R0 Ŵ | Ψ(0) -R0 Ĝ0 | Ψ(1) , (A22) where Ĝ0 = 2 drdr ′ n(r)| Ψ(0) δ 2 F [n (0) ] δn(r)δn(r ′ ) Ψ(0) |n(r ′ ).

TABLE II :

 II Reduced parameters of the calculated MP2, CCSD(T) and RSH+MP2 (µ =0.5) potential energy curves obtained by the aug-cc-pVTZ (AVTZ), aug-cc-pV5Z (AV5Z) and d-aug-cc-pV5Z (d-AV5Z) basis sets. Reduced experimental parameters are listed in the first line for each dimer. Absolute reference values are given in TableI.

	Method	MP2		CCSD(T)	RSH+MP2
	System	r * m	U * m	r * m	U * m	r * m	U * m
	He2
	AVTZ	-0.007 -0.125	-0.008 -0.121	-0.000 -0.063
	AV5Z	-0.004 -0.048	-0.002 -0.037	-0.001 -0.015
	d-AV5Z -0.008 -0.160	-0.007 -0.113	-0.001 -0.031
	Ne2
	AVTZ	-0.049 -0.547	-0.035 -0.674	-0.031 -0.335
	AV5Z	-0.011 -0.150	-0.006 -0.148	-0.001 -0.025
	Ar2
	AVTZ	-0.020 -0.263	-0.023 -0.239	-0.007 -0.101
	AV5Z	-0.005 -0.138	-0.004 -0.103	-0.002 -0.022
	Kr2
	AVTZ	-0.013 -0.191	-0.017 -0.174	-0.007 -0.126
	AV5Z	-0.003 -0.073	-0.002 -0.049	-0.002 -0.039
	TABLE III: BSSE correction for the reduced parameters r * m and U * m .


Julien Toulouse

Simple model of the static exchange-correlation kernel of a uniform electron gas with long-range electron-electron interaction

Keywords: numbers: 71, 10, Ca,71, 45, Gm

Simple model of the static exchange-correlation kernel of a uniform electron gas with long-range electron-electron interaction

Julien Toulouse

I. INTRODUCTION

Density functional theory (DFT) 1 applied within the Kohn-Sham (KS) 2 scheme is nowadays a widely-used method for electronic calculations in condensed-matter physics and quantum chemistry. The design of better approximations to the exchange-correlation energy functional, the central quantity of the theory, constitutes an important topic of research.

Following early ideas [3][4][START_REF] Stoll | Density Functional Method in Physics[END_REF] , the possibility of improving the description of the exchange and/or correlation energy functionals by treating separately its long-range and short-range components has recently gained growing interest [START_REF] Savin | [END_REF][START_REF] Savin | Recent Advances in Density Functional Theory[END_REF][START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Leininger | [END_REF][10][11][12][13][14][15][16][17][18][19][20][21][22] . In this approach, the Coulomb electronelectron interaction is decomposed as

1 r = v µ (r) + vµ (r), (1)
where v µ (r) is a long-range interaction, vµ (r) is the complement short-range interaction and µ is a parameter controlling the range of the separation. For example, a convenient long-range interaction which has been often used is the so-called erf interaction 8-12,14-21,23

v µ erf (r) = erf(µr) r , (2)
which vanishes for µ = 0 and reduces to the Coulomb interaction for µ → ∞. The Coulombic exchangecorrelation energy functional E xc,coul [n] can be in turn decomposed as [START_REF] Savin | [END_REF][START_REF] Savin | Recent Advances in Density Functional Theory[END_REF][START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Leininger | [END_REF][10]12,20,21

E xc,coul [n] = E µ xc,erf [n] + Ēµ xc,erf [n], (3)
where E µ xc,erf [n] is the long-range component of exchangecorrelation energy, associated to the interaction v µ erf (r), and Ēµ xc,erf [n] is the complement short-range part. The long-range exchange-correlation energy, E µ xc,erf , can be efficiently calculated by configuration interaction (CI) [START_REF] Leininger | [END_REF]10 or multi-configurational self-consistent-field (MCSCF) 24 methods enabling to handle near degeneracy in atoms and molecules, and by second-order perturbation theory 25 or adiabatic connection fluctuationdissipation (ACFD) 13,26 approaches for describing van der Waals interactions. The short-range exchangecorrelation energy Ēµ xc,erf can be well approximated by simple (semi)local density functional approximations 20,21 based a uniform electron gas model with a modified electron-electron interaction [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF]27 (see also Ref. 28).

Knowledge of the properties of the uniform electron gas is always useful to gain more insight into the physics of electronic correlation of inhomogeneous systems and ultimately to construct better approximations. In particular, in standard KS-DFT, the static exchange-correlation kernel of the Coulombic uniform electron gas, f xc,coul , is of interest and several accurate parametrizations are available (see, e.g., Refs. 29-31). Indeed, f xc,coul has been recognized as a potential ingredient for density functional approximations (see, e.g., Refs. 1,2,32-37). f xc,coul can also be used in principle in ACFD approaches (see, e.g., Ref. 38 for an example using a simplified kernel) or time-dependent density functional theory (TDDFT) for density-density response calculations, although in these cases clear strategies to use f xc,coul are still under investigations. As a matter of fact, the high non-locality of f xc,coul questions its transferability to inhomogeneous systems and limits its use.

In the framework of a long-range/short-range energy decomposition in DFT, knowledge of the static exchangecorrelation kernel f µ xc,erf of a uniform electron gas with the long-range erf interaction v µ erf (r), as well as its complement analogue f µ xc,erf = f xc,coul -f µ xc,erf associated to the short-range part of the Coulomb interaction, is also of interest. In particular, f µ xc,erf can be used in approximations for the long-range exchange-correlation energy E µ xc,erf within ACDF approaches and f µ xc,erf can be used to construct density functional approximations for the short-range exchange-correlation energy Ēµ xc,erf , fol-lowing the same ideas as for the Coulombic functional E xc,coul . In comparison to f xc,coul , the non-locality character of f µ xc,erf and f µ xc,erf tends to be attenuated by the reduction of the interaction and the modified kernels may thus be more transferable to inhomogeneous systems. In principle, these kernels can also be useful in the context of TDDFT.

The main purpose of this work is therefore the determination of the static exchange-correlation kernel of the uniform electron gas with the long-range erf interaction

f µ xc,erf , the complement kernel f µ xc,erf = f xc,coul -f µ xc,erf
being directly deducible from it. A simple approximation is constructed for f µ xc,erf from its exact asymptotic behaviors which in turn require, among other things, information from the momentum distribution associated to the modified interaction. Therefore, the momentum distribution is preliminary calculated by many-body perturbation theory (MBPT), in particular using the GW approximation 39 . All the calculations are also performed with the Coulomb interaction which enables a check by comparison to other available data. The production of data of the same quality for f xc,coul and f µ xc,erf is also important to construct the complement kernel f µ xc,erf . Let's describe now the system under consideration. The uniform electron gas is described as N electrons in a box of volume Ω with a uniform neutralizing background, studied in the thermodynamic limit (i.e. N → ∞ and Ω → ∞ such that the density n = N/Ω remains constant). In second quantization, its Hamiltonian writes

Ĥ = T + V , (4)
where T is the kinetic energy operator

T = kσ ε(k)c † kσ c kσ , (5)
with ε(k) = k 2 /2, and V is the electron-electron interaction

V = 1 2Ω q =0 kσ k ′ σ ′ v(q)c † k+qσ c † k ′ -qσ ′ c k ′ σ ′ c kσ , (6)
where the constant term q = 0 has been omitted since it cancels out with the electron-background and background-background interactions. In Eq. ( 6), v(q) is the Fourier transform of an arbitrary electron-electron interaction. The standard case corresponds to the Coulomb interaction

v coul (q) = 4π q 2 , (7)
while for the long-range erf interaction the Fourier transform writes

v µ erf (q) = 4π q 2 e -q 2 /(4µ 2 ) . ( 8
)
The paper is organized as follows. In Sec. II, we discuss the calculation of the momentum distribution and associated properties. In Sec. III, the static exchangecorrelation kernel is obtained by interpolation from its limiting behaviors. In Sec. IV, as a basic example of the use of this kernel, the long-range correlation energy of the uniform electron gas is calculated. Sec. V summarizes and concludes this work. Details of derivations and analytical parametrizations can be found in Appendices A to D.

Unless otherwise stated, atomic units is assumed throughout this work.

II. MOMENTUM DISTRIBUTION

One of the most reliable calculations of the momentum distribution of the Coulombic electron gas for a wide range of densities are that of Takada and Yasuhara 40 using the effective-potential expansion method, and the quantum Monte Carlo (QMC) calculation of Ortiz and Ballone 41 . Here, we will consider the use of the more traditional methods of MBPT.

In MBPT, the momentum distribution n(k) is expressed as

n(k) = +∞ -∞ dω 2πi e iω0 + G(k, ω), (9)
where G(k, ω) is the one-particle Green function, calculated in a given approximation. Daniel and Vosko 42 have calculated the momentum distribution of the Coulombic electron gas in the RPA. The fact that the RPA probably constitutes the simplest MBPT approximation for the electron gas and becomes exact in the high-density limit (r s → 0 where r s = (4πn/3) -1/3 ) makes it valuable for comparison to more elaborate approximations. Later, additional exchange terms have been included by Lam 43 , but the improvement over the RPA remains modest. More recently, more accurate MBPT calculations of the momentum distribution of the Coulombic electron gas have been reported using the GW approximation 39 at different self-consistency levels [44][45][46][47] . In its semi-selfconsistent, GW 0 , or fully self-consistent, GW , versions this method is known to be particle-conserving 45,48,49 , i.e. fulfilling the exact condition

2 dk (2π) 3 n(k) = n, (10)
while the non-self-consistent G 0 W 0 approximation is not 50 . However, even in this case, the violation of Eq. ( 10) is numerically small (2% at most for r s ≤ 10 according to our own calculations). Moreover, imposition of the self-consistency does not necessary improve the momentum distribution; in particular the self-consistency increases the quasiparticle renormalization factor at the Fermi surface Z F [Eq. ( 11)] 45 which is suspected to be already too large at the G 0 W 0 level (see below). We therefore choose to compute the momentum distribution of the uniform electron gas with the modified erf interaction in the G 0 W 0 approximation (see Refs. 44,46,51 for details). For comparison purpose, calculations of the momentum distribution in the RPA (see Refs. 42,51 for details) have also been performed. The range of densities explored goes from r s = 0.1 to r s = 10 and that of interaction parameters from µ = 0.1 to µ = 30.

In order to assess the quality of the obtained momentum distributions, we look especially at two quantities that characterizes well the effect of correlations on the momentum distribution: the quasiparticle renormalization factor at the Fermi surface

Z F = n(k → k + F ) -n(k → k - F ), (11)
where k F = (3π 2 n) 1/3 is the Fermi wave-vector, and the fraction of electrons in the correlation tail (i.e., k > k F )

∆N N = 2 n ∞ kF dk (2π) 3 4πk 2 n(k). (12)
We also compute some moments of the momentum distribution: the second reduced moment δ 2 related to the kinetic energy, and the fourth reduced moment δ 4 . The dimensionless quantity δ 2 is defined as

δ 2 = ÊK -ÊK 0 ÊK 0 , (13)
where ÊK = (1/N ) kσ ε(k)c † kσ c kσ is the kinetic energy per particle operator, • • • 0 and • • • mean expectation values between the non-interacting and interacting ground states, respectively. Similarly, δ 4 is defined as

δ 4 = Ê2 K -Ê2 K 0 Ê2 K 0 , (14)
where

Ê2 K = (1/N ) kσ ε(k) 2 c † kσ c kσ .
It is these moments that appear in the asymptotic expansion of the static exchange-correlation kernel (see Sec. III).

We first present results for the Coulombic uniform electron gas and compare to available results in the literature in order to assess the accuracy of the G 0 W 0 approximation.

Fig. 1 reports the momentum distribution n coul (k) calculated at the RPA and G 0 W 0 levels for r s = 5. The parametrization of Gori-Giorgi and Ziesche (GZ) 52 using accurate data from Takada and Yasuhara 40 and known theoretical constraints is also reported. One sees that the G 0 W 0 calculation greatly improves the RPA and is very closed to the GZ parametrization, especially for the correlation tail.

The renormalization factor Z F,coul is plotted in Fig. 2. Z F,coul given by the RPA is largely underestimated and even negative at very low densities (r s > 7). As announced before, G 0 W 0 gives a renormalization factor slightly too large compared to that of the GZ parametrization which uses supposedly more accurate data of Takada and Yasuhara 40 for Z F,coul . The fraction of electrons in the correlation tail (∆N/N ) coul is reported in Fig. 3. The RPA largely overestimates (∆N/N ) coul . The G 0 W 0 approximation gives (∆N/N ) coul that is overall close to that given by the GZ parametrization, slightly deviating for large r s (r s > 8). As an independent check, we also report in Fig. 3, the fraction of electrons in the correlation tail extracted from a Coupled-Cluster calculation with double excitations (CCD) 53 . In this calculation, the energy of a uniform electron with unoccupied non-interacting levels shifted by a common gap is calculated and (∆N/N ) coul is deduced as the derivative of the energy at zero gap according to the Hellmann-Feynman theorem. The CCD calculation is close to both the GZ and G 0 W 0 results.

We now discuss results for the moments of the momentum distribution.

Using the virial theorem, the reduced moment δ 2,coul for the Coulomb interaction can be expressed exactly in term of the correlation energy of the uniform electron gas (see Appendix A). This enables us to perform another test of the accuracy of the G 0 W 0 calculation. Fig. 4 compares δ 2,coul obtained with the virial theorem [Eq. (A5)] using the usual Vosko-Wilk-Nusair (VWN) parametrization of the correlation energy 54 , with the GZ parametrization, and with the RPA and G 0 W 0 calculations. The RPA calculation gives very poor results for this quantity, especially for large r s . The GZ parametrization has been designed to reproduce the correct kinetic energy and consequently gives δ 2,coul identical to the formula from the virial theorem. The G 0 W 0 approximation agrees well with this result for all r s . Fig. 5 shows the reduced moment δ 4,coul given by the GZ parametrization and by the RPA and G 0 W 0 calculations. For comparison, a parametrization of δ 4,coul due to Farid, Heine, Engel and Robertson (FHER) 55 is also reported. The G 0 W 0 calculation agrees with the GZ parametrization for r s ≤ 5 but deviates from above for larger values of r s . The FHER parametrization is constantly below the GZ result. These results reflect the fact that δ 4,coul is difficult to compute accurately due to its important sensitivity to the correlation tail of the momentum distribution. In comparison to the RPA calculation, the G 0 W 0 approximation, as well as the GZ and FHER parametrizations, give "reasonable estimates" of δ 4,coul . For later use, δ 4,coul given by the GZ parametriza- Having assessed the precision of the G 0 W 0 approximation on the Coulombic electron gas, we now present the results for the long-range erf interaction. Fig. 6 reports the momentum distribution n µ erf in the RPA and G 0 W 0 approximations for r s = 5 and µ = 0.3. For comparison with Fig. 1, the G 0 W 0 calculation for the Coulomb interaction is also shown. The reduction of the interaction naturally brings the momentum distribution closer to the non-interacting momentum distribution n 0 (k) = θ(k F -k), even if the modifications are small at the scale of the plot, especially for the correlation tail. Also, the difference between RPA and G 0 W 0 is reduced compared to the Coulombic case at the same density.

The variation of the renormalization factor Z µ F,erf with the interaction parameter µ is reported in Fig. 7 for the RPA and G 0 W 0 approximation with r s = 5. One sees that Z µ F,erf is very sensitive with respect to µ near µ = 0. Thus, even a very small interaction introduces significant incoherence into the evolution of an electron in the system. The fraction of electrons in the correlation tail (∆N/N ) µ erf is plotted in Fig. 8 for the RPA and G 0 W 0 calculations. A CCD calculation similar of that of Ref. 53 but with the erf interaction is also reported. As for the Coulomb case, the G 0 W 0 and CCD agree well in all the range of µ.

The reduced moment δ µ 2,erf associated to the long-range erf interaction can also be expressed in term of the correlation energy via a generalization of the virial theorem to this modified interaction (see Appendix A). Fig. 9 compares δ µ 2,erf derived from the virial theorem [Eq. (A6)] using the parametrization of the correlation energy of Ref. 27 with the RPA and G 0 W 0 results. The G 0 W 0 calculation is in overall agreement with the formula from the virial theorem. The slight inaccuracy of G 0 W 0 in the Coulombic limit (µ → ∞) is inherent to the method which neglects part of the short-range interactions. On the contrary, the inaccuracy for small µ is of numerical origins: when the interaction is too small, the magnitude of the correlation tail is lower than the precision of the method (which does not exceed 10 -4 for the momentum distribution in the present implementation). Finally, the reduced moment δ µ 4,erf calculated in the RPA and G 0 W 0 approximations are reported in Fig. 10 for r s = 5. Analytical parametrizations of these results are given in Appendix B. The parametrization of δ µ 4,erf in the G 0 W 0 approximation will be used in the next section for the determination of the exchange-correlation kernel.

III. STATIC EXCHANGE-CORRELATION KERNEL

The dynamical exchange-correlation kernel f xc (k, ω), or equivalently the dynamical local-field factor G(k, ω) = -v(k)f xc (k, ω), is related to the linear dynamical (interacting) response function χ(k, ω) via the Dyson-like screening equation

χ(k, ω) = χ 0 (k, ω) + χ 0 (k, ω) [v(k) + f xc (k, ω)] χ(k, ω), (15)
where χ 0 (k, ω) is the non-interacting dynamical response function. The static exchange-correlation kernel is

f xc (k) = f xc (k, ω = 0).
The study of the local-field factor of the Coulom-bic electron gas has generated an abundant literature (see, e.g., Refs. 29,30,55-71). The QMC simulations 30,69 give probably one of the most reliable static exchangecorrelation kernel f xc,coul (k). A remarkable feature of the QMC results is that f xc,coul (k) does not have much structure, and can be essentially reproduced by combining its small-k constant value for k 2k F and its large-k twoterm asymptotic expansion for k 2k F 30 . For the case of the long-range erf interaction, accurate data are not available to check this property. However, it seems reasonable to assume that no additional structure appears in the exchange-correlation kernel when the interaction is reduced. As a first approximation, we will thus construct f µ xc,erf (k) from its limiting behaviors for the case of the erf interaction too.

The limiting behaviors of f xc (k) for both the Coulomb and long-range erf interactions can be cast into the same form. For k → 0, the static exchange-correlation kernel has the limit

f xc (0) = - 4π k 2 F A, (16)
where the dimensionless coefficient A has been introduced to retain the notation of the literature 30 . When k → ∞, f xc (k) has the following asymptotic expansion (up to k -2 order) (see Appendix C)

f ∞ xc (k) = - 4π k 2 F C - 4π k 2 B, (17)
where C and B are also dimensionless quantities.

The coefficients A, B, C naturally depend on the interaction chosen. A is given by the compressibility sum rule (see, e.g., Refs. 72,73) which is valid for any electronelectron interaction. C is related to the reduced second moment δ 2 of the momentum distribution. The expression of B involves, among other things, the second and fourth moments δ 2 and δ 4 .

For the Coulomb interaction, A coul is given by

A coul = - k 2 F 4π ∂ 2 (nε xc,coul ) ∂n 2 , (18)
with the exchange-correlation energy per particle ε xc,coul taken from the usual VWN parametrization 54 . The coefficients C coul and B coul have been calculated by Holas [START_REF] Holas | Strongly Coupled Plasma Physics[END_REF][START_REF] Holas | Electronic Structure of Solids[END_REF] [see Appendix C, Eq. (C25)] and read

C coul = k 4 F δ 2,coul 5ω 2 p , (19)
and

B coul = 2 3 (1 -g coul (0)) + 12k 4 F δ 4,coul 35ω 2 p - 4k 4 F (2δ 2,coul + (δ 2,coul ) 2 ) 25ω 2 p , (20)
where ω p = √ 4πn is the plasma frequency, g coul (0) is the on-top pair-distribution function taken from Ref. 76, and δ 2,coul is given by the virial formula of Eq. (A5). The different coefficients B coul resulting from the GZ, RPA and G 0 W 0 parametrizations of δ 4,coul [Eq. (B1)] are compared in Fig. 11. The parametrization of Moroni, Ceperley and Senatore (MCS) 30 is also reported. The G 0 W 0 calculation gives a reasonable estimation of B coul , in agreement with the GZ and MCS parametrizations.

The coefficients A coul , B coul C coul have been calculated (using the G 0 W 0 parametrization of δ 4,coul ) and fitted to simple analytical formulas for convenience in Appendix D.

For the long-range erf interaction, the coefficient A erf is

A µ erf = - k 2 F 4π ∂ 2 (nε µ xc,erf ) ∂n 2 , ( 21
)
where ε µ xc,erf is the exchange-correlation energy per particle associated to the erf interaction taken from Ref. 27. The calculation of Holas is generalized in Appendix C for the coefficients C erf and B erf [Eq. (C26)]; they write

C µ erf = k 4 F δ µ 2,erf 5ω 2 p , (22)
and

B µ erf = - 2µ 3 3πn ∂ε µ xc,erf ∂µ + 12k 4 F δ µ 4,erf 35ω 2 p - 4k 4 F (2δ µ 2,erf + (δ µ 2,erf ) 2 ) 25ω 2 p , (23)
where δ µ 2,erf is given by the virial formula of Eq. (A6). The three coefficients have the correct Coulombic limit when µ → ∞:

A µ→∞ erf = A coul , B µ→∞ erf = B coul and C µ→∞ erf = C coul .
The coefficients A µ erf , B µ erf and C µ erf have been calculated (using the G 0 W 0 parametrization for δ µ 4,erf [Eq. (B2)]) and fitted to analytical formulas in Appendix D. For both the Coulomb and long-range erf interactions, f xc (k) is then approximated over the whole range of k by the simple following stepwise interpolation between the two limiting behaviors at k → 0 and k → ∞

f xc (k) = f xc (0), f xc (0) > f ∞ xc (k) f ∞ xc (k), f xc (0) < f ∞ xc (k). ( 24
)
The junction point between the short and long wave vector regions where

f xc (0) = f ∞ xc (k) is located at k = √ γk F with γ = B/(A -C),
and Eq. ( 24) can be re-written in the more compact form

f xc (k) = - 4πA k 2 F -4π B k 2 + C -A k 2 F θ (k - √ γk F ) . (25)
In the case of the Coulomb interaction and for r s = 5, Fig. 12 compares interpolation of Eq. ( 25) with the parametrization of MCS [Eq. (7) of Ref. 30 with n = 8] and the parametrization of Corradini, Del Sole, Onida and Palummo (CDOP) 31 based on the same QMC data. The three curves are in overall agreement, the main discrepancy arising around k ≈ 2k F where the presence of a singularity or not is still a matter of debate [29][30][31]55 . Of course, our simple model cannot be trusted for k ≈ 2k F .

Taking the inverse Fourier transform of Eq. ( 25), we obtain the static exchange-correlation kernel for both the Coulomb and long-range erf interactions in real space

f xc (r) = - 4πC k 2 F δ(r) - B r 1 - 2 π Si( √ γk F r) + 2(C -A) πr 3 k 2 F [sin( √ γk F r) - √ γk F r cos( √ γk F r)] , (26)
where Si(x) =

x 0 sin(t)/t dt is the sine integral function. In Eq. ( 26), the second term gives the major contribution to f xc (r) for r = 0, the third term being only a small correction. Fig. 13 shows the static exchange-correlation kernel in reciprocal and real spaces for r s = 5 with the erf interaction for a series of interaction parameters µ. f µ xc,erf is naturally reduced when the interaction decreases. More precisely, in reciprocal space, f µ xc,erf (k) is flatten and the junction point is slightly shifted toward larger values of k, while, in real space, the spatial extension of f µ xc,erf (r) decreases. The kernel f µ xc,erf thus becomes more and more local when µ decreases. Symmetrically, the complement kernel f µ xc,erf = f xc,coul -f µ xc,erf also becomes more and more local when µ increases (not shown).

In the limit of the Coulomb interaction, the kernel of Eq. ( 26) displays Friedel-like long-range oscillations in real space of quasiperiodicity √ γk F ≈ 2k F coming mainly from the second term in Eq. ( 26). These oscillations connected with the behavior of the kernel in reciprocal space around k ≈ 2k F constitute an usual feature of models for the exchange-correlation kernel of the uniform electron gas. In our model, these oscillations are likely to be exaggerated due to the sharp cut-off in Eq. ( 25) and their relevance for inhomogeneous systems is questionable. However, the oscillations in our model for f µ xc,erf (r) rapidly spread out and vanish when µ decreases, i.e. when the interaction is reduced, (see Fig. 13) so that this problem disappears. Symmetrically, the oscillations in the the complement kernel f µ xc,erf (r) also disappears when µ increases. This results are consistent with the common intuition that the decrease of the interaction reduces the structure of the exchange-correlation kernel. Note that in the uniform electron gas the interaction can also be reduced by decreasing r s so that an equivalent way to look at f µ xc,erf (r) in the limit of a small interaction, µ → 0 for a fixed r s , is to look at the high-density limit, r s → 0, for a fixed µ.

IV. CORRELATION ENERGY FROM STATIC EXCHANGE-CORRELATION KERNEL

As an example of the usefulness of the static exchangecorrelation kernel, we follow Lein, Gross and Perdew 77 and compute the correlation energy of the uniform electron gas from it.

Indeed, the correlation energy per particle can be exactly deduced in principle from the dynamical exchangecorrelation kernel using the ACFD approach. For comparison purposes, we begin by the standard case of the Coulomb interaction for which the expression is (see, e.g., Ref. 77) 25)] and direct [Eq. ( 26)] spaces for rs = 5 with the erf interaction (dashed curves) for µ = 0.2, 0.3, 0.5 and 1.0. For comparison, the case of the Coulomb interaction (solid curve) is also shown.

ε ACFD c,coul = -1 2πn 1 0 dλ dq (2π) 3 v coul (q) × ∞ 0 dω v λ coul (q) + f λ xc,coul (q, iω) χ 0 (q, iω) 2 1 -v λ coul (q) + f λ xc,coul (q, iω) χ 0 (q, iω) , (27)
with the scaled interaction v λ coul (q) = λv coul (q) and the associated exchange-correlation kernel f λ xc,coul (q, ω) related to f xc,coul (q, ω) by a simple scaling relation 77 28) We can thus use Eq. ( 27) to estimate the relevance of our simple model of the static exchange-correlation kernel to obtain correlation energies, even though in theory the frequency-dependence is also needed to obtain the exact correlation energy. Fig. 14 27), the correlation energy is greatly improved compared to the RPA. The exchange-correlation kernel of Eq. ( 25) gives results very close to those given by the MCS and CDOP parametrizations, validating our simple interpolation.

f λ xc,coul [n](q, ω) = λ -1 f xc,coul [n/λ 3 ](q/λ, ω/λ 2 ). (
We now generalize the ACFD approach to the case of the long-range erf interaction. In this case, the long-range correlation energy per particle is given by

ε µ,ACFD c,erf = -1 2πn 1 0 dλ dq (2π) 3 v µ erf (q) × ∞ 0 dω v µ,λ erf (q) + f µ,λ xc,erf (q, iω) χ 0 (q, iω) 2 1 -v µ,λ erf (q) + f µ,λ xc,erf (q, iω) χ 0 (q, iω) , (29)
with the scaled interaction v µ,λ erf (q) = λv µ erf (q) and the associated exchange-correlation kernel f µ,λ xc,erf (q, ω). The scaling relation of Eq. ( 28) is easily generalized to the erf interaction (see also Ref. 81)

f µ,λ xc,erf [n](q, ω) = λ -1 f µ/λ xc,erf [n/λ 3 ](q/λ, ω/λ 2 ). (30)
The error in the long-range correlation energy obtained with Eq. ( 29), ∆ε µ c,erf = ε µ,ACFD c,erf -ε µ c,erf , in the RPA (f µ xc,erf = 0) and with the exchange-correlation kernel of Eq. ( 25) are represented in Fig. 15 for µ = 1. The reduction of the interaction mainly acts at small r s decreasing the error given by the RPA. As for the Coulomb case, the exchange-correlation kernel of Eq. ( 25) enables to correct the RPA for larger values of r s .

V. SUMMARY AND CONCLUSION

We have given a simple approximate expression in real and reciprocal spaces for static exchange-correlation kernel f µ xc,erf of a uniform electron gas interacting with the long-range part only of the Coulomb interaction [Eq. (25) and Eq. (26]. This expression interpolates between the exact asymptotic behaviors of the kernel at small and large wave vectors which in turn requires, among other thing, information from the momentum distribution of the uniform electron gas with the same interaction that we have calculated in the G 0 W 0 approximation. In the limit of the Coulomb interaction (µ → ∞), the proposed exchange-correlation kernel f µ→∞ xc,erf = f xc,coul is closed to other parametrizations of the literature based on QMC data. As a matter of course, the complement kernel f µ xc,erf = f xc,coul -f µ xc,erf associated to the short-range part of the Coulomb interaction is also obtained in this work.

In the same way that the static exchange-correlation kernel of the Coulombic uniform electron gas is an essential ingredient in ACFD and DFT approaches of the electronic correlation problem, the modified kernels f µ xc,erf and f µ xc,erf can be used in alternative ACFD and DFT approaches based on a separate treatment of long-range and short-range correlations. As a simple illustration, the long-range correlation energy of the uniform electron gas has been calculated from f µ xc,erf . In comparison to the Coulombic kernel f xc,coul , the modified kernels f µ xc,erf and f µ xc,erf are more local which may improve their transferability to inhomogeneous systems and therefore facilitate their use in applications. We are currently investigating such an issue.

Coulomb interaction

For the Coulomb interaction, δ 4,coul is parametrized as

δ 4,coul = 6 i=3 d i x i , (B1) where x

erf interaction

For the erf interaction, we take the parametrization In this appendix, following the procedure proposed by Holas [START_REF] Holas | Strongly Coupled Plasma Physics[END_REF][START_REF] Holas | Electronic Structure of Solids[END_REF] , we derive the asymptotic expansion of the static exchange-correlation kernel f xc (k) of the uniform electron gas for k → ∞ for the erf interaction. For comparison, results for the Coulomb interaction are also given. Contrary to the rest of this work and except for Eqs. (C25) and (C26), momentum are expressed in units of k F , energies and frequencies in units of k 2 F throughout this appendix.

δ µ 4,erf = 6 i=3 e i x i , (

Frequency moments of the linear response function

Knowledge of the frequency moments of the frequencydependent linear response function χ(k, ω) is useful to study the asymptotic behavior of the static response function for k → ∞. These moments are defined as

M l (k) = 2 π ∞ 0 Im χ(k, ω)ω l dω. (C1)
The zeroth, first and third moments are wellknown [START_REF] Holas | Strongly Coupled Plasma Physics[END_REF][START_REF] Holas | Electronic Structure of Solids[END_REF]84 M 0 (k) = -2nS(k), (C2)

M 1 (k) = -nk 2 , (C3) M 3 (k) = -nk 2 k 4 2 + 2k 2 ÊK + nk 2 v(k) -4πnG PV (k) , (C4)
where S(k) is the static structure factor, ÊK = (1/N ) kσ ε(k)c † kσ c kσ is the kinetic energy per particle operator, • • • means expectation value between the interacting ground state, and G PV (k) is the "Pathak-Vashishta" 84 (or "Niklasson" 85 ) local-field factor

G PV (k) = G PV a (k) + G PV b (k), (C5)
with

G PV a (k) = 1 n dp (2π) 3 (k • p) 2 k 2 v(p) 4π (S(p) -1) , (C6)
and

G PV b (k) = - 1 n dp (2π) 3 (k • p) 2 k 2 v(p) 4π (S(|k -p|) -1) .
(C7) The asymptotic expansion of M 0 (k) as k → ∞ is given by that of S(k) 86

M 0 (k) = -2n 1 + C k 4 + o 1 k 4 . ( C8
)
For the Coulomb interaction, C is directly related to the cusp condition: C coul = -2ω 2 p g ′ (0) = -2ω 2 p g(0) where ω p = √ 4πn is the plasma frequency, g(0) and g ′ (0) are the on-top values of the pair-distribution function and of its spherical-average derivative. For the cusp-less erf interaction, C vanishes: C erf = 0.

The asymptotic expansion of M 3 (k) writes

M 3 (k) = - nk 6 4 1 + 8 k 2 ÊK + 4 k 4 nv(k) - 16 k 4 πnG PV (k → ∞) + o 1 k 4 , (C9)
where the limit G PV (k → ∞) is determined as follows.

G PV a (k) is actually independent of k and writes

G PV a = 1 3 1 n dp (2π) 3 p 2 v(p) 4π (S(p) -1) , (C10)
which gives, for the Coulomb interaction,

G PV a,coul = 1 3 1 n dp (2π) 3 (S coul (p) -1) = 1 3 (g coul (0) -1) , (C11)
and, for the erf interaction,

G PV,µ a,erf = 1 3 1 n dp (2π) 3 e -p 2 /(4µ 2 ) (S µ erf (p) -1) = µ 3 6πn dp (2π) 3 ∂v µ erf (p) ∂µ (S µ erf (p) -1) = µ 3 3πn ∂ε µ xc,erf ∂µ , (C12)
where the Hellmann-Feynman theorem has been used for the last line, and ε µ xc,erf is the exchange-correlation energy per particle with the erf interaction.

After a trivial variable transformation, the limit of G PV b as k → ∞ is easily seen to be

G PV b (k → ∞) = - 1 n dp (2π) 3 (k -p) 2 v(|k -p|) 4π × (S(p) -1) , (C13)
which gives, for the Coulomb interaction,

G PV b,coul (k → ∞) = - 1 n dp (2π) 3 (S coul (p) -1) = -(g coul (0) -1) , (C14)
and, for the erf interaction,

G PV,µ b,erf (k → ∞) = - 1 n dp (2π) 3 e -(k-p) 2 /(4µ 2 ) × (S µ erf (p) -1) = -e -k 2 /(4µ 2 ) (g µ erf (0) -1) , ( C15
) where g µ erf (0) is the on-top pair-distribution function associated to the erf interaction.

Static linear response function

The static linear response function χ(k, 0) can be expressed by the spectral representation

χ(k, 0) = 2 π ∞ 0 Im χ(k, ω) ω dω, (C16)
i.e., we simply have χ(k, 0) = M -1 (k). Using this result, Holas [START_REF] Holas | Strongly Coupled Plasma Physics[END_REF][START_REF] Holas | Electronic Structure of Solids[END_REF] showed that the asymptotic expansions for large k of M 0 (k), M 1 (k) and M 3 (k) [Eqs. (C8), (C3) and (C9)] are sufficient to determine the asymptotic expansion of χ(k, 0) up to k -4 order

χ(k, 0) = -4n k 2 1 + 8 3 ÊK 1 k 2 + o 1 k 2 . (C17)
He then determined the following terms in the asymptotic expansion of χ(k, 0) by (not rigorously) inferring them as follows. The following term in k -3 in the square-bracket of Eq. ( C17) is set to zero by a simple argument involving the continuity and the likely non-oscillating behavior of the moments M l (k) with respect to l. The next term in k -4 is inferred from the asymptotic expansion of the static response function calculated to first-order with respect to the electron-electron interaction. We follow the same procedure. The (zeroth-order) Lindhard static response function χ 0 (k, 0) is well-known. Its asymptotic expansion for large k writes

χ 0 (k, 0) = -4n k 2 1 + 8 3 ÊK 0 1 k 2 + 64 5 Ê2 K 0 1 k 4 + o 1 k 4 , (C18)
where Ê2 K = (1/N ) kσ ε(k) 2 c † kσ c kσ and • • • 0 means expectation value between the non-interacting ground state.

The asymptotic expansion of the first-order correction to the static response function χ 1 (k, 0) has been calculated for the Coulomb interaction by Geldart and Taylor 59,87 (see also Ref. 63). This result is conveniently generalized to an arbitrary interaction as

χ 1 (k, 0) = -4n k 2 -4ω 2 p 2G PV a,HF + k 2 v(k) 4π 1 k 4 +o 1 k 4 , (C19)
with

G PV a,HF = 1 3 1 n dp (2π) 3 p 2 v(p) 4π (S HF (p) -1) , (C20)
and the Hartree-Fock (HF) static structure factor

S HF (p) = 1 - 2 n dk (2π) 3 n 0 (k)n 0 (|k + p|). (C21)
For the Coulomb interaction, we simply have G PV a,HF,coul = (g HF,coul (0) -1)/3 where g HF,coul (0) = 1/2 is the HF on-top pair-distribution function. For the erf interaction, G PV,µ a,HF,erf = µ 3 /(3πn)∂ε µ x,erf /∂µ where ε µ x,erf is the exchange energy per particle associated to this modified interaction. Following Holas [START_REF] Holas | Strongly Coupled Plasma Physics[END_REF][START_REF] Holas | Electronic Structure of Solids[END_REF] , we then infer the following term in k -4 in the square-bracket of Eq. (C17) from the corresponding term in the expansion of χ 0 (k, 0)+χ 1 (k, 0) with the substitutions Ê2 K 0 → Ê2 K and G PV a,HF → G PV a . We thus finally obtain

χ(k, 0) = -4n k 2 1 + 8 3 ÊK 1 k 2 + 64 5 Ê2 K -4ω 2 p 2G PV a + k 2 v(k) 4π 1 k 4 + o 1 k 4 . (C22)
Observe that the form of the asymptotic expansion of χ(k, 0) [Eq. (C22)] does depend explicitly on the interaction v(k) and is therefore different for the Coulomb and erf interactions.

Static exchange-correlation kernel

The static exchange-correlation kernel f xc (k) writes [see Eq. ( 15)] with ÊK 0 = 3/10 and Ê2 K 0 = 3/28, we obtain (with k in units of k F and ω p in units of k 2 F ) To facilitate the evaluation of the static exchangecorrelation kernel for the Coulomb and erf interaction, we give analytical parametrizations of the coefficients A, B and C determining its limiting behaviors for small and large k (see Sec. III).

f xc (k) = χ 0 (k, 0) -1 -χ(k, 0) -1 -v(k), (C23)
f xc (k) = -4πδ 2 5ω 2 p - 4π k 2 -2G PV a + 12δ 4 35ω 2 p - 4(2δ 2 + δ 2 2 ) 25ω 2 p +o 1 k 2 . ( C24

Coulomb interaction

For the Coulomb interaction, the coefficient A [Eq. (18)] is parametrized as follows

A coul = 6 i=0 a i x i , (D1)
with x = √ r s and the fitted parameters a 0 = 0.250019, a 1 = -0.000162, a 2 = 0.013441, a 3 = -0.003591, a 4 = 0.000380, a 5 = 0.000002, a 6 = -0.000003.

For the coefficient B (Eq. 20), we adopt the analytical form of Moroni, Ceperley and Senatore 30 incorporating the correct limit B(r s → 0) = 1/3 (19)] is parametrized as

B coul = 1 + b 1 x + b 2 x 3 3 + b 3 x + b 4 x 3 , (
C coul = 6 i=1 c i x i 1 + 4 i=1 f i x i , ( D3
)
with c 1 = 0.002127, c 2 = 0.169597, c 3 = 0.450771, c 4 = -0.023265, c 5 = 0.001855, c 6 = -0.000069, f 1 = 7.062604, f 2 = 8.589773, f 3 = 2.747407, f 4 = 0.648920.

erf interaction

For the erf interaction, A µ erf [Eq. (21)] is represented by the analytical parametrization

A µ erf = 6 i=2 g i x i 1 + h 4 x 4 , (D4)
with g i = g i1 µ + g i2 µ 2 and h 4 = g 41 µ + h 42 µ 2 . The fitted parameters are g 21 = -0.029315, g 22 = -0.000927, g 31 = 0.061867, g 32 = 0.010970, g 41 = -0.053761, g 42 = 0.078580, g 51 = 0.012970, g 52 = 0.014669, g 61 = -0.001232, g 62 = -0.000891, h 41 = -0.025963, h 42 = 0.389673. The coefficient B µ erf [Eq. ( 23)] is parametrized as

B µ erf = 3 i=0 j i x i 1 + k 1 x + k 3 x 3 , ( D5
)
where j i = (j i1 µ + j i2 µ 2 )/(1 + j i3 µ 2 ) and k i = (k i1 µ + k i2 µ 2 )/(1+k i3 µ 2 ). Using the parametrization of Eq. (B2) for δ µ 4 in the G 0 W 0 approximation, the fitted parameters obtained are j 01 = 0.010533, j 02 = -0.002640, j 03 = 0.314403, j 11 = -0.143455, j 12 = 0.046302, j 13 = 0.014315, j 21 = -0.415043, j 22 = 0.194149, j 23 = 0.078849, j 31 = 0.164085, j 32 = 0.925083, j 33 = 0.491612, k 11 = 0.291633, k 12 = 0.102905, k 13 = 0.013551, k 31 = 0.014324, k 32 = 0.714449, k 33 = 0.473130.

Finally, the coefficient C µ erf [Eq. (22)] is parametrized as

C µ erf = 6 i=1 l i x i 1 + m 4 x 4 , (

FIG. 1 :FIG. 2 :

 12 FIG.1: Momentum distribution with the Coulomb interaction for rs = 5: GZ parametrization 52 (solid curve), RPA (short-dashed curve) and G0W0 (long-dashed curve) calculations.

FIG. 3 :FIG. 4 :

 34 FIG.3: Fraction of electrons in the correlation tail of the moment distribution with the Coulomb interaction: GZ parametrization52 (solid curve), RPA (short-dashed curve), G0W0 (long-dashed curve) and CCD53 (dotted curve) calculations.

FIG. 5 : 3 FIG. 6 :

 536 FIG.5: δ 4,coul with the Coulomb interaction: GZ 52 (solid curve) and FHER55 (dotted curve) parametrizations, RPA (short-dashed curve) and G0W0 (long-dashed curve) calculations.

FIG. 7 :

 7 FIG.7: Renormalization factor with the erf interaction for rs = 5: RPA (short-dashed curve) and G0W0 (long-dashed curve) calculations.

FIG. 8 :

 8 FIG.8: Fraction of electrons in the correlation tail of the momentum distribution with the erf interaction: RPA (shortdashed curve), G0W0 (long-dashed curve) and CCD (dotted curve) calculations.

FIG. 9

 9 FIG. 9: δ µ 2,erf with the erf interaction for rs = 5: formula from virial theorem [thick solid curve, Eq. (A6)], RPA (shortdashed curve) and G0W0 (long-dashed curve) calculations.

FIG. 10

 10 FIG. 10: δ µ4,erf with the erf interaction for rs = 5: RPA (shortdashed curve) and G0W0 (long-dashed curve) calculations.

FIG. 11 :

 11 FIG.11: Coefficient B coul with the Coulomb interaction: GZ (solid curve) and MCS30 (dotted curve) parametrizations, RPA (short-dashed curve) and G0W0 (long-dashed curve) calculations.

FIG. 12 :

 12 FIG. 12: Static exchange-correlation kernel f xc,coul (k) with the Coulomb interaction for rs = 5: simple interpolation formula [Eq. (25), dashed curve], MCS 30 (solid curve) and CDOP 31 (dotted curve) parametrizations.

FIG. 13 :

 13 FIG. 13: Static exchange-correlation kernel f µ xc,erf in reciprocal [Eq. (25)] and direct [Eq. (26)] spaces for rs = 5 with the erf interaction (dashed curves) for µ = 0.2, 0.3, 0.5 and 1.0. For comparison, the case of the Coulomb interaction (solid curve) is also shown.

1 + 1 =

 11 reports the error in the correlation obtained with Eq. (27), ∆ε c,coul = ε ACFD c,coul -ε c,coul , using several approximations for f xc,coul : RPA (f xc,coul = 0), MCS and CDOP parametrizations and the interpolation formula of Eq. (25). The expansion of the exact correlation energy for r s → 0 (see, e.g., Refs. 78-80) is ε c,coul = C 0 ln r s + C 1 + O(r s ln r s ) with C 0 = (1 -ln 2)/π 2 and C 1 = -0.0469205 while in the RPA it is ε RPA c,coul = C 0 ln r s + C RPA O(r s ln r s ) with C RPA -0.071100. The error in the RPA correlation energy at r s → 0 is due to the difference C X 1 = C 1 -C RPA 1 = 0.024179 corresponding to the contribution of exchange diagrams. When a static exchange-correlation kernel is used in Eq. (

FIG. 14 :

 14 FIG. 14: Error in the correlation energy ∆ε c,coul = ε ACFD c,coulε c,coul for the Coulomb interaction [Eq. (27)] with different exchange-correlation kernels: RPA [fxc = 0, shortdashed curve], MCS 30 (solid curve) and CDOP 31 (dotted curve) parametrizations and the simple interpolation formula [Eq. (25), long-dashed curve].

FIG. 15 :

 15 FIG. 15: Error in the correlation energy ∆ε µ c,erf = ε µ,ACFD c,erf ε µ c,erf for the erf interaction [Eq. (29)] with different exchangecorrelation kernels: RPA [fxc = 0, short-dashed curve], the simple interpolation formula [Eq. (25), long-dashed curve], for µ = 1.

  = √ r s . A least-square fit of the result of the GZ parametrization gives d 3 = 0.271191, d 4 = -0.009998, d 5 = -0.036383, d 6 = 0.006706. For the RPA calculation we obtain d 3 = 0.093623, d 4 = 0.194288, d 5 = 0.051445, d 6 = 0.005449, and for the G 0 W 0 calculation, d 3 = 0.126362, d 4 = 0.001428, d 5 = 0.014278, d 6 = -0.004522.

  B2) with e i = (e i1 µ + e i2 µ 2 )/(1 + e i3 µ 2 ). The optimal parameters reproducing the RPA calculation are e 31 = -0.690727, e 32 = 0.393525, e 33 = 3.398631, e 41 = 0.231529, e 42 = 0.885698, e 43 = 5.999882, e 51 = 0.001233, e 52 = -0.061940, e 53 = 2.887302, e 61 = 0.001083, e 62 = -0.000060, e 63 = 0.391712. Those for the G 0 W 0 calculation are e 31 = -0.525472, e 32 = 0.417720, e 33 = 5.102281, e 41 = 0.223039, e 42 = 0.316492, e 43 = 5.984490, e 51 = 0.002486, e 52 = -0.000770, e 53 = 0.265086, e 61 = -0.004998, e 62 = -0.003175, e 63 = 1.191797.

  Appendix C: Asymptotic expansion of the static exchange-correlation kernel

  and its asymptotic expansions for large k is therefore determined from the expansions of χ 0 (k, 0) and χ(k, 0) [Eqs. (C18) and (C22)]. Introducing the quantities δ 2 = ( ÊK -ÊK 0 )/ ÊK 0 and δ 4 = ( Ê2 K -Ê2 K 0 )/ Ê2 K 0

)

  The term explicitly depending on v(k) in the asymptotic expansion of χ(k, 0) [Eq. (C22)] exactly cancels with the Hartree kernel v(k) in the expression of f xc (k) [Eq. (C23)]. Therefore, the asymptotic expansion of f xc (k) [Eq. (C24)] do has the same form for the Coulomb and erf interactions. Note however that the local field factor G(k) = -v(k)f xc (k) involves explicitly the interaction and has therefore different asymptotic behaviors for the Coulomb and erf interactions.Coming back to atomic units now and specializing to the Coulomb interaction, Eq. (C24) becomesf xc,coul (k) = -4πk

D2) but the parameters are refitted with the parametrization of δ 4 [

 4 Eq. (B1)] with the G 0 W 0 data: b 1 = 0.721543, b 2 = 0.317320, b 3 = -0.133379, b 4 = 0.269494. The coefficient C [Eq.

  D6) where l i = (l i1 µ + l i2 µ 2 )/(1 + l i3 µ 2 ) and m 4 = (m 41 µ + m 42 µ 2 )/(1 + m 43 µ 2 ). The fitted parameters are l 11 = -0.039662, l 12 = 0.002346, l 13 = 0.405492, l 21 = 0.187782, l 22 = 0.066673, l 23 = 0.759503, l 31 = -0.270823, l 32 = -0.083482, l 33 = 0.511806, l 41 = 0.133824, l 42 = 0.061658, l 43 = 0.348085, l 51 = -0.027392, l 52 = 0.024384, l 53 = 0.261739, l 61 = 0.001877, l 62 = -0.005601, l 63 = 0.351272, m 41 = -0.090916, m 42 = 1.238974, m 43 = 0.251286.

Acknowledgments

I thank A. Savin (Paris, France) who initiated this work, P. Gori-Giorgi (Paris, France) for useful advices and in particular for Eqs. (C12) and (C15), F. Colonna (Paris, France), J. Ángyán (Nancy, France) and J. Dobson (Griffith, Australia) for interesting discussions. I am also grateful to M. Vogt and R. Needs (Cambridge, UK) who provided me with technical details on GW calculations.

Appendix A: Kinetic energy from virial theorem In this appendix, using the virial theorem, we give the expression of δ 2 [Eq. (13)] for the Coulomb and erf interactions.

The virial theorem

The virial theorem for the uniform electron gas with the Coulomb interaction writes 82,83

where t coul (r s ) = t 0 (r s ) + t c,coul (r s ) is the interacting kinetic energy per particle which can be decomposed into a non-interacting contribution t 0 (r s ) and a correlation contribution t c,coul (r s ), v coul (r s ) is the electronelectron repulsion energy per particle and

is the total energy per particle which can be decomposed into the non-interacting kinetic energy t 0 (r s ) and the exchangecorrelation energy per particle ε xc,coul (r s ). Eliminating v coul (r s ) for ε xc,coul (r s ) in Eq. (A1) and using the virial theorem for the non-interacting electron gas, 2t 0 (r s ) = -r s dt 0 (r s )/dr s , and the virial relation for exchange only, ε x,coul (r s ) = -r s dε x,coul (r s )/dr s , we arrive at the virial theorem for the correlation part (see also Ref. 55)

For the erf interaction, the virial theorem of Eq. (A1) is generalized as (see Ref. 81)

Following the same steps as for the Coulomb interaction, it can be shown that the virial theorem for the correlation contributions reads now

.

(A4)

Expression of δ2

Using t c,coul (r s ) = -d(r s ε c,coul (r s ))/dr s coming from Eq. (A2) and t 0 (r s ) = 3k 2 F /10, the expression of δ 2,coul for the Coulomb interaction writes

For the erf interaction, Eq. (A4) gives t µ c,erf (r s ) = -∂(r s ε µ c,erf (r s ))/∂r s + µ∂ε µ c,erf (r s )/∂µ and the expression of δ µ 2,erf is

+ µ ∂ε µ c,erf (r s ) ∂µ .

(A6)

Appendix B: Analytical parametrization of δ4

In this appendix, we give analytical parametrizations for the reduced moment δ 4 of the momentum distribution [Eq. ( 14)] calculated by several methods for the Coulomb and erf interactions.