Calculating potential energy curves with fixed-node diffusion Monte Carlo: CO and N₂

Quantum Monte Carlo (QMC) methods have shown promise in performance and scalings as an approach to quantum mechanical calculations.^1–5 These methods have been applied not only to the electronic structure, especially of solids^6–10 and medium-sized molecules,¹¹ but also to atoms and small molecules^12–28 presenting an alternative to traditional high-accuracy quantum chemistry methods such as configuration interaction (CI)²⁹ and coupled-cluster (CC).^30,31 Some typical limitations of standard ab initio methods are as follows: (1) high-order dynamic electron correlation may be neglected or is prohibitively costly to compute; (2) scaling with the number of electrons is poor;^32,33 (3) some errors may be introduced by internal contraction;³⁴ and (4) many algorithms are not yet efficient for large scale parallelization.

Key advantages for QMC methods are favorable scaling with the number of electrons and efficient parallelization (scaling with the number of computing cores). CI and CC methods can scale as n⁷ (for n electrons) or worse³² and thus rapidly become prohibitively expensive with the increasing system size, whereas QMC methods scale as n³ making them especially favorable for larger systems. In practice for small systems, despite the impressive n³ scaling, QMC tends to have a large cost-prefactor making it relatively expensive compared with traditional quantum methods. Nevertheless, the vastly better scaling means that there exists a crossover point in system size beyond which QMC is favored. In addition, QMC algorithms can be very efficiently parallelized, scaling essentially linearly with the number of cores.³⁵ QMC methods can take full advantage of massively parallel machines, and are thus well-suited for next-generation computer architectures with millions of cores.^3,5

During the development of QMC methods over the past few decades, there have been numerous reported studies of first-row atoms and (mostly) homonuclear diatomics^12–28,36 as well as hydrides.^37,38 The majority of those studies have focused on a single equilibrium geometry for each species while treating the total binding energy as well as components of the energy as method development benchmarks. It is well established that QMC methods can capture large fractions of both the strong and dynamic correlation energies, illustrated for example, by its impressive performance for the challenging Be₂ system.^17,25 Methods to compute forces have also been developed³⁹ and estimates of anharmonic force constants based on a few near-equilibrium points have been reported.¹⁵ 

Some PECs calculated with QMC have been reported^40–42 including the most recent work of Giner et al.⁴⁰ Those studies focus on technical aspects of the wavefunction construction and sampling as well as accuracy performance, but with less emphasis on assessments of cost and routine feasibility. They did achieve high-accuracy in comparison with reference curves, and reported that the fixed-node error (due to the fixed-node approximation)⁴³ was reduced as an increasing number of determinants were used. Some of the studies did not pursue the highest degree of accuracy, limiting the size of basis set in the trial wavefunction.⁴⁰ The largest bond dissociation distances were the least well converged relative to reference PECs,^37,40 and this was ascribed to limitations in the number of determinants.

For PECs, single-reference methods (with typical levels of truncation of the excitation operators) often do not produce the correct behavior over the entire bond-distance range, because they cannot easily account for the evolution toward other configurations as the bond is stretched toward dissociation. The problem is common when breaking multiple bonds, a well-known example being N₂.⁴⁴ Even though CCSD(T), i.e., coupled-cluster with single, double, and perturbative triple excitations, is considered the “gold standard” of quantum chemistry, for N₂ it cannot accurately calculate the region of the PEC corresponding to dissociation. CCSD(T) breaks down at about twice the equilibrium bond length r_e and produces an artificial hump.^44,45

Multi-reference methods are often preferred as a globally applicable approach with correct dissociative behavior and have been found necessary in previous QMC studies.³⁷ The single determinant fixed-node (FN-DMC) atomization benchmarks reported by Grossman in 2002 found a 2.9 kcal/mol average absolute deviation over a 55 molecule test set, with deviations of 3.0 and 4.1 kcal/mol for CO and N₂, respectively.¹⁴ The efficient use of multideterminants in QMC codes is not trivial but has been successfully addressed by an algorithm known as the table method^46,47 implemented in the QMCPACK⁴⁸ code package. The work reported here focuses on assessing the accuracy, cost, and routine feasibility of using QMC as implemented in two freely available packages, calculating ground state potential energy curves (PECs) for two test systems (CO and N₂). The multi-determinant trial wave functions used in this study were generated from orbitals and determinants using the multi-configurational self-consistent field (MCSCF) and configuration interaction (CI) methods from GAMESS (U.S.).⁴⁹ The CASINO⁵⁰ package was used for most of the QMC calculations reported here, along with some timing comparisons conducted using the QMCPACK⁴⁸ code.

