X marks the spot: Accurate energies from intersecting extrapolations of continuum quantum Monte Carlo data Open Access

The nominal computational cost of cQMC calculations scales polynomially with system size N, typically as N²–N⁴, and the quality of the resulting energies depends on the accuracy of the trial wave function for VMC and on the accurate location of its nodes for DMC. The cQMC methods excel at describing explicit dynamic and long-ranged correlations, but the error incurred by the fixed-node approximation is often significant. By contrast, FCIQMC is formally an exponentially scaling method that trivially captures static correlations but requires a very large number of walkers to provide a good description of dynamic correlations. The complementary nature of the strengths of cQMC and FCIQMC makes combining these methods highly desirable. Several ways of combining cQMC and FCIQMC have been presented in the literature, such as using DMC to assist in the extrapolation to the thermodynamic limit of FCIQMC energies of the electron gas¹¹ or using VMC-optimized Jastrow factors in FCIQMC with the transcorrelated method.^12,13 Here, we shall focus on the use of selected CI wave functions generated with FCIQMC to construct multideterminantal trial wave functions for cQMC calculations.

Multideterminant expansions have been used for decades in cQMC calculations of atomic and molecular systems, including ground-state energy calculations,^14–19 excitation energies,^20–23 and geometry optimizations.^22,24 The use of truncated CI expansions in cQMC presents the problem that no reliable criteria exist to truncate wave functions for different systems in a consistent manner, resulting in energy differences of questionable accuracy. One possible approach is to use extremely large multideterminantal wave functions^19–27 under the expectation that the fixed-node error in the total energies will become smaller than the target error. While algorithmic developments have vastly reduced the computational cost associated with the use of multideterminantal wave functions in cQMC,^25,28–31 this remains an expensive choice. Using trial wave functions without a Jastrow factor reduces the nominal computational burden^20,21,27 at the cost of losing the accurate, compact description of dynamic correlation afforded by fully optimized trial wave functions. By including explicit correlations, in the present paper we are able to explore the use of relatively small multideterminantal wave functions to perform an extrapolation of the cQMC total energy to the full-CI, complete orbital-basis limit. We test our method on a variety of molecular systems, obtaining total and relative energies within uncertainty of benchmark-quality results from the literature.

The rest of this paper is structured as follows: In Sec. II, we present the methodological details of our extrapolation method, which we illustrate with calculations of the carbon dimer and the water molecule. We then apply our method to several atomic and molecular systems, and we report the results in Sec. III. Our conclusions and outlook are presented in Sec. IV. Hartree atomic units (ℏ = |e| = m_e = 4πϵ₀ = 1) are used throughout; the uncertainties and error bars we report refer to standard 68.3% (one-sigma) confidence intervals except when explicitly noted otherwise.

In what follows, we develop the methodology to enable the application of the xspot method in practice using as test systems the C, N, and O atoms, the ground-state C₂, N₂, H₂O, and CO₂ molecules, and the C₂ molecule in its lowest-lying singlet electronic excited state, which we refer to simply as $C 2 *$⁠. These atoms and molecules are simulated as all-electron, both in the sense that no effective-core potentials are used and that excitations from “core” orbitals are allowed in the CI wave function. In Table I, we give the states and geometries we have used for these systems.

The extrapolation shown in Fig. 1 might seem simplistic from a quantum chemical perspective, given that all calculations involved have been performed with the same, finite orbital basis, so one would expect an orbital-basis dependent result, which should itself be extrapolated to the complete-basis limit. For instance, the FCIQMC energy tends to a basis-set dependent FCI limit as the number of walkers tends to infinity, and this must in turn be extrapolated to the basis-set limit in order to obtain the exact energy of the system.

In what follows, we will conceptually combine the choice of molecular orbitals (e.g., Hartree–Fock, natural orbitals, …) with the choice of basis set (e.g., cc-pCVDZ, cc-aug-pVTZ, …), so we shall discuss the completeness of the (molecular) orbital basis instead of that of the basis set alone to emphasize this point.

