A Equivalence of the density matrix and logarithmic expressions

Dielectric matrix formulation of correlation energies in the Random Phase Approximation (RPA): inclusion of exchange effects

Abstract

Starting from the general expression for the ground state correlation energy in the adiabatic connection fluctuation dissipation theorem (ACFDT) framework, it is shown that the dielectric matrix formulation, which is usually applied to calculate the direct random phase approximation (dRPA) correlation energy, can be used for alternative RPA expressions including exchange effects. Within this famework, the ACFDT analog of the second order screened exchange (SOSEX) approximation leads to a logarithmic formula for the correlation energy similar to the direct RPA expression. Alternatively, the contribution of the exchange can be included in the kernel used to evaluate the response functions. In this case the use of an approximate kernel is crucial to simplify the formalism and to obtain a correlation energy in logarithmic form. Technical details of the implementation of these methods are discussed and it is shown that one can take advantage of density fitting or Cholesky decomposition techniques to improve the computational efficiency; a discussion on the numerical quadrature made on the frequency variable is also provided. A series of test calculations on atomic correlation energies and molecular reaction energies shows that exchange effects are instrumental to improve over direct RPA results.

I Introduction

In the past few years, random phase approximation (RPA) approaches, which belong to the methods on the 5th rung of Jacob’s ladder Perdew and Schmidt (2001), have been on the way of becoming a practical tool to construct correlation energy functionals Furche (2001); Fuchs and Gonze (2002); Miyake et al. (2002); Furche and van Voorhis (2005); Fuchs et al. (2005); Furche (2008); Harl and Kresse (2008); Lu et al. (2009); Li et al. (2010); Grüneis et al. (2009); Nguyen and de Gironcoli (2009); Nguyen and Galli (2010); Paier et al. (2010); Hesselmann and Görling (2010, 2011); Eshuis et al. (2012); Hesselmann (2012); Ren et al. (2012); Toulouse et al. (2009); Janesko et al. (2009a, b); Rocca (2014). Specifically, within the framework of the adiabatic-connection fluctuation-dissipation-theorem (ACFDT), the electronic correlation energy is expressed in terms of dynamical (frequency-dependent) linear response functions and the electron-electron interaction is gradually switched on from the independent particle reference to the fully interacting state using an adiabatic connection parameter. Therefore the correlation energy takes the form of a multiple integral involving both the frequency and the adiabatic connection parameters. Instead of the exact linear response function, convenient approximations, like the RPA, are used.

Among the numerous possible ways to express the RPA ground state correlation energy, depending on the order in which the analytical and/or numerical frequency- and interaction strength-integrations are performed, one may cite (1) the density matrix formulation Furche (2001), which leads to expressions involving numerical integration over the interaction strength and (2) the dielectric matrix formulation which involves a numerical integral over frequency of a logarithmic expression involving the dynamical dielectric function. This approach has mainly been used by solid state theorists Miyake et al. (2002); Harl and Kresse (2009); Lu et al. (2009); Nguyen and de Gironcoli (2009) and recently adapted also in a density fitting framework by Furche et al. Eshuis et al. (2010). Within this formulation the direct calculation of the lowest eigenvalues and eigenvectors of the dielectric matrix Wilson et al. (2008) can be used to achieve a more compact representation of dielectric functions Lu et al. (2009); Nguyen and de Gironcoli (2009); Nguyen and Galli (2010); Li et al. (2010); Rocca (2014); Kaoui and Rocca (2016). Finally, there are approaches which avoid numerical integration altogether, like (3) the plasmon formula, obtained after a double analytical integration on both the frequency and the interaction strength Furche (2008). An elegant way to obtain the plasmon expression consists in solving the algebraic Riccati-equations of the RPA problem Sanderson (1965); Szabo and Ostlund (1977). This method has been shown to be strictly equivalent to a coupled cluster doubles approach in the ring-approximation (rCCD) Toulouse et al. (2011); Scuseria et al. (2008). Modifications and approximations to the rCCD equations and energy expression has led to a whole class of additional RPA-based approaches for the ground state correlation energy Szabo and Ostlund (1977); Lotrich and Bartlett (2011); Klopper et al. (2011); Hesselmann (2012).

The RPA problem can be set up either with the inclusion of nonlocal exchange effects (leading to a class of approximations denoted RPAx) or by restricting the coupling between independent particle responses to the direct Hartree interaction (leading to the class denoted dRPA). Additionally, according to the classification and nomenclature proposed in Ref. Ángyán et al., 2011, RPA methods for the ground-state correlation energy can be sought in two main flavors: either by the complete neglect of the exchange integrals, i.e. by taking the contraction of the dRPA response function with non-antisymmetrized two-electron integrals (dRPA-I), or by a full consideration of non-local exchange, when antisymmetrized two-electron integrals are contracted with the RPAx response function matrix elements (RPAx-II). This latter approach has been followed by some early works in quantum chemistry, based on Hartree-Fock orbitals Oddershede (1978); Szabo and Ostlund (1977); Ostlund and Karplus (1971); McLachlan and Ball (1964), while dRPA is usually the method of choice in the density functional context, starting from Kohn-Sham orbitals. Partial inclusion/omission of nonlocal exchange leads to “mixed” methodologies, like dRPA-II and RPAx-I.

