Prof Werner

Andreas Heÿelmann

Andreas Görling

Random phase approximation correlation energies with exact Kohn-Sham exchange

Keywords: random phase approximation, exact exchange, Kohn-Sham, orbitaldependent functional, exact exchange kernel

The random phase approximation (RPA) correlation energy is expressed in terms of the exact local Kohn-Sham (KS) exchange potential and corresponding adiabatic and nonadiabatic exchange kernels for density-functional reference determinants. The approach naturally extends the RPA method in which, conventionally, only Coulomb interactions are included. By comparison with the coupled cluster singles doubles with perturbative triples method it is shown for a set of small molecules that the new RPA method based on KS exchange yields correlation energies more accurate than RPA on the basis of Hartree-Fock exchange.

Introduction

The description of electron correlation eects has been shown to be essential for an accurate description of molecular energetics, structures and properties. Because of this, in recent years, there has been much eort to improve the methods that are in use to determine correlated wave functions and energies and also the underlying algorithms that are used for their computer implementation. Concering the latter point the most signicant developments probably were density-tting methods [START_REF] Bartlett | ManyBody Methods in Chemistry and Physics: MBPT and CoupledCluster Theory[END_REF] that reduce the computational cost with respect to the basis set size and secondly methods that treat electron correlation using localised orbitals [59].

The nowadays most widely used approaches to describe electron correlation effects are Møller-Plesset (MP) perturbation theory [1014], coupled-cluster (CC) methods [START_REF] Helgaker | Molecular ElectronicStructure Theory[END_REF][START_REF] Piela | Ideas of quantum chemistry[END_REF][START_REF] Bartlett | ManyBody Methods in Chemistry and Physics: MBPT and CoupledCluster Theory[END_REF] and density-functional theory (DFT) [START_REF] Piela | Ideas of quantum chemistry[END_REF]1618]. Principally all three approaches aim at describing electron correlation exactly, i.e. they offer the possibility to obtain the exact solutions of the nonrelativistic Schrödinger equation. In praxis, however, approximations have to be introduced: in MP and CC theory the wave function is usually restricted to include at most triple excitations from the given reference determinant while in the case of DFT the yet unknown exchange-correlation (xc) functional has to be approximated.

The most prominent approximations of the xc functional are the local density approximation (LDA) [1618] and generalised gradient approximation (GGA) functionals [1620] that describe the xc energy in terms of the electron density and its gradient. More recently these methods were further improved by introducing functionals depending on the laplacian of the density or the kinetic energy density (meta GGA's) [2124] or that include fractions of nonlocal exchange (hyper GGA's) [2528]. Unfortunately it was found that GGA's and many of their extensions suffer from a number of shortcomings, most importantly they are unable to describe long-range correlation eects or dispersion interactions [START_REF] Koch | A Chemist's Guide to Density Functional Theory[END_REF]29].

An important step forward in the development of density functionals was there- fore to introduce functionals that are dependent both on occupied and unoccupied molecular orbitals and thus only implicitly dependent on the electron density [30 47]. While such orbital-dependent functionals on the one hand have the advantage that exchange interactions can be treated exactly [3537,41,48,49] and thus the self-interaction error of the Coulomb interaction term is cancelled exactly, on the other hand, within the framework of orbital-dependent functionals, accurate expressions for the correlation functional can be obtained from many-body perturbation theory [START_REF] Görling | [END_REF]33,34,45,5055] or coupled-cluster theory [42]. Indeed it has been demonstrated that orbital-dependent functionals derived from second-order perturbation theory are capable to describe dispersion interactions [56] and yield improved total energies and molecular properties compared with standard GGA functionals [39,40]. There are, however, a number of points that limit the use of correlation functionals derived from perturbation theory: rst the standard secondorder correlation energy expression is variationally unstable [52,57] and so has to be modied in the energy denominators [40,45], second the use of perturbation theory methods is problematic for extended and periodic systems which have a small or vanishing band gap. Finally it was shown in case of Møller-Plesset perturbation theory that the perturbation expansion of the correlation energy or of molecular properties often diverges [5861].

F
Because of this there is the need to derive new orbital-dependent functionals that sum certain types of perturbation diagrams up to innity as is for example done in coupled-cluster theory. One such method is the random phase approximation (RPA) [START_REF] Harris | Algebraic and diagrammatic methods in many fermion theory[END_REF][START_REF] Mcweeny | Methods of molecular quantum mechanics[END_REF]6270] which can be seen as an approximate coupled-cluster doubles (CCD) approach in which the doubles amplitudes are restricted to those portions that lead to ring diagrams [START_REF] Harris | Algebraic and diagrammatic methods in many fermion theory[END_REF][START_REF] Scuseria | [END_REF].

The name random phase approximation stems from a classical mechanical treatment of the collective properties of the electron gas by Bohm and Pines [72,73]. Later on Nozieres and Pines demonstrated that in a quantum-mechanical framework RPA is equivalent to the addition of singly excited determinants to an independent particle determinant [74] and Ehrenreich et al. have then shown the equivalence with the perturbed self-consistent eld method [75]. Historically, however, the earliest RPA calculations were not based on antisymmetrised wave functions and RPA was originally associated with ring sums for the electron gas in which the interactions were not antisymmetrised. This term of the RPA is still used and usually refered to in solid state physics (see e.g. Ref. [76]) while in quantum chemistry the term RPA is often used as another term for time-dependent Hartree-Fock (TDHF) because RPA here is used mainly for the determination of excitation energies and response properties [START_REF] Mcweeny | Methods of molecular quantum mechanics[END_REF]7779]. In order to avoid the confusion that arises due to these dierent denitions it has become common to term RPA methods including exchange interactions either 'RPA with exchange' or 'full RPA' while RPA methods excluding exchange terms are termed just RPA or 'direct RPA'.

The RPA method without exchange has been shown to yield the exact correlation energy of the electron gas in the high density limit [START_REF] Harris | Algebraic and diagrammatic methods in many fermion theory[END_REF][START_REF] Mattuck | A guide to Feynman diagrams in the manybody problem[END_REF] and thus may be very useful for larger and particularly extended electronic systems. However, from a diagrammatical point of view 'direct RPA' has the disadvantage that it contains many exclusion principle violating (EPV) diagrams whose cancelling counterparts are absent. This deciency is most severe for small systems with small basis sets and this might explain that it did not play a signicant role in quantum chemical treatments. This is dierent for the 'RPA with exchange' which apart from its use for the calculation of response properties and excitation energies has also been used to directly determine non-multipole expanded dispersion energies [START_REF] Mcweeny | Methods of molecular quantum mechanics[END_REF]8082] correlation energies [8386]. In 1977 Szabo and Ostlund analysed the full RPA correlation contribution to the interaction energy of two remote molecules and showed that the dispersive part of the interaction is described on the coupled Hartree-Fock level (identical with TDHF) [84]. The 'full RPA' method should therefore be more suited to describe intermolecular interactions than second-order Møller-Plesset theory which describes the dispersion energy contribution only on an uncoupled Hartree-Fock level [84,87]. However, the full RPA method nevertheless has not been widely used for the determination of correlation energies and this may stem from the fact that coupled-cluster methods turned out to be more accurate while not much more expensive.

