Assessment of correlation energies based on the random-phase approximation

The random-phase approximation to the ground state correlation energy (RPA) in combination with exact exchange (EX) has brought Kohn-Sham (KS) density functional theory one step closer towards a universal,"general purpose first principles method". In an effort to systematically assess the influence of several correlation energy contributions beyond RPA, this work presents dissociation energies of small molecules and solids, activation energies for hydrogen transfer and non-hydrogen transfer reactions, as well as reaction energies for a number of common test sets. We benchmark EX+RPA and several flavors of energy functionals going beyond it: second-order screened exchange (SOSEX), single excitation (SE) corrections, renormalized single excitation (rSE) corrections, as well as their combinations. Both the single excitation correction as well as the SOSEX contribution to the correlation energy significantly improve upon the notorious tendency of EX+RPA to underbind. Surprisingly, activation energies obtained using EX+RPA based on a KS reference alone are remarkably accurate. RPA+SOSEX+rSE provides an equal level of accuracy for reaction as well as activation energies and overall gives the most balanced performance, which makes it applicable to a wide range of systems and chemical reactions.


I. INTRODUCTION
In the context of first principles electronic structure theory "exact-exchange plus correlation in the randomphase approximation (EX+RPA)" 1,2 has recently generated renewed and widespread interest.  In practice, the RPA calculations are most often carried out in a nonself-consistent manner where the exchange-correlation (xc) energy contributions are evaluated with input orbitals corresponding to an approximate, usually semilocal xc energy functional. The great interest in EX+RPA is largely due to its three attractive features: (i) The exact-exchange energy (EX) cancels the spurious selfinteraction error present in the Hartree energy, (ii) the RPA correlation energy is fully non-local and includes long-range van der Waals (vdW) interactions automatically and seamlessly, 31 and (iii) EX+RPA is applicable to small-gap or metallic systems by summing up the sequence of "ring" diagrams to infinite order. The latter is in contrast to order-by-order perturbation theories [e.g. 2nd-order Møller-Plesset (MP2) 32 ] which break down for systems with zero gap. Moreover one can interpret the RPA as an approach that screens the nonlocal exchange resulting in a frequency dependent nonlocal screened exchange interaction, as opposed to conventional or global hybrid functionals where the parameters that reduce or "screen" the exact-exchange contribution are fixed and system independent. [33][34][35] Such a system independent "screening" is expected to be unreliable for metals or wide gap insulators, where non-local exchange is almost entirely screened (metals) or prevails to a large extent (insulators).
While a critical assessment of EX+RPA is emerging 5-10,12-26,28 some shortcomings have been known for a while. Total energies are typically significantly overestimated, 3,9,11,28,36,37 which is caused by an overestimation of the correlation energy at the short range. Binding energies, on the other hand, show a tendency to be underestimated. 3,7,8,12,13,16,17,25,38 Moreover, the RPA correlation energy is not self-correlation free. 9,37,39 It has been demonstrated that the overestimation of the absolute correlation energy can be almost entirely removed by adding a second-order screened exchange (SOSEX) term. 36,37,39 For one-electron systems the self-interaction error in EX+RPA is exactly canceled by adding this term, 37,39 however, for systems with more than one electron, a many-electron self-interaction error 40,41 prevails. 39 The SOSEX can also be interpreted as a correction to the RPA correlation energy that can be included to approximately restore the antisymmetry of the many electron description. 39 Furthermore, SOSEX improves binding energies, although a sizeable underestimation persists. 9,36,37,39,42,43 The underbinding problem can also be alleviated, in particular for weakly interacting systems, by adding a correction deriving from single excitations (SE) 25 to EX+RPA built on a reference state obtained from Kohn-Sham (KS) density-functional theory (DFT). This suggests that RPA(+SOSEX) yields good estimations for the correlation energy but errors in the exchange energy are sizeable if Kohn-Sham orbitals are used to evaluate the exact exchange.
In light of these observations it is timely to extend the critical assessment of EX+RPA to a wider class of systems and to consider combinations of the corrections suggested before. In this paper we will address this objective by performing benchmark calculations for atomization energies on an appreciable test set of archetypal insulating solids and small molecules, [44][45][46][47] as well as reaction and activation energies for hydrogen and nonhydrogen transfer reactions. 48,49 The schemes we include are EX+RPA based on KS-DFT reference states, and those beyond EX+RPA by adding corrections from SE or SOSEX individually, or both of them. In addition we will also assess the hybrid-type schemes 25 where one replaces the total energy at the EX level evaluated with KS-DFT orbitals by that evaluated with Hartree-Fock (HF) orbitals, as an effective way to approximate the SE contribution. 50 The second-order single excitation correction can diverge when the gap between occupied and virtual states closes, with detrimental effects for the description of the transition states in chemical reactions. As briefly discussed in Ref. 25, including higher-order terms in the spirit of RPA permits a resummation of the SE correction, as will be demonstrated in Section II C. This so called renormalized SE (rSE) is well behaved and is included in our benchmark tests.
The paper is organized as follows: Sections II and III briefly summarize the important aspects of the underlying theory and the computational parameters of our work. Results on molecular and solid-state atomization energies as well as reaction energies and barrier heights are presented in Section IV before we draw conclusions in Section V.

