János G Angyán
email: angyan@crm2.uhp-nancy.fr

Ru-Fen Liu

Julien Toulouse
email: julien.toulouse@upmc.fr

Georg Jansen
email: georg.jansen@uni-due.de

János G Ángyán

Correlation energy expressions from the adiabatic-connection fluctuation-dissipation theorem approach

come

Correlation energy expressions from the adiabatic-connection fluctuation-dissipation theorem approach

Introduction

There is a recent revival of interest in the random phase approximation (RPA) to obtain groundstate correlation energies of electronic systems. The RPA is considered as a promising first approximation to obtain non-perturbative, exact-exchange-compatible, post-Kohn-Sham correlation energy corrections in density-functional theory. In particular, the RPA is thought of as a remedy for the bad description of London dispersion forces by conventional local and semi-local density-functional approximations. However, it is widely admitted that while RPA is well adapted to long-range electron-electron interactions, for small interelectronic distances its performance is even poorer than that of semi-local density functionals. 42,43 An efficient way to make an optimal use of RPA is to apply it in a range-separated approach, 44,45 where the short-range interactions are described by an exchange-correlation density functional, and long-range exchange and correlation are treated by Hartree-Fock (HF) and RPA, respectively. Computational schemes following these principles have been recently proposed and applied mainly to van der Waals complexes. [15][16][17]28,31,33,46 Several formulations of RPA have been developed. Perhaps, the most well-known approach to RPA is the one based on the adiabatic-connection fluctuation-dissipation theorem (ACFDT). 47,48 In this approach, the correlation energy expression involves integrations over both the frequency and the interaction strength, which can be performed either numerically or analytically. Obviously, an expression which has already been integrated analytically along at least one or both of these variables is more advantageous than the repeated use of numerical quadratures. If an analytical integration over the frequency is performed first, followed by a numerical integration over the interaction strength, one obtains an expression that is of the form of an interaction-strength-averaged twoparticle density matrix contracted with the two-electron integrals. This is the adiabatic-connection formulation. An analytical integration over the interaction strength followed by a numerical integration along the frequency leads to an expression involving the dynamic dielectric matrix. This is the dielectric-matrix formulation. With a second analytical integration (either along the interaction strength starting from the adiabatic-connection expression, or along the frequency starting from the dielectric-matrix expression) both of these intermediate forms can be converted to a common expression, which consists in a sum of the shifts of electronic excitation energies when passing from an independent-particle to the RPA description of the excited states. This is the plasmon formulation. The plasmon expression can be further converted to an equivalent expression involving coupled-cluster doubles (CCD) amplitudes calculated in the ring-diagram approximation. 14 This is the ring CCD formulation. The relationship between the adiabatic-connection and ring CCD formulations of RPA has been recently discussed in Ref. 34.

In this work, we study different variants of RPA within the adiabatic-connection formulation, which differ in the way the exchange interactions are handled. If the exchange interactions are neglected in the density matrix, we obtain the direct RPA (dRPA) approach (also called timedependent Hartree), while inclusion of the non-local HF exchange response kernel leads to the RPAx approach (also called time-dependent Hartree-Fock, or full RPA). A third possibility, not discussed here, consists in including an exact exchange response kernel from a local exact exchange potential. 27 If the dRPA density matrix is contracted with non-antisymmetrized two-electron integrals, one obtains what we call the dRPA-I variant, while if it is contracted with antisymmetrized two-electron integrals, one obtains the dRPA-II variant. Similarly, if the RPAx density matrix is contracted with non-antisymmetrized two-electron integrals, the RPAx-I variant is obtained, while if it is contracted with antisymmetrized two-electron integrals, one obtains the RPAx-II variant.

The dRPA-I variant is just the commonly called "RPA" of the density-functional/material-science community. The dRPA-II variant, which is similar to the second-order screened exchange (SO-SEX) expression introduced by Grüneis et al. 23 in the ring CCD formulation, is original to this work. In contrast to SOSEX, it involves higher-order screened exchange effects. The RPAx-II variant was first introduced by McLachlan and Ball, 49 but here we derive a new adiabatic-connection expression for it. Finally, the RPAx-I variant has been recently introduced by Toulouse et al. 15,33 When possible, for the case of dRPA-I and RPAx-II, we also compare with the equivalent plasmon formulation, and clarify the origin of the prefactor of 1/4 in the plasmon formula of RPAx-II in place of the prefactor of 1/2 appearing for dRPA-I. We remind the reader that in spite of the very different formulations, the dRPA-I variant is the same as the direct ring-CCD method, while the RPAx-II approach is identical to ring-CCD. 34,46 For the sake of simplicity, we give all the expressions without range separation, but it is straightforward to generalize them for the case of range separation, as done in Ref. 33. The paper is organized as follows. In Sec. II, we first provide an overview of the adiabatic-connection RPA correlation energy variants. In Sec. III, we review how the two-particle density matrix is obtained from the RPA polarization propagator. In Sec. IV, we derive the expressions of RPA correlation energy variants in spin-orbital basis. In Sec. V, we derive the corresponding spin-adapted expressions for closed-shell systems. In Sec. VI, we perform numerical comparisons of different variants on a few atomic and molecular systems with and without range separation. Finally, Sec. VII contains our conclusions. The analysis of the second-order limit in the electron-electron interaction of each variant is given in Appendix.

Overview of RPA correlation energy variants in the adiabaticconnection formulation

In the adiabatic-connection formalism, the correlation energy in a spin-orbital basis can be expressed as

E c = 1 2 1 0 dαTr V P c,α = 1 2 1 0 dα pq,rs rq|sp P c,α pq,rs , (1)
where V sr,qp = rq|sp are the two-electron integrals, P c,α is the correlation part of the two-particle density matrix at interaction strength α, and Tr denotes the trace (sum over the indices rs and pq). Using the antisymmetry of P c,α with respect to the permutation of the indices p and s, the correlation energy can also be expressed as

E c = 1 4 1 0 dαTr W P c,α = 1 4 1 0 dα pq,rs rq||sp P c,α pq,rs , (2)
where W sr,qp = rq||sp = rq|sp -rq|ps are the antisymmetrized two-electron integrals. In RPAtype approximations, P c,α is obtained via the fluctuation-dissipation theorem

P RPA c,α = - ∞ -∞ dω 2πi e iω0 + Π RPA α (ω) -Π 0 (ω) , (3)
where Π RPA α (ω) is the four-index matrix representation of dynamic polarization propagator at interaction strength α and frequency ω, and Π 0 (ω) is the corresponding non-interacting (Hartree-Fock or Kohn-Sham) polarization propagator. In the dRPA variant (or time-dependent Hartree) the polarization propagator is obtained from the response equation with the Hartree kernel V

Π dRPA α (ω) -1 = Π 0 (ω) -1 -αV, (4)
whereas in the RPAx variant (or time-dependent Hartree-Fock) the polarization propagator is obtained using the Hartree-Fock kernel W

Π RPAx α (ω) -1 = Π 0 (ω) -1 -αW . (5)
The obtained dRPA and RPAx correlation density matrices P dRPA The dRPA-I variant is obtained by inserting the dRPA correlation density matrix in 1,