Here, we give a brief summary of some important aspects of QMC. For more technical details, the reader is referred to Refs. 50 and 51.

Variational Monte Carlo (VMC), using an approximate trial wave function $ΨT$⁠, calculates the expectation value of a Hamiltonian H, with the integrals being performed by a Monte Carlo method.^52–54 The variational energy $EVMC$ is an upper bound to the exact ground state energy $E0$ and is mathematically defined as

where R is a 3N-dimensional vector of the coordinates (r₁,r₂,…,r_N) of the N particles in the system (in this case electrons).⁵⁰ The expectation value of the Hamiltonian H can be rewritten with respect to the trial wave function $ΨT$ as

with the local energy $Eloc=ΨT−1HΨT$⁠. By sampling the points $𝑹i$ according to the distribution $|ΨT|2$ with $M$ configurations, the variational energy can then be computed from a set of local energies,

Statistical uncertainty in a Monte Carlo method decreases as $1/M$⁠, where M is the number of samples. The primary purpose of the VMC method in this application is to optimize the parameters of a trial wave function, such as the Jastrow factor,^55–59 for subsequent use in the more accurate DMC (Diffusion Monte Carlo) method.

In this work, to help describe dynamic electron correlation, a three-body Jastrow factor was used,⁵² which includes electron-electron u terms, electron-nucleus χ terms, and electron-electron-nucleus f terms. The Jastrow factor makes the trial wave function depend explicitly on particle separations and is symmetric under the interchange of identical particles.⁶⁰ Note that in a similar application to the F₂ molecule, Giner et al. justified not employing a Jastrow factor, primarily in order to avoid the costly optimization of parameters.⁴⁰ They state that non-linear optimizations of Jastrow factors may make it more difficult to obtain smooth PECs. The main drawback to lacking a Jastrow factor conceded by Giner et al. is a greatly increased variance in the trial wavefunction and corresponding increased simulation times.⁴⁰ For a detailed description of the form of the Jastrow factor that is used in CASINO, see Ref. 54.

The Diffusion Monte Carlo (DMC) method^61,62 uses the importance-sampled imaginary-time Schrödinger equation to evolve an ensemble of electronic configurations toward the ground state. The Schrödinger equation can be written in integral form as

where R is a point in the configuration space of an N-particle system. Here, $f(𝑹,t)$ is a mixed distribution dependent on some trial wave function $ΨT$⁠, and written mathematically as $f(𝑹,t)=Ψ(𝑹,t)ΨT(𝑹)$⁠. The Green’s function $G(𝑹←𝑹′,t−t′)$ gives the probability of a transition from R to $𝑹′$ in the time interval $t−t′$⁠, and satisfies the initial condition $G(𝑹←𝑹′,0)=δ(𝑹−𝑹′)$⁠.

Due to the fermion sign problem,⁵¹ the DMC method in CASINO adopts the fixed-node approximation, in which the nodes of f are constrained to be the same as those of the trial wave function $ΨT$⁠. The DMC method then produces the lowest energy possible for this nodal surface and can be considered variational with respect to it.⁴¹ Interestingly, for both atoms and diatomic molecules, Giner et al. reported a systematic decrease in the fixed-node error with respect to both the number of determinants retained and the size of the one-particle basis.⁴⁰ 

For all of the QMC calculations, multi-determinant Slater-Jastrow (MD-SJ) trial wave functions were used, which can be written mathematically as

where $ΨT$ is the trial wave function, $eJ$ is the Jastrow factor, $ck$ are the determinant coefficients for the multi-determinant expansion describing static correlation,⁵¹ and $Dk↑$ and $Dk↓$ are the spin-up and spin-down Slater determinants, respectively. For both CO and N₂, the aug-cc-pwCV5Z basis sets⁶³ were truncated to 11s10p8d, i.e., angular momentum functions $≥f$ were removed. The orbitals in the trial wave functions were cusp corrected using the scheme of Ma et al.⁶⁰ 