A. Basics on RPA
In order to properly position the methods applied in the present work within the formal framework of DFT, we briefly recapitulate essential equations and outline the structure of the functionals used. Currently, for total energy calculations, RPA-based functionals usually use either KS-DFT reference states, i.e. single-particle wavefunctions and eigenvalues or generalized KS (GKS) 51 reference states to compute the nonlocal EX energy as well as the nonlocal correlation energy. 3,5,7 Within this context, the total energy is defined as (1) where the terms deriving from the potential contributions in the Hamiltonian, E H , the electrostatic Hartree or Coulomb energy and E ext , the (external) electron-ion interaction depend on the local density, whereas the last two terms, EX energy E x and correlation energy E c , are nonlocal contributions. Note that the KS kinetic energy, in analogy to the exact-exchange energy, is not an explicit functional of the density, but rather of the KS orbitals. The nonlocality in E x is due to the nonlocal exchange operator acting on each (occupied) orbital φ iσ (r) associated to spin σ and its well known dependence on the (nonlocal) reduced one-particle density matrix ρ σ (r, r ′ ) = occ j φ jσ (r)φ * jσ (r ′ ) reads In contrast to the optimized effective potential (OEP) method, [52][53][54] in HF theory the exchange operator is fully nonlocal, and the action of the exchange operator on a single-particle wavefunction (i.e. orbital) depends on the value of that very orbital throughout the entire space (see Ref. 55). Note that the correlation energy E c is a functional of both occupied as well as unoccupied eigenstates and requires knowledge of the associated eigenenergies as well (see below). However, both, E x and E c are implicit functionals of the electron density n (see e.g. Ref. 56). Recent work pursuing the construction of a local RPA correlation potential are presented in Refs. 57−62. Work in this direction is of great value, since it ultimately enables calculations of a self-consistent RPA correlation energies staying rigorously within the KS-DFT picture.
The RPA correlation energy can be conveniently derived from (i) perturbation theory or (ii) from the adiabatic-connection fluctuation-dissipation (ACFD) theorem. [63][64][65] Fundamental to the formalism is the adiabatic connection between the HamiltonianĤ of an interacting many-electron system and the corresponding noninteracting KS HamiltonianĤ KS . Formally, both systems may be simultaneously described by a coupling constant dependent HamiltonianĤ(λ) with λ being the coupling-constant or the scaling factor in the electronelectron interaction, v λ = λv(r − r ′ ). The electrons move in a λ-dependent external potential v λ ext (r). Note that the ground-state density ofĤ(λ) for all λ ∈ [0, 1] is constant and equals the physical ground-state density n(r), i.e. the ground-state density of the real system. H(λ = 1) is the physical many-electron Hamiltonian with v λ=1 (r) = v ext (r), andĤ(λ = 0) is the KS Hamiltonian is the electrostatic Hartree potential and v xc (r) is the xc potential. Within ACFD, the exact KS correlation energy can be written as Here ν(r−r ′ ) = 1/|r−r ′ | is the bare Coulomb interaction kernel, and χ 0 is the KS independent-particle response function at imaginary frequencies iu, where c.c. denotes "complex conjugate" and the prefactor 2 acounts for the spin-degeneracy in closed-shell systems. In Eq. (3), χ λ is the density-density response function of the "intermediately" interacting many-electron system employing a scaled Coulomb potential ν λ . We adhere to the commonly used notation of i, j . . . being occupied, i.e. hole KS states and a, b . . . being unoccupied or virtual, particle states. In principle, a Dyson-type integral equation 66 has to be solved for χ λ , with f λ xc as the xc kernel, i.e. the functional derivative of the exchange-correlation potential with respect to the density. Within RPA, f xc = 0, i.e. using many-body terminology, 67 so-called vertex corrections are not included in the response function χ or equivalently in the screening of the Coulomb interaction. Solving Eq. (5) for χ λ with f λ xc = 0 corresponds to the diagrammatic resummation of ring graphs 36,68 to infinite order. In passing we note that, working within RPA, Eq. (5) can be rearranged to reflecting the above mentioned summation of the (screened) Coulomb interaction up to infinite order in χ 0 ν λ .
As will be seen later, Eq. (6) resembles the coupled-cluster amplitude equations where socalled particle-particle, particle-hole, and hole-hole ladder terms have been removed (see Eq. 17). Starting from Eq. (6), the λ-integral is readily done and the final expression for the RPA correlation energy reads Tr{ln(1 − χ 0 (iu)ν) + χ 0 (iu)ν} . (7) B. From coupled-cluster theory to RPA and RPA+SOSEX From a DFT purist's point of view, the previously outlined ACFD terminology for the RPA is certainly the most consistent way to classify "RPA" as a correlation energy functional to the many-electron groundstate. An alternative formulation of the RPA may be motivated starting from many-body theory. Many-body or equivalently field-theoretical diagrammatic techniques originally developed in quantum electrodynamics and nuclear physics 69 have been applied to the homogeneous electron gas as well as finite systems like atoms and molecules for several decades already. For systems that are not strongly correlated, the most successful diagrammatic, partial summation technique (see Refs. 56 and 70) is the coupled cluster (CC) expansion of the manyelectron wavefunction. The CC expansion to the homogeneous electron gas has been applied by Freeman, 36 Kümmel, Lührmann, and Zabolitzky, 71 as well as Bishop and Lührmann. 72,73 The same CC expansion techniques are indispensable ingredients for highly accurate molecular calculations. Here, pioneers have beenČížek, 74,75 Paldus et al., 76 and Bartlett and Purvis 77 to name a few. A more complete list of references may be found in the recent review article by Bartlett and Musia l. 70 The CC expansion relies on the ansatz for the manyelectron wavefunction, |Ψ , to generate the exact ground state from the ground state |Φ of the reference system commonly within the HF approximation. Note thatT may be represented by a sum of single, double, and higher-order excitation operators, generating in a similar way to configuration interaction (CI) techniques, singly, doubly substituted determinants based on the HF reference wavefunction |Φ . However, the CC expansion is distinct from CI by virtue of the exponential ansatz used in CC expansions (Eq. 8) for the wavefunction |Ψ , with introducing so-called disconnected products of excitations responsible for the size-extensivity of the coupled cluster correlation energy. 78 In coupled-cluster doubles theory (CCD) the excitation operator corresponds to a double excitation operator only, whereT ≡T 2 with (10) The amplitudes t ab ij are obtained from solving a set of so-called doubles amplitude equations reading Φ ab ij |e −TĤ eT |Φ = 0.
Solving Eq. (12) self-consistently for t ab ij results in a resummation of infinitely many diagrams of a certain type. Removing all terms from the above amplitude equation that do not correspond to so-called ring-diagrams defines the so-called ring-CCD.
Recently, the equivalence between direct, i.e. "Coulomb term only" ring-CCD (drCCD) and RPA as considered by Freeman, 36 reexamined by Grüneis and Kresse 43 and Scuseria et al., 6 has been demonstrated.
Scuseria et al. algebraically showed that the CCD approximation to the many-electron wavefunction contains the ring-approximation, i.e. the RPA to the ground-state correlation energy, but also includes selected higher-order exchange and ladder diagrams. 36,72,73 In other words, RPA equals drCCD and therefore corresponds to a subset of CCD diagrams.
Within the framework of CC expansions, the RPA and RPA+SOSEX correlation energies may be calculated using drCCD amplitudes {t ab ij } by employing the respective equations, 6,36,37 The matrices B ia,jb and K ia,jb are of rank N occ × N virt , and they are defined by two-electron integrals B ia,jb = ij | ab and K ia,jb = ij | ab − ij | ba , respectively, with x={r, σ}. The amplitudes {t ab ij } are obtained from solving a set of nonlinear Riccati equations, closely related to to the time-dependent HF or more precisely the time-dependent Hartree method, The previous equation can be rewritten in a more compact form, 6 with and T ia,jb = t ab ij , underlining the quadratic order in the amplitudes' matrix T.
Freeman has evaluated the RPA correlation energy of the unpolarized electron gas for various electron densities 36 using the drCCD equations and compared them to Hedin's RPA results (see Table II in Ref. 79) following an approach suggested by Nozières and Pines. 80 Both agree to within the numerical accuracy employed in the calculations. Moreover, Freeman has gone beyond RPA via inclusion of the second-order screened exchange (SOSEX) diagram. He found that SOSEX reduces the correlation energy by about 30%. Monkhorst and Oddershede came to similar conclusions employing RPA and RPA+SOSEX to metallic hydrogen, 81 and Grüneis observed a similar reduction of the correlation energy for small atoms 37 finding good agreement with highly accurate coupled cluster correlation energies only after inclusion of SOSEX. Finally we note that until recently the formulation of SOSEX within an ACFD framework has not been entirely clear, but has lately been shown by Jansen et al.. 82