Very recently a revival of RPA methods in quantum chemistry has occured as RPA was found to be an attractive starting point for the development of new orbitaldependent density functionals [44,46,[START_REF] Scuseria | [END_REF]8894]. Compared with the second-order perturbation theory functionals described above, RPA based functionals do not possess a variational instability, can be applied also to extended systems with small band gaps and they describe dispersion interactions on a coupled TDDFT level which has been found to give very accurate dispersion energies in the framework of the DFT-SAPT (symmetry-adapted intermolecular perturbation theory) method [95]. Moreover, Furche has shown [90] that the correlation energy of the direct RPA method can be computed with an ecient algorithm that scales only as N 5 with the molecular size and therefore the computational cost does not exceed the cost from perturbation theory based functionals. However, while for bulk properties of solids the direct RPA method has been shown to give a good agreement with experimental results [76], for molecular systems its performance for, e.g., the prediction of atomization energies was found to be not much better than with standard GGA functionals [89]. Because of this the exchange(-correlation) interactions omitted in direct RPA have to be accounted for and a straightforward approach to do so is the inclusion of local exchange-correlation kernels derived from the functional derivatives of GGA functionals [44]. It has been shown that in conjunction with hybrid functionals this extension to the RPA, termed uctuation-dissipation theorem DFT (FDT-DFT), yields atomization energies and intermolecular interactions that outperform direct RPA results [44]. However, it has been shown that the use of local xc kernels in the framework of FDT-DFT leads to unphysical pair density functions that diverge for small interelectronic distances and this apparently leads to a slow basis set convergence of correlation energies [44].

One solution to this problem consists in a separation of the electron interaction into a short-range and a long-range part and to treat the RPA correlation for the long-range part only while the short-range interactions are described using GGA functionals. Such short-range-long-range functionals have very recently been developed by Scuseria et al. [92,93] and Toulouse et al. [94]. Both types of functionals were especially tested for the description of van-der-Waals interactions and it was found that they perform exceedingly better than direct RPA [93,94]. It should be noted here that the long-range RPA part is dierent in both functionals: the functional from Scuseria et al. uses direct RPA only while the functional from Toulouse et al. in addition contains the exchange kernel according to the TDHF formalism. Another extension to the RPA correlation functional was proposed by Jiang and Engel that was based on adding the second-order exchange energy perturbation expression to the direct RPA energy [46]. This approach, however, leads to the same problematic properties of the functional as described above for the pure second-order perturbation theory functionals.

A natural extension to direct RPA functionals are full RPA functionals that account for exchange interactions by using directly the exact exchange kernel of DFT [9698]. The exact exchange kernel, that is the functional derivative of the exact KS-exchange potential with respect to the electron density in contrast to standard GGA kernels is both nonlocal and frequency-dependent and therefore does not lead to an unphysical behaviour of the pair density for small interelectronic distances. The exact KS exchange potential is determined and used in EXX (exactexchange) methods [35,36,41] within the framework of the optimised eective potential method [START_REF] Görling | [END_REF]99,100] that yields derivatives of orbital-dependent functionals with respect to the electron density. First implementations of the adiabatic and nonadiabatic KS exchange kernel [101,102], however, relying on expansions of the Kohn-Sham response function in auxiliary basis sets suered from numerical instabilities and are therefore impractical for the accurate determination of response properties and energies. We have recently implemented the exact-exchange kernel using a modied TDDFT method (Ref. [103], see also section 2.3) and this new approach has been used to calculate excitation energies of weakly bound dimers, demonstrating that the time-dependent EXX method is capable to describe chargetransfer excitations. Independently a time-dependent EXX method for atomic systems using cubic splines radial basis functions has been implemented by Hellgren and von Barth and it was shown that this approach yields excitation energies and dispersion coecients that are close to the corresponding TDHF results [91,104]. Furthermore Hellgren and von Barth have also calculated correlation energies using RPA with the exact-exchange kernel for a few closed-shell atomic systems and report a signicantly good agreement with accurate conguration interaction (CI) correlation energies [91].

In this work a RPA method with exact KS exchange (RPA(EXX)) is presented that is applicable to general molecular systems. For this the TDDFT formalism for orbital-dependent xc kernels from Ref. [103] will be used that avoids the numerically problematic inversion of the Kohn-Sham response matrices. It will be shown in section 2.3 that from the full RPA(EXX) method two further approximations can be devised, namely an RPA method including only the adiabatic exact-exchange kernel (RPA(AEXX)) and an RPA method with an exchange kernel obtained by scaling the nonadiabatic part of the full exact-exchange kernel with one half (RPA(EXXh)). Section 4 shows and discusses the results for total energies, correlation energies and reaction energies for a set of small organic molecules. Finally section 5 summarises and concludes.

Theory

Electron correlation energies from the uctuation-dissipation theorem

The electronic energy of a many-body interacting system is usually split into the reference energy of a single-determinant wave function and a remainder termed the correlation energy. The reference energy is given by:

E 0 = Φ 0 | Ĥ|Φ 0 = occ i 2 φ i | - 1 2 ∇ 2 + v ext |φ i + occ ij 2(ii|jj) -(ij|ij)
where Φ 0 is the reference determinant, Ĥ is the electronic hamilton operator, φ i denotes an occupied orbital, v ext is the external potential (usually containing electronnucleus interactions) and (ij|kl) is a two-electron integral in chemists notation. Note that, in the following, for convenience closed-shell formalism will be considered and spin will be taken into account by appropriate prefactors.

E c = Ψ| Ĥ|Ψ -Φ 0 | Ĥ|Φ 0 (1)
depends on the reference determinant Φ 0 which a priory is not specied (Ψ may be the exact or an approximation to the exact wave function). If Φ 0 is chosen as the HF determinant the denition of the correlation energy results that is usually refered to in wave-function based quantum chemistry methods. In density-functional theory Φ 0 is the KS determinant. This leads to another denition of the correlation energy.

In this work the reference determinant will be the EXX determinant with orbitals obtained from an exact-exchange method [41] which does not yield Hartree-Fock but EXX orbitals.

Using the uctuation-dissipation theorem [105,106] the correlation energy of the interacting many-body system can be obtained from a coupling strength integration over the Coulomb-type integral of the correlated part of the pair correlation function P c α [44,105,106]:

E c = 1 2 1 0 dα dr 1 dr 2 P c α (r 1 , r 2 ) r 12 = - 1 2π 1 0 dα dr 1 dr 2 1 r 12 ∞ 0 dω χ α (r 1 , r 2 , ω) -χ 0 (r 1 , r 2 , ω) ( 2
)
where α is the coupling strength, χ 0 is the uncoupled response function and χ α is the coupled response function at coupling strength α. The integration over ω goes over imaginary frequencies, i.e. here and in the following ω is dened to be an imaginary quantity. The uncoupled response function is given in terms of the orbitals and eigenvalues of the noninteracting KS system:

χ 0 (r 1 , r 2 , ω) = ia 4ε ia ε 2 ia + ω 2 φ ia (r 1 )φ ia (r 2 )
where φ ia (r) = φ i (r)φ a (r) denotes an occupied-virtual orbital product and the convention ε ia = ε a -ε i is used. The coupled response function χ α can be obtained from the Dyson-type equation:

χ α (r 1 , r 2 , ω) = χ 0 (r 1 , r 2 , ω)+ dr 3 dr 4 χ 0 (r 1 , r 3 , ω) α r 34 + f α xc (r 3 , r 4 , ω) χ α (r 4 , r 2 , ω) (3)
which is rewritten symbolically by

χ α = χ 0 + χ 0 W χ α ( 4
)
where W = 1/r 12 + f xc is the interaction operator comprising Coulomb-, exchangeand correlation-eects. It can be seen that Eq.( 4) can be solved iteratively obtaining an n-order expansion of the coupled response function: and in innite order the interacting response function is given by:

χ (n) α = χ 0 + χ 0 W χ 0 + χ 0 W χ 0 W χ 0 + . . . ( 5
χ α = (1 -χ 0 W ) -1 χ 0 = (χ -1 0 -W ) -1 (6)
The response function can also be written as (we assume in the following an implict dependency on the coupling strength α)

χ α (r 1 , r 2 , ω) = p ia,jb 4ω p ω 2 p + ω 2 U p,ia U p,jb φ ia (r 1 )φ jb (r 2 )
where ω p are the excitation energies and U p are the corresponding eigenvectors of the TDDFT eigenvalue equation:

ε 2 -4ε 1/2 K(ω)ε 1/2 U = ω 2 U (7)
with

K ia,jb (ω) = dr 1 dr 2 φ ia (r 1 )W (r 1 , r 2 , ω)φ jb (r 2 )
and ε is a diagonal matrix with the elements ε ia,jb = δ ia,jb (ε a -ε i ).

It is now assumed that the coupling matrix K and thus the eigenvectors U are frequency-independent. Then the ω-integration in Eq.( 2) can be done analytically. If in addition also the Coulomb-type integrals over r 1 and r 2 are performed the correlation energy can be rewritten as:

E c = - 1 0 dα p ia,jb U p,ia U p,jb (ia|jb) - ia (ia|ia) = - 1 0 dα p ia,jb U p,ia U p,jb (ia|jb) + 1 2 ia (ia|ia) (8)
If however the xc-kernel is frequency dependent the correlation energy would have to be obtained from:

E c = - 1 2π 1 0 dα ∞ 0 dω   p ia,jb 4ω p ω 2 p + ω 2 U p,ia (ω)U p,jb (ω)(ia|jb)   + 1 2 ia (ia|ia)
Alternatively the computation of the correlation energy can also be done without solving the eigenvalue equation (7). For this we dene: (compare with Eq.( 6)) so that E c can be written as:

λ ia,jb (ω) = δ ia,jb 4ε ai ε 2 ai + ω 2 K ia,jb (ω) = ia 1 r 12 + f xc (ω) jb χ α = (1 -αλK) -1 λ (9) F o r P e e r
E c = - 1 2π 1 0 dα ∞ 0 dω ia,jb (χ α ) ia,jb -λ ia,jb (ia|jb) (10)
Note that of course again the uncoupled part involving an αand ω-integration of the λ-term can be separated as in Eq.( 8) and again yields the term 1 2 ia (ia|ia).

Coupling strength integration

If the exchange-correlation kernel f α xc is approximated by the exchange-kernel

f α x = αf 1 x = αf x with f x = f α=1 x then the electron interaction operator W α = α r12 + f α xc turns into W α = αW 1 = αW = α( 1 r12 + f x ) with W = W α=1
and depends linearly on the coupling strength. The response function at coupling strength α then is given by:

χ α = χ 0 + αχ 0 W χ 0 + α 2 χ 0 W χ 0 W χ 0 + . . .
and in the occupied-virtual orbital product space (Eq. ( 9)):

χ α = λ + αλKλ + α 2 λKλKλ + . . . ( 11
)
so that the coupling strength integration gives:

1 0 dα χ α = λ + 1 2 λKλ + 1 3 λKλKλ + 1 4 λKλKλKλ + . . .
It is now convenient to exploit the fact that λ is diagonal and its square root can easily be taken. Using K = λ 1/2 Kλ 1/2 the coupling strength integrated response matrix can be rewritten:

1 0 dα χ α = λ 1/2 1 + 1 2 K + 1 3 K K + 1 4 K K K + . . . λ 1/2 .
Compared with the power series:

-ln(1 -x) = x + 1 2 x 2 + 1 3 x 3 + 1 4 x 4 + . . .
the coupling strength integrated response function can be expressed by

1 0 dα χ α = λ 1/2 -ln(1 -K) K-1 λ 1/2 (12)
The logarithm of the symmetric matrix(1 -K) can be evaluated analytically as:

ln(1 -K) = U T ln(Λ) U (13)
where the matrix U contains the eigenvectors and Λ is a diagonal matrix containing the eigenvalues of the symmetric matrix (1 -K). Note that the precondition for applying Eq. ( 13) is of course that the matrix (1 -K) is positive denite. This was the case for all systems that were studied in this work and so the logarithm according to Eq. ( 13) could always be taken. Finally note also that apart from an eigenvalue decomposition of (1 -K) other numerical techniques exist to compute logarithms of matrices that might be more ecient for larger matrices [107109]. Using Eq. ( 13) the correlation energy can now be computed as:

E c = - 1 2π ∞ 0 dω ia,jb λ 1/2 -U T ln(Λ) U K-1 λ 1/2 -λ ia,jb
(ia|jb) .

FD-DFT using the exact-exchange kernel

The frequency dependent exact-exchange kernel can generally be written as [9698]:

f x (r 1 , r 2 , ω) = dr 3 dr 4 χ -1 0 (r 1 , r 3 , ω)h x (r 3 , r 4 , ω)χ -1 0 (r 4 , r 2 , ω) (14)
with

h x = h [1]
x + h [2] x given by (note that ω is assumed to be imaginary, so for real frequencies the signs in front of the ω 2 contributions would have to be changed):

h [1] x (r

1 , r 2 , ω) = ij,ab -4ε ia ε jb + 4ω 2 (ε 2 ia + ω 2 )(ε 2 jb + ω 2 ) (ij|ab) + δ ab i|v NL x -vx |j -δ ij a|v NL x -vx |b + + -4ε ia ε jb -4ω 2 (ε 2 ia + ω 2 )(ε 2 jb + ω 2 ) (ib|ja) φ ia (r 1 )φ jb (r 2 ) (15)
h [2] x (r

1 , r 2 , ω) = ij,a -4ε ja ε ia (ε 2 ja + ω 2 ) i|v NL x -vx |a φ ij (r 1 )φ ja (r 2 ) + φ ja (r 1 )φ ij (r 2 ) + + i,ab -4ε ib ε ia (ε 2 ib + ω 2 ) i|v NL x -vx |a φ ab (r 1 )φ ib (r 2 ) + φ ib (r 1 )φ ab (r 2 ) (16)
with i|v NL x -vx |j giving the dierences of the matrix elements of the nonlocal (v NL

x ) and local (v x ) exchange potentials. Unfortunately the computation of the exact-exchange kernel via Eq. ( 14) using nite basis sets is problematic due to the occurence of the inverses of the uncoupled response functions χ 0 . Because of this in Ref. [103] the TDDFT response equations were reformulated yielding the response to the eective Kohn-Sham potential and not, as usually, the electron density. With this a close approximation of the exact-exchange kernel is obtained that does not involve the inverses of the KS response functions anymore and that would be the exact functional derivative of the exact local KS exchange potential if the occupiedvirtual orbital products were linearly independent. As in case of the exact-exchange kernel dened in Eqs. [START_REF] Piela | Ideas of quantum chemistry[END_REF][START_REF] Bartlett | ManyBody Methods in Chemistry and Physics: MBPT and CoupledCluster Theory[END_REF][START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF] it can be split into two parts and the rst part can be written as [103]:

[K [1] x ] ia,jb = 1 -

ω 2 ε ia ε jb [-(ij|ab) + δ ij a|v NL x -vx |b -δ ab i|v NL x -vx |j ] -1 + ω 2 ε ia ε jb (ib|ja) (17)
while the second part, corresponding to Eq. ( 16), contains the occupied-unoccupied matrix elements i|v NL x -vx |a . While its computation is more demanding than the x (see Ref. [103] for details) it has been found that its inclusion has a negligible eect for excitation energies or correlation energies of small molecules. This is concordant with the fact that the EXX orbitals can approximately be obtained from an occupied-occupied virtual-virtual unitary transformation of the Hartree-Fock orbitals which is explicitly exploited in approximate exact-exchange density functional theories [37,49]. Therefore the K [2] x contribution has not been included in the calculations of this work.

Since the exchange-kernel of Eq. ( 17) is frequency-dependent we can dene two additional approximations that will be used in the calculations: adiabatic approximation (AEXX): ω-dependent terms in Eq. ( 17) are omitted EXXh: ω-dependent terms in Eq. ( 17) are scaled with one half It will now be shown that the correlation energy dened in terms of the exchangekernel of Eq. ( 17) is related to the leading, i.e., second order, term in perturbation theory along the adiabatic connection which is the analogue of Møller-Plesset perturbation theory within the KS formalism [START_REF] Görling | [END_REF]100]. In order to show this the power series of Eq. ( 11) is truncated after the second term and inserted in Eq. ( 10) yielding:

E (2) c = - 1 2π 1 0 dα ∞ 0 dω ia,jb α(λKλ) ia,jb (ia|jb) = - 1 4π ∞ 0 dω ia,jb (λKλ) ia,jb (ia|jb)
The coupling matrix K can be split into a Coulomb contribution and an exchange contribution:

K = K C + K x
with K x given by Eq. ( 17) and [K C ] ia,jb = (ia|jb). Accordingly the second-order correlation energy is split into a Coulomb-and an exchange-part. For the Coulombpart the frequency integration yields (see also Ref. [89]):

E (2) c,C = -2(ia|jb)(ia|jb) ε ia + ε jb (18)
which is the direct part of the second-order correlation energy. For the exchangepart one can make use of the integrals

∞ 0 dω xy -ω 2 (x 2 + ω 2 )(y 2 + ω 2 ) = 0 ∀ x 2 , y 2 , xy > 0 ∞ 0 dω xy + ω 2 (x 2 + ω 2 )(y 2 + ω 2 ) = π x + y ∀ x 2 , y 2 , xy > 0
and obtains:

E (2) c,x = (ia|jb)(ib|ja) ε ia + ε jb ( 19
)
which shows that the rst part of the exchange-kernel of Eq. ( 17) vanishes upon ω-integration and does not contribute to the second-order correlation energy. The c,x dened in Eqns. ( 18) and ( 19) yields the correlation energy through second-order in the interelectronic interaction as given in Refs. [START_REF] Görling | [END_REF]100] except for terms proportional to | i|v NL

x -vx |a | 2 that are small in magnitude and are neglected here. The neglected terms have their origin in the fact that products of occupied and unoccupied orbitals are linear dependent, which is also the reason why HF and exact exchange-only KS determinants are dierent [110,111].

Finally it is noted that while the full RPA approach with exact KS exchange described in this section scales as N 6 with the molecular size we have recently implemented the method exploiting density-tting of occupied-virtual orbital products using a similar approach as the one used in Ref. [112] to obtain the response functions for the TDHF or hybrid-DFT methods. The density-tting RPA(EXX) method only scales as N 5 with the molecular size and is therefore applicable to relatively large molecular systems. This will be demonstrated for weakly bound van-der-Waals complexes in a subsequent publication [START_REF] Heÿelmann | ????[END_REF].

Computational Details

Total energies and correlation energies were calculated with the approach described in section 2 for a set of 21 small organic molecules. The geometries used in this work were optimised at the MP2 (second-order Møller-Plesset perturbation theory) level with the aug-cc-pVTZ basis set [START_REF] Kendall | [END_REF] and are the same as those used in Ref. [115].

The exact-exchange Kohn-Sham (EXX) orbitals and eigenvalues were obtained by rst performing numerically stable EXX calculations with the method from Ref. [41] for all 21 molecules in twos steps. In the rst step the exact local KS exchange potential was obtained with the balanced uncontracted triple-zeta orbital-and auxiliary basis sets of Ref. [41]. In the second step EXX orbitals and eigenvalues were calculated by using the exact local exchange-potential represented in the auxiliary basis set obtained in the rst step in a subsequent EXX calculation with the smaller contracted aug-cc-pVTZ orbital basis set [START_REF] Kendall | [END_REF]. In the EXX calculation of the second step only the Coulomb-potential was optimised self-consistently. Correlation energies using the approach described in section 2 were then calculated using the EXX orbitals and eigenvalues as input, i.e., the calculations were performed in a post Kohn-Sham (post selfconsistent-eld) way.

While the coupling-strength integration required for the computation of the correlation energies (see e.g. Eq. ( 2)) could be done analytically using the approach described in section 2.2, the ω-integration over imaginary frequencies can not be done analytically for frequency-dependent exchange-correlation kernels. In this work we have used a 12-point Gauss-Legendre quadrature [116] for the frequency integration which was found to be accurate to approximately 10 -5 hartree.

For comparison the molecular energies were also calculated using Hartree-Fock (HF), MP2, coupled cluster with singles and doubles (CCSD) and additionally also perturbative triples (CCSD(T)), the Becke-Lee-Yang-Parr 3-parameter hybrid density functional (B3LYP) [25,26], and the random-phase approximation with exact Hartree-Fock exchange (RPA(HF)). The correlation energy for the latter approach can be obtained from Eq. ( 8) if the response eigenvectors U are given by the solutions of the time-dependent Hartree-Fock equations [START_REF] Mcweeny | Methods of molecular quantum mechanics[END_REF]79].

All calculations were done using the aug-cc-pVTZbasis set [START_REF] Kendall | [END_REF] and all electrons were correlated. The method described in section 2 has been implemented in the development version of the Molpro quantum chemistry package [START_REF] Werner | MOLPRO, version 2008.2, a package of ab initio programs[END_REF]. The total molecular energies calculated for the set of 21 molecules are displayed in table 1. The last three lines of the table show the root-mean square errors, the mean absolute errors and the relative percentual deviations from the CCSD(T) reference values of the fth column. It can be observed that while the MP2 method generally captures over 90 percent of the correlation eect, with only few exceptions the CCSD energies are closer to the CCSD(T) reference results than the MP2 ones, as expected. In contrast to this the B3LYP method almost always yields total energies that are about 0.1 -0.2 hartree larger in magnitude than the respective CCSD(T) values. On the other hand the RPA method with Hartree-Fock exchange (RPA(HF)) tends to underestimate the CCSD(T) energies slightly. It is interesting to note here that while RPA(HF) is superior to MP2 from the theoretical point of view (it sums up certain types of diagrams to innity including all second-order diagrams), on average it does not give more accurate total energies than the MP2 method as can be obeserved from table 1. For example the mean absolute errors from RPA(HF) are about twice as large as the MAE's of the MP2 method for the given set of molecules. In contrast to this the RPA method with KS exchange (RPA(EXX)) gives total energies that are almost always closer to CCSD(T) than the MP2 energies and actually the RPA(EXX) energies are very close to the CCSD energies, exhibiting similar RMS errors and MAE. It can also be seen that the RPA(EXX) energies are always larger than their RPA(HF) counterparts. This is likely due to the dierent eigenvalue spectrum of the HF and EXX method, more precisely the EXX method yields bound unoccupied orbitals and thus smaller energy gaps between the occupied and lowest unoccupied molecular orbitals. In contrast to the RPA(EXX) energies the RPA energies from the adiabatic exchange kernel (RPA(AEXX)) always overestimate the CCSD(T) energies and this lead us to devise the RPA(EXXh) approach with an EXX kernel in which the frequency dependent part is scaled by one half. The total energies of this method are shown in the last column of table 1. It can be observed that the correspondence with the CCSD(T) reference values is excellent, the mean absolute deviations being only about 5 millihartree for the studied systems. The diagram in gure 1 shows the deviations of the correlation energies of the respective methods from the CCSD(T) reference correlation energies. Note that the correlation energy displayed in Fig. 1 in all cases is dened as the dierence of the total energies and the Hartree-Fock energies of the second column of table 1 in order to consider comparable values. Fig. 1 also shows that RPA(EXX) correlation energies are close to CCSD ones while being always larger than their RPA(HF) counterparts, again likely due to the dierent eigenvalue spectrum of HF and EXX. The diagram in Fig. 1 also shows clearly that the correlation energies of the RPA(EXXh) method interpolate between the RPA(AEXX) and RPA(EXX) energies and as a consequence are very close to the CCSD(T) reference correlation energies. This is also summarised in the diagram in Fig. 2 which displays the absolute deviations from the CCSD(T) correlation energies for all methods averaged over the 21 molecules studied in this work. Fig. 2 shows that the RPA(AEXX) and RPA(EXX) method clearly improve the accuracy of the correlation energy if compared with the corresponding RPA(HF) method, but they have about the same average deviation of 25 millihartree as obtained with the MP2 and CCSD methods. On the other hand the RPA(EXXh) correlation energies deviate from the CCSD(T) reference energies only by 5 millihartree on average.