All DMC calculations were performed with the electrons moving one at a time (electron-by-electron). The time steps for all the systems were chosen so that the acceptance ratio of the proposed moves would be ∼99.5%.

For both systems, the Jastrow factor was defined using an expansion of order 8 for the u terms (⁠$Nu=8$⁠), an expansion also of order 8 for the $χ$ terms (⁠$Nχ=8$⁠), and an order of 4 for the f terms (⁠$Nfen=Nfee=4$⁠), resulting in a total of 149 variable parameters for N₂ and 281 parameters for CO in the two respective Jastrow factors. For all trial wave functions, the truncation order-parameter C, which determines the behavior at the cutoff lengths, was set to a value of 3, providing a local energy that is continuous in configuration space.⁶⁰ For the electron-electron u terms and electron-electron-nucleus f terms, different parameter values were used for the parallel- and the antiparallel-spin electron-pairs, and, for the electron-nucleus χ terms, different parameters were used for spin-up and spin-down electrons.

The parameters of the Jastrow factor were optimized with un-reweighted variance minimization. For each initial optimization cycle, the default cutoff lengths for the u, χ, and f terms were used and the initial Jastrow parameters were set to zero. To see how uncertainties in the optimized Jastrow factors would carry over to the subsequent DMC calculations, two procedures with vastly different costs were tested. Procedure 1 is designed to seek high precision during the optimization stage, thus incurring a high computational cost. Procedure 2 explores what compromises in accuracy are suffered when a less rigorous and time consuming optimization is performed. In Procedure 1, the Jastrow parameters were optimized using 5.0 × 10⁵ configurations for one initial cycle, followed by a system-dependent number of additional optimization cycles. For N₂, the number of additional cycles was 6 and, for CO, the number of additional cycles was 9. In Procedure 2, the Jastrow parameters were optimized using 1.0 × 10⁵ configurations (five times fewer than procedure 1) for 10 cycles with the cost of the optimizations being roughly an order of magnitude fewer CPU-hours than Procedure 1. It is common in studies using multiconfigurational trial wavefunctions to also re-optimize the CI coefficients.²⁵ In this study, this was tested for each system at a few points. Following Jastrow optimization using the VMC variance minimization procedure, additional cycles were performed to re-optimize both the Jastrow parameters and the CI coefficients, using the VMC energy minimization method. For these systems, further optimization using energy minimization did not yield significantly lower VMC energies and also resulted in larger variances and uncertainties. The energy minimization optimizations were less robust (particularly at stretched geometries), requiring additional sampling, thus making the procedure less straightforward and further increasing costs. We therefore elected to use the variance minimized Jastrow parameters and not re-optimize the CI coefficients. In both procedures, the subsequent DMC calculations were performed with target populations of 2500, a minimum of 1.5 × 10⁵ sample points, and time steps of 0.002 a.u. for all geometries. For comparison with DMC, multireference configuration interaction (MRCI) calculations were also done using the MOLPRO⁶⁴ package.

When the ground states of carbon and oxygen atoms combine, the number of molecular states of CO, resolved into C_2v symmetry, is⁶⁵ 

which represents a total of nine states of each of three spin-multiplicities (singlet, triplet, and quintet). Since there are nine singlet-states, dynamic weighting⁶⁶ was used to facilitate robust convergence near the asymptote for the full-valence (10e,8o), state-averaged multi-configurational self-consistent field (SA-MCSCF) calculations. As in the case of N₂, a high-spin (quintet) restricted-open Hartree-Fock (ROHF) calculation was performed before the subsequent singlet DW-SA-MCSCF calculations.

To control the cost of the QMC calculations, only a limited number of determinants were retained (specified by two contending criteria). First, a coefficient cutoff of 0.002 was used at each point, such that all determinants with an absolute weight coefficient value >0.002 were retained. This resulted in a varying number of retained determinants at each $rCO$ distance. Second, for comparison, a fixed number of 250 determinants (those with the largest coefficients) were retained throughout the coordinate range.