C. Single excitations and their renormalization
As alluded to above, in most practical calculations, RPA and SOSEX correlation energies are evaluated using KS orbitals from local or semilocal density functionals, 3,12 or generalized KS orbitals 7,9,17 from range-separated density functionals. This way, both RPA and SOSEX can be interpreted in terms of many-body perturbation theory (MBPT) based on a (generalized) KS reference state, where only a selected type of diagrams are summed up to infinite order. If one performs a simple Rayleigh-Schrödinger perturbation theory (RSPT) starting from an (approximate) KS-DFT reference, and examines the perturbation series at second order, one can identify a term arising from single excitations (SE), that is not included in RPA or SOSEX correlation energies. In terms of single-particle orbitals, this term can be expressed as where v HF is the self-consistent HF potential, and v eff is the effective single-particle potential that defines the noninteracting reference Hamiltonian h eff giving rise to the single-particle orbitals |i and |a in the above expression.
(See the supplemental material of Ref. 25 for a detailed derivation.) As is obvious from Eq. (18), E SE c trivially vanishes for the HF reference, i.e., when v eff = v HF , but is nonzero otherwise. It has been shown that adding this term to RPA improves the description of weak interactions significantly. 25 Note that the choice of v eff in Eq. (18) is slightly different in RSPT from that in the 2nd-order Görling-Levy perturbation theory (GL2). 83 In the latter case, v eff = v EXX-OEP , with v EXX-OEP being the exact-exchange OEP 52-54 potential. The difference of the two perturbation theories lies in the choice of the adiabatic-connection path (λ-integral) -in GL2 the electron density is kept fixed along the path way and the perturbative Hamiltonian has a non-linear dependence on λ, whereas in RSPT the λ-dependence of the perturbative Hamiltonian is linear, but the electron density varies along the λ-integral. Eq. (18) in RSPT is more efficient and practically useful in the sense that there is no need to solve the computationally intensive and sometimes numerically problematic EXX-OEP equation and more flexible in the sense that it can be matched to any suitable reference state. The price one has to pay is that the theory, strictly speaking, is not KS-DFT formulated within the ACFD framework.
The SE contribution at second order as given by Eq. (18) may become ill-behaved when the single-particle gap closes. To deal with this problem, in Ref. 25 a sequence of higher-order terms involving SE processes have been identified and summed up in the spirit of RPA. This leads to a "renormalized" SE (rSE) contribution to the correlation energy, where ∆v = v HF − v eff . The additional term i|∆v|i − a|∆v|a in the denominator of Eq. (19) is negative definite, and prevents the possible divergence of the ex- pression even when the KS gap closes. The rSE correction is therefore expected to have a more general applicability, while preserving the good performance of the 2nd-order SE for wide-gap molecules and insulators. In deriving Eq. (19), however, the "non-diagonal" elements in the higher-order SE diagrams have been neglected for simplicity. Such an approach lacks invariance with respect to unitary transformations (orbital rotations) within the occupied and/or unoccupied subspaces. The orbital-rotation-invariance can be restored by including also the "non-diagonal" elements. This can be achieved by first semi-diagonalizing the Fock Hamiltonian f = h eff + v HF − v eff separately within the occupied and unoccupied subspaces of h eff and utilizing the resultant (so-called semi-canonical ) orbitals and orbital energies in Eq. (18). A detailed description of this procedure will be presented in a forthcoming paper. However, we emphasize that results presented in this work are based on Eq. (19), but despite the lack of rotational invariance in the orbitals of this approach, numerical results are only very little affected.
As also demonstrated in Ref. 25, the SE contributions to the correlation energy can be effectively accounted for to a large extent by replacing the non-self-consistent HF total energy computed using KS orbitals by its selfconsistent counterpart. In this so-called hybrid-RPA scheme the RPA correlation energy is still evaluated using KS orbitals, whereas the EX term is evaluated using HF orbitals. The same strategy can be applied to "RPA+SOSEX" calculations. In this work, we will benchmark the influence of SE contributions on the performance of RPA and SOSEX both by explicitly including the (r)SE corrections and in terms of the hybrid scheme.
As outlined in Ref. 25 by Ren et al., rendering the energy functional stationary with respect to variations in the orbitals implies a zero correlation energy contribution stemming from SEs. This is well known as Brillouin's theorem. It will be demonstrated in this work, that SE effects represent a non-negligible contribution to the correlation energy and consequently affect results on thermochemistry and kinetics. In the field of quantum chemistry effects induced by SEs are known as orbitalrelaxation effects. 84,85 Besides MBPT discussed above, the SE terms are present in the CC theory as well. In this context, Scuseria and Schaefer have shown that CCD employing optimized-orbitals (see Ref. 86) gives results very close to CCSD. On the other hand, optimizing orbitals for CCSD calculations does not lead to significant improvements in the wavefunction. In other words, changes in the correlation energy induced upon inclusion of SEs may be effectively incorporated by means of a unitary transformation, i.e. rotation of the orbitals, as given in Eq. (6) of Ref. 86. We close this section by presenting Table I, which summarizes the acronyms of the various methods applied in this work. For the KS single-determinant reference wave function we use the Perdew, Burke, and Ernzerhof (PBE) 87 generalized gradient approximation (GGA). We adopt the notation introduced by Ren et al. in Ref. 17, hence "@PBE" means "evaluated using PBE orbitals and orbital energies". This particular choice of orbitals is mainly driven by the following arguments: (i) PBE contains no empirically adjusted parameters, (ii) performs slightly better than LDA (see e.g. Ref. 12), and (iii) it is computationally less expensive to calculate the orbitals using semilocal functionals instead of e.g. hybrid functionals. 17 In addition, once one restricts the input orbitals to KS orbitals, results have shown to be virtually identical to those obtained using PBE orbitals. 24,88

III. COMPUTATIONAL DETAILS
Computational results of the present work are based on calculations using (i) the Vienna ab initio simulation package vasp, 89-91 (ii) a development version of the gaussian 92 suite of programs, and (iii) FHI-aims. 93,94 All of the software packages used have the RPA and RPA+SOSEX functionals available since recently. 7,9,12,17 vasp uses periodic boundary conditions and projector augmented plane waves as a basis set, which makes it ideally suited for extended, crystalline systems. gaussian is based on local, analytic Gaussian type (GT) basis functions using open boundary conditions and the linear combination of atomic orbitals to expand the molecular orbitals. FHI-aims primarily uses numeric, atom-centered basis functions, but GT orbitals can be employed as well. In both cases, all the required integrals are evaluated numerically on an overlapping atom-centered grid. 93 The resolution-of-identity approximation is used to handle the four-centered Coulomb repulsion integrals and the KS response function (details of the implementation have been presented in Ref. 94). In this work GT orbitals are used in FHI-aims calculations to facilitate a direct comparison with gaussian and the extrapolations to the complete basis set (CBS) limit.
In this work we present statistical errors for the G2-1 set, 44-47 as well as for BH6, 95 HTBH38/04, and NHTBH38/04 sets of 38 hydrogen transfer and 38 nonhydrogen transfer barrier heights after Zhao et al. 48,49 Results for the molecular test sets use a two-point extrapolation procedure on the correlation energies to attain the complete basis set (CBS) limit. [96][97][98] The chosen ansatz is motivated by an atomic partial wave expansion of the two-particle many-body wavefunction, 97 where the E X corr are correlation energies corresponding to the cc-pVXZ basis sets. For G2-1, CBS calculations are based on Dunning's correlation-consistent cc-pVQZ and cc-pV5Z basis sets. 99,100 Note that throughout this work CBS extrapolation will be denoted by, e.g., cc-pV(Q,5)Z.
Moreover, G2-1 calculations employ the Boys-Bernardi counterpoise correction 101 to correct for basis set superposition errors (BSSE) within a particular basis set. Therefore, we emphasize that the CBS procedure uses BSSE free correlation energies. In order to avoid inaccuracies from numerical quadrature of xc energy contributions, gaussian calculations use a grid of 400 radial shells and 770 angular points in each shell to converge the KS orbitals. gaussian employs a root-mean-square convergence criterion for the density matrix in the SCF iteration of 0.1 µHartree, which implies an energy convergence no worse than at least 0.01 µHartree (gaussian keyword: SCF=tight). In FHI-aims the grid setting "tight" together with "radial multiplier=6" has been used to achieve convergence within one µHartree. Results on barrier heights in BH6, HTBH38/04, and NHTBH38/04 use a cc-pV(T,Q)Z CBS extrapolation of the correlation energies and do not employ counterpoise corrections. To test for the errors incurred, we compare with benchmark results obtained using RPA and RPA+SOSEX given in Ref. 9. The statistical errors in barrier heights deviate from the aforementioned benchmark values by at most 1 kJ/mol. Hence, the errors incurred using smaller basis sets are minute and consequently are not expected to bias the conclusions. The test set on atomization energies for crystalline solids includes 11 archetypal semiconductors and insulators. Specifically it comprises C, Si, SiC, BN, BP, AlN, AlP, LiH, LiF, LiCl, and MgO. The projector augmented wave (PAW) pseudopotentials (technical details in Tab. II) and kinetic energy cutoffs employed in the present calculations are identical to the ones used in Ref. 102. Table III summarizes the lattice constants used in "post-RPA" calculations. Moreover, we specify plane wave cutoffs for the overlap charge densities described in Refs. 38 and 102. The SOSEX correlation energy was calculated using a (3 × 3 × 3) Γ-centered k mesh, except for BN and BP due to a slower k-point convergence of the energy. For these systems a (4 × 4 × 4) mesh was used. RPA correlation energies are taken from the literature (see Ref. 38). In vasp, atoms are calculated using a supercell approach. The dimension of the supercells has been chosen as (9 × 9 × 9)Å 3 in size. To reduce the computational cost of the "RPA+SOSEX" calculations for isolated atoms, natural orbitals obtained using second order perturbation theory have been employed. As outlined in Ref. 107, natural orbitals substantially improve convergence of the correlation energy with respect to the number of virtual orbitals.
To assess the codes used in this work, we compare numerical results obtained using the "RPA" and "RPA+SOSEX" implementations of gaussian and FHIaims. Table IV shows correlation energies for the He atom obtained using the cc-pV5Z basis set. In order to avoid errors caused by numerical integration, we decided to use (restricted, i.e. spin-unpolarized) HF orbitals and eigenvalues for the calculation of RPA and RPA+SOSEX. The agreement found is close to perfect. Differences between results are within a sub-micro-Hartree error margin. In passing we mention that FHIaims employs the resolution-of-identity (RI) technique, 94 which (i) reduces the computational workload significantly and (ii), as shown in Tab. IV, does not sacrifice accuracy. For the molecular test sets, we always crosscheck the "RPA" and "RPA+SOSEX" results obtained with the gaussian suite of program and FHI-aims to make sure that the results presented in this work are not affected by the actual implementations. "SE" and "rSE" have so far only been implemented in FHI-aims and we use these results throughout.

IV. RESULTS AND DISCUSSION
Central findings of this work are summarized in Tab. V presenting binding energies in molecules (G2-1) and solids, HT activation energies or barrier heights (BH6, HTBH38) as well as NHT barrier heights (NHTBH38). Whenever results are compared to experiment or to the best theoretical estimates, we use mean error (ME) and mean unsigned error (MUE) as statistical measures to assess the accuracy of individual methods employed. Note  Table VI, are calculated from barrier heights in HTBH38 and NHTBH38, respectively.

A. Atomization energies of small molecules and solids
The notorious underbinding of (EX+RPA)@PBE in molecules and solids has already been demonstrated in many studies. 3, 9,12,13,17,37,38 Table V presents MEs and MUEs in binding (atomization) energies obtained using (EX+RPA)@PBE for insulating solids (see Sec. III) as well as for the molecules contained in the G2-1 set. Following the suggestion of Ren et al., 25 effects incurred by replacing the EX@PBE reference energy by the HF total energy have been checked for both molecules and solids. Indeed, HF+RPA@PBE improves binding energies of molecules and solids by almost 50% compared to (EX+RPA)@PBE. Fig. 1 presents mean unsigned relative errors (MURE) in molecular (full bars) as well as solid state (squared bars) binding energies. Overall differences in MUREs are rather small. For the commonly applied (EX+RPA)@PBE method, the MURE is found to be approximately 6%. Using HF at the exact exchange level reduces the MURE by more than 1%. It appears that the aforementioned improvements are less pronounced at the relative scale, and the error reduction is apparently bigger for solids than for molecules. The explicit inclusion of the SE contribution to the correlation energy "SE@PBE" obtained using Eq. (18) has been evaluated for molecules only. Adding "SE@PBE" to (EX+RPA)@PBE results in an ME of approximately −14 kJ/mol (see Tab. V) and an MUE of approximately 23 kJ/mol, clearly outperforming HF+RPA@PBE. Relative unsigned errors in G2-1 collected in Fig. 2 further corroborate the improvements of (EX+RPA+SE)@PBE over HF+RPA@PBE. Overall, these results confirm the findings presented by Ren et al. in Ref. 25. However, "Renormalization" of the SE contributions, as required for systems with a small one-electron band gap in PBE (see activation energies discussed in Sec. IV B) brings the atomization energies in the G2-1 test set back into almost perfect agreement with HF+RPA@PBE. Therefore, the good agreement with experiment for the G2-1 test set on the level of (EX+RPA+SE)@PBE is most likely to some extent fortuitous.
In summary, single excitation diagrams improve (EX+RPA)@PBE atomization energies of small molecules at virtually zero additional computational cost. However, as we will see below, this method fails when the one-electron band gaps in PBE become small. The better founded rSE does not perform equally well for atomization energies when comined with RPA. Combined with RPA+SOSEX it yields impressive atomization energies that are also in almost perfect agreement with the "hybrid variants" e.g. the (self-consistent) HF total energy together with "(RPA+SOSEX)@PBE". Overall this indicates that (EX+RPA+SOSEX+rSE)@PBE or HF+(RPA+SOSEX)@PBE are the methods of choice for atomization energies. The ability to accurately describe the topology of multidimensional potential energy surfaces spanned by the internal molecular degrees of freedom, i.e. the reaction coordinates, in the course of a chemical reaction, is central to first principles electronic structure methods. Calculating the energy difference between reactants and transition states (see Fig. 3) is a stringent test for the accuracy of density functionals. As mentioned in Sec. III, the HTBH38 and NHTBH38 test sets established by Truhlar and coworkers 48,49 will be used here to test the RPA-based functionals considered.
Our findings on barrier heights, i.e. activation energies (Fig. 3), are summarized in Table V Fig. 4. Note that legends given in Fig. 4 follow the color code used in Fig. 2. To establish a connection to Ref. 9, Tab. V also shows MEs and MUEs for the BH6 test set, 95 which has been introduced as a computationally less intensive, but statistically representative subset of HT/NHTBH38. However, we do not present a detailed discussion on BH6 here, but stress that errors in BH6 essentially follow the trends found for HT/NHTBH38. One of the main findings of this work is the astonishingly good performance of the conventional (EX+RPA)@PBE scheme for activation energies. To be more specific, (EX+RPA)@PBE performs significantly better for the transfer of hydrogen atoms than for reac-tions involving heavier atoms. For HTBH38, the ME obtained using (EX+RPA)@PBE amounts to −0.8 kJ/mol and the associated MUE amounts to 7.1 kJ/mol. These error margins are similar to those of some of the rangeseparated, generalized KS-DFT functionals like e.g. LC-ωPBE. 110 The latter performs very well for chemical reaction barriers (see also Sec. IV D). However, for (EX+RPA)@PBE, the MUE increases by more than 50% when elements heavier than H, like e.g. F or Cl, are transferred. The MUE in NHTBH38 obtained using (EX+RPA)@PBE amounts to 12.1 kJ/mol together with a rather distinct underestimation of the barriers by −10.5 kJ/mol (ME).
On a relative scale, the MURE for HT reactions obtained using (EX+RPA)@PBE (see Fig. 4) amounts to approximately 20%, but increases to a value approximately twice as large for NHT reactions [panel (b)]. Note that reaction 7 in NHTBH38 has a barrier height of only −1.42 kJ/mol. For this reaction MUREs are extraordinarily large leading to an increase, which is seven to eight times as large as the corresponding value in HT reactions. The statistics would be drastically biased by such a case being very likely compensated by significantly extending the test set. Therefore, we decided to exclude reaction number 7 from the test set when calculating the MURE in NHTBH38.
Both HF+RPA@PBE and (EX+RPA+SE)@PBE show a strong underestimation of barriers with maximal errors as large as 50 kJ/mol. As mentioned above, the reason for this behavior has been attributed to small HOMO-LUMO differences found for some of the transition state structures, which are not correctly described by the simple (EX+RPA+SE)@PBE scheme. Indeed, the renormalization of SE alleviates the problem, and the corresponding ME and MUE in HTBH38 obtained using (EX+RPA+rSE)@PBE are improved by almost 60% compared to (EX+RPA+SE)@PBE. Note that numerical results given in Tab. V nicely reflect the trend induced by incorporation of SE effects in the correlation energy contribution, i.e. it partially takes care of the lack of self-consistency in the EX@PBE energy. However, the rSE corrects for the strong underestimation of barriers seen in HF+RPA@PBE and (EX+RPA+SE)@PBE, but qualitatively reflects the same trend compared to (EX+RPA)@PBE.
The performance of (EX+RPA+SOSEX)@PBE for barrier heights has already been tested by Paier et al. for the BH6 test set. 9 This work extends the findings of Ref. 9 by discriminating HT and NHT reactions.
(EX+RPA+SOSEX)@PBE is less accurate for HT barriers than (EX+RPA)@PBE as indicated by an MUE of about 22 kJ/mol compared to 7 kJ/mol. Quantitatively, (EX+RPA+SOSEX)@PBE on average overestimates barrier heights for HTBH38 by the aforementioned 22 kJ/mol. This is in perfect agreement with the errors found for the BH6 test set. 9 On the other hand, (EX+RPA+SOSEX)@PBE performs substantially better for NHT barrier heights, where ME and MUE are found to be close to the ones obtained using (EX+RPA)@PBE. On average, (EX+RPA+SOSEX)@PBE overestimates NHT barriers by approximately 13 kJ/mol, whereas (EX+RPA)@PBE underestimates them by 11 kJ/mol.
To summarize this section, SOSEX and rSE tend to overestimate and underestimate reaction barrier heights, respectively. Thus it appears advantageous to combine the correction schemes in order to achieve a partial error compensation. This mechanism works most efficiently for HT reactions and somewhat less so for NHT reactions.
Taking the excellent performance of (EX+RPA+SOSEX+rSE)@PBE for binding energies (see previous Section) into account, this functional offers the most balanced description in terms of binding energies as well as activation energies.

C. Reaction energies in HTBH38 and NHTBH38
As shown in Fig. 3, knowing both forward (V = f ) and reverse (V = r ) reaction barrier heights, corresponding reaction energies ∆E are readily calculated using Note that 17 out of the 38 reactions contained in HTBH38 lead to a nonzero ∆E, whereas NHTBH38 comprises 13 reactions with a forward barrier different from the reverse barrier. The corresponding MEs and MUEs of the RPA-based functionals are compiled in Tab. VI, and the MUREs are depicted in Fig. 5 Similar to the trends found for atomization energies, HT reaction energies are significantly improved upon in- clusion of (SOSEX)@PBE to (EX+RPA)@PBE as reflected in the MUEs. For (EX+RPA)@PBE the MUE in HT reactions amounts to 18.2 kJ/mol and drops down to 4.6 kJ/mol employing (EX+RPA+SOSEX)@PBE.
Hence, it appears that eliminating the one-electron selfcorrelation error contained in RPA@PBE is beneficial for HT reaction energies. This is not entirely surprising, since the aforementioned error will be largest for breaking and creating covalent hydrogen bonds. For reactions involving heavier atoms, as exemplified by the reaction energies in NHTBH38, the correction due to (SO-SEX)@PBE appears to perform less favorably. This can be seen by inspection of Fig. 5 presenting MUREs in HT [panel (a)] as well as NHT reaction energies [panel (b)]. For (EX+RPA)@PBE the MUE in NHTBH38 amounts to 9.7 kJ/mol, which is rather low, whereas for NHT reaction energies obtained using (EX+RPA+SOSEX)@PBE, the MUE increases to 20.5 kJ/mol.
We now turn to a discussion of results obtained using the "hybrid variants", which employ the HF energy as the reference energy on the EX level. Specifically for (EX+RPA)@PBE, HT reaction energies are substantially improved upon replacement of EX@PBE through HF. As can be seen from Tab. VI, the MUE is reduced by approximately 6 kJ/mol, which translate into an improvement of the MURE by approximately 50%. HT reaction energies obtained using (EX+RPA+SOSEX)@PBE, which are fairly accurate, are hardly affected by changing to the corresponding hybrid scheme.
Employing HF+(RPA+SOSEX)@PBE, however, the MUE in NHT reaction energies is reduced by 5 kJ/mol. In addition, the ME amounts to −12 kJ/mol, which compares very favorably to the ME of −18 kJ/mol obtained using (EX+RPA+SOSEX)@PBE. In terms of performance, the combined scheme (EX+RPA+SOSEX+rSE)@PBE is on par with HF+(RPA+SOSEX)@PBE for both HTBH38 and NHTBH38 reaction energies. (EX+RPA+SOSEX+rSE)@PBE has two apparent favorable features: (i) it substantially improves HT reaction energies obtained using (EX+RPA)@PBE, and (ii) it performs approximately similarly well for all of the test sets investigated in this work.
In other words, the overall variation in error margins for atomization energies, barrier heights, and reaction energies is smallest for (EX+RPA+SOSEX+rSE)@PBE lending the functional robustness. Among the functionals discussed in this work, (EX+RPA)@PBE performs best for NHT reaction energies. Nevertheless, (EX+RPA+SOSEX+rSE)@PBE performs only slightly worse, but given the better HT reaction barrier heights and the significantly better reaction energies in HTBH38, (EX+RPA+SOSEX+rSE)@PBE is among the RPA-based functionals tested in this work, the functional of broadest applicability.

D. Comparing RPA to semilocal and hybrid functionals
To close the discussion on the performance of the RPA-and RPA+SOSEX-based functionals, we briefly compare molecular atomization and activation energies to results obtained using commonly applied semilocal as well as HF/DFT hybrid functionals. To render direct comparisons easier, Table VII repeats MUEs for G2-1, BH6, HTBH38, and NHTBH38 for three of the RPA-based functionals, which perform best, namely (EX+RPA)@PBE, HF+(RPA+SOSEX)@PBE, and (EX+RPA+SOSEX+rSE)@PBE. The above mentioned statistical errors are compared to PBE-GGA, BLYP-GGA 113,114 as well as the PBE0 111,115 and B3LYP 116 global hybrid functionals. In addition, we also present results obtained using the above mentioned LC-ωPBE range-separated hybrid functional. 110 LC-ωPBE mixes a fraction of EX at the long-range of the Coulomb interaction (for definitions, see Ref. 110), but uses only one parameter (0.40 bohr −1 ) for defining a universal range separation. It is remarkable that LC-ωPBE describes reaction barriers and atomization energies extremely accurately representing a landmark among hybrids for thermochemistry and kinetics. Admittedly, for extended systems admixture of EX on the long-range is detrimental and leads to e.g. strongly overestimated band gaps. 117 Returning to RPA, activation energies obtained using (EX+RPA)@PBE are de facto on par with LC-ωPBE (Tab. VII). Trends for GGA and global hybrid functionals like PBE0 or B3LYP are rather general, hence other GGA-type and global hybrid functionals perform quite similarly (for other functionals, see e.g. Ref. 112). Although, HF+(RPA+SOSEX)@PBE does not perform as well as (EX+RPA)@PBE for activation energies of non-hydrogen transfer reactions (corresponding MUE is almost twice as large), it performs certainly better than PBE and BLYP. HF+(RPA+SOSEX)@PBE is only slightly outperformed by B3LYP for the aforementioned activation energies in NHTBH38. According to this synopsis, (EX+RPA+SOSEX+rSE)@PBE certainly shows the most balanced description of molecular binding and barrier heights. It performs similarly well as hybrid functionals in terms of atomization energies, outperforms both GGA and hybrid functionals in terms of hydrogentransfer barrier heights, and performs equivalently well for non-hydrogen barrier heights as aforementioned hybrids do. Although, this work is not devoted to weak, van-der-Waals-type of interactions, it should be emphasized that all of the RPA-based functionals presented here do include them seamlessly as already mentioned in the introduction. It is well known that neither GGA nor hybrid functionals do show the correct van der Waals asymptote.

V. CONCLUSIONS
In summary, we have reported an extensive assessment of several exact-exchange based post-KS density functionals involving RPA correlation energies and beyond. Correlation energies have been assessed for solids as well as for small molecules. Specifically we calculated atomization energies of solids and molecules using (EX+RPA)@PBE, (EX+RPA+SOSEX)@PBE as well as HF+RPA@PBE and HF+(RPA+SOSEX)@PBE, where the latter approach gives binding energies improved by approximately 50% compared to "conventional" (EX+RPA)@PBE. Furthermore, we investigated the performance of individual functionals for chemical reaction barrier heights or activation energies em-ploying large and well established test sets. Generally, we found that it is rather difficult to systematically improve on (EX+RPA)@PBE reaction barrier heights, although modest improvements using (EX+RPA+SOSEX+rSE)@PBE were achieved for HT barriers. Importantly, the favorable impact of the correlation energy contribution stemming from SE effects on binding energies does not translate to reaction barriers. This has been explained by divergent correlation energy contributions in systems with small HOMO-LUMO gaps. Therefore, application of "SE" to systems with small oneelectron band gaps is not possible, but a renormalization of "SE" helps to alleviate the problem. Surprisingly, (EX+RPA)@PBE yields reaction energies of high accuracy for reactions involving non-hydrogen atoms. Good and robust performance of a novel RPA-based functional (EX+RPA+SOSEX+rSE)@PBE is a central finding of this work. It improves on binding or atomization energies compared to (EX+RPA)@PBE, improves on HT barrier heights as well as reaction energies.