Efficient Treatment of Correlation Energies at the Basis-Set Limit by Monte Carlo Summation of Continuum States

The calculation of electron correlation is vital for the description of atomistic phenomena in physics, chemistry, and biology. However, accurate wavefunction-based methods exhibit steep scaling and often sluggish convergence with respect to the basis set at hand. Because of their delocalization and ease of extrapolation to the basis-set limit, plane waves would be ideally suited for the calculation of basis-set limit correlation energies. However, the routine use of correlated wavefunction approaches in a plane-wave basis set is hampered by prohibitive scaling due to a large number of virtual continuum states and has not been feasible for all but the smallest systems, even if substantial computational resources are available and methods with comparably beneficial scaling, such as the Møller–Plesset perturbation theory to second order (MP2), are used. Here, we introduce a stochastic sampling of the MP2 integrand based on Monte Carlo summation over continuum orbitals, which allows for speedups of up to a factor of 1000. Given a fixed number of sampling points, the resulting algorithm is dominated by a flat scaling of . Absolute correlation energies are accurate to <0.1 kcal/mol with respect to conventional calculations for several hundreds of electrons. This allows for the calculation of unbiased basis-set limit correlation energies for systems containing hundreds of electrons with unprecedented efficiency gains based on a straightforward treatment of continuum contributions.