E dRPA-I c = 1 2 1 0 dαTr V P dRPA c,α . ( 6
)
This variant is commonly called "RPA" in the density-functional/material-science community. It corresponds to the first RPA correlation energy approximation historically developed and is still widely used. Since the dRPA response equation involves the mere Hartree kernel, only the screening effect of the bare Coulomb interaction is taken into account in the polarization propagator and all exchange-correlation screening effects are neglected. The resulting correlation energies tend to be too strongly negative. At second order in the electron-electron interaction, the dRPA-I correlation energy does not reduce to the standard second-order Møller-Plesset (MP2) correlation energy, but instead to a "direct MP2" expression, i.e. without the MP2 exchange term. 2,50 The dRPA-II variant is obtained by contracting the dRPA correlation density matrix with the antisymmetrized two-electron integrals W ,

E dRPA-II c = 1 2 1 0 dα Tr W P dRPA c,α , (7)
which re-establishes the correct second-order MP2 limit. 2 could have suggested to use a factor of 1/4 instead of 1/2 in 7, but in fact the correct MP2 limit is only recovered with the factor 1/2. This variant can also be obtained from 6 by antisymmetrizing the correlation density matrix with respect to the permutation of p and s: (P dRPA ) sq,r p . As far as we know, the dRPA-II variant has never been described before. It is similar to the second-order screened exchange (SOSEX) expression introduced by Grüneis et al. 23 but the latter does not involve integration over the adiabatic connection and treats exchange effects only at the lowest order of perturbation.

The RPAx-I variant is obtained by inserting the RPAx correlation density matrix in 1,

E RPAx-I c = 1 2 1 0 dα Tr V P RPAx c,α , (8)
and has been introduced recently by Toulouse et al. 15,33 In this variant, the exchange screening effects are taken into account in the polarization propagator. The matrix P RPAx c,α is properly antisymmetric at first order, and therefore the RPAx-I correlation energy has the correct MP2 limit. At propagator at interaction strength α writes

Π RPA α (ω) = -(Λ α -ω ∆) -1 , (12)
where the supermatrix Λ α is calculated with the Hartree kernel V in the case of dRPA,

Λ dRPA α = Λ 0 + αV, (13)
and with the Hartree-Fock kernel W in the case of RPAx,

Λ RPAx α = Λ 0 + αW . ( 14
)
From now on, we will consider real-valued orbitals. In this case, the Hartree kernel is made of four identical blocks,

V =            K K K K            , (15)
where K ia, jb = ab|i j are non-antisymmetrized two-electron integrals. Similarly, the Hartree-Fock kernel writes

W =            A ′ B B A ′            , (16)
with the antisymmetrized two-electron integrals

A ′ ia, jb = ib||a j = ib|a j -ib| ja = K ia, jb -J ia, jb , (17)
and

B ia, jb = ab||i j = ab|i j -ab| ji = K ia, jb -K ′ ia, jb . (18)
Let us consider now the generalized non-hermitian RPA eigenvalue equation

Λ α C α,n = ω α,n ∆ C α,n , (19)
whose solutions come in pairs: positive excitation energies ω α,n with eigenvectors C α,n = x α,n , y α,n and negative excitation energies ω α,-n = -ω α,n with eigenvectors C α,-n = y α,n , x α,n . The spectral representation of Π RPA α (ω) then writes

Π RPA α (ω) = n        C α,n C T α,n ω -ω α,n + i0 + - C α,-n C T α,-n ω -ω α,-n -i0 +        , ( 20
)
where the sum is over eigenvectors n with positive excitation energies ω α,n > 0. The fluctuationdissipation theorem [3] leads to the supermatrix representation of the correlation density matrix P RPA c,α (using contour integration in the upper half of the complex plane)

P RPA c,α = - ∞ -∞ dω 2πi e iω0 + [Π RPA α (ω) -Π 0 (ω)] = n C α,-n C T α,-n -C 0,-n C T 0,-n , (21)
with the non-interacting eigenvectors C 0,-n = y 0,n , x 0,n with y 0,n = 0 and x 0,n = 1 n (where 1 n is the vector whose n th component is 1 and all other components are zero). The explicit supermatrix expression of the RPA correlation density matrix is thus

P RPA c,α =            Y α Y T α Y α X T α X α Y T α X α X T α            -            0 0 0 I            , ( 22
)
where X α and Y α are the matrices whose columns contain the eigenvectors x α,n and y α,n . The dRPA and RPAx correlation density matrices have the same form in terms of the eigenvector matrices X α and Y α , although the eigenvectors are of course different for dRPA and RPAx.

Correlation energy expressions in spin-orbital basis

We give here the expressions in a spin-orbital basis for calculating the different RPA correlation energy variants. We first consider the dRPA-I and RPAx-II variants which have similar expressions.

In both cases the integration over the adiabatic connection can be done analytically, leading to plasmon formulae. We then examine the dRPA-II and RPAx-I variants. They have in common that they are mixing the non-antisymmetrized integrals V and the antisymmetrized integrals W , which makes it impossible to do the integration over the adiabatic connection analytically. Although the dRPA-I variant is well-documented in the literature after the work of Furche and coworkers, 2,13,32 the review that we give here is useful to define our notations and for comparisons with other variants. The RPAx-I variant has been discussed in detail in the context of range separation by Toulouse et al. 15,31,33 The RPAx-II variant is much less documented and the dRPA-II is new, so most of the expressions that we give for them are original to this work.

dRPA-I correlation energy

There are several equivalent expressions for the dRPA-I correlation energy.

Adiabatic-connection formula

The dRPA-I correlation energy of 6 can be expressed with the eigenvectors of the dRPA polarization propagator according to the general prescription to form the correlation density matrix, 21,

E dRPA-I c = 1 2 1 0 dα n Tr V C α,-n C T α,-n -V C 0,-n C T 0,-n , (23)
or, using the explicit expressions in terms of the block matrix components [15 and 22],

E dRPA-I c = 1 2 1 0 dα tr (X α + Y α ) (X α + Y α ) T -I K , (24)
where tr refers to the trace now applied to the block matrices (which are half the size of the supermatrices). As shown by Furche, 2 one does not need to calculate explicitly the eigenvector matrices X α and Y α to get the correlation energy; it is sufficient to form the matrix

Q α = (X α + Y α ) (X α + Y α ) T , (25)
which can be obtained directly from the matrices involved in the RPA response equation. In the case of dRPA, it simply reads

Q dRPA α = ε ε ε 1/2 M dRPA α -1/2 ε ε ε 1/2 , (26)
with

M dRPA α = ε ε ε 1/2 (ε ε ε + 2 α K) ε ε ε 1/2 . ( 27
)
The adiabatic-connection formula for the dRPA-I correlation energy is then finally

E dRPA-I c = 1 2 1 0 dα tr Q dRPA α -I K . (28)
In previous papers, this equation was written with the matrix

P dRPA α = Q dRPA α -I.

Plasmon formula

The plasmon formula for the dRPA-I correlation energy is found by starting from an equivalent form of 23,

E dRPA-I c = 1 2 1 0 dα n Tr C T α,-n V C α,-n -C T 0,-n V C 0,-n , (29)
obtained by a cyclic permutation of the matrices in the trace. Since the positive excitation energies can be written as 13,49

ω dRPA α,n = C T α,-n Λ dRPA α C α,-n , (30)
the derivative of ω α,n with respect to α gives

