Accelerating the convergence of auxiliary-field quantum Monte Carlo in solids with optimized Gaussian basis sets

The accurate theoretical description of complex materials is one of the grand challenges of chemistry and physics. Mean field methods, such as density functional theory^1,2 (DFT), are by far the most widely used approaches, offering often satisfactory accuracy relative to their computational cost. Despite this, DFT suffers from a number of well-known challenges, including the choice of approximate exchange–correlation functional³ and the somewhat ad hoc extensions necessary for the treatment of strongly correlated systems.⁴ Motivated by these limitations, new methodological developments, and the ever increasing computer power available,⁵ there has been a renewed interest in applying wavefunction based quantum chemistry methods in solids.⁶ Although typically orders of magnitude more expensive than DFT, wavefunction methods are appealing as they form a systematic hierarchy of approximate approaches that yield a higher accuracy albeit at an ever increasing cost.

Perhaps, the most widely developed solid state quantum chemistry methods to date are the second-order Møller–Plesset perturbation theory^7–14 (MP2) and the Random Phase Approximation (RPA).^15–18 Coupled cluster methods have also been recently developed and applied to solid state problems^19–21 (CCSD). In parallel with this, orbital-based statistical approaches such as full configuration interaction Monte Carlo^22,23 and, in particular, auxiliary-field quantum Monte Carlo (AFQMC)^24,25 have begun to be widely applied in solids.^6,26–34 Despite this rapid development, the widespread adoption of these approaches has in part been hindered by the lack of compact basis sets available. This is critically important if experimentally accurate relative energies are desired as all of these approaches scale with a high power of the basis set size.

Many different approaches have been developed over the years to accelerate the convergence of wavefunction methods in solids with respect to the single particle basis set and system size. These include constructing “downfolded” Hamiltonians,^31,32 using natural orbitals from MP2,³⁵ employing explicitly correlated approaches,^19,36,37 developing mixed Gaussian-plane-wave approaches,³⁸ the use of MP2 to find optimal twist angles to reduce size effects,³⁹ and using correction schemes from lower order theories.^40,41 Despite these developments, no one approach is ideal for every method and convergence to the complete basis set limit remains slow. In addition to the slow convergence with basis set size, approaches based on GTOs optimized for the atomic environment are also known to suffer from linear dependency issues when transferred to the solid.^42–45 This is particularly a problem when approaching the complete basis set limit in general and when investigating high pressure phases of materials.

Recently, it has been suggested that custom built GTO basis sets tailored for the solid environment often outperform standard basis set libraries for DFT calculations.⁴⁶ In this paper, we explore a similar approach but instead targeted toward generating correlation-consistent basis sets for correlated methods. In essence, we propose optimizing the exponents through the minimization of the MP2 correlation energy computed in the unit cell of the material. The resulting basis sets are then used for this system in different conditions, e.g., different supercell sizes, similar structures, and different volumes. Although our approach is far from ideal, as a new basis set needs to be generated for each new application, we argue that this is a price worth paying, given that highly accurate approaches are only ever going to be applied in targeted, relatively simple applications for the foreseeable future. Moreover, the cost of the subsequent correlated calculation almost always dwarfs the cost of generating the basis set itself.

In this paper, we focus on deploying these basis sets in AFQMC calculations, but we stress that these basis sets can be used with other highly accurate many-body methods. By studying a range of insulating solids, we show that these optimized GTOs outperform “downfolded” PBE⁴⁷ orbitals, atomic GTOs, and yield comparable correlation energies to those found using MP2 natural orbitals³⁵ at a fraction of the cost. We also suggest the use of a simple MP2-based basis set correction scheme to further accelerate convergence that is independent of the underlying wavefunction method in contrast to that suggested in Ref. 41. With these developments, we are able to reliably produce AFQMC results for the energy in the complete basis set and thermodynamic limit for LiH, Si, and MgO. Our results for the bulk modulus, lattice constant, and cohesion all agree excellently with experiment.

The rest of this article is outlined as follows: In Sec. II, we outline our method for generating optimized basis sets for correlated calculations in solids. In Sec. III, we describe our MP2-based basis set correction and analyze its performance in a range of settings. In Sec. IV, we compare the performance of the various basis sets and present results for the equilibrium properties of LiH, Si, and MgO. Finally, we close in Sec. V with some discussion on the remaining challenges, including the adaptation of our approach to metals and the development of better finite-size corrections for orbital-based approaches.