It has been suggested that certain shortcomings of the dRPA correlation energy can be remedied by including nonlocal exchange interactions in a perturbative way, i.e. with the dRPA polarization propagator being contracted with a list of fully antisymmetrized two-electron integrals Grüneis et al. (2009); Paier et al. (2010). We can mention, in this aspect, the SOSEX (second order screened exchange) corrections, which have been formulated originally within the rCCD formalism. It has been pointed out that “it is difficult to motivate this approximation in the framework of ACFDT” Grüneis et al. (2009). This situation seems to be somewhat paradoxical, since the plasmon formula, identical to the rCCD correlation energy expression, can be derived from ACFDT in a straightforward manner Furche (2008), and there is no fundamental reason to think that an analogous derivation is impossible for SOSEX. In fact, it has been shown that a very similar perturbative screened exchange formula, which has been designated by the acronym dRPA-IIa, can be obtained within the density matrix formulation of RPA Ángyán et al. (2011). Although this expression gives correlation energies numerically very close to the rCCD-based SOSEX, it has been proven that they are not strictly identical Jansen et al. (2010). Recently, in a similar but different fashion, the AXK introduced by Bates et al. reduces the self-interaction error and improves the description of static correlation over dRPA.

In this work we discuss different approximations to efficiently include exchange effects within the dielectric matrix formulation of the RPA correlation energy. In a quantum chemical context these methodologies provide an alternative to the ring-CCD-based RPA formalism. Additionally, these methods are of significant interest for the solid state physics community, where the dielectric matrix approach is almost exclusively used. As it will be shown in the following, the dRPA-IIa approximation can be derived also in the dielectric matrix formulation of RPA, leading to an alternative algorithm to calculate the SOSEX-like dRPA-IIa energy. At the same time we have access to the conventional MP2 energy, which corresponds to the lowest non-vanishing term in the series expansion of the dRPA-IIa correlation energy. We note that approximations in the RPAx class lead to other correlation energy estimates, which can be also adapted to the logarithmic formulation. Specifically, the RPAx-Ia approximation is derived, which includes only the particle-hole contribution to the exchange kernel. It will be shown that this last approach provides the most promising results in the numerical tests considered in this work.

In Sec. II.1, the working equations for the well-known dielectric matrix formulation of the direct RPA method will be derived from the ACFDT. This version is particularly simple because of the complete neglect of the exchange effects. Incorporation of exchange leads to a more complicated expression and one is constrained to apply various additional approximations to obtain practical expressions. As discussed in Sec. II.2, one possibility is a SOSEX-like correlation energy expression, which is the dielectric matrix formulation of the dRPA-IIa variant. In Sec. II.3 we will discuss an additional variant, which includes exchange effects in the response function. The possible practical implementations of these methodologies are discussed in some details in Sec. III. Several atomic and molecular systems will be used for numerical illustration of the formalism, in Sec. IV. Finally, Sec. V contains our conclusions.

Ii RPA versions from ACFDT and their dielectric matrix formulation

In the adiabatic-connection fluctuation-dissipation theorem (ACFDT) approach the correlation energy reads as Langreth and Perdew (1975, 1977):

(1)

where is the matrix representation of the interaction which contains two-electron integrals (see Eqs. (7) and (16)). In Eq. 1 the four-index matrix representation of the dynamical polarization propagator (also called density matrix response function) at interaction strength is obtained from a Dyson-like equation,

(2)

where is the interaction kernel matrix (see Eqs. (7) and (24)). In the two previous equations, the polarization propagator of the non-interacting reference system, , is given by:

(3)

with

(4)

In the random phase approximation (RPA), the elements of the matrix are the independent one-particle excitation energies: , where , is the energy of a virtual orbital and the energy of an occupied orbital. Here and in the following we assume that a finite-dimensional basis set is used for the representation of occupied and virtual orbitals.

The polarization propagator of the non-interacting reference system Oddershede (1978) then reads:

(5)

which introduces compact notations for the diagonal blocks, and .

The dimensions of the matrices appearing in Eq. (1) are , where is the number of the products between occupied and virtual orbitals. However, since we need only the trace for the correlation energy, the same result can be obtained using matrices with a dimension reduced to , as it will be shown below.

In several earlier works Harl (2008); Hesselmann and Görling (2010) the derivation of the logarithmic energy expression have been based on a Taylor expansion of matrix functions (namely the inverse and the logarithmic functions). Such a procedure raises some problems of general validity, since the conditions for an absolute convergence of the series expansions cannot be always satisfied. Below, we propose an alternative derivation using matrix functions, which requires only the weaker condition that the matrix be diagonalizable. The general definition of a matrix function is

(6)

where is the eigendecomposition of the matrix . Hence, all along the manuscript, the matrices are diagonal matrices of eigenvalues and are eigenvectors.

ii.1 Direct RPA

In the direct RPA (dRPA), the nonlocal Hartree-Fock exchange is not taken into account in the interaction kernel of the Dyson-like equation in Eq. (2). Furthermore, we will first consider an interaction matrix constituted of non-antisymmetrized two-electron integrals in Eq. (1). This yields:

(7)

where are the non-antisymmetrized two-electron integrals (physicist’s notation) in spin-orbitals. In this situation, one only has to deal with the matrix product (see Eqs. (1) and  (2)). Considering the block-structure of ,

(8)

application of the following unitary transformation:

(9)

to the integrand of Eq. (1) yields:

(10)

where we use the notation . The notation emphasizes that we have reduced the matrix dimensions as compared to (the same applies later to / and to /). Defining the function that satisfies the dimension-reduced Dyson-like equation,

(11)

