Plane wave basis set correction methods for RPA correlation energies

Electronic correlation energies from the random-phase approximation converge slowly with respect to the plane wave basis set size. We study the conditions, under which a short-range local density functional can be used to account for the basis set incompleteness error. Furthermore, we propose a one-shot extrapolation scheme based on the Lindhard response function of the homogeneous electron gas. The different basis set correction methods are used to calculate equilibrium lattice constants for prototypical solids of different bonding types.


I. INTRODUCTION
The random-phase approximation (RPA) for the electronic correlation energy was developed during the beginnings of many-body perturbation theory in the 1950s. The detailed analytical studies of the homogeneous electron gas (HEG) in the RPA have contributed greatly to the understanding of many-body effects in solid state physics such as plasmons, screening and electronic correlation in general. 1 However, large scale practical applications of the RPA have been possible only in the last one or two decades due to increasing computer power.
In the 1970s, the RPA was formulated in the framework of density functional theory (DFT) 2,3 by Langreth and Perdew via the adiabatic connection formalism. 4,5 It can be categorized as a fifth rung functional (i.e. includes unoccupied orbitals) on Jacob's ladder. 6 In the early 2000s, renewed interest was drawn to the RPA following the work of Yan, Perdew, and Kurth 7 and Furche. 8 Most importantly, the RPA -unlike lower rung functionals -describes long-range effects like van der Waalsbonding and dispersion very well. Furthermore, the RPA is compatible with exact exchange and behaves well in the metallic limit, where other methods such as Møller-Plesset perturbation theory diverge. 9,10 Finally, the adiabatic connection formalism provides a natural way to construct beyond-RPA theories, see e.g. Ref. 11. The recent development has been summarized in the review of Ren et al., 12 which also includes derivations of the standard RPA expressions.
One limiting factor of the RPA, like for other methods based on many-body perturbation theory, is the slow convergence of the correlation energy with respect to the basis set size. In a canonical plane wave implementation, the RPA scales at least as N 3 PW with respect to the number of plane waves N PW . [13][14][15] On the other hand, the basis set incompleteness error decays only as 1/N PW . 9,16,17 The reason for this slow convergence is connected to the a) Author to whom correspondence should be addressed: stefan.riemelmoser@univie.ac.at electronic cusp condition, 18,19 which states that the exact many-body wave function has a kink at electron coalescence. This (integrable) UV -divergence is caused by the 1/r -singularity of the Coulomb potential, or equivalently by its 1/q 2 -high momentum behavior. The cusp-related convergence problem is not limited to plane waves, but applies also to local basis sets. 8 The most straightforward strategy to reduce the computational cost is to perform a basis set extrapolation. 9,17 Alternatively, range-separated DFT 20 circumvents the cusp condition directly by replacing the short-range part of the Coulomb interaction, which includes the singularity, by DFT. Range-separation has been applied to the RPA by Toulouse et al., 21,22 Janesko et al. [23][24][25] and Bruneval,26 who have all reported a significant reduction of computational cost.
Range-separation has possibly further advantages. The RPA has well known shortcomings in describing short-and mid-range correlation effects, such as the prediction of an unphysical bump in the dissociation curve of Be 2 . 21 Unsurprisingly, it is thus possible to improve on full-range RPA by creating range-separated hybrid functionals. The price to pay is the introduction of empirical parameters into the theory. Inverse range-separation, which aims to cure IR divergences, has been applied extensively to the exchange energy, e.g. in the popular HSE hybrid functional. 27 Range-separation schemes up to now have been based mostly on the error function. In this paper, we drop this constraint and investigate the performance of alternative long-range potentials. First, we recall the wave vector decomposition method of old 28 to study range-separated RPA based on momentum cutoffs. Second, we propose an optimized long-range potential ("squeezed Coulomb kernel") that performs an implicit basis set extrapolation.
In Sec. II, we discuss our range-separation scheme, and study exacts limits for range-separation based on wave vector decomposition. Then, we use the results for the low-density HEG to construct the squeezed Coulomb kernel. In Sec. III, we present numerical studies for the HEG and discuss the analytical representations of our local density functionals. In Sec. IV, we show test cal-arXiv:2001.08124v2 [cond-mat.mtrl-sci] 16 Jul 2023 culations for a small set of prototypical materials and compare the different basis set correction methods. Conclusions are drawn in Sec. V. The discussion is restricted to the case of non-spin-polarized electrons. Unless stated otherwise, Hartree units are used throughout the work.