dω dRPA α,n dα = C T α,-n dΛ dRPA α dα C α,-n = C T α,-n VC α,-n , (31)
which allows one to perform the integral over α in 29 analytically, leading to the plasmon formula

E dRPA-I c = 1 2 n ω dRPA 1,n -ω 0,n -C T 0,-n VC 0,-n = 1 2 n ω dRPA 1,n -ω dTDA n , (32)
where

n ω dTDA n = n C T 0,-n Λ dRPA 1 C 0,-n = n ω 0,n + C T 0,-n VC 0,-n
is the sum of the (positive) excitation energies in the direct Tamm-Dancoff approximation (dTDA). The sum of the dTDA excitation energies can also be expressed as n ω dTDA n = tr{ε ε ε + K}.

Alternative plasmon formula

An alternative form of the plasmon formula can be found by rewriting 32 as

E dRPA-I c = 1 2 n Tr Λ dRPA 1 C 1,-n C T 1,-n -Λ dRPA 1 C 0,-n C T 0,-n , (33)
where the cyclic invariance of the trace has again been used. Using then 22 and recalling that the diagonal blocks of Λ dRPA 1 are ε ε ε + K and the off-diagonal blocks are K, the correlation energy becomes

E dRPA-I c = 1 2 tr Y 1 Y T 1 + X 1 X T 1 -I (ε ε ε + K) + Y 1 X T 1 + X 1 Y T 1 K . ( 34
)
Introducing now the inverse of the Q α matrix, 43

Q -1 α = (X α -Y α ) (X α -Y α ) T , (35)
which in the case of dRPA can be written as

Q dRPA α -1 = ε ε ε -1/2 M dRPA α 1/2 ε ε ε -1/2 , (36)
the correlation energy can be expressed as

E dRPA-I c = 1 2 tr 1 2 Q dRPA 1 + Q dRPA 1 -1 -I (ε ε ε + K) + 1 2 Q dRPA 1 -Q dRPA 1 -1 K , (37)
or, equivalently,

E dRPA-I c = 1 2 tr Q dRPA 1 -I K + 1 2 Q dRPA 1 + Q dRPA 1 -1 -2I ε ε ε , (38)
or, rearranged in a different way

E dRPA-I c = 1 2 tr 1 2 Q dRPA 1 (ε ε ε + 2K) + 1 2 Q dRPA 1 -1 ε ε ε -(ε ε ε + K) . (39)
Using the expressions of Q dRPA 1 [26], (Q dRPA 1 ) -1 [36], and M dRPA 1 [27], and the cyclic invariance of the trace, we finally arrive at the alternative form of the plasmon formula for the dRPA-I correlation energy . Indeed, the potential correlation energy is just the value of the integrand in 28 at α = 1, i.e.

E dRPA-I c = 1 2 tr M dRPA 1 1/2 -(ε ε ε + K) . (40)
U dRPA-I c = 1 2 tr Q dRPA 1 -I K , (41)
and thus, by subtraction, according to 38, the kinetic correlation energy is

T dRPA-I c = 1 4 tr Q dRPA 1 + Q dRPA 1 -1 -2I ε ε ε . ( 42
)
In the limit of a system with orbitals that are all degenerate, i.e. with static correlation only, then ε ε ε = 0 and the kinetic correlation energy vanishes as it should. This is in agreement with the statement that dRPA-I correctly describes left-right static correlation in bond dissociations. 7,51

RPAx-II correlation energy

We now derive several equivalent RPAx-II correlation energy expressions by proceeding in an analogous way to the case of dRPA-I.

Adiabatic-connection formula

The RPAx-II correlation energy of 9 can be written in terms of the eigenvectors of the RPAx polarization propagator

E RPAx-II c = 1 4 1 0 dα n Tr W C α,-n C T α,-n -W C 0,-n C T 0,-n , (43)
or, using the block structure of W [16],

E RPAx-II c = 1 4 1 0 dα tr Y α Y T α + X α X T α -I A ′ + Y α X T α + X α Y T α B . (44)
Using the matrix Q α which in the case of RPAx is given by

Q RPAx α = ε ε ε + αA ′ -αB 1/2 M RPAx α -1/2 ε ε ε + αA ′ -αB 1/2 , (45)
with

M RPAx α = ε ε ε + αA ′ -αB 1/2 ε ε ε + αA ′ + αB ε ε ε + αA ′ -αB 1/2 , (46)
and the inverse

Q -1 α Q RPAx α -1 = ε ε ε + αA ′ -αB -1/2 M RPAx α 1/2 ε ε ε + αA ′ -αB -1/2 , ( 47
)
we arrive at the adiabatic-connection formula for the RPAx-II correlation energy

E RPAx-II c = 1 4 1 0 dα tr 1 2 Q RPAx α A ′ + B + 1 2 Q RPAx α -1 A ′ -B -A ′ . (48)
Since Q α = I+P α , if P α is small, we can consider the approximation

Q -1 α = (I+P α ) -1 ≈ I-P α = 2I -Q α ,
which leads to the following approximation for the RPAx-II correlation energy

E RPAx-IIa c = 1 4 1 0 dα tr 1 2 Q RPAx α A ′ + B + 1 2 2I -Q RPAx α A ′ -B -A ′ = 1 4 1 0 dα tr Q RPAx α -I B . (49)
So, we have the interesting result that this approximate correlation energy expression is analogous to the dRPA-I correlation energy expression of 28, the only differences being that the matrix Q α is now obtained from the RPAx response equation and that it is contracted with the antisymmetrized two-electron integrals B, along with the corresponding change of the prefactor from 1/2 to 1/4.

Plasmon formula

As in the case of dRPA-I, the plasmon formula for the RPAx-II correlation energy is found by taking profit of the cyclic invariance of the trace to rewrite 43 as

E RPAx-II c = 1 4 1 0 dα n Tr C T α,-n W C α,-n -C T 0,-n W C 0,-n , (50)
and then using

dω RPAx α,n /dα = C T α,-n (dΛ RPAx α /dα)C α,-n = C T α,-n W C α,-n to integrate analytically over α E RPAx-II c = 1 4 n ω RPAx 1,n -ω 0,n -C T 0,-n W C 0,-n = 1 4 n ω RPAx 1,n -ω TDAx n , (51)
where

n ω TDAx n = n C T 0,-n Λ RPAx 1 C 0,-n = n ω 0,n + C T 0,-n W C 0,-n
is the sum of the (positive) excitation energies in the Tamm-Dancoff approximation with exchange (TDAx) or configuration interaction singles (CIS). The sum of the TDAx excitation energies can also be expressed as n ω TDAx n = tr {ε ε ε + A ′ }. This plasmon formula was first presented by McLachlan and Ball. 49 The presence of a factor of 1/4 in 51 and not a factor of 1/2 like in 32 has been debated in the literature. 52 The present exposition makes it clear that this factor of 1/4 is due to the use of the antisymmetrized two-electron integrals W .

Alternative plasmon formula

As in the case of dRPA-I, the alternative plasmon formula is found by rewriting 51 as

E RPAx-II c = 1 4 n Tr Λ RPAx 1 C 1,-n C T 1,-n -Λ RPAx 1 C 0,-n C T 0,-n , (52)
and inserting the diagonal blocks of Λ RPAx 1 which are ε ε ε + A ′ and the off-diagonal blocks which are B,

