Correlation energy within exact-exchange ACFD theory: systematic development and simple approximations

Correlation energy within exact-exchange ACFD theory: systematic
development and simple approximations

Nicola Colonna International School for Advanced Studies (SISSA), Via Bonomea 265, I-34136 Trieste, Italy    Maria Hellgren International School for Advanced Studies (SISSA), Via Bonomea 265, I-34136 Trieste, Italy    Stefano de Gironcoli International School for Advanced Studies (SISSA), Via Bonomea 265, I-34136 Trieste, Italy CRS Democritos, CNR-IOM Democritos, Via Bonomea 265, I-34136 Trieste, Italy
July 12, 2019
Abstract

We have calculated the correlation energy of the homogeneous electron gas (HEG) and the dissociation energy curves of molecules with covalent bonds from a novel implementation of the adiabatic connection fluctuation dissipation (ACFD) expression including the exact exchange (EXX) kernel. The EXX kernel is defined from first order perturbation theory and used in the Dyson equation of time-dependent density functional theory. Within this approximation (RPAx), the correlation energies of the HEG are significantly improved with respect to the RPA up to densities of the order of . However, beyond this value, the RPAx response function exhibits an unphysical divergence and the approximation breaks down. Total energies of molecules at equilibrium are also highly accurate but we find a similar instability at stretched geometries. Staying within an exact first order approximation to the response function we use an alternative resummation of the higher order terms. This slight redefinition of RPAx fixes the instability in total energy calculations without compromising the overall accuracy of the approach.

I Introduction

Kohn-Sham (KS) methods that treat the exchange and correlation energy on the basis of the Adiabatic Connection Fluctuation-Dissipation (ACFD) theorem Langreth (); Langreth2 () have raised considerable interest in recent years PhysRevB.65.235109 (); Fuchs (); Jiang (); Hellgren3 (); Hellgren1 (); Nguyen (); range-separated (); PhysRevLett.103.056401 (); PhysRevB.81.115126 (); jcp.132.044109 (); Hellgren4 (); Hebelmann (); SE (); review-RPA (); Thygesen () mainly because they provide a route to overcome the shortcomings of standard local-density-approximation/generalized-gradient-approximation density-functional theory (LDA/GGA DFT). In particular i) an exact expression for the exchange-correlation (xc) energy in term of density-density response function can be derived from the ACFD theorem providing a promising way to develop a systematic improvement for the xc functional; ii) all ACFD methods treat the exchange energy exactly thus canceling out the spurious self-interaction error present in Hartree energy; moreover iii) the correlation energy is fully non-local and automatically includes van der Waals interactions.

The ACFD method is computationally very demanding and most often it is limited to a post self-consistent correction where the xc energy is computed from the charge density obtained from a self-consistent calculation performed with a more traditional xc functional. The basic ingredients needed to compute the correlation energy within the ACFD formalism are the density-density response function of the non-interacting KS system and the density-density response function of a system where the electron-electron interaction is scaled by a coupling constant. While for the former an explicit expression exists, the latter is usually calculated from the Dyson equation of time-dependent density functional theoryLRTDDFT () containing the xc kernel, , that needs to be approximated.

The random phase approximation (RPA) is the simplest approximation; the xc kernel is simply neglected and only the frequency-independent Coulomb or Hartree kernel is taken into account. While correctly describing van der Waals interactions Dobson (); Topics () and static correlation Fuchs (); Hebelmann (); caruso_bond_2013 (), as seen for instance when studying H dissociation, RPA is known to overestimate the correlation energies and thus to poorly describe total energies. Hellgren3 (); Jiang ()

In this respect various approaches have been developed in order to correct the RPA RPA+ (); range-separated (); SE (). A systematic possibility to address the shortcomings of RPA is to include all terms up to a given power of the interaction strength in the kernel. To linear order this implies including not only the Coulomb kernel, defining RPA, but also an exchange contribution. The frequency-dependent exact-exchange kernel, , has been derived by Görling Gorling2 (); Gorling1 (); kim-gorling () from the time-dependent optimized effective potential (TDOEP) method and by Hellgren and von Barth Hellgren1 (); Hellgren2 (); Hellgren3 () from a variational formulation of many-body perturbation theory (MBPT). The corresponding approximation for the density-density response function, named RPAx, is obtained by solving the Dyson equation setting and has been successfully used in the ACFD formula to compute correlation energies of atoms Hellgren1 (); Hellgren4 () and molecules. hesselman_random_2010 (); Hebelmann (); bleiziffer_RI_2012 ()