II. THEORY
A. Range-separated density functional theory Since we are interested in a simple basis set correction method for post-DFT RPA calculations, we adopt the range-separation scheme of Bruneval. 26 Generally, the Coulomb kernel is decomposed in a long-range and a complementary short-range part where µ is a tunable range-separation parameter and cuts of the Coulomb potential at a cutoff radius r cut ≈ 1/µ. The long-range version of the RPA correlation energy E RPA,LR c is obtained by replacing V by V LR in the text-book expression. The short-range part is then formally defined as and is, in practice, approximated by a density functional This scheme provides a generalized adiabatic connection 20,29 between the full-range RPA, which is obtained in the limit µ → ∞ and DFT in the limit µ → 0. As in standard DFT, the transferability of the density functional is key to the success of the method. Near the full-range limit, short-range effects are essentially localized. Thus, we can expect that the short-range part is described exactly by the local density approximation (LDA). 3 A more concrete theorem will be given further below.
In the LDA, the short-range correction required to recover the full RPA correlation energy is given by where n(r) is the local electronic density and ε RPA,SR c,HEG is given by the difference between the full-range and the long-range RPA correlation energy per particle for the homogeneous electron gas ε RPA,SR c,HEG (n) = ε RPA c,HEG (n) − ε RPA,LR c,HEG (n).
Naturally, it is also possible to use more sophisticated density functionals. For example, various gradient correction schemes were discussed by Toulouse, Colonna, and Savin. 30 However, in this work we will stay at the LDA level of DFT, which was also the choice of Bruneval. 26 Nowadays, the most common way to separate the Coulomb potential in long-and short-range parts is based on the error function Besides the error function, other separation methods have been suggested. Toulouse, Colonna, and Savin 20 have used a scheme based on the "erfgau interaction" where c = (1 + 6 √ 3) 1/2 is some scaling constant introduced in order to achieve similar cutoff radii as the error function. The erfgau interaction provides a much sharper separation between long and short-range interactions. They showed that this is manifest in the DFT limit, where both E SR x and E SR c were flat as a function of the range-separation parameter for the erfgau interaction, but not for the error function. 20 This means that the erfgau interaction does a better job at separating out long-range effects.