E RPAx-II c = 1 4 tr 1 2 Q RPAx 1 ε ε ε + A ′ + B + 1 2 Q RPAx 1 -1 ε ε ε + A ′ -B -(ε ε ε + A ′ ) . ( 53
)
Using the expressions of Q RPAx 1 [45], Q RPAx 1 -1 [47], and M RPAx 1 [46], and the cyclic invariance of the trace, we arrive at the alternative plasmon formula for the RPAx-II correlation energy

E RPAx-II c = 1 4 tr M RPAx 1 1/2 -ε ε ε + A ′ . (54)
Finally, just as for dRPA-I, comparison of 48 and 53 provides us with a decomposition of the correlation energy into the potential energy contribution to the correlation energy

U RPAx-II c = 1 4 tr 1 2 Q RPAx 1 A ′ + B + 1 2 Q RPAx 1 -1 A ′ -B -A ′ , (55)
and the kinetic correlation energy

T RPAx-II c = 1 8 tr Q RPAx 1 + Q RPAx 1 -1 -2I ε ε ε . ( 56
)
The RPAx-II kinetic correlation energy vanishes in the limit where ε ε ε = 0 as for dRPA-I.

dRPA-II correlation energy

The dRPA-II correlation energy of 7 writes in terms of the eigenvectors of the dRPA polarization propagator

E dRPA-II c = 1 2 1 0 dα n Tr W C α,-n C T α,-n -W C 0,-n C T 0,-n , (57)
leading to

E dRPA-II c = 1 2 1 0 dα tr 1 2 Q dRPA α A ′ + B + 1 2 Q dRPA α -1 A ′ -B -A ′ . (58)
Equation ( 58) is similar to 48, with Q dRPA α instead of Q RPAx α and a factor 1/2 instead of 1/4.

The approximation Q -1 α ≈ 2 I -Q α leads to the following approximate dRPA-II correlation energy

E dRPA-IIa c = 1 2 1 0 dα tr Q dRPA α -I B , (59)
which is in close analogy (but usually not equal) to the SOSEX correlation energy in the ring-CCD formulation. The analytic relationship of this "adiabatic-connection SOSEX" (AC-SOSEX)

variant with the original SOSEX has been discussed in detail in Ref. 34.

RPAx-I correlation energy

Finally, the RPAx-I correlation energy of 8 writes in terms of the eigenvectors of the RPAx polarization propagator

E RPAx-I c = 1 2 1 0 dα n Tr V C α,-n C T α,-n -V C 0,-n C T 0,-n , (60)
leading to

E RPAx-I c = 1 2 1 0 dα tr Q RPAx α -I K , (61)
which has the same form than 28 but with the RPAx matrix Q RPAx α . This last variant has been discussed in detail and applied in the context of range-separated density-functional theory. 15,31,33 5 Correlation energy expressions in spatial-orbital basis for closed-shell systems

For spin-restricted closed-shell calculations, all the matrices in the spin-orbital excitation basis occurring in the RPA equations have the following spin block structure

C =                              C ↑↑,↑↑ C ↑↑,↓↓ 0 0 C ↓↓,↑↑ C ↓↓,↓↓ 0 0 0 0 C ↑↓,↑↓ C ↑↓,↓↑ 0 0 C ↓↑,↑↓ C ↓↑,↓↑                              . ( 62
)
This structure is a consequence of the fact that the two-electron integrals can be non-zero only for pairs of identical spins. The orthogonal transformation

U = 1 √ 2                              1 1 0 0 1 -1 0 0 0 0 1 1 0 0 1 -1                              , (63)
leads to a spin-adapted matrix C = U T C U, which in the case of the matrices involved in RPA simplifies into a block-diagonal form with a spin-singlet excitation block 1 C and three spin-triplet excitation blocks 3,0 C, 3,1 C, and 3,-

1 C C =                              1 C 0 0 0 0 3,0 C 0 0 0 0 3,1 C 0 0 0 0 3,-1 C                              , (64)
with the matrix elements (i, j and a, b referring now to occupied and virtual spatial orbitals, respectively)

1 C ia, jb = 1 2 (C i↑a↑ j↑b↑ + C i↑a↑ j↓b↓ + C i↓a↓ j↑b↑ + C i↓a↓ j↓b↓ ), ( 65a
) 3,0 C ia, jb = 1 2 (C i↑a↑ j↑b↑ -C i↑a↑ j↓b↓ -C i↓a↓ j↑b↑ + C i↓a↓ j↓b↓ ), ( 65b
) 3,±1 C ia, jb = 1 2 (C i↑a↓ j↑b↓ ± C i↑a↓ j↓b↑ ± C i↓a↑ j↑b↓ + C i↓a↑ j↓b↑ ). ( 65c
)
Let us start with dRPA. Spin-adaptation of the non-antisymmetrized integrals matrix K gives only a contribution from the singlet excitations where 1 A ′ ia, jb = 2 ib|a j -ib| ja , 3 A ′ ia, jb = -ib| ja , 1 B ia, jb = 2 ab|i j -ab| ji , and 3 B ia, jb = -ab| ji . Notice the minus sign for the last triplet block in the B matrix which makes spin-adaptation less trivial for RPAx. By 46, it leads to the following spin-adaptation for the matrix

M RPAx α MRPAx α =                              1 M RPAx α 0 0 0 0 3 M RPAx α 0 0 0 0 3 M RPAx α 0 0 0 0 3 N RPAx α                              , (70)
with the expected spin-adapted blocks

1 M RPAx α = ε ε ε + α 1 A ′ -α 1 B 1/2 ε ε ε + α 1 A ′ + α 1 B ε ε ε + α 1 A ′ -α 1 B 1/2 and 3 M RPAx α = ε ε ε + α 3 A ′ -α 3 B 1/2 ε ε ε + α 3 A ′ + α 3 B ε ε ε + α 3 A ′ -α 3 B 1/2
, along with the less expected last triplet block with opposite signs for 3 B,

3 N RPAx α = ε ε ε + α 3 A ′ + α 3 B 1/2 ε ε ε + α 3 A ′ -α 3 B ε ε ε + α 3 A ′ + α 3 B 1/2
. By 45, it gives the following spin-adaptation for the matrix

Q RPAx α QRPAx α =                              1 Q RPAx α 0 0 0 0 3 Q RPAx α 0 0 0 0 3 Q RPAx α 0 0 0 0 3 Q RPAx α -1                              , (71)
with the spin-adapted blocks 1

Q RPAx α = ε ε ε + α 1 A ′ -α 1 B 1/2 1 M RPAx α -1/2 ε ε ε + α 1 A ′ -α 1 B 1/2 and 3 Q RPAx α = ε ε ε + α 3 A ′ -α 3 B 1/2 3 M RPAx α -1/2 ε ε ε + α 3 A ′ -α 3 B 1/2 .
The last triplet block turns out to be the inverse 3

Q RPAx α -1 = ε ε ε + α 3 A ′ + α 3 B 1/2 3 N RPAx α -1/2 ε ε ε + α 3 A ′ + α 3 B
1/2 since according to 25 and 35 one goes from Q α to Q -1 α by changing the sign of Y α which is equivalent to changing the sign of B.

The spin-adapted correlation energy expressions can be easily obtained by using the invariance of the trace under the transformation C → U T CU. The spin-adapted adiabatic-connection formula for the dRPA-I correlation energy is thus