In order to assess the presented RPA approach with KS exchange for the prediction of reaction energies we have considered the 14 chemical reactions listed 2. Note that compared to the prediction of total energies this is an even more dicult test for electron correlation methods as the energy dierences are typically an order of magnitude smaller than total (correlation) energies and for the reactions studied here they lie in the range of 8 to 80 millihartree. The input energies used for the calculation of the reaction energies all were taken from table 1.

Figure 3 shows the statistical error estimates for all methods averaged over the 14 reactions with CCSD(T) as reference. Interestingly the B3LYP method that yields large deviations from CCSD(T) total energies gives a comparable accuracy for the chemical reactions of table 2 as MP2. The RPA(EXXh) method on the other hand yields mean absolute errors for the reaction energies of about 1.7 kcal/mol which is about 0.7 kcal/mol larger than the deviations obtained with the CCSD method, despite the rather accurate total energies obtained with the RPA(EXXh) method. It can however be observed from Fig. 3 that both, the RPA(EXX) and RPA(EXXh) methods improve upon the accuracy of RPA(HF) for the chemical reactions while RPA(AEXX) and RPA(HF) itself give reaction enthalpies comparable to MP2. While the comparison with CCSD(T) reaction energies is certainly the ultimate test for any lower level electron correlation method, it may be more equitable to use the CCSD method as the reference approach since both, MP2 and the RPA methods account only for double excitations and thus are rather approximations to the CCSD method. Therefore deviations to CCSD for the chemical reactions are displayed in Fig. 4. Here it can now be seen that the RPA(EXX) and RPA(EXXh) methods perform clearly better than MP2 and even than the related RPA(HF) method. For example the root mean square errors are only 0.9 and 1.7 kcal/mol while the MP2 method has a root mean square error of nearly 12 kcal/mol for the chemical reactions. As expected the RPA(AEXX) does not perform as well as the approaches with frequency dependent kernels, but anyhow it performs as well as RPA(HF) and also improves the MP2 method.

Summary

Random phase approximation methods for the correlation energy including exchange interations have been developed for exact-exchange Kohn-Sham reference determinants using the recently developed TDDFT formalism for orbital-dependent xc functionals that avoids the numerically dicult inversion of the Kohn-Sham response matrix. The exchange kernel derived from this approach, in contrast to standard GGA kernels, is both nonlocal and frequency dependent and therefore improves upon the unphysical behaviour of the pair density at small interelectronic distances occuring in the case of local GGA kernels. Besides RPA with the full exact-exchange kernel (RPA(EXX)) two approximate methods were derived, one in which the frequency dependent part is completely neglected (adiabatic approximation, RPA(AEXX)) and one in which the frequency dependent part is scaled with one half (RPA(EXXh)).

It has been shown for a set of 21 small organic molecules that the RPA(EXX) method underestimates the correlation energy compared to coupled cluster singles doubles with perturbative triples (CCSD(T)) correlation energies while the RPA(AEXX)) overestimates the CCSD(T) correlation energies to about the same amount. Accordingly the RPA(EXXh) method yields correlation energies that are close to the CCSD(T) reference energies with an average deviation of only 0.005 hartree for the systems that were studied. All three RPA approaches with exact KS exchange give correlation energies that are on average in a much better agreement with CCSD(T) correlation energies than the corresponding RPA method including Hartree-Fock exact-exchange (RPA(HF)) and the RPA(EXXh) method clearly outperforms second-oder Møller-Plesset (MP2) and CCSD in this respect. However, it was found that the reaction energies for 14 organic reactions obtained by all three introduced RPA methods are only slightly better than corresponding MP2 and RPA(HF) reaction energies if compared to CCSD(T) reference energies and gave worse results than CCSD. However, if compared to CCSD as reference method for the reaction energies the RPA(EXX) and RPA(EXXh) method performed clearly better than the MP2 and RPA(HF) method for the reactions that were studied. This nding is signicant from the point of view that the computational cost of full RPA(EXX) calculations should be not much higher than that of MP2 if density-tting methods are exploited, as is shown in Ref. [START_REF] Heÿelmann | ????[END_REF], since both methods then scale as N 5 with respect to the molecular size. The presented RPA(EXX) method is therefore an orbital-dependent functional that can be applied to large molecules and may become an attractive method for quantum chemistry applications in the future.

F

3.4 • 10 -1 1.2 • 10 -3 5.3 • 10 -4 2.0 • 10 -2 5.7 • 10 -3 1.2 • 10 -3 7.5 • 10 -4 4.1 • 10 -5 MAE 5.4 • 10 -1 3.2 • 10 -2 2.1 • 10 -2 1.3 • 10 -1 7.0 • 10 -2 3.2 • 10 -2 2.4 • 10 -2 5.3 • 10 -3 ME 5.4 • 10 -1 3.2 • 10 -2 2.1 • 10 -2 -1.3 • 10 -1 7.0 • 10 -2 -3.2 • 10 -2 2.

H

Figure 1 .

 1 Figure1. Deviations of correlation energies (dened in all cases as dierence of the total energy with the HF total energy) from CCSD(T) correlation energies, aug-cc-pVTZ basis set.

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Mean absolute deviations from CCSD(T) correlation energies (dened in all cases as dierence of the total energy with the HF total energy) of 21 molecules, aug-cc-pVTZ basis set.

Table 1 .

 1 Total energies in hartree, aug-cc-pVTZ basis set. The last four lines display the root mean squared error (RMS), the mean absolute error (MAE), the mean error (ME) and the total relative deviation (|∆|) to the CCSD(T) reference values for each method.

	eptember 30, 2009	18:26	Molecular Physics	rpa_exx
			o r
			P
				e
				e r
				R
				e
				v i e
				w
				O n l
				y

URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 2 .

 2 Reactions C 2 H 2 +H 2 → C 2 H 4 C 2 H 4 +H 2 → C 2 H 6 C 2 H 6 +H 2 → 2CH 4 CO+H 2 → H 2 CO H 2 CO+H 2 → CH 3 OH H 2 O 2 +H 2 → 2H 2 O C 2 H 2 +H 2 O → CH 3 CHO C 2 H 4 +H 2 O → C 2 H 5 OH CH 3 CHO+H 2 → C 2 H 5 OH CO+NH 3 → HCONH 2 CO+H 2 O → CO 2 +H 2 HNCO+NH 3 → NH 2 CONH 2 CH 3 OH+CO → HCOOCH 3 CO+H 2 O 2 → CO 2 +H 2 O

	-4.6 • 10 -3
	4 • 10 -2

URL: http://mc.manuscriptcentral.com/tandf/tmph

Acknowledgements

The authors gratefully acknowledge the funding of the German Research Council (DFG), which, within the framework of its 'Excellence Initiative', supports the Cluster of Excellence 'Engineering of Advanced Materials' (www.eam.uni-erlangen.de) at the University of Erlangen-Nuremberg.


Hongjun Luo

Wolfgang Hackbusch

Heinz-Jürgen Flad

Quantum Monte Carlo study of the transcorrelated method for correlation factors

Keywords: Transcorrelated method, Jastrow factor, Quantum Monte Carlo

We have studied the performance of Boys and Handy's transcorrelated equation for Jastrow factors, without further approximations on the integrals involved, using quantum Monte Carlo methods. Contrary to previous statements in the literature, we observed no violation of variational bounds. Our results agree to very high accuracy with variational calculations for various atoms and small molecules. It should be mentioned that our findings are not in direct contradiction with the results of Boys and Handy obtained forty years ago. Instead their results were slightly obscured by obvious limitations of computational resources and resulting approximations for complicated integrals. According to our results and related work, the transcorrelated method of Boys and Handy might be an interesting alternative to conventional F12 methods where the correlation factor is not the subject of optimisation.

Introduction

A prevalent prejudice in quantum chemistry is that correlated wave functions of many-electron systems are unwieldy objects which require an enormous number of degrees of freedom for an accurate representation. This point of view originated from large scale configuration interaction (CI) calculations where up to billions of Slater determinants were used in order to achieve accurate results for rather small molecules. A more economical representation of wave functions can be achieved via the exponential ansatz which became popular through its use in coupled cluster (CC) theory. This ansatz incorporates essential physical insights and guarantees size consistency of the method. Another important insight from Hylleraas and collaborators was the astonishing accuracy of expansions in terms of inter-particle distances. A combination of both worlds is the Jastrow ansatz for the wave function

Ψ(R) = e τ (R) Φ(R), R = (r 1 , r 2 , • • • , r N ), (1.1)
where R denotes the combined coordinates of all N particles. The correlation function τ represents a symmetric function of the electron coordinates and the reference wave function Φ(R) can be chosen as a single Slater determinant, e.g. from Hartree-Fock, or an appropriate linear combination. Generally, a given ansatz for the reference wave function has an inherent error in the nodal structure which cannot be compensated by the Jastrow factor. Therefore, the Jastrow ansatz should be considered only as an approximation to the exact wave function. The remaining, so-called fixed-node error is typically small, i.e., less that 5% of the correlation energy, but hard to improve in a systematic manner. Jastrow factors have been extensively studied in condensed matter and nuclear physics [START_REF] Clark | Progress in Nuclear and Particle Physics[END_REF] but rarely considered in quantum chemistry. Remarkable exceptions are the transcorrelated (TC) method of Boys and Handy [START_REF] Boys | [END_REF]3,4,5,6,7,8,9,10] developed forty years ago and more recent developments in quantum Monte Carlo (QMC) methods [11,12,13,14]. Because of the exponential ansatz, the optimisation of the correlation function is a highly nonlinear problem. This can be handled in a variational manner using either QMC [11,15,16,17] or highly sophisticated diagrammatic techniques, like the Fermi hypernetted chain method [18,19]. An alternative nonvariational approach is the TC method of Boys and Handy [4]. The basic idea is to transform the formal Schrödinger equation for the Jastrow ansatz

( Ĥ -E)e τ Φ = 0, (1.2)
into the form (e -τ Ĥe τ -E)Φ = 0, (

where the similarity transformed Hamiltonian gives only rise to quadratic non-linear terms in the correlation function. According to our previous discussion, Eqs. (1.2) and (1.3) are generally not exactly solvable for a given ansatz of the reference wave function Φ. Instead an approximate solution can be envisaged which satisfies

R 3N d 3N R ∂e -τ Φ ∂c α ( Ĥ -E)e τ Φ = 0, (1.4)
for any free parameter c α of the correlation function τ or reference wave function Φ. This transcorrelated formulation strongly resembles the projection equations of CC theory. The error of this approach in the wave function and energy has been characterized by Boys and Handy [START_REF] Boys | [END_REF]4] via the least squares errors of the basis {∂e τ Φ/∂c α } and corresponding projection basis {∂e -τ Φ/∂c α } with respect to the exact wave function. For actual calculations, the original Eq. (1.4) poses severe technical difficulties concerning the computation of the integrals. Three different approaches have been studied by Boys and Handy to circumvent this problem. The first approach is to replace the projection basis by a computationally simpler CI-type basis consisting of a set of single and double excited Slater determinants. This has been tested for the Ne atom [5] and gave a total energy which is more than 20 mHartree below the "exact" energy. Another approach is to fix the leading order term which gives rise to Kato's cusp condition and to choose a separable ansatz for the optimised part of the correlation factor. The difficult to handle leading order term therefore does not show up as a derivative on the left side of the Hamiltonian in (1.4) which greatly simplifies the computations. Applications to LiH and the Be atom are reported in Refs. [6] and [7], respectively. The third approach uses a Gaussian type geminal basis set for the correlation factor which enables the computation of all the integrals involved. Interestingly, it has been observed by Handy that the latter approach behaves similar to a variational calculation [9] in contrast to the other two approaches mentioned before. It is fair to say that the TC calculations of Boys and Handy were of amazing accuracy for those days. However, due to obvious limitations of the computing facilities sometimes very small basis sets were taken which might obscure, together with a mixture of different computational approaches, the actual behaviour of the TC method. It is the purpose of the present work to provide unambiguous results which shed light on the performance of the TC method for Jastrow factor optimisation. The lack of a variational lower bound, generally attributed to the non-hermitian TC Hamiltonian, is usually considered as a severe problem of the TC approach [20]. Therefore it was largely abandoned until the revival of the method by the work of Ten-no and collaborators. In their work [21,22,23,24,25], the correlation factor was kept fixed, e.g., using Gaussian geminals, and the transcorrelated Hamiltonian was treated using perturbation theory or the coupled electron-pair approximation. In this sense Tenno's approach resembles to R12 or F12 methods [20] which improve convergence with the basis set for conventional many-particle methods. Another work along this line is the approach of Umezawa et al. [26,27,28,29] to solve the orbital equations in a self-consistent manner for a single Slater determinant, where the correlation factor is fixed up to a single parameter. Like in the present work, QMC techniques were used to compute the integrals.

Within the present work we want to reconsider the TC method in its original formulation (1.4) using the Hylleraas-like correlation basis introduced by Boys expressions in the inter-particle distances. The complicated integrals, which prevented straightforward computations forty years ago, have been calculated with QMC. For comparison we performed variational Monte Carlo (VMC) calculations with TC and energy optimised Jastrow factors. Originally, the TC equations were designed to optimise the reference function, i.e., the orbitals which constitute the Slater determinant, and the correlation function simultaneously. We refrain from optimising the orbitals and consider only the correlation function. With this no indeterminacy appears concerning the one-electron part of the correlation function which required an additional contraction equation in the work of Boys and Handy [3].

