Approaching the basis-set limit of the dRPA correlation energy with explicitly correlated and Projector Augmented-wave methods

The direct random-phase approximation (dRPA) is used to calculate and compare atomization energies for the HEAT set and 10 selected molecules of the G2-1 set using both plane waves and Gaussian-type orbitals. We describe detailed procedures to obtain highly accurate and well converged results for the projector augmented-wave (PAW) method as implemented in the Vienna Ab-initio Simulation Package (VASP) as well as the explicitly correlated dRPA-F12 method as implemented in the TURBOMOLE package. The two approaches agree within chemical accuracy (1 kcal/mol) for the atomization energies of all considered molecules, both for the exact exchange as well as for the dRPA. The root mean-square deviation is 0.41 kcal/mol for the exact exchange (evaluated using density functional theory orbitals) and 0.33 kcal/mol for exact exchange plus the random-phase approximation.


I. INTRODUCTION
Density functional theory (DFT) is without question one of the most versatile approaches to approximate the many-body Schrödinger equation. 1 A great advantage is that calculations of the ground-state energy require only the calculation of the occupied manifold of the one-electron orbitals since unoccupied orbitals do not enter by construction in the Kohn-Sham (KS) energy functional. This makes the method not only fast but also makes it easy to construct suitable and compact basis sets for KS ground-state energy calculations. Similarly, reliable and transferable pseudopotentials can be constructed by choosing a small core radius and matching the scattering properties of atoms at a small set of energies. The advent of the projector augmented wave (PAW) method made the construction of potentials arguably even simpler, since the PAW method allows to reconstruct the all-electron orbitals, and hence by construction takes into account the exact shape of the wavefunctions. 2,3 Pseudopotentials, on the other hand, make approximations for the shape. Using the PAW method, constructing potentials for the entire periodic table is indeed now as simple as matching the scattering properties of atoms at a sufficiently small core radius, at a set of energies for sufficiently many angular quantum numbers. It has repeatedly been demonstrated that such PAW potentials yield excellent agreement with other all-electron methods. For instance, Paier and Kresse showed an in-depth comparison between the results for molecules, obtained from the PAW method and Gaussian-type orbitals, 4 finding agreement within 1 kcal/mol. Studies comparing many different all-electron and pseudopotential codes also confirmed that the PAW method is capable of predicting reference type quality results for solids. 5 However, present density functionals are by no means infallible. In fact, in the last decades, we have seen an everincreasing number of density functionals, addressing one or the other pitfall. It goes well beyond the scope of this introduction to attempt an even cursory review, but, the most important recent developments include machine-learned (ML) density functionals that are trained on highly accurate reference data. Attempts to construct ML density functionals have been particularly successful for small molecules where fairly large training data sets are available or can be computed with reasonable effort. [6][7][8][9][10][11] This underlines the need for highly accurate reference data, in particular for condensed matter systems, where such calculations are largely missing.
Performing reference type calculations for solids is still an extraordinarily difficult task, for the following reasons. First, correlated wavefunction methods, which are penultimately required to accurately solve the many-body Schrödinger equation, are traditionally implemented in Gaussian-type orbital codes. Now, unfortunately, many of the standard correlation consistent basis sets 12 can not be used for solid-state calculations, since these basis sets involve diffuse orbitals leading to almost singular overlap matrices that can not be inverted. Thus, reaching basis-set convergence is not as simple as for molecular systems. Second, many of the correlated wavefunction methods have adverse scaling both concerning the number of occupied states as well as the basis-set size. For instance, the coupled-cluster singles doubles with perturbative triples method 13 (CCSD(T)) scales with the 7 th power of the size of the molecular system (number of electrons and orbitals). It is hence far from trivial to perform converged allelectron calculations for solid-state systems, even if the basisset problem had been solved. Finally, a method implemented for molecules might not be readily available for calculations using periodic boundary conditions. Often the four-center Coulomb integrals become very long ranged and storing all of them might not be possible for solids. Furthermore, special attention is required in dealing with the tail of the 1/r Coulomb potential, which, even in the simplest approximation (the random phase approximation), leads to long-range dipole interactions that fall off like 1/r 6 . Hence, simple truncations, or hierarchical methods that treat long-range interactions using a lower-level theory are potentially unreliable. 14,15 Ultimately, plane waves are much simpler to use for solids and can be easily made complete by increasing the number of plane waves. But as already eluded to above, plane waves necessarily need to be used in unison with pseudopotentials or other types of projection methods, such as the PAW method. Now it might seem that one could apply pseudopotentials constructed for ground-state KS calculations also in correlated wavefunction calculations. Unfortunately, this hope turns out to be unjustified. For instance, Klimes, Kaltak and Kresse have demonstrated that the pseudized orbitals should match closely in their norm to the all-electron orbitals; an aspect usually not considered when PAW potentials are constructed. 16 If they do not match, the energy contributions related to high energy dipole fluctuations can be underestimated. Also, as for Gaussian basis set codes, it is important to describe not only the ground-state orbitals accurately but the inter-electron cusp needs to be captured as well. There is no a priori telling whether this will be the case for standard PAW potentials, although, in principle, plane waves should be sufficiently flexible to model the cusp. Therefore, validation of the PAW potentials against accurate reference calculations is very desirable.
Ideally, we would like to do this validation for a computationally "simple" method that allows one to reach accurate convergence using Gaussian basis sets and plane waves. Specifically, convergence with box size and plane-wave cutoff must be attainable. The validation method that we have chosen in the present work is the random-phase approximation (RPA) to the correlation energy. Bohm and Pines [17][18][19] developed the RPA in a series of papers between 1951 and 1953. In 1957 Gell-Mann and Brueckner 20 pointed out that the RPA correlation energy corresponds to a partial summation of the many-body perturbation theory expansion. Combining the adiabatic connection framework [21][22][23] with the fluctuationdissipation theorem 24 also recovers the RPA as a low-order approximation. For simplicity, we use the RPA in a postprocessing manner on top of Kohn-Sham Perdew-Burke-Ernzerhof 25 (PBE) orbitals. This not only simplifies the setup but the PBE functional is also consistently implemented in most codes, including the codes in the present study: VASP and TURBOMOLE.
Although the RPA is not an exceedingly accurate theory for the correlation energy it exposes many of the problems also present in more involved methods. For instance, Scuseria et al 26 established a connection between the RPA and coupled-cluster theory (direct ring coupled-cluster doubles theory). Furthermore, van-der-Waals-(vdW)-like interactions between fluctuating dynamic dipoles are fully accounted for by the RPA. These kinds of high energy fluctuations are potentially difficult to capture using pseudopotentials. 16 From a computational point of view, we note that the RPA is implemented in many different plane wave, atomic orbital, and Gaussian-type orbital codes: Furche 27 was the first to apply the RPA to molecules in a post-Kohn-Sham context. This marks the beginning of a renaissance of the RPA in quantum chemistry with a focus on improved lower scaling RPA algorithms. 28,29 In the following years further studies, applying the RPA to different molecular systems including the HEAT 30 and G2-I AE 31,32 sets, were published using plane waves, [33][34][35][36] Gaussian-type orbitals, 27,28,37 or numeric atomcentered-orbital basis sets. 38,39 A somewhat troublesome observation is that the atomization energies reported in these studies markedly disagree even for simple diatomics: for example, the N 2 RPA atomization energy was reported to be 224 , 223 and 220 kcal/mol by Olsen and Thygesen, 35 Furche and Van Voorhis 40 and Ren et al. 39 , respectively. This calls for a careful evaluation and publication of reference type results that can help to validate other codes and other pseudopotentials. Our ultimate goal is, of course, to use the validated PAW potentials for solid-state calculations and to apply them with methods beyond the RPA. [41][42][43][44][45][46][47] This implies that we seek to reach chemical precision, that is, at least 1 kcal/mol. We show conclusively that this goal can be reached.
The remainder of the paper is organized as follows. In Sec. II the adiabatic-connection fluctuation-dissipation theory leading to the random-phase approximation is presented. The computational details are discussed in Sec. III followed by the results in Sec. IV. Finally, Sec. V is devoted to the conclusion.

II. THEORY
Hartree atomic unitsh = m e = e = 1 are used throughout this paper. The adiabatic-connection formalism 21-23 links the non-interacting Hamiltonian adiabatically to the fully interacting Hamiltonian by introducing a coupling-strength parameter λ that scales the electron-electron interaction at fixed electronic density n(r). The exchange-correlation energy functional is then formally given as where the exchange-correlation hole n λ xc (r, r ′ ) is linked to the response properties of the system by invoking the fluctuationdissipation theorem 24 The linear response function χ λ of the λ -interacting system fulfills a Dyson-type equation 48,49 χ λ (r, r ′ , iω) =χ 0 (r, r ′ , iω) + dr 1 dr 2 χ 0 (r, r 1 , iω) Here f xc denotes the exchange-correlation kernel at coupling strength λ and χ 0 (r, r ′ , iω) is the independent particle response function of the KS system with the KS orbitals φ i , their energies ε i and occupation number f i . Together, Eqs.
(1)-(4) define the exchange-correlation energy within the adiabatic-connection fluctuation-dissipation theory (ACFDT). These equations are exact if the KS groundstate density agrees with the true many-body ground-state density, however, the evaluation of the response function of the interacting system (Eq. (3)) requires some approximation for the exchange-correlation kernel. The simplest one-the RPA-is obtained by setting f xc to zero. This allows to write the exchange-correlation energy as E RPA x is the Hartree-Fock (HF) exchange energy functional evaluated with KS orbitals φ KS . Together with the kinetic energy and the Hartree term this corresponds to the HF energy functional evaluated using KS orbitals, which we will refer to as Hartree plus exact-exchange (HXX) energy throughout this paper. This name originates in DFT where usually only approximate exchange functionals are used. The direct RPA correlation energy E RPA c is given by with Tr[AB] := drdr ′ A(r, r ′ )B(r, r ′ ) and ν abbreviates the Coulomb kernel.

III. COMPUTATIONAL DETAILS
The ACFDT-RPA atomization energies are assessed for the 26 molecules included in the HEAT set 30 and for 10 additional molecules (BeH, C 2 H 4 , C 2 H 6 , H 2 CO, CH 4 , CH 3 OH, Li 2 , LiH, LiF and N 2 H 4 ) from the G2-1AE set. 31,32 The resulting test set contains molecules involving first and second row atoms (i.e., H to F). The HEAT geometries are found in a paper by Harding et al. 50 while the G2-1AE geometries have been taken from the supplementary material of Ref. 51.

A. Plane-wave basis-set calculations
In a periodic code, an isolated atom or molecule is simulated by placing it in a sufficiently large unit cell to suppress Shown are the orbitals that are treated as valence orbitals and the label of the respective pseudopotential in the VASP database. The radial cutoffs for each angular momentum quantum number are specified as n × r cut , where n specifies the number of projectors and r cut is the radial cutoff in atomic units for the specific angular quantum number. The local potential corresponds to the all-electron potential which is replaced by a local pseudopotential below the radius r core (a.u. interactions between its periodic replica. To obtain reliable atomization energies, we kept the unit cell and basis-set size consistent between the atomic and molecular calculations. Symmetry broken unit cells were used to avoid fractional occupancies. We use the Vienna Ab initio Simulation Package (VASP) that relies on the PAW method 2,3 and thus employs pseudopotentials. We found the set of PAW potentials listed in Table I to yield the most accurate correlation energies. Note that the semi-core 1s states of Li and Be were treated as valence states. The total ground-state energy in the RPA is given by the HF energy functional plus the ACFDT-RPA correlation energy. The RPA is employed in a post-processing manner and thus the necessary inputs -the KS orbitals and eigenvalues -were obtained using the GGA by Perdew, Burke, and Ernzerhof (PBE). 25 Note that the calculation of E RPA c involves the evaluation of the independent particle response function and as such requires summation over many virtual orbitals (see Eq. (4)). The evaluation of E RPA c implemented in VASP scales cubically with the system size 29 and is thus expensive compared to the evaluation of E HF [{φ KS }]. Hence, the HXX contribution to the total ACFDT energy was converged separately and subsequently added to the extrapolated RPA correlation energy. In the following, we discuss technical details concerning the HXX and ACFDT-RPA calculation. By setting LASPH =.TRUE. in all DFT calculations, we take the aspherical contributions inside the PAW sphere exactly into account for the electrostatic energy as well as the exchange and correlation energy.

HXX energy
In the evaluation of the Fock operator, the Coulomb kernel is given by and is thus ill defined at G + q = 0, where V is the volume of the unit cell, q are the crystal momentum points chosen to sample the Brillouin zone, and we denote the reciprocal lattice vectors with G. A similar divergence issue occurs for the Hartree term but by combining the Hartree term with the electron-ion and ion-ion term, a well defined expression for the total electrostatic energy can be obtained. 52 In order to deal with the G + q = 0 component in the exchange operator, the Coulomb kernel is truncated in real space 53 to a sphere with unit-cell volume V and radius R c . This yields the following modified Coulomb kernel This effectively cuts off inter-molecular interactions in the exchange energy. The HXX energy has been evaluated using 24 Å boxes and a 1000 eV energy cutoff. The resulting atomization energies are converged to 1 meV indicating that intermolecular static dipole-interactions are well suppressed by the large cell and thus no volume extrapolation is required. For a detailed discussion on cell size extrapolation for isolated molecules in VASP, we also refer to the supplementary material of Ref. 54.

ACFDT-RPA correlation energy
Similar to the exchange operator, the RPA correlation energy, Eq. (5), is also ill-defined at G = 0 due to the Coulomb kernel, Eq. (6). Unfortunately, the truncation procedure outlined above is not applicable, since it results in an oscillatory convergence of the correlation energy related to the cosine modulations in the Coulomb kernel. Hence, alternative procedures need to be found to deal with this issue. In a standard RPA calculation in VASP, all components G + q = 0 are excluded in the calculation of the trace in Eq. (5), that is all rows and columns for G + q = 0 are suppressed. This introduces an error that is proportional to the volume of the unit cell, since important long-range contributions related to the correlation energy are not included. We correct for this by calculating the Taylor expansion of the polarizability χ 0 G,G ′ (q, ω) for |q| → 0. The wings are obtained by setting |G ′ | = 0 or |G| = 0 and approaching |q| → 0 from three different Cartesian directions. The head is a 3 × 3 tensor describing the polarizability of the molecule at G ′ = 0 and G = 0 for infinitesimal small wave vectors |q|. For the three Cartesian directions, the head and wings are added back to the polarizability matrix and a wave vector dependent correlation energy (i.e., dependent on the direction ofq) is calculated. The final correlation energy is the average of the values obtained for the three Cartesian directions. This specific treatment is selected by specifying LOPTICS=.TRUE. during the RPA calculations in VASP, and it significantly improves the cell size convergence. Note that this correction has not always been used in the past for calculations on small molecules. If the correction is not applied, it is mandatory to correct for the resulting error, which is proportional to the inverse of the volume of the unit cell. 54 Even with this correction, a residual undesirable contribution that is proportional to 1/V 2 prevails. This contribution relates to the vdW interaction between the repeated atoms or molecules (see second next paragraph).
The calculation of the independent particle response function χ 0 requires a large number of virtual orbitals and their respective eigenvalues. To avoid any basis-set induced errors, all virtual orbitals spanned by the plane-wave basis set are calculated by exact diagonalization of the KS Hamiltonian, and hence as many orbitals as plane waves in the basis set are included in the calculation of χ 0 . Note that only plane waves with a kinetic energy smaller than E cut are included in the basis set. The RPA correlation energy converges slowly with respect to the dimension of the response function χ 0 governed by the energy cutoff E χ cut . Specifically, the response function includes only components up tō For sufficiently large values of E χ cut , the correlation energy shows the following functional behaviour 34 allowing to extrapolate to infinite cutoffs. This procedure is done automatically for RPA calculations in VASP based on eight different cutoff values E χ cut . The key trick is to evaluate the independent particle response function for the largest value of these cutoffs E χ cut , then progressively remove rows and columns from the matrix and calculate the correlation energy for each of these truncated response function arrays. If the largest cutoff E χ cut equals two-thirds of the plane-wave cutoff for the orbitals, E χ cut = 2/3E cut , then the automatic extrapolation is accurate to few meV. Specifically, we systematically increased the plane-wave cutoff E cut and E χ cut = 2/3E cut and found that the correlation energy contribution to the atomization energy is converged to 8 meV (0.2 kcal/mol) for all molecules using an 800 eV plane-wave cutoff. For these convergence tests, we used (8 Å) 3 unit cells for both the atoms and molecules .
As already eluded to above, the RPA correlation energy also needs to be converged with respect to the cell size. Note that, this is made easier in the present work since we only desire correlation energy differences. Fixing the energy cutoff now to 800 eV and performing RPA calculations using (8 Å) 3 , (9 Å) 3 and (10 Å) 3 cells we can estimate the errors related to the interaction between repeated images. The correlation contribution to the atomization energy ∆E c obeys a V −2 volume dependence (see Fig. 1) due to the dynamic dipole interactions between periodically repeated entities. These are taken into account by extrapolating ∆E c to the infinite cell size limit using In Fig. 1, the extrapolation of ∆E c is shown graphically for LiH and N 2 . For instance, the LiH correlation energy difference ( Fig. 1 (a)) changes significantly between the (10 Å) 3 cell (0.916 eV) and the infinite cell size limit (0.893 eV). Similar behavior is observed for the other lithium-compounds (Li 2 and LiF) as well as for BeH suggesting that, for (10 Å) 3 cells, there are still sizable long-range interactions present between the repeated molecules and atoms. On the other hand, N 2 ( Fig. 1 (b)) converges quickly with respect to the unit cell volume. Note the different scales of the y-axis in Fig. 1. Hence for N 2 , only a minor difference of 2 meV is observed between the (10 Å) 3 cell result and the infinite cell size limit. For the majority of the molecules, long-range vdW interactions are sufficiently suppressed in the (10 Å) 3 cells resulting in converged correlation energy differences. However, for systems with spatially diffuse highest molecular orbitals such as the considered Li-compounds as well as BeH an extrapolation to the infinite cell size limit is necessary. In order to be consistent, all correlation energies were extrapolated to the infinite cell size limit using (8 Å) 3 , (9 Å) 3 and (10 Å) 3 results at an 800 eV plane-wave cutoff.

B. Gaussian basis-set calculations
The Gaussian basis-set calculations were carried out using the TURBOMOLE program package. 55 The "direct" RPA approach (dRPA) 56 was applied using converged PBE orbitals. The density-functional computations were performed with regular two-electron integrals and employing TURBO-MOLE's quadrature grid 5a. 57 To ensure rapid convergence to the basis-set limit, the corresponding dRPA computations were carried out in terms of explicitly correlated direct-ringcoupled-cluster-doubles (drCCD-F12), which is equivalent to dRPA-F12 calculations. 58,59 In order to monitor the convergence to the basis-set limit, the correlation-consistent basis set series aug-cc-pCVXZ (X= D, T, Q, 5, 6) 12,60-62 were employed in a deciontracted fashion in all computations for the elements H, C, N, O, and F. Due to the fact that decontraction of the aug-cc-pCVXZ basis leads to (numerical) linear dependencies, the sp-sets were replaced by those of the corresponding primitive sp sets of the aug-cc-pVXZ 12,60,62 basis sets.
Besides the regular one-electron basis set, dRPA-F12 calculations require three auxiliary basis sets. These are a basis set for the density-fitting approximation to certain twoelectron integrals (MP2-fitting basis; CBAS), a basis set for the density-fitting approximation to the matrix representation of the exchange operator of HF theory (exchange-fitting basis; JKBAS), and a complementary auxiliary basis set (CABS) for the dRPA-F12 calculations. In order to minimize effects due to the incompleteness of the MP2-fitting basis, large MP2-fitting basis sets optimized for the augcc-pCV6Z basis were used in all computations for the elements H, C, N, O, and F. The optimization was done as described in Ref. 63, and resulted in a 14s13p12d10f9g7h6i2k set for H and 22s21p20d19f18g16h15i8k sets for C, N, O, and F. Similar considerations led to an optimized large JKBAS for the elements H, C, N, O, and F, which resulted in 13s12p5d4f1g, 24s18p13d12f12g9h5i, 27s18p15d14f10g7h8i, and 28s18p15d16f11g9h6i for H, C, N, O, and F, respectively. In the case of Li and Be, the JKBAS was taken from Ref. 64.
A restricted reference determinant was employed for closed-shell systems, all calculations for open-shell systems were performed using an unrestricted reference determinant. The HF self-consistent field (HF-SCF) energy was improved by taking into account single excitations into the CABS (CABS singles correction). The dRPA-F12 calculations were carried out using ansatz 2A [T+V] and the fixed amplitudes (sp) approximation with respect to spin-flipped geminals as proposed in Ref. 65. In all cases, a Slater-type geminal with an exponent of γ = 1.4 a −1 0 was used. All convergence criteria have been set to at least 10 −11 .

IV. RESULTS
In the following two subsections, we report the calculated non-relativistic HXX and RPA atomization energies for the 26 molecules included in the HEAT set and the 10 selected molecules of the G2-1 set. We compare plane-wave and Gaussian basis-set calculations and provide scalar relativistic results in the Supplementary material. The difference between the molecular energy ε(M) and the energies of the individual atoms ε(X) yields the atomization energy D(M) of a molecule where n X labels how many atoms of the species X are included in the molecule.
A. HXX atomization energies Table  II summarizes the HXX atomization energies obtained using plane waves (PW) and Gaussian-type orbitals (GTO). The maximal discrepancy between the two approaches is 0.85 kcal/mol, and the root-mean-square deviation amounts to 0.41 kcal/mol. We note that the plane-wave results are systematically smaller than the GTO calculations (i.e., underbinding in the PAW case). Fig. 2 visualizes the difference between the atomization energies predicted by the two methods ∆E HXX = E HXX PW − E HXX GTO . Paier et al. 4 observed a similar behavior when comparing PBE and PBE0 atomization energies for plane waves and GTOs. They suggested that these discrepancies originate either from basis-set superposition errors (BSSE) in the case of the GTO calculations or from the frozen-core approximation used in the PAW method.  30 and selected molecules of the G2-1 test set. 31,32 The energies were calculated using plane-wave (PW) and Gaussian (GTO) basis sets both without relativistic corrections. Energies are given in kcal/mol. The differences between PW and GTO are displayed in the column ∆ PW−GTO . Additionally, the root-meansquare deviation (RMSD) and mean signed deviation (MSD) are given in the final rows. In the present work, we have taken great care to remove any remaining uncertainties from the GTO results with respect to the basis-set size. Thus, the remaining uncertainties are well below 0.04 kcal/mol. To estimate the effect of longrange interactions between adjacent molecules on the planewave results we performed calculations for a series of unit cells, as explained before. We found that long-range effects are sufficiently suppressed by using (24 Å) 3 unit cells. Furthermore, the plane-wave results were also converged with respect to the energy cutoff to attain about the same accuracy of 0.04 kcal/mol. As such, we do not expect BSSE or finite-size errors to affect the reported numbers.

Molecule
Hence, we need to conclude that the remaining error is due to the PAW approximation. Specifically, the core electrons are kept frozen in the PAW calculation at the level of the PBE reference atom calculations. This is, in fact, a somewhat uncontrolled approximation, but one that is very difficult to remove in the PAW methodology. Specifically, the PAW approximation is so fast and numerically robust because the core electrons are kept frozen, which removes a huge energy contribution from the total energy. We speculate that this causes an underbinding because the core relaxation effects are likely to be larger for the lower symmetry molecules than for atoms: keeping the core rigid and frozen results in an upper bound for the total energy, possibly with a larger underestimation of the energy for molecules. We, furthermore, modified the energies at which the partial waves are constructed in the PAW method, but found no systematic way to improve the atomization energies.
The assumption that freezing the core electrons is the remaining source of the error, is also corroborated by the observation that H 2 , LiH, Li 2 as well as BeH show a remarkably good agreement between both methods. In these cases, the 1s semi-core states are treated as valence states in the plane-wave calculations (and, as matter of fact, also in the GTO calculations). The remaining discrepancies are certainly within the error bars of the present calculations.
We conclude that the majority of the deviations between plane-wave and GTO calculations at the HXX level are due to the frozen-core approximation used in the PAW method. Hence, we estimate the residual error of the used GW PAW potentials (see Table I) to be around 0.5 kcal/mol for Hartree-Fock type calculations. The agreement between both methods is certainly within the desired chemical precision (1 kcal/mol) and satisfactory for KS, Hartree-Fock, and most likely any other mean-field calculation.

B. RPA atomization energies
We obtain the plane-wave total RPA atomization energies by adding the HXX atomization energies (see Table II) and the RPA correlation energies (RPAc). Table III summarizes the non-relativistic RPA results using plane waves as well as GTOs. We refer to the supplementary material for the scalar relativistic RPA atomization energies. Comparing the plane-wave RPA atomization energies to the GTO results, we find the largest deviation of −0.8 kcal/mol for H 2 O 2 while the root-mean-square deviation amounts to 0.33 kcal/mol. The uncertainty of the GTO calculations with respect to the complete-basis-set limit is at most 0.2 kcal/mol. Regarding the plane-wave results, we estimate the errors due to the finite basis-set size and the extrapolation to infinite cell size to be around 0.2 kcal/mol as well.  30 and selected molecules of the G2-1 test set. 31,32 The energies were calculated using plane-wave (PW) and Gaussian (GTO) basis sets without relativistic corrections. Energies are given in kcal/mol. The differences between PW and GTO are displayed in the column ∆ PW−GTO . Additionally, the root-mean-square deviation (RMSD) and mean signed deviation (MSD) are given in the final rows.  Fig. 2 shows the deviation of the plane-wave atomization energies from the GTO calculations for the HXX and total RPA atomization energies as well as the RPA correlation energy (RPAc) contribution to the atomization energy. We observe that the plane-wave RPAc energies tend to over-correlate in comparison to the GTO results. Nevertheless, the reported RPAc energies are within the desired 1 kcal/mol accuracyexcept for the molecules C 2 H 4 and C 2 H 6 whose correlation energies deviate by 1.04 and 1.05 kcal/mol respectively. In contrast to the HXX atomization energies, we observe that the RPAc energies are systematically less biased.
Furthermore, Fig. 2 shows that there is some degree of error cancellation between the exchange and correlation energy contributions. One notes that the dips in the yellow line (∆E RPAc ) are mirrored by corresponding dips with an alternate sign in the blue line (∆E HXX ). Thus, whenever the exchange energy is more strongly underestimated, the correlation energy RPAc is more strongly overestimated. The error in the combined exchange-correlation energy is therefore on average reduced compared to the individual errors in HXX and RPAc. On average there is no clear tendency for over-nor underbinding in the PAW method for the final atomization energies, which is also supported by the comparatively small mean signed deviation (MSD) in Table III. Of particular interest are again those cases where we treat all electrons exactly in the PAW case: H 2 , Li 2 , LiH and BeH. It is clear that for these molecules, both the HXX and RPAc errors are small, and agreement for GTO and PAW is well within the combined estimated error bars of 0.3 kcal/mol. As to why there is some error cancellation we can only speculate: it is clear that the RPA screens the exchange interaction between the electrons, that is, in essence the bare exchange is screened by an electron cloud related to the polarizability of the electron gas. This might imply that errors in the description of the bare exchange are also partly canceled (screened) by the response of the electron gas. This compensation is, however, not exact; a residual error is clearly observed. Note that we have ordered the molecules according to their residual difference between plane waves and GTOs, but inspection of the order of the molecules does not suggest any systematics for this residual error. Possible sources are i) the fact that we keep the 1s orbital frozen at the level of the PBE atoms in the PAW case, whereas they are fully relaxed using the PBE functional for the GTO case, ii) residual shape approximations used when evaluating the PAW RPAc correlation energy, iii) prevailing technical errors such as the extrapolation to infinite unit cell volume and basis-set extrapolation errors in the plane-wave calculations, but also remaining small basisset truncation errors for the GTO case. The latter two errors might well still amount to 0.2-0.3 kcal/mol. To be certain that some of the prevailing errors are also related to i) and ii), we spot-checked one of the outliers, H 2 O 2 , very carefully. For this specific case, increasing the plane-wave energy cutoff or the box size does not significantly improve the agreement between the PAW and GTO methods. Also, different PAW potentials yield consistent results within 0.2 kcal/mol. We conclude that the VASP PAW results can hardly be improved for H 2 O 2 and that some of the remaining errors are likely related to the PAW technique, either the freezing of the core orbitals at the level of the oxygen atom in the PBE, or other small numerical inaccuracies introduced during the fairly complex construction of the PAW potentials.
To investigate the impact of including the 1s core states in the correlated calculations, we also evaluated the GTO RPA correlation energies when explicitly correlating the 1s orbitals for the C, N, O, and F atoms. The comparison of the all-electron GTO calculations to the plane-wave results (still treating the 1s electrons as core states) is provided in the Supplementary material. We find that including the core states in the GTO calculations leads to a sizable increase of the atomization energies of 0. 8 (Table II and  ual error is larger than the differences between the PAW and GTO calculations when the 1s core states are consistently uncorrelated. We conclude that some of the residual PAW errors in the HXX and RPAc energies compensate each other, resulting in a very good agreement of the total RPA atomization energies. While BSSE and finite-size errors (volume and plane-wave basis-set size errors) hardly contribute to the discrepancies observed for the HXX energies, this can not be entirely ruled out for the RPA correlation energies. The observed discrepancy might originate from an interplay of the BSSE, errors in the extrapolation to infinite volume and infinite basis-set size, as well as residual PAW errors. Nevertheless, basis-set extrapolation and unit-cell-size extrapolation are necessary in order to obtain accurate RPA atomization energies. If done carefully, both approaches agree well within the desired accuracy of 1 kcal/mol for all investigated molecules. This confirms that consistent RPA correlation energies between different basis sets can be obtained by a careful convergence study and suitable PAW choices in the plane-wave calculation.

V. CONCLUSION
The present work is motivated by the need of validation of PAW potentials for correlated wavefunction calculations. PAW potentials are usually constructed for KS ground-state calculations and their accuracy is rarely assessed for such correlated calculations. We decided to evaluate the potentials for the atomization energies of molecules, since comparison to high-quality GTO calculations is at least in principle straightforward. The random phase approximation was chosen as a "simple" but prototypical correlated method, since its favourable scaling allows to achieve accurate and technically converged results using both types of basis sets, but the method also exposes many of the problems present in more involved methods. To this end, reference type RPA atomization energies for the HEAT set as well as 10 selected molecules of the G2-1 set were calculated using plane waves and the PAW method, and GTOs. All calculations were performed on top of PBE ground-state calculations, using PBE orbitals and oneelectron energies.
After reviewing the basic theory, we discussed the technical details and strategies required to obtain highly accurate RPA atomization energies using plane waves in the PAW method as well as using GTOs. Regarding the plane-wave calculations, special care must be devoted to the treatment of the G = 0 component in the exchange and correlation energy. Furthermore, while obtaining cell and basis-set converged HXX energies is straightforward by increasing the cell size and plane-wave cutoffs, this is computationally still rather challenging for the RPA correlation energy. Convergence of the RPA correlation energy with respect to the plane-wave basis-set size was achieved by an extrapolation with respect to the dimension of the response matrix. Furthermore, longrange vdW-like interactions between repeated images were removed via a cell size extrapolation. This was found to be particularly important for systems with spatially diffuse molecular orbitals, such as the considered Li-containing molecules and BeH. Likewise, the Gaussian basis-set correlation energies were evaluated at the complete basis-set limit using an explicitly correlated dRPA-F12 method. We expect the uncertainties of both methods with respect to the complete-basis-set limit to be smaller than 0.2 kcal/mol. Using these setups, we compared the plane-wave and GTO results. The two methods are found to agree within chemical accuracy (1 kcal/mol) for the atomization energies of all considered molecules. We report root-mean-square deviations of 0.41 and 0.33 kcal/mol for the HXX and RPA energies, respectively.
For the HXX atomization energies, we observed a systematic underbinding using the plane-wave results. Since we do not expect BSSE or finite-size errors to be present in these calculations, we attribute this difference to the frozen core approximation employed in the PAW method. More specifically, the 1s core states of the atoms (except for H, Li, and Be) were kept frozen at the level of the PBE reference atom calculations and were not relaxed in the preparatory PBE mean-field calculations. This means that our plane-wave HXX energies are upper bounds, and relaxation effects that are expected to be larger for the molecules will increase the binding energies in the plane-wave calculations potentially bringing them in line with the GTO HXX results.
For the RPA correlation energy, we found a reasonably systematic trend for over-correlation, but the relative error is not as consistent among different molecules as it was for HXX. We gather that besides a residual PAW error and the just discussed froze core error, also BSSE and extrapolation errors contribute to the discrepancies between the two codes (the combined error is likely to be 0.3 kcal/mol). Nevertheless, the opposing systematic trends for the HXX and RPA correlation energies results in an error cancellation leading to the very satisfactory agreement of the final RPA atomization energies (RMSD 0.33 kcal/mol), with differences being below the threshold of chemical accuracy.
In summary, the present work reports reference-type RPA atomization energies for the HEAT set and selected molecules of the G2-1 set for GTOs and plane waves. Our assessment convincingly demonstrates that the employed PAW potentials are applicable to correlated wavefunction calculations and that the residual error due to the PAW approximation is below the threshold of chemical accuracy. Going forward, this validates the use of the PAW potentials in solid-state correlated wavefunction calculations where reference calculations are not yet easily and systematically doable.

SUPPLEMENTARY MATERIAL
See the supplementary material for scalar relativistic HXX and RPA atomization energies using plane-wave and Gaussian basis sets as well as the comparison of all-electron GTO RPA atomization energies to the plane-wave results.