E dRPA-I c = 1 2 1 0 dα tr 1 Q dRPA α -I 1 K , (72)
i.e. only singlet excitations contribute. Similarly, the corresponding plasmon formula contains only singlet excitation energies

E dRPA-I c = 1 2 n 1 ω dRPA 1,n -1 ω dTDA n . ( 73
)
The triplet term vanishes since both 3 ω dRPA 1,n and 3 ω dTDA n are equal to the one-particle excitation energies ε aε i . Finally, the spin-adapted alternative plasmon formula is

E dRPA-I c = 1 2 tr 1 M dRPA 1 1/2 -ε ε ε + 1 K . ( 74
)
Both singlet and triplet excitations contribute the RPAx-II correlation energy. The spin-adapted adiabatic-connection formula for RPAx-II is

E RPAx-II c = 1 4 1 0 dα tr 1 2 1 Q RPAx α 1 A ′ + 1 B + 1 2 1 Q RPAx α -1 1 A ′ -1 B -1 A ′ + 3 4 1 0 dα tr 1 2 3 Q RPAx α 3 A ′ + 3 B + 1 2 3 Q RPAx α -1 3 A ′ -3 B -3 A ′ . ( 75
)
The last triplet term gives a contribution identical to the other two triplet terms because the expression is invariant under the replacements

Q α → Q -1 α and B → -B.
The spin-adaptation of the approximate RPAx-II correlation energy of 49 is

E RPAx-IIa c = 1 4 1 0 dα tr 1 Q RPAx α -I 1 B + 2 4 1 0 dα tr 3 Q RPAx α -I 3 B - 1 4 1 0 dα tr 3 Q RPAx α -1 -I 3 B , (76)
where now the last triplet term is not identical to the other two triplet terms. If we make the

additional approximation 3 Q RPAx α -1 ≈ 2I -3 Q RPAx α
we arrive at the following expression

E RPAx-IIb c = 1 4 1 0 dα tr 1 Q RPAx α -I 1 B + 3 4 1 0 dα tr 3 Q RPAx α -I 3 B , (77)
which could also have been obtained by starting from the spin-adapted formula of 75 and making the approximation Q -1 α ≈ 2I -Q α in both the singlet and the triplet terms. The RPAx-II plasmon formula decomposes into sums over singlet and triplet excitation energies

E RPAx-II c = 1 4 n 1 ω RPAx 1,n -1 ω TDAx n + 3 4 n 3 ω RPAx 1,n -3 ω TDAx n , (78)
and similarly for the alternative plasmon formula

E RPAx-II c = 1 4 tr 1 M RPAx 1 1/2 -ε ε ε + 1 A ′ + 3 4 tr 3 M RPAx 1 1/2 -ε ε ε + 3 A ′ . ( 79
)
The last triplet term is identical to the other two because 3 N RPAx The spin-adapted dRPA-II correlation energy involves only singlet excitations

E dRPA-II c = 1 2 1 0 dα tr 1 2 1 Q dRPA α 1 A ′ + 1 B + 1 2 1 Q dRPA α -1 1 A ′ -1 B -1 A ′ , (80)
since for the triplet blocks 3 Q dRPA α = I and the contribution vanishes. Likewise, the spin-adaptation of the approximate dRPA-II correlation energy of 59 is simply

E dRPA-IIa c = 1 2 1 0 dα tr 1 Q dRPA α -I 1 B . (81)
Finally, the spin-adapted RPAx-I correlation energy expression is

E RPAx-I c = 1 2 1 0 dα tr 1 Q RPAx α -I 1 K , (82)
where only single excitations contribute since the triplet blocks of the matrix K are zero.

Numerical illustrations

The above-described spin-adapted RPA correlation energy variants based on the adiabatic-connection formula have been implemented in the development version of the MOLPRO quantum chemistry package. [START_REF] Werner | version 2010.2, a package of ab initio programs[END_REF] The numerical equality of the alternative but equivalent expressions has been carefully tested and has been confirmed within the usual accuracy of quantum chemical calculations. In each case, we start by doing a usual Kohn-Sham (KS) calculation with some approximate density functional, and evaluate the RPA correlation energy with the KS orbitals. The total RPA energy is calculated as

E RPA tot = E EXX + E RPA c , (83)
where E EXX is the exact exchange (EXX) energy expression evaluated with the same KS orbitals.

This exchange energy is Hartree-Fock type, and it is not to be confused with the optimized effective potential (OEP) type local exchange, often denoted by the same acronym. For comparison, we also perform range-separated calculations, in which we start from a range-separated hybrid (RSH), 45 using the short-range PBE exchange-correlation functional of Ref. 54, and add the long-range RPA correlation energy evaluated with RSH orbitals

E RSH+RPA tot = E RSH + E lr,RPA c . ( 84
)
The long-range RPA correlation energy E lr,RPA c is calculated by replacing the Coulombic twoelectron integrals by the two-electron integrals with the long-range interaction erf(µr)/r, just as in Refs. 15,31,33. We use a fixed value of the range-separation parameter of µ = 0.5 bohr -1 .

This value corresponds to a reasonable global compromise, as it has been shown previously 55 by a study of thermochemical properties, and as it has been confirmed later by using alternative criteria leading to similar estimates of the µ parameter (see, e.g. Ref. 56). In all cases, the adiabaticconnection integration is performed by a 8-point Gauss-Legendre quadrature.

The RPA correlation energies are extrapolated to the complete basis set (CBS) limit by the usual 1/X 3 formula 57 for a series of Dunning basis sets. In contrast to the usual two-point extrapolation procedure 58,59 all the available correlation energies calculated by at least triple zeta basis set are used. The single-determinant reference energies are evaluated with a large basis set so that they can be considered as converged.

Atomic correlation energies

As a first test, we have calculated correlation energies for a series of atoms and atomic cations and compared with full configuration interaction (FCI) quality correlation energies as estimated by Davidson and coworkers. 60,61 In order to make a direct comparison with the FCI-quality correlation energies which are defined with respect to the HF energies, we redefine RPA correlation energies as the difference between the total RPA energies and the regular HF energies. The singledeterminant reference energies are calculated with a large even-tempered basis set. With this basis 60,61 with and without range separation. All the correlation energies have been extrapolated to the CBS limit. The RPA correlation energies E * c (RPA) are redefined here as the difference between the total RPA energies and the regular HF energies.

