Beyond the Random Phase Approximation for the Electron Correlation Energy: The Importance of Single Excitations

The random phase approximation (RPA) for the electron correlation energy, combined with the exact-exchange energy, represents the state-of-the-art exchange-correlation functional within density-functional theory (DFT). However, the standard RPA practice -- evaluating both the exact-exchange and the RPA correlation energy using local or semilocal Kohn-Sham (KS) orbitals -- leads to a systematic underbinding of molecules and solids. Here we demonstrate that this behavior is largely corrected by adding a"single excitation"(SE) contribution, so far not included in the standard RPA scheme. A similar improvement can also be achieved by replacing the non-self-consistent exact-exchange total energy by the corresponding self-consistent Hartree-Fock total energy, while retaining the RPA correlation energy evaluated using Kohn-Sham orbitals. Both schemes achieve chemical accuracy for a standard benchmark set of non-covalent intermolecular interactions.

The concept of cRPA dates back to the many-body treatment of the uniform electron gas in the 1950's [1,2], and was later formulated [20] within the context of density-functional theory (DFT) [19]. Recent years have witnessed a revived interest in EX+cRPA and its variants in quantum chemistry [3][4][5][6][7][8][9], solid state physics [10][11][12], and surface science [13][14][15]. Within the framework of Kohn-Sham (KS) DFT, EX+cRPA embodies an orbital-dependent functional that can in principle be solved self-consistently via the optimized effective potential approach [21]. This is however numerically very demanding, and practical EX+cRPA calculations are commonly performed in a post-processing fashion, where single-particle orbitals from a self-consistent DFT calculation in the local-density approximation (LDA), generalized gradient approximations (GGAs), or alike, are used to evaluate both the EX and cRPA terms. Alternatively, one can formulate cRPA in terms of many-body perturbation theory (MBPT) based on a Hartree-Fock (HF) reference.
Throughout this Letter we will adopt the following nomenclature: E F @SC is the total energy of the functional F , evaluated with the orbitals of a self-consistent (SC) scheme, e.g., HF, or the Perdew-Burke-Ernzerhof (PBE) [22] GGA. The corresponding theoretical scheme is then labelled as F @SC. We also use the letter "x" or "c" in front of F or as a subscript of E F to refer to the exchange or correlation part of the scheme explicitly. The functional F can be exact-exchange (EX), or additionally contain the RPA correlation (EX+cRPA), etc. For instance, E EX @HF is the self-consistent Hartree-Fock energy, whereas the conventional RPA scheme based on PBE orbitals is referred to as (EX+cRPA)@PBE.
The original (EX+cRPA)@PBE and (EX+cRPA)@HF schemes both exhibit systematic underbinding for a large variety of systems, including covalent molecules [3], weakly bonded molecules [7,8], solids [11], and molecules adsorbed on surfaces [13][14][15]. Several attempts have been made to improve the accuracy of EX+cRPA. The earliest is the so-called RPA+ scheme [23] where a local correction at the LDA/GGA level is added to cRPA. More recent attempts add second-order screened exchange (SO-SEX) [9,24]) to make the entire approach self-correlation free, or invoke cRPA in a range-separated framework where only the long-range part of cRPA is incorporated [7,8]. Among these, RPA+ improves total correlation energies considerably [25], but not binding energies [3]. The SOSEX correction performs well [9,24] with considerable additional numerical effort. Range-separated RPA schemes also improve upon the standard EX+cRPA scheme [7,8,16], however, at the price of introducing empirical parameters in the approach.
In this Letter, we offer a new perspective, based on MBPT, for going beyond cRPA, and show that a simple modification of the standard EX+cRPA scheme leads to a significant accuracy increase for molecular binding energies. We first illustrate our key idea us-  [26]. Panel (b): Decomposition of the (EX+cRPA)@HF ((EX+cRPA)@PBE) binding energy of Ar2 into individual contributions: EX@HF (EX@PBE) and cRPA@HF (cRPA@PBE). The difference between EX@HF and EX@PBE, and the SE@PBE term are also plotted. The vertical dashed line marks the equilibrium distance. Calculations are done using FHI-aims [27,28] and Dunning's aug-cc-pV6Z basis [29]. The basis set superposition error (BSSE) is corrected here and in the following.
ing the example of Ar 2 . The (EX+cRPA)@PBE and (EX+cRPA)@HF binding energy curves for Ar 2 are plotted in Fig. 1(a). Both schemes show a significant underbinding behavior compared to the reference curve modeled by Tang and Toennies [26] based on experimental data. To gain more insight into the origin of the underbinding, the EX+cRPA binding energies are decomposed into two contributions in Fig. 1(b): the exchangeonly part and the remaining cRPA part. Inspection of the individual components reveals that E cRPA c @HF is (much) more repulsive than E cRPA c @PBE, whereas at the EX level E EX @PBE is (much) more repulsive than E EX @HF. The fact that E cRPA c @PBE is more attractive than E cRPA c @HF is easy to rationalize by inspecting the corresponding frequency-dependent polarizabilities. Extensive benchmark calculations for 1225 molecular pairs [30] show that asymptotic C 6 dispersion coefficients derived from E cRPA c @HF are systematically too small by approximately 40% [28], while this error is only ∼ 10% for E cRPA c @PBE. Adding ∆vdW corrections in an attempt to reduce the remaining error in cRPA@PBE [31] only leads to minor changes in the binding energy at the equilibrium distance. What is more striking, however, is the considerable difference in binding energies at the EX level --E HF @HF−E EX @PBE (plotted also in Fig. 1(b) (red stars)). It amounts to ∼6 meV at the equilibrium distance and is thus close to the deviation of the (EX+cRPA)@PBE binding energy from the reference value.
From the viewpoint of Rayleigh-Schrödinger perturbation theory (RSPT), E EX @HF and E EX @PBE correspond to the sum of the zeroth and first-order terms in the perturbative expansions based on HF and PBE reference state respectively [33]. The difference between E EX @HF and E EX @PBE must therefore be compensated by higher-order terms in the perturbation series since the final result should be independent of the reference state, if all terms were summed up. The next term in the series is the 2nd-order correlation energy E (2) c , to which only single and double excitation configurations contribute. Here we particularly examine the contribution of single excitations (SE) to E (2) c , which can be expressed [33] as Here ψ i and ǫ i are the single particle orbitals and orbital energies of the reference state, andf is the singleparticle HF Hamiltonian -the Fock operator. A more detailed derivation of Eq. (1) is given in the supplementary material (where we simply follow the RSPT instead of the Görling-Levy PT [32]). As a consequence of the Brillouin theorem [33], E SE c trivially vanishes for HF orbitals, but is in general non-zero for KS orbitals. The contribution of E SE c evaluated with PBE orbitals (referred to as SE@PBE) to the binding energy of Ar 2 is plotted in Fig. 1(b) (violet crosses). It amounts to 50% of the binding energy at the equilibrium distance, and is close in magnitude to the contribution from E EX @HF−E EX @PBE, and to the amount of underbinding in the original (EX+cRPA)@PBE scheme. We therefore propose a new scheme by adding E SE c to E EX+cRPA (subsequently referred to as EX+cRPA+SE). In Fig. 1(a) the resultant (EX+cRPA+SE)@PBE binding energy curve is also plotted, which improves considerably over the (EX+cRPA)@PBE results, and is in close agreement with the Tang-Toennies reference curve.
It appears that the quantitative agreement between E SE c defined in Eq. (1) and E EX @HF−E EX @PBE is a general feature. We found for a set of 50 atoms and molecules that the agreement typically ranges between 70% and 100%, suggesting that replacing E EX @PBE by E EX @HF is an effective way to account for the SE contributions. This leads to a "hybrid-RPA" scheme, whose total energy is given by as an alternative to boost the accuracy of RPA. Fig. 1(a) shows that the resultant binding energy curve is in almost perfect agreement with the reference curve. At this point, it is illustrative to take a closer look at the individual contributions to E EX @HF−E EX @PBE. In Fig. 2  binding energies into their kinetic (T s ), electrostatic (E elec , external potential energy and Hartree energy combined), and EX components (E EX x ) for Ar 2 . All three energy components behave quite differently for HF and PBE orbitals. The HF kinetic energy is purely repulsive, whereas the PBE one exhibits spurious attraction at intermediate and large distances. The HF electrostatic and exact-exchange energies, on the other hand, are purely attractive and decay to zero from below, while the corresponding PBE ones become repulsive in the intermediate range and decay to zero from above at large distances. Since the PBE orbitals are less localized than their HF counterparts all three energy components decay much slower in PBE than in HF. The overall effect is that E EX @PBE becomes significantly more repulsive than E EX @HF, resulting in the underbinding behavior of (EX+cRPA)@PBE. The more physical behavior of EX@HF than EX@PBE at the EX level provides a sound basis for the systematic improvement from (EX+cRPA)@PBE to hybrid-RPA.
Indeed, the exceptional performance of the hybrid-RPA and (EX+cRPA+SE)@PBE schemes for rare-gas dimers carries over to many other molecular systems. As a second example we show results for the N 2 molecule adsorbed on benzene (N 2 @benzene), which is an important model system for studying molecular adsorption on graphene and graphite surfaces [34]. We consider two possible configurations: N 2 placed parallel or perpendicular to the benzene plane. A successful theoretical approach for this system must be able to describe the delicate balance between electrostatic and dispersion interactions. We use FHI-aims [27,28] numeric atom-centered orbital basis (6s5p4d3f 2g for C, O, N, and 5s3p2d1f for H ) augmented with gaussian diffuse functions from aug-cc-pV5Z to achieve convergence of the binding energy to within 1 meV. The results shown in Fig. 3 are very similar to the rare-gas dimers: (EX+cRPA)@HF and (EX+cRPA)@PBE underbind significantly at the equilibrium distance, while hybrid-RPA and (EX+cRPA+SE)@PBE bring the bind-  ing energy into much closer agreement with the reference curve computed with the coupled cluster method including single, double and perturbative triple excitations (CCSD(T)) [34]. In contrast, the traditional MP2 method vastly overbinds the system. Finally we examine the performance of hybrid-RPA and (EX+cRPA+SE)@PBE for the S22 database of Jurečka et al. [35], which represents a balanced benchmark set for non-covalent interactions. The molecular dimers in this database can be divided into three groups of different bonding types: hydrogen-bonded, dispersionbonded, and mixed complexes. We note that RPA in a range-separated framework has been applied to the S22 database very recently [16]. In Fig. 4 we plot the deviation from the CCSD(T) reference values [36] for the binding energies in the S22 database [35] for four RPA-based approaches and MP2. The basis set type and quality is the same as for N 2 @benzene. A detailed error analysis is presented in Table I.
We observe that the standard (EX+cRPA)@PBE scheme systematically underbinds all complexes. (EX+cRPA)@HF performs even worse for dispersion and mixed bonding, but better for hydrogen bonding. The latter case can be explained by the fact that the better performance of EX@HF dominates over the bad  [36] for the binding energies of the S22 database [35] for RPA-based approaches as well as MP2. Positive errors correspond to overbinding and negative ones to underbinding.
performance of cRPA@HF for hydrogen bonded systems. Again hybrid-RPA and (EX+cRPA+SE)@PBE correct the underbinding behavior of the standard EX+cRPA scheme, and improve the accuracy considerably. The hybrid-RPA scheme yields a mean absolute error (MAE) of 14 meV. The performance of (EX+cRPA+SE)@PBE is very similar to hybrid-RPA for dispersion and mixed bonding, albeit somewhat worse for hydrogen bonding. However, the mean absolute percentage error for hydrogen bonding (6%) is still quite small. The accuracies achieved here compare favorably to the recently developed vdW functional (vdW-DF) [37], where the MAE for the PBE-based vdW-DF results for S22 [38] is 54 meV. We also note that for covalent molecules the accuracies in the atomization energies are improved considerably by the two schemes. For instance, the MAE of the atomization energies of the G2-I set is reduced from 10.5 kcal/mol to 6.2 kcal/mol by (EX+cRPA+SE)@PBE and 6.3 kcal/mol by hybrid-RPA.
To summarize, we have unraveled the origin of the underbinding that plagues the standard (EX+cRPA)@PBE scheme, which is mostly due to the too-repulsive nature of E EX @PBE rather than the (slightly) underestimation of the long-range dispersion force by E cRPA c @PBE. This problem can be largely solved either by replacing E EX @PBE by the self-consistent HF energy E EX @HF, or by adding a SE correction to the standard (EX+cRPA)@PBE approach. Particularly (EX+cRPA+SE)@PBE is a well-defined parameter-free scheme in which the SE term does not add any significant computational cost to the approach. In addition, the SE correction is compatible with other beyond-RPA schemes like RPA+ or SOSEX. We also like to emphasize that in both schemes the cRPA evaluated with KS orbitals is retained, which is essential for producing quantitatively correct asymptotics for vdW bonded systems. Despite its success for describing vdW and covalently bonded molecules, one obvious deficiency of the 2nd-order SE as given by Eq. (1), however, is that it is not well-behaved for systems with vanishing gaps. In such cases, we propose to "renormalize" the SE contribution via a resummation of a geometrical series of higher-order diagrams involving single excitations (in the spirit of cRPA). This leads to additional terms in the denorminator of Eq. (1) which prevent the possible divergence even when the KS gap closes. A brief derivation of this renormalized SE (RSE) scheme is presented in the supplementary material. Further details and benckmark calculations will be published elsewhere [28].