INTRODUCTION
Electron correlation lies at the heart of a wide range of fundamental physical and chemical phenomena, which range from the structural diversity and dynamics of water 1 over the dissociation of liquid hydrogen at high pressure 2 and the stability and mobility of point defects in semiconductors 3 to the barrier height of chemical reactions. Wavefunction-based methods allow for a conceptually simple and convenient treatment of electron correlation 4−6 and have found widespread and long-lasting use in theoretical chemistry. Correlated wavefunction methods have been widely applied as a benchmarking tool 7,8 in the development of computationally more expedient methods such as the Kohn−Sham density functional theory (KS-DFT). 9,10 More recently, their scope has been enlarged by rigorous hybridization schemes that combine KS-DFT with correlated wavefunction approaches, 11−15 giving rise to some of the most accurate density functional approximations available to date. 12,13,16−19 In particular, while it has been pointed out that many recently developed density functional approximations fail to yield correct densities and energies at the same time, 20 double hybrid (DH) functionals have been shown to be able to overcome this fundamental problem. 21 Recently, modern machine learning techniques have considerably increased timescales and system sizes that can be sampled on conventional infrastructure, but the generation of reliable input data for the training of such methods still relies on the computational feasibility of reference calculations of sufficient accuracy. In this perspective, the importance and scope of wavefunction-based firstprinciples techniques applied to condensed matter systems can therefore only be expected to grow further, both as a standalone method and in combination with DFT.
To this day, wavefunction-based correlation methods are hampered by a scaling that is polynomial at best and that is associated with a considerable prefactor. This implies that for larger systems, trade-offs have to be made between the accuracy of the basis set employed and the number of electrons that can be treated with reasonable computational resources. Moreover, correlated wavefunction approaches have only scarcely been applied in the condensed phase, which is due to additional difficulties encountered in periodic systems. 22−40 These difficulties can be further exacerbated by the large basis sets needed to obtain basis-set limit reference energies. 41 This precludes the routine use of wavefunction-based methods for large condensed phase systems; at the same time, benchmarking possibilities, for example, against newly developed density functionals, remain restricted to comparably few atoms and small supercells 27,31,34,35,38 or have to be based on basis sets which are far from completeness. In benchmarking, this can be particularly problematic in combination with the erroneous convergence behavior of certain density functionals, which obfuscates any comparison that is not explicitly made at the complete basis-set limit. 42,43 The availability of basis-set limit results is therefore necessary not only for formal reasons, but is also of great importance for the assessment of the physical accuracy of existing models and approximations, representing an important guideline in the development of new techniques and approximations that are able to reach far beyond current system sizes and limitations.
In principle, plane wave (PW) basis sets would constitute an ideal choice for the calculation of basis-set limit correlation energies since they do not introduce any localization bias and allow for a controlled, simple, and well-defined extrapolation to the complete basis-set limit 44 without the linear dependency issues commonly encountered in large atom-centered bases. 28,35,45,46 In particular, since a single PW is the solution of the Schrodinger equation of a free electron, its use enables the description of continuum states, which have been shown to play a crucial role in the complete description of electron correlation. 47 In the following, we will refer to continuum states in finite systems as those virtual states that resemble a free electron; it is also this resemblance that lies at the heart of simple extrapolation to basis-set limit values. 44 By virtue of their very nature, conventional atom-centered bases such as Gaussian functions or combined Gaussian/PW (GPW) 48 bases are unable to describe such continuum states, which also account for the absence of physical models that would allow for a simple extrapolation to the basis-set limit. Instead, they usually rely on specifically constructed basis sets that allow for certain extrapolation models to be applied; this, however, does not commonly hold for density functionals. 42,49 PWs have not been reported to suffer from this drawback. 43 In addition, PWs are equally suitable for the treatment of both periodic and nonperiodic systems with either wavefunction-based methods, density functional techniques, or hybridizations thereof. These advantages, however, come at a price: The presence of a large number of continuum states in PW setups exacerbates the steep scaling of correlated wavefunction methods, making them computationally intractable for all but the smallest systems, which on their own will already require substantial resources on conventional high-performance compute clusters. In the following, we will show that PW-based correlated wavefunction calculations can be sped up by a factor of up to 1000 by stochastically sampling continuum state contributions via Monte Carlo summation. The error introduced by this stochastic approach remains below 0.01 kcal/mol per electron. This enables correlated wavefunction calculations in PWs for unprecedented system sizes on conventional computational infrastructure, making unbiased basis-set limit values routinely accessible for systems with up to hundreds of electrons. The same reflections hold for hybrid wavefunction/DFT methods, such as DH 11 density functionals.
1.1. Møller−Plesset Perturbation Theory. Among the correlated wavefunction methods, second-order Møller− Plesset perturbation theory (MP2) 50 exhibits a comparably flat scaling of N ( ) 5 with the number of electrons or basis functions N, making it one of the flattest scaling correlated, wavefunction-based approaches available, second only to the random phase approximation (RPA) and the direct RPA (dRPA), respectively. 51,52 In general, MP2 has been found to provide a good first estimate of the dynamic correlation energy. 6,53 Conceptually simple, the MP2 correlation energy E c MP2 is obtained by a perturbative treatment that includes up to doubly excited determinants and summation over pairs of all N occ occupied and N vir virtual orbitals. For the spin-restricted case, 4 where i, j and a, b denote spatial occupied and virtual orbitals ψ, respectively, that are eigenstates of the Hartree−Fock operator with eigenvalues ε. The MP2 matrix element E̅ (i,j,a,b) is expressed in terms of four-electron Coulomb energies, which can be cast into a positive-definite form where ⟨ij|ab⟩ are two-electron matrix elements. In DH density functionals, the second-order integrals in eq 1 are evaluated using a ground-state KS determinant, rather than the Hartree− Fock solution, 11 and will only contribute to a fraction of the total correlation energy, the remainder being treated by pure density functional methods. If the advantages of DHs are to be made routinely available in condensed matter applications, they too require an efficient treatment of the terms in eq 1, and any improvements in the calculation of the MP2 term will directly benefit DH calculations.
Historically, E c MP2 has been evaluated using atom-centered basis sets or mixed GPW 54 implementations. In particular, in periodic setups, use of localized basis functions has been reported to be susceptible to basis-set convergence issues, 28,35,45,46 whereas problematic basis-set superposition effects and possible linear dependencies are absent in a PW representation. Applications in solid-state physics have also been scarce, which is in part due to the divergence of the MP2 integrand in zero-band gap systems, but a Thomas−Fermi screening of the MP2 amplitudes can resolve this issue. 55 The presence of continuum (or continuum-like) states allows for simple extrapolation to the complete basis-set limit: Alavi and co-workers have demonstrated that in a PW basis at large N vir , E c MP2 decays as ∝ε a −3/2 . 44 PW-based methods therefore offer the unique advantage of a simple evaluation of basis-set limit values in both periodic and isolated systems, making them a potentially invaluable tool for basis-set bias-free calculations. PW MP2 calculations are few and have only recently been reported. 34,35,56,57 This is in part due to their extensive memory requirements. The presence of a number of virtual states close to the number of PWs themselves (up to 10 9 ) can further complicate the calculations. However, given sufficient memory, the integrals of eq 2 are easily evaluated in the reciprocal space: Exchange-like matrix elements are easily obtained by solving the Poisson equation for a set of pair densities ρ ia (r) = ψ i *(r)ψ a (r) in reciprocal space, where the Coulomb operator is diagonal. The real-space pair densities can simply be subjected to a fast Fourier transform (FFT), ρ ia (G) = FFT[ρ ia (r)]. Then, at the Γ-point, the reciprocalspace equivalent of a matrix element reads 58 where † stands for the complex conjugate and index permutation. Ω is the volume of the system, G are reciprocal-space vectors, and Φ(G) is a suitably generalized form of the reciprocal-space Coulomb potential. 59,60 Attempts to reduce the overall cost of the method have included mapping the virtual states onto a localized subspace, 37,61−65 the use of stochastic orbitals 66−69 and (real-space 70,71 or graphbased 57,72 ) sampling approaches, respectively, as well as exploiting Laplace transforms 53,56,73−76 to enhance parallel efficiency. In an alternative strategy, one seeks to accelerate convergence of the correlation energy by improving the description of the electron−electron cusp. To this end, explicitly correlated basis sets can be used. In an ansatz commonly referred to as the F12 theory, 77,78 the description of the electron−electron cusp is improved by combining a conventional, atom-centered Gaussian basis with a set of strongly orthogonal geminals. This approach has recently been extended to PWs. 36 However, localization procedures fail for continuum states; previously reported stochastic techniques either introduce system-dependent errors of up to 2 kcal/mol per occupied orbital 69 or carry a prefactor too large for practical applications, 70 and approaches based on the Laplace transform 56 require repeated calculations at different quadrature points, increasing the overall operation count for the sake of parallel efficiency, with modest reported speedups of around 4 to 5. A combined stochastic/Laplace-transform approach that results in appreciable computational speedups leads to errors of at least 20%. 57 Similarly, due to the presence of nonfactorizable many-electron integrals, the cost of evaluating the MP2 integrand in PW-F12 is higher than for a pure PW basis set, 36 and the efficiency of a graph-based approach 72 was hindered by the absence of an optimized weighting scheme for graph generation. In a more general scope, a recent diagrammatic decomposition of the coupled cluster pair correlation function has allowed for the introduction of a basis-set correction in PWs that results in speedups of 2 orders of magnitude. 79 Alternatively, in the context of Full Configuration Interaction Quantum Monte Carlo calculations, 80 use of an effective, transcorrelated Hamiltonian 81 has been shown to substantially improve convergence of the correlation energy of the uniform electron gas. In the GW theory, 82 stochastic sampling schemes have successfully been applied with competitive accuracy and favorable timings with respect to deterministic calculations. 83 In the following, we shall demonstrate how the presence of continuum states can be exploited to drastically reduce the computational cost of MP2 calculations without impacting their accuracy. The approach is based on a simple stochastic sampling of the integrands of eq 1 and can be implemented with little effort, representing a sleek and clean approach to tackle the issues arising in the continuum. This opens the path to routine applications of PW MP2-based approaches in both isolated and periodic systems with up to hundreds of occupied orbitals, making it possible to obtain basis-set limit DH or MP2 correlation energies on conventional computational infrastructure within a reasonable time.