set, the HF energies agree within all significant digits with Davidson's reference data. Core excitations are included in the calculation of the RPA correlation energies and are extrapolated from the series of aug-cc-pCVXZ basis sets for He up to X=6, for B + , Al + , Ne, Ar up to X=5 and for Li + , Na + , Be, Mg up to X=Q. 1 (a)-(c) show the ratios of the correlation energies for each full-range RPA variant (dRPA-I, dRPA-II, dRPA-IIa, RPAx-I, RPAx-II) to the FCI-quality correlation energies, using orbitals obtained with the local-density approximation (LDA), 62 the Perdew-Burke-Ernzerhof (PBE), 63 and the Zhao-Morrison-Parr (ZMP) 64 exchange-correlation potentials. The ZMP potentials have been constructed from high-quality ab initio wave functions (Brueckner coupled cluster doubles). [START_REF] Boese | Private communication. The ZMP potentials have been prepared by the methodology described in 68,69 using densities calculated with Brueckner CCD method[END_REF] It appears that the correlation energies are only marginally dependent on the quality of the KS orbitals, at least for this series of atomic systems. The full-range RPAx-I and RPAx-II variants suffer from instabilities in the RPAx response equation for the Be, B + , Mg, and Al + systems, and additionally Ar in the case of RPAx-II with the ZMP orbitals. In fact, the strongly overestimated RPAx-II correlation energies of Ar obtained with the LDA and PBE orbitals indicate a situation close to an instability. More generally, the presence of near instabilities may be considered as being at the origin of the relatively strong overestimation (usually more than 150 %) of the correlation energy in RPAx-II. In view of the poor performance of RPAx-II, we did not test the approximate versions of 76 and 77. The RPAx-I variant only involves singlet excitations, and thus is not subject to triplet instabilities. It gives quite reasonable correlation energies (maximum 25% of overestimation) for He, Li + , Ne, Na + , and even for Ar. However, RPAx-I is subject to singlet instabilities which appear for the rest of the systems. The dRPA-I variant is free of any instability problems, since the dRPA response matrix is positive definite by construction, but has nevertheless a tendency for overestimating correlation energies by a factor of 1.5 to 2. This systematic error be easily corrected by including exchange in the energy expression. In fact, the dRPA-II variant and especially its approximate form dRPA-IIa (AC-SOSEX) lead to a very good reproduction of the reference correlation energies. Similar effects could be observed recently in the direct ring-CCD (dRPA-I) and SOSEX calculations of correlation energies by Klopper et al., 40 performed with a much smaller basis set.

As mentioned previously, dRPA-IIa (or AC-SOSEX) and the ring-CCD-based SOSEX correlation energies are expected to be quite close to each other. Numerical results (not shown in the figures) confirm this expectation. For two-electron systems (He, Li + ) the dRPA-IIa and SOSEX correlation energies are identical, while for the rest of the systems the relative difference is less than 0.15 %. The largest absolute difference, 0.82 mHartree, has been found in full-range calculations on the Ne atom. It is interesting to note that the ring-CCD based SOSEX correlation energies are always lower than the dRPA-IIa values. This fact cannot be interpreted simply by the comparison of the third order energy expressions, reported in Ref. 34 1 (d) shows the same total correlation energies obtained with range separation, i.e. the sum of the short-range PBE correlation energy and the long-range RPA correlation energy. The situation is quite different from the full-range RPA calculations. First, we do not encounter any instability problems anymore. Second, all the range-separated RPA variants give essentially identical correlation energies. Third, the correlation energies are systematically underestimated, for most of the systems with less than 20% of error, but with the notable exceptions of Li + , Be, and B + , for which the correlation energies are underestimated by as much as 50%. These findings may be due to the fact the systems considered here have very compact densities, and for the value of the range separation used here, µ = 0.5 bohr -1 , the major part of correlation is assigned to the short-range density functional rather than to the long-range RPA calculation. Improvement over these results would require either increasing the value of µ or using a more accurate short-range density-functional approximation.

Bond lengths and harmonic vibrational frequencies

London dispersion interactions

3 shows the interaction energy curves of the three rare-gas dimers He 2 , Ne 2 , and Ar 2 , calculated with the full-range and range-separated RPA variants. The full-range RPA calculations are done with PBE orbitals, while the range-separated RPA calculations are done with the short-range PBE density functional. All RPA calculations are without core excitations, and extrapolated to the CBS limit with the series of basis sets aug-cc-pVXZ with X=T,Q,5,6. The single-determinant reference energies are obtained with the aug-cc-pV6Z basis set. We note that when using LDA orbitals (not shown), instabilities are found for Ne 2 and Ar 2 in a wide range of interatomic distances. In contrast, no instabilities are encountered in the case of PBE, neither with nor without range separation.

The continuous curves without points represent on the one hand the accurate reference curves according the analytical potential energy expression of He 2 dimer without range separation, for which the EXX energy is too strongly repulsive. The prerequisite of the good performance of the range-separated calculations is obviously the excellent accuracy of the RSH energy, which, for He 2 , is in almost perfect agreement with the reference repulsive curve.

The full-range RPAx-II variant overestimates systematically the binding energy by a factor of 3 or more. The dRPA-I method largely underestimates the interaction energies and for He 2 it does not provide any minimum at all, although the non-binding character is mostly due to the bad single-determinant energy. The dRPA-II variant systematically gives more binding than dRPA-I but also tends to underestimate the interaction energies. The approximate dRPA-IIa variant gives results that are always very close to those of dRPA-I. This is not surprising since the dRPA-I and the dRPA-IIa methods differ only by the presence of exponentially decaying exchange integrals in the interaction matrix which become quite rapidly negligible to the interaction energy in van der Waals complexes. This behavior is analogous to that of the SOSEX method which gives dispersion interaction energies also very close to those of dRPA-I. 46 The best full-range method for these rare gas dimers is RPAx-I which is in quite good agreement with the reference curves for Ne 2 and Ar 2 .

With range separation, all the RPA variants give much closer interaction energy curves to each other, but the same trends are found. Range-separated dRPA-I, dRPA-II, and dRPA-IIa methods systematically underestimate interaction energies, the range-separated RPAx-II significantly overbinds Ar 2 , and the range-separated RPAx-I globally gives the most accurate interaction energies.

Conclusions

We have analyzed several RPA correlation energy variants based on the adiabatic-connection formula: dRPA-I , dRPA-II, RPAx-I, and RPAx-II. These variants have the generic form of an interaction-strength-averaged two-particle density matrix contracted with two-electron integrals.

They differ in the way the exchange interactions are treated. The dRPA-I variant is just the usual separated RPAx-I is the best variant for dispersion interaction energies.

We hope that the overview of the RPA correlation energy variants provided in this work will be useful for a better understanding of RPA methods and can serve as a starting point for the design of improved approximations.

all virtual (a and b) orbitals are real. In the following the index pairs ia and jb will be replaced with simple indices m and n. Note that the matrices are symmetric: A ′ n,m = A ′ m,n and B n,m = B m,n . The solutions of 85 come in pairs, i.e., if C α,n = x α,n , y α,n is an eigenvector with a positive eigenvalue ω n,α > 0, then C α,-n = y α,n , x α,n is an eigenvector with the negative eigenvalue ω α,-n = -ω α,n .

In the following we will use positive integer indices to denote solutions which connect to positive eigenvalues in the limit of a vanishing coupling parameter α, i.e., to ω 0,n > 0. Note that we also suppose a non-vanishing HOMO-LUMO gap.

The positive energy solutions of 85 for α = 0 are trivially given by ω 0,n = ε n , x 0,n = 1 n and y 0,n = 0, where 1 n denotes the n-th unit vector, i.e., a vector with vanishing components except for the n-th component which is equal to one. We now wish to find the first-order correction C (1) n to the eigenvector employing the power-series Ansatz ω α,n = ω 0,n + αω (1) n + . . . ,

C α,n = C 0,n + αC (1) n + . . . . (87)
Plugging this into 85 one sees that the first-order corrections are obtained from solving Λ 0 C (1) n + W C 0,n = ω 0,n ∆C (1) n + ω (1) n ∆C 0,n .