Here we set the RPAx within the context of a general scheme which allows to formally define a power expansion of the xc kernel combining the general ACFD theory with a many-body approach, specifically the Görling-Levy Perturbation Theory Gorling0 () (GLPT), along the adiabatic-connection path. To first order this reduces to the RPAx for which a novel and efficient implementation based on an eigenvalue decomposition of the interacting time-dependent density response function in the limit of vanishing electron-electron interaction, is proposed.

The performance of the RPAx has in this work been tested on the homogeneous electron gas (HEG) at different values of as well as on the dissociation of diatomic molecules with covalent bonds such as H and N. The results give further support to the accuracy of the RPAx but also reveal an instability or pathological behavior in the low density regime of the HEG and N, which leads to a break-down of the approximation. This breakdown points to the need for including correlation or a screening of the exchange kernel. However, we here show that such a procedure is not always necessary and, in particular, if the aim is to calculate total energies. Instead we reduce the effect of the “bare” particle-hole interaction by omitting all higher order particle-hole scatterings. This can be achieved by expanding the RPA response function in the irreducible polarizability, approximated to first order. In this way we are able to fix the instability and at the same time keep the overall accuracy of the RPAx.

Ii Systematic improvement of the correlation energy

Within the ACFD framework a formally exact expression for the exchange-correlation energy of an electronic system can be derived: Langreth (); Langreth2 ()

(1)

where is the density response function at imaginary frequency, , of a system whose electrons interact with a scaled Coulomb interaction, , and move in a local potential chosen in such a way to keep the electronic density fixed to the ground state density of the fully interacting system (). At , the local potential is equal to the external potential (usually the nuclear potential) of the fully interacting system and coincides with the fully interacting Hamiltonian while at the local potential coincides with the KS potential and is the KS Hamiltonian. For intermediate values of the Hamiltonian of the system is Gorling0 ()

(2)

where is the Hartree potential, is the local exchange potential and is the correlation contribution to the potential. The exchange potential is defined as the functional derivative of the exact KS exchange energy

(3)

that has the same expression as the Hartree-Fock exchange energy but is evaluated with the KS orbitals . It is easy to verify that it can be derived from Eq. (1) replacing with the non-interacting density response function

(4)

where , and are the KS eigenvalues, KS orbitals and occupation numbers, respectively. Subtracting the KS exchange energy from Eq. (1) the correlation energy is obtained in terms of linear density responses:

(5)

where is the Coulomb kernel and is the density response function of the non-interacting KS system. For the interacting density response function can be related to the non-interacting one via a Dyson equation obtained from time-dependent density functional theory (TDDFT):

(6)

where is the scaled frequency-dependent xc kernel. Spatial coordinates dependence is implicit in the matrix notation. When the xc kernel is specified one can thus determine a corresponding correlation energy via Eq. (5).

In the following we will describe a general scheme which allows us to compute the xc kernel to a given order, thus establishing a link between the TDDFT expression for the response function in Eq. (6) and the power expansion of in the interaction strength, which can be obtained resorting to the well established GLPT Gorling0 () along the adiabatic-connection path.

Considering the power expansion for the xc kernel and explicitly expanding the Dyson equation (6) in power of the interaction strength

(7)

it can be seen that the first order kernel, , is intimately related to the first order variation of with respect to and similarly higher order correlation contributions to the kernel are related to the corresponding power in the expansion.

Therefore i) we can define an arbitrarily accurate approximation to the density-density response function considering the expansion of the kernel up to a desired order in :

(8)

where ii) the kernel up to order can be exactly determined by comparing with the expansion of from GLPT and iii) the solution of the Dyson equation for leads to a density-density response function which is exact to order but also contains higher-order terms.

In order to solve the Many Body Hamiltonian in Eq. (2) so as to obtain the xc kernel to a given order in , the xc potential, and hence the xc energy, must be known up to the same level. This apparent circular dependence does not actually hinder the application of the procedure since, thanks to the coupling constant integration involved in Eq. (5), the knowledge of the xc energy, and therefore its functional derivatives, up to order only depends on the xc kernel up to order Our strategy can thus be applied in the sequential way

(9)

showing that to order, i.e. replacing with its non-interacting counterpart , the exact-exchange KS energy is obtained; moving to the next step, the exact-exchange kernel can be derived from first order GLPT and we recover the so-called RPAx approximation for the response function, i.e. , and for the correlation energy, i.e. . Notice that the RPAx correlation energy is exact to order but also contains, although in an approximate way, all higher-order terms, and should not be confused with the perturbative correction to the correlation energy in the Görling-Levy perturbation theory Gorling0 ().