DISTRIBUTION OF CONTINUUM STATES AND
STOCHASTIC SAMPLING The possibility of introducing stochastic sampling is rooted in the behavior of the integrand at large ε a , where continuum states arise. In this regime, E c MP2 grows as ε a −3/2 , 44 and the overlap elements ⟨ij|ab⟩ must therefore be of low magnitude. It is obvious from eq 3 that overlaps will only be non-negligible whenever the symmetries of the continuum states match, but explicit symmetry determination for all states would be prohibitively expensive. Instead, given two high-lying virtual states, the statistical distribution of non-negligible overlaps is expected to be similar between truncations of eq 1 at a and some subsequent a + δ with arbitrary δ. For a spatially infinite supercell at infinite PW cutoff, one can define a cutoff energy ε c with orbital index c, from where on all orbitals a ≥ c are part of the continuum. Then, the correlation energy can be separated with the continuum contribution η c (a) accounting for the incremental change in correlation energy when adding an additional continuum orbital to the system: where ∑ b′ is due to pairs of continuum (a) and noncontinuum (b′) virtual orbitals and we have used that, at the Γ-point, ψ a = ψ a *. Note that the last term in eq 4 contains only one explicit sum over virtual orbitals N vir , with the orbital pairs themselves being formed by the triple sum in eq 5. The simplest possible stochastic treatment of eq 5 is given by a uniform sampling of the summand, but this calls for a regular distribution of the overlap values obtained over all tuples ijab, which are the arguments of η c (a); that is, the highlying virtual orbitals of a finite system at finite PW cutoff need to reasonably approximate free, continuum electrons. Figure 1 shows the distribution of these arguments for a finite periodic box with ε c − ε HOMO = 50 eV and a total of N vir = 3000 virtual orbitals. With the positive-definite definition of E̅ adopted in eq 2, the resulting distribution is indeed regular and smooth, indicating that η c (a) could be predestinated to be treated by uniform Monte Carlo sampling. This is the approach we will privilege in the following. Note that in the limit of an infinite system, this is equal to Monte Carlo integration over the continuum; such integration techniques have been used as Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article early as 1957 to calculate the high-density correlation energy of the uniform electron gas. 84 Let p s ∈ [0, 1) be a predefined sampling probability (i.e., p s is larger than or equal to zero but always smaller than 1). We define the function p as Symmetry at the Γ-point, which gives rise to a factor of 2 in continuum/continuum and mixed continuum/noncontinuum terms of eq 5, can be accounted for by introducing a Kronecker delta. Then, the sums in η c (a) can be simplified, and the stochastic expression for η c (a) is where x ijab ∼ U([0, 1)) is a sequence of random numbers drawn from a uniform distribution U and b covers both b and b′ of eq 5. In the limit of a continuum, the error of the sampling is expected to decrease as N 1/ MC , with N MC being the number of tuples ijab that are explicitly sampled. The combination of eqs 4 and 7 amounts to a stochastic summation over all orbital pairs that contribute to an increase in E c MP2 , once a continuum orbital a is added to a system already containing a − 1 virtual orbitals.
For any sampling probability p s < 1, only elements of η c (a) with p(x ijab , p s ) = 1 have to be evaluated. In a finite cell with a finite PW cutoff and a sufficiently large number of virtual orbitals N vir , estimates for η c (a) are then obtained for every given continuum orbital a ∈ [c, N vir ] as follows: rather than randomly generating a tuple of indices for a predetermined number of Monte Carlo moves, a random number p(x ijab , p s ) is drawn for every element of the summand in eq 7, determining whether a particular tuple ijab will enter into the estimator of η c (a). p s therefore defines a target value N MC ∝ p s of tuples that are expected to be sampled for every continuum orbital a. The advantages of such an algorithm are twofold: On the one hand, it avoids under-or oversampling of the subspace associated with orbital a; on the other hand, it enables efficient extrapolation to the basis-set limit in one single calculation, directly providing E c MP2 as a function of the highest virtual orbital included in eq 4.
In the following, we will adopt an orbital-dependent sampling probability p s (a) = N MC /N card (a), where N card (a) = (2a − 1)N occ (N occ + 1)/2 is the product of the cardinalities of the sums in eq 7. N card (a) therefore explicitly depends on the virtual orbital a that is added to eq 4. With this choice of a continuum-orbital-dependent p s (a), the density of the sampling decreases as a increases. This allows for N MC to remain a fixed, system-independent input quantity. We shall later show that conservative estimates for N MC and ε c can be regarded as system-independent.
For orbitals that are part of the continuum, the resulting algorithm scales formally as N N N ( ) PW vir MC , where N PW is the number of PWs in the basis set. This follows from eqs 4 and 7 since the cost of evaluation of the triple sum over N occ 2 N vir in eq 7 is reduced to elements with p(x ijab , p s (a)) = 1, which in turn is proportionate to N MC . For virtual orbitals with ε a < ε c , the conventional N N N ( ) PW occ 2 vir 2 scaling applies (cf. eqs 1 and 3). Further on in the text, we will show that in practice, the number of terms due to eigenvalues ε a < ε c does not dominate scaling. Once N MC can be made both independent of the orbital index a and the system at hand, the scaling of the resulting method reduces to N N ( ) PW vir integral evaluations for all orbitals with eigenvalue ε > ε c .