Multiplication of this equation from the left with C T 0,n and using C T 0,n Λ 0 C (1) n = ω 0,n C T 0,n ∆C (1) n along with the normalization condition C T 0,n ∆C 0,n = 1 gives the first-order correction to the eigenvalue

ω (1) n = C T 0,n W C 0,n = A ′ n,n . (90)
Multiplying 89 from the left with C T 0,m for m n, using C T 0,m Λ 0 C (1) n = ω 0,m C T 0,m ∆C (1) n , and employing the orthogonalization condition C T 0,m ∆C 0,n = 0 leads to

C T 0,m ∆C (1) n = - C T 0,m W C 0,n ω 0,m -ω 0,n , (91)
provided that the zeroth-order eigenvalues are non-degenerate, i.e., that no two occupied-virtual orbital energy differences match. Repeating the same operations for C T 0,-m one arrives at

C T 0,-m ∆C (1) n = C T 0,-m W C 0,n ω 0,m + ω 0,n , (92)
where ω 0,-m = -ω 0,m has been used. Using the resolution of identity, 1 = m C 0,m C T 0,m + -m C 0,-m C T 0,-m , the orthogonality of C (1) n to the zeroth-order eigenvector, i.e. C T 0,n ∆C (1) n = 0, and ∆ 2 = 1, we find the expansion of the first-order correction to the positive-energy eigenvectors

C (1) n = - m n C T 0,m W C 0,n ω 0,m -ω 0,n ∆C 0,m + m C T 0,-m W C 0,n ω 0,m + ω 0,n ∆C 0,-m . (93)
From 93 it follows that the first-order corrections read more explicitly

x (1) n = - m n A ′ m,n ε m -ε n 1 m , (94a)
y (1) n = - m B m,n ε m + ε n 1 m . (94b)
The first-order corrections to the negative-energy solutions are simply:

ω (1) -n = -ω (1) n , x (1)
-n = y (1) n , and y (1) -n = x (1) n . We can obtain the first-order expansion of the matrix

Q RPAx α Q RPAx α = n (x α,n + y α,n )(x α,n + y α,n ) T = n 1 n 1 T n + α n x (1) n 1 T n + 1 n x (1) n T + y (1) n 1 T n + 1 n y (1) n T + O(α 2 ), (95)
where the sum over n refers to positive-energy eigenvectors only. The first term is simply the identity matrix

n 1 n 1 T n = I. (96)
Using 94a, one can show that the term depending on x (1) n vanishes n

x (1)

n 1 T n + 1 n x (1) n T = - n m n A ′ m,n ε m -ε n 1 m 1 T n - n m n A ′ m,n ε m -ε n 1 n 1 T m = 0. (97)
This is seen by swapping n and m in the last term and noting that A ′ m,n /(ε mε n ) is antisymmetric when exchanging m and n. Finally, using 94b, the term depending on y (1) n gives n y (1)

n 1 T n + 1 n y (1) n T = - n m B m,n ε m + ε n 1 m 1 T n - n m B m,n ε m + ε n 1 n 1 T m = -2B, ( 98
)
where B is the matrix with elements B m,n = B m,n /(ε m + ε n ) or, more explicitly, B ia, jb = B ia, jb /(ε a + ε bε iε j ). Therefore, we have

Q RPAx α = I -2αB + O(α 2 ), (99)
and, similarly, the first-order expansion of the inverse matrix (Q RPAx α ) -1 = n (x α,ny α,n )(x α,n -

y α,n ) T yields Q RPAx α -1 = I + 2αB + O(α 2 ). ( 100
)
99 and 100 show that the approximation Q α + Q -1 α ≈ 2I, which lead to the definitions of E RPAx-IIa c [49] and E dRPA-IIa c [59], is correct up to first order in α.

All the above considerations remain valid for the dRPA case, except for the replacements A ′ → K and B → K, with the obvious results

Q dRPA α = I -2αK + O(α 2 ), (101)
and

Q dRPA α -1 = I + 2αK + O(α 2 ), (102)
where the matrix elements of K are given by K m,n = K m,n /(ε m + ε n ) or, more explicitly, K ia, jb = K ia, jb /(ε a + ε bε iε j ).

We can give now the second-order limits of the RPA correlation energy variants. Using 101, we find the second-order limit of the dRPA correlation energy variant of 28

E dRPA-I c ≈ 1 2 1 0 dα tr -2αK K = - 1 2 tr K K , (103)
which is not the normal MP2 correlation energy, but a MP2-like correlation energy without exchange, also called direct MP2 or JMP2. 50 In a similar way, 101 and 100 give the second-order limit of the RPAx-II correlation energy variant of 48, which is the same for its approximation of 49,

E RPAx-II c ≈ E RPAx-IIa c ≈ 1 4 1 0 dα tr -2αB B = - 1 4 tr B B , (104)
which is exactly the MP2 correlation energy expression (except for the possible replacement of Hartree-Fock orbitals and orbital energies with corresponding Kohn-Sham quantities). The secondorder limit of the dRPA-II correlation energy variant of 58 and its approximation of 59 are found with 101 and 102

E dRPA-II c ≈ E dRPA-IIa c ≈ 1 2 1 0 dα tr -2αK B = - 1 2 tr K B . ( 105
)
Using the antisymmetry of B and observing the prefactor of 1/2, it can easily be seen that this is another way to write the usual MP2 correlation energy expression. Finally, the RPAx-I correlation energy variant of 61 has the following second-order limit

E RPAx-I c ≈ 1 2 1 0 dα tr -2αB K = - 1 2 tr B K , (106)
which again exactly corresponds to the usual MP2 correlation energy expression.

Let us now consider the case of a closed-shell system. In this case, there is (at least) a fourfold degeneracy in the ε ε ε block of Λ 0 since ε i↑ = ε i↓ and ε a↑ = ε a↓ . As a consequence, the condition of non-degeneracy of zeroth-order excitation energies ω 0,n = ǫ ia leading to 91 and 94a is violated.

Even if the final results for the second-order correlation energies do not contain differences of excitation energies anymore, a different derivation is needed. This may be achieved by first spinadapting the RPA response equation (for the details, see, e.g., Ref. 33), and only subsequently making the perturbation expansion on the spin-adapted energy expressions of 5. Assuming the absence of further degeneracies between orbital energy differences (zeroth-order excitation energies), one obtains formally identical expansions for the singlet and triplet blocks. For example, the spin-adapted matrices 1 Q α = n ( 1 x α,n + 1 y α,n )( 1 x α,n + 1 y α,n ) T and 3 Q α = n ( 3 x α,n + 3 y α,n )( 3 x α,n + 3 y α,n ) T , where ( 1 x α,n , 1 y α,n ) and ( 3 x α,n , 3 y α,n ) are the singlet and triplet eigenvectors, and the corresponding inverse matrices ( 1 Q α ) -1 and ( 3 Q α ) -1 have the following expansions in the case of

  basis of occupied-virtual orbital products, i.e. pq = ia or ai and rs = jb or b j where i, j refer to occupied orbitals and a, b to virtual orbitals. Neither P dRPA c,α nor P RPAx c,α are properly antisymmetric. As a consequence, the two correlation energy expressions, 1 and 2, are no longer equivalent in dRPA or RPAx. This leads to at least four RPA variants for calculating correlation energies, denoted here by dRPA-I, dRPA-II, RPAx-I, and RPAx-II, depending whether the correlation density matrix is contracted with the non-antisymmetrized two-electron integrals V (variants I) or the antisymmetrized two-electron integrals W (variants II).