The paper is organized as follows: In Section 2 we give a detailed description of the TC equation and discuss its relations with the variational approach. Test calculations for atoms and small molecules are presented in Section 3 where we compare the performance of the TC and variational approach. Finally, in Section 4 we give some concluding remarks.

Transcorrelated versus variational treatment of electron correlation

For a given reference function Φ, the Jastrow ansatz (1.1) usually does not provide an exact solution for the stationary Schrödinger equation ĤΨ = EΨ (2.1)

with nonrelativistic Hamiltonian Ĥ = - 1 2 i ∇ 2 i + i v ext (r i ) + i>j 1 r ij . (2.2)
in the Born-Oppenheimer approximation. Atomic units have been used throughout the paper. Instead of the exact solution one can search for the correlation function τ (R) which minimises the expectation value of the energy

E = min τ R 3N d 3N R e τ Φ(R) Ĥe τ Φ(R) R 3N d 3N R e τ Φ(R) e τ Φ(R) . (2.3)
Here and in the following we assume all wave functions to be real valued. According to the standard procedure we introduce a basis set for the correlation function

U (1) α (R) = i u (1) α (r i ), U (2)
β (R) = 1 2 i =j u (2) β (r i , r j ), . . . (2.4)
where

U (k)
α denotes symmetric k-particle basis functions. For notational simplicity we skip the index k in the following discussion. Inserting the expansion of the correlation function

τ (R) = α c α U α (R), (2.5) into (2.
3) yields a nonlinear eigenvalue problem

U α e τ ( Ĥ -E)e τ = 0, (2.6)
where denotes the expectation value of an operator Ô with respect to the normalized reference wave function.

