*Ab initio* quantum Monte Carlo (QMC) methods are a state-of-the-art computational approach to obtaining highly accurate many-body wave functions. Although QMC methods are widely used in physics and chemistry to compute ground-state energies, calculation of atomic forces is still under technical/algorithmic development. Very recently, force evaluation has started to become of paramount importance for the generation of machine-learning force-field potentials. Nevertheless, there is no consensus regarding whether an efficient algorithm is available for the QMC force evaluation, namely, one that scales well with the number of electrons and the atomic numbers. In this study, we benchmark the accuracy of all-electron variational Monte Carlo (VMC) and lattice-regularized diffusion Monte Carlo (LRDMC) forces for various mono- and heteronuclear dimers (1 ≤ *Z* ≤ 35, where *Z* is the atomic number). The VMC and LRDMC forces were calculated with and without the so-called space-warp coordinate transformation (SWCT) and appropriate regularization techniques to remove the infinite variance problem. The LRDMC forces were computed with the Reynolds (RE) and variational-drift (VD) approximations. The potential energy surfaces obtained from the LRDMC energies give equilibrium bond lengths (*r*_{eq}) and harmonic frequencies (*ω*) very close to the experimental values for all dimers, improving the corresponding VMC results. The LRDMC forces with the RE approximation improve the VMC forces, implying that it is worth computing the DMC forces beyond VMC despite the higher computational cost. The LRDMC forces with the VD approximations also show improvement, which unfortunately comes at a much higher computational cost in all-electron calculations. We find that the ratio of computational costs between QMC energy and forces scales as *Z*^{∼2.5} without the SWCT. In contrast, the application of the SWCT makes the ratio *independent* of *Z*. As such, the accessible QMC system size is not affected by the evaluation of ionic forces but governed by the same scaling as the total energy one.

## I. INTRODUCTION

The quantum Monte Carlo (QMC) technique^{1} is one of the state-of-the-art numerical methods for evaluating the expectation values of many-body wave functions, and it provides extremely accurate energetics. There are two major real-space QMC frameworks, specifically the variational Monte Carlo (VMC) and fixed-node diffusion Monte Carlo (FNDMC) methods.^{1} In the VMC framework, expectation values such as energy are calculated using a given trial wave function and a Metropolis algorithm.^{2,3} Once the energy has been calculated for a given trial function, one can readily optimize the variational parameters included in the trial wave function using an appropriate method to minimize the energy and obtain the best possible representation of the ground-state wave function. One of the drawbacks of the VMC method is its strong dependence on the quality of the trial wave function. The FNDMC framework is a stochastic method that filters out the ground state from a given trial wave function using the so-called projection technique or the power method. Since electrons are fermionic, to avoid the well-known sign-problem instability, the projection to the ground state is restricted to leave the nodal surface of the initial trial wave function unchanged, and this is usually given by the best available VMC variational guess. Despite this approximation, the FNDMC method provides results that are usually of much higher quality than the VMC results. TurboRVB, the QMC package used in this work, implements the lattice-regularized version of the FNDMC method,^{4,5} which is hereafter referred to as the lattice-regularized diffusion Monte Carlo (LRDMC) method.

To date, QMC methods have been widely applied to various materials that density functional theory (DFT) fails to cope with, such as molecular crystals,^{6,7} two-dimensional materials,^{8–10} superconductors,^{11} and materials at extreme pressures.^{12–14} After several successful QMC applications, the field is still undergoing rapid development, and several improvements are expected from the methodological perspective. One of the major drawbacks that hinders QMC applications is the lack of a method to compute exact atomic forces. Here, “exact” refers to a force that is consistent with the derivative of the potential energy surface (PES), which is not always true in the QMC framework, as will be discussed later. In the VMC framework, the modern energy minimization^{15,16} and regularization^{12} techniques provide exact force evaluations within a reasonable computational effort. In contrast, the situation is not satisfactory in the DMC framework.

Many researchers have proposed schemes for exact and approximate DMC force calculations, including Reynolds *et al.*,^{17} Badinski and Needs,^{18} Assaraf and Caffarel,^{19,20} Chiesa *et al.*,^{21} Filippi and Umrigar,^{22} and Moroni *et al.*^{23,24} Some of these approaches guarantee an exact force evaluation (full consistency between the PES and the forces) but become impractical as the system size increases. Other schemes are only approximately consistent but allow for the computation of forces with an affordable computational cost. Therefore, to exploit QMC forces in realistic simulations, it is very important to systematically study the reliability of the approximated force evaluations. In this study, we computed the PESs and atomic forces of six dimers—H_{2}, LiH, Li_{2}, CO, MgO, and SiH—at the VMC and LRDMC levels and investigated their corresponding consistency. The LRDMC forces were computed using the Reynolds (RE) approximation^{17} and the variational-drift–diffusion (VD) approximation.^{23}

Another important concern in QMC calculations is the scaling of the computational cost with the system size (*N*) and the atomic number (*Z*). On the one hand, one of the advantages of the QMC framework over other accurate quantum chemistry methods is the scaling with respect to *N*. For example, the computational cost of the coupled-cluster method with single, double, and perturbative triple excitations, or CCSD(T), which is considered to be the gold standard in quantum chemistry, grows as *N*^{7}, while that of DMC methods grows as *N*^{4} at most. On the other hand, in the QMC framework, we should consider the fact that the computational cost of QMC with respect to *Z* is much worse than to *N*. The accelerated Metropolis algorithm^{25} and the real-space double-grid method^{26} have been devised to alleviate this problem. The scaling with *Z* is, however, still higher than that with *N*; for example, it is *Z*^{5.0–6.5} for DMC energy calculations.^{26–29} Last but not least, Tiihonen *et al.*^{30} showed that the scaling of QMC atomic force calculations with pseudopotentials is worse than that of energy calculations, remarkably by orders of magnitude. In detail, they claimed that, at constant computational cost and constant statistical resolution, the accessible system size scales as $Zeff\u22122$, where *Z*_{eff} is the effective valence charge. This result is particularly disappointing in the QMC community since one of the most important properties in computational materials science, the atomic force, is not available for heavier atoms such as transition metals, where QMC methods have proven to be very useful in energetics.^{31–45}