The mathematical complexity of this sequential procedure increases very rapidly and makes extremely hard its application already at the second order; nevertheless the prescription is in principle given. The functional derivative of with respect to the density defines the exact correction to the Hamiltonian in Eq. (2) and allows to apply the GLPT to second order and hence to have access to corresponding second-order contribution to the xc kernel. Solving the Dyson equation with the improved kernel defines a new approximation for the response function which is exact up to second order. Plugging into the ACFD formula (5), leads to a new approximation for the correlation energy, , which is exact to order but also contains, although in an approximate way, all higher-order terms.

Essentially this scheme can be regarded as a revised version of the standard GLPT Gorling0 (); Gorling0-td (); Gorling1 () with the additional step provided by the solution of the Dyson equation for the response function and the calculation of a non-perturbative correlation energy (all order in the coupling-constant appears in and following approximation to ) from the ACFD formula in Eq. (5). In this way we expect this approach to be applicable also to small gap or metallic systems where finite-order many-body perturbation theories break down Gell-Mann (); Giuliani-Vignale ().

Having introduce the general framework, we apply our strategy to first order in the coupling strength, hence we focus on the frequency-dependent exact-exchange kernel and on the calculation of the contribution to the correlation energy (previously denoted as RPAx Hellgren1 (); Hellgren4 () or EXXRPA bleiziffer_RI_2012 (); Hebelmann ()) for which we propose a novel and efficient implementation.

Iii Efficient calculation of RPAx correlation energy

Our implementation for computing the RPAx correlation energy is based on an eigenvalue decomposition of the time-dependent response function in the limit of vanishing coupling constant. The scheme described below is a generalization of the implementation proposed by Nguyen and de Gironcoli Nguyen () for computing RPA correlation energies.

iii.1 RPAx Correlation Energy

Let us start by defining the following generalized eigenvalue problem:

(10)

where the eigenpairs and all the operators implicitly depend on the imaginary frequency . Once the solution of the generalized eigenvalue problem (10) is available, the trace in Eq. (5) is simply given by

(11)

and the integration over the coupling constant can be calculated analytically, leading to the final expression

(12)

Notice that Eq. (10) and (12) demonstrate that knowledge of is sufficient for computing the RPAx correlation energy and the exact-exchange kernel alone is not needed.

iii.2 Exact-Exchange kernel

The exact expression for in term of the KS eigenvalues and eigenfunctions has been derived by Görling starting from the time-dependent optimized potential method equationGorling2 () and by Hellgren and von Barth starting from the variational formulation of many-body perturbation theoryHellgren1 (); Hellgren3 (). Here we propose an alternative derivation staying within the general scheme described in the previous section.

In section II it has been shown that is the first order correction to the non-interacting response function due to the switching on of the perturbation . Moreover in the previous sub-section it has been shown that the eigenvalues and eigenvectors of are sufficient for computing RPAx correlation energies. In what follows we derive the exact expression for the matrix elements of in term of the KS eigenvalues and eigenfunctions and their first order corrections only, and show how they can be efficiently computed resorting to the linear-response techniques of density functional perturbation theory DFPT ().

Let us start by considering the matrix element of on two arbitrary, and , time-dependent perturbing potentials at imaginary frequency

(13)

For a non degenerate ground state the linear response density at imaginary frequency can be written as

(14)

where are the first order corrections to the KS wavefunction due to the perturbation and satisfy the linearized time-dependent KS equations

(15)

Eq. (13) becomes and if the (static) perturbation is turned on, the first order correction to , i.e. , in the coupling constant can be computed:

(16)

where is obtained by taking the linear variation of Eq (15)

(17)

while the static correction vector satisfies the linearized time-independent Schrödinger equation

(18)

with .

With a simple manipulation it’s easy to show that depends only on the GS wavefunction and its first order corrections (and not on the second order correction ). Taking the h.c. of Eq. (15) and multiplying it on the right by and Eq. (17) on the left by and subtracting the two identities so obtained, an expression for is obtained where the second order corrections cancel out.

The final expression for becomes:

(19)
(20)
(21)

Eq. (LABEL:eq2.12) together with Eq. (15) and Eq. (18) defines the matrix elements as a function of the KS many-body ground-state wavefunctions and its first order corrections and . Introducing their definitions in terms of the single particle KS orbitals, ’s, and their first order variations, ’s and ’s, Eq. (LABEL:eq2.12) becomes:

(23)

where the sums run over the occupied single-particle KS state only and and are the (conduction-band projected) variations of the occupied single-particle state. They can be efficiently computed resorting to the linear-response techniques of density functional perturbation theory DFPT ():

(24)

where is the non-local exchange operator identical to the Hartree-Fock one but constructed from KS orbitals, is the projector on the occupied manyfold and is a positive constant larger than the valence bandwidth in order to ensure that the linear system is not singular even in the limit for .

Inserting the formal solutions for and from Eq. (24) into Eq. (23), removing the trial perturbing potential and sending , the expression for previously derived in Ref. Gorling2, ; bleiziffer_RI_2012, and in Ref. Hellgren2, is recovered.

The scheme described above has been implemented in the Quantum ESPRESSO distribution QE (). The basic operations involved in the calculation of the matrix elements are the same required for the calculation of RPA energy and potential in the implementations proposed by Nguyen and de Gironcoli Nguyen () and Nguyen et al. Linh () respectively, meaning that our RPAx calculation has a computational cost comparable to their RPA implementations and maintains their favorable scaling.

Iv Homogeneous electron gas

As a test of the accuracy of the new approximation we choose the simple homogeneous electron gas. The homogeneous electron gas is an idealized system of electrons moving in a uniform neutralizing background. At zero temperature it is characterized by two parameters only, i.e. the number density , or equivalently the Wigner-Seitz radius , and the spin polarization , where is the density of spin up (down) electrons and . Despite its simplicity, (i) the HEG model represents the first approximation to metals where the valence electrons are weakly bound to the ionic cores, (ii) the system is found to display a complex phase diagram including transition to the Wigner crystal and in addition (iii) it provides the basic ingredient of any practical density functional calculation. The most widely used approximations for the unknown xc-energy functional are based on properties of the HEG.

iv.1 Unpolarized HEG

We begin by studying the unpolarized HEG. While the solution of Dyson equation is demanding in general, it becomes trivial in the case of the HEG; the response functions and the kernels are all diagonal in momentum space and the RPAx Dyson equation can be easily solved as

(25)

where and is the exchange kernel at a given momentum and frequency.

The correlation energy per electron follows from Eq. (5) where the trace has been replaced by an integral over momentum and the integration over has been done analytically

(26)

Here has been defined as

(27)

While the Lindhard function at imaginary frequency is known exactly Giuliani-Vignale (), the function can be directly derived from the general expression given in Eq. (23) and is given by a six-fold integral over crystal momenta. Its static values was computed first numerically by several author Geldart (); Antoniewicz (); Chevary () and later analytically by Engel and Vosko EV ().The frequency dependence of has been calculated by Bronsens, Lemmens and Devreese Bronsens1 (); Bronsens2 () for real frequencies and by Richardson and Ashcroft Richardson () for imaginary frequencies. Following Bronsens et al. four integrations can be done analytically using cylindrical coordinates; we used numerical quadrature for the two remaining integrations. Our numerical integration is able to recover the analytic results of Engel and Vosko EV () in the limit . Finally the integration over momentum and imaginary frequency in Eq. (26) has been computed numerically. The results are listed in Table 1 and Fig. 1. RPA can be easily obtained from Eq. (26) and Eq. (27) with and can be seen to seriously overestimate the correlation energy at all densities. Including the exact exchange kernel greatly improves over simple RPA and the RPAx correlation energy per particle is close to the accurate Quantum Monte Carlo (QMC) results Ceperly ().

Figure 1: Correlation energy per particle in the homogeneous electron gas as a function of the Wigner-Seitz radius evaluated with different kernels: RPA (red squares), RPAx (green diamonds) and QMC calculation (black circles).

As expected RPAx works well for small values of and becomes less accurate when increases. According to our calculation, within RPAx for there is a charge density instability with wave-vector . In Fig. 2 the critical behavior of the static density-density RPAx response function is shown for the full interacting system (Eq. (25), ). When the density decreases a pronounced peak appears at indicating the instability with respect to charge modulations with this wave-vector. As can be seen from the inset in Fig. (2), for sufficient large values or , approaches unity and the denominator in Eq. (25) tends to vanish leading to the appearance of the peak. Beyond , exceeds unity and RPAx approximation breaks down as the density-density response function is not anymore negative definite.

This instability resembles the charge density wave instability, already observed at the Hartree-Fock level by Overhauser Overhauser (); Giuliani-Vignale () and is an artifact of the truncation of the kernel expansion to first order in the interacting strength. A full treatment of correlation in the QMC calculations moves the density instability toward the Wigner crystal to much smaller densities corresponding to Ceperly ()

RPA RPAx QMC
0.5 -0.194 -0.154 -0.153
1.0 -0.157 -0.121 -0.119
3.0 -0.105 -0.077 -0.074
5.0 -0.085 -0.060 -0.056
8.0 -0.068 -0.047 -0.043
10.0 -0.061 -0.042 -0.037
11.0 -0.058 -0.035
Table 1: Corelation Energy per particle with different kernels: RPA (), RPAx (), and Quantum Monte Carlo calculation Ceperly (). All energies are in Rydberg.
Figure 2: Critical behavior of the static density-density RPAx response function when the density decreases. For the system becomes unstable respect to charge modulation with wavevector .

iv.2 Alternative RPAx resummations

In Sect. II we have established a strategy for a systematic improvement of the xc kernel. However, because of the complexity of the procedure, rather than proceeding along this way we propose here two simple modifications to the original RPAx approximation which are able to fix the instability problem and, at the same time, to give correlation energies on the same level of accuracy as RPAx.

Introducing the irreducible polarizability , it is possible to write the interacting response function asGiuliani-Vignale ()

(28)

where . Neglecting and summing up to infinite order leads again to the RPAx approximation defined above. If we instead replace with only its first order expansion we can define a new approximation, here named tRPAx, which contains only a subset of the original RPAx expansion:

(29)

with . In this way we are only including terms which contain first order particle-hole interactions.

Figure 3: Correlation energies per particle as a function of evaluated from the RPA, original and modified RPAx response functions and compared to accurate QMC calculations.

A similar idea has been proposed in Ref. RPAx-F, where the authors suggest to expand the TDDFT response function in a power series of the RPA response function (instead of the non-interacting one), and then to keep only the first order. This amounts in a alternative re-summation, here named tRPAx, for the interacting response function:

(30)

We notice that tRPAx and t’RPAx both only require to be defined. Both approximations thus neglect all higher order particle-hole scatterings which in the original RPAx are simulated by the kernel.

Up to first order the alternative RPAx response functions coincide with the original one, while they have different power expansions starting from the term, meaning that only contributions already approximated at the RPAx level are affected by these different re-summations.

Fig. 3 shows the correlation energies per particle obtained starting from the alternative RPAx approximations of the response function. As expected, for high density electron gases (small values of ) the correlation energies are essentially identical to the original one, since the underlying response functions are the same in the limit for . At the same time they are well behaved also where the original RPAx approximation breaks down.

In Fig. 4 we compare the corresponding static density response functions (calculated at full interaction strength ) with the exact one, obtained from QMC calculation moroni_static_1995 (), for a density corresponding to . The difference between RPA and QMC results reveals that exchange and correlation effects in the kernel are important already at this density; including the exact-exchange kernel (original RPAx) overcorrects the RPA deficiency, in particular between and , while both the alternative RPAx approximations give a much better agreement with accurate QMC calculations. Thus despite the fact that the RPAx energy is better at this value of the static response function is worse suggesting that the RPAx results are subjected to a cancellation of errors when integrated over the frequency.

Figure 4: Approximate static response functions for the HEG at as compared to the exact QMC calculations moroni_static_1995 (). The alternative RPAx approximations show a better agreement with QMC results.

In the range of densities analyzed, tRPAx and tRPAx response functions do not show any critical behavior; moreover when the density decreases a trend opposite to the one found for the RPAx response function is observed with a reduction (instead of the enhancement shown in Fig. 2) of the height of the peak near , suggesting no divergence would appear even for smaller densities.

iv.3 Spin-polarized HEG

We continue our analysis of the HEG at the RPAx level by studying the spin magnetization dependence of the correlation energy of the system. We start noticing that for the non-interacting system the spin-up and spin-down components of the gas are independent so that a simple scaling relation between the non-interacting density-density response functions of the polarized and unpolarized gas can be derived:

(31)

while .

The spin-up and spin-down components behave as independent constituents of the system at the exchange level too and a scaling relation similar to Eq (31) holds true also for the exchange energy exchange-spin () and, accordingly, for the exchange potential and kernel:

(32)

while . Thus at the RPAx level the interaction between the spin-up and spin-down components of the system is only mediated by the Coulomb kernel .

Although more involved than for the unpolarized case, the solution of the RPAx Dyson equation for the polarized gas is anyway straightforward and using the definitions in Eq. (31) and Eq (32), the RPAx response function of the polarized HEG can be written as:

(33)

where and are the same functions already used for the unpolarized case but evaluated at density or .

Figure 5: Spin polarization function for from RPA (red line), RPAX (green squares), Perdew-Zunger parametrization Perdew-Zunger () (blue solid line), Perdew-Wang parametrization Perdew-Wang () (brown dashed line) and Quantum Mote Carlo calculation Ortiz () (open circles).

Integrating Eq. (5) with the new definition of in Eq. (33), gives the correlation energy per particle, , as a function of and or, equivalently, as a function of and . At the RPA level the dependence of the correlation energy on the spin magnetization has been already calculated long time ago by Von Barth and Hedin vBH () and more recently by Vosko, Wilk, and Nusair VWN (). Our RPA results, simply obtained by setting in Eq. (33), are, within the numerical accuracy, in perfect agreement with both the above mentioned calculations. Fig. 5 shows the spin-polarization function defined as

(34)

for the case evaluated at the RPA and RPAx level, and compares it with the exchange-only dependence that is the one assumed in the Perdew-Zunger parametrization Perdew-Zunger () of Local Spin Density Approximation (LSDA) and with the Perdew-Wang parametrization Perdew-Wang (), that is based on the more physically motivated spin-interpolation expression proposed by Vosko, Wilk, and Nusair VWN (). While within RPAx the correlation energy significantly improves with respect to RPA results, there is essentially no difference between the RPA and RPAx spin-polarization functions. For this value of , calculations done with the alternative re-summations (tRPAx and tRPAx) give essentially the same results as the original RPAx and are not shown in Fig. (5). Thus for this property of the system RPA and all the RPAx (original and alternative) approximations give results in very good agreement with accurate Quantum Monte Carlo calculationsOrtiz () performing much better than the Perdew-Zunger parametrization and slightly better than the more sophisticated Perdew-Wang parametrization.

V Bond dissociation of dimers

As a second test for the RPAx approximation we studied the dissociation curve of the hydrogen and nitrogen molecules.

Figure 6: Dissociation curve of H molecule. PBE, RPA and RPAx (original and alternative) results are compared with accurate calculations kolos_potential_1965 ().

Within standard Density Functional Approximations (DFAs) the proper (singlet) KS ground state of these molecules at large interatomic separations has too high total energy (as illustrated later in Figs. 6 and 7). A better agreement with the experimental potential energy curve can be achieved resorting to a spin-polarized calculation which give good energies, however at the price of a qualitatively wrong spin-density. In a spin-unrestricted calculation, beyond a certain value of the interatomic separation the two spin components, defining the total electron density are no longer equal leading to solution which is no more a singlet as it should.

The H and N dissociation curves at the RPA level has been previously studied Fuchs (); caruso_bond_2013 (); furche_molecular_2001 (). In Refs. Fuchs, ; rohr, the authors have shown RPA to be size-consistent, and thus to correctly describe the dissociation without resorting to any artificial spin-symmetry breaking. However the total energy is far too negative because of the well know overestimation of the correlation energy RPA+ (). Moreover an erroneous repulsion bump appears in the dissociation curves at intermediate distances.

Recently Heßelmann et al. reported the H dissociation curve within the RPAx approximation showing very good result for the total energy both around the equilibrium position and at dissociation but still the problem of the unphysical bump at intermediate bond-lengths remains. Görling and coworkers have also computed RPAx total energies for a set of 21 molecules but always in their equilibrium geometries hesselman_random_2010 (); bleiziffer_RI_2012 ().

Figure 7: Dissociation curve of N molecule. The plot compares results from PBE, RPA and RPAx (original and alternative) calculations.

Here we would like to asses the performance of the RPAx (original and alternative) approximations for molecules beyond their equilibrium geometries studying the dissociation curves of H and N.

The dimers and the corresponding isolated atoms were simulated using a simple-cubic super cell with a size length and bohr, respectively. A kinetic energy cut-off of Ry were used for both systems and up to 200 lowest-lying eigenpairs of the generalized-eigenvalue problem in Eq. (10) were used to compute the RPA and RPAx correlation energies. All the calculations have been done starting from well converged PBE orbitals.

PBE RPA RPAx tRPAx tRPAx Refs.
H
(Å) 0.755 0.740 0.738 0.742 0.738 0.741
(eV) 6.78 4.85 4.41 4.48 4.45 4.75
(cm) 4219 4520 4560 4506 4406 4529
N
(Å) 1.102 1.100 1.085 1.090 1.086 1.095
(eV) 16.86 9.92 9.22 9.07 9.91
(cm) 2274 2322 2569 2430 2482 2383
Table 2: Equilibrium properties of hydrogen and nitrogen dimers computed within different functionals: PBE, RPA and all the RPAx. Accurate values extracted from dissociation curves from Ref. kolos_potential_1965, for H and from Ref. gdanitz_accurately_1998, for N are also given. Equilibrium bond length () in Å  binding energy () in meV , and vibrational frequency () in cm

In Figs. 6 we report our results for the dissociation curves of H and in Tab 2 the structural parameters extracted from them. Comparison with accurate calculations kolos_potential_1965 () illustrates the aforementioned deficiencies of PBE and RPA dissociation curves: standard DFAs give too high total energy in the dissociation limit while RPA overestimates the correlation energy leading to a curve well below the reference one. Including the exact-exchange kernel leads to a sensible improvement in the total energy description; as can be seen from the inset in Fig. 6 the RPAx total energies around the equilibrium position is in very good agreement with accurate quantum chemistry calculations. The alternative re-summations while essentially giving the same energy as the original RPAx in the minimum region, have a positive effect on the dissociation curve at intermediate distances reducing the height of the repulsive hump. We notice that at large interatomic separations all the RPAx approximation drops below the exact dissociation limit of Ry in agreement with the analysis reported in Ref. Fuchs, .

With the simple H example in mind we can turn to analyze the more interesting case of the N molecule. In Fig. 7 we report our results for the dissociation curve and in Tab. 2 the structural parameters obtained from them. As already observed for the H dimer also in this case the whole RPA dissociation curve lies far below all the other curves. Nevertheless the structural parameters at the RPA level are in very good agreement with results from accurate quantum chemistry calculations gdanitz_accurately_1998 (). Including the exact-exchange contribution to the kernel corrects for the RPA overestimation of the correlation energy shifting the RPAx dissociation curve upward. At the same time, the good performance for the equilibrium bond length and the vibrational frequency already obtained at the RPA level is maintained. However, unlike what happens for the H molecule, in this case the original RPAx approximation breaks-down when the nitrogen atoms are separated. For bond lengths greater than Å the RPAx response function is no more negative-definite leading to an instability which is very similar the one observed for the low-density homogeneous electron gas and, ultimately, causes the break-down of the approximation. The alternative re-summations proposed to fix the pathological behavior of the RPAx response function in the HEG, turn out to be effective also in this very different situation. The tRPAx and tRPAx dissociation curves are close to the RPAx one in the equilibrium region (see the inset in Fig. 7) but they are well-behaved also for bond lengths greater that Å overcoming, also in this case, what appears to be an intrinsic inadequacy of the original RPAx approximation.

Vi Conclusions

In this work, we have setted the RPAx approximation for the correlation energy within a general scheme that combines the general framework of the ACFD theory with a systematic many-body approach along the adiabatic-connection path and allows in principle to improve the xc kernel for the purpose of calculating increasingly more accurate correlation energy. We have shown that, in a perturbative approach, RPA is an “incomplete” approximation and that the exact-exchange kernel has to be taken into account for a consistent description to first order in the interaction strength. An efficient method for the calculation of the RPAx correlation energy has been proposed, based on an eigenvalue decomposition of the time-dependent response function of the many body system in the limit of vanishing coupling constant.

The accuracy of the RPAx approximation has been tested on the homogeneous electron gas revealing a great improvement over RPA results and a very good agreement with accurate QMC calculations. The spin magnetization dependency of the RPA and RPAx correlation energies has been calculated as well, showing a big improvement if compared to standard parametrization and a nearly perfect agreement with QMC calculation.

These encouraging results are however disturbed by the break-down of the procedure for large values of where the RPAx density-density response function unphysically changes sign thus indicating that correlation contributions to the kernel are needed to obtain accurate results for the HEG at low densities. Staying within an exact first order approximation to the particle-hole interaction we have suggested two simple and inexpensive modifications of the RPAx approximation which lead to a good description of the correlation energy of the system even in the limit of small densities.

We then examine molecular dissociation of H and N within the RPAx approximation, discovering the same virtues and vices already observed in the HEG case. A sensible improvement of the total energy description is disturbed by a pathological behavior of the response function which ultimately poses doubts on the broad applicability of the RPAx approximation. The alternative re-summations, tRPAx and tRPAx, proposed here, have been shown to be able to fix the RPAx inadequacy without compromising its virtues. Although more tests are needed in order to completely characterize them, tRPAx and tRPAx emerge as promising and stable alternatives to the original RPAx approximation.

Acknowledgements.
We thank CINECA (Bologna, Italy) and SISSA (Trieste, Italy) for the availability of high performance computing resources and support.

References

  • (1) D. Langreth and J. Perdew, Solid State Communications 17, 1425 (1975).
  • (2) D. C. Langreth and J. P. Perdew, Phys. Rev. B 15, 2884 (1977).
  • (3) M. Fuchs and X. Gonze, Phys. Rev. B 65, 235109 (2002).
  • (4) M. Fuchs, Y.-M. Niquet, X. Gonze, and K. Burke, The Journal of Chemical Physics 122, (2005).
  • (5) H. Jiang and E. Engel, The Journal of Chemical Physics 127, (2007).
  • (6) M. Hellgren and U. von Barth, Phys. Rev. B 76, 075107 (2007).
  • (7) M. Hellgren and U. von Barth, Phys. Rev. B 78, 115107 (2008).
  • (8) H.-V. Nguyen and S. de Gironcoli, Phys. Rev. B 79, 205114 (2009).
  • (9) J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. Ángyán, Phys. Rev. Lett. 102, 096404 (2009).
  • (10) J. Harl and G. Kresse, Phys. Rev. Lett. 103, 056401 (2009).
  • (11) J. Harl, L. Schimka, and G. Kresse, Phys. Rev. B 81, 115126 (2010).
  • (12) H.-V. Nguyen and G. Galli, The Journal of Chemical Physics 132, (2010).
  • (13) M. Hellgren and U. von Barth, The Journal of Chemical Physics 132, (2010).
  • (14) A. Heßelmann and A. Görling, Phys. Rev. Lett. 106, 093001 (2011).
  • (15) X. Ren, A. Tkatchenko, P. Rinke, and M. Scheffler, Phys. Rev. Lett. 106, 153003 (2011).
  • (16) X. Ren, P. Rinke, C. Joas, and M. Scheffler, Journal of Materials Science 47, 7447 (2012).
  • (17) T. Olsen and K. S. Thygesen, Phys. Rev. Lett. 112, 203001 (2014).
  • (18) E. K. U. Gross and W. Kohn, Phys. Rev. Lett. 55, 2850 (1985).
  • (19) J. Dobson, Topics in Condensed Matter Physics, Nova, 1994.
  • (20) E. Gross, J. Dobson, and M. Petersilka, Topics in Current Chemistry, volume 181, Springer-Verlag, 1996.
  • (21) F. Caruso, D. R. Rohr, M. Hellgren, X. Ren, P. Rinke, A. Rubio, and M. Scheffler, Phys. Rev. Lett. 110, 146403 (2013).
  • (22) S. Kurth and J. P. Perdew, Phys. Rev. B 59, 10461 (1999).
  • (23) A. Görling, International Journal of Quantum Chemistry 69, 265 (1998).
  • (24) A. Görling, Phys. Rev. A 57, 3433 (1998).
  • (25) Y.-H. Kim and A. Görling, Phys. Rev. B 66, 035114 (2002).
  • (26) M. Hellgren and U. von Barth, The Journal of Chemical Physics 131, 044110 (2009).
  • (27) A. Heßelmann and A. Görling, Mol. Phys. 108, 359 (2010).
  • (28) P. Bleiziffer, A. Heßelmann, and A. Görling, The Journal of Chemical Physics 136, (2012).
  • (29) A. Görling and M. Levy, Phys. Rev. B 47, 13105 (1993).
  • (30) A. Görling, Phys. Rev. A 55, 2630 (1997).
  • (31) M. Gell-Mann and K. A. Brueckner, Phys. Rev. 106, 364 (1957).
  • (32) G. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid, Cambridge University Press, 1971.
  • (33) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  • (34) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (19pp) (2009).
  • (35) N. L. Nguyen, N. Colonna, and S. de Gironcoli, Phys. Rev. B 90, 045138 (2014).
  • (36) D. Geldart and R. Taylor, Can. J. Phys. 48, 155 (1970).
  • (37) P. R. Antoniewicz and L. Kleinman, Phys. Rev. B 31, 6779 (1985).
  • (38) J. A. Chevary and S. H. Vosko, Phys. Rev. B 42, 5320 (1990).
  • (39) E. Engel and S. H. Vosko, Phys. Rev. B 42, 4940 (1990).
  • (40) F. Bronsens, L. Lemmens, and J. Devreese, Phys. Stat. Sol. (b) 74, 45 (1976).
  • (41) F. Bronsens, L. Lemmens, and J. Devreese, Phys. Stat. Sol. (b) 80, 99 (1977).
  • (42) C. F. Richardson and N. W. Ashcroft, Phys. Rev. B 50, 8170 (1994).
  • (43) D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
  • (44) A. W. Overhauser, Phys. Rev. 128, 1437 (1962).
  • (45) J. E. Bates and F. Furche, The Journal of Chemical Physics 139, (2013).
  • (46) S. Moroni, D. M. Ceperley, and G. Senatore, Phys. Rev. Lett. 75, 689 (1995).
  • (47) G. L. Oliver and J. P. Perdew, Phys. Rev. A 20, 397 (1979).
  • (48) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  • (49) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  • (50) G. Ortiz and P. Ballone, Phys. Rev. B 50, 1391 (1994).
  • (51) U. von Barth and L. Hedin, Journal of Physics C: Solid State Physics 5, 1629 (1972).
  • (52) S. H. Vosko, L. Wilk, and M. Nusair, Canadian Journal of Physics 58, 1200 (1980).
  • (53) W. Kolos and L. Wolniewicz, The Journal of Chemical Physics 43, 2429 (1965).
  • (54) F. Furche, Phys. Rev. B 64, 195120 (2001).
  • (55) M. Hellgren, D. R. Rohr, and E. Gross, J. Chem. Phys 136, 034106 (2012).
  • (56) R. J. Gdanitz, Chemical Physics Letters 283, 253 (1998).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
114368
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description