Singles correlation energy contributions in solids

The random phase approximation to the correlation energy often yields highly accurate results for condensed matter systems. However, ways how to improve its accuracy are being sought and here we explore the relevance of singles contributions for prototypical solid state systems. We set out with a derivation of the random phase approximation using the adiabatic connection and fluctuation dissipation theorem, but contrary to the most commonly used derivation, the density is allowed to vary along the coupling constant integral. This yields results closely paralleling standard perturbation theory. We re-derive the standard singles of G\"orling-Levy perturbation theory [G\"orling and Levy, Phys. Rev. A {\bf 50}, 196 (1994)], highlight the analogy of our expression to the renormalized singles introduced by Ren and coworkers [Ren, Tkatchenko, Rinke, and Scheffler, Phys. Rev. Lett. {\bf 106}, 153003 (2011)], and introduce a new approximation for the singles using the density matrix in the random phase approximation. We discuss the physical relevance and importance of singles alongside illustrative examples of simple weakly bonded systems, including rare gas solids (Ne, Ar, Xe), ice, adsorption of water on NaCl, and solid benzene. The effect of singles on covalently and metallically bonded systems is also discussed.


I. INTRODUCTION
In the last decade the interest in many body perturbation theory has risen significantly. This is to some extent related to the enormous increase in the available computer performance, but it is also driven by the realization that many of the presently available density functionals have limited predictive accuracy. Improving density functionals is a very active field of research in itself. In fact, tremendous progress has been made for semiconductor and insulator modelling by the inclusion of exact exchange, 1-5 as well as for weakly bonded systems by including either atom centered dispersion corrections [6][7][8][9] or non-local van der Waals corrections regarding the density at two points in space. [10][11][12][13][14] However, a unified comprehensive, accurate and predictive framework for metallic, covalently as well as dispersion driven interactions is hard to attain using present density functionals: most available density functionals require careful evaluation against more accurate methods before one can trust them to predict accurate numbers for a specific material.
In the quantum chemistry community, such a concise hierarchy of methods for evaluating and benchmarking more approximate methods, such as density functional theory, is well established. The highest rung of this hierarchy is made up by the full configuration interaction method, followed by a variety of more approximate methods. For instance, if the material is dominated by a single Slater determinant, the methods of choice are coupled cluster methods, 15 as well as Møller-Plesset perturbation theory 16 for large band gap systems. Only recently these methods have become available for solids. [17][18][19] For solids, calculations using coupled cluster methods or full configuration interaction methods are, however, exceedingly demanding approaching several 100.000 CPU hours for a single material with a few atoms in the unit cell. The approach taken in solid state systems is therefore often "bottom up", i.e. starting with an approximate scheme such as density functional theory and improving upon the description until results compare reasonably well with experiment. The random phase approximation (RPA) is one promising approach to achieve this goal. 20 In fact, the RPA yields a balanced description of most bonding types, including metallic bonding, covalent, ionic, as well as van der Waals (vdW) bonding. Initial applications were limited to small molecules. 21 Although first applications to bulk materials were not encouraging, 22 for many materials results nowadays surpass those from semi-local functionals. [23][24][25][26] The studies now span a wide range of applications, including molecular reactions, 27 rare gas solids, 24 properties of covalent, metallic and ionic solids, [28][29][30][31][32] dispersion forces in graphite and between graphene and surfaces, [33][34][35] layered compounds, 36 adsorption of molecules on surfaces, 37 bulk ice properties, 38 and many more applications are emerging. Recent advances include highly efficient implementations scaling only cubically with system size and linearly in the number of k-points, 39,40 as well as implementations scalable to massively parallel computers. 41 There is no denying that RPA is not perfect. Among all the possible many body diagrams, direct RPA exclusively sums the bubble diagrams. Attempts to include other kinds of diagrams, for instance higher or-der exchange interactions, 42,43 the contribution of single excitations [44][45][46][47][48] or approximate exchange and correlation contributions inspired by density functional theory 49 are currently vigorously explored research directions. Also better starting points than standard density functionals are explored. 50,51 Likewise, forces are yet only implemented in two molecular codes, 52,53 and they are not available in solid state codes.
The present paper mainly concentrates on the already mentioned singles, contributions that arise from determinants where one electron is moved from a state occupied in the initial Slater determinant to an originally unoccupied state. If the initial determinant is the Hartree-Fock determinant, then the singles are zero in lowest order perturbation theory because of Brillouin's theorem. 54 If one sets out from the Kohn-Sham determinant, the singles however contribute even in lowest order perturbation theory. In the adiabatic-connection framework this contribution was first derived by Görling and Levy. 55 In the present work, for reasons of consistency we will rederive the same term, albeit doing so using a concise Green's function notation, making the derivation algebraically simpler (see Sec. II C). The main point of our work is, however, that we give up the assumption that the charge density is kept constant along the coupling constant integration (see Sec. II A), an approximation made in most derivations of the adiabatic-connection fluctuation-dissipation theorem (AC-FDT). 20 This approximation was originally adopted, since density functionals were considered to be very accurate in predicting the groundstate density, and because AC-FDT was used as a theoretical pathway to derive improved density functionals. However, since AC-FDT is now a computational framework, and since density functionals do not always yield accurate densities, we feel that this point urgently needs to be revised. We are aware of at least two papers where this assumption has been given up as well, albeit in the first one the point was only made in passing 56 and the second one only considers small prototypical molecules. 57 Our present formal derivation (see Sec. II D) yields results similar to the singles and re-normalized singles proposed by Ren et al. 44,45,48 Ren and coworkers, however, formally gave up the adiabatic-connection framework and used standard second order perturbation theory to motivate their singles and re-normalized singles. We feel that a concise perturbational framework (here the adiabatic-connection framework) is helpful to better understand the underlying approximations. We will discuss that the singles account for the changes of the mean field density matrix along the coupling constant integral. This is exactly analogous to coupled cluster theory or, since coupled cluster theory is just a re-summation of certain Goldstone diagrams, standard many body perturbation theory. This insight explains why the inclusion of singles increases the bonding between weakly interacting fragments, i.e. atoms and molecules. When the singles are included the charge density contracts compared to the original density functional, and this results in a de-crease of the Pauli repulsion. We demonstrate this effect here for rare gas solids, as well as, the cubic phase of ice, the benzene crystal and water adsorption on NaCl. For these systems, the predicted cohesive energies and lattice constants are, after inclusion of the singles, in very good agreement with experiment (see Sec. IV A,  IV B, IV C, IV D). We also present results for covalently bonded as well as metallic systems. Here no improvements are discernible, or rather, the corrections of the singles are mostly tiny and do not worsen the already excellent agreement with experiment (see Sec. IV E).