we retrieve the analog of Eq. (1) with dimension-reduced matrices

(12)

Equation (12) can be further transformed either by an analytical frequency-integration, which leads to the density matrix expression of the dRPA correlation energy (this has already been discussed in a previous publication Ángyán et al. (2011), and is briefly recalled in Appendix A), or by an analytical integration over the adiabatic connection parameter. This will be done by considering the eigenvalue decomposition of . Using the dimension-reduced Dyson-like equation (Eq. 11), we express seen in Eq. (12) as a matrix function of :

(13)

where, as stated before, is the diagonal matrix of the eigenvalues of and contains its eigenvectors. Using the cyclic invariance of the trace and denoting by the -th diagonal element of , Eq. (12) becomes

(14)

We recognize here the matrix representation of the dielectric function, , and this expression may be called the dielectric matrix formulation of the dRPA-I correlation energy.

Although it might not be obvious at first sight, a further analytical integration on the frequency leads directly to the plasmon formula. Such a relationship has already been mentioned by McLachlan et al. in the sixties McLachlan et al. (1963). Inversely, the dielectric matrix formulation can be derived from the plasmon formula, as it has been shown recently by Eshuis et al. Eshuis et al. (2010). In Appendix B, we present an alternative derivation of the plasmon formula for the dRPA-I correlation energy from the dielectric formulation of Eq. (12). This seals the strict equivalence of all formalisms in the case of dRPA-I. Notice that the plasmon formula applies only to the dRPA-I (and to the RPAx-II, which is not discussed here) energy expressions.

Moreover, doing a series expansion of the matrix logarithm in Eq. (14) and retaining only the first non-vanishing term, we obtain:

(15)

which corresponds to a direct-MP2 method, that is to say a version of the MP2 energy without exchange terms Furche (2001) (in its range-separated variant, JMP2 Janesko and Scuseria (2009)).

ii.2 Rpa+sosex

Considering the fact that, in the ACFDT expression, the polarization propagator carries information about the screening properties in the system, a SOSEX-like Ansatz consists in replacing the direct interaction matrix by the antisymmetrized expression ,

(16)

where the antisymmetrized two-electron integral matrices have the elements and . This particular form of the interaction matrix can be derived from the time dependent Hartree-Fock equations using the fact that the independent particle wave function, providing the bare, unscreened response, satisfies the Brillouin theorem. We note that one can argue for other choices of the matrix (vide infra).

The resulting correlation energy expression corresponds to the dRPA-II variant Ángyán et al. (2011),

(17)

which can be designated also as screened exchange dRPA (dRPA-SX), since the dRPA response function screens the full, Coulomb plus exchange interaction represented by the matrix .

As a consequence of the more complicated block structure of the interaction matrix, it is less obvious to reduce the dimensions of the problem in the dRPA-II correlation energy as compared to the dRPA-I case. It is reasonable to decompose into a “major” contribution, involving four identical blocks and a “minor” correction, which consists in the difference of the diagonal blocks:

(18)

The major contribution to the integrand can be brought to a convenient form by applying the unitary transformation (Eq. 9) to the matrix products under the trace in Eq. (17). After evaluating the inverse of the blocked matrix and subsequent matrix multiplications, one gets the trace of the major contribution in dimension-reduced matrices. The minor contribution to the energy expression of Eq. (17) can be neglected, as is shown in Sec. 2 of the Supporting Information, and this leads to the following approximation of the screened exchange dRPA, denoted in our previous publication Ángyán et al. (2011) as dRPA-IIa:

(19)

As in the dRPA-I case, after analytical frequency integration we get the same density matrix formulation expression as the one we have obtained in a quite different manner in Ref. Ángyán et al., 2011.

Before exploring another alternative, which consists in an analytical integration according to the interaction strength parameter, the reader must bear in mind the close analogy of dRPA-IIa and the so-called SOSEX (second order screened exchange) approximation Grüneis et al. (2009); Freeman (1977), usually defined in the framework of a drCCD (direct ring coupled cluster doubles) theory. Although the above ACFDT-based expression is not strictly equal to the drCCD-based SOSEX, their difference appears only at the third order of perturbation, as was demonstrated in Ref. Jansen et al., 2010. The difference of those two variants has been found numerically small for all of the systems studied up to now Ángyán et al. (2011). The approximation, which leads from dRPA-II to dRPA-IIa shows in a self-explanatory manner the perturbation character of the SOSEX (dRPA-IIa), i.e. “second-order screened exchange” compared to the fully screened exchange (SX, i.e. dRPA-II) version.

In order to perform the analytical integration along the adiabatic connection path, we take advantage of the matrix function formalism (with now obvious notations for and ):

(20)

Particular care has to be taken to define the inverse of in Eq. (II.2), since this matrix might be singular or close to singular in certain representations. Since the (close to) zero eigenvalues do not contribute to the sum over in Eq. (II.2), as can be seen from an expansion of the logarithm, it is implicit that is defined in terms of a pseudoinverse where for finite and for (or below a certain small threshold).

By decomposing the antisymmetrized two-electron integrals as , where the matrix has the elements , the correlation energy can be separated to a dRPA-I contribution and a SOSEX correction:

(21)
(22)

Furthermore, it is easy to verify that the second order approximation to the dRPA-IIa correlation energy is exactly the usual MP2 correlation energy (again by doing a series expansion of the matrix logarithm):

(23)