To assess the quality of the QMC results for CO in comparison to typical high-level conventional electronic structure methods, an accurate MRCI-based reference PEC by Dawes et al.⁶⁷ was used. It was constructed from Davidson-corrected MRCI data, based on a dynamically weighted state-averaged CASSCF reference (DW-SA-CASSCF), with a full-valence active space. The MRCI(Q) data were extrapolated to the complete basis set (CBS) limit with the aug-cc-pwCVnZ (n = 3-5) bases and all electrons correlated. The vibrational levels on the PEC are of spectroscopic accuracy. The MRCI-based PEC also includes small spin-orbit (SO) and scalar-relativistic (SR) corrections. The small SO and SR correction terms were removed from the MRCI-based PES for the comparisons with DMC presented here. This permits a more direct comparison since the DMC Hamiltonian does not include those terms.

When two ground state N(⁴S_u) nitrogen atoms combine, the total number of molecular states of N₂, resolved into D_2h symmetry, is ^7,5,3,1Ag (one state of each of four spin-multiplicities).⁶⁵ 

For the calculations in GAMESS, to ensure convergence to the ground state of N₂, a high-spin (in this case septet) restricted-open Hartree-Fock (ROHF) calculation was performed at the largest N–N distance, followed by a 1-state MCSCF calculation for the desired singlet ground state with the full-valence (10e,8o) active space. In the QMC calculations, based on the results for CO discussed above, the strategy of retaining a variable number of determinants based on a coefficient threshold was abandoned, and a fixed (generous) number of 396 determinants were retained. The QMC data obtained along the coordinate range were assessed by comparison with an empirical spectroscopic PEC by Le Roy et al.⁶⁸ 

For CO, initially a coefficient weight cutoff of 0.002 was used to restrict the number of determinants retained in the trial wavefunctions, generally resulting in a different number of determinants employed at each CO bond distance. The DMC energy data and the number of determinants are given in Tables I and II for Procedure 1 and Procedure 2 (respectively), which differ by roughly an order of magnitude in the computational expense devoted to optimizing the Jastrow factor (the parameters are better converged by Procedure 1). The DMC energy data were overlaid with the reference PEC using weighted least-squares to account for the different uncertainties at each data point. The root mean-squared deviation (RMSD) of the DMC energy data with respect to the reference PEC is 315 cm⁻¹. As shown in Figure 1, the points obtained via both (Jastrow optimization) procedures follow the reference PEC quite well near the equilibrium geometry. At bond distances >2.0 Å, where fewer determinants are retained, the DMC points tend to deviate slightly from the reference PEC (to higher energies). To see if this behavior might be related to the reduced numbers of determinants, trial wave functions with a fixed number of determinants (250) were optimized via both procedures and then recalculated with DMC. The results are given in Tables III and IV and in Figure 2. The RMSD of the DMC energy data computed using a fixed number of determinants with respect to the reference PEC is 292 cm⁻¹. This is only slightly less than that of the variable determinants data. Note that corresponding to the number of DMC samples, the uncertainties at each data point are in the range of 40-90 cm⁻¹ (see Tables). It is not clear from these results that retaining a fixed number of determinants produces significantly more consistent results than the lower-cost coefficient weight cutoff strategy (which results in a variable number of retained determinants). Regarding the Jastrow optimization, since the results obtained by procedures 1 and 2 are of similar quality and are essentially interchangeable, it appears that the investment of considerably more CPU time in Procedure 1 was not warranted.

To test whether sensible spectroscopic parameters could be obtained from the data, a fit to the Morse function was performed (it is straightforward to convert Morse parameters into anharmonic parameters).¹⁵ First, since the most accurate PEC will not be precisely Morse-like, parameters obtained by fitting the reference PEC were compared to the experimental CO parameters of r_e = 1.128 Å, $ωe=2170.2$ cm⁻¹, and $ωexe=13.46$ cm⁻¹.⁶⁵ The Morse form could not accurately accommodate the entire coordinate range of the reference PEC and produced a fitted value for the harmonic constant that is significantly too high (⁠$ωe=2215.1$ cm⁻¹). Fitting data in a more limited range of r_CO = [0.9, 1.6] Å produced more reasonable values of r_e = 1.128 Å, $ωe=2174.9$ cm⁻¹, and $ωexe=13.05$ cm⁻¹. Fitting the DMC data over the same coordinate range yields parameters of r_e = 1.129 Å, $ωe=2187.0$ cm⁻¹, and $ωexe=13.28$ cm⁻¹, all of which are quite useful estimates. Since the spectroscopically accurate reference PEC and the DMC data are globally consistent, it seems likely that further reducing the DMC uncertainties (through additional sampling) and a denser grid would yield even better constants.