c)

  pq,rs → (P dRPA c ) pq,rs -(P dRPA c

Recently, 40 1 are 1 ) 1 / 2 }

 401112 have been used by Eshuis et al.32 as the starting point for developing a computationally efficient algorithm for calculating the dRPA-I correlation energy. Note that expression 40 could have also been found by noting that the eigenvalues of M dRPA . However, working with Q -1 α will be useful for the other variants. Also, comparison of 28 and 38 provides us with a decomposition of the correlation energy into kinetic and potential contributions, E dRPA-I c = T dRPA-I c + U dRPA-I c

1 and 3 M 1 ) 1 / 2 } 1 ) 1 / 2

 13112112 RPAx 1 have the same eigen-values and thus tr{( 3 N RPAx = tr{( 3 M RPAx }.

Figure 1 :

 1 Figure1: Ratios between various RPA correlation energy variants and the FCI-quality correlation energy as estimated by Davidson and coworkers,60,61 with and without range separation. All the correlation energies have been extrapolated to the CBS limit. The RPA correlation energies E * c (RPA) are redefined here as the difference between the total RPA energies and the regular HF energies.

2Figure 2 :

 2 Figure 2: Errors in the equilibrium bond lengths and harmonic vibrational frequencies for simple diatomic molecules, calculated with the full-range and range-separated RPA variants and compared to experimental reference values. All the correlation energies have been extrapolated to the CBS limit. The experimental reference values are (in bohr and cm -1 ) H 2 R e =1.40112, ω e = 4401.21; HF R e = 1.73250, ω e = 4138.32; N 2 R e = 2.07431, ω e = 2358.57. 66

Figure 3 :

 3 Figure 3: Interaction energy curves of He 2 , Ne 2 , and Ar 2 , calculated with the full-range and rangeseparated RPA variants. All the correlation energies have been extrapolated to the CBS limit.

RPAx 1 , 3 Q

 13 RPAx α ±1 = I ∓ 2α 1,3 B + O(α 2 ). (107) with 1 B m,n = 1 B m,n /(ε m + ε n ) and 3 B m,n = 3 B m,n /(ε m + ε n ).Using these results, one can easily check that all the spin-adapted correlation expressions of 5 correctly reduce to MP2 at second order, except for the dRPA-I variant which reduces to direct MP2.

  Tang and Toennies, 67 and on the other hand the repulsive (exponential) component of the same potential. These latter curves serve as useful guides to estimate the accuracy of the single-determinant reference energies, i.e. EXX energies with PBE orbitals or RSH energies. It is quite clear that the quality of the results depends strongly on the quality of the repulsive curve. The poorest interaction energy curves are obtained for the

		60					Accurate EXX(PBE)			30					Accurate RSH
							PBE + dRPA-I							RSH + dRPA-I
		40					PBE + dRPA-II PBE + dRPA-IIa		20					RSH + dRPA-II RSH + dRPA-IIa
							PBE + RPAx-I							RSH + RPAx-I
		20					PBE + RPAx-II		10					RSH + RPAx-II
	E int [ µhartree]	-20 0								E int [ µhartree]	-10 0
		-40									-20
		-60									-30
		5	6		7	8	9	10	11	5	6	7		8	9		10	11
					R [bohr]								R [bohr]
		(a) He without range separation			(b) He with range separation
							Accurate								Accurate
		200					EXX(PBE) PBE + dRPA-I		100					RSH RSH + dRPA-I
							PBE + dRPA-II							RSH + dRPA-II
							PBE + dRPA-IIa							RSH + dRPA-IIa
		100					PBE + RPAx-I PBE + RPAx-II		50					RSH + RPAx-I RSH + RPAx-II
	E int [ µhartree]	0								E int [ µhartree]	0
		-100									-50
		-200									-100
		5	6	7	8	9	10	11 12	13	5	6	7	8	9	10	11 12	13
					R [bohr]								R [bohr]
		(c) Ne without range separation			(d) Ne with range separation
		600									600
							Accurate								Accurate
							EXX(PBE)								RSH
		400					PBE + dRPA-I		400					RSH + dRPA-I
							PBE + dRPA-II							RSH + dRPA-II
							PBE + dRPA-IIa							RSH + dRPA-IIa
		200					PBE + RPAx-I PBE + RPAx-II		200					RSH + RPAx-I RSH + RPAx-II
	E int [ µhartree]	0								E int [ µhartree]	0
		-200									-200
		-400									-400
		-600									-600
		6	7	8	9	10	11	12	13	14	6	7	8	9	10	11	12	13	14
					R [bohr]								R [bohr]
		(e) Ar without range separation			(f) Ar with range separation

Acknowledgments

We thank A. Savin (Paris) for discussions. This work was supported by ANR (Agence National de Recherche) via contract number ANR-07-BLAN-0272 (Wademecom). G. J. thanks Nancy University for a visiting professorship during which part of this work has been carried out.

A Second-order approximations to the RPA correlation energy expressions In this appendix, we explicitly derive the approximations at second order in the electron-electron interaction of the RPA correlation energy variants.

We will deal with the more general RPAx response equation and obtain dRPA as a special case.

We thus start from the response equation

RPA of the density-functional/material-science community and neglects all exchange interactions. The dRPA-II variant uses a density matrix without exchange but contracted with antisymmetrized two-electron integrals. It is original to this work, although it resembles the SOSEX method, 23 especially in its approximate form named dRPA-IIa. The RPAx-I uses a density matrix with exchange but contracted with non-antisymmetrized two-electron integrals. It has previously been discussed in the context of range separated density-functional theory. 15,33 The RPAx-II variant uses a density matrix with exchange and contracted with antisymmetrized two-electron integrals. The RPAx-II method itself is obviously not new, 49 but we have derived several new expressions for it. Contracting the density matrix with either non-antisymmetrized or antisymmetrized two-electron integrals is not equivalent because of the breaking of the antisymmetry of the density matrix in RPA. For the dRPA-I and RPAx-II variants, we have made the connection with the plasmon formulation, and clarify the origin of the factor of 1/4 in the plasmon formula for RPAx-II instead of the factor of 1/2 for dRPA-I. We have carefully studied the second-order limit in the electron-electron interaction, and showed that all the correlation energy variants except for dRPA-I correctly reduce to the MP2 correlation energy (see Appendix). Finally, we have derived the spin-adapted forms of all these methods for closed-shell systems, and implemented and tested them with and without range separation of the electron-electron interaction.

The numerical examples on atomic and molecular systems show that the RPAx variants without range separation frequently suffer from instabilities in the RPAx response equation which make it impossible to extract a meaningful correlation energy in these cases. However, no instabilities are encountered with range separation, and the RPAx variants can be thus safely applied. The tests performed do not allow us to identify an RPA variant which would be uniformly better than the others. Without range-separation, dRPA-II performs well for atomic correlation energies and equilibrium molecular properties, but significantly underestimates London dispersion interaction energies for which RPAx-I is more accurate. With range separation, all the RPA variants tend to give more accurate results, and they also become much more similar to each other. Rangeseparated RPAx-II appears as the best variant for equilibrium molecular properties and range-