An identical result has been obtained in the Appendix of Ref. Ángyán et al., 2011. It is important to mention that in this context the MP2 correlation energy represents the first non-vanishing contribution to a converging series only if all the eigenvalues of are smaller than 1.

ii.3 Approximate RPA with exchange

The RPAx variants are based on a response function which fully takes into account the nonlocal exchange effects by a Hartree-Fock type kernel. In other words, the response function satisfies the following Dyson-like equation, where in Eq. (2):

(24)

In the RPAx-I case, that is to say with in Eq. (1), the ACFDT expression becomes:

(25)

Let us write again as a sum of two terms,

(26)

and use the transformation on the matrix product under the trace. Neglecting the minor correction contribution leads to the RPAx-Ia energy:

(27)

where we define which satisfies the dimension-reduced Dyson-like equation, obtained from Eq. (24) after taking an approximate kernel including only the major contribution from Eq. (26):

(28)

The notation RPAX (with capital X) refers to an analogy to an approximate RPA variant with exchange, proposed by Hesselmann Hesselmann (2012). The relationships to Hesselmann’s approach will be discussed elsewhere. Comparison to Eqs. (19) and (11) shows that, formally, the roles of the matrices and is merely exchanged with respect to the dRPA-IIa energy expression. Hence, the corresponding logarithmic formula obtained after -integration is:

(29)

again defining in terms of a pseudoinverse involving the eigenvalues of .

The second order approximation yields the MP2 energy in the exact same manner as the dRPA-IIa case (consider the exchange of the matrices and in Eq. (II.2)).

As shown in Sec. IV, the approximate RPAx-Ia correlation energy expression leads to the most accurate results for the systems considered in this work. Additionally, while the RPAx-I approximation (Eq. II.3) is known to suffer from instabilities when the initial states are approximated from semi-local functionals Ángyán et al. (2011), the RPAx-Ia is numerically stable.

Iii Computational realization

For the practical implementation, it is worth noting that spin-adaptation of all the equations is trivial: in spacial orbitals, the two-electron integrals read as , , and .

In the following, we present an orbital-based implementation of the equations derived in Sec. II. An implementation using density fitting is shown in the Appendix. C.

iii.1 Orbital-based implementation

The computational realization of equations (14), (II.2) and (29) may proceed following a common scheme. In fact, Eq. (14) can be considered as a special case of Eq. (II.2), where , and Eq. (29) can be obtained by interchanging the roles of and in Eq. (II.2). Therefore we focus our attention to the case of Eq. (II.2), i.e. the dRPA-IIa energy expression.

We can rewrite Eq. (II.2) in terms of the symmetric matrices and as

(30)

where we took advantage of the cyclic invariance of the trace and used the symmetry of with respect to the replacement of by to restrain the integral to the positive imaginary axis. The derivation of Eq. (30) also requires the following identity for a generic function of a matrix :

(31)

that can be easily derived from the matrix function formalism and from the fact that implies (the matrix defines the unitary transformation that diagonalizes the Hermitian matrix ).

By expressing Eq. (30) in terms of the (diagonal) eigenvalue matrix we obtain

(32)

with . It should be emphasized that the matrix has been transformed with the eigenvectors of and therefore is not diagonal. Nevertheless in order to calculate the trace of its product with the diagonal matrix we need only its diagonal elements, which will be designated by the shorthand notation .

The presence of potentially very small eigenvalues may lead to numerical instabilities. This can be avoided by using the following power series expansion of the logarithm under a given threshold of the eigenvalue (in practice, the summation is carried up to and the threshold is 0.0001):

(33)

At second order, , we obtain the MP2 energy, which can be designated as “Casimir-Polder transform MP2”

(34)

and which offers an interesting alternative to calculate conventional MP2 energies, especially in solids Marsman et al. (2009). Higher order terms of the series expansion do not correspond exactly to the MP3 and MP4 energies. A vague analogy can be noted between this expression and the Laplace-transform method to obtain the MP2 energy Eshuis et al. (2010).

iii.2 Numerical frequency integration

In principle, the numerical frequency integration is expected to be “fairly unproblematic” Harl et al. (2010), since the integrand is expected to have a smoothly decaying behavior. While it seems to be really the case for solids, where a mapping of the Gauss-Legendre quadrature to the interval (truncated at about 30 a.u.) with an exponentially decaying weighting function ensures a reasonable convergence (0.05 mH) with only 16 quadrature points Harl et al. (2010), an accuracy of 0.2 mH (about 0.15 kcal/mol) claimed by the same authors is clearly insufficient for atoms and molecules. The situation for atomic and molecular systems is not better for more sophisticated weighting function models either.

Lu et al. have performed the frequency integration by a 10-point Gauss-Legendre quadrature in the range of , where , with  a.u. This quadrature was found rather insensitive to the choice of the parameter Lu et al. (2009).

We have implemented the method proposed by Eshuis, Yarkony and Furche Eshuis et al. (2010), based on the Clenshaw-Curtis quadrature Boyd (1987), which consists in writing the integral of a function as

(35)

where with and . This quadrature scheme has a single scaling parameter , which can be adjusted to each individual system by requiring that the following analytically solvable expression, based on a diagonal model of the dRPA dielectric matrix formulation of Eq. (14) (remember that ) :

(36)

be reproduced at best by the numerical integral,

(37)

This tuning method of the parameter can be extended for the dRPA-IIa case, using the analytically solvable dRPA-IIa diagonal approximation:

(38)