For N₂, a fixed number of 396 determinants were employed for all geometries. The data for Jastrow optimizations via Procedure 1 and Procedure 2 are given in Tables V and VI, respectively, and are plotted in Figure 3. The data are compared with an empirical spectroscopic PEC by Le Roy et al.⁶⁸ The DMC data are generally consistent with the empirical PEC and are much closer to it than moderately high-level MRCI/AVTZ data shown for additional comparison. The RMSD of the DMC data with respect to the empirical PEC is 598 cm⁻¹, which is significantly larger than was found for CO discussed above. Since the empirical PEC for N₂ was obtained via direct fit to spectroscopic data, it implicitly includes small effects such as relativistic corrections not included in the DMC Hamiltonian (that were removed for this reason from the CO PEC for the previous comparison discussed above). However, this is not likely to be a significant source of discrepancy as these sorts of corrections are relatively small for N₂. It is noteworthy that the DMC data point at r = 0.8 Å (listed in Tables V and VI) was excluded from the comparison as it is an outlier roughly 5000 cm⁻¹ more stable than the value given by the reference PEC. The empirical PEC is expected to be unreliable for repulsive geometries far beyond the turning points of the contributing spectroscopic data. Indeed, the empirical PEC was confirmed to be much more repulsive at short range than a high quality ab initio PEC by Gdanitz which is consistent with the DMC value.⁶⁹ Again, as was found for CO, the additional cost of a stricter Jastrow optimization procedure did not yield obviously improved results.

As was done for CO, the reference and DMC-based PECs were each fit to a Morse function to extract spectroscopic parameters. The experimental values used for comparison are r_e = 1.094 Å, $ωe=2359.6$ cm⁻¹, and $ωexe=14.46$ cm⁻¹.⁶⁵ Fitting the reference PEC over the range r = [0.9, 1.3] Å produced values of r_e = 1.098 Å, $ωe=2360.6$ cm⁻¹, and $ωexe=14.70$ cm⁻¹. Fitting the DMC data over the same coordinate range yields parameters of r_e = 1.097 Å, $ωe=2373.4$ cm⁻¹, and $ωexe=14.45$ cm⁻¹, all of which are again quite accurate.

The cost of generating high quality energies via DMC for the two 14-electron systems (CO and N₂) was assessed. A significant fraction of the total cost relates to the optimization of a Jastrow factor that was employed in this study. As mentioned previously, some advantages and disadvantages of using a Jastrow have been noted by Giner et al.⁴⁰ The initial Jastrow optimizations for N₂ with a fixed number of determinants using cheaper Procedure 2 took ∼800-1300 CPU-h per point. With CO, for Jastrow optimizations via Procedure 2 with a varying number of determinants (using a coefficient cutoff strategy), the CPU-hours depend significantly on the number of retained determinants. The lowest cost was the bond distance of 3.5 Å, which includes 87 determinants and took about 300 CPU-h. The highest cost was 1.5 Å, which includes 317 determinants and took about 1860 CPU-h. For Procedure 1, the cost for CO increased to ∼5000-15 000 CPU-h per point, and the cost for N₂ increased to ∼18 000-25 000 CPU-h per point.

Once the Jastrow has been optimized by VMC via either Procedure 1 or Procedure 2, the subsequent DMC energy calculations add a significant additional cost which depends on the number of DMC samples which in turn determines the final uncertainties (related as $1/M$⁠, where $M$ is the number of samples). For N₂, to reach uncertainties on the order of 100 cm⁻¹ (see Tables V and VI), the DMC cost was ∼1800-2200 CPU-h per point. Specifically, after Jastrow optimization via Procedure 1, the cost of the DMC sampling for the N₂ points located at r = 0.8 Å, 1.1 Å, and 2.5 Å took 2170 CPU-h, 1850 CPU-h, and 2010 CPU-h, respectively, to reach an average uncertainty in the three points of ∼86 cm⁻¹. Similarly, after Procedure 2 optimization, the DMC cost for the same points was about 1880 CPU-h, 2035 CPU-h, and 2030 CPU-h, respectively, for an average uncertainty of ∼93 cm⁻¹.