The xspot extrapolation procedure can be easily justified in the hypothetical case of using an infinite, “complete” orbital basis. The FCIQMC energy with this “complete” orbital basis would tend to the exact total energy of the system E₀ in the infinite walker-number limit, and the sum of the squared CI coefficients would also tend to that of the exact wave function, w₀. The exact wave function has no dynamic correlation left to recover, so the Jastrow factor and backflow displacement in the cQMC trial wave function would optimize to zero, and the VMC and DMC energies would both coincide with E₀. At finite expansion sizes, w < w₀, the VMC and DMC methods yield variational energies satisfying E_VMC ≥ E_DMC ≥ E₀, which, assuming these to be smooth functions of w, validates the xspot method with the “complete” orbital basis.

We note that in a truncated CI wave function, the infinite, “complete” orbital basis is effectively finite since a finite number of determinants can only contain a finite number of distinct orbitals. Conversely, a sufficiently small selected CI wave function with a finite orbital basis is indistinguishable from a CI wave function of the same size with the “complete” orbital basis—assuming the finite basis contains the first few orbitals in the “complete” orbital basis.

As the orbitals in a finite basis get used up, the cQMC energies can be expected to plateau as a function of w as they tend to their orbital-basis dependent limit. We refer to this phenomenon as “orbital-basis exhaustion” and to the onset of this plateau as the exhaustion limit w_exh.. Note that orbital bases such as natural orbitals can be constructed so as to compactly describe the system with fewer orbitals, which has the side effect of reducing the value of w_exh.. We discuss this aspect further in Sec. II C.

As a proxy for the degree of orbital-basis exhaustion, in Table II, we show the fraction of orbitals used in CI wave functions of the same size for Hartree–Fock orbitals expanded in four different basis sets in the cc-pVxZ and cc-pCVxZ families^32,33 for the all-electron carbon dimer. Based on these numbers, we use the cc-pCVTZ basis throughout this paper to ensure we have enough leeway to increase the multideterminant wave function size before hitting the exhaustion limit. We provide an a posteriori assessment of this choice in Sec. III.

Finite-orbital-basis FCIQMC and cQMC calculations performed at w < w_exh. behave as if one were using the “complete” orbital basis. Therefore, it is legitimate to expect that the extrapolation of quadratic fits to these VMC and DMC data intersect at w = w₀ and E = E₀, provided that the VMC and DMC energies are smooth functions of w representable by a second-order polynomial for w_h.o. < w < w_exh., where w_h.o. is a threshold below which higher-order polynomials would be needed.

Note that the initial FCIQMC wave function with M_gen determinants is not required to be below the exhaustion limit since it simply serves to construct selected CI wave functions of size M ≪ M_gen, which are required to be below the exhaustion limit, and to define the arbitrary point at which w = 1 in the plots; w = 1 has no special significance in this method.

In our calculations, we choose CI wave function sizes so that the points are more or less evenly spaced on the w axis, and we make sure that different points correspond to wave functions containing a different number of distinct spatial orbitals so as to capture the effect of simultaneously growing the CI expansion and the orbital basis.

In order to account for the statistical uncertainty and quasi-random fluctuations in the xspot method, we use a Monte Carlo resampling technique in which we generate 100 000 instances of each VMC and DMC dataset in which a random amount proportional to $α 2 + Δ E i 2$ is added to the original energy values. We then perform fits on these shifted data, find the intersection point for each such instance, and obtain the final result by averaging over instances; see Fig. 3 for an illustration of this process. This procedure provides meaningful uncertainties on the intersection energies, which account for both the cQMC statistical uncertainty and quasi-random deviations from the smooth trend.

We demonstrate the full statistical procedure of the xspot method on multideterminant-Jastrow data for the carbon dimer in Fig. 4. Notice that the distribution of intersection points shown in the inset of Fig. 4 has a tail extending toward low E and large w. These tails become more problematic the more parallel the two intersecting curves are, eventually preventing the evaluation of an intersection point at all. It is, therefore, important to try to apply the xspot method to curves that are as close to perpendicular as possible.

In Fig. 5, we include additional VMC and DMC data using an inhomogeneous backflow transformation (“bVMC” and “bDMC”) for the carbon dimer. We list the intersections between pairs of curves in Table III, all of which are within uncertainty of each other. The VMC and bDMC curves provide the best-resolved results, which is to be expected since these curves intersect at the widest angle among all pairs of curves, as shown in Fig. 5. By contrast, the DMC and bVMC curves intersect at a narrow angle and incur a small but non-zero fraction of “missed” intersections, i.e., random instances of the data whose fits fail to intersect at w > 1, which signals the presence of heavy tails in the intersection distribution, resulting in a large uncertainty on the intersection energy.