COMPUTATIONAL METHODS
3.1. General Setup. Hard pseudopotentials of the Goedecker−Teter−Hutter (GTH) form 85 parametrized for Hartree−Fock calculations 31 have been used for all calculations. PW MP2 and stochastic MPs2 energies were calculated using a modified version of the CPMD code. 86 The convergence threshold on the residual gradient on occupied orbitals was set to 10 −7 a.u., whereas a threshold of 10 −5 a.u. was used for the virtual space. For isolated systems, periodic images were decoupled using the Poisson solver by Martyna and Tuckerman. 87 The wavefunction cutoff energies E cut ψ were set to 150 Ry for all systems but for the ethylene crystal, where a value of 140 Ry was used. A density cutoff E cut ρ = 4E cut ψ was adopted for all systems, while cutoff energies for MP2 pair densities were set to E cut ψ without impacting accuracy. The calculation of the MP2 term is based on straightforward evaluation of eq 3 as in ref 58. Since no derivatives with respect to E c MP2 have to be evaluated, reciprocal-space orbital pairs ρ ia (G) are stored in memory, allowing for substantial speedups with respect to an on-the-fly evaluation of every pair density. The size of the periodic super-or unit-cells, the number of electronic states as well as molecular geometries, and PW energy cutoffs are given in the Supporting Information.
3.2. Extrapolation of Correlation Energies. Correlation energies are obtained from single-point extrapolation using the ε N vir −3/2 dependency (equivalent to 1/N vir ) as described in ref 44. Such a leading-order approximation requires a sufficient number of high-energy orbitals in order to gather sufficient statistics and reliable extrapolated values. In this work, we adopt a fitting scheme with various windows (ranges of points to fit) moving along the curve. Those vary in size with respect to the number of orbitals included, with a shift of ∼200 orbitals between them. The maximum window size is taken as the one that retains a fitting error comparable to smaller windows when finishing fitted ranges at N vir max , the highest orbital available from a Davidson diagonalization. The smallest window is given by the smallest possible window size that provides a stable fit.
MP(s)2 curves are fitted according to The series of α(N vir ), obtained from all windows along the curves that terminate at N vir up to N vir max , converge to an estimate of the MP(s)2 energies at the basis-set limit, E c MP(s)2 (ε N vir → ∞) Averaging over different windows allows to account for sensitivity and variance of extrapolated values, which are given in the Supporting Information along with figures that illustrate the extrapolation procedure and error determination.