For CO, the DMC sampling cost after Jastrow optimization was ∼775-1900 CPU-h per point to reach similar uncertainties. The cost of the DMC calculations for CO and N₂ (both 14-electron systems) was approximately the same, but the Jastrow optimizations via Procedure 1 took about an order of magnitude longer than those of Procedure 2. The quality of the final energies is similar for the two procedures, indicating that a more conservative Jastrow optimization is reasonable since it appears that diminishing returns are realized with respect to further optimization. Trial wave functions for CO with 250 (fixed) determinants were also optimized using Procedure 2 with a cost of ∼1000-1500 CPU-h per point. To achieve an average uncertainty ∼33 cm⁻¹ in the DMC, an additional cost of ∼7000-10 000 CPU-h per point was required.

The QMCPACK code was used to test the improved efficiency that is expected by the use of the table method algorithm for multideterminant calculations.^46,47 For these comparisons, the same numbers of determinants were retained and the same Jastrow factors were employed. The QMCPACK code was found to reach similar uncertainties for our two test systems in about 30% fewer CPU-hours. For much larger numbers of determinants (up to 16 000), Clark et al. reported much more significant speedups in the range of factors of 15-40.⁴⁷ They noted smaller speedup factors for smaller numbers of determinants. The speedup that one might expect from the QMCPACK algorithm will depend on the size of the system and the number of determinants as well as hardware limitations such as memory and cache.⁷⁰ In our study of two small systems with only 14 electrons and modest numbers of determinants, the speedup is already significant indicating that this method should be preferred in future larger scale multideterminant applications.

It was determined that straightforward application of QMC methods implemented in two freely available codes (CASINO and QMCPACK) could robustly compute electronic energies along dissociation coordinates of small molecules that are comparable in accuracy to high level traditional quantum chemistry. A QMC tutorial aimed at graduate students who have some familiarity with traditional quantum chemistry, but no experience with QMC, is provided as supplementary material.

Points along the potential energy curves of the ground states of CO and N₂ were generated with multideterminant fixed-node diffusion Monte Carlo methods and were found to be in close agreement with spectroscopically accurate curves. The spectroscopic constants obtained by fitting the data are in close agreement with the experiment. In particular, the equilibrium bond distance parameters obtained by fitting to only a few points are within 1-3 pm of experiment.

For the two 14-electron test systems, generating comparably high-quality electronic structure data by conventional methods such as MRCI takes less than 0.5 CPU-h per point, compared with at least 4000 CPU-h for the employed DMC method, as implemented in CASINO, depending on the desired final uncertainty. The QMCPACK code is known to be more efficient for multideterminant trial wavefunctions. For the small test systems and modest numbers of determinants employed in this study, the QMCPACK code was only slightly faster (∼30%), but is indicated for larger scale applications where more significant speedups have been reported.⁴⁷ The favorable n³ scaling of DMC does ensure a crossover point in system size beyond which it becomes cheaper than conventional high-accuracy electronic structure methods (scaling as n⁷ or worse). The large cost pre-factor of DMC seems to preclude it from the routine use in the construction of global PESs (which for 3-5 atom systems typically require thousands of points) at this time. However, in addition to the favorable scaling with system size, QMC methods scale nearly linearly with the number of cores, which could lead to short time-to-solution using next generation architectures with millions of cores.

Even now, given the high accuracy that is achievable via QMC methods, we also see it as a possible arbiter in difficult cases where high-level conventional methods might disagree about the presence or height of a rate determining reaction barrier.^71,72

It is anticipated that QMC methods will become increasingly relevant in the near future.

See supplementary material for a guide illustrating the use of CASINO to perform VMC and DMC calculations with trial wavefunctions generated in GAMESS. The guide is intended to be accessible to readers with some experience in traditional quantum chemistry.

This work was supported by the U.S. Department of Energy [Grant No. DE-SC0010616]. R.D. thanks Paul Kent for useful discussions. A.P. thanks Albert DeFusco and Mike Towler for help with interfacing GAMESS with CASINO.