Orbital-based quantum many-body methods are typically formulated starting from the Hamiltonian in second quantization,

where $ĉiσ†$ and $ĉiσ′$ represent the fermionic creation and annihilation operators in a given single particle orbital (SPO) set and h_ij and v_ijkl are the one- and two-electron matrix elements of the Hamiltonian.

The development of appropriate basis sets for electronic structure calculations has a long and rich history in both quantum chemistry and condensed matter communities. This is particularly true in the case of molecular calculations where the development of optimized basis sets has been extended to almost all aspects of the electronic structure calculation including electron correlation,^48–50 density fitting,^51–53 F12 and other explicitly correlated approaches,^54–58 core-valence correlation,⁵⁹ molecular properties,⁶⁰ and relativistic effects,⁶¹ among many others. On the other hand, for periodic calculations in condensed matter systems, the work has mainly focused so far on the development of good basis sets for mean field methods such as DFT. Due to the low cost of DFT calculations in comparison with typical many-body methods such as MP2, RPA, CCSD, or AFQMC, the development of DFT basis sets is typically focused on the efficiency of computation, robustness in convergence, and generic applicability. Typical choices for basis sets in this case include plane waves,^62,63 periodic versions of the linear combination of atomic orbitals approach,^64–67 and wavelets,⁶⁸ among others.

While the above approaches can work out of the box for any correlated wavefunction approach, they all suffer from different drawbacks. Ideally, we would employ a plane-wave approach, which is perhaps the most natural basis for simulating periodic systems and would allow us to build on decades of experience and developments. In addition, using plane waves would give us access to robust and efficient codes^69–72 as well as highly developed databases of pseudopotentials⁷³ and projector augmented wave potentials.^74,75 Unfortunately, despite offering a robust and straightforward way to converge the total energy with respect to the basis set size, a pure plane-wave approach would be computationally demanding, given the large plane-wave cutoffs required to converge the electron–electron cusps.⁷⁶ This is crucially important as AFQMC calculations scales like $O[M3]−O[M4]$ (with a large prefactor), where M is the number of single particle basis functions. One method to overcome this is “downfolding”^31,32 Kohn–Sham or Hartree–Fock states from a prior plane-wave calculation, i.e., working in an active space of eigenstates of the prior mean field calculation. As we will see, this approach typically converges slowly with the number of Kohn–Sham states and relies on good error cancellation, which may be difficult to achieve uniformly as conditions are changed, for example, by changing volumes in an equation of state calculation.

The use of localized basis sets for periodic calculations also has a long history,⁶⁴ which has mostly been focused on mean field approaches such as Hartree–Fock and DFT.^43,65,67 In this case, typically, one employs custom-built basis sets for solid state DFT simulations⁶⁶ or periodized atomic GTO basis sets.^43,44 Unfortunately, basis sets optimized for DFT simulations are not designed for correlated calculations. In particular, they do not follow the usual angular momentum hierarchy developed for atoms and molecules and thus cannot be used to reliably extrapolate to the Complete Basis Set (CBS) limit.³⁴ Although it is possible to modify these basis sets with basis functions from atomic Dunning basis sets,^14,34 this procedure is error prone and suboptimal. On the other hand, directly using all-electron correlation-consistent basis sets developed for atoms and molecules is also problematic, given the existence of diffuse functions. These basis functions typically lead to linear dependency issues, making the initial mean field simulation difficult to converge. The common strategy of removing basis functions with small exponents often helps; however, it can lead to subsequent difficulties with basis set extrapolation. These issues are well known and a number of alternative approaches have been developed, including using local-correlation methods⁴⁵ or natural orbitals from MP2³⁵ (MP2-NO) or RPA⁷⁷ (RPA-NO).

Here, we seek a practical middle ground that allows us to use existing plane-wave methodologies while wedding these with more compact localized basis functions designed from the outset for correlated solid state simulations. Our guiding principles are as follows: we desire a scheme that (1) has general applicability, (2) can be applied independently of the underlying basis used in the mean field calculations, (3) maximizes the amount of correlation energy for a give basis set size, and (4) lends itself to corrections from lower levels of theory such as MP2 or the RPA. We build on the hybrid PW-PGTO approach recently developed by Booth et al.³⁸ by not relying on a pre-existing Gaussian basis set for the virtual orbitals; instead, we optimize the exponents in the basis (and possibly the contraction coefficients if desired) by minimizing the MP2 energy of a small representation of the periodic solid, typically a single primitive unit cell. Similar to the approach developed by Booth et al., the current approach allows us to use a fully converged set of occupied orbitals in the reference calculation, which improves the transferability of the basis across multiple environments and also allows us to use any representation of the electron–ion interaction without the need to have the PGTO basis involved in the reference calculation. As we will see, the minimization of the MP2 energy leads to a set of virtual states that efficiently captures electron correlation while at the same time avoiding most of the problems associated with the use of atomic basis sets in solids, such as the ubiquity of small exponents that result in catastrophic linear dependency problems. We note that while the MP2 energy is non-variational and also typically overestimates the correlation energy, it still provides an appropriate optimization cost function for the basis set of virtual orbitals as we show below.

Our approach to generate a basis set is as follows:

1. We first perform a converged PW DFT or HF calculation to compute the occupied states.

2. We take the valence states from an existing correlation-consistent atomic basis set for the atomic species under consideration and generate the appropriate PGTO orbitals on a FFT grid.

3. We orthogonalize this PGTO basis against the occupied orbitals. This set of occupied and virtual orbitals now makes up our basis set.

4. We next perform a HF calculation within this new set of orbitals.

5. We compute the MP2 energy from the resulting HF eigenvalues and eigenvectors.

6. Iterate steps 2–5 by updating all the exponents for the valence states until the MP2 energy is minimized.

In order to avoid developing our own optimization routines, we created a Python interface for Quantum Espresso and relied on the SciPy package to perform the optimization of exponents in step no. 6. We used a finite difference scheme to calculate energy gradients. For the smaller basis sets, we observe a significant dependence of the optimized exponents on initial conditions, pointing to the existence of multiple local minima in the MP2 energy as a function of basis set exponents. As a result, we had to resort to global optimization algorithms in order to obtain a smooth and optimal improvement of the energy with increased basis set size. In particular, we used the basin hopping global optimizer from SciPy and found it to be sufficient. As the basis set size grows (e.g., pV5Z and beyond), the problem with multiple minima becomes less severe.

As we are looking for a relatively efficient and computationally inexpensive approach, we only reoptimize a basis set for an atom when the environment changes significantly. This means that in calculations within similar phases, we reuse the basis set generated at one volume for all volumes and system sizes. In the work presented here, all basis sets were optimized for a single primitive unit cell at the equilibrium volume of the material and reused everywhere. We also only considered single Gaussians without contraction for the virtuals. In a future publication, we will investigate the transferability of the resulting basis for different phases, different combination of elements, and other optimization strategies including the use of contracted Gaussians. For now, we are interested in comparing this approach with two other popular and general approaches developed so far, namely, downfolded Kohn–Sham states and MP2-NOs. We do not compare against RPA-NOs, but we expect similar performance in the materials studied here.

As an example of the utility of our approach, in Fig. 1, we show the error in the MP2 correlation energy, measured with respect to an MP2 calculation using all the plane-wave orbitals within a 70 Ry cutoff, as a function of the number of single particle orbitals of a two atom primitive cell of LiH in the B1 structure at the equilibrium volume. The results are shown for PBE orbitals, MP2-NOs, and the MP2-optimized hybrid PW-PGTO basis (from now on referred to as the optimized GTO basis (OPTM-GTO) for simplicity). As observed by others in the case of LiH,^35,38 the convergence of the MP2 energy with mean-field orbitals, either KS states such as PBE or HF orbitals (not shown here), is extremely slow compared to the convergence on the MP2 natural orbital basis. While convergence on the latter basis set is always expected to be faster, the different rate of convergence on LiH is quite dramatic. The convergence of the optimized GTO basis is comparable to the MP2 natural orbital basis, only showing a larger error in the correlation energy for the cc-pVDZ basis. In this case, it appears that either option represents a reduction of at least a factor of four in the number of states needed to reach a given accuracy, which would imply reductions of over an order of magnitude in AFQMC calculations at the same accuracy.

In addition to fast convergence, another important component of an efficient basis set for correlation is its capacity to provide a framework for basis set corrections, either through extrapolations or in combination with less expensive methods. While basis set extrapolations are very successful in molecular calculations, its utility in condensed matter systems is still unclear. In fact, without properly developed correlation-consistent basis sets for solids, this is in general not possible.³⁴ Instead, we will rely on less expensive approaches that offer a basis set dependence that is compatible with AFQMC. In this work, we explore the use of MP2 for such a purpose. As is the case with basis set development, the concept of basis set corrections based on MP2 for more expensive many-body methods such as CCSD(T) has a long history in the quantum chemistry community,^78–81 and its extension to condensed systems is straightforward.

The simplest MP2 correction adds the difference between the MP2 energy at the complete basis set limit and its value at the same basis used in the AFQMC calculation, namely, $ΔEMP2(M)=EMP2CBS−EMP2(M)$⁠, where M is the number of single particle basis functions. This correction should be reasonable as long as MP2 is applicable to the system since MP2 is known to diverge for metallic systems and will be inaccurate for systems with significant static correlation. Overall, we can expect this correction to overestimate the basis set incompleteness error of AFQMC since MP2 tends to overestimate the correlation energy in general.

An improved correction can be obtained following the recent work of Irmler and Gruneis^40,41 on MP2 basis set corrections for coupled cluster theory. In their work, they showed that the basis set incompleteness error in CCSD theory is dominated by the MP2 correlation energy and the particle–particle ladder term in the diagrammatic expansion of the CCSD equations, which led them to propose the following correction for CCSD:⁴⁰ 

where E^ppl is the energetic term associated with the particle–particle ladder (ppl) diagram.⁸² Even though the equation is approximate and based on somewhat empirical observations, Irmler and Grüneis subsequently showed that it can be used successfully to significantly reduce basis set incompleteness errors in CCSD.⁴¹ Looking at Eq. (3) more closely, we note that two approximations are made: (1) the MP2 and ppl terms have the same dependence with basis set size and (2) the rest of the CCSD energy is independent of the basis set.

Following similar approximations in this case for the AFQMC energy, we can write

where α is a system dependent constant that accounts for the part of the energy that scales like E_MP2(M) and E₀ represents those terms that we assume to be insensitive to the basis set due to their weak dependence. With this in mind, we can rewrite Eq. (3) as

which depends only on α and can be seen as a scaled version of the original MP2 correction. If we compare this expression with Eq. (3) in the case of CCSD, we can identify $α=1+EpplEMP2$⁠, which can be directly evaluated at any given point. In AFQMC, unfortunately, we do not have a clear way to calculate it directly. Instead, we use Eq. (4) and solve for α and E₀, given the value of the AFQMC and MP2 energies for two different vales of M. In the rest of this paper, we refer to Eq. (5) as the “scaled MP2 correction,” where α is calculated using results from pVTZ and pVQZ basis sets (or the equivalent number of states with PBE and MP2-NO orbitals) unless otherwise stated.

As an example of this, Fig. 2 shows the error in the AFQMC correlation energy of MgO as a function of the number of orbitals in the OPTM-GTO basis for a single primitive cell of the solid in the B1 phase. The CBS limit of the AFQMC energy was obtained from an extrapolation using the results from the cc-pV5Z and cc-pV6Z basis sets. The figure shows the results for the uncorrected AFQMC energies as well as those with the two corrections discussed above, the unscaled and scaled MP2 corrections. A single value of α was used, calculated using energies from pVTZ and pVQZ basis sets. It is clear that at least for this material, the scaled correction is superior to the bare correction, providing results within error bars of the extrapolated CBS energy already with a cc-pVTZ basis. While this is a particularly successful application of the correction, below we show evidence that the scaled version is in general a superior correction for basis set incompleteness error.

In this section, we compare the performance of PBE, MP2-NO, and OPTM-GTO orbitals in the calculation of the equation of state of LiH, Si, and MgO near equilibrium. We used a locally modified version of Quantum Espresso⁶⁹ to perform HF, DFT, and MP2 calculations (including the generation of MP2 natural orbitals), as well as to generate the one- and two-electron Hamiltonian matrix elements for AFQMC. PySCF^83,84 was used to evaluate the PGTO orbitals in the FFT grid. We used optimized norm-conserving pseudopotentials^85,86 for all atoms. As described in Sec. II, we generated MP2-optimized PGTO basis sets for all atoms by minimizing the MP2 energy of a single primitive unit cell of each material at the equilibrium volume, and this basis set was then reused for all volumes and system sizes presented here. The optimized basis sets have been included in the supplementary material. We used the phaseless AFQMC²⁵ method as implemented in QMCPACK.^87,88 We used the hybrid method for weight updates²⁶ and used the pair-branch population control algorithm.⁸⁹ All AFQMC simulations used a single determinant restricted Hartree–Fock trial wavefunction. The AFQMC input and output data are available at Ref. 90.

Figure 3 shows the AFQMC correlation energy for LiH, Si, and MgO, for a primitive unit cell of the material at the equilibrium volume, as a function of the number of orbitals in the single particle basis set. The results are shown for PBE, MP2-NO, and OPTM-GTO orbitals in all cases. The dependence of the correlation energy on the basis set size will show weak dependence on system size (e.g., as a function of the number of orbitals per atom), which makes calculations on a single primitive cell a useful tool in determining the quality and convergence rate of a given basis. In the case of LiH and MgO, both MP2-NO and OPTM-GTO offer a far superior alternative to downfolded PBE orbitals. We observe a similar performance in these two cases from HF orbitals and from DFT orbitals with other typical exchange–correlation functionals such as local-density approximation (LDA) and PBE0, implying that orbitals obtained from mean-field SCF calculations offer a poor choice on these materials. In the case of silicon, all three orbitals offer an acceptable convergence rate, with slightly superior results from MP2-NO. We believe that the somewhat slower convergence rate of the OPTM-GTO orbitals is the result of using uncontracted periodic Gaussians in the basis. We are currently exploring the use of an optimized contracted basis to further improve performance, which we will report in a future publication.

Even though MP2-NO and OPTM-GTO orbitals seem to provide a fast convergence of the correlation energy, the high cost of correlated calculations in solids combined with the need to account for finite-size effects, a highly expensive computational task, makes the development of basis set correction schemes a high priority. Figure 4 shows the error in the AFQMC correlation energy for the same systems shown in Fig. 3 but with the inclusion of the ΔE_MP2 and $ΔEMP2scaled$ corrections described in Sec. III. In all cases, the scaling factor in $ΔEMP2scaled$ was obtained with the results of pVTZ and pVQZ basis sets (or the equivalent number of orbitals in PBE and MP2-NO). The results are plotted with respect to extrapolations to the complete basis set limit from either MP2-NO results (LiH and Si assuming a 1/M dependence) or from OPTM-GTO results (for MgO assuming a typical 1/L³ dependence). As expected from the discussion in Sec. III, MP2’s tendency to overestimate the correlation energy will generally lead to a correction that is too large. The scaled MP2 correction shows a better performance, leading to errors on the order of a few mHa already for pVTZ basis sets. All corrections vanish in the complete basis set limit, leading to the same result at large enough basis set size (this is not always clear with PBE orbitals due to the slow convergence). Nonetheless, both corrections offer a significant reduction in the basis set incompleteness error, particularly for the smaller basis sets where the scale of errors can be reduced by over an order of magnitude. In addition, both corrections are able to reduce the error in calculations with PBE orbitals to magnitudes comparable to the other two orbital sets. We, nonetheless, still prefer orbitals that directly capture a larger amount of correlation since with PBE orbitals the correction will be larger than the AFMQC correlation energy with smaller orbital sets.

The comparison of these basis sets at a single volume shows similar performance between MP2-NO and OPTM-GTO orbitals, which are both typically superior to mean-field SCF orbitals. If we are interested in calculations at a single volume, both options seem perfectly acceptable. When we are interested in thermodynamic properties of the material that depend on volume, such as lattice constants, bulk moduli, or the equation of state in general, the more relevant question is not the rate of convergence of the energy at a given volume but rather the cancellation of errors between different volumes. We are interested in a basis set that can approach the correct volume dependence of the correlation energy with as few terms as possible. In the quantum chemistry community, these are known as nonparallelity errors in the case of dissociation curves of molecular complexes.

Figure 5 shows the AFQMC correlation energy of MgO as a function of the lattice constant (the experimental lattice constant, corrected for zero-point effects, is 7.916 bohr) for all basis sets and varying number of orbitals. In order to better visualize the nonparallelity errors, we plot the error in the correlation energy, measured with respect to the complete basis set extrapolated value at every volume, with respect to its value at the lowest volume, namely, ΔE(a) = (E(a) − E^CBS(a)) − (E(a_min) − E^CBS(a_min)). As currently presented, the plotted value is a combination of the errors at the current and reference volumes. Nonetheless, it clearly shows the magnitude of nonparallelity errors and allows for a clear comparison between various orbital sets. In addition to bare AFQMC energies in the first row, the figure also shows results with the unscaled and scaled MP2 corrections in the second and third rows, respectively. The shaded region in each plot represents the error bar of the basis set extrapolation and contains contributions from both the current lattice constant and the subtracted reference value. Similar figures for LiH and Si can be found in the supplementary material.

Several points are immediately clear from Fig. 5. On the one hand, PBE orbitals not only have significant nonparallelity errors, but they also decrease very slowly in magnitude with basis set size. Even though we decided not to pursue convergence of relative energy differences with PBE further in the case of MgO, it is apparent that over 1000 orbitals will be needed to reach convergence in the mHa range. OPTM-GTO orbitals, on the other hand, show much more reasonable errors across the volume range, and even small basis sets already produce acceptable results depending on the application. In the case of MP2-NO orbitals, we observe a somewhat erratic behavior in MgO, which is not observed in either LiH or Si. In the latter cases, both MP2-NO and OPTM-GTO orbitals produce comparable results and are far superior to PBE in the case of LiH. At this point, it is not clear what is the cause of the unusual behavior of MP2-NO orbitals in MgO as a function of volume, but it is clear that their performance is quite poor in this material. When MP2 corrections are applied to all three orbital sets, the magnitude of the errors is significantly reduced in all cases, particularly in the case of PBE orbitals. In this case, the scaled correction does not necessarily outperform the unscaled version, even though it is usually superior on total energies. Note that the scaling factor is calculated for every volume and system size independently. This implies that while the scaled correction leads to a more accurate estimate of the basis set incompleteness error at each volume, it does not mean that its volume dependence is better. Nonetheless, both corrections seem to always provide a reduction in errors, of significant magnitude in some cases. The picture in LiH and Si is consistent with the discussion above (see the supplementary material) and with the results in Fig. 3.

The results presented in the Secs. IV A and IV B indicate that the performance of each orbital set is heavily dependent on the system. PBE orbitals, for example, perform quite poorly on LiH and MgO but are a good choice in Si. It is hard, in practice, to determine which orbital set is appropriate for a given system without more experience and testing. On the other hand, mean-field SCF orbitals such as PBE have the smallest overhead associated with them and can be obtained from many widely available software packages. OPTM-GTO orbitals provide a good performance in all the systems we have considered to date but have the additional burden of requiring an optimization framework and careful testing. Once the optimized orbital sets have been developed though, they represent a negligible overhead compared to the use of mean-field orbitals since the evaluation of the periodic GTO orbitals on the FFT grid is a small overhead compared to the cost of evaluating the second-quantized Hamiltonian necessary for AFQMC and many other accurate orbital-based quantum many-body approaches. Out of the three orbital sets considered in this work, MP2-NO are the only ones that represent a significant overhead compared to the generation of the second-quantized Hamiltonian. Our current implementation of the MP2 energy scales roughly as $O[n2m2Nk3Ng]$⁠, where n is the number of electrons, m is the size of the basis set, N_k is the total number of k-points, and N_g is the number of states in the auxiliary basis used to represent the electron repulsion integrals (typically given by the size of the FFT grid within the plane-wave cutoff for the density). The calculation of the MP2-NO is more expensive by roughly a factor of $O[mNk]$⁠, which makes the evaluation of natural orbitals considerably more expensive than the evaluation of the energy. Even the evaluation of approximate natural orbitals for MP2, presented in Ref. 35 as a reasonable alternative to full natural orbitals, represents an increase in the cost of $O[m/n]$ compared to energy evaluations in our implementation. As the system size grows, this becomes significantly more expensive to the point that it becomes slower than the AFQMC calculation itself. This is mainly due to the fact that, in MP2 calculations, we must keep all the states in the PW basis in order to produce high quality natural orbitals, so m is typically an order of magnitude larger in MP2-NO calculations compared to the number of states used in the subsequent AFQMC calculation. Taking this into consideration, we believe that the OPTM-GTO orbitals offer the best balance between performance and computational cost.

All results presented up to this point have been obtained from calculations on the primitive unit cell of the material. In this section, we study the convergence of the equation of state properties, including the lattice constant, bulk modulus, and atomization energies, on the basis set and system size. We only present results for the OPTM-GTO basis in this section since it is the main focus of this work and since the calculation of MP2 natural orbitals represented a big computational cost in larger system sizes (more expensive than the AFQMC calculations themselves for systems with k-point grids beyond 2 × 2 × 2).

Table I shows a summary of the equation of state properties of MgO for various OPTM-GTO basis sets for a gamma-centered 2 × 2 × 2 k-point grid. We show results for uncorrected AFQMC energies as well as for the two MP2-based corrections described above. Properties were obtained from a fit to the Vinet equation of state using the results at seven volumes near equilibrium. Only the correlation energy was evaluated on the OPTM-GTO basis, the exchange energy was calculated separately with Quantum Espresso on a 10 × 10 × 10 k-point grid. Both basis set correction schemes give reasonable results at the pVTZ level already. We observe similar behavior in all three materials. See the supplementary material for a similar table with additional system sizes and for results on the other materials.

We performed calculations on 2 × 2 × 2, 3 × 3 × 3, and 4 × 4 × 4 k-point grids and various OPTM-GTO basis sets. We used the results of the 3 × 3 × 3 and 4 × 4 × 4 k-point grids to perform an extrapolation of the energies to the thermodynamic limit, assuming a 1/N_k dependence at each volume independently. We then combined these finite-size corrections with the scaled MP2-correction for basis set incompleteness to estimate the equation of state in the thermodynamic and complete basis set limits. Table II shows a summary of equation of state properties for all materials. We present results for uncorrected AFQMC on a 3 × 3 × 3 k-point grid with a pVTZ basis, the same calculations with scaled MP2 basis set and finite-size corrections, as well as experimental results corrected for zero-point vibrational effects. Overall, we obtain excellent agreement with experimental results for almost all properties. Even uncorrected AFQMC results are in reasonable agreement with experiment in most cases, which shows the excellent performance of the basis set generation scheme in these materials. We note that the only significant discrepancy with experiment comes from the atomization energy of Si. While several effects can be contributing to this error, including the phaseless approximation in AFQMC, pseudo-potential approximations, and remaining finite-size error, a large fraction of the error most likely comes from an overestimation of the basis set incompleteness error when using the scaled MP2 correction in Si. From Fig. 3, this overestimation is on the order of 0.14 eV in the primitive unit cell for the OPTM-GTO basis set. We are currently investigating alternative basis set correction schemes for systems where MP2 is not appropriate, including RPA-based methods. The results of this study will be presented in a future publication.

We present an efficient approach for the generation of basis sets for correlated calculations in the solid state. By minimizing the MP2 energy of a small representation of the solid, we are able to generate a basis set that offers a convergence of the correlation energy as a function of system size that is comparable to more sophisticated approaches such as the use of MP2 natural orbitals but at a fraction of the cost. The resulting basis sets offer excellent transferability as a function of volume and lead to reasonably converged equation of state properties with modest sizes and much reduced computational costs. We also presented a simple but powerful basis set correction scheme based on MP2, inspired by recent work on coupled cluster theory.⁴⁰ The correction is able to provide a reasonable estimate of the basis set incompleteness error in AFQMC calculations and helps extend the range of systems within the range of the method. We are currently exploring extensions of the work presented here, including the use of contracted basis sets, the use of RPA as the correlated method instead of MP2, and the transferability of the resulting basis sets with different environments and varying chemical composition. We believe that this work can pave the way for more robust and efficient basis set generation schemes in correlated solid state calculations, reducing one of the main hurdles in widespread application of quantum chemistry methods in solids.

See the supplementary material for additional figures showing the non-parallelity errors in Si and LiH, including MP2 and scaled MP2 corrections for both materials. We also provide equation of state properties for MgO, LiH, and Si for additional system sizes and basis sets.

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials (CPSFM). This work was performed under the auspices of the U.S. DOE by LLNL under Contract No. DE-AC52-07NA27344. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. Some AFQMC calculations received computing support from the LLNL Institutional Computing Grand Challenge program.

The data that support the findings of this study are openly available in the Materials Data Facility (MDF).⁹⁰