Note that the diagonal approximation to the dRPA-IIa correlation energy is half the diagonal approximation to the dRPA-I correlation energy seen in Eq. III.2 (remember that in spatial orbitals we have and ) and there is no need to re-optimize the free parameter of the numerical quadrature.

Iv Numerical results

We have implemented the dielectric matrix based dRPA-I (Eq. 14), dRPA-IIa (Eq. II.2), and RPAx-Ia (Eq. 29) energy formulae in an occupied-virtual basis set representation within the development version of the MOLPRO quantum chemistry package Werner et al. (2012a, b). Although this algorithm is expected to be less efficient than the density fitting approach (in particular for larger systems), we considered this implementation useful to produce benchmark results exempt of density fitting uncertainties.

Below we present electron correlation energies for some atoms and ions Davidson et al. (1991); Chakravorty et al. (1993) and a test set of reaction energies Hesselmann (2012). The correlation energy from the different RPA approximations has been calculated starting from a self-consistent DFT calculation based on the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional. The total RPA energy is then evaluated from the expression

(39)

where is the Hartree-Fock energy computed from PBE orbitals.

For comparison purposes, we also provide results obtained from the rCCD-based SOSEX approximation Grüneis et al. (2009); Freeman (1977) and the dRPA-II approximation as derived within the adiabatic connection approach Ángyán et al. (2011). As already discussed in Sec. I and Sec. II.2 the dRPA-IIa (Eq. II.2) and the rCCD-based SOSEX are analogous although not strictly equivalent; the numerical calculations below help to better quantify this statement. The dRPA-II approximation as derived within the adiabatic connection approach is equivalent to the dynamical polarizability expression in Eq. 17, that, however, cannot be conveniently implemented (both frequency and coupling constant integration are necessary and the dimensionality of the polarizability cannot be reduced). In this case the comparison between the dielectric matrix-based dRPA-IIa and the adiabatic connection-based dRPA-II is useful to understand the effect of the approximation of the kernel introduced by discarding the minor contribution to Eq. 18. To avoid confusion the dRPA-I, dRPA-IIa, and RPAx-Ia approximations derived within the dielectric matrix formulation will be denoted with an additional “DIEL”, the SOSEX approximation with an additional “rCCD”, and the dRPA-II approximation within the adiabatic connection approach with an additional “AC”.

iv.1 Frequency quadrature

Figure 1: Convergence with respect to the number of quadrature points of the Gauss-Chebyshev, Gauss-Legendre and Clenshaw-Curtis schemes for Li (up) and Mg (down). The dRPA-I DIEL energy is compared to the dRPA-I PLASMON energy, and the dRPA-IIa DIEL and RPAx-Ia DIEL are compared against their AC analog.

The frequency integrations in Eqs. 14II.2 and  29 are carried out by a quadrature. We tested the convergence of the Gauss-ChebyshevRijks and Wormer (1988), Gauss-LegendreFurche (2001) and Clenshaw-Curtis (described in Sec. III.2) schemes with respect to the number of quadrature points, for calculations on several atoms and ions. The RPA correlation energies were computed using the aug-cc-pCVXZ basis sets (X=6 for He, X=5 for Ne, Ar, B, and Al, and X=Q for Li, Na, Be, and Mg). The results of the convergence study can be found in Table 1 and Figure 1, where we use as reference the dRPA-I PLASMON (for which both the frequency and coupling constant integrations are analytical, see Appendix B), and the dRPA-IIa AC and RPAx-Ia AC energies (for which no PLASMON analog exists). We see that the Gauss-Chebyshev and Gauss-Legendre schemes yield, in this case, unsatisfying results: the Gauss-Chebyshev scheme has a slow convergence with respect to the number of quadrature points in all the cases studied, and while the Gauss-Legendre quadrature performs well for Be, He and Li, it exhibits a slow convergence in the cases of B, Mg and Na. On the other hand, the Clenshaw-Curtis quadrature was found to converge rapidly for all studied atoms and ions. As a result, we choose in the following to perform the frequency integration with a 48-points Clenshaw-Curtis quadrature.

B+ Mg Na+ Be He Li+
Gauss-Chebyshev
24 1.51 0.07 0.90 1.63 0.51 1.81
48 1.09 1.42 1.05 0.84 0.16 0.51
72 0.78 1.45 0.81 0.44 0.04 0.15
96 0.58 1.16 0.66 0.22 0.01 0.05
Gauss-Legendre
24 -2.22 -3.38 -2.23 -1.06 -0.13 -0.46
48 -0.78 -1.17 -0.78 -0.23 -0.02 -0.09
72 -0.37 -0.52 -0.35 -0.09 -0.01 -0.03
96 -0.21 -0.27 -0.19 -0.04 0.00 -0.01
Clenshaw-Curtis
24 -0.35 -0.77 0.02 -0.23 0.00 0.00
48 -0.02 -0.04 0.00 -0.01 0.00 0.00
72 0.00 0.00 0.00 0.00 0.00 0.00
96 0.00 0.00 0.00 0.00 0.00 0.00
Table 1: Errors made in the calculation of the dRPA-I DIEL energy using 24-, 48-, 72- and 96-points Gauss-Chebyshev, Gauss-Legendre and Clenshaw-Curtis quadrature schemes (in percentage with respect to the dRPA-I PLASMON result).

iv.2 Atomic correlation energies

Figure 2: Percentage error of total correlation energies compared to the FCI-quality estimates by Davidson and his coworkers Davidson et al. (1991); Chakravorty et al. (1993) for simple closed shell atoms and ions calculated with different RPA variants. Correlation energies have been extrapolated to the CBS limit.