B. Wave vector decomposition
When working with a plane wave basis, it is natural to seek a range-separation scheme in Fourier space. In this section, we review the available literature on what is commonly referred to as "wave vector decomposition". We comment first on the limit of small wave vectors corresponding to the DFT limit (which is less relevant to the present work but included here for completeness) and then on the more relevant case of large wave vectors, corresponding to the full-range limit.
The error function exhibits Gaussian decay over the Coulomb potential whereas the wave vector decomposition method is based on a hard momentum cutoff The latter scheme was used by Nozières and Pines 28 (NP) to investigate the RPA for the HEG. For the small wave vector limit, NP developed a series expansion for ε RPA,LR c in terms of the cutoff momentum where α = (4/9π) 1/3 and r s = (αk F ) −1 is the Wigner-Seitz radius. NP studied corrections to the RPA from second and higher order exchange diagrams and found that they contribute to the expansion above only at order O(Q 4 cut ). Thus, the contribution of small wave vectors to the correlation energy is exactly described by the RPA ε LR,RPA c,HEG → ε LR c,HEG for Q cut → 0.
The series expansion (10) is valid, if the relevant momentum transfers are small compared to the Thomas-Fermi wave vector k s = 4k F /π. For large Q cut , NP used second order perturbation theory for ε SR c,HEG (Q cut ) and then interpolated between this and the small Q cut expansion (10). This gave the NP interpolation formula for the electronic correlation energy ε NP c,HEG (r s ) = (0.031 ln(r s ) − 0.115) Rydberg.
A similar interpolation method was used by Langreth and Perdew 4,5 to describe the exchange-correlation energies of metallic surfaces. However, they used local density functional theory for large momentum transfers. Hence, their method closely resembles our own wave vector decomposition scheme.
C. Approaching the full-range limit The short-range RPA correlation energy per particle for the HEG, ε RPA,SR c,HEG , is an important quantity in our range-separation scheme [see Eq. (4)]. One requires to know ε RPA,SR c,HEG as a function of the range-separation parameter Q cut . As Q cut → ∞, V LR approaches the full Coulomb potential and E RPA,SR c vanishes. In Sec. II D, we will show that the large Q cut limit yields important information on the plane wave basis set incompleteness error. Furthermore, knowledge of the exact limits is useful for finding practical representations of ε RPA,SR c,HEG . Paziani et al. 31 have conducted studies along these lines for the error function, here we concentrate on the hard cutoff.
Numerically, ε RPA,SR c,HEG can be evaluated with relative ease, since the Lindhard response function for the HEG in terms of imaginary frequencies χ 0,HEG (q, iω) is well known [see for example Eq. (4.8) in Ref. 32]. For a given long-range potential and value of r s , we first obtain ε RPA,LR c,HEG (r s ) via (see Appendix B) and then ε RPA,SR c,HEG (r s ) using Eq. (5). We calculate ε RPA,SR c,HEG according to this procedure and show results in Fig. 1 for two different cutoffs Q cut .
In the low density or large Q cut limit, ε RPA,SR c,HEG behaves as (see Ref. 17 To estimate the validity range of the low density expansion, we have to compare Q cut to some natural momentum scale. The relevant scale is the Fermi wave vector k F , which is manifest in Eq. (14) (k F ∝ 1/r s ). Fig.  1 shows that the low density approximation is accurate for Q cut 2k F . Furthermore, Fig. 1 shows that ε RPA,SR c,HEG is nearly scale invariant when plotted versus Q cut /2k F . For Q cut → 0, ε RPA,SR c,HEG converges towards the full range RPA value ε RPA c,HEG , and for large Q cut , the scale invariance follows from the asymptotics [Eq. (14)]. However, the near scale invariance at intermediate densities is not an obvious observation.
Since we will also apply range-separation to the exchange energy, we require expressions for the respective short-range LDA correction. Analogous to Eq. (5), the short-range exchange energy per particle for the HEG is defined as where ε x,HEG = −3/(4παr s ) is the familiar Dirac expression for the full-range exchange energy per particle and ε LR x,HEG is given by (see Appendix A) The constant term is connected to the normalization of the exchange hole, see Ref. 29. As pointed out by NP, the term linear in r s cancels the respective correlation term. This cancellation applies for the error function as well, but with linear terms ±3α/(2π)µ 2 r s . 31

D. Plane wave basis set incompleteness error
The large Q cut limit of ε RPA,SR c,HEG relates to an equation previously derived by Klimeš, Kaltak, and Kresse, 16 as we briefly discuss in this section. To show this, one inserts the leading order term in Eq. (14) into Eq. (4) and uses r s = (3/(4πn)) 1/3 . This means we construct a local density functional approximation for the short-range part, yielding By Fourier transforming the density to reciprocal space, we obtain where Ω is the system volume. This is exactly the expression that has been derived by Klimeš This shows that at sufficiently high cutoffs the plane wave basis set incompleteness error is described exactly by the short-range LDA. It is important to point out that this statement contains additional information over the simple fact that E RPA,LR c (Q cut ) approaches the fullrange value as Q cut → ∞, since it precisely predicts the leading order correction.
The key assumption in the derivation of Klimeš, Kaltak, and Kresse 16 was that for high energies the unoccupied orbitals can be approximated by plane waves. This is always guaranteed as the kinetic energy then dominates the Hamiltonian. Clearly, using the HEG as a reference system makes a similar assumption. An alternative proof of this theorem using the coupling constant formalism was given previously by Burke, Perdew, and Langreth. 33 Toulouse, Colonna, and Savin 20 have derived similar theorems for the error function.
Finally, we briefly discuss how these results can be extended beyond the RPA. In the RPA+SOSEX (second order screened exchange), 34,35 Eq. (17) is merely reduced by a factor of two for non-spin-polarized systems. This can be easily understood, as the SOSEX correction simply restores the Pauli principle in the large Q cut limit, i.e. the self-interaction between electrons of same spins is exactly cancelled. 28,36 However, the exact (rather than RPA) short-range correlation energy is not simply a function of the local density alone, but also involves the on-top pair density. 20,33 Although the LDA does not generally describe the latter exactly, it is still an accurate approximation, if the ground state wave function in the weak coupling limit is well described by a single Slater determinant. 37

E. Plane wave basis set extrapolation
Prior to the work of Klimeš et al., 16 Harl and Kresse 9 provided numerical evidence that the plane wave basis set incompleteness error falls off as 1/Q 3 cut and suggested a basis set extrapolation method based on higher order corrections to Eq. (14) and suggested the extrapolation where A 5 and A 7 are further fit parameters. In Appendix B, we provide an analytical evaluation of A 3 -A 7 for the HEG, yielding where n is the electronic density of the HEG [note that our coefficients differ from the ones given in Ref. 17, see discussion after Eq. (B3)]. The terms proportional to n 2 stem from the direct MP2 diagram, compare Fig. 2.
The terms proportional to n 3 stem from the third order ring diagram and so forth. As was pointed out by Gulans, these terms can be interpreted neatly in terms of the electronic cusp. They represent probabilities that there are two, three, ..., electrons at the same place. This interpretation suggests a beyond-LDA short-range functional based on an expansion involving the pair density. However, we do not attempt to construct such a functional here and leave it up to future work.

F. One-shot extrapolation method
The basis set extrapolation schemes described above require repeated evaluations of RPA correlation energies for a set of cutoff energies Q cut . Since this is typically one of the bottlenecks in modern RPA implementations, the basis set extrapolation increases the computational cost. We now propose a one-shot extrapolation scheme that accounts for basis set incompleteness a priori via an optimized long-range potential. The Coulomb interaction is enhanced at intermediate momentum transfers q ≈ Q cut in order to make up for the cutoff at q = Q cut + ∆Q We base the construction of this "squeezed Coulomb kernel" (SCK) on the fact that at large momentum transfers, the RPA is dominated by the direct MP2 term (see Fig. 2) where χ 0,HEG (q) is the frequency integrated Lindhard response function, see Appendix B. In the same limit, χ 0,HEG (q) falls off as 1/q, which yields the leading order term in Eq. (21) The SCK is constructed to make up for the loss of correlation energy by enhancing the Coulomb potential at intermediate q assuming Details on the choice of f SCK (q) and ∆Q will be given in section III, where we also show that ∆Q plays only the role of a window parameter. This means that the effective cutoff momentum is located at q = Q cut , i.e. at the center of the window, rather than at q = Q cut − ∆Q. This approach is in spirit similar to standard ionpseudo-potential methods, which have been employed to handle the nuclear cusp. 39 The standard pseudopotentials are constructed such that they reproduce the electronic properties of an atomic reference. In our case, it is the low-density HEG that plays the role of the reference system.
We do not combine the SCK with a short-range LDA correction, because the SCK already describes low densities well enough. This is exactly the region where the LDA works best, while for high densities it transfers spurious long-range effects to inhomogeneous systems. 30,40 In fact, we argue in the following that neglecting the LDA correction can be seen as an attempted effective gradient correction. In the local interaction parameter effective gradient correction of Toulouse, Colonna, and Savin, 30 ε SR c,HEG is reduced for small r s to prevent the over-correction of correlation. This was done by choosing the range-separation parameter locally as where µ l (r) is some typical correlation length, for instance µ l (r) = αk s (r). This reduces the short-range correction for high densities and hence mimics a gradient correction, as We return now to the SCK and assume that ε RPA,SR c,HEG vanishes exactly for Q cut > µ l . We will show in section III that this is an excellent approximation for µ l = 2k F . Then, an effective gradient correction is obtained by simply neglecting the short-range LDA correction, as the assumption above implies

III. COMPUTATIONAL DETAILS
In the following we present results from various numerical studies of the HEG. We argue that for the purpose of basis set correction, one can identify µ ≈ Q cut for a fair comparison between the error function and the plane wave cutoff schemes. Since its structure is visually more revealing, we show the long-range RPA correlation energy per particle ε RPA,LR c,HEG rather than ε RPA,SR c,HEG throughout. The latter can be easily obtained by subtracting ε RPA,LR c,HEG from the full-range reference.

A. Smooth momentum cutoff
The discontinuous cutoff in Eq. (9) causes technical problems such as shell-filling effects. We follow Ref. 10 and replace the step at q = Q cut by a smooth cosine window where f cos (q) = 1 2 and ∆Q controls the "smoothness" of the window. For ∆Q → 0, Eq. (9) is recovered. Even a narrow window with ∆Q = 0.1Q cut remedies the technical issues, and changes the short-range functional very little. Most notable is a slight reduction of ε RPA,LR c,HEG at densities where r s ≈ 4/Q cut [see Fig. 3 (solid lines)]. The choice of the window parameter ∆Q represents a trade-off between smoothness and plane wave basis set convergence, since the latter is now tied to a cutoff energy (Q cut + ∆Q) 2 /2.

B. Form of the squeezed Coulomb kernel
In the construction of the SCK, we impose the following constraints: (i) it should join onto the bare Coulomb kernel at Q cut − ∆Q, (ii) it should vanish at Q cut + ∆Q, and (iii) it should be positive definite. A convenient form that satisfies the constraints (i)-(iii) as well as Eq. (25) is Similar to the smooth momentum cutoff above, ∆Q controls the width of the window [see Fig. 3 (dashed lines)]. As the function f SCK (q) contains more structure than the cosine window, we use a broader window with ∆Q = 0.2Q cut . The form of the SCK is displayed in Fig. 4. As variations in the window parameter ∆Q do not change ε RPA,LR c,HEG significantly, the effective cutoff momentum is located at the center of the window, i.e. at q = Q cut , as is the case for the cosine window as well.
C. Analytical representation of the short-range local density functional  10) and (14). For each long-range potential, there is a distinct minimum at intermediate densities, which indicates a shift from the high density to the low density regime. The relevant limit for basis set correction is that of large range-separation parameters, corresponding to the low density regime. Judging from the numerical data, we identify this regime roughly at densities given by r s 4/Q cut and r s 4/µ, which corresponds to Q cut 2k F and µ 2k F .
For the error function, ε RPA,LR c,HEG approaches the fullrange value only slowly as r s → ∞. Identifying Q cut ≈ µ, the cosine window describes low densities better than the error function, but high densities worse. This means that the cosine window separates long-and short-range effects better than the error function. The SCK is exact in the low density limit per construction, and ε RPA,SR c,HEG = The a i are determined by a non-linear least square fit for a given range-separation parameter and the constant A enforces the high density limit of Gell-Mann and Brueckner 41 ε RPA,SR c,HEG (r s ) ≈ A ln(r s ) for r s → 0 This is the appropriate behavior for any reasonable LDA functional, since (i) ε RPA,LR c,HEG vanishes for r s → 0, and (ii) the RPA becomes exact in this limit. The fit parameters a i are given in Appendix C for selected range-separation parameters.