RESULTS AND DISCUSSION
4.1. Accuracy of Stochastic Summation. Our stochastic sampling scheme was tested both on isolated systems as well as in periodic, condensed-phase setups. Test systems in the Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article condensed phase include solid-state ethylene and benzene molecular crystals as well as a liquid consisting of a hydronium ion solvated in 32 water molecules. Isolated (gas-phase) systems are represented by a benzene monomer and a dimer in sandwich configuration. All calculations were carried out with a modified version of the CPMD code 86 using GTH pseudopotentials generated for Hartree−Fock calculations. 31,85 We first investigate the dependency of the accuracy of our stochastic sampling of eqs 4 and 7, called MPs2, on the number of orbitals in the occupied subspace and the choice of the continuum cutoff value ε c gap . This quantity is given with respect to the highest occupied molecular orbital (HOMO) in order to be independent of the system setup or reference frame: ε c gap = ε c − ε HOMO . Figure 2 shows the absolute difference between MP2 and MPs2 calculations for a benzene crystal and a benzene monomer for ε c gap ranging from 20 to 180 eV. Differences rapidly decrease by increasing ε c gap . Errors averaged over independent stochastic runs remain <0.01% for continuum cutoffs ≥120 eV. Reducing this value by half to 60 eV results in doubling of the relative error, which can still be acceptable. Further lowering of ε c gap , however, leads to rapidly increasing errors. The standard deviation of the sampling error depends much more strongly on ε c gap than on N occ . From ε c gap = 120 eV on, standard deviations become negligible and virtually identical for both systems. Notably, relative sampling errors are lower for the crystal at N occ = 60 than for the monomer (graphs are provided in the Supporting Information). This strongly supports that accuracy is mainly influenced by ε c gap , whereas at constant N MC , the error is independent of the partitioning of occupied and virtual states in eq 7.
In the following and in line with the values of Figure 2, we will adopt N MC = 12,000 and ε c gap = 120 eV to investigate the accuracy of stochastic sampling for different systems. Figure 3 shows the correlation energy as a function of the eigenvalue of the highest virtual orbital included (ε N vir ) for two exemplary systems: one periodic (ethylene crystal) and one isolated (benzene monomer) setup, calculated both with a full summation according to eq 1 as well as with our stochastic sampling. The corresponding extrapolated basis-set limit values are shown in Table 1. The largest errors of the extrapolated, absolute MP2 correlation energies lie between 0.02 and 0.1 kcal/mol. Errors in the binding energy of the benzene sandwich are of comparable magnitude, which is far below chemical accuracy. Differences per electron, ΔE/e − , do not exceed 0.1 meV. These values compare well to an expected stochastic error of ≤0.01%. Table 1 also shows the number of matrix element evaluations for all setups. With our systemindependent choice of N MC , stochastic sampling reduces this number by 1−2 orders of magnitude compared to conventional calculations.
For both absolute energies and energy differences, the observed deviations are several orders of magnitude lower than the reported error of other stochastic or Laplace transformbased schemes. 57,66−72 For comparison, values of correlation energies obtained with atom-centered all-electron 88 and GPW codes 89 are listed in the Supporting Information; all values are consistently higher than those reported for our PW calculations. In particular, we note that differences between different basis sets tend to be substantially larger than the stochastic sampling error, further confirming the viability of uniform, stochastic sampling of the continuum space.