Ô := R 3N d 3N R Φ(R) ÔΦ(R) (2.
The variational equation (2.6) can be interpreted as the Schrödinger equation for the Jastrow ansatz projected onto the subspace spanned by {U α e τ Φ}. This means that Schrödinger's equation is only satisfied in a weak sense with respect to the set of test functions {U α e τ Φ} and not on the whole Hilbert space. In contrast to the variational approach, the TC method of Boys and Handy for the correlation factor leads to a weak formulation of the Schrödinger equation of the form

U α e -τ ( Ĥ -E)e τ =
U α ( Ĥ + [ Ĥ, τ ] - 1 2 i (∇ i τ ) 2 -E) = 0, (2.9)
and the TC energy is

E = Ĥ - 1 2 i (∇ i τ ) 2 .
(2.10)

By constant shifts of the basis functions

Ũα = U α -U α , (2.11)
Eq. (2.9) becomes Ũα ( Ĥ + [ Ĥ, τ ] -

1 2 i (∇ i τ ) 2 ) = 0, (2.12)
which is decoupled from the transcorrelated energy (2.10). The non-linear character of Eq. (2.12) requires an iterative solution method. Following Boys and Handy [4], we solve the corresponding linear defect equation Ũα

[ Ĥ, δτ ] - i Ũα (∇ i τ ) • (∇ i δτ ) = -Ũα ( Ĥ + [ Ĥ, τ ] - 1 2 i (∇ i τ ) 2 ) , (2.13)
and update the correlation function τ + δτ → τ of the Jastrow factor in an iterative manner until convergence has been achieved. Some of the integrals in Eq. (2.12) are notoriously difficult to compute and must be treated by elaborate numerical quadrature schemes. In order to simplify the problem Boys and Handy devised an alternative approach where the set of test functions {U α e -τ Φ} is replaced by {e -τ Φ ′′ r } which is more convenient from a computational point of view. Here Φ ′′ r denotes Slater determinants generated via double excitations from the reference determinant Φ. By construction, these test functions, however, have different nodal structures and their span differs from the subspace spanned by {U α e -τ Φ}. Therefore no basis set limit exits, in the sense discussed in the previous paragraph, such that the TC and variational approaches become equivalent. Furthermore, for the new set of test functions, Eq. (2.12) cannot be rigorously satisfied anymore and has to be solved in the sense of a least squares approximation [5].

Even for present day computing facilities, not to mention the situation faced by Boys and Handy forty years ago, TC calculations with the original set of test functions {U α e -τ Φ} are rather challenging. On the one hand basis functions U α are required which provide compact and accurate representations for correlation functions and on the other hand these functions must be also convenient from the computational point of view. In our previous work we studied tensor product wavelets within the framework of diagramatic multiresolution analysis [30] in order to meet both requirements. Quite recently we have also considered canonical tensor product approximations for correlation functions [31] as an alternative approach to reduce the computational complexity of the corresponding integrals. In order to be applicable to the TC method both approaches require extensive software development. The main purpose of the present work is to test the accuracy of the original TC equation (2.12) and to perform a direct comparison with variational calculations. Therefore we have chosen QMC techniques for the computation of the integrals. This enables us to use well established basis functions for the Jastrow factor. The expectation value of a general operator Ô can be rewritten as

Ô = R 3N d 3N R O L (R)ρ(R), (2.14)
with

O L (R) = ÔΦ(R) Φ(R) (2.15)
and measure

ρ(R) = |Φ(R)| 2 R 3N d 3N R |Φ(R)| 2 .
(2.16)

Generating a set of random samples {R m } m=1,M according to the distribution function ρ from a Metropolis algorithm, the expectation value of Ô can be estimated by

Ô ≈ 1 M M m=1 O L (R m ).
(2.17)

The statistical error of this approximation converges with O(M -1 2 ) independent on the number of particles.

Transcorrelated QMC calculations

We have performed some test calculations for the ground states of the atoms Be, O, Ne and the molecules CH 4 , H 2 O. In these calculations, we have used the polynomial basis originally suggested by Boys and Handy [5] which became the standard basis in QMC after its revival by the work of Schmidt and Moskowitz [12]. For the representation of the pair-correlation function this basis is given in terms of polynomials of rational functions in the inter-partical distances, i.e.,

u (2) α (r i , r j ) = (r m iI rn jJ + rm jI rn iJ )r o ij , (3.1)
with riI = r iI 1 + r iI , rij = r ij 1 + r ij , (3.2)
where m, n, o are taken to be non-negative integers and I, J are indices of nuclei. Such basis functions are characterized by a multi-index (m, n, o, I, J) which is simply denoted by α here and in the following. For atoms, the largest basis used in our calculations contains 49 polynomials. It is simply constructed by taking into account all polynomials (3.1) with m + n + o ≤ 6. Among these 49 polynomials, there are 6 one body terms (with m = 0, n = o = 0) and 43 two body terms to which spin dependent coefficients have been assigned except for He atom. Therefore our atomic basis for the pair-correlation function has in total 92 spin dependent degrees of freedom. Wave functions with spin dependent Jastrow factors may not be eigenstates of the total spin operator Ŝ2 . A study by Huang et al. [32], however, revealed that possible spin contaminations are small. For molecules, the multi-index of the basis functions also depends on indices for the nuclei (I, J) which increases their number considerably and requires a more adaptive construction of the basis set. In the case of H . degrees of freedom have been employed in our calculations, respectively. Within the present work, we consider only the TC equation for the correlation factor, and the reference wave functions were taken from Hartree-Fock calculations. For atoms, Slater-type orbitals were taken from Clementi and Roetti [33] and for the molecules Gaussian-type orbitals have been obtained from MOLPRO [34].

As can be expected, the polynomial basis (3.1) leads to ill-conditioned matrices, a problem which is enhanced by statistical fluctuations from the Monte Carlo algorithm. To keep the equations tractable a singular value decomposition (SVD) has been applied to get rid of statistical fluctuations and almost linear dependencies in the basis. The details of our procedure are as follows: usually we take an ensemble of 100 random walkers which undergo random walks of 2 × 10 5 Monte Carlo steps. The matrix elements for the linear TC equation (2.13), in the following briefly denoted A x = b, are calculated for each of the walkers separately, so we have an ensemble of matrices A (i) and vectors b (i) where i refers to the index of the walker. Taking an ensemble average yields the matrix A, and vector b. The SVD applied to the matrix A gives à = U T AV,

where U and V are orthogonal matrices and à is a diagonal matrix with non-negative elements (singular values). Thus the linear equation is diagonalized as

Ãx = b, x = V T x, b = U T b. (3.4)
In order to get an estimate of the statistical error of the matrix à and vector b, we performed the transformations (3. The first test calculations were carried out for the Ne atom. In order to check the dependence of the TC correlation energy on the choice of the basis functions {U α }, we have performed calculations for three different basis sets, i.e., the 8 term basis of Boys and Handy (cf. Table 1. in Ref [5]), the 17 term basis of Schmidt and Moskowitz (cf. Table VI. in Ref [12]) and the 92 term spin dependent basis described before. The convergence of the iterative solution of the transcorrelated equation (2.12) using (2.13) is shown in Fig. 3.1, where the initial correlation function was set to τ = 0. It turned out that the iteration scheme converged very fast, i.e., within two to three iterations for all three basis sets. For comparison, we have also calculated energy optimised correlation functions τ opt for the same basis sets using first-order iterative Jastrow perturbation theory [17]. The corresponding variational correlation energies, which are obtained from VMC calculations of the expectation value of the energy, are indicated in Fig. which can be considered as benchmarks for variational energies in the complete basis set limit. We obtained very close agreement, i.e., within 2% for all three types of correlation energies. In particular, we never observed that (2.10) gives an energy significantly lower than the corresponding variational value. The good agreement between transcorrelated and variationally optimised correlation energies suggests that also the corresponding Jastrow factors are of comparable quality.

Our calculations show a similar qualitative behaviour as the Gaussian-type geminal approach of Handy [9] which can be expected since this approach uses the same TC equation (2.12) and differs only in the choice of the basis set for the correlation function. Based on their calculations for the Be atom and LiH [6,7], it has been argued by Handy [9] that a better convergence behaviour can be expected if the leading order term is not fixed by Kato's cusp condition. We have tested this argument in our calculations and obtained almost the same results independent of whether or not this term was fixed.

  and Handy in terms of polynomials of rational 2

4 Page

 4

Figure 3 . 1 :

 31 Figure 3.1: Convergence of the TC correlation energy (Hartree) of the Ne atom for different basis sets: ( ) 8 term basis of Boys and Handy, ( ) 17 term basis of Schmidt and Moskowitz, (♦) 92 term spin dependent basis. For comparison we have shown the corresponding variational correlation energies for energy optimised correlation functions τ opt from first-order iterative Jastrow perturbation theory (dotted lines). .

  3) and (3.4) for each A (i) and b (i) separately, and calculated the standard error of the matrix elements of à and b based on the fluctuation of Ã(i) and b(i) , respectively. By denoting the diagonal matrix elements of à and their standard errors as Ãα and δ Ãα , respectively, the solution of the linear equation simply becomes xα = bα / Ãα . To prevent the fluctuations being amplified by small singular values, we set xα = 0 if |δ bα | > bα or |δ bα | > Ãα or |δ Ãα | > Ãα .

Figure 3 . 2 :

 32 Figure 3.2: Comparison of TC and variational correlation energies. Ratios with respect to the corresponding fixed-node DMC correlation energies are shown. The variational energies are calculated with TC correlation factors as well as for energy optimised correlation factors from first-order iterative Jastrow perturbation theory. For atoms the 92 term basis has been used whereas for H 2 O and CH 4 the 103 and 125 term basis have been employed, respectively.

3 . 1

 31 as dotted lines. It can be seen that for all three basis sets, the TC and variational correlation energies are in very close agreement to each other. For our other test systems Be, O, H 2 O and CH 4 , calculations have been only performed with the large spin dependent bases discussed before. It is interesting to calculate besides the TC correlation energies E (T C) corr from Eq. (2.10) also the variational correlation energy E (var) corr (τ T C ) for the TC Jastrow factors. The results for TC Jastrow factors have been compared in Table 3.1 with correlation energies for energy optimised Jastrow factors E (var) corr (τ opt ). In Fig. 3.2 these correlation energies are shown for our test systems relative to the fixed-node diffusion quantum Monte Carlo (DMC) results E (DM C) corr

7

 7

  0,(2.8) where a set of test functions {U α e -τ Φ} has been applied. In fact Eqs. (2.6) and (2.8) are closely related since the test functions U α e τ Φ and U α e -τ Φ share the same nodal structure. If we consider the limiting case where {U α } forms a complete basis for n-body permutationally symmetric functions, the spaces spanned by {U α e τ Φ} and {U α e -τ Φ} are identical which means that also Eqs. (2.6) and (2.8) become equivalent.It is only the kinetic energy part of the Hamiltonian which does not commute with the Jastrow factor. Therefore, Eq. (2.8) can be rewritten as

6 Page 6 of 9 URL: http://mc.manuscriptcentral.com/tandf/tmph

				Molecular Physics
			1
			0.98				E corr (TC) (τ TC ) E corr (var) (τ TC )
	(DMC)	corr	0.96				E corr (var) (τ opt )
	corr (τ)/E	0.94
	E
			0.92
			0.9	Be	O	Ne	H 2 O	CH 4
	F o
	r
	P
			e e r
				R e v i e w
							O n
							l y

Table 3 .

 3 1: Transcorrelated and variational correlation energies (mHartree). The variational energies are calculated with VMC for both TC (τ T C ) and energy optimised correlation functions (τ opt ). For comparison, the fixed-node DMC and exact nonrelativistic correlation energies are given. Ref.[35], b Ref.[36], c Ref.[37] 4 ConclusionTest calculations on several atomic and molecular systems demonstrate that the performance of the original TC equation (2.12) is very close to variational calculations for the determination of correlation functions. This important fact was somewhat obscured by the calculations of Boys and Handy owing to technical limitations and related further approximations. Our work together with the results of Ten-no et al. and Umezawa et al. suggests that the TC method is an interesting alternative to conventional F12 methods where the correlation function is kept fixed throughout the calculation. QMC offers an efficient approach to compute complicated integrals which are otherwise difficult to access.

		Be	O		Ne		H 2 O	CH 4
	N term	92	92	8	17	92	103	125
	TC	-0.0769(2)	-0.2253(2)	-0.3370(4) -0.3431(4) -0.3567(3) -0.3290(5) -0.2634(3)
	VMC(τ T C ) -0.0772(2)	-0.2260(5)	-0.3367(5) -0.3446(3) -0.3581(2) -0.3310(4) -0.2657(2)
	VMC(τ opt ) -0.0778(2)	-0.2268(4)	-0.3369(5) -0.3458(4) -0.3599(5) -0.331(1)	-0.2678(5)
	DMC	-0.08415(7) a -0.2422(5) a		-0.3752(3) a		-0.352(1) c -0.283(1) c
	Exact	-0.09434 b	-0.2579 b	-0.3905 b			-0.3719 c	-0.296 c
	a

URL: http://mc.manuscriptcentral.com/tandf/tmph

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics