Accurate ground state- and quasiparticle energies: beyond the RPA and GW methods with adiabatic exchange-correlation kernels

We review the theory and application of adiabatic exchange-correlation (xc-) kernels for ab initio calculations of ground state energies and quasiparticle excitations within the frameworks of the adiabatic connection fluctuation dissipation theorem and Hedin's equations, respectively. Various different xc-kernels, which are all rooted in the homogeneous electron gas, are introduced but hereafter we focus on the specific class of renormalized adiabatic kernels, in particular the rALDA and rAPBE. The kernels drastically improve the description of short-range correlations as compared to the random phase approximation (RPA), resulting in significantly better correlation energies. This effect greatly reduces the reliance on error cancellations, which is essential in RPA, and systematically improves covalent bond energies while preserving the good performance of the RPA for dispersive interactions. For quasiparticle energies, the xc-kernels account for vertex corrections that are missing in the GW self-energy. In this context, we show that the short-range correlations mainly correct the absolute band positions while the band gap is less affected in agreement with the known good performance of GW for the latter. The renormalized xc-kernels offer a rigorous extension of the RPA and GW methods with clear improvements in terms of accuracy at little extra computational cost.


Introduction
For decades density functional theory (DFT) has been the workhorse of first-principles materials science. Immense efforts have gone into the development of improved exchange-correlation (xc)-functionals and today hundreds of different types exists, including the generalized gradient approximations (GGA), meta GGAs, (screened) hybrid functional, Hubbard corrected functionals (LDA/GGA+U), and the non-local van der Waals density functionals. Typically, these contain several parameters that have been optimized for a particular type of problem or class of materials. Moreover they rely on fortuitous and poorly understood error cancellations. This limits the universality and predictive power of commonly applied xc-functionals and the accuracy is often highly system dependent.
At the highest rung of the current hierarchy of xc-functionals, lie those based directly on the adiabatic-connection fluctuation-dissipation theorem (ACFDT). The ACFDT provides an exact expression for the electronic correlation energy in terms of the interacting density response function [1,2]. An attractive feature of the ACFDT is that it provides the pure correlation energy, which should then be combined with the exact exchange energy. This clear division removes the reliance on error cancellation between the exchange and correlation terms, which is significant (and uncontrolled) in the lower rung xc-functionals. A further advantage of the ACFDT is that, even in its simplest form, it captures dispersive interactions very accurately through the non-locality of the response function.
The simplest approximation to the response function beyond the non-interacting one is the random phase approximation (RPA). The RPA generally provides an excellent account of longrange screening and it cures the pathological divergence of second-order perturbation theory for the homogeneous electron gas. However, an important shortcoming of the RPA response function is that the local (r close to r ) correlation hole derived from it, via the ACFDT, is much too deep leading to a drastic overestimation of the absolute correlation energy by several tenths of an eV per electron. This key observation is responsible for most of the failures of the RPA and GW schemes to be discussed in this review. It occurs because the RPA response function only accounts for the Hartree component of the induced potentials. The neglected xc-component of the induced potential is short range in nature and therefore mainly influences the local shape of the correlation hole.
Early work [3][4][5] applied the ACFDT-RPA to compute the dissociation energies of small molecules finding a systematic tendency of the RPA to underbind and generally lower accuracy than the generalized gradient approximations (GGA). It was also demonstrated that RPA accounts well for strong static correlation and correctly describes the dissociation curve of the N 2 molecule. Around ten years later, the RPA was applied to calculate cohesive energies of solids [6] again finding that RPA performs significantly worse than GGA with a systematic tendency to underbind. In contrast, RPA was found to produce excellent results for structural parameters of solids [7,8] as well as bond energies in van der Waals systems like graphite [9] and noble gas solids [10], which are poorly described by semi-local approximations. In addition, for the case of graphene adsorbed on metal surfaces, where dispersive and covalent interactions are equally important, the RPA seems to be the only non-empirical method capable of providing correct potential energy curves [11,12]. While the RPA method has many attractive features, it is clear that its poor description of short-range correlations, which results in overestimation of absolute correlation energies, and underestimation of covalent bond strengths, disqualifies it as the highly desired fully ab initio and universally accurate total energy method. One approach to improving the RPA total energy is based on the idea of correcting the RPA self-correlation energy by including higher order exchange terms in many-body perturbation theory (MBPT) and is referred to as second-order screened exchange (SOSEX). The SOSEX correlation energy vanishes for one-electron systems and improves the accuracy of covalent bonds slightly compared to RPA. However, it deteriorates the good description of static correlation and barrier heights in the RPA [13,14]. In addition, SOSEX scales as N 5 with system size and therefore comprises a significant computational challenge compared to RPA, which scales as N 4 .

The adiabatic connection fluctuation-dissipation theorem
In the Kohn-Sham (KS) scheme, a non-interacting Hamiltonian is constructed such that it has a ground state Slater determinant |ϕ 0 , which yields the same ground state density as the true ground state wavefunction |ψ 0 . The adiabatic connection denotes a generalization of this scheme where the Coulomb interaction v c is rescaled by λ, such that the ground state wavefunction |ψ λ 0 reproduces the electronic density of |ψ 0 . The procedure can be accomplished by modifying the external potential v λ KS (r) and we have that |ψ λ=1 0 = |ψ 0 and |ψ λ=0 0 = |ϕ 0 . The adiabatic connection allows one to obtain a highly useful expression for the correlation energy. To begin with the Hartree-exchange-correlation energy can be written as [39] where E tot is the total electronic ground state energy, T KS is Kohn-Sham kinetic energy, and v ext is the expectation value of the external potential. In the last quality we used the Hellmann-Feynman theorem and the fact that ψ λ is defined as the state that minimizes the expectation value of T + λv c . Inserting the second quantized form of the Coulomb interaction the expression becomes wheren(r) = Ψ † (r)Ψ(r) and we have introduced the notation . . . λ = ψ λ 0 | . . . |ψ λ 0 . Since the last term is independent of λ, we can get rid of it by subtracting the Hartree-exchange energy E Hx , which is given by a similar expression with |ψ λ 0 replaced by |ϕ 0 . We then have E c = 1 2 1 0 dλ drdr |r − r | n(r)n(r ) λ − n(r)n(r ) 0 The density-density correlation function is closely related to the density-density response function. The retarded response at vanishing temperature is defined by χ λ (r, r ; t, t ) = −iθ(t − t ) [n(r, t),n(r , t )] λ , where the expectation value is with respect to the ground state. In the frequency domain it becomes χ λ (r, r ; ω) = m =0 n λ 0m (r)n λ m0 (r ) ω − E m0 + iη − n λ 0m (r )n λ m0 (r) where n 0m (r) = ψ λ 0 |n(r)|ψ λ m , E λ m0 = E λ m − E λ 0 are the eigenvalue differences, and η is a positive infinitesimal. It is then clear that −1 π ∞ 0 dωImχ λ (r, r ; ω) = m =0 n λ 0m (r)n λ m0 (r ) = n(r)n(r ) λ − n(r)n(r ) = δn(r)δn(r ) λ , with δn(r) ≡n(r) − n(r). The equality is an example of a fluctuation-dissipation theorem, since it relates the imaginary (dissipative) part of the density response to the correlation between density fluctuations. The retarded response function only has poles in the negative imaginary half-plane and its frequency integral on a closed loop in the upper right quarter of the complex plane vanishes since χ ∼ 1/|ω| 2 for |ω| → ∞. We can thus switch the integration path to the positive imaginary axis where the frequency dependence is smooth. Noting that χ * (r, r ; iω) = χ(r , r; iω) we obtain where ⟪. . .⟫ indicates the trace of the two-point functions involved in the adiabatic-connection integrand.
The problem of calculating the correlation energy has thus been rephrased into finding a good approximation for the density-density response function. The simplest non-trivial approximation is the Random Phase Approximation (RPA), which can be obtained from many-body perturbation theory by assuming a non-interacting irreducible response function. Alternatively, the full (reducible) response function can be obtained from Time-Dependent Density Functional Theory (TDDFT), where it can be shown to satisfy the Dyson equation where all quantities are functions of r and r and integration of repeated variables is implied. f xc (ω) is the temporal Fourier transform of the exchange-correlation kernel and any approximation to f xc thus implies an approximation for the ground state correlation energy in the framework of the adiabatic-connection combined with the fluctuation dissipation theorem. In the context of TDDFT, the RPA is simply obtained by neglecting the xc-kernel when solving Eq. (8).
In order to calculate correlation energies from Eqs. (7) and (8) it is necessary to generalize the kernel to an arbitrary coupling strength λ. In Ref. [15] it was shown that f λ xc can be obtained from f xc by the rescaling f λ xc (n, q, ω) = λ −1 f xc (n/λ 3 , q/λ, ω/λ 2 ).
In particular it is straightforward to show that any bare exchange kernel satisfies f λ x = λf x . 6

RPA renormalization
While the Dyson equation (8) provides an exact representation of χ for a given kernel, the solution of the equation may exhibit pathological behavior related to electronic instabilities [16,40,41]. The simplest examples are for the low-density homogeneous electron gas and stretched  diatomics, where Colonna et al. [41] demonstrated that the exact-exchange kernel leads to a divergence in the density-density response function computed via Eq. (8). To avoid this problem, an exact refactorization of Eq. (8) was introduced by Bates and Furche in 2013 [42], The series expansion of χ(ω) in powers of χ KS (ω) generates an unscreened perturbation theory that is equivalent to Görling-Levy Perturbation theory [43]. Since the bare Kohn-Sham orbital energy differences appear in the denominator of the non-interacting response function, the unscreened perturbation series diverges for small-gap or metallic systems. [44] This divergence can be eliminated by expanding in powers of the RPA response function, which leads to the following series In addition to eliminating the divergences related to the non-interacting response function, this expansion also eliminates the electronic instabilities resulting from the kernel since inversions are never needed directly involving the xc-kernel. [41,42,45] The decomposition in Eq. (11) also naturally leads to a simple partition of the correlation energy into two pieces The beyond-RPA (bRPA) piece incorporates all of the terms in Eq. (14) beyond the "bare" RPA response function, which can be collected and exactly expressed as [45] By truncating χ λ to a low-order in χ RP A λ , one hopes to faithfully reproduce the infinite-order correlation energy while avoiding the need to invert a function that directly contains f xc . We stress that this approximation scheme can never exceed the accuracy of the infinite-order approach for energy differences and material properties, but it does guarantee the stability of the scheme to compute the correlation energy. [46] Furthermore, this division naturally separates the long-range and short-ranged contributions to the correlation energy, enabling approximations for ∆E bRP A c to be added directly on top of the already robust random phase approximation.
The first-order approximation derived from RPA renormalization, RPAr1, recovers a significant part (∼90%) of the total bRPA correlation energy for a given kernel [45,46] ∆E RP Ar1 This approximation has several key features: it recovers the exact, second-order correlation energy given the exact-exchange kernel, the coupling strength integral can be performed analytically for exchange-like kernels leading to efficient implementations [42,46,47], and it reduces properly to RPA for stretched bonds unlike other second-order schemes such as SOSEX [48]. In fact, within the adiabatic-connection framework, SOSEX can be obtained directly from RPAr1 [42,45,46] through the replacement of one χ RP A with χ KS in Eq. (17) This approximation was shown to be less consistent than RPAr1 due to the reintroduction of the KS response function for molecular energy differences [42] and structural properties of simple solids [46]. To recover the remaining ∼10% of the bRPA correlation energy, corrections beyond RPAr1 to the response function can be systematically added order-by-order until convergence to Eq. (8). Rather than compute these terms exactly, a simple approximation can be introduced to eliminate the coupling-strength integration and utilize information from second-order to estimate third and higher order terms in the RPA renormalized expansion. This approximation method was termed the Higher-Order Terms (HOT) approximation [49] and is obtained through a rescaling of the second-order RPAr correction at λ = 1. The HOT approximation usually reproduces the total correlation energy to within 1-2%, and, consequently, accurately reproduces the performance of a given kernel for chemical or physical properties of molecules and materials.

xc-kernels from the HEG
Although Eq. (9) provides a definition of the kernel f xc , the absence of an exact expression for the exchange-correlation potential v xc requires that an approximate form of f xc must be used in practical calculations. The homogeneous electron gas (HEG) provides a valuable testing ground for approximations of f xc and allows the kernel's limiting behaviour to be studied. The analogue of Eq. (8) for the HEG is where the dependence on the wavevector q has now been made explicit. χ 0 is the textbook Lindhard function [50], which coincides with χ KS when Eq. (8) is applied to the HEG [51]. A quantity commonly found in the HEG literature is the local field factor G(q, ω), which is closely related to f HEG xc as f HEG xc (q, ω) = −v c (q)G(q, ω). Theoretical work on G (and thus f HEG xc ) can be traced at least as far back as Hubbard (see Sec. III C of Ref. [52] for a review), and exact limits have been derived for a number of cases.
First, the long wavelength and static limit (q → 0, ω = 0) actually corresponds to the adiabatic local density approximation (ALDA) commonly employed in TDDFT, k F = (3π 2 n) 1/3 is the Fermi wavevector for the HEG of density n, and E C is the correlation energy per electron. The two terms in Eq. (20) correspond to exchange and correlation contributions. Eq. (19) is intuitive in stating that the ALDA should be exact in describing the HEG response to a uniform, static perturbation [53], and is more formally derived from the compressibility sum rule [52]. Next, the short wavelength and static limit has the form [54,55] f HEG while the long wavelength, high frequency limit has the form [56] f HEG The parameters A, B, C and D depend on the HEG density, which in turn can be written in terms of the Fermi wavevector or Wigner radius r s = (3/4πn) 1/3 . Practically, A, C and D can be obtained from a parameterization of the correlation energy E C for a HEG of density n, wherase B requires additional knowledge of the momentum distributio of the HEG [57]. For intermediate q values it is necessary to turn to diffusion Monte Carlo calculations. The study of Ref. [58] investigated the q-dependence of the static kernel f HEG xc (q, ω = 0) for a range of densities. A key conclusion of that work was that for wavevectors q ≤ 2k F , f HEG xc (q, ω = 0) remains close to its q = 0 value, (=f ALDA xc , Eq. (19)), while for q > 2k F , the kernel can be reasonably well described by the short wavelength limit (Eq. (21)).
In the context of approximations to f HEG xc , it is worth stressing a point discussed in Ref. [51]: there is no particular reason why approximate, frequency-independent kernels should display the same limiting behaviour as the exact, frequency-dependent kernel evaluated at ω = 0. Indeed, having a frequency-independent kernel which is finite at large q (obeying Eq. (21)) will in fact lead to a pair-distribution function which is singular at the origin [59][60][61].
In Fig. 1 we show some approximate forms for f HEG xc , which have been proposed in the literature [17,[62][63][64]. The "rALDA" kernel will be discussed in some detail in the following sections. Briefly describing the other kernels, "CDOP" refers to the frequency-independent kernel proposed by Corradini, Del Sole, Onida and Palummo [62], which has the same limiting behaviour as the exact static kernel (Eqs. (19) and (21)). "CDOPs" refers to the kernel introduced in Ref. [63], which modifies CDOP so that it vanishes at large q. "CPd" refers to the dynamical kernel proposed by Constantin and Pitarke [64], which satisfies the long wavelength static and high frequency limits (Eqs. (19) and (22)). The frequency-independent "CP" kernel corresponds to the CPd kernel at ω = 0. Although the HEG carries the advantage of being a very well-studied system, it is worth remembering that fundamentally it is metallic. The exchange-correlation kernel of a periodic insulator is known to display different limiting behaviour to that of a metal, diverging as 1/q 2 in the q → 0 limit [65,66]. This aspect is especially important in TDDFT calculations of optical spectra including excitonic effects [67][68][69]. This consideration led to the development of the frequency-independent jellium-with-gap model kernel, which has the 1/q 2 divergence [67]. The slightly simpler "JGMs" kernel shown in Fig. 1 is described in Ref. [70]. Here the band gap E g enters parametrically. In the limit E g → ∞ the correlation energy disappears, while E g → 0 the metallic CP kernel is recovered.
Ref. [70] provides a more thorough discussion of all of the kernels shown in Fig. 1, including the expressions used to evaluate them and their forms in real space. In the current work we focus our attention on the renormalized adiabatic kernels rALDA and rAPBE, although a comparison of all the different xc-kernels for the structural parameters of solids will be presented in Sec. 4.2.

The renormalized adiabatic LDA kernel
Ideally one should aim at obtaining a general approximation to f xc that can reproduce various physical quantities such as optical absorption spectra and ground state electronic correlation energies. However, finding good approximations for f xc is highly challenging and it is often necessary to limit the approximation to a given application. As mentioned previously, it is crucial to use a form of f xc that has the correct 1/q 2 behavior in the long wavelength limit in order to capture excitonic effects in absorption spectra. On the other hand, ground state correlation energies involve q-space integrals making it extremely important to obtain a good approximation at large values of q, whereas the long wavelength limit is less important. In the following we will focus on obtaining an approximation that provides accurate ground state correlation energies.
The correlation energy per electron is directly related to the integral of the coupling constant wheren c (q) is the Fourier transform ofn c (r). A parametrization of the exact n c (q) has been provided by Perdew and Wang [71] based on quantum Monte Carlo simulations of the homogeneous electron gas at various densities. Approximations for n c (q) can be obtained from Eqs. (7) and (8) using the Lindhard function for χ 0 (q, ω). The correlation hole in q-space is shown in Fig. 2 calculated with RPA and ALDA for two different densities. Compared to the exact parametrization it is clear that RPA severely overestimates the magnitude of the correlation hole and the RPA will predict a correlation energy that is ∼ 0.5 eV too low per electron for a wide range of densities. The ALDA on the other hand straddles the exact parametrization for a wide range of q-values but decays too slowly at large q compared to the exact results. This is a consequence of the locality of the approximation, which translates into an independence of q. At large q the xc-kernel will thus dominate the Coulomb kernel and fail to reproduce the exact limit (21). Since the total energy involves a q-space integral over all space the slow decay of the correlation hole introduces significant errors and overestimates the correlation energy by ∼ 0.3 eV per electron.
The ALDA x kernel provides a good approximation to the exact one for both low r s = 1 and high r s = 10 densities for q < 2k F , where the correlation hole has a zero point in q-space. However, for q > 2k F the exact correlation hole largely vanishes and we expect to obtain a better approximation for the correlation energy if we simply truncate the q-integration at 2k F when evaluating (23) using the ALDA x approximation. We will refer to this scheme as renormalized ALDA x (rALDA), since the truncation preserves the integral of the correlation hole in real space. The correlation energy per electron evaluated in this scheme is shown in Fig. 3 obtained with RPA, ALDA x , and rALDA. Evidently, the errors in the correlation energy obtained with rALDA are much smaller than both RPA and ALDA x .
For the homogeneous electron gas, the truncation is equivalent to using the Hartree-exchange- Fourier transforming this expression yields Since k F is related to the density, one can attempt to generalize this scheme to inhomogeneous systems. We then take r → |r − r | and k F → (3π 2 n(r, r )) 1/3 , but there is no unique way to define the two-point density n(r, r ). A natural choice is to take [15] n(r, r ) = (n(r) + n(r ))/2, but other choices are possible, which will be discussed in Sec. 3. The simple truncation procedure has thus led to a non-local rALDA kernel that does not contain any free parameters and significantly improves correlation energies for homogeneous systems. We have split the Hartree-exchange-correlation kernel into a renormalized exchangecorrelation part f rALDA xc [n] and a renormalized Coulomb part v r [n]. However, both terms depend on the density and contain exchange-correlation effects. The f rALDA xc part can be regarded as an ALDA x kernel where the delta function has acquired a density dependent broadening, whereas v r is the Coulomb interaction reduced by a density and distance dependent factor that approaches unity for large densities or distances. In fact, at large separation f rALDA Hxc reduces to the pure Coulomb kernel and it is expected to retain the accurate description of long range van der Waals interactions within RPA. For example, in a jellium with r s = 2.0 two points separated by 5Å gives a renormalized Coulomb interaction v r [r s = 2](r) = 0.97v(r) and the magnitude is a factor of 30 times larger than f rALDA xc . Interestingly, both v r and f rALDA xc becomes finite at the origin which implies f rALDA

Hxc
[n](r = 0) = 0. This property is related to the fact that the position weighted correlation hole entering the first integral in Eq. (23) vanishes at the origin [18] and is highly convenient for numerical real-space evaluation of the kernel.
It is often more convenient to separate the Hxc-kernel into the exact Coulomb kernel and an xc-kernel and one is then led to define This expression is typically more useful for applications to periodic systems since and v(r) are long ranged.

Generalized truncation scheme
The truncation scheme defined above is easily generalized to any adiabatic semi-local local kernel. The correlation hole of the homogeneous electron gas calculated with ALDA x has a zero point exactly at 2k F . This is not true in general, but for any adiabatic kernel the correlation hole becomes zero at the point where the Hartree-exchange-correlation kernel vanishes. This leads to the zero-point wavevector where f A xc [n] is the spatial Fourier transform of the adiabatic kernel corresponding to a particular semi-local approximation. Renormalized kernels for any semi-local approximation for the exchange-correlation functional can then be defined by replacing to 2k F by q 0 in Eq. (25) and the kernel is again generalized to inhomogeneous systems by taking r → |r − r | and n → n(r, r ) in addition to a scheme that defines n(r, r ) in terms of n(r). For generalized gradient corrected functionals, q 0 will depend on the gradient of the density as well, which may lead to positive values of f A xc at certain points. At those points we set q 0 [n](r, r ) = 0 in order to maintain a well defined kernel. Below we will only consider rALDA and rAPBE.

Spin
The inclusion of spin degrees of freedom in RPA is almost trivial since the correlation energy involves the sum over all spin components χ σσ , which is obtanined by the simple substitution . This is due to the fact that F Hxc is independent of spin in RPA. In general, however, it is not straightforward to generalize a kernell for spin-paired systems to the spin-polarized case.
In the case of exchange one can resort to the spin dependence of the exchange energy. In particular one has which yields It is possible to enforce this condition on the rALDA kernel as well, but we have found that it makes the correlation energy difficult to converge. The reason is that the off-diagonal (in spin) componenets of the Hxc-kernel involves a bare Coulomb interaction, whereas the diagonal components lack a long-range cancellation between v(r) and v r [n]. This failure is clearly a limitation of the rALDA scheme and an additional approximation is required to maintain the the accuracy of rALDA for spin-polarized systems. To this end we start with the Dyson equation with explicit spin dependence where it was used that χ KS is diagonal in spin. For the spin-paired case one has that which will always hold if Eq. (32) if Eq. (32) is satisfied. To reintroduce a balanced expression for the renormalized kernel in each spin component we relax Eq. Eq. (32) and use where n = n σ + n σ . Eq. (34) is now satisfied, but Eq. (32) is not. This choice is not unique though and another choice is comprised by f rALDA − v, which was used in Ref. [17]. However, Eq. (35) appears to yield better results for atomization energies where spin-polarized isolated atoms are used as a reference.

Hedin's equations and vertex corrections
So far we have discussed the use of xc-kernels in the context of the ACFDT formula for the ground state correlation energy. However, it is possible, and in fact quite effective, to apply the same xc-kernels to describe the effect of vertex corrections in the electron self-energy. In 1965 Lars Hedin introduced a set of coupled equations relating the single-particle Green's function G, the electron self-energy, Σ, to the polarization, P , the screened electron-electron interaction, W , and the 3-point vertex function Γ [72], where we employed the notation (1) = (r 1 , t 1 , σ 1 ) and G H is the Hartree Green's function.
The well known and widely used GW approximation is obtained by iterating Hedin's equations once starting from Σ = 0, i.e. the Hartree approximation. This produces the trivial vertex function Γ = δ(1, 2)δ(1, 3), which corresponds to the time-dependent Hartree approximation for the polarization P , which is the approximation refeered to as RPA in the present review. There are basically two issues with this approach. First of all, it starts from G H , which is known to be a poor approximation. Secondly, it neglects vertex corrections completely. In practice, the latter issue is rarely dealt with because of the complex nature of Γ, while the first is overcome by following a "best G, best W " philosophy [21]. Within the popular G 0 W 0 method, one evaluates the self-energy from a non-interacting G 0 obtained from a DFT calculation while W is obtained within RPA using the polarization P 0 = G 0 G 0 . Today, the G 0 W 0 method remains the stateof-the-art for calculation of QP band structures of inorganic solids [23][24][25] and nano-structures including two-dimensional materials [73][74][75].
Another question related to Hedin's equation is the role of self-consistency. In principle, the five equations should be solved self-consistently. However, while self-consistency improves the description of energy levels in molecules [76,77] and is essential for systems out of equilibrium [78][79][80][81], it does not in general improve the the band structure and spectral functions of solids when vertex corrections are neglected [82,83]. The role of self-consistency will not be further discussed in this review where we instead concentrate on the problems of vertex corrections.
Rather than starting the iterative solution of Hedin's equations with Σ = 0 (which leads to the GW approximation at first iteration), it is of course possible to start with a local approximation, Σ 0 (1, 2) = δ(1, 2)v xc (1). As shown by Del Sole et al. [28] this leads to a self-energy of the form where and is the adiabatic xc-kernel. By inspection it becomes clear that W (1, 2) is the screened effective potential generated at point 2 by a charge at point 1. This potential consists of the bare Coulomb potential plus the induced Hartree and xc-potential. Consequently, it represents the potential felt by an electron. In contrast, the potential felt by a classical test charge is the bare potential screened only by the induced Hartree potential: In Eq. (41) the replacement of W by W thus corresponds to including the vertex in the polarisability, P , but neglecting it in the self-energy. In the following we refer to these two alternative schemes as G 0 W 0 Γ 0 and G 0 W 0 P 0 , respectively. As usual the subscripts indicate that the quantities are evaluated non self-consistently starting from DFT. In addition to the fact that the vertex correction accounts for the change in the xc-potential and therefore should be more accurate, an attractive feature of the G 0 W 0 Γ 0 scheme is that DFT becomes the consistent starting point for non-selfconsistent calculations when the relation (43) is satisfied. This is in stark contrast to G 0 W 0 for which the Hartree approximation is the consistent starting point. In Sec. 4.11 we show that when the rALDA xc-kernel is used to include vertex corrections through Eq. (42) the improved description of short-range correlations lead to a significant upshift of QP energies by 0.3-0.5 eV in agreement with experiments [20]. Since both occupied and unoccupied states are shifted up, the band gap is not affected as much but a small increase is generally observed again in agreement with experiments.

Evaluating non-local kernels for inhomogeneous densities
Kernels like the rALDA which were derived from the HEG (a uniform system) have the form f m xc [n](q, ω) in reciprocal space or f m xc [n](|r − r |, ω) in real space. As mentioned in Sec. 2.4.1, this nonlocality |r − r | leads to a question regarding the treatment of the density argument when calculating the correlation energy of inhomogeneous systems [n = n(r)]. To illustrate this point more clearly, we consider the plane wave representation of the kernel, Here V is the volume of the crystal, G and G are reciprocal lattice vectors and q lies within the first Brillouin zone. In the case that the system under investigation is homogeneous [n(r) = n 0 ] then the kernel is diagonal, On the other hand, if the kernel is fully local (independent of q, e.g. the ALDA) it is natural to use the local density to evaluate the kernel, obtaining where Ω is the unit cell volume. However, for the general case of a inhomogeneous system and non-local kernel there is no unique way of constructing f xc [n](r, r ) from the knowledge of f xc [n](r) in a homogeneous system. The problem is that the density is a one-point function and it is not clear how to treat the dependence of r and r . One important constraint is that the resulting kernel must be symmetric in r and r [84], and in the following we assume the form where S is a functional of the density symmetric in r and r and we have restricted ourselves to frequency-independent kernels.

Density symmetrization
The density symmetrization scheme used in the rALDA/rAPBE calculations in Refs. [17,19,85] employed a two-point average, but more elaborate functionals are possible [86,87]. A kernel satisfying Eq. (47) with a general two-point density is only invariant under simultaneous lattice translation in r and r . Its plane wave representation can then be written in the form and R ij = R i − R j . f (q; r, r ) is periodic in r and r and Eq. (49) must be converged by unit cell sampling, which should typically match the k-point sampling in periodic systems.

Kernel symmetrization
A second approach symmetrizes the kernel itself [63]. Starting from a non-symmetric kernel, and inserting into equation 45 gives It is now possible to obtain asymmetric kernel by taking the average f NS,GG xc (q, ω) and its Hermitian conjugate, which can be seen equivalently as inserting the two-point average 1/2[f NS xc (r, r , ω)+f NS xc (r , r, ω)] into Eq. (45) [63]. Compared to density symmetrization, Eq. (52) has the advantage that the integral is performed over one unit cell only, and that only the density has to be represented on a real space grid, while the kernel is defined by its plane wave representation.

Wavevector symmetrization
The third approach we consider is that employed in Ref. [67], which retains the computational advantages of the kernel symmetrization scheme. Here the wavevector |q + G| entering Eq. (52) is replaced by the symmetrized quantity |q + G||q + G |, such that A kernel constructed using Eq. (54) will automatically satisfy the symmetry requirement of Eq. (47). Furthermore, for the specific case of kernels based on the jellium-with-gap model, the head and wings of the matrix in G and G have the correct 1/q 2 and 1/q divergences, respectively [67].
The main drawback of Eq. (54) is that, compared to a two-point average of the density or the kernel, the procedure of symmetrizing the wavevector is not physically transparent. Of course, two-point schemes also suffer from limitations (e.g. the kernel has no knowledge of the medium between r and r ). The fact that we have to invoke any averaging system at all is an undesirable consequence of using kernels derived from the HEG to describe inhomogeneous systems. In what follows we present calculations using both the density and wavevector symmetrization schemes.

Computational details and convergence
The kernels described in previous sections has been implemented in the DFT code GPAW [88,89], which uses the projector augmented wave (PAW) method. [90] The calculation of correlation energies energies in the framework of the ACFD are performed in four steps. 1) A standard LDA of PBE calculation is carried out in a plane wave basis. 2) The full plane wave Kohn-Sham Hamiltonian is diagonalized to obtain all unoccupied electronic states and eigenvalues.
3) A plane wave cutoff energy is chosen and the Kohn-Sham response function [91] is calculated by putting number of unoccupied bands equal to the number of plane waves defined by the cutoff. 4) The correlation energy is evaluated according to Eqs. (7) and (8). The calculated correlation energies are finally added to non-selfconsistent Hartree-Fock energies evaluated on the same orbitals as the correlation energy. The coupling constant integration is evaluated using 8 Gauss-Legendre points and the frequency integration is performed with 16 Gauss-Legendre points with the highest point situated at 800 eV .
The main convergence parameter for these calculations is thus plane wave cutoff energy (E cut ) used for the response function and kernel. In the case of RPA calculations it has been shown that for sufficiently high cutoff energies the correlation energy scales as [10,85] and it is thus possible to perform accurate extrapolation to the converged results from a few calculations at low cutoff energies. When the ACFD method is used with a kernel the extrapolation (55) is less accurate, but the calculations often converge much faster than RPA such that extrapolation is either not needed at all or only introduce small errors. As an example we show the correlation energy of bulk Na in Fig. 4 calculated with RPA, ALDA and rALDA. It is expected that the correlation energy should resemble that of a HEG with the average valence density of Na due to the delocalized valence electrons in Na [92]. The rALDA calculations are rapidly converged with respect to unit cell sampling (two nearest unit cells are sufficient) and the result are shown in Fig. 4 as a function of plane wave cutoff energy along with the RPA and ALDA results. Similarly to the HEG we find that RPA significantly underestimates the correlation energy while ALDA X overestimates it. We also note the very slow convergence of the ALDA calculation with respect to plane wave cutoff due to the q-independent kernel. In a plane wave representation the non-local kernels considered in the present work takes the form of Eq. (49). While the response function is calculated within the full PAW framework, it is not trivial to obtain the PAW corrections for a non-local functional. However, since the ALDA kernel vanishes for large densities, the non-local kernels considered in the present work tends to be small in the vicinity of the nuclei where it is usually difficult to represent the density accurately. As a consequence the kernels are rather smooth -even at the points where the density is non-analytical -and the kernel can be evaluated from the all-electron density represented on a uniform real-space grid using Eqs. (49) and (50). This is illustrated in Fig. 5 where the correlation energy of an N 2 molecule is shown as a function of grid spacing. The energy difference (contribution to the atomization energy) converges rapidly and is accurate to within 10 meV at 0.17Å. For the calculations in the present work the grid spacing was determined by the plane wave cutoff as h = π/ √ 4E cut and a plane wave cutoff of 600 eV for the initial DFT

Absolute correlation energies
We have already shown that the RPA underestimates the correlation energy of the homogeneous electron gas by 0.6-0.3 eV per electron, whereas the ALDA overestimates the correlation energy by 0.3 eV per electron compared to the RPA. This is significantly improved by the rALDA functional, which gives an error of less than 0.05 eV (See Fig. 3). In Tab. 1, we show that this trend remains true for simple atoms and molecules. For the H atom RPA gives a correlation energy of -13 kcal/mol (-0.56 eV) whereas ALDA gives 6 kcal/mol (0.26 eV). rALDA on the other hand gives -2 kcal/mol, which is a factor of three better than ALDA and a factor 6 better than RPA. A similiar picture emerges from the correlation energy of H 2 and the He atom.  Figure 6: Absolute correlation energy per valence electron calculated with different xckernels [70]. Also shown is the correlation energy for Si obtained from diffusion Monte Carlo (DMC) calculations in Ref. [97].
Although the RPA is free from self-interaction (cancellation of Hartree and exchange for single electron systems), it still contains a large amount of self-correlation, which originates from the fact that the Hartree kernel is not balanced by an exchange kernel in the ACFDT formalism. The self-correlation can be cancelled exactly by including second order screened exchange (SOSEX) [14,94] or an exact exchange kernel [95,96]. However, these approaches are far more computationally demanding than the ACFDT formalism. In contrast the rALDA kernel has a similar computational cost as RPA and reduces the self-correlation error to less than 0.1 eV for a H atom. The remaining error in rALDA is largely due to the choice of the LDA functional as the starting point and choosing the rAPBE kernel instead reduces the error to less than 1 meV.
In Figure 6 we show the correlation energy per valence electron calculated for a number of solids using various xc-kernels, and the RPA [70]. As for the molecules and homogeneous electron gas, the RPA correlation energy is consistently larger (by a few tenths of an eV/electron) than that calculated using the xc-kernels. However, the differences in correlation energy between the kernels themselves are much smaller. For instance, including the correlation contribution of the ALDA in the rALDA kernel ("rALDAc") increases the correlation energy by around just 1% (∼0.01 eV/electron) compared to the standard rALDA. Fig. 6 also shows the result of diffusion Monte Carlo (DMC) calculations of the correlation energy of Si [97], which we can tentatively compare to our own results. The DMC correlation energy lies among the values calculated from the kernels. Ref. [63] also found a correlation energy close to the DMC value using the CDOP kernel and a different averaging scheme. This result is of course reassuring, but we note that care must always be exercised when making such comparisons. First, one must expect the calculated value to have some dependence on the treatment of the core-valence interaction (e.g. PAW, pseudopotentials, all-electron). More generally, the concept of the correlation energy "per valence electron" becomes less well-defined if the calculated correlation energy includes the contribution of semicore states. For instance, comparing Figs. 6 and 4 reveals an apparent contradiction, that the valence energy per electron of Na in Fig. 6 is apparently approximately half its value in Fig. 4. This discrepancy is resolved, however, by noting that in the calculations of Fig. 6[70], the entire 2s 2 2p 6 shell of Na was included as valence in addition to the 3s electron, unlike in Fig. 4. Therefore the correlation energy of Na reported in Fig. 6 represents an average over the free-electron-like 3s electron and the more localized 2s 2 2p 6 shell, which cannot be straightforwardly compared to the free electron gas as in Fig. 4.
It is remarkable that the amount of self-correlation introduced by RPA is similar for widely different systems and it indicates that there will be a large energy cancellation when considering energy differences. As we will see below this is true to some extent, but the cancellation is far from perfect and RPA gives rise to systematic errors in cohesive energies of solids and atomization energies of molecules.

Effect of kernel averaging scheme
As discussed in Sec. 3.1, there is a choice in how one constructs the XC kernel for an inhomogeneous system. For instance, the correlation energies displayed in Table 1 were calculated using a two-point average of the density. If we repeat the rALDA calculations using the symmetrized wavevector scheme (equation 54) we obtain correlation energies of 6, -24 and -23 kcal/mol for H, H 2 and He, i.e. a difference of +8, +4 and +4 kcal/mol compared to the values of Table 1. The small differences in H and H 2 cancel when calculating the atomization energy [70]. We find it encouraging that the symmetrized wavevector approach agrees very well with the more intuitive two-point density average when calculating the atomization energy.

Structural parameters
Having established that the renormalized kernels yield greatly improved absolute correlation energies, we now consider physical observables, starting with lattice constants and bulk moduli of a test set of 10 crystalline solids. In these calculations, the Kohn-Sham states and energies obtained self-consistently within the LDA were used to calculate the noninteracting response function and exact exchange contribution to the total energy. A number of different XC kernels (including the rALDA) were used to calculate the response function and correlation energy through equations 7 and 8. The wavevector symmetrization scheme (equation 54) was used to construct the kernel. The lattice constants and bulk moduli were extracted by calculating the total energy as a function of lattice spacing and fitting the results to a Birch-Murnaghan equation of state. More computational details are given in Ref. [70]. Figure 7 shows the percentage deviation between the calculated structural parameters and the experimental data listed in Ref. [6]. Calculations performed within ground-state DFT within the LDA or GGA are also shown for comparison, which show the well-known tendency for the LDA/GGA to over/underbind, respectively. Using exact exchange and the ACFD correlation energy (with any kernel, or the RPA) systematically improves the agreement with experiment, going from a mean absolute error of 1.3%/7% for PBE to ≤0.7%/4% for lattice constants/bulk moduli, respectively.  Figure 7: Lattice constants (left) and bulk moduli (right) calculated using different approximations for the XC kernel for a test set of 10 crystalline solids, compared to the experimental values listed in Ref. [6]. The experimental lattice constants were corrected for zero-point motion [6]. Note the CP and JGMs kernels coincide for metallic systems [70].
The difference between the various kernels, and even the RPA, is rather small. In particular, the rALDA, rALDAc, CDOPs and CP kernels yield very similar results. The very close agreement between rALDA and rALDAc supports the use of the simpler rALDA kernel, which uses only the exchange part of the ALDA (Eq. (25)). One attractive property of the rALDA is that the calculated values display the fastest convergence with respect to the number of plane waves used to construct the response function, allowing a saving in computational time.
In terms of the other kernels, the strongest outlier is the JGMs (jellium-with-gap) kernel, particularly for the ionic solid LiCl. The agreement with experiment for the JGMs lattice constants can be improved even further by replacing the experimental optical gap E g that appears in the kernel definition with an effective gap inspired by excitonic calculations involving a long-range (LRC) attractive kernel [68,70]. One can also see that the CDOP kernel produces respectable structural parameters, despite it actually having a divergent pair-distribution function, while the variation between the CP and CPd kernels illustrate the potential importance of dynamical effects.
However, the overall differences between all of the kernels are rather small, and based on these calculations it is hard to argue that there are particular benefits in going beyond a simple, static kernel which tends to a density-dependent constant at small q and decays as 1/q 2 at large q. The rALDA satisfies these properties and also carries the particular advantage of scaling simply with the coupling constant λ. Finally, a significant strength of the rALDA is that, unlike the other kernels derived from the HEG, it has a spin-dependent generalization (Sec. 2.5), which is essential when calculating molecular atomization energies.

Atomization energies of molecules
In Fig. 8 we compare the performance of LDA, PBE, RPA@LDA, RPA@PBE, rALDA and rAPBE for the atomization energies of 14 small molecules using the experimental atomic positions. All numbers are shown relative to experimental values corrected for non-adiabatic effects [98]. RPA is seen to systematically underestimate the binding energies with RPA@LDA being slightly worse than RPA@PBE. In contrast LDA and PBE overestimate binding energies and the performance of PBE is similar to RPA whereas LDA is much worse. Both rALDA and rAPBE provides a significant improvement over RPA. rAPBE performs slightly better than rALDA, which is most likely due to the poor description of the ground state within LDA compared to PBE. In Fig. 9 we show the mean absolute percentage error (MAPE) of RPA, rALDA and rAPBE compared with that obtained with PBE0 as well as SOSEX and rP2T, which constitute two other beyond RPA methods [99,100]. rALDA and rAPBE are a factor of three more accurate than RPA@LDA and RPA@PBE respectively. Moreover, the rAPBE MAPE is less than 1.5 % and outperform both PBE0 and r2PT on this small test set.

Cohesive energies of solids
While hybrid functionals may provide a computationally cheap way of obtaining accurate ground state energies for atoms and molecules, they typically fail dramatically for solids. Moreover, quantum chemistry methods are prohibitively demanding for solid state systems and DFT and DFT-based methods currently seem to be the only possible choice when dealing with solids. In Ref. [6], it was shown that RPA performs somewhat worse than PBE for the cohesive energies of solids although it does provide significantly better results than LDA. This is in contrast to the case of molecules where RPA yields slightly better results than PBE. The reason is that the accuracy of RPA for atomization energies crucially depends on error cancellation of the ubiquitous self-correlation in RPA. The cancellation of errors is likely to work better when comparing similar systems, but for the cohesive energy of solids one has to consider ground state energies of atoms with ground state energies of solids, so the error-cancellation can become more inaccurate. Since the rALDA and rAPBE functionals to a large extent eliminate the selfcorrelation error of RPA, it is expected that these approaches should perform significantly better than RPA. In Fig. 10   RPA. The mean absolute percentage error is shown in Fig. 11 and the rAPBE deviation from experiment is less than 2 %, whereas PBE, RPA@LDA and RPA@PBE give errors of 4 %, 9 % and 7 % respectively. PBE0 gives a mean error of 7.5 %, which is four times worse than rAPBE.

Formation energies of metal oxides
In the previous two sections we considered the problems of calculating the atomization energies of molecules and solids. This problem gauges the ability of a method to describe the absolute energy cost of breaking a chemical bond. In most practical situations, however, it is often more relevant to consider the material's formation energy, i.e. its energy relative to the standard states of its constituent elements rather than the isolated atoms. The calculation of formation energies thus gauges the ability of a method to describe the energy of one type of chemical bond relative to another. Predicting the heat of formation of metal oxides has proven to be particularly challenging for a wide range of commonly applied xc-functionals. The RPA has previously been shown to significantly improve the accuracy of calculated formation energies of group I and II metal oxides as compared to semi-local functionals [102]. In the following we briefly assess the performance of the rAPBE method and compare it PBE, RPA and the BEEF-vdW functional, which contains non-local correlation to account for van der Waals interactions [103]. The formation energy per oxygen atoms was obtained from the computed total energies as ] are the total energies of the oxide, the bulk metal and the O 2 molecule in the gas phase respectively. Zero-point energy contributions were not included in the present study as previous work has shown that they affect the formation energies of oxides by less than 0.01 eV [102]. The formation energies computed with PBE, BEEF-vdW, EXX, RPA, and rAPBE are summarized in Fig. 12. For the latter three methods, single-particle wave functions and energies were obtained from a self-consistent PBE calculation. All structures were optimized with PBE. The BEEF-vdW was included here to compare the performance of RPA and rAPBE methods to a semi-empirical method that explicitly includes dispersive interactions. The mean signed error (MSE), mean absolute error (MAE) and mean absolute percentage error (MAPE) with respect to experiment are shown Table 2. Comparing the MSE and MAE shows that formation energies from PBE, BEEF-vdW, EXX and RPA have a clear systematic tendency to overstimate formation energies; that is oxides are predicted less stable than found in experiments (see Ref. [102] for references for the experimental data). In contrast, while rAPBE shows the same tendency   of destabilizing the metal oxide, it is less pronounced, and CaO 2 , KO 2 , CsO 2 and RuO 2 are in fact predicted to be more stable than experiment. The rAPBE method is in better agreement with experiments than all the other methods with a MAE of only 0.21 eV compared to 0.38 eV for RPA and 0.40 for BEEF-vdW. We can partly attribute the failure of the RPA to the lack of error cancellation between the correlation energy of the oxide and the bulk metal and oxygen molecule, which are all separately underestimated by the RPA (see Ref. [104]). The errors of the DFT xc-functionals and the RPA are to some extent systematic and can be ascribed to a bad description of the O 2 molecule. In fact, treating the energy of the O 2 reference as a fitting parameter, the MAE for all the methods become comparable and lie in the range 0.15-0.2 eV/O [104].

Surface-and adsorption energies
For applications of DFT to problems in surface science, in particular heterogeneous catalysis and electrocatalysis, the ability to predict stability and reactivity of metal surfaces is of crucial importance. It has been established that the RPA yields very good results for surface energies and chemisorption energies of atoms and small molecules on transition metal metal surfaces and greatly improves the accuracy of the xc-functionals [105][106][107][108][109]. As shown below, the good performance of RPA for surface-and adsorption energies is preserved and probably even improved by the renormalized kernels. The difference between RPA and rALDA for surface reaction energies is on the order 5-10%, which is comparable to the difference found for atomization energies of solids and molecules. This suggests that the better description of short range correlations by the kernel, which was found to improve atomization energies in Secs. 4.3 and 4.4, carries over to metal-molecule bonding. However, due to the significantly smaller magnitude of such bond energies compared to covalent bonds in solids/molecules and the lack of experimental data with sub 100 meV accuracy, it is not possible to confirm this hypothesis directly. Figure 13 shows the adsorption energy of CO on Pt(111) for a coverage of 1/4 plotted against the surface energy. Results obtained with RPA, rAPBE, and a range of xc-functionals are shown together with the experimental result states in Ref. [105]. It is evident that the RPA and rAPBE methods are able to break the incorrect correlation between adsorption energy and surface energy exhibited by the xc-functionals [105] and thereby yield an accurate description of both quantities simultaneously. Similar conclusions have recently been demonstrated for other adsorbates and surfaces [110]. Table 3 reports reaction energies in the full coverage limit for the set of benchmark surface reactions introduced in Ref. [111] over the early 3d transition metal surfaces. The absolute and relative deviation from the rALDA result is shown in the last two rows. It is interesting to note that the xc-functional yielding the best agreement with the rALDA (and RPA) is the RPBE, which is a revised version of the PBE fitted to experimental data on surface reactions like the ones considered here [112]. The rALDA and RPA results are in good agreement, with the largest deviation occurring for OH adsorption (difference of 0.11 eV). This agrees well with results of Fig. 13 for CO adsorption on Pt (111) and for graphene adsorbed on Ni(111) (see Fig. 16). Although the absolute deviation between rALDA and RPA for the surface reactions is small (MAE of 40 meV), the relative deviation is in fact similar to that obtained for the atomization energies of molecules and solids. However, it is interesting to note that RPA shows a small but systematic tendency to overbind the adsorbates relative to rALDA. This is opposite to the trend observed for the atomization energies of molecules and solids where RPA underestimates the bond strength by around 0.3-0.5 eV/atom relative to rALDA (see Figs. 8 and 10). For metals, the cohesive energy is also underestimated by RPA, but by somewhat smaller degree (MSE of -0.15 eV relative to rALDA). The average deviation between RPA and rALDA for the bond energy per atom is shown in Table 4 for the set of materials in Figs. 8 and 10 grouped according to material type.
The fact that molecule-metal bonding shows an opposite trend compared to atomization energies can be explained as follows: for reactions involving the breaking of covalent bonds in the initial adsorbate molecule, e.g., for reactions of the type 1/2A2 -A/metal (reactions 1, 2, 6-8 in Table 3), RPA will overestimate the reaction energy because the A-A bond strength is underestimated more than the A-metal bond (the latter being weaker than the former). For pure adsorption reactions of the form A2 -A2/metal (reactions R3-R5), RPA will overestimate the adsorption energy because the reduction of the internal A-A bond upon adsorption is underestimated more than the A2/metal bond. In both cases the reason for the (slight) overestimation of the reaction energy can thus be traced to a larger underestimation by the RPA in describing pure covalent bonds compared to bonds with partial metallic character.

Static correlation
High highly attractive feature of the RPA is the correct description of bond dissociation in molecular dimers, which cannot be captured by ()restricted semi-local functionals. The bond dissociation curves in RPA are, however, only accurate of the reference energy of the isolated atoms are subtracted. In the case of H 2 , for example, the dissociation curve is 1.1 eV below the exact result due to the self-correlation energy of the H atom. As we have seen, the rALDA  Table 3: Adsorption energies (in eV) for a few reactions at full coverage calculated with rALDA, RPA and seven different DFT xc functionals. Data taken from Ref. [110].   large eliminates the self-correlation and one can obtain accurate dissociation curves of molecular dimers become without subtracting reference energies for the isolated atoms. This is shown in Fig. 14, where we also compare with LDA, PBE and Hartree-Fock, which all fail to yield the correct disociation energy. We note, however, that rALDA as well as rALDA exhibits a spurious maximum at intermediate bond lengths, which is not present in the exact result. The exact monotonous increase in energy signals a failure of both approximations. We note that is has previously been shown thatr the correct dissociation curve can be obtained from the ACDFT formalism if the interacting response function χ is evaluation from the Bethe-Salpater equation [113].

MAE
Finally it is worth mentioning that both RPA and rALDA fails dramatically for the dissociation of the H + 2 dimer. In contrast, the SOSEX method yields exact result for this problem but cannot dissociate H 2 correctly.

Dispersive interactions
One of the main qualities of the RPA is its ability to account for van der Waals interactions in weakly interacting systems. Specifically, the RPA has been shown to provide an excellent description of the equilibrium geometry in hBN [114] and graphite [9] as well as graphite adsorbed on metal surfaces [11,12,85]. Conserving the accurate description of dispersive interactions is thus a major success criterion for any beyond-RPA method for ground state energies.

Bilayer graphene
In Fig. 15 we show the binding energy as a function of interlayer distance for a bilayer of graphene calculated with LDA, PBE, the van der Waals functional of Dion et al. [115] (vdW) along with the results of RPA and rALDA. The RPA and rALDA is seen to yield nearly identical equilibrium distances of 3.4Å, but rALDA gives a slightly smaller binding energy (22 meV per atom) compared to RPA (25 meV per atom). The semi-local functionals are not expected to capture dispersive interactions and PBE predicts a very shallow minimum at 4.4Å. However, LDA yields an equilibrium distance of 3.3Å, which is close to the RPA/rALDA value, but only half the binding energy. It is remarkable that LDA provides a seemingly qualitative correct description of the binding energy in this system. Similar results for LDA has previously been obtained for graphite [9] and graphene on metal surfaces [85], but the results are likely to be fortuitous, and the fact that LDA and various GGAs produce qualitatively different results renders the result dubious. The vdW functional predicts an equilibrium distance of 3.7Åand a binding energy of 22 meV per atom. In contrast to LDA, this functional binds the two layers for the right reasons, but the binding distance is quantitatively wrong. It is also noteworthy that the tails of the potential energy surfaces of RPA and rALDA coincide as expected, but deviate from the prediction of the vdW functional. It should be noted that there is no experimental values for the binding energy and binding distance of bilayer graphene, but it is expected that the binding distance should be close to the value of 3.34Å of graphite.

Graphene on metals
The description of graphene adsorbed on metal surfaces has proven a highly challenging task for first principles methods due to equal contributions of weak covalent and van der Waals bonding in these systems [11,12]. In Fig. 16 we compare the binding energy of graphene on a Ni (111) surface calculated with LDA, PBE, RPBE, vdW, RPA and rAPBE [19]. The RPA shows two distinct minima at 2.2Å (chemisorbed) and 3.3Å (physisorbed), which originates from weak covalent interactions and dispersive forces respectively. The chemisorption minimum is a few meV lower than the physisorption minimum and agrees well with the experimentally determined value of 2.1Å for this system. The rAPBE functional closely mimics the RPA result, but gives a slightly larger energy difference between the two local minima. The vdW functional completely misses the weak covalent interactions responsible for the chemisorption minimum and only yields a physisorbed state at 3.8Å.

C 6 coefficients
Long-range interactions are typically dominated by Coulomb interactions and accuracy of dispersive interactions are thus determined by the response functions of the individual systems. For well-separated atoms the binding energy asymptotically becomes E B (r) = C 6 /r 6 , where the C 6 coeffients are only depend on the polarizabilities of the isolated atoms. In particular the C ij 6 coefficient relating atoms i and j can be calculated from the Casimir-Polder formula where is the polarizability of atom i and χ i is the interacting response function, which can be calculated form the Dyson equation with a given approximation for the xc-kernel. The C 6 coefficients for eight different atoms (i = j) is displayed in Tab. 5 calculated with LDA (χ = χ KS ), RPA and rALDA. We observe that RPA performs significantly better than LDA, but rALDA is a factor of three better than RPA on average. The performance is, however, very dependent on the type of atom and rALDA performs better for the noble gas atoms, except for He, which is more accurately described in RPA. For Li and Na RPA fails completely, whereas LDA provides rather accurate predictions (better than rALDA for Li and slightly worse tahn rALDA for Na).

Structural phase transitions
Bulk solids usually exist in various polymorphic forms. Under changes in pressure or temperature, one structure may transition into another. The structural phase transition of solids has large theoretical and practical importance. With the external influence, the space-group symmetry and associated internal structural parameters change from one crystal structure to another. Temperature or pressure induced structural phase transitions can also change the electronic structures of the corresponding materials, such as from insulator to metal and vice versa, resulting in changes of the band-gap or conductance [116]. Structural phase transitions also sometimes lead to different magnetic states [117]. As a consequence, structural phase transitions offer an opportunity to tune a material toward particular applications in electronics, optics and other relevant fields [118][119][120]. Structural phase transitions have largely remained a challenge for electronic structure methods. When computing phase transitions, a robust theory must overcome the dilemma of simulta-  Table 5: C 6 coefficients between identical atoms (i = j in Eq. (57)) calculated with LDA, RPA, and rALDA. All values are in atomic units. We also show the mean absolute relative error (MARE).
neously predicting equilibrium structural properties and the phase transition parameters. Since experiments often fail to precisely measure the coexistence temperatures or pressures of two different structural phases of a solid, there is a high demand for a robust theoretical method. The failure of semilocal density functional approximations for structural phase transitions was earlier attributed to the underestimation of the band-gap with these methods [116,[121][122][123].
This assumption was later questioned by the results from the HSE (Heyd-Scuseria-Ernzerhof) approximation [124,125]. HSE is nonlocal in the exchange and predicts more realistic fundamental gaps. HSE is better than semilocal functionals for the transition pressures of Si and SiO 2 , but seriously overestimates the transition pressure in Zr.
In combination with the rAPBE and rALDA kernels, the RPAr approximation was investigated for the structural phase transitions of a small but representative group of materials [126]. The examples were chosen to incorporate several changes in band structure, including from semiconductor to semiconductor, semiconductor to metal, and metal to metal transitions. The assessment includes the phase transitions of the diamond phase of Si to the metallic beta-tin form, the zinc-blende (ZB) to rocksalt transition of SiC, the ZB to Cmcm phase transition of GaAs, the quartz to stishovite transition of SiO 2 , the transition from fcc to hcp structure of Pb, and finally the phase transition of the hexagonal to cubic structures of BN.
The transition pressure can be found as the negative slope of the common tangent line between the two phases. At the transition pressure, the difference in Gibbs free energies for the two phases should be equal to zero. At a finite, but constant, temperature, the pressure can be calculated as the negative derivative of the free energy with respect to the volume At zero temperature the equivalent condition is that the enthalpy difference of the two phases is zero, and the pressure can be found directly as the derivative of the electronic energy.
To include the thermal effects, the equilibrium parameters were obtained from fitting the  Table 6: Finite temperature transition pressures, in GPa, predicted at 300 K using different levels of density functional approximations [126]. The rAPBE kernel was used in combination with RPAr to obtain the bRPA results. RPAr1 is sufficient at capturing the needed bRPA correlation for predicting the phase transitions, and the HOT correctio introduces a small shift.
third-order Birch-Murnaghan equation of state (EOS) including the thermal corrections from the vibrational degrees of freedom. The results of the fitting are F 0T , the minimum of the Helmholtz energy at temperature T and (V 0T , B 0T , B 0T ), the equilibrium parameters such as the equilibrium volume, bulk modulus and the derivative of the bulk modulus at that same temperature. The transition pressures were then obtained using the isothermal EOS fitting parameters obtained from Eq. (9) in Ref. [126]. Compared to zero Kelvin, the addition of thermal corrections introduces a rigid shift in the performance of all the methods for most of the materials in the assessment, and the agreement typically improves with experiment when the thermal corrections are included. The results of different functionals for predicting the phase transitions generally follows the same trend as the structural parameters, Figure. (1) in Ref. [126]. In general, all beyond-RPA methods deliver an overestimation of the equilibrium and transition volumes by about the same amount as the local density approximation (LDA) underestimates them. Beyond-RPA approximations yielded an overall improvement compared to the bare RPA and semilocal functionals for the transition pressures, Table 6, although the behavior is less systematic than for RPA.
For Si and Ge the RPAr1 approximation delivers a reasonable result, however, RPA shows closer agreement with the experimental transition pressure. For SiC, neither RPA nor beyond-RPA show agreement with the experimental transition pressure, indicating that there is a significant difference between equilibrium and non-equilibrium transition pressures [128]. For SiO 2 , and GaAs the transition pressure predicted with an xc kernel are more accurate compared to experiment than the bare RPA. Adding bRPA correlation from rAPBE at any level of RPAr reduces the transition pressure of RPA and comes quite close to the experiment. Ref. [129] attributes this deviation of RPA to its poor performance for some molecular-like solids where there is less cancellation of error between dissimilar phases. It is remarkable that RPA fails to predict the correct phase ordering for BN, while the rAPBE kernel in conjunction with RPAr brings the transition pressure close to the experimental value.  [126]. The negative slope of the common tangent line corresponds to the transition pressure. The kernel corrections for SiO 2 increase the equilibrium energy difference between phases and correct the large underestimation of the transition pressure by RPA as a result. The kernel-corrected curves have been rigidly shifted up in energy by 0.05 eV compared to RPA for visual clarity. more prominent role in the phase transitions of carbon and BN, and indicate further investigations are necessary. RPA yields too low a transition pressure for the phase transition of carbon without thermal corrections, and reverses the phase stability of BN without a finite-temperature correction. All beyond-RPA approximations largely overestimate the transition pressure in C with the inclusion of thermal corrections, while the prediction for BN is significantly improved compard to RPA. The inclusion of the higher-order terms with the rAPBE kernel leads to a better accuracy against the experimental transition pressure in BN. The unexpected inaccuracies in materials such as Pb, C and BN, can be explained by near degeneracies.
To illustrate the benefits of including the exchange-correlation kernel, Fig. 17 shows the EOS data and common tangent results for SiO 2 predicted with RPA and the HOT approximation in combination with rAPBE. In this case, RPA underestimates the energy difference between the phases resulting in an underestimation of the transition pressure. The root cause of this can be understood if the transition pressure is thought of as ∆E/∆V evaluated at the transition volumes for each phase. The energy underestimation is more severe than the transition volume errors, and so RPA underestimates the pressure (as with most semilocal functionals). With the addition of the xc-kernel, the energy difference between phases is increased by an appropriate ammount to bring the rAPBE transition pressure within the experimental range.

Cesium halide stability
The difficulty in predicting the energy difference between similar phases of a material is a more general problem than phase transitions alone. The stability of different phases in alkali-halides is also a strong probe of various electronic structure methods and the correlation effects they incorporate. Among the alkali halides, those formed from cesium show an interesting behavior. CsF is experimentally stable in the B1 (rocksalt or NaCl) structure, while the Cl, Br, and I materials exist experimentally in the B2 (CsCl) structure. For all of the other alkali halides, the stable structure is B1. The stability of the B2 phase for certain cesium halides is a direct consequence of weak van der Waals bonding [130][131][132][133][134][135]. All ACFD-based methods naturally account for long-range van der Waals forces [136], but an accurate treatment of their structure requires both short and long-range interactions to be accounted for. Semilocal density functional approximations miss the long-range van der Waals forces. Ref [134] indicated the relevance of long-range interactions through the PBE+D2 method compared to the bare PBE-GGA for the phase ordering of these cesium halides. Since D2 is a semi-empirical method [137], the energy differences reported in Ref. [134] are not necessarily usable as a benchmark.
To go beyond the semilocal level, the energy differences between the B1 and B2 phases were explored using ACFD-based methods. The performance of the previously discussed ACFDTbased approximations was assessed in detail in Ref [138] for the structural parameters and cohesive energies of cesium-halides. The rALDA kernel was used primarily for this study. Compared to semilocal density functionals, RPA yields superior structural parameters for all of the stable cesium-halides. Beyond-RPA approximations in combination with the rALDA kernel are in general even more accurate for predicting the lattice constants and bulk moduli. The predicted cohesive energies for all of these ACFDT methods are also more accurate compared to the PBE-GGA due to the incorporation of van der Waals interactions. RPA predicts the proper phase stability of all Cs halides whereas PBE only predicts the correct order for the fluoride, Fig. 18.
By definition the cohesive energy includes the energy of the bulk and the constituent atoms. Describing accurately the energy of the bulk and free atoms simultaneously is a challenge for most electronic structure methods, and a method biased towards one paradigm will lead to inconsistent predictions of the cohesive energy. Depending on the kernel and level of RPAr approximation, the resulting methods yield different descriptions for the bulk and free atoms. Overall the cohesive energies of the stable Cs halides are more accurate with RPA for the F and Cl compounds than for the heavier halides. For Br and I, the rALDA kernel within RPAr improves the predicted cohesive energies compared to RPA. A possible explanation is that RPA is applied in a non-self-consistent manner with PBE input orbitals, which is not necessarily the best starting point for the more ionic X-F and Cl bonds. The addition of a kernel tends to improve the short-range correlation of the ACFDT and therefore compensates the error of the PBE orbitals that is more prominent in the cohesive energies for RPA for the larger Cs halides.  [134]. Positive ∆E coh corresponds to the B1 phase being preferred as the ground state,whereas negative values indicate the preferred stability of the B2 phase. PBE predicts all ground state cesium halides to be in the B1 phase whereas all other methods favor the B2 structure except in CsF. S+D3 [139] and S+rVV10 [140] correspond to the SCAN [141] semilocal results plus the dispersion method specified.
In order to correctly predict the splitting between the B1 and B2 phases a delicate balance of short and long-range correlation is required, Fig. 19. The large overestimation of the stability difference between these two phases with the PBE+D2 method indicates the incompleteness of this approximation. RPA in contrast correctly incorporates long-range correlation, but is incomplete for the short-range part and so underestimates the energy difference. The energy difference is increased by any RPAr approximation using the rALDA kernel. rAPBE could be naturally expected to be the ideal kernel for evaluating the bRPA corrections on top of PBE orbitals, but interestingly this kernel struggles to predict the correct phase ordering in these materials. The other kernels tested, including rALDA, CP07 [142] and CDOP [143], all predict a consistent phase stability ordering, regardless of the RPAr approximation, further indicating the error lies in the rAPBE kernel itself [138]. This seems to be an isolated case, however, as the rest of the results in this review clearly demonstrate the utility and success of the rADFT kernels.

Vertex corrected quasiparticle energies
This section presents some results for quasiparticle energies obtained using the G 0 W 0 , GW 0 , G 0 W 0 P 0 , and G 0 W 0 Γ 0 self-energy methods outlined in Sec. 2.6. For the latter two, the rALDA x xc-kernel was used and all calculations used an LDA starting point. The GW 0 refers to eigenvalue self-consistency in the Green function. All structures were relaxed using the PBE xc-functional and QP calculations were based on norm-conserving PAW potentials and spin-orbit coupling was included for the band structures. Further details on the calculations can be found in Ref.
HOT ACSOSEX Figure 19: Bar graph representing the difference in cohesive energies, as in Fig. 18 [20]. Table 7 reports the band gaps obtained with the different self-energy methods together with their mean absolute error (MAE) and mean signed error (MSE) relative to the experimental reference values for eight bulk crystals and three transition metal dichalcogenides in monolayer form. As expected G 0 W 0 @LDA underestimates the experimental band gaps. In agreement with previous findings, iterating to self-consistency in the Green's function (the GW 0 method) improves the situation somewhat, but leads to a small systematic overestimation of the gaps. This overestimation becomes even larger in fully self-consistent GW (not shown) where the MAE and MSE increase to 0.7 eV (see Ref. [20]). Including the rALDA vertex correction in a nonselfconsistent G 0 W 0 Γ 0 calculation, reduces the systematic underestimation of the gap somewhat (the MSE is reduced from -0.19 eV to -0.12 eV). Summarizing, we find that G 0 W 0 Γ 0 and GW 0 perform the best for the band gaps closely followed by G 0 W 0 . Including the vertex only in the polarizability (G 0 W 0 P 0 ) closes the gap because of the increased screening (as previously reported in Refs. [144,145]) and results in a significant underestimation of the gaps. Fig. 20 show the absolute band positions for the valence band maximum (VBM) and conduction band minimum (CBM) relative to vacuum. For the bulk materials the band positions were determined by aligning the Hartree potential of a bulk calculation with the potential in the center of a slab. The slab thickness was 10-24 atomic layers depending on the material and the surface orientation and reconstruction was taken from available experimental studies. The inclusion of the vertex has a striking effect of blue shifting the band edges by around 0.5 eV. Remarkably, this upshift yields a better overall agreement with experimental values (dashed black lines). Interestingly, this shift of band energies is not observed when the vertex is only included only in the polarizability. In addition, no systematic shift of the band edges is observed for the self-consistent GW approximations (without vertex corrections). We are thus led to conclude that the blue shift of band energies originates from the inclusion of vertex corrections in the self-energy.  The physical origin of these effects can be traced to the improved description of short range correlations provided by the rALDA vertex. Indeed, the induced potential from Eq. (44) is the Hartree potential generated by the induced density (the screening cloud of the quasiparticle). This Hartree potential is too deep and results in too deep-lying QP energies. This is the same reason underlying the systematic underestimation of the correlation energy by RPA. Adding the induced xc-potential by Eq. (42) reduces the size of the screening potential and thus shifts the QP energies up. As discussed in Ref. [20], the band gap size is governed by long range correlations, which are well described by the RPA, while the absolute band energies also depend on the short range correlations whose proper description require the vertex function.
The observed upshift in QP energies by around 0.5 eV due to the vertex correction can also be understood by noting that the ionization potential and electron affinity, i.e. the VBM and CBM relative to vacuum, obtained from G 0 W 0 (G 0 W 0 Γ 0 ) can be related to total energy differences between N and N ±1 ground states evaluated from the ACFDT formula employing the RPA (rALDA) with a "frozen orbitals" assumption, i.e. a generalized Koopman's theorem [146]. Indeed, as we have shown in Sec. 4.1 the RPA underestimates the absolute correlation energy by around 0.3-0.6 eV/electron and this error is largely repaired by rALDA due to the improved Energy relative to vacuum (eV) description of the short range correlations (see Fig. 1). The incorrect description of absolute correlation energies by RPA largely cancels for energy differences when the states in question contain the same number of electrons. However, for QP energies where the initial and final states differ by ±1 electron, such errors are directly revealed. Finally, it should be mentioned that the suppression of the large q components in the selfenergy resulting from the rALDA kernel not only improves the description of local correlations but also leads to faster convergence with respect to the number of plane waves and unoccupied states, as compared to standard GW calculations [20]. The situation is very similar to that reported in Fig. 4 for the ground state correlation energy in rALDA versus RPA and ALDA.

Conclusions and outlook
We have reviewed the theory of exchange-correlation (xc-)kernels derived from the homogeneous electron gas (HEG), and illustrated how they can be used to obtain ground state correlation energies and quasiparticle band structures beyond the RPA and GW methods, respectively. While several xc-kernels have been introduced, we have mainly focused on the renormalized adiabatic kernels rALDA and rAPBE, which are obtained from the (semi)local LDA and PBE xc-functionals by a simple renormalization procedure. The renormalization procedure consists of a truncation of the large q-components of Fourier transformed xc-kernel, which renders it non-local in real space -an essential property to avoid a divergent on-top correlation hole.
The main observation is that the xc-kernels greatly improve the description of the shortrange correlations relative to RPA (which correspond to setting f xc = 0)). In particular, the coupling constant averaged correlation hole, obtained from the density response function via the fluctuation dissipation theorem, becomes almost exact for the HEG and reduces the error on the correlation energy from about 0.5 eV/electron in the RPA to below 0.05 eV/electron. This improvement was shown to carry over from the extreme limit of delocalized electrons in the HEG to limit of localized states of simple atoms and molecules. For example, the spurious self correlation energy of the hydrogen atoms is reduced from 0.56 eV (RPA) to 0.04 eV (rALDA). The much better reproduction of absolute correlation energies reduces the reliance on error cancellation effects when computing energy differences. For RPA, error cancellation is very significant on the order of 0.5 eV/electron. For this reason RPA performs well for isoelectronic problems, problems where the electronic structure of the initial and final states are similar, e.g. for the calculation of structural parameters or the breaking/formation of weak (dispersive) bonds. When the xc-kernels are invoked the need for error cancellation is reduced, and as a consequence covalent bond energies, for which short-range correlations play an important role, are much more accurately described. This was explicitly demonstrated for atomization energies of molecules, cohesive energies of solids, formation energies of metal oxides, surface energies, and chemisorption energies of atoms and molecules on metal surfaces. In all these cases, the xckernels cure the systematic underbinding by the RPA and significantly improve the agreement with experimental data. Importantly, because the xc-kernels are of short-range nature they do not affect the shape of the correlation hole at longer distances, and consequently the excellent performance of the RPA for van der Waals interactions is preserved or even improved. For example the mean absolute error on the C 6 coefficients of noble gas and alkali elements are reduced from 0.3 eV (RPA) to 0.1 eV (rALDA). Similar conclusions apply to bonding with intermediate correlation ranges such as the mixed covalent-dispersive forces governing organicmetal interfaces as exemplified here by the prototypical graphene/metal interfaces. Furthermore, the positive effects that the xc-kernels bring can be adequately captured using RPA renormalized perturbation theory to first order, demonstrating that the impact of higher-order correlation effects are less important for correcting the shortcomings of RPA for structural and energetic properties.
In the context of quasiparticle band structure calculations within the formal framework of Hedin's equations, we have shown that the renormalized xc-kernels can be used to include vertex corrections in the electron self-energy. When the four-point kernel of many-body perturbation theory is approximated by the two-point kernel of time-dependent DFT, the effect of the vertex corrections attains a simple and physically transparent form. Namely, the effect of the vertex is to include the xc-potential into the dynamically screened Coulomb potential. We have shown that this is crucial in order to obtain a correct description of absolute band energies relative to the vacuum level. The effect of the xc-potential is to reduce the magnitude of the attractive QP screening cloud (in GW only the Hartree potential of the screening cloud is accounted for) and this lead to a significant up-shift of 0.3-0.5 eV for all energy levels leaving the band gap unchanged to within 0.2 eV. These effects should be important to include in order to correctly describe band alignment at surfaces and interfaces.
While the adiabatic xc-kernels discussed in this review provide an improved description of the short range correlations leading to the many positive derived effects described above, they still present shortcomings. Most importantly, they do not provide any improvements for strongly correlated systems such as Mott insulators. For such systems it might be necessary to include frequency dependence or develop kernels with a more sophisticated density dependence than the HEG based kernels considered here. Another short-coming is the failure to reproduce the correct wavevector-dependence in the limit of long wavelengths. This is not an important issue for total energy calculations, but for optical absorption spectre it is crucial to have the right limiting behaviour in order to capture excitonic effects. The deficiency stems from the fact that the renormalization scheme is based on the correlation hole of the homogeneous electron gas and a proper treatment of the long wavelength limit in insulators is likely to require a scheme, which is not solely based on the density. One possible route towards this is to note that the truncation factor θ(2k F − q) is the Fourier transform of the density matrix of a homogeneous electron gas with Fermi wavevector 2k F . In contrast to the bare density (defined by k F ), the density matrix provides a qualitative distinction between insulators and metals, but a direct truncation scheme based on the density matrix is not straightforward. Finally, the issue of self-consistency deserves a comment. In principle it is possible to define an optimized local effective potential (OEP) from the ACDFT with renormalized kernels, which may then form the basis of self-consistent solutions in the framework of the Kohn-Sham equations. However, the non-selfconsistent calculations are already rather computationally demanding and the OEP method does not seem to be a viable route to follow. In addition, the adiabatic kernels considered in the present work cannot be derived from any known ground state energy functional, which implies that it is somewhat inconsistent to evaluate the rALDA energy based on LDA orbitals for example. On the other hand, the situation is certainly improved compared to RPA where the Hartree orbitals comprises the only natural choice, but it would be highly desirable to have a class of xc-energy functionals that yield the renormalized kernels as their second functional derivative with respect to the density. A possible step in this direction is to define an LDA functional where the input density is replaced by a density convoluted with the Fourier transform of the truncation factor θ(2k F −q). This yield a ground an xc-energy functional, which leads directly to rALDA kernel if the density dependence of k F is neglected [19]. Although such a scheme is approximate, the construction itself provides an intriguing new route to the development of novel xc-energy functionals that incorporate non-local effects through weighted density averaging.