As a first test, we applied the dielectric matrix formalism to compute correlation energies for several atoms and ions. To verify the accuracy of the RPA-based approximations we compared the obtained results with the full configuration interaction (FCI) quality correlation energy estimates by Davidson and collaborators Davidson et al. (1991); Chakravorty et al. (1993). For a meaningful comparison it is necessary to keep into account that the FCI-quality correlation energies have been obtained with respect to a Hartree-Fock reference point. For this reason we redefined the RPA correlation energies as the difference of RPA total energies and regular Hartree-Fock energies. This procedure was already used in Ref. Ángyán et al., 2011. Core excitations have been included in the calculations. The RPA-based correlation energies were computed in the complete basis set (CBS) limit by the usual 1/X formula Kutzelnigg and Morgan III (1992) considering aug-cc-pCVXZ basis sets up to X=6 for He, up to X=5 for Ne, Ar, B, and Al, and up to X=Q for Li, Na, Be, and Mg. The ratios between the RPA-based correlation energies and reference FCI correlation energies are shown in Fig. 2. As expected the dRPA-I DIEL approximation strongly overestimates the absolute value of the correlation energy (by 50 to 100). On the other hand, the dRPA-IIa DIEL and RPAx-Ia DIEL approximations, which include exchange contributions, significantly improve over the direct RPA results and lead to a percentage of the correlation energy close to . On the scale of Fig. 2 the dRPA-IIa DIEL and SOSEX rCCD results are basically indistinguishable (the largest difference is found for Al, where SOSEX rCCD gives 95.79 of the FCI correlation energy and dRPA-IIa 95.57). The dRPA-II AC, which does not rely on the neglect of the minor contribution to the kernel discussed in Eq. 18, provides the worst results among the approximations including exchange effects. Additionally, in order to quantify the one-electron self-interaction error, we computed correlation energies for the hydrogen atom. Also in this case the exchange contribution significantly improves the results: The dRPAI-I correlation energy of -0.08 Ha decreases to about -0.04 Ha for the dRPA-II, dRPA-IIa, SOSEX, and RPAx-Ia approximations. As expected, the dRPA-IIa correlation energy was found to be exactly half of the dRPA-I value.

iv.3 Application to reaction energies

In this section we discuss a series of results for the reaction energy test set proposed by Hesselmann Hesselmann (2012). In this case the accuracy is evaluated with respect to CCSD(T) benchmark results; CCSD values are also provided for comparison purposes. Similarly to Ref. Hesselmann, 2012, the correlation energy in the complete basis set limit is obtained by extrapolating the values obtained with the aug-cc-pVTZ and aug-cc-pVQZ basis sets. Reaction energy results are detailed in Table 2 for the different approximations considered in this work. The simple dRPA-I DIEL approximation gives substantial mean error (ME) and mean absolute error (MAE): 2.18 kcal/mol and 2.32 kcal/mol, respectively. The dRPA-IIa DIEL approach, with a MAE of 1.99 kcal/mol and a ME of -1.43 kcal/mol, improves to a certain extent the results of dRPA-I DIEL. Table 2 also shows that dRPA-IIa and SOSEX produce similar results, with a difference of at most 0.2 kcal/mol. The use of the full exchange kernel (Eq. 18), as done in dRPA-II AC, deteriorates the accuracy and leads to a MAE that is even larger than in the case of the simple dRPA-I approach. We finally consider the RPAx-Ia DIEL approximation. The corresponding MAE and ME are significantly decreased with respect to the other methods and close to CCSD values. Considering these results for reaction energies and the previous results for atoms and ions in Sec. IV.2, we can notice that the RPAx-Ia approach is the most promising within the dielectric matrix approximations studied in this work. Additionally, RPAx-Ia is more stable than RPAx-I (Eq. II.3) when using PBE as starting point Ángyán et al. (2011). It is important to notice that the PBE and the non-self-consistent PBE0 (post-PBE) approximations lead to the significant MAEs of 4.68 and 4.57 kcal/mol, respectively.