DERIVATION OF THE SINGLE EXCITATION CONTRIBUTION TO THE 2ND ORDER CORRELATION ENERGY
In this section we derive Eq. (1) that is presented in the main part of this Letter -the single excitation contribution to the 2nd-order correlation energy -from Rayleigh-Schrödinger perturbation theory (RSPT). The interacting N -electron system at hand is governed by the Hamiltonian wherev ext (r) is a local, multiplicative external potential.
In RSPT,Ĥ is partitioned into a non-interacting meanfield HamiltonianĤ 0 and an interacting perturbationĤ ′ , Herev MF is some mean-field-type single-particle potential, which can be non-local, as in the case of Hartree-Fock (HF) theory, or local, as in the case of Kohn-Sham (KS) theory. Suppose the solution of the single-particle Hamiltonian h 0 is knownĥ then the solution of the non-interacting many-body Hamiltonian H 0 followŝ The |Φ n are Slater determinants formed from N of the spin orbitals |p = |ψ p determined in Eq. (1). These Slater determinants can be distinguished according to their excitation level: the ground-state configuration |Φ 0 , singly excited configurations |Φ a i , doubly excited configurations |Φ ab ij , etc., where i, j, . . . denotes occupied orbitals and a, b, . . . unoccupied ones. Following standard perturbation theory, the single-excitation (SE) contribution to the 2nd-order correlation energy is given by where we have used the fact E To proceed, the numerator of Eq. (2) needs to be evaluated. This can be most easily done using secondquantization formulation where p, q, r, s are arbitrary spin-orbitals from Eq. (1), c † p and c q , etc. are the electron creation and annihilation operators, and pq|rs the two-electron Coulomb integrals The expectation value of the two-particle Coulomb operator between the ground-state configuration Φ 0 and the single excitation Φ a i is given by where v HF is the HF single-particle potential. The expectation value of the mean-field single-particle operatorv MF , on the other hand, is given by Combining Eqs. (2), (3), and (4), one gets where ∆v ia is the matrix element of the difference between the HF potentialv HF and the single-particle meanfield potentialv MF in question.
Observing that the ψ's are eigenstates ofĥ 0 = − 1 2 ∇ 2 + v ext + v MF , and hence all non-diagonal elements ψ i |ĥ 0 |ψ a are zero, one can alternatively express Eq. (5) as wheref is the single-particle HF Hamiltonian, or simply Fock operator. Thus Eq. (1) in the main paper is derived.  For the HF reference state, i.e., whenv MF =v HF , the ψ's are eigenstates of the Fock operator, and hence Eq. (5) (or (6)) is zero. For any other reference state, e.g., the KS reference state, the ψ's are no longer eigenstates of the Fock operator, and Eq. (5) is in general not zero. This gives rise to a finite SE contribution to the second-order correlation energy. In the following we assumev MF =v KS .