Performance and Speedups.
With the error of the stochastic sampling scheme being considerably lower than the errors documented for other PW implementations, effective speedups remain to be determined. Figure 4 shows the resulting cumulative execution times and speedups of the stochastic sampling compared to the direct implementation of eqs 1 and 3 for the ethylene crystal. Timings are reported as a function of the highest orbital index included in the expansion, N vir . For the ethylene crystal, at N vir = 10,000, speedups of up to 957 can be reached with stochastic summation, making the calculation about 3 orders of magnitude faster compared to a conventional implementation. All the while, the error introduced by uniform stochastic summation with respect to a full calculation is only about 10 −2 kcal/mol for this system (cf. Table 1). This has to be compared to maximum speedups of about 5 documented for Laplace transform-based schemes that allow for similar accuracy to be retained 56 and 20% errors in correlation energies for algorithms that allow for larger speedups 57 versus errors around 0.01% reported here. Figure 4 also includes a fit of the CPU time of our stochastic scheme to a function O(N vir ) = a 0 + a 1 N vir , demonstrating that the scaling is described well by N N ( ) PW vir , assuming that N MC = cnst. N ( ) 2 constitutes a formal improvement over the scaling achieved with a stochastic graph-based approach in atomcentered basis sets, 72 which was reported to be of N ( ) 2.6 . Together with the high accuracy of the method, this considerable gain in efficiency allows for the treatment of systems that would be intractable when treated with conventional algorithms. One typical usage example of accurate wavefunction-based theories or computationally demanding  , are plotted on a secondary y-axis. Extrapolated basis-set limit values for all systems described here are found in Table 1.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article high-quality DFT methods such as DH functionals lies in postprocessing of simulation data, for example, in a posteriori calculations of potential energy surfaces or reaction paths generated using lower-level methods in the context of molecular dynamics or Monte Carlo simulations. Recent developments have even made it possible to directly sample 90 high-quality potential energy surfaces by virtue of multiple time step (MTS) propagators, 91 which allow for an important decrease in computational cost and substantial improvements in tractability by permitting less-frequent evaluation of the full, high-quality Hamiltonian during a first-principles molecular dynamics run. A liquid constituted of 1 hydronium ion solvated in 32 water molecules (264 electrons) will serve as an example of the performance of stochastic sampling in a typical setup encountered in first-principles (MTS-)molecular dynamics. Using the stochastic summation scheme reported here, calculating the basis-set limit correlation energy of this system, shown in Figure 5, is feasible in about 15 h on 25 16-core compute nodes with 128 GB of RAM, using the same, conservative estimates for ε c gap reported in Table 1, which yielded accurate results for all test systems considered so far. Additionally, calculations using ε c gap = 60 eV and ε c gap = 90 eV have been carried out for the sake of comparison. Already at ε c gap = 90 eV, the execution time is almost halved to about 8 h and can be further reduced by using ε c gap = 60 eV, at which the basis-set limit correlation energy can be calculated in a mere 5 h. In particular, this drastic reduction in execution time is not accompanied by a considerable loss of accuracy. Extrapolated correlation energies are given in Table 2. It should be noted that postprocessing protocols can be applied when training machine learning algorithms with high-quality data, based on a coarse sampling of configuration space with a lower-level method (cf., e.g. ref 1). Speedups in the calculation of correlation energies with MP2 or DH functionals will therefore be directly reflected in less time-consuming training procedures, thus considerably increasing throughput.
4.3. Generalization to RPA. The promising performance of the stochastic summation scheme described here also opens the possibility of its application to similar approaches in which continuum states can play a role. In DFT, the exact-exchange plus RPA approach has emerged as a promising method capable of more accurately predicting van der Waals binding energies, adsorption energies on surfaces, or lattice constants in molecular solids. 44,54,92−94 Based on the similarity of the MP2 energy expression and the RPA, one can expect transferability of the stochastic sampling approach to the evaluation of the and ΔE/e − denote absolute and per-electron energy differences between stochastic sampling and conventional calculations, respectively. N ijab denotes the number of matrix elements (ijab-tuples) sampled in a conventional calculation, and p s N ijab is the number of effectively sampled matrix elements in a stochastic MPs2 calculation, with both numbers rounded to three digits. The threshold for stochastic sampling expressed with respect to the HOMO was identical for all systems, ε c gap = ε c − ε HOMO = 120 eV. The same holds for N MC = 12,000, the number of terms sampled per virtual contribution η c (a), as described in eq 7. For details on the systems used, cf. the Supporting Information.
Tr ln 1 ( , ) ( , ) where χG ,G′ (q,iω) are elements of the full density response function, including the Coulomb interaction. At the Γ-point with q = 0, the diagonal elements in a stochastically sampled RPA scheme become where the variables c, x jb , and p s are used analogously to eq 7.
In this limit, we investigate the RPA integrand of eq 9 for an ethylene crystal at different values of ω. Figure 6 shows the values of E c RPA at varying iω for N vir = 11,040 as well as an estimate of the overall RPA correlation energy based on trapezoidal integration. Stochastic sampling with p s = 1/3 introduces a maximum error of about 1% in the integrand. Overall, the stochastic sampling introduces a final error of less than 0.03 kcal/mol, which is comparable to the error obtained in MP2 calculations for the same system. These results demonstrate that the stochastic sampling of continuum states can also be applied for methods other than MP2 that include a substantial continuum-state contribution.