A. Applied settings
We use the RPA implementation in the Vienna Ab Initio Simulation Package (VASP) 44 as presented in Ref. 10, which also describes the calculation of the exact exchange energy. The orbitals of the underlying DFT calculations are obtained self-consistently using the Perdew-Burke-Ernzerhof (PBE) functional. 45 The equilibrium lattice constants are found by a Murnaghan equation of state fitted to energies evaluated at seven lattice volumes centered at the experimental value, spanning a window of ±15% . If the lattice constants deviate from experiment by more than 3%, additional points are considered. In the present work we consider six prototypical materials, which are summarized in Table I. Notably, we use a large plane wave cutoff energy E max for the orbital basis set (ENCUT in VASP). This ensures sufficiently converged lattice constants for all long-range potentials. The relative errors in the lattice constants are estimated to be under 0.1%. For the response func- tion χ 0 , we use a smaller cutoff energy E χ max = 2/3 E max (ENCUTGW in VASP). Fig. 6 compares the convergence behavior of the longrange RPA energies with respect to E χ max for C and Pd at the experimental lattice constant, with E max fixed to the values stated in Table I. For both the cosine window and the SCK, E RPA,LR is fully converged with respect to E χ max once the potentials vanish, i.e. beyond the energies (Q cut +∆Q) 2 /2. These cutoff energies are given in Table  II for selected range-separation parameters.

B. Convergence behavior for total energies
This makes it possible to automatically adapt the range-separation parameter Q cut to a given plane wave cutoff E χ max or vice versa. For the error function, however, E RPA,LR c is in principle fully converged only at infinite cutoff energies, and this simple automatic adaption is not possible. In practice, one would have to perform additional convergence tests for each value of µ or resort to more elaborate schemes. 46 For the sake of simplicity, we omit these convergence tests throughout and use the large cutoffs E χ max = 2/3 E max for all values of µ, compare Table I. Fig. 6 shows that replacing E RPA c → E RPA,LR c generally underestimates the correlation energy, whereas including the short-range LDA correction (dashed curves) overcorrelates. Fig. 6 furthermore demonstrates the basis set extrapolation according to Eq. (19) (dotted curves, 1/N PW ∝ 1/Q 3 cut ), which provides full-range reference values in the following.
For C, the full-range data closely follows the 1/N PW behavior that is expected from the HEG model even at smaller values of E χ max . In contrast, the full-range curve for Pd deviates strongly from the 1/N PW behavior except for very large values of E χ max . This is due to the properties of the highly localized 4d -electrons of Pd. Thus, for our purposes C can be considered to be an "easy" material, and Pd a "difficult" one. Heuristically, this can also be related to the fact that for the same range-separation parameter the relative amount of short-range correlation is smaller for C than for Pd. Likewise, total correlation energies are better reproduced for C than for Pd.
Due to systematic error cancellation, relative energies are in general easier to converge than total energies. To obtain equilibrium lattice constants, we have to compare systems of slightly different lattice volumes. We can expect good cancellation, if the orbitals of these systems are similar. A more detailed analysis will be given in the TABLE II. Cutoff energies (Qcut + ∆Q) 2 /2 for the cosine window (∆Q = 0.1Qcut) and the SCK (∆Q = 0.2Qcut). The range-separation parameters are given in a −1 B , the associated cutoff energies are converted to eV.