RENORMALIZED SINGLE EXCITATION (RSE)
CORRECTION SCHEME The SE contribution as given by Eq. (5) corresponds to a correction to the correlation energy at 2nd-order, which can be represented by the Goldstone diagram [1] shown in Fig. 1(a). Similar to 2nd-order Møller-Plesset perturbation theory, the 2nd-order SE is ill-behaved when the single-particle gap closes. An important lesson from diagrammatic perturbation theory to deal with problems of this kind is to resum higher order diagrams of the same type to infinite order. The random phase approximation itself is a good example of this. In case of SE, one can also identify a series of higher order diagrams up to infinite order, as illustrated in terms of Goldstone diagrams in Fig. 1. Here we only pick the "diagonal" terms where a particle state a or a hole state i is, in the intermediate process, scattered into the same state by the perturbative potential ∆v =v HF −v KS . These terms dominate over the "non-diagonal" ones [2] and facilitate a simple mathematical treatment, because they form a geometrical series that can easily be summed up. Here we call the summation of this series renormalized single excitations (RSE). Following the textbook rule [3] for evaluating the Goldstone diagrams, one obtains The additonal term ∆v ii − ∆v aa in the denominater of Eq. (7), which appears to be negative definite in practical calculations, ensures that E RSE c is finite and prevents possible divergence problems even when the single-particle KS gap closes. In this context, we note in passing that a similar partial resummation of the so-called Epstein-Nesbet series [3] of diagrams to "renormalize" the 2ndorder correlation energy has been performed by Jiang and Engel [4] in the KS many-body perturbation theory.