In this study, we revisited and estimated the computational costs of the QMC energy and force calculations with and without the space-warp coordinate transformation (SWCT)^{46} using eight homonuclear atoms: H_{2}, Li_{2}, N_{2}, F_{2}, P_{2}, S_{2}, Cl_{2}, and Br_{2}. We found that the ratio of the computational costs between the QMC energies and forces scales as *Z*^{∼2.5} without the SWCT, while the application of the SWCT makes the ratio *independent* of *Z*. Indeed, the accessible system size is not affected by QMC force calculations when the SWCT variance-reduction technique is applied.

## II. IMPLEMENTATION

As noted above, TurboRVB implements two types of real-space QMC methods: VMC and LRDMC. In both cases, the expectation value of the energy evaluated for a given trial wave function Ψ and a complementary sampling wave function Φ is formally written as

where $x=r1\sigma 1,r2\sigma 2,\u2026,rN\sigma N$ here and henceforth represents the *N* electron coordinates and their spins, while

are the so-called local energy and the probability of configuration **x**, respectively. In VMC, Ψ = Φ, while in DMC, Φ is the projected fixed-node ground state, which is obtained through a stochastic implementation of the power method; thus, it is no longer explicitly known in an explicit form.

Let **R**_{α} be the atomic position of a nucleus *α*. The atomic force acting on *α* is represented as the negative gradient of the energy with respect to **R**_{α},

### A. Variational Monte Carlo forces

As mentioned already, the complementary sampling wave function Φ is equal to Ψ in the VMC framework. One usually ignores Eq. (2c) when evaluating a force, resulting in

This is the conventionally adopted expression for the VMC forces. Neglecting Eq. (2c) is justified only when the system is at a variational minimum for all the parameters (i.e., $\u2202E\u2202ci=0$, *∀i*) or when the variational parameters, which are *implicitly* dependent on the atomic position, do not vary with the atomic position (i.e., $dcidR\alpha =0$; *∀i*); otherwise, a VMC force obtained by using Eq. (3) can be biased. Indeed, the exact force obtained from the slope of the PES is not consistent with the atomic force obtained by the above expression. This error is referred to as the “self-consistency error.”