Reaction dRPA-I dRPA-IIa SOSEX dRPA-II RPAx-Ia PBE PBE0 CCSD CCSD(T)
DIEL DIEL rCCD AC DIEL post-PBE
CH+H CH -48.03 -51.97 -51.99 -48.29 -51.57 -52.15 -55.32 -50.46 -49.44
CH+H CH -37.54 -42.41 -42.23 -38.31 -40.26 -40.66 -43.96 -40.49 -39.47
CH+H 2CH -17.50 -19.00 -18.98 -18.15 -18.45 -18.61 -19.01 -18.76 -18.18
CO+H HCHO -3.07 -4.94 -4.98 -1.39 -4.84 -12.45 -12.62 -5.61 -5.47
HCHO+H CHOH -26.88 -32.75 -32.55 -27.77 -30.96 -29.66 -33.44 -30.78 -29.70
HO+H 2HO -82.53 -92.66 -92.48 -86.07 -89.62 -81.81 -86.78 -89.55 -87.63
CH+HO CHCHO -38.01 -39.68 -39.73 -37.33 -39.28 -44.25 -45.49 -38.53 -38.28
CH+HO CHOH -13.15 -15.73 -15.58 -12.77 -14.35 -15.67 -17.85 -14.43 -14.12
CHCHO+H CHOH -23.17 -28.02 -27.85 -23.73 -26.64 -23.57 -27.67 -26.36 -25.28
CO+NH HCONH -6.61 -9.77 -9.83 -5.06 -9.35 -21.04 -19.80 -9.35 -10.26
CO+HO CO+H -5.04 -3.81 -3.85 -1.93 -3.80 -17.62 -12.86 -3.76 -6.18
HNCO+NH NHCONH -17.34 -22.89 -22.81 -18.71 -22.02 -18.45 -23.03 -22.04 -20.70
CO+CHOH HCOOCH -10.07 -12.49 -12.56 -8.54 -12.14 -22.44 -20.85 -12.12 -13.59
HCOOH+NH HCONH+HO -0.79 -1.76 -1.76 -1.17 -1.60 -1.76 -2.03 -1.24 -1.16
CO+HO CO+HO -87.57 -96.47 -96.33 -87.99 -93.42 -99.43 -99.64 -93.31 -93.81
HCCO+HCHO CHO+CO -4.99 -5.67 -5.47 -5.44 -5.26 5.06 0.52 -4.93 -3.83
ME 2.18 -1.43 -1.37 2.15 -0.40 -2.34 -3.92 -0.29
MAE 2.32 1.99 1.90 2.36 1.12 4.68 4.57 0.95
Table 2: Reaction energies for 16 chemical reactions (in kcal/mol). The mean error (ME) and mean absolute error (MAE) are computed with respect to CCSD(T)/CBS. The results with the smallest deviation from the benchmark value are marked in bold face.

The reaction test set in Table 2 involves only energy differences between molecular energies. To better understand the performance of each method it is interesting to also discuss the total energy of each molecule. Detailed results are summarized in Table 3, where the deviation of the total energy of each molecule from CCSD(T) reference values is presented. Since an error cancellation is expected when energy differences are computed, it is not surprising that the approximations in Table 3 lead to MEs and MAEs with respect to CCSD(T) that are substantially larger than what previously seen in Table 2. For total energies the dRPA-I approximation leads to the largest deviation with a MAE of about 168 kcal/mol. The dRPA-IIa DIEL, the SOSEX rCCD, and the dRPA-II AC approximations all have a considerably lower MAE (about 30 kcal/mol). However, it is important to notice that dRPA-II AC tends to overestimate the total energy while dRPA-IIa DIEL and SOSEX rCCD underestimate it. Considering Table 3, the method that gives the most accurate results is RPAx-Ia DIEL, which provides a MAE of less than 9 kcal/mol and outperforms CCSD.

Figure 3: Relative error of the total energy computed with different RPA methods with respect to CCSD(T)/CBS results.
Figure 4: Dissociation curve of H
Molecule
DIEL DIEL rCCD AC DIEL
H -24.67 0.77 0.77 -6.41 -1.92 0.00
CH -92.61 9.34 9.14 -21.57 -1.22 4.60
NH -93.92 12.96 12.70 -18.89 1.56 5.90
HO -94.95 15.86 15.58 -15.22 4.05 6.25
CH -115.23 22.65 22.09 -26.26 6.50 11.82
CH -138.49 20.89 20.30 -31.52 2.45 10.80
CH -161.23 18.73 18.31 -36.76 -0.26 9.78
CO -120.23 26.80 26.14 -23.32 9.69 12.52
HCHO -142.51 28.09 27.39 -25.65 8.39 12.38
CHOH -164.36 25.82 25.31 -30.13 5.21 11.30
HO -170.33 35.98 35.23 -25.59 12.01 14.42
HCCO -186.86 37.15 36.17 -35.69 11.62 19.21
CHO -210.29 36.61 35.79 -39.63 8.91 17.96
CHCHO -209.91 37.11 36.21 -40.53 9.55 17.82
CHOH -232.47 35.14 34.42 -45.38 6.27 16.74
HNCO -188.46 41.28 40.25 -32.07 15.23 21.17
HCONH -210.50 40.25 39.26 -37.01 12.16 19.34
CO -189.37 44.27 43.27 -27.87 18.04 21.19
HCOOH -211.90 43.75 42.74 -33.33 15.09 19.76
NHCONH -279.03 52.05 50.85 -48.97 15.47 25.73
HCOOCH -281.07 53.72 52.47 -48.40 16.35 25.29
ME -167.54 30.44 29.73 -30.96 8.34 14.48
MAE 167.54 30.44 29.73 30.96 8.66 14.48
Table 3: Differences between total energies computed within different approximations and CCSD(T)/CBS reference values (in kcal/mol). Mean error (ME) and mean absolute error (MAE) of total energies are included in the table. The results are presented for the 21 molecules involved in the reactions considered in Table 2 (in kcal/mol). The molecules are sorted in increasing order of the absolute values of CCSD(T) total energies.

To help the visualization of the results in Table 3, we show in Fig. 3 the relative error of total energies obtained within different RPA approximations with respect to CCSD(T) results. The dRPA-I DIEL approximation leads to substantially larger errors and, in order to make the graph more readable, has not been included in the figure. The trend of the curves shows that the relative error tends to be approximately constant by increasing the molecular size. Since CCSD/CBS is exact for two-electron systems it is not surprising that the corresponding relative error for H is zero. For all the other methods H is instead the most problematic case, which is related to the fact that the reference determinant is constructed from PBE orbitals. The dRPA-IIa DIEL and SOSEX rCCD approximations give almost identical curves for all the molecules considered here (the difference is at most 1.2 kcal/mol, for NHCONH and HCOOCH). With the exception of H, the RPAx-Ia relative error is systematically lower than in the CCSD case. However, as seen in Table 2, the accuracy of RPAx-Ia slightly worsens with respect to CCSD in the case of reaction energies, whose evaluation takes advantage of an error cancellation effect.

In Fig. 4, we show dissociation curves for the singlet state of the H molecule, for which CCSD is exact. In the literature we can find several studies of H dissociation in the context of RPA. The influence of the reference orbitals has often been studiedFuchs et al. (2005); Dahlen et al. (2006); Janesko et al. (2009a). The use of unrestricted Kohn-Sham orbitals is thought to allow for some inclusion of static correlation in the description of the dissociationFuchs et al. (2005); Dahlen et al. (2006). Note first that we observed that the Clenshaw-Curtis frequency quadrature is suitable for a certain range of distances only, and that one will need to switch to other quadrature schemes a large distances. Indeed, for large interatomic distances H is a multireference problem characterized by quasi-degenerate determinants and dielectric matrix diverges for , hence the currently implemented integration schemes fail. This is an interesting subject for future work. The dRPA-II AC diverges at larger distances, which can also be explained by the fact that the system is dominated here by non-dynamical correlation. Hence, the coupling-constant integrand is not smooth enough for numerical AC quadratures, as can be inferred from a study of the integrand itself (see for example Fuchs et al. (2005)). The SOSEX curve shows no bump and is consistent with earlier studies in the literatureHenderson and Scuseria (2010); Paier et al. (2010); Ren et al. (2012); Bates and Furche (2013). The dRPA-IIa AC and dRPA-IIa DIEL, although not theoretically equivalent to SOSEX, are both very close to the SOSEX curve. Note that the dRPA-IIa DIEL starts to show a deviation from the dRPA-IIa AC results: this is related to the frequency quadrature problem. The RPAx-Ia-DIEL energy, calculated from a broken-symmetry SCF is very close to the CCSD curve, despite the fact that the fractional deviation between RPAx-Ia-DIEL and CCSD is the worst for H among all the molecules shown on Fig. 3. Note, however, that the absolute deviation is not large in the equilibrium structure and at dissociation it gets even smaller.

V Conclusions, perspectives

The main result of the present work is to show that in contrast to a widely accepted view it is possible to conveniently include exchange effects in the dielectric matrix formulation of RPA correlation energy. Two particular cases have been derived and numerically implemented: the SOSEX-like dRPA-IIa correlation energy and an approximate variant of the RPAx-I method, named RPAx-Ia. Our derivation of the dielectric matrix formulation of the RPA proceeds in a “conventional” or “traditional” way, i.e. from the full adiabatic connection formula at an RPA level. This is in contrast to the work of Eshuis and Furche, who started from the plasmon formula of the dRPA-I energy and derived directly the density fitting expressions. Our approach is strictly equivalent to theirs for dRPA-I but offers in addition a well-defined route to get alternative variants including exchange in a similar form. Furthermore, the connections with the previously studied density matrix formulation are straightforwardly established and therefore our previous results concerning the relationship between the adiabatic connection and rCCD RPA can be simply transferred.

We think that the main interest of the dielectric matrix formalism of RPA in a quantum chemical (more precisely LCAO based) context is mostly conceptual. As demonstrated by the work of Eshuis and Furche, this formalism is well-adapted for density fitting implementations, which is computationally advantageous, in particular for the dRPA-I case, which is an method. The exchange-including variants are expected to show a somewhat less advantageous behavior which is comparable to the scaling of DF-MP2 Werner et al. (2003) or RI-MP2 Weigend and Häser (1997); Weigend (2002). The brute-force orbital-based implementations are of scaling.

It is important to notice the new approximations developed within this work may be of significant interest for the condensed matter physics community, that mostly rely on the dielectric matrix formalism to compute RPA correlation energies Fuchs and Gonze (2002); Fuchs et al. (2005); Harl and Kresse (2008); Lu et al. (2009); Nguyen and de Gironcoli (2009); Rocca (2014).

As shown in this work, the relatively poor performance of the dRPA-I is improved to a certain extent by the dRPA-IIa (SOSEX-like) method. The approximate RPAx-Ia approach is surprisingly good for the systems considered in this paper. In the future more test cases will certainly be necessary to fully establish the accuracy of the RPAx-Ia method. Additionally, this methodology is suitable for implementation in plane-wave codes Rocca (2014); Kaoui and Rocca (2016) and applications to solids will be soon possible.

Vi Aknowledgement

D.R. acknowledges financial support from Agence Nationale de la Recherche under grant number ANR-15-CE29-0003-01. The authors thank Andreas Hesselmann for useful discussions and for providing data on the molecular reaction test set. Computer time was provided by GENCI-CCRT/CINES under grant x2015-085106.

Appendix A Equivalence of the density matrix and logarithmic expressions

In this appendix we show the equivalence of the density matrix and logarithmic expressions of the dRPA-I correlation energies described in this paper. The derivation holds for dRPA-IIa, by replacing with . Let us recall Eq. (12) as:

(40)

Since

(41)

from the “dimension-reduced” Dyson equation , we have

(42)

where

(43)

The symmetric matrix can be brought to a diagonal form using its eigenvectors and eigenvalues:

(44)

Thus we find:

(45)

and furthermore:

(46)

Noting that

(47)

where labels the eigenvalues of , i.e. the diagonal elements of , we get

(48)

Furthermore, since , we have

(49)

which finally leads to