C. Lattice constants
Harl, Schimka, and Kresse 10 have studied lattice constants for solids of various bonding types and found the following general behavior: LDA overbinds in comparison to the experimental values, PBE underbinds, and the RPA (with standard PAW potentials) underbinds, but less than PBE. This trend also applies to noble gas solids, though the (semi-)local functionals fail quite dramatically at describing van der Waals bonds with PBE yielding lattice constants that are significantly too large. 9 Table III confirms these results for the materials considered here. Slight differences to the full-range RPA lattice constants reported in the recent study of Klimeš et al. 47 (compare values listed as "RPA std-PAW") are due to higher cutoff energies used in the present work. Fig. 7 shows that replacing E RPA c → E RPA,LR c typically leads to larger lattice constants (dashed lines). The lattice constants tend to be overcorrected when E RPA,SR c is included (solid lines). In this section, the quality of the short-range LDA correction is only judged by whether or not it yields lattice constants that are closer to the full-range RPA base line than their respective long-range only counterparts (dashed lines), whereas in Sec. IV D, we will compare with experimental results.
The performance of the short-range LDA correction is closely linked to the performance of the standard LDA, compare Table III. The short-range correction works well for C and Si, where the difference between standard LDA and full-range RPA lattice constants is also smallest. Between those two materials, convergence with respect to the range-separation parameter is faster for Si, which is attributable to a size effect, i.e. 2k F is smaller for Si than for C. For the other materials, however, the short-range LDA correction strongly overcorrelates. The most extreme cases are MgO, Pd and Kr in combination with the cosine window. At Q cut = 2 a −1 B , i.e. (Q cut + ∆Q) 2 /2 ≈ 66 eV, most of the relevant contributions to the correlation energy are simply disregarded in the RPA. In these cases, the LDA correction results in a quite severe underestimation of the lattice constant. The upshot is that an adequate cutoff must be at least around Q cut = 3 a −1 B . Turning briefly to range-separation using the error function, we note that it is generally better behaved than the hard cosine cutoff, but shares similar issues as the cosine window. The uncorrected lattice constants are too large, whereas the corrected lattice constants are a little bit too small. But "stark" outliners as for the cosine window are missing. We relate this to the fact that the error function cuts off the Coulomb kernel very slowly so that for momentum transfers q ≈ µ a sizeable fraction of the Coulomb kernel is still present. This also means that one must use fairly large plane wave cutoffs for the response function to recover truly converged results for the correlation energy. Hence, computational gains are small when the error function is used for range-separation.
To judge the plane wave basis set extrapolation via the SCK, we compare to the non-extrapolated plane wave cutoff scheme, i.e. cosine windows using the same effective cutoff momentum Q cut . Fig. 8 shows that the SCK shares the good description for C and Si with the shortrange LDA correction, as well as the poor description for Kr and Pd. This can be attributed to the success or failure of the HEG model, which is underlying both methods: if the correlation part that is not handled explicitly corresponds to plane wave-like contributions, both methods are reliable, otherwise the errors are sizable. Hence, as for the cosine window, results for Q cut = 2 a −1 B are fairly unreliable.
As discussed earlier, we expect that the SCK mimics some aspects of a short-range effective gradient correction, if the long-range effects are sufficiently well described. In fact, the SCK underbinds C and Si, but overbinds the other materials. Noticeably, the MgO lattice constant is almost converged at Q cut = 3 a −1 B , and long-range correlation effects play an important role for MgO. However, the good description of MgO is maybe somewhat coincidental, as it again relies on fortuitous error cancellation, whereas the good description of C and Si also extends to total correlation energies (not shown here).

D. Semi-empirical RPA-LDA hybrid functional
With decreasing range-separation parameter, the results for the error function change more slowly than those for the cosine window. This is because the error function mixes long-and short-range effects and does not abruptly cut off the Coulomb kernel. As discussed in Sec. II, we expect similar results for both methods for Q cut = 2µ in the DFT limit [compare the discussion following Eq. (16)]. In contrast, we identified that both methods work similar for Q cut = µ near the full-range limit, hence the slower change for the error function. The slow change of the results with respect to µ for the error function is in fact advantageous for the construction of semi-empirical RPA-LDA hybrid functionals, where the optimal µ for one system should be transferable to other systems.
For every material considered here, there is an optimal range-separation parameter µ or Q cut that reproduces the experimental lattice constant. This is to be expected, since (i) the range-separated functional switches between the full-range RPA and LDA. (ii) the overbinding of standard LDA is generally due to the correlation part, whereas the exchange part underbinds. (iii) LDA functionals based on the RPA instead of exact (QMC) data overbind even more than standard LDA, which can be deduced from the fact that the RPA+ correction typically increases the RPA lattice constants by 0.2-0.3%. 10 Thus, by choosing the range-separation parameter empirically, we can improve the lattice constants upon fullrange RPA. From our limited data set, we observe that the error function with µ 2k F gives the best lattice constants with respect to experiment [compare Fig.  7(b)]. Kr seems off, though the experimental value used here is possibly about 0.5% too small. 43 For the cosine window, however, such a simple approximation does not work, since the results vary too rapidly as the cutoff changes.
In summary, we find that the error function is better suited for the construction of semi-empirical RPA-LDA hybrid functionals. Even though the mixing of long-and short-range effects is not ideal in terms of removing the cusp-related UV divergence, it leads to improved transferability with respect to the optimal range-separation parameter.

E. Range-separated exchange
To study the influence of short-range exchange, we now replace the short-range part of the exchange energy by LDA similar to the treatment of correlation as described in section II, The treatment of E LR x in VASP is technically elaborate, as one-center PAW terms are explicitly evaluated for The error cancellation between exchange and correlation has played a central role in the success of standard LDA. This is manifest in the cancellation of the HEG terms linear in r s , as discussed earlier. The importance of joint treatment of exchange and correlation was also emphasized by Langreth and Perdew 4 for their wave vector decomposition method. Joint treatment of rangeseparated exchange has also been employed by most of the authors investigating range-separated RPA based on the error function, 21-25 though none of these studies have commented on the role of short-range exchange on its own. Fig. 9 shows that entirely neglecting short-range exchange deteriorates the results even for large µ. However, the short-range corrected lattice constants differ only slightly from the full-range exchange counterparts, compare Fig. 7. Thus, we can conclude that the LDA describes short-range exchange very well, in fact much better than short-range correlation.
In the following, we will comment on the observation that the short-range correction works better for exchange than for correlation near the full-range limit. This finding is not entirely surprising, since exchange is essentially a long-range effect. The exchange hole is quadratic in r 12 for small electronic distances r 12 ("no cusp for exchange"). 20 This is drastically manifest in the wave vector analysis of the HEG, where only momentum transfers q ≤ 2k F contribute [compare Eq. (16)]. For the error function, the separation between long-and shortrange effects is not as extreme. As shown in Appendix A, there is always at least a small short-range contribution [compare Eq. (A3)]. For large µ the short-range LDA exchange correction still becomes exact, as was shown originally by Gill, Adamson, and Pople, 51 This is analogous to our theorem (17) on short-range RPA correlation, see also Appendix A. But even beyond that, it seems plausible that the short-range LDA correction for exchange should work well near the full-range limit. There is no cusp for exchange, and so the shortrange exchange contribution can be plane wave-like. In contrast, for correlation a large amount of plane waves is always needed to properly resolve the cusp. It may still seem non-intuitive that the joint treatment of local exchange and local correlation is not beneficial here. However, we have focused on large µ, i.e. on shortrange effects, whereas the cancellation is valid for small µ, compare Eqs. (10) and (16). Most clearly perhaps, this is seen in a real space picture. The cancellation of the HEG exchange and correlation holes pertains to the long-range parts, making the combined exchangecorrelation hole more local. As for the short-range parts, we reiterate that there is a cusp for correlation but not for exchange.
Finally, the SCK method, which is designed to remove the electronic cusp explicitly, is not applicable to exchange. However, there exist similar methods for inverse range-separation, such as the cutoff scheme of Rozzi et al. 52

V. CONCLUSION AND OUTLOOK
We have shed new light on the plane wave basis set incompleteness error in the RPA by relating it to rangeseparated density functional theory. In a gist, the plane wave basis set incompleteness error of the RPA correlation energies can be described by a complementary shortrange local density functional. In the limit of large cutoff energies, this partitioning becomes even exact. Furthermore, we have introduced a one-shot basis set correction method based on an optimized long-range potential. The method modifies the Coulomb kernel at intermediate wave vectors such that the truncated and "squeezed" Coulomb kernel reproduces the RPA correlation energy of the low-density homogeneous electron gas exactly.
In practice, we find that the success of the two plane wave basis set correction methods used here is tied to that of the underlying HEG model. If the states that are removed from an explicit correlation treatment are plane wave-like, then the correction methods are successful. If the removed states are not plane wave-like, the correction methods tend to fail or are at least less accurate. Difficult materials are obviously transition metals (here Pd), oxides (here MgO), but also in van der Waals bonded solids (Kr) it is not a simple matter to develop accurate methods for correcting basis set errors.
An interesting aspect that emerged from the present study is that since the RPA tends to underbind (underestimates relative binding energies, overestimates lattice constants) and the local density approximation overbinds (overestimates relative binding energies, underestimates lattice constants), there is always one specific range-separation parameter that reproduces the experiment. Although the few data we have inspected suggest that µ 2k F yields the best results, also a fixed µ 3 a −1 B seems to work remarkably well for lattice constants. More extensive tests will be required in order to tell whether such an approach will work equally well for other properties.
Somewhat surprisingly, we have found that the joint treatment of local exchange and local correlation does not improve the predicted lattice constants. As the short-range exchange correction is essentially exact for large µ, the errors introduced in the correlation cannot be compensated by errors in the exchange. This can be interpreted by the fact that there is no cusp for exchange.
Finally, we note that the extrapolation via the squeezed Coulomb kernel can be further improved by including higher order terms in the low density expansion. Moreover, it would be desirable to generalize the method towards non-locality in order to describe more strongly correlated systems better. Climbing this Jacob's ladder for the optimized long-range potential should lead to faster convergence with respect to the plane wave basis set size. where x = µ/2k F . For large µ, this can be expanded as By inserting this into equation (34), we obtain an expression for the LDA approximation to the short-range exchange energy Gill, Adamson, and Pople 51 have shown that, up to leading order, this is also the expression for the exact shortrange exchange, compare Eq. (35).
Appendix B: Short-range RPA energy for the HEG: Approaching the full-range limit Our derivation of the high Q cut -limit of ε RPA,SR c,HEG (Q cut ) follows the study of Gulans. 17 The text-book equation for full-range expression reads 12 where χ 0,HEG (q, iω) is the Lindhard polarizibility for the HEG in terms of imaginary frequencies. The long-range version is obtained by replacing V → V LR (Q cut ) [see Eq. (9)], and ε RPA,SR c,HEG (Q cut ) is given by the difference between those. It is convenient to express χ 0,HEG in the closed form (see Ref. 54, chap. 4) where Ψ(z) is defined as Note that Gulans missed a factor of 2 in his expression for χ 0,HEG , that accounts for summation over spin. This also affects his result for the GW basis set incompleteness error [compare Eq. (28) in Ref. 16]. We introduce the dimensionless variables y = q/2k F and t = 2ω/q 2 and use the series expansion Ψ(z) Hence, for large momentum transfers q we obtain χ 0,HEG (y, t) = − k F 2π 2 y 2 2 3 In this limit, we may also use the Mercator expansion to decompose the RPA into its ring diagram components, compare Fig. 2. We proceed to evaluate the frequency integrations, starting with the direct MP2 part (B7) After momentum integration, we find that the second order contribution to ε RPA,SR c,HEG is given by For contributions from the third order ring diagram, we evaluate ∞ 0 dω 2π and obtain for the short-range correlation energy per particle Contributions of higher order diagrams are O(1/Q 9 cut ) as well. The last equation shows that the low-density scaling behavior starts to break down at O(1/Q 7 cut ) through the influence of higher order diagrams. (C1) In Tables IV and V, the fit parameters are given for selected range-separation parameters. Figs. 10 and 11 demonstrate the fits for the cosine window and the error function respectively.