For the four curves in Fig. 5, the estimated amplitude of the quasi-random fluctuations ranges from α = 0.17 to 0.41 mHa. Throughout this paper, we converge the cQMC calculations so that the uncertainties satisfy $⟨ ( Δ E i ) 2 ⟩ < α 2$⁠, i.e., so that quasi-random fluctuations represent the main contribution to the uncertainty on the fits and intersections; see the supplementary material for the list of values of α obtained. The statistical uncertainties on the cQMC energies can thus be neglected for all practical purposes.

As alluded to in Sec. II A, the choice of molecular orbitals plays a crucial role in the behavior of the cQMC energies as a function of w, modifying the point w_exh. at which the effects of exhaustion start to become noticeable. In Fig. 6, we demonstrate this for the H₂O molecule by comparing the cQMC energies obtained using Hartree–Fock orbitals expanded in the cc-pCVTZ basis set and natural orbitals expanded in the same basis set constructed so as to diagonalize the one-body density matrix in coupled cluster singles and doubles (CCSD). While CCSD natural orbitals produce lower cQMC energies and correspond to larger values of w at fixed expansion sizes, the cQMC energies we obtain using Hartree–Fock orbitals follow the quadratic trend throughout the whole w range considered, while those obtained with natural orbitals plateau very early on, preventing their meaningful extrapolation.

In our calculations, we use Hartree–Fock orbitals expanded in the cc-pCVTZ Gaussian basis set^32,33 obtained using molpro.³⁶ We perform a small-scale FCIQMC calculation using the neci package³⁷ with configuration state functions (CSFs) instead of determinants as walker sites,³⁸ which reduces the number of FCIQMC walkers required to accurately represent the wave function; note that the use of CSFs is not a requirement of the xspot method. The FCIQMC population is grown to 10⁶ walkers and equilibrated, and the coefficients of the M_gen = 10 000 most-occupied CSFs are recorded from a population snapshot. From this information, we build CI expansions with the M CSFs with the largest absolute coefficients, where M ≤ 1500 ≪ M_gen. In our cQMC calculations, the orbitals are corrected to obey the electron–nucleus cusp condition.³⁹ 

The CSF coefficients are reoptimized in the presence of a Jastrow factor of the Drummond–Towler–Needs form^5,6 and an optional inhomogeneous backflow transformation including electron–electron, electron–nucleus, and electron–electron–nucleus terms.⁷ We do not optimize any of the parameters in the molecular orbitals, which provide degrees of freedom that overlap significantly with those in the backflow transformation. Note that even though CSFs are used, the presence of the Jastrow factor and the backflow transformation prevents the cQMC trial wave function from formally being an exact spin state.⁴⁰ We optimize our wave function parameters using linear least-squares energy minimization^3,4 with 10⁶ statistically independent VMC-generated electronic configurations, a number large enough that the optimized cQMC energy can be assumed to lie reasonably close to its variational minimum; note that any remaining optimization error can be considered to be absorbed into the quasi-random error.

The resulting trial wave function is then used to run two DMC calculations with time steps 0.001 and 0.004 a.u. and target populations of 2048 and 512 configurations, respectively, except for bDMC runs on CO₂ at 500 and 1000 CSFs, for which we use 65 536 and 16 384 configurations. These energies are then linearly extrapolated to the zero time-step, infinite-population limit.^2,41

We use the casino package² to run the cQMC calculations and use multi-determinant compression³⁰ to reduce the computational expense of evaluating the trial wave function. We perform the fits to the data and find their intersections using our custom polyfit tool.⁴² The cQMC energies obtained for all systems can be found in the supplementary material.

In this section, we test the xspot method on all eight systems under consideration to assess the different aspects discussed in Sec. II and determine the broader applicability of the method. The VMC and bDMC energies and fits we obtained for the eight systems are shown in Fig. 7; additional plots containing the bVMC and DMC energies can be found in the supplementary material.

All the curves in Fig. 7 are relatively smooth and provide a well-defined intersection. The apparent non-monotonicity of the bDMC curve for the carbon atom is an artifact of the use of a fit function that formally allows non-monotonic behavior, and should be interpreted accordingly. E(w) can be regarded as approaching the intersection with negligible slope, and the region w > w₀ should be ignored since E(w) does not have a physical meaning there. All of the other fits appear to be monotonic in the range shown.

The fraction of orbitals used, a proxy for the degree of orbital-basis exhaustion, is plotted for each of the systems in Fig. 8. We do not use up all of the orbitals in the basis in any of our calculations, and the curves in Fig. 7 do not seem to exhibit symptoms of orbital-basis exhaustion. The bVMC energies do seem to plateau somewhat, which we discuss briefly in the supplementary material; note that we do not use the bVMC data to obtain our final results.

In Table IV, we compare the total energies obtained from applying the xspot method to VMC and bDMC energy data with benchmark-quality estimates of the exact nonrelativistic energies of the systems from the literature, along with our best bDMC result and prior cQMC results for reference. The atomization energies of the ground-state molecules are shown in Table V. For excited-state $C 2 *$⁠, we compare the vertical excitation energy with that calculated with internally contracted multi-reference coupled cluster (ic-MRCC) theory;^44–46 we have computed the total energy of $C 2 *$ shown in Table IV by adding the ic-MRCC excitation energy to the estimated ground-state energy of C₂ from Ref. 34.

All of the total energies reported in Table IV are within statistical uncertainty of their corresponding benchmark values. An important observation is that our individual cQMC energies are not lower than those from prior cQMC calculations, implying that we incur a lower computational cost, but our xspot results are in general closer to the benchmarks than cQMC results from prior studies. Our xspot energies are on average 1.1 standard errors above the benchmark, with a root-mean-square deviation of 1.9 standard errors. These results are compatible with the xspot method being exact when the method’s assumptions are satisfied. The relative energies are likewise in agreement with the benchmark values.

We find that the magnitude α of the quasi-random fluctuations of the cQMC energies is up to 0.7 mHa. These fluctuations are particularly visible in the VMC data for N, O, and H₂O in Fig. 7, for example; in the supplementary material, we give the values of α we have obtained for each of the curves. The magnitude of the quasi-random fluctuations does not seem to increase too rapidly with system size, but their effect on the extrapolated energy becomes more pronounced the further the cQMC data are from the intersection in the plots. This increasing uncertainty on the xspot total energies, reaching 4 mHa for the CO₂ molecule, hints at a limitation of the methodology: the cQMC energies and values of w obtained using modest-sized multideterminant expansions with a fixed basis set to move away from the intersection point with increasing system size, which in turn exacerbates the effects of quasi-random noise on the uncertainty of the xspot energy; one would have to use bigger basis sets and larger multideterminantal expansions to get data closer to the intersection in order to reduce this uncertainty, increasing the computational cost of the approach.

We have presented an empirical extrapolation strategy for cQMC energies as a function of the sum of the squared multideterminant coefficients in the initially selected CI wave function from which the trial wave function is constructed. This approach is made possible by the smoothness of the energies as a function of the CI expansion size, and we have presented a simple statistical procedure to handle the quasi-random non-smoothness in the data, which we show to work very well in practice. We find that Hartree–Fock orbitals expanded in standard basis sets provide the type of gradual convergence required for the xspot method to work well. The results from the tests we have conducted are compatible with the xspot method being capable of obtaining exact total energies, with the caveat that trial wave function complexity must increase with system size in order to control the uncertainty on the results.

See the supplementary material for the cQMC data used in this paper, a table of the magnitude of the quasi-random fluctuations encountered, and a discussion of connected extrapolation approaches. The supplementary material additionally cites Refs. 47–51.

P.L.R. and A.A. acknowledge the support from the European Center of Excellence in Exascale Computing, TREX, funded by the Horizon 2020 program of the European Union under Grant No. 952165. The views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or the European Research Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

The authors have no conflicts to disclose.

Seyed Mohammadreza Hosseini: Investigation (equal); Software (equal); Writing – original draft (equal). Ali Alavi: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Pablo López Ríos: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.