To obtain a statistically meaningful value of the VMC force with finite variance, the so-called “reweighting” technique is needed because the HF and Pulay terms may diverge when the minimum electron–nucleus distance vanishes and when an electronic configuration is close to the nodal surface. The first type of divergence can be cured by using the method described in Refs. 19 and 48 or by applying the SWCT algorithm,^{46,49,50} whereas the second type of divergence can be alleviated by modifying the VMC sampling distribution using a modified trial wave function that differs from the original trial wave function only in the vicinity of the nodal surface.^{12} By defining *d* as the distance from the nodal region, the local energy (*E*_{L}) and the logarithmic derivative of the wave function ($\u2202\u2061log\Psi \u2202R\alpha $) diverge as ∼1/*d*, leading to an infinite variance problem in Eq. (3). The reweighting method proposed by Attaccalite and Sorella,^{12} which we dub the AS regularization, can cancel out the divergence without introducing bias, realizing finite variance evaluation of VMC forces.

### B. Diffusion Monte Carlo forces

The situation is more complicated in the DMC framework because of the complementary sampling wave function Φ. Indeed, within the DMC framework, there is no straightforward way to evaluate the logarithmic derivative of the complementary sampling wave function ($\u2202\u2061log\Phi \u2202R\alpha $). Therefore, to overcome this difficulty, Reynolds *et al.*^{17} proposed the approximation $\u2202\u2061log\Phi \u2202R\alpha \u2248\u2202\u2061log\Psi \u2202R\alpha $, leading to

The so-called “hybrid estimator” technique can be further used to correct the bias in the Reynolds force due to the above approximation,^{48}

which dubbed the “hybrid force.”

Later on, a more sophisticated approximation, the VD approximation, was proposed by Filippi and Umrigar,^{22} and this was recently revisited by Moroni *et al.*^{23} This approximation replaces the logarithmic derivative of the drift–diffusion part of the Green function along the random walk with the logarithmic derivative of the trial wave function at the current configuration and the sum of the derivative of the local energy along with the previous *k* steps. Indeed, the VD force is evaluated as

where *τ* is the time step defining the discretized times *t*_{n} = *nτ* when the local energy $ELxn$ is evaluated. This formulation was originally proposed within the standard DMC framework,^{23} but it can also be used in the LRDMC framework. Note also that, due to time discretization, a bias proportional to *O*(*τ*^{2}) is also present within LRDMC, which otherwise has only finite lattice mesh error. We have chosen in all forthcoming calculations a time step *τ* = 0.1 a.u., small enough to have negligible errors of this type. Moreover, *k* is a hyperparameter, whose convergence was checked for each dimer.

Here, we note that there are two sources of the self-consistency error in LRDMC calculations: One is the same one as discussed in the VMC calculation and the other comes from the aforementioned approximation (e.g., Reynolds). In principle, it is not possible to discuss the two errors separately because a better trial WF considerably alleviates both errors as shown later in the Li_{2} dimer case.

The AS regularization has also been used for DMC force calculations,^{23,51} but it is not an optimal regularization for this purpose because it enforces a finite density of walkers on the nodal surface.^{24} In this study, we employed the novel regularization technique recently proposed by Pathak and Wagner^{52} for both the RE and VD forces, which we shall refer to as the PW regularization. The PW regularization was originally proposed for computing the derivative of the wave function in the VMC optimization,^{52} but it is also applicable to DMC force calculations, as proposed in Ref. 24. The PW regularization multiplies the terms in the brackets of Eqs. (4) and (6) by the polynomial reweighting function *f*_{ϵ} = 7|*d*/*ϵ*|^{6} − 15|*d*/*ϵ*|^{4} + 9|*d*/*ϵ*|^{2} when |*d*/*ϵ*| < 1, where *d* is the distance from the nodal surface, defined by $d=\Psi (x\u20d7)/\Vert \u2207\Psi (x\u20d7)\Vert $.

### C. Space-warp coordinate transformation

We employed the SWCT^{46,49} in both the VMC and LRDMC force calculations. The SWCT, which was originally proposed by Umrigar^{46} in 1989, is used to mimic the displacement of charges around the nucleus,

where $\kappa r$ is a function that decays for large *r* and *ω* → 0 as a consequence of this property. The original^{46} simple choice $\kappa r=1/r4$ is adopted in TurboRVB. The atomic force is evaluated with the SWCT,^{49}

for the VMC (Φ = Ψ_{T}) and RE-LRDMC forces and

for the VD-LRDMC force, where the brackets indicate a Monte Carlo average over the trial wave function and $rij$ are the electronic coordinates corresponding to **x**_{j}. All the terms above can be very efficiently evaluated by the adjoint algorithmic differentiation^{49} in TurboRVB.^{53} Some of the adjoint routines in TurboRVB were generated from the corresponding forward ones using TAPENADE.^{54}

## III. COMPUTATIONAL DETAILS

The all-electron VMC and LRDMC calculations for the selected dimers (H_{2}, LiH, Li_{2}, N_{2}, CO, F_{2}, MgO, SiH, P_{2}, S_{2}, Cl_{2}, and Br_{2}) were performed using TurboRVB,^{55} and this was combined with the Python package Turbo-Genius^{55} to avoid any human error and to achieve efficient calculations. TurboRVB employs the Jastrow antisymmetrized geminal power (JAGP)^{56} *Ansatz*. The *Ansatz* is composed of a Jastrow and an antisymmetric part (Ψ = *J* · Φ_{AGP}). The singlet antisymmetric part is denoted as the singlet antisymmetrized geminal power (AGPs), which reads

where $A\u0302$ is the antisymmetrization operator and $\Phi r\u2191,r\u2193$ is called the pairing function. The spatial part of the geminal function is expanded over the Gaussian-type atomic orbitals (GTOs),

where *ψ*_{a,l} and *ψ*_{b,m} are primitive GTOs, their indices *l* and m refer to different orbitals centered on atoms *a* and *b*, and *i* and *j* are coordinates of spin-up and spin-down electrons, respectively. When the JAGPs is expanded over *p* = *N*_{el}/2 molecular orbitals, the JAGPs coincides with the Jastrow–Slater determinant (JSD) *Ansatz*.^{5,57} The JSD with DFT orbitals is referred to as JDFT hereafter. Indeed, the coefficients in the pairing functions (i.e., variational parameters in the antisymmetric part) are obtained by the built-in DFT code, named prep, and the coefficients are fixed during the VMC optimization step. This is the simplest choice, but it is reasonable in the majority of QMC applications. The coefficients are further variationally optimized with the JSD and JAGPs *Ansätze* at the VMC level, leading to better nodal surfaces and variational energies at the DMC level.^{58} In this study, we employed the JDFT *Ansatz* except for the case of the Li_{2} dimer. The JDFT, JSD, and JAGPs *Ansätze* were employed for the Li_{2} dimer to further investigate the effect of orbital optimization on the QMC energies and forces. The Jastrow term is composed of one-body, two-body, and three-/four-body factors (*J* = *J*_{1}*J*_{2}*J*_{3/4}). The one-body and two-body factors are used to satisfy the electron–ion and electron–electron cusp conditions, respectively, and the three-/four-body factor is adopted to take the further electron–electron correlation into consideration. The one-body Jastrow factor reads as follows:

where

**r**_{i} are the electron positions, **R**_{I} are the atomic positions with the corresponding atomic number *Z*_{I}, *l* runs over atomic orbitals $\chi I,lJ$ centered on atom *I*, and $ur$ is a short-range function containing a variational parameter *b*,

The two-body Jastrow factor is defined by a long-range function as

where $vr$ is given as

and *F* is a variational parameter. The three- and four-body Jastrow factors are

and

where the indices *l* and m again indicate different orbitals centered on corresponding atoms *a* and *b*. Only three-body (*a* = *b*) factors were used in this study.

We employed *modified* cc-pVTZ basis sets taken from the EMSL basis set library^{59} for the antisymmetric parts. Here, “modified” refers to the modification of the original basis set such that the s orbitals whose exponents are larger than 8 · *Z*^{2} (where *Z* is the atomic number) are disregarded, and they are implicitly compensated by the homogeneous one-body Jastrow part [$J\u03031r$] to fulfill the electron–ion cusp condition explicitly, even at the DFT level.^{58} Indeed, a single-particle orbital is modified as

where $J\u03031r$ is the same as in Eq. (14) and the parameter *b* in Eq. (15) is optimized by direct minimization with the chosen DFT functional. In this way, each element of the modified basis set satisfies the electron–ion cusp conditions even at the DFT level. For the inhomogeneous one-body and three-body Jastrow parts, we employed the cc-pVDZ basis set taken from the EMSL basis set library^{59} and neglected the large exponents according to the same criteria as in the antisymmetric part. The variational parameters in the antisymmetric and Jastrow parts were optimized by the modified linear method^{15,55} implemented in TurboRVB. The exponents of the Jastrow and determinant basis sets were also optimized at the VMC level with the JSD and JAGPs *Ansätze*. The LRDMC calculations for computing the PESs were performed by the single-grid scheme^{4} with lattice spaces *a* = 0.30, 0.10, 0.10, 0.05, 0.05, and 0.03 bohr for H_{2}, LiH, Li_{2}, CO, MgO, and SiH, respectively. The double-grid scheme is more efficient in all-electron LRDMC calculations,^{4,26} but we chose the single-grid scheme in this work to make it easier to compare our computational scaling with other benchmark tests performed within the standard DMC framework. We chose 0.001, 0.001, 0.001, 0.01, 0.02, and 0.01 bohr for H_{2}, LiH, Li_{2}, CO, MgO, and SiH, respectively, for the PW regularization parameter (*ϵ*). These are compatible with the values proposed in the original paper^{52} (∼10^{−2}). It turned out that the choice of *ϵ* did not introduce bias (less than 0.5% error in the DMC force evaluations). The LRDMC calculations for investigating the computational costs of the QMC energies and forces were also performed by the single-grid scheme^{4} with lattice spaces *a* = 1/(2.0 · *Z*) bohr for H_{2}, Li_{2}, N_{2}, F_{2}, P_{2}, S_{2}, Cl_{2}, and Br_{2}. The regularization parameter *ϵ* was set to 0.001 bohr for all the dimers.

## IV. RESULTS AND DISCUSSION

### A. PES calculations and self-consistency errors with JDFT *Ansatz*

Figure 1 shows the PESs of the CO dimer obtained from the VMC calculations and the LRDMC calculations with the RE approximation, wherein the JDFT *Ansatz* was employed. The derivatives of the PESs (i.e., exact forces) and the obtained atomic forces are also shown in Fig. 1. The broken vertical lines in Fig. 1 indicate the equilibrium bond lengths obtained from the PESs, the atomic forces, and experiments.^{60} A comparison of the LRDMC calculations with the VD approximation is shown in the supplementary material. Figure 2 focuses on the region near the experimental equilibrium distance. We can see at a glance that the left (C) force is exactly opposite to the right (O) force at both the VMC and LRDMC levels. We can readily establish that this is because the SWCT should give the same values of forces with the opposite signs in a dimer as a consequence of the following rigorous and straightforward derivation. The local energy and the logarithm of the wave function, which are explicitly dependent on the atomic and electron positions, satisfy the following translational symmetries:

and

They are equivalent to, e.g., by differentiating the right-hand sides of Eqs. (21) and (22) with respect to Δ,

and

Then, as a consequence of $\u2211\alpha \omega \alpha ri=1$, the sum of the HF forces in Eq. (9) over all the atoms in a system is

and the sum of the Pulay forces is

These relations clearly also hold in Eq. (10). Thus, the total force in a system is always zero regardless of the existence of statistical errors. Indeed, atomic forces in an isolated atom are always zero, and if the SWCT is employed correctly, the right and left forces in a dimer should have the same values with opposite signs, regardless of the presence of statistical errors. We also stress that since the application of the SWCT does not change the mean values but only reduces variances at the VMC level, the left and right VMC forces should be consistent within their statistical uncertainties even if the SWCT is not applied, as shown in Tables I and II.

Li_{2}
. | Left force (hatree/Å) . | Right force (hatree/Å) . |
---|---|---|

w/SWCT | −0.027 278(24) | +0.027 278(24) |

wo/SWCT | −0.027 243(58) | +0.027 254(62) |

Absolute difference | 0.000 035(63) | 0.000 024(66) |

Li_{2}
. | Left force (hatree/Å) . | Right force (hatree/Å) . |
---|---|---|

w/SWCT | −0.027 278(24) | +0.027 278(24) |

wo/SWCT | −0.027 243(58) | +0.027 254(62) |

Absolute difference | 0.000 035(63) | 0.000 024(66) |

SiH . | Left force (hartree/Å) . | Right force (hartree/Å) . |
---|---|---|

w/SWCT | −0.224 59(16) | +0.224 59(16) |

wo/SWCT | −0.218 1(53) | +0.224 757(92) |

Absolute difference | 0.006 5(53) | 0.000 17(19) |

SiH . | Left force (hartree/Å) . | Right force (hartree/Å) . |
---|---|---|

w/SWCT | −0.224 59(16) | +0.224 59(16) |

wo/SWCT | −0.218 1(53) | +0.224 757(92) |

Absolute difference | 0.006 5(53) | 0.000 17(19) |

In Fig. 1, the *r*_{eq} value obtained from the LRDMC PES is the closest to the experimental value. This is also true in the other dimers, as will be shown later. Figure 1 also indicates that self-consistency errors exist at both the VMC and LRDMC levels since the orbitals are not variationally optimized. In the CO dimer, the VD approximation alleviates the self-consistency error more than the RE approximation. This is not true in all six dimers. Figures 3 and 4 and Table III show summaries of the equilibrium distances (*r*_{eq}) and harmonic frequencies (*ω*) estimated from the PESs and forces for the six dimers. The mean absolute errors (MAEs) of *r*_{eq} and *ω* are given in Table IV. The MAEs of the LRDMC PESs with respect to the experimental values are 0.009 13(49) Å and 14.1(6.4) cm^{−1} for *r*_{eq} and *ω*, respectively. These are much better than the corresponding VMC values, which are 0.021 50(92) Å and 27(10) cm^{−1} for *r*_{eq} and *ω*, respectively. This is, indeed, expected, and it is consistent with the consensus of the QMC community that DMC calculations usually give much better results than VMC ones.

Molecule . | Method . | r_{eq,E} (Å)
. | r_{eq,F} (Å)
. | ω_{E} (cm^{−1})
. | ω_{F} (cm^{−1})
. |
---|---|---|---|---|---|

H_{2} | VMC | 0.741 355(89) | 0.741 322(36) | 4398.3(5.1) | 4410.6(1.4) |

LRDMC(RE) | 0.741 328(92) | 0.741 319(71) | 4404.3(6.8) | 4405.1(2.6) | |

LRDMC(VD) | ⋯ | 0.741 373(72) | ⋯ | 4398.9(2.8) | |

LRDMC(Hyb) | ⋯ | 0.741 32(15) | ⋯ | 4399.7(5.7) | |

Expt.^{a} | 0.741 44 | ⋯ | 4401.21 | ⋯ | |

LiH | VMC | 1.603 1(11) | 1.607 89(28) | 1384(22) | 1388.2(2.9) |

LRDMC(RE) | 1.598 3(17) | 1.599 36(58) | 1421(33) | 1401.8(6.3) | |

LRDMC(VD) | ⋯ | 1.592 41(62) | ⋯ | 1408.3(6.3) | |

LRDMC(Hyb) | ⋯ | 1.590 9(12) | ⋯ | 1418(12) | |

Expt.^{a} | 1.595 7 | ⋯ | 1405.65 | ⋯ | |

Li_{2} | VMC | 2.762 5(41) | 2.754 55(37) | 335.3(8.3) | 333.04(60) |

LRDMC(RE) | 2.713 5(21) | 2.732 09(40) | 349.1(5.6) | 341.18(61) | |

LRDMC(VD) | ⋯ | 2.697 41(29) | ⋯ | 352.18(52) | |

LRDMC(Hyb) | ⋯ | 2.710 75(81) | ⋯ | 349.6(1.3) | |

Expt.^{a} | 2.672 9 | ⋯ | 351.43 | ⋯ | |

CO | VMC | 1.120 98(39) | 1.114 521(85) | 2232(13) | 2277.8(2.0) |

LRDMC(RE) | 1.124 12(15) | 1.116 556(57) | 2206.1(4.9) | 2257.0(1.4) | |

LRDMC(VD) | ⋯ | 1.122 9(32) | ⋯ | 2193(70) | |

LRDMC(Hyb) | ⋯ | 1.118 63(14) | ⋯ | 2235.7(3.5) | |

Expt.^{a} | 1.128 323 | ⋯ | 2169.813 58 | ⋯ | |

MgO | VMC | 1.728 0(26) | 1.728 98(19) | 830(24) | 863.4(1.8) |

LRDMC(RE) | 1.753 34(84) | 1.728 91(16) | 793.8(4.9) | 848.6(1.1) | |

LRDMC(VD) | ⋯ | 1.728 0(32) | ⋯ | 896(32) | |

LRDMC(Hyb) | ⋯ | 1.728 84(39) | ⋯ | 833.5(3.0) | |

Expt.^{a} | 1.749 | ⋯ | 785.0 | ⋯ | |

SiH | VMC | 1.523 6(25) | 1.532 81(39) | 2030(48) | 2008.9(4.7) |

LRDMC(RE) | 1.523 03(80) | 1.524 32(36) | 2061(16) | 2052.8(4.3) | |

LRDMC(VD) | ⋯ | 1.505 7(35) | ⋯ | 2156(46) | |

LRDMC(Hyb) | ⋯ | 1.516 18(81) | ⋯ | 2089.2(9.2) | |

Expt.^{a} | 1.520 1 | ⋯ | 2041.80 | ⋯ |

Molecule . | Method . | r_{eq,E} (Å)
. | r_{eq,F} (Å)
. | ω_{E} (cm^{−1})
. | ω_{F} (cm^{−1})
. |
---|---|---|---|---|---|

H_{2} | VMC | 0.741 355(89) | 0.741 322(36) | 4398.3(5.1) | 4410.6(1.4) |

LRDMC(RE) | 0.741 328(92) | 0.741 319(71) | 4404.3(6.8) | 4405.1(2.6) | |

LRDMC(VD) | ⋯ | 0.741 373(72) | ⋯ | 4398.9(2.8) | |

LRDMC(Hyb) | ⋯ | 0.741 32(15) | ⋯ | 4399.7(5.7) | |

Expt.^{a} | 0.741 44 | ⋯ | 4401.21 | ⋯ | |

LiH | VMC | 1.603 1(11) | 1.607 89(28) | 1384(22) | 1388.2(2.9) |

LRDMC(RE) | 1.598 3(17) | 1.599 36(58) | 1421(33) | 1401.8(6.3) | |

LRDMC(VD) | ⋯ | 1.592 41(62) | ⋯ | 1408.3(6.3) | |

LRDMC(Hyb) | ⋯ | 1.590 9(12) | ⋯ | 1418(12) | |

Expt.^{a} | 1.595 7 | ⋯ | 1405.65 | ⋯ | |

Li_{2} | VMC | 2.762 5(41) | 2.754 55(37) | 335.3(8.3) | 333.04(60) |

LRDMC(RE) | 2.713 5(21) | 2.732 09(40) | 349.1(5.6) | 341.18(61) | |

LRDMC(VD) | ⋯ | 2.697 41(29) | ⋯ | 352.18(52) | |

LRDMC(Hyb) | ⋯ | 2.710 75(81) | ⋯ | 349.6(1.3) | |

Expt.^{a} | 2.672 9 | ⋯ | 351.43 | ⋯ | |

CO | VMC | 1.120 98(39) | 1.114 521(85) | 2232(13) | 2277.8(2.0) |

LRDMC(RE) | 1.124 12(15) | 1.116 556(57) | 2206.1(4.9) | 2257.0(1.4) | |

LRDMC(VD) | ⋯ | 1.122 9(32) | ⋯ | 2193(70) | |

LRDMC(Hyb) | ⋯ | 1.118 63(14) | ⋯ | 2235.7(3.5) | |

Expt.^{a} | 1.128 323 | ⋯ | 2169.813 58 | ⋯ | |

MgO | VMC | 1.728 0(26) | 1.728 98(19) | 830(24) | 863.4(1.8) |

LRDMC(RE) | 1.753 34(84) | 1.728 91(16) | 793.8(4.9) | 848.6(1.1) | |

LRDMC(VD) | ⋯ | 1.728 0(32) | ⋯ | 896(32) | |

LRDMC(Hyb) | ⋯ | 1.728 84(39) | ⋯ | 833.5(3.0) | |

Expt.^{a} | 1.749 | ⋯ | 785.0 | ⋯ | |

SiH | VMC | 1.523 6(25) | 1.532 81(39) | 2030(48) | 2008.9(4.7) |

LRDMC(RE) | 1.523 03(80) | 1.524 32(36) | 2061(16) | 2052.8(4.3) | |

LRDMC(VD) | ⋯ | 1.505 7(35) | ⋯ | 2156(46) | |

LRDMC(Hyb) | ⋯ | 1.516 18(81) | ⋯ | 2089.2(9.2) | |

Expt.^{a} | 1.520 1 | ⋯ | 2041.80 | ⋯ |

^{a}

These values are taken from Ref. 60.

A
. | B
. | MAE (r_{eq,A} − r_{eq,B}) (Å)
. | MAE (ω_{A} − ω_{B}) (cm^{−1})
. |
---|---|---|---|

Expt. | E_{VMC} | 0.021 50(92) | 27(10) |

Expt. | E_{DMC} | 0.009 13(49) | 14.1(6.4) |

E_{VMC} | F_{VMC} | 0.004 90(93) | 20(10) |

E_{LRDMC} | F_{VMC} | 0.015 73(51) | 41.3(6.5) |

E_{LRDMC} | F_{RE-LRDMC} | 0.008 82(51) | 23.5(6.6) |

E_{LRDMC} | F_{VD-LRDMC} | 0.011 0(11) | 38(16) |

E_{LRDMC} | F_{Hyb-LRDMC} | 0.007 83(57) | 17.6(7.0) |

A
. | B
. | MAE (r_{eq,A} − r_{eq,B}) (Å)
. | MAE (ω_{A} − ω_{B}) (cm^{−1})
. |
---|---|---|---|

Expt. | E_{VMC} | 0.021 50(92) | 27(10) |

Expt. | E_{DMC} | 0.009 13(49) | 14.1(6.4) |

E_{VMC} | F_{VMC} | 0.004 90(93) | 20(10) |

E_{LRDMC} | F_{VMC} | 0.015 73(51) | 41.3(6.5) |

E_{LRDMC} | F_{RE-LRDMC} | 0.008 82(51) | 23.5(6.6) |

E_{LRDMC} | F_{VD-LRDMC} | 0.011 0(11) | 38(16) |

E_{LRDMC} | F_{Hyb-LRDMC} | 0.007 83(57) | 17.6(7.0) |

The self-consistency error is one of the main concerns in QMC force calculations when wave functions are not variationally optimized (i.e., with the JDFT *Ansatz*). The MAEs between a (VMC or LRDMC) PES and the corresponding forces are a measure of the self-consistency error. The MAEs between VMC energies and VMC forces are very small, 0.004 90(93) Å and 20(10) cm^{−1} for *r*_{eq} and *ω*, respectively, while the MAEs between DMC energies and RE-DMC forces are larger than the VMC ones, 0.008 82(51) Å and 23.5(6.6) cm^{−1} for *r*_{eq} and *ω*, respectively, implying that the self-consistency error is more severe in the DMC force calculations. The MAEs of the DMC energies and VD-DMC forces are as large as the RE-DMC ones, 0.011 0(11) Å and 38(16) cm^{−1} for *r*_{eq} and *ω*, respectively. Indeed, the VD forces do not show systematic improvement when compared with the Reynolds forces. The hybrid forces [MAEs: 0.007 83(57) Å and 17.6(7.0) cm^{−1} for *r*_{eq} and *ω*] slightly improve the corresponding RE forces [MAEs: 0.008 82(51) Å and 23.5(6.6) cm^{−1} for *r*_{eq} and *ω*], but not as drastically as in the original study,^{48} which is consistent with the recent DMC-force benchmark test.^{30}

We see that the MAEs of the VD-DMC forces are slightly worse than RE ones overall, but the quality of the VD force depends on each dimer. For instance, as shown in Fig. 3, the VD approximation corrects the forces in the right directions for the CO and MgO dimers, but in the cases of LiH, Li_{2}, and SiH, the corrections overshoot. In an earlier VD paper,^{23} it was mentioned that the quality of the trial function crucially affects the bias introduced within the VD approximation. This is why the correction depends on the dimers in this study. We also see that the VD forces show much larger error bars than the RE forces, especially in heavier atoms. For instance, in the SiH dimer at *d* = 1.200 Å, the ratio between the variance of force and that of energy ($\sigma F2/\sigma E2$) obtained from the same run is ∼510 with the VD approximation, while that with the RE approximation is just ∼2. The earlier VD paper also mentioned that the additional Pulay terms {$\u2211n\u2212kn\u22121\u2207\alpha \u22c5[ELxi+1\u2212ELxi]$} increased the variance,^{23} but the ratio between the variance of the force and that of the energy calculated in the same run was just ∼4^{2} = 16 at most. This is certainly because pseudopotentials were used in the VD work,^{23} whereas our calculations are all-electron. Our result implies that the alleviation of the nodal surface divergence by the PW regularization is not sufficient to significantly reduce the error bars due to the presence of the additional Pulay term {$\u2211n\u2212kn\u22121\u2207\alpha \u22c5[ELxi+1\u2212ELxi]$}, especially in all-electron calculations. A new regularization technique that can remove the divergence could alleviate the large-variance problem and make all-electron VD force calculations more practical.

Another concern is whether or not DMC forces are better than VMC forces (i.e., closer to the exact derivative obtained from a DMC PES). If DMC forces do not improve on the VMC ones, it is not worth computing the DMC force with its associated higher computational cost, for instance, to perform structural optimizations or molecular dynamics. The present study shows that the MAEs of VMC and RE-DMC forces with respect to LRDMC energies for *r*_{eq} are 0.015 73(51) Å and 0.008 82(51) Å, respectively, and those for *ω* are 41.3(6.5) cm^{−1} and 23.5(6.6) cm^{−1}, respectively. This implies that RE-DMC forces improve VMC ones, i.e., structural optimizations or molecular dynamics based on RE-DMC forces are better than those based on VMC forces at the JDFT level.

### B. Impact of WF optimization on self-consistency errors: Li_{2} dimer

The self-consistency error discussed so far can be removed by optimizing the variational parameters of the wave functions at the VMC level, and this also usually decreases the self-consistency error at the LRDMC level. For example, if a given trial wave function is exact, the RE and VD expressions are consistent with the exact force expression.^{23} Table V shows the *r*_{eq} and *ω* values of the Li_{2} dimer obtained from the PESs and atomic forces with the JDFT, JSD, and JAGPs *Ansätze*. Comparisons with experimental values are shown in Fig. 5, and the corresponding PESs are shown in the supplementary material. This result shows that, as expected, the optimized JSD and JAGPs wave functions removed the self-consistency error at the VMC level, and they also alleviated the error at the LRDMC level. Both the VMC and the LRDMC PESs show big discrepancies from the experimental *r*_{eq} and *ω* values with the JDFT *Ansatz* for the Li_{2} dimer (Fig. 5), while the optimized JAGPs *Ansatz* almost completely removes these discrepancies. The LRDMC PES with the JAGPs *Ansatz* gives 2.673 78(58) Å, 352.4(1.3) cm^{−1}, and 1.054 25(79) eV for *r*_{eq}, *ω*, and the binding energy, respectively. These are in very good agreement with the experimental values of 2.6729 Å, 351.43 cm^{−1}, and 1.055 944(67) eV and also with the CCSD(T) values, which are 2.672 Å, 355.10 cm^{−1}, and 1.057 eV.

Method . | Ansatz
. | r_{eq,E} (Å)
. | r_{eq,F} (Å)
. | ω_{E} (cm^{−1})
. | ω_{F} (cm^{−1})
. | Binding energy (eV)^{a}
. |
---|---|---|---|---|---|---|

VMC | JDFT | 2.762 5(41) | 2.754 55(37) | 335.3(8.3) | 333.04(60) | 0.776 6(14) |

JSD | 2.754 8(14) | 2.751 78(14) | 337.0(2.8) | 339.18(23) | 0.797 80(90) | |

JAGPs | 2.679 41(82) | 2.677 49(15) | 351.3(1.8) | 354.28(16) | 1.037 63(66) | |

LRDMC-RE | JDFT | 2.713 5(21) | 2.732 09(40) | 349.1(5.6) | 341.18(61) | 0.975 9(13) |

JSD | 2.719 9(15) | 2.726 52(21) | 345.3(3.8) | 346.16(36) | 0.972 58(96) | |

JAGPs | 2.673 78(58) | 2.674 35(15) | 352.4(1.3) | 354.05(21) | 1.054 25(79) | |

LRDMC-VD | JDFT | ⋯ | 2.697 41(29) | ⋯ | 352.18(52) | ⋯ |

JSD | ⋯ | 2.693 94(24) | ⋯ | 356.68(38) | ⋯ | |

JAGPs | ⋯ | 2.671 82(19) | ⋯ | 353.57(27) | ⋯ | |

CCSD(T)^{b} | ⋯ | 2.672 | ⋯ | 355.10 | ⋯ | 1.057 |

Expt. | ⋯ | 2.672 9^{c} | ⋯ | 351.43^{c} | ⋯ | 1.055 944(67)^{d} |

Method . | Ansatz
. | r_{eq,E} (Å)
. | r_{eq,F} (Å)
. | ω_{E} (cm^{−1})
. | ω_{F} (cm^{−1})
. | Binding energy (eV)^{a}
. |
---|---|---|---|---|---|---|

VMC | JDFT | 2.762 5(41) | 2.754 55(37) | 335.3(8.3) | 333.04(60) | 0.776 6(14) |

JSD | 2.754 8(14) | 2.751 78(14) | 337.0(2.8) | 339.18(23) | 0.797 80(90) | |

JAGPs | 2.679 41(82) | 2.677 49(15) | 351.3(1.8) | 354.28(16) | 1.037 63(66) | |

LRDMC-RE | JDFT | 2.713 5(21) | 2.732 09(40) | 349.1(5.6) | 341.18(61) | 0.975 9(13) |

JSD | 2.719 9(15) | 2.726 52(21) | 345.3(3.8) | 346.16(36) | 0.972 58(96) | |

JAGPs | 2.673 78(58) | 2.674 35(15) | 352.4(1.3) | 354.05(21) | 1.054 25(79) | |

LRDMC-VD | JDFT | ⋯ | 2.697 41(29) | ⋯ | 352.18(52) | ⋯ |

JSD | ⋯ | 2.693 94(24) | ⋯ | 356.68(38) | ⋯ | |

JAGPs | ⋯ | 2.671 82(19) | ⋯ | 353.57(27) | ⋯ | |

CCSD(T)^{b} | ⋯ | 2.672 | ⋯ | 355.10 | ⋯ | 1.057 |

Expt. | ⋯ | 2.672 9^{c} | ⋯ | 351.43^{c} | ⋯ | 1.055 944(67)^{d} |

### C. Computational scaling of QMC atomic force calculations

As mentioned in the Introduction, another important concern in QMC calculations is the scaling of the computational cost with the atomic number (*Z*). A recent VMC and DMC benchmark test carried out by Tiihonen *et al.*^{30} suggests that the scaling of the computational cost of forces with respect to the atomic number is much worse than that of energy even if the SWCT is applied. They claimed that the accessible system size scales approximately as *Z*^{−2}. However, within our implementation, i.e., VMC with the AS regularization and LRDMC with the RE approximation and the PW regularization combined with the SWCT, we did not observe any degradation in the scaling of forces at either the VMC or LRDMC levels for the above all-electron PES calculations.

To study the scalings more systematically, we calculated the energies and forces of homonuclear dimers—H_{2}, Li_{2}, N_{2}, F_{2}, P_{2}, S_{2}, Cl_{2}, and Br_{2}—and computed energy variances ($\sigma E2$), force variances ($\sigma F2$), and the ratios between the force and energy variances ($\sigma F2/\sigma E2$), wherein *σ*_{E} is the error bar of the total energy of the system and *σ*_{F} is that of the atomic force acting on an atom. The energy and force error bars were evaluated with the same Monte Carlo sampling. The LRDMC forces refer to those with the RE approximation. The variances of energy and force were obtained by applying the reblocking technique to remove the auto-correlations. Our findings are that the computational costs to achieve a constant statistical error in energy evaluation scale with *Z*^{∼4.5} and *Z*^{∼6.1} for VMC and LRDMC calculations, respectively. They are consistent with earlier works. For example, Hammond *et al.*^{28} and Ma *et al.*^{29} showed that the computational cost of all-electron DMC energy calculations scales with *Z*^{∼6.5} and *Z*^{∼6.0}, respectively. On the other hand, our results show that the corresponding costs in the force evaluations scale with *Z*^{∼4.4} and *Z*^{∼5.9} for VMC and LRDMC with SWCT, respectively, while they scale with *Z*^{∼7.2} and *Z*^{∼8.4} for VMC and LRDMC without SWCT, respectively. Figure 6 shows the ratios between the force and energy variances evaluated for each dimer. A summary of the fittings is shown in Table VI. We found that the ratio scales as *Z*^{2.65} and *Z*^{2.42} without the SWCT in the VMC and LRDMC calculations, respectively, reproducing the previous report.^{30} On the other hand, the ratio is *independent* of *Z* with the SWCT in both the VMC and LRDMC calculations, indicating that the computational cost of QMC forces with respect to *Z* is no worse than that of energy. Indeed, the accessible system size is not affected by QMC force calculations when the SWCT variance-reduction technique is applied. We note that the application of the SWCT is also drastic in pseudopotential calculations. To confirm this, we computed the ratio of Cl_{2} (*Z* = 17, *Z*_{eff} = 15) with the He-core ccECP pseudopotential^{64} at the VMC level. We obtained $(\sigma F/\sigma E)2\u223c1$ with the SWCT, while $(\sigma F/\sigma E)2\u223c317$ without the SWCT, indicating that the SWCT is also effective for pseudopotential calculations.

Method . | SWCT . | α
. | β
. |
---|---|---|---|

VMC | No | 1.04 | 2.65 |

Yes | 1.88 | −0.15 | |

LRDMC | No | 1.15 | 2.32 |

Yes | 1.97 | −0.24 |

Method . | SWCT . | α
. | β
. |
---|---|---|---|

VMC | No | 1.04 | 2.65 |

Yes | 1.88 | −0.15 | |

LRDMC | No | 1.15 | 2.32 |

Yes | 1.97 | −0.24 |

The following argument clarifies, in our opinion, why the SWCT is so effective for reducing the variance for large atomic numbers. First of all, the SWCT provides zero forces with zero variances for the calculations of forces of isolated atoms by definition, as already derived in this paper and also pointed out in earlier works.^{46,49} Although the result is trivial in this case, evaluating the zero force without the SWCT becomes problematic in QMC because the derivative of the local energy (required for the HF term) and the derivative of the logarithm of the wave function (required for the Pulay term) scale with ≃*Z* in the vicinity of nuclei because these functions have length scales of order 1/*Z*, implying large fluctuations of the mentioned derivatives in a QMC sampling. Instead, the application of the SWCT is able to suppress the fluctuations because electron coordinates follow the displacement of the atomic position in the vicinity of the nucleus. This argument clearly holds in complex systems containing several atoms because the damping function (*ω*_{α}) makes electron coordinations follow an atomic displacement in the vicinity of the nucleus, and the mentioned differentiations are, therefore, smooth, as in the isolated case. In this way, one can suppress the large fluctuations originating from small electron–ion distances.

## V. CONCLUDING REMARKS

In this study, we benchmarked the accuracy of all-electron VMC and LRDMC forces for various mono- and heteronuclear dimers (1 ≤ *Z* ≤ 35). Although all our calculations were carried out by considering all the electrons of the corresponding atomic species, we do not expect—nor do we have evidence that—the main conclusions of our work will depend on the presence of pseudopotentials, which are often used in QMC to remove the chemically inert core electrons. The VMC and LRDMC forces were calculated in combination with the SWCT and appropriate regularization techniques to remove the infinite variance problem. The equilibrium distances and the harmonic frequencies of the six dimers (H_{2}, LiH, Li_{2}, CO, MgO, and SiH) obtained from the DMC PESs are very close to the experimental values, but the RE and VD forces are less close as a result of residual self-consistency errors. Self-consistency errors in forces are also seen at the VMC level, but they are much smaller than those in the DMC forces. As shown in the Li_{2} dimer case, the self-consistency error with the JDFT *Ansatz* can be eliminated by carefully optimizing the determinant part, while until now, there has been no practical way to remove it at the DMC level. Since VMC calculations are sometimes not accurate enough to describe peculiar molecular interactions such as water clusters, developing a practical DMC force calculation without suffering from self-consistency errors is an important challenge for QMC applications. We also found that the ratio of the computational costs between QMC energies and forces scales as *Z*^{∼2.5} without the SWCT, while the application of the SWCT makes this ratio *independent* of *Z*. Indeed, the accessible system size is not affected by QMC force calculations when the SWCT variance-reduction technique is employed. This is quite a promising conclusion from the viewpoint of the QMC computation of forces.

## SUPPLEMENTARY MATERIAL

See the supplementary material for PESs of the six dimers (H_{2}, LiH, Li_{2}, CO, MgO, and SiH).

## ACKNOWLEDGMENTS

The computations in this work were mainly performed using the supercomputer Fugaku, which was provided by RIKEN through the HPCI System Research Project (Project No. hp210038). This work was supported by the European Centre of Excellence in Exascale Computing TREX—Targeting Real Chemical Accuracy at the Exascale. This project received funding from the European Union’s Horizon 2020 research and innovation program (Grant No. 952165). K.N. and A.R. acknowledge computational resources from the facilities of the Research Center for Advanced Computing Infrastructure at Japan Advanced Institute of Science and Technology (JAIST). K.N. acknowledges support from the JSPS Overseas Research Fellowships, from the Grant-in-Aid for Early Career Scientists (Grant No. JP21K17752), and from the Grant-in-Aid for Scientific Research (Grant No. JP21K03400). A.R. acknowledges the financial support from the MEXT (Monbukagakusho) scholarship. S.S. acknowledges financial support from PRIN 2017BZPKSZ. The authors would like to thank C. Filippi (University of Twente) and S. Moroni (SISSA) for valuable comments through the TREX collaboration. The authors are also grateful to A. Zen (University of Naples Federico II) for fruitful discussion and for providing them with Orca input files for CCSD(T) calculations.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material and the cited references. TurboRVB is available from S.S.’s website (https://people.sissa.it/∼sorella/) upon reasonable request.

## REFERENCES

The Hellmann–Feynman (HF) term is just conventionally used when the SWCT is applied. This is because the mean value of the term in the presence of the SWCT does not correspond to the expectation value of $\u2202RH\u0302$.

In TurboRVB, the memory requirement for the adjoint algorithm differentiation is kept to a reasonable small value (less than the one required by the direct algorithm) by recomputing intermediate variables of the direct algorithm on the fly, with an almost negligible computational overhead.