CONCLUSIONS AND OUTLOOK
Continuum states have been shown to be an important contributor to the overall electron correlation energy. 47 Among the basis sets commonly used in solid-state physics, quantum chemical calculations, and first-principles molecular dynamics, PWs stand out as the only choice that can effectively account for continuum contributions, which are also the base of a simple basis-set limit extrapolation technique. 44 Here, we have introduced a uniform stochastic sampling approach to treat continuum states arising in correlation energy calculations, where contributions due to states with orbital eigenvalues beyond some threshold ε c are added by stochastic summation. This algorithm has been applied to the calculation of secondorder perturbation energies which occur in both MP2 and DH density functionals. We have shown that stochastic summation over the continuum orbitals allows for the calculation of MP2 correlation energies with speedups of up to 3 orders of magnitude at remarkably low errors. This significant increase in efficiency enables calculations with several hundreds of electrons at a relatively low computational cost, making it possible to standardly apply MP2 and DH methods in a PW basis, which has so far been intractably expensive even on highperformance compute clusters. Importantly, the results presented here also enable straightforward DH calculations in the condensed phase, thus extending the availability of one of the most accurate density functional methods available to date for condensed matter systems. We have also shown that stochastic sampling of the continuum orbitals can easily be extended to other approaches, demonstrating the generality of the ansatz employed here. Calculations carried out within a stochastic RPA suggest errors comparable to the stochastic MP2 scheme. The stochastic sampling scheme itself is straightforward, easy to implement, and based on simple physical concepts.
Stochastic sampling of continuum states permits to easily obtain basis-set limit values, be it for periodic or isolated systems, with maximum errors of only 0.1 meV per electron. This accuracy makes basis-set limit values for correlation energies available using reasonable computational resources and execution times. This will allow for thorough benchmarking of new computational methods without basis-set bias, for routine postprocessing of potential energy landscapes generated using lower-level methods, as well as on-the-fly generation of high-accuracy first-principles molecular dynamics trajectories using multiple time stepping schemes. This data, in turn, can be used to feed high-throughput methods based on artificial intelligence. Overall, the techniques presented in this text pave the road to routinely apply accurate MP2 and DH calculations at the basis-set limit in condensed matter systems, ultimately extending the use of a method well-established for isolated systems to the condensed phase.  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article