II. THEORY
A. Adiabatic switching with density change In this section, we give a brief derivation of the adiabatic switching, where the correlation energy is obtained by switching from the single determinant reference system to the many body system of interest. We consider the usual adiabatic connection, where the density is kept constant along the pathway, as well as adiabatic switching allowing for a change of the density.
In the adiabatic-connection framework, it is common practice to switch from the purely local, in real space multiplicative potentialV λ to the exact many electron HamiltonianĤ λ=1 : Herev is a two-electron operator, typically the Coulomb kernel 1/|r − r ′ |, andĥ is an arbitrary one electron operator, typically the kinetic energy operator ∆ and the potential of the ionic cores or some other external oneelectron potentialV ext where we have assumed atomic (Hartree) units. The correlation energy E tot c is defined as the difference between a single Slater determinant Ψ 0 evaluated using the orbitals at zero coupling (λ = 0) and the exact many-electron wave function Ψ 1 evaluated at full coupling (λ = 1): At zero coupling, the Hamiltonian is chosen to be the usual Kohn-Sham Hamiltonian, whereV KS 0 is the standard Kohn-Sham potential (we keep the subscript 0 to indicate zero coupling). Note that the Kohn-Sham potential includes exchange, correlation, and the Hartree contributions. In the following, we will also use the shorthand |n = |φ n , and the indices i, j, k, ... always refer to occupied one-electron orbitals, whereas the indices a, b, c, ... refer to virtual (unoccupied) one-electron orbitals of the Kohn-Sham Hamiltonian.
Subtracting and addingV KS 0 to the second term on the r.h.s of (2) one can rewrite the correlation energy as In going from the first to the second equation, the "Hellman-Feynman" theorem has been used, i.e. it is assumed that Ψ λ is the groundstate wave function of the HamiltonianĤ λ , and terms involving the derivative of the wave function d Ψ λ /d λ thus vanish. The first term in the second line is exactly the Hartree and exchange energy evaluated for Kohn-Sham orbitals. Inserting the derivative of the Hamiltonian d H λ /d λ immediately yields (compare Equ. (1)) We note that a similar term is given in the appendix of Ref. 56 without a clear derivation. The most common way to realize the coupling integral is to chose the poten-tialV λ such that the density remains constant along the coupling constant integral. 20 Then the expectation value of the density operator yields the (constant) groundstate density Ψ λ |n(r)|Ψ λ = n(r), and the first term in the second line can be integrated yielding canceling the last term in the second line of Equ. (5).
Since the additive potentialV λ must be zero at full cou-plingV 1 = 0, one obtains the standard equation for the AC-FDT: which is only valid if the density is constant along the integration pathway, whereas the full version is obviously given by Equ. (5).
To obtain a simplified version of the full equation, one can switch off the Kohn-Sham potential linearly i.e.
Linear switching was first considered by Harris and Jones, 58 and subsequently discussed in Refs. 56 and 57.
Linear switching is also exact (as long as no phase transition is encountered) and yields results identical to standard Rayleigh-Schrödinger perturbation theory where the perturbation is also switched on linearly. Linear switching yields for the second line in Equ. (5) a simple correction to the correlation energy Here n λ is the charge density at coupling λ. The total correlation energy is then given by the sum of Equ. (6) and Equ. (8).

B. Fluctuation-Dissipation theorem
It is common to rephrase the correlation energy in Equ. (6) using the fluctuation-dissipation theorem. The derivation is sketched in appendix A, however, similar derivations can be found elsewhere. 20,56,57,59 The final result becomes: Here we have introduced short-hands In the equations above, χ λ is the reducible polarizability (or density fluctuation response function) of the many electron system at coupling λ: T is the time ordering operator, and δn the density fluctuation operator δn =n − n λ . The polarizabilities are evaluated at a small negative infinitesimal imaginary time τ → 0 − . Note that throughout the paper, τ is the imaginary part of an imaginary time t = −iτ . The polarizability in imaginary time can be related to the polarizability in imaginary frequency iω by a Fourier transformation: At zero coupling, χ 0 becomes the independent particle polarizability P 0 of the Kohn-Sham system, and can be written in terms of the one particle Green's function 56,60 We define the independent particle polarizability at coupling λ generally as where G λ is the one-particle Green's function of the λ interacting systemĤ λ . At zero coupling, i.e. for the Kohn-Sham case, the Green's function in imaginary time G 0 are given by (µ is the chemical potential of the electrons, i are occupied and a are unoccupied one electron orbitals): Before continuing, we note that at small λ and neglecting fluctuations, the term in the integral on the second line in Equ. (9) can be also written as an integral of the exchange energy where γ λ (r, r ′ ) is the one-particle density matrix at coupling λ. This follows from expanding the δ function in one-electron orbitals δ(r, r ′ ) = m∈all r|m m|r ′ , and inserting the independent particle approximation for the polarizability (compare Equs. 13 and 14). Hence, in Equ. (9) the first line represents the change of the Hartree energy along the coupling constant integration, whereas the second line accounts for the change of the exchange energy (for uncorrelated wave functions). Returning now to Equ. (13), we note that replacing χ λ by P λ in Equ. (9) neglects important many body effects related to changes in the Hartree or exchange potential. In many body perturbation theory, these can be included exactly by solving the Bethe-Salpether equation for the polarization propagator. [61][62][63][64] The simplest approximation to the BSE equation is the common RPA approximation which only includes the Hartree-related ring or "bubble" diagrams. In time dependent DFT, the related equation isχ where f xc (λ) accounts for all correlation effects, including the change of the independent particle polarizability along the coupling path. If the density is kept constant (n λ = n 0 ) and one approximates P λ (ω) = P 0 (ω) in Equ. (16), one obtains for Equ. (9) the direct RPA correlation energy: 65 Here and in the following, we often drop position indices and implicitly assume integration over the missing spatial coordinates: C. Singles contribution in Görling-Levy perturbation theory Before considering density changes along the coupling constant integral, we will briefly derive the singles expression in the AC-FDT for a fixed charge density in order to compare with this equation later.
In standard RPA, one neglects that the independent particle polarizability changes along the coupling constant integration, as the one-particle Green's function changes as λ changes (see Equ. (19) below). Oddly this term is not often considered in the AC-FDT framework, although Görling and Levy have already highlighted its relevance (albeit not in a Green's function formalism). A partially analogous discussion using the Green's function formalism can be found in e.g. Ref. 56 and 66. At small coupling λ, the change of the one-particle Green's function is described by the first order term in the Dyson equation (see e.g. Ref. 67): where ∆V KS (r)δ(r, r ′ ) is the change of the local Kohn-Sham potential, andV x (r, r ′ ) is the exact non-local exchange potential evaluated using the DFT orbitals at coupling constant zero: this term originates from switching on the exact many body potential, which in lowest order is equivalent to switching on the non-local exchange. If one performs the coupling constant integral keeping the density fixed, then the change of the Kohn-Sham potential must be chosen such that the density remains exactly constant to the original density at λ = 0. This requirement can be written as: Since the Kohn-Sham potential is local, one can factor outP 0 (ω = 0) = dτ G 0 (τ )G 0 (−τ ) (see Equ. (12)) and obtain: This is just the local exact-exchange optimized-effective potential (EXX-OEP) V EXX . In summary, if the density is supposed to remain constant along the coupling constant, the first order change of the local Kohn-Sham potential is exactly given by the exact-exchange OEP potential. This is not surprising, since the defining property of the Sham-Schlüter equation is that the density from the non-local exchange potential and the effective local potential must equal each other.
One can now calculate the change of the correlation energy by inserting the first order change of the Green's function into the expression for the correlation energy in the AC-FDT. Approximating χ λ → P λ = G λ G λ , inserting (13) and (19) into the second line of Equ. (9), and rewriting the δ function as a sum over states, or alternatively, starting from Equ. (15) and identifying γ λ = G λ (0 − ) yields a change of the correlation energy of (yet omitting the integration over λ): Because of the trace, the second term yields just the complex conjugate of the first term. Furthermore, the integral over the positive half-plane of τ gives the same value as that over the negative half-plane (τ < 0). Finally, the term −G 0 (0 − )v corresponds to the non-local exchange potentialV x . Inserting the defining equation for the one-particle Green's function (14) and performing the integration over τ and λ yields what is commonly referred to as the singles contribution: Note that one needs to add Equ. (21) to Eq. (23) multiplied byV EXX (r) and integrated over r to derive this convenient equation. The corresponding (time-ordered) Goldstone diagrams are also shown schematically in Fig.  1. A few comments are in place here. The term has been first derived by Görling and Levy. 55 In their derivation it is, however, not obvious that this term describes the changes of the correlation energy from the one-particle density matrix (which equals the one particle Green's function at τ = 0 − ) along the coupling constant integral. In fact, most standard AC-FDT calculations neglect the term. Only the RPAX (RPA including eXchange) following Görling and coworkers accounts for this term (often referred to as EXX-RPA). 69 In these formulations, however, the change of the Green's function is accounted for by recasting it (as well as the particle-hole ladder diagrams) into an effective exchange kernel f xc (λ) ≈ λf x for the polarizability (compare Equ. 17). 66,69 We conclude, the singles account for the change of the density matrix along the coupling constant integral. For a constant density, obviously only changes of the exchange and kinetic energy can be included. In the next step, we will also allow for changes of the charge density along the coupling constant integral.

D. Singles contribution with density changes
We now derive the singles contribution to the correlation energy for the case when the density does not stay constant during the coupling constant integration. Changes of the charge density are most easily accounted for by linearly switching off the Kohn-Sham potential (see Equ. (7)) and linearly switching on the exact many body interaction. In principle, this makes the derivation even less involved, since the determination of the local exact exchange potential is no longer required. We first again derive the expression in lowest order, where the Green's function is now given by In the lowest order, the change of the potential is now the difference between the Hartree-Fock potentialV HF = V H +V x , the sum of the Hartree and exchange potential, and the original Kohn-Sham potential, which is adiabatically switched off. As before, the change of the density matrix is given by the change of the Green's function G(τ → 0 − ) (compare Equ. (20)). In this case, the first and second line of Equ. (9) can be combined to yield: The first line of Equ. (9) yields the Hartree-potential times the change of the density, whereas the second line yields the exchange potential times the change in the density matrix, in sum the change of the density matrix timesV HF .
After performing the coupling constant integral this yields This term needs to be combined with the term ∆E c given in Equ. (8), which can be calculated by inserting Equ. (19) and performing the τ and λ integration. Both contributions combined yield a simple term which exactly corresponds to the singles suggested by Ren et al. 44,45 Here, we have performed the derivation concisely within the AC-FDT framework instead of Rayleigh-Schrödinger perturbation theory, and, as it must be, both are exactly equivalent. As in the previous paragraph, the singles account for the change of the mean field exchange energy. However, now they also include the change of the mean field Hartree energy along the coupling constant integration. Here and in the following, we define the mean field as the contributions that stem from the one-particle Green's function and the related density matrix.
In the renormalized singles of Ren and coworkers 48 also higher order contributions are accounted for. However, it is not entirely straightforward to generalize our results to include higher order terms in λ, and to continue, we make one crucial approximation. Let us introduce this approximation for the density term, which can be written as (compare Equ. (9)) In the second line, we assume that the differences between n λ and n 0 are small so that we can approximate the sum by 2n 0 . Analogous manipulation is possible for the term involving the polarizability, if the full polarizability is approximated by the independent particle approximation Equ.
After collecting all terms, adding Equ. (8), and noting that n(r) = G(0 − , r, r), one obtains the approximate renormalized singles correlation energy It is fairly straightforward to backtrack this into a simple total energy relation (essentially reversing the steps introduced in Sec. II A). We first note that the difference between the Hartree-Fock potential and the Kohn-Sham potential equals the difference in the corresponding oneelectron Hamiltoniansĥ HF −ĥ KS . Next, the constant term G 0 (0 − ) is integrated over so that we obtain As the one-electron Green's function G λ (0 − ) is the exact Green's function of the one-electron Hamilton operator h λ =ĥ KS + λ(ĥ HF −ĥ KS ), one can rewrite the first line as (Hellman-Feynman theorem): Combining this with the second line in Equ. (31) Here γ HF = G HF (0 − ) is the Hartree-Fock density matrix, determined for the Hartree-Fock Hamiltonianĥ HF , where the Hartree-Fock potential is set up with DFT-orbitals. This is exactly the "single-shot" Hartree-Fock energy: it can be calculated by diagonalizing the HF Hamiltonian (set up with DFT-orbitals), summing the eigenvalues of the occupied states and subtracting the original diagonal part of the HF Hamiltonian for the occupied manifold evaluated using the original DFT orbitals (i DFT ): This prescription is one central result of the present paper.
In all practical calculations, we found the value of the single shot HF energy to be exceedingly close to the renormalized singles introduced by Ren and coworkers. 48 For diamond the difference is on the order of 0.3 meV (with the singles being of the order of 0.3 eV). Even for small band gap systems, such as metallic Pd or Al, differences hardly ever exceed 1 % and are entirely irrelevant when evaluating relative energies or lattice constants. The relation of our equation to the singles of Ren et al. is fairly straightforward to see. Ren et al. essentially renormalizes the propagator G in the occupied-occupied block as well as the unoccupied-unoccupied block, by diagonalization of these sub-blocks using the HF Hamiltonian h HF . This step is crucial, since the one-electron eigenvalues are renormalized to the HF-eigenvalues; if LDA eigenvalues were used in the evaluation of the singles in Equ. (28), the response of the system to the change of the potential from Kohn-Sham to HF would be strongly overestimated. Ren then calculates the change of the mean field energy in second order with the DFT eigenvalues in Equ. (28) replaced by the renormalized HF eigenvalues. As opposed to this prescription, we also "renormalize" the propagator in the occupied-unoccupied block. With the present derivation at hand, there is no obvious reason why not to chose the simpler prescription of the present work. Since results obtained using the renormalized singles (rSE) are in practice indistinguishable from results obtained using the single-shot HF energy change, all calculations here use the single-shot HF energy change, but are nevertheless labeled as "rSE".

E. Singles contribution beyond the Hartree-Fock description
An obvious extension to the renormalized singles approach is to use the full RPA density matrix instead of the HF density matrix to estimate the change of the mean field energy: where the γ RPA is the RPA density matrix. In the singles derived in the previous section, it was assumed that the density matrix of the interacting system is well approximated by the Hartree-Fock case, which seems a crude approximation. Given the excellent performance of the RPA for total energies and band gaps, the RPA density matrix, however, should approximate the density matrix of the real interacting system very accurately. At this point, we have, however, no entirely concise derivation for the term, although the physical interpretation is clear: it should account approximately for the change of the mean-field kinetic, Hartree and exchange energy along the coupling constant integral.
To evaluate the RPA density matrix, we calculate the RPA Green's function Here Σ is the self-energy in the GW or random phase approximation (the two approximations are synonymous) and G 0 and W are the Kohn-Sham Green's functions (14) and the screened potential evaluated using Kohn-Sham polarizabilities. 60V KS xc is the exchange-correlation contribution to the Kohn-Sham potential.
To evaluate the RPA density matrix numerically, we determine the interacting Green's function at full coupling G RPA (iω) using Equ. (35), transform it to the imaginary time at a small negative infinitesimal τ → 0 − to obtain the one-particle density matrix for the RPA, γ RPA (r, r ′ ) = G RPA (0 − , r, r ′ ), and finally diagonalize the density matrix to obtain the natural orbitals In the present work, we use one more crucial approximation: instead of using the occupancies of the actual density matrix, we keep the original occupancies evaluated on the level of density functional theory (1 and 0 for insulators). This has two reasons: first the density matrix evaluated from the Green's function (35) is not particle conserving, i.e. the number of electrons deviates from the original number. 67 Only if the Green's function in (35) and (36) is iterated to self-consistency, the particle number is conserved. To test our present code, we have also iterated the Green's function in both equations to selfconsistency, and in that case, the electron number is indeed exactly conserved. However, such calculations are fairly expensive and difficult to apply routinely. They would also most likely require us to combine it with a different treatment of the fluctuation terms beyond the standard RPA treatment as used here. 67 The second reason is based on the definition of singles in quantum chemistry. A density matrix with occupancies 0 and 1 corresponds to a single Slater determinant. We, therefore, approximate the density matrix by the "best" single Slater determinant approximating the correlated RPA density matrix. This is consistent with Brueckner coupled cluster orbitals, 70 which are obtained by performing a rotation between the occupied and unoccupied manifold to determine an optimal reference single Slater determinant. The rotation is chosen to remove all "singles" contributions from the correlation energy. In quantum chemistry, the fluctuations that cause the fractional occupancies in the density matrix are, in fact, not included in the singles. For instance, in the coupled cluster theory, fluctuations are accounted for by double excitation operators (doubles). Per definition, singles are diagrams ending in a single excited determinant with one hole in an orbital i and an additional electron in an orbital a (compare Fig. 1). We essentially follow the quantum chemistry convention in partitioning the correlation energy into a "singles" term and fluctuation terms described entirely by Equ. (18). We finally note that the RPA orbitals constructed in this manner are similar but not identical with the Brueckner RPA direct ring coupled-cluster orbitals suggested by Moussa. 71

III. COMPUTATIONAL SETUP
The present calculations were performed using the VASP code. 72 We used the new implementation of the RPA routines as documented in Refs. 39 and 40. The code has been extended to allow for (self-consistent) calculations of the one particle Green's function in the random phase approximation i.e.
solving Equ. (35) and (36). Here only single shot calculations are performed, by calculating the self-energy once at full coupling and determining the RPA Green's function and RPA density matrix (Equ. (35)). The singles contributions are evaluated either using the single shot HF density matrix (Equ. (32)) or the single shot RPA density matrix (Equ. (34)). We term the two results either rSE (renormalized singles) or GWSE, respectively. For the fluctuation term we use the standard random phase approximation as defined in Equ. (18). Furthermore, we

A. Rare gas solids
Rare gas solids constitute a prototypical test case for van der Waals bonded solid state systems. Although, the standard RPA performs reasonably well for rare gas solids, 24 one observes that the binding energy is quite significantly underestimated, in particular for He. Inclusion of singles has remedied this issue for the rare gas dimers, and one would expect that this also applies to solids. 44,48 We first note that quantum chemical coupled cluster calculations at the level of CCSD(T) (coupled cluster with singles and doubles and perturbational triples) using the method of increments and including up to four body interactions yield essentially exact results within a few µHartree. 75,76 To compare with experiment we have corrected the experimental data for zero point vibration energies in both the cohesive energy, as well as in the lattice constants. 75 The present calculations have been performed using relatively high plane wave energy cutoffs of 1000 eV, 600 eV and 500 eV for Ne, Ar and Kr in order to obtain smooth energy-volume curves. 4×4×4 k-points were used, although already 3 × 3 × 3 k-points yield practically identical results. The PAW potentials are approximately norm-conserving to avoid errors in the vdW contributions from excitations into high lying unoccupied orbitals. 77 As already mentioned, standard RPA combined with exact exchange yields at best modest agreement with experiment (see Tab. II). Specifically, the equilibrium volumes are overestimated and the binding energies are, as already noted above, about 20 % too small, with the errors being particularly large for Ne. Inclusion of the singles contributions from HF (rSE) and using the RPA density matrix (GWSE) yields a clearly improved agreement with experiment in particular for Ar and Kr. Ne is less satisfactory. For Ne, the binding energy significantly overshoots if the singles are evaluated on the Hartree-Fock level (rSE). This improves, if the RPA density matrix is used, however, including GW-singles worsens the volume even compared to RPA. We believe that this is a result of the PBE approximation being particularly inaccurate for large band gap system such as Ne. This is for instance exemplified by the fact that the binding energy almost doubles, when the exact exchange evaluated using DFT orbitals is replaced by the exact exchange evaluated using Hartree-Fock orbitals (HF+RPA). For Ar and Kr the changes are typically only 25 %. Also other metrics indicate that the error in the orbitals is significant using the PBE functional for this case. Hence, single shot corrections, be it RPA+rSE or RPA+GWSE, yield less reliable results then for Ar and Kr. For Ar and Kr, the performance of RPA+rSE and RPA+GWSE is remarkably good, approaching that of state of the art quantum chemistry methods. Furthermore, the differences between rSE and GWSE are small (except for Ne), which is related to the fact that all these systems screen relatively weakly. Therefore, the RPA density matrix is very close to the HF density matrix (see below).

B. Ice
Ice is another system, where singles are expected to have a significant impact. As for rare gas solids, the PBE charge density is too spread out and replacing the exact exchange evaluated using PBE orbitals by the exact Hartree-Fock exchange increases the binding energy by 100 meV. 38 The predicted volumes using EXX+RPA and HF+RPA almost exactly bracket the experimental values.
Here we only concentrate on one ice phase, the lowest energy cubic phase of ice, I c (a), with a ferroelectric order. In our previous study, we found this phase to be practically iso-energetic with the ferroelectrically ordered ice XI in the Cmc2 1 space group. Ice XI is a proton ordered variant of the common proton disordered phase of hexagonal ice. In the present study, we used PBE relaxed structures and identical potentials as in our previous study. 38 However, the cutoff was increased from 450 eV to 600 eV. This generally results in smoother energyvolume curves and changes the calculated volumes by about 0.5 %: because of the increased cutoff and reduced noise in the calculated data, the present calculations are slightly more accurate.
As already observed in our recent work, the EXX+RPA underbinds, whereas combining the RPA with Hartree-Fock energies overestimates the binding energies (see Tab. III). Including the singles in the HF approximation improves the description already significantly, although the results are still too close to the HF case. In this case, the singles evaluated using the RPA-density matrix (GWSE) yield results very close to experiment and practically identical to the diffusion Monte Carlo data, which gives an equilibrium volume per molecule of 31.7Å 3 and a binding energy of 605 meV for hexagonal ice. 78 Since hexagonal ice has a 0.2Å 3 smaller volume and a 5 meV reduced binding energy due to disorder (compare Table III in Ref. 38), the results are virtually on top of the DMC data at a tiny fraction of the compute cost.

C. Benzene
Benzene is the simplest aromatic molecule and therefore it represents an interesting reference point for the study of more complex molecular solids. The delocalized character of π bonds in benzene makes an accurate calculation of interactions or cohesive energies difficult, both, in the gas phase, as well as in the condensed phases. Simple methods, such as MP2, lead to overestimated binding energies and more involved schemes need to be used to obtain accurate results. For example, for crystalline benzene Wen and Beran computed a post-MP2 correction of 10.4 kJ/mol (108 meV), 80 reducing the MP2 cohesive energy considerably. Li et al. estimated the RPA cohesive energy to be 47 kJ/mol or 487 meV, 81 significantly below the recently revised experimental estimate of −55.3 ± 2.2 kJ/mol (about 573 meV). 82 We have used the geometry of the crystal and monomer optimised with the optB88-vdW functional and performed both the extrapolation to infinite k-point mesh and infinite basis sets. To obtain the reference energy of the monomer, we also performed calculations at different volumes and extrapolated to the infinite cell volume. The RPA and EXX calculations used up-to 3×3×3 and 6×4×6 k-point meshes, respectively.
The data are summarised in Table IV  ing energy increases by 40 meV for the HF based singles and by about 55 meV when we use the GW singles. This means that the final error is halved compared to the original RPA calculation. The hybrid scheme, where self-consistent Hartree-Fock is used for the mean field part, gives results in a better agreement with the reference data, underestimating it by about 43 meV. We note that such a performance is in agreement with the results of Ren et al., 47 who found that RPA underestimates interaction energies for benzene dimers both in the stacked and in the T-shaped geometries. Moreover, including single excitations did not improve the results substantially. This clearly shows that singles have limitations. We believe that the main error in this case stems from an inaccurate description of the delocalized π electrons in PBE.

D. Water on NaCl
To explore the accuracy of the RPA and the singles corrections for adsorption on surfaces, we have calculated the adsorption energy of a water molecule on the NaCl surface. This is a prototypical system and one for which an estimate of the adsorption energy of water has been obtained. 83 However, this reference data was calculated using an embedded cluster approach to obtain the correlation energy. This makes a direct comparison to our calculations, that are necessarily a finite coverage, difficult. Therefore, we opted not to use this reference data and instead obtained an estimate from MP2 calculations for a small system. To reduce the computational cost of MP2, we, furthermore, restrict the study to a small p(1×1) supercell with a single water molecule and two layers of NaCl and use that consistently for both MP2 and RPA. This corresponds to a high coverage, with one molecule per one surface sodium atom. Moreover, the reference surface and isolated molecule have the same geometry as in the adsorbed structure, and the same simulation cell is used for all cases to allow for efficient error compensation.
The interaction energy depends only weakly on the cutoffs chosen for orbital and auxiliary plane-wave basis sets. For MP2 we set the cutoff for the orbital basis and auxiliary polarizability basis to 350 eV and 450 eV, respectively. The interaction energy depends also only weakly on the k-point sampling, we have used up to 3×3×1 k-points (the cutoff for the auxiliary basis in this case was 250 eV). After accounting for basis set convergence and k-point convergence, we obtained an estimate of the molecule-surface interaction energy of E int = −420 meV. Corrections beyond MP2 are expected to be small. In the work of Li et al. 83 , CCSD(T) calculations lead to a 10 meV stronger binding energy than MP2. Further corrections will arise from the use of hard PAW potentials, for example, the RPA binding increases by 15 meV when small core potentials are used, but since we compare MP2 and RPA with similar setups this is irrelevant in the present case.
Our results are summarized in Table V. As a reference, we use our MP2 value corrected with the post-MP2 correction of Li et al., which yields, in total, a binding energy of −430 meV. As expected, RPA underestimates this value, by about 50 meV. Adding the singles correction (rSE) leads, in this case, to a perfect agreement with the estimated reference data. This is also in agreement with the calculations of Ren et al. performed on the S22 test set. 47 They found a very good agreement between the reference data and RPA+rSE for those dimers in S22 that are bound by hydrogen bonds. The hybrid HF+RPA overestimates the reference interaction energy, which is also in agreement with the findings of Ren et al. In this case, the singles in the GW approximation (GWSE) do not quite work as well, but the agreement with the reference data is still reasonable.

E. Covalent and metallic solids
To evaluate the influence of the singles on the lattice constants of covalent and metallic solid state systems, we show the equilibrium lattice constants for selected materials in Tab. VI. We also use the present opportunity to evaluate whether improved PAW potentials have an effect on the equilibrium lattice constants. As shown in one of our recent work, 77 quasiparticle energies as evaluated in the GW approximation can have large errors, since the PAW projectors possibly do not span the unoccupied orbital space accurately. As a remedy to this problem, we have suggested to use PAW potentials with norm-conserving partial waves. These unfortunately increase the computational cost, sometimes, even quite significantly. In Tab. VI, the first column reports the lattice constants evaluated using such norm-conserving PAW potentials. In the present calculations, to attain the highest possible accuracy, the entire lower lying core shell was included in the correlated calculations, except for oxygen, carbon, nitrogen, fluorine (2p elements) as well as phosphorus and chlorine. For instance for Si and Al, the 2s and 2p states were treated as valence states, for In and Sb the 4d, 4s and 4p states were fully included. The calculations were performed using 6 × 6 × 6 k-points and 8 × 8 × 8 k-points for gapped systems and metals, respectively. Increasing the k-point set for selected semiconductors and insulators from 6 × 6 × 6 to 8 × 8 × 8 changed the lattice constants by less than 0.1 % (C, Si, SiC, LiCl). For metals, an increase of the k-point set to 10 × 10 × 10 k-points changed the lattice constants also only by typically 0.1 % (the results for the transition metals and std-PAW are reported for 10 × 10 × 10 k-points). This suggests that the lattice constants are kpoint converged to about 0.1 %. Errors incurred by the finite plane-wave basis set are of the same order, so that we estimate the accuracy of the present calculations to be about 0.2-0.3 % in the lattice constants (or 1 % in the volumes).
For the RPA the mean relative error with respect to the zero-point corrected experimental lattice constants is just 0.06 % in the present calculations, whereas the mean absolute error is about 0.3 %. We note that this is within the estimated error bars of our calculations. It is therefore futile to seek for any systematic errors: the RPA seems to be able to predict lattice constants in almost perfect agreement with experiment. Only for the Na metal, the lattice constant is obviously significantly underestimated (excluding Na from the calculations, the MRE and MARE drop to 0.01 % and 0.25 % respectively). There are very few density functionals that yield a similar accuracy. In fact, the present results slightly surpass those for the HSEsol functional. 84  The most commonly used functional, the PBE functional, overestimates the lattice constants by about 1.3 %, and even the PBEsol functional yield a mean absolute relative error of 0.46 % for a slightly larger set. 84 The use of not norm conserving PAW potentials increases the lattice constants, on average by 0.3 %. Also the mean absolute relative error increases slightly from 0.3 % to 0.35 %. For most elements, the differences between standard and NC PAW potentials are small. However, they can approach up to 0.4 % for elements with 3d and 4d semi-core states and (AlAs and GaAs) and up to 1 % for ionic compounds with strongly polarizable cores (Na). Note that in ionic solids, vdW interactions are sizable, since the Na 2s and 2p core electrons are hardly screened and interact via van der Waals interactions with the neighboring halide atoms.
The origin for the increase in the lattice constants from the NC PAW potentials to the standard PAW potentials is that the standard potentials underestimate the polarizability of the core and this in turn yields too large lattice constants. For most applications, this small error of the lattice constants should be acceptable, however. We finally note that the present values for the standard potentials are also in good agreement with our first publication, 24 although the potentials have been improved since our previous calculations published in 2008. Specifically, the present set of standard potentials preserves the norm of the pseudo-orbitals better (although not perfectly, as the NC PAW potentials do), which in turn increases the core polarizability and decreases the lattice constants compared to the original values in Ref. 24.
Because the compute cost for the NC PAW potentials is very high, we have evaluated the singles only for standard potentials. The important result is that the singles hardly change the lattice constants, except for some transition metals where an increase by 1 to 2 % is observed for Cu and Pd. Fig. 2 shows that singles indeed shift significantly the equilibrium volume to larger values, clearly worsening the agreement with experiment. GWSE rectifies this, and yields almost identical values to the RPA for all considered elements. As we will discuss in the next paragraph, the RPA density matrix for metals is seemingly very different from the HF density matrix, and most likely close to the DFT density matrix. This suggests that HF singles are not adequate for metals, whereas, we expect GW singles to be accurate across all systems. Fig. 3 indicates how close the GW density matrix is to the HF density matrix. In order to measure this, one has to introduce a metric to sensibly determine the difference. An obvious choice is the total energy difference between the Hartree-Fock energy evaluated using the HF density matrix (rSE) and the Hartree-Fock energy evaluated using the RPA/GW density matrix To present the values in a concise manner, we divide this by where γ DFT is the DFT density matrix and finally take the square-root. The reason for including the squareroot is that the functional is variational and therefore quadratic around γ HF , since γ HF is the groundstate density matrix ofĥ HF . If the value is 0, the RPA density matrix coincides with the HF density matrix, whereas for 1 it is closer to the DFT density matrix. One clearly recognizes that the RPA density matrix is generally quite close to the HF density matrix for light elements and insulators. However, for metallic materials Na, Pd, Cu, Rh, and Ag, as well as for small gap semiconductors and semiconductors with heavy elements (Ge, Ga and In compounds) this is not the case. The variational properties, discussed above, suggest that an evaluation of the mean field contribution using HF orbitals will be generally superior to an evaluation using DFT orbitals. In particular, for large gap systems, such as rare gas solids, ice, but also C, Si and SiC, LiF and MgO, the differences between rSE and GWSE mean field contributions are tiny and only of the order of 10 % (essentially the square of the distance shown in Fig. 3). For metals, this is, however, clearly not the case, and already shown in the previous paragraph, erroneously increases the lattice constants compared to RPA or RPA+GWSE.

V. DISCUSSION AND CONCLUSIONS
The present work is devoted to the performance of the random phase approximation for extended systems if singles contributions are taken into account. The first part of the paper focuses on the derivation of the singles within the adiabatic-connection fluctuation-dissipation framework. Not unexpectedly, this derivation yields exactly the same contributions as the singles originally suggested by Ren and coworkers, 44 because standard Rayleigh-Schrödinger perturbation theory and coupling-constant integration are identical. The coupling-constant integration and the formulation used here has, however, the advantage that it gives a very clear picture of what the singles describe. They account for the "mean field" energy changes from the non-interacting Kohn-Sham reference system to the interacting system, where we define mean field as those contributions arising from a changes of the one-electron density matrix. This makes it very clear why cohesive energies are increased when going from the DFT mean field description to the HF mean field description: the latter contracts the orbitals and thus reduces the Pauli repulsion between the atoms or molecules.
Renormalized singles can be also derived in the present framework and the final equation for them is particularly revealing (compare Equ. (32)). The "renormalized" singles describe the energy difference between the Hartree-Fock eigenvalues and the diagonal of the HF matrix evaluated using DFT orbitals (compare Equ. (33)). This is not exactly identical to the renormalized singles suggested by Ren et al. 48 although we found, in practice, that our simpler equation gives virtually the same results as Ren's renormalized singles and, as a bonus, it is trivial to implement and most likely already available in most codes.
Inspired by the simple physically transparent form of the singles, we have also suggested an alternative form for the singles that relies on the RPA-density matrix instead of the HF-density matrix (compare Equ. (34)). We have termed this correction GWSE, singles in the GW approximation. Clearly this description should be superior to the simple HF description, as the RPA density matrix should be fairly close to the exact groundstate one-particle density matrix. In practice, for large band gap systems, such as rare gas solids, ice and adsorption of water on NaCl, differences between the rSE and GWSE are small. This suggests that the HF density matrix is often astonishingly close to the RPA density matrix for systems with light atoms and large band gaps (compare Fig. 3). In such cases, the rSE can be used instead of the GWSE with little loss of accuracy. We believe this explains why the rSE approximation was so successful for molecules. For metals and heavier elements, the approximation becomes increasingly worse and an erroneousness increase of the lattice constants is observed for some transition metals for the rSE approximation, which is rectified by the GWSE.