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 (req) 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.

The quantum Monte Carlo (QMC) technique1 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 minimization15,16 and regularization12 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—H2, LiH, Li2, CO, MgO, and SiH—at the VMC and LRDMC levels and investigated their corresponding consistency. The LRDMC forces were computed using the Reynolds (RE) approximation17 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 N7, while that of DMC methods grows as N4 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 algorithm25 and the real-space double-grid method26 have been devised to alleviate this problem. The scaling with Z is, however, still higher than that with N; for example, it is Z5.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 Zeff2, where Zeff 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: H2, Li2, N2, F2, P2, S2, Cl2, and Br2. 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.

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

(1)

where x=r1σ1,r2σ2,,rNσ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α,

(2a)
(2b)
(2c)

wherein Eqs. (2a)(2c) are called the Hellmann–Feynman (HF),47 Pulay, and variational terms, respectively. In the variational term, the set of {c1,,cNc} are the variational parameters included in the wave function.

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

(3)

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., Eci=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α=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 (EL) and the logarithmic derivative of the wave function (logΨRα) 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.

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 (logΦRα). Therefore, to overcome this difficulty, Reynolds et al.17 proposed the approximation logΦRαlogΨRα, leading to

(4)

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

(5)

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

(6)

where τ is the time step defining the discretized times tn = 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 Li2 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 Wagner52 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=Ψ(x)/Ψ(x).

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

(7)
(8)

where κr is a function that decays for large r and ω → 0 as a consequence of this property. The original46 simple choice κr=1/r4 is adopted in TurboRVB. The atomic force is evaluated with the SWCT,49 

(9)

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

(10)

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 xj. All the terms above can be very efficiently evaluated by the adjoint algorithmic differentiation49 in TurboRVB.53 Some of the adjoint routines in TurboRVB were generated from the corresponding forward ones using TAPENADE.54 

The all-electron VMC and LRDMC calculations for the selected dimers (H2, LiH, Li2, N2, CO, F2, MgO, SiH, P2, S2, Cl2, and Br2) were performed using TurboRVB,55 and this was combined with the Python package Turbo-Genius55 to avoid any human error and to achieve efficient calculations. TurboRVB employs the Jastrow antisymmetrized geminal power (JAGP)56Ansatz. 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

(11)

where  is the antisymmetrization operator and Φr,r is called the pairing function. The spatial part of the geminal function is expanded over the Gaussian-type atomic orbitals (GTOs),

(12)

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 = Nel/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 Li2 dimer. The JDFT, JSD, and JAGPs Ansätze were employed for the Li2 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 = J1J2J3/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:

(13)

where

(14)

ri are the electron positions, RI are the atomic positions with the corresponding atomic number ZI, l runs over atomic orbitals χI,lJ centered on atom I, and ur is a short-range function containing a variational parameter b,

(15)

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

(16)

where vr is given as

(17)

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

(18)

and

(19)

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 library59 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 · Z2 (where Z is the atomic number) are disregarded, and they are implicitly compensated by the homogeneous one-body Jastrow part [J̃1r] to fulfill the electron–ion cusp condition explicitly, even at the DFT level.58 Indeed, a single-particle orbital is modified as

(20)

where J̃1r 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 library59 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 method15,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 scheme4 with lattice spaces a = 0.30, 0.10, 0.10, 0.05, 0.05, and 0.03 bohr for H2, LiH, Li2, 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 H2, LiH, Li2, CO, MgO, and SiH, respectively, for the PW regularization parameter (ϵ). These are compatible with the values proposed in the original paper52 (∼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 scheme4 with lattice spaces a = 1/(2.0 · Z) bohr for H2, Li2, N2, F2, P2, S2, Cl2, and Br2. The regularization parameter ϵ was set to 0.001 bohr for all the dimers.

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:

(21)

and

(22)

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

(23)

and

(24)

Then, as a consequence of αωαri=1, the sum of the HF forces in Eq. (9) over all the atoms in a system is

(25)

and the sum of the Pulay forces is

(26)
FIG. 1.

Potential energy curves of the CO dimer obtained by (a) VMC calculations with the AS regularization and (b) LRDMC calculations with the PW regularization. The JDFT Ansatz was employed in these calculations. The derivatives and atomic forces acting on the left (C) and right (O) atoms are also shown. The vertical broken lines show the equilibrium bond lengths obtained from previous experiments60 and from these calculations.

FIG. 1.

Potential energy curves of the CO dimer obtained by (a) VMC calculations with the AS regularization and (b) LRDMC calculations with the PW regularization. The JDFT Ansatz was employed in these calculations. The derivatives and atomic forces acting on the left (C) and right (O) atoms are also shown. The vertical broken lines show the equilibrium bond lengths obtained from previous experiments60 and from these calculations.

Close modal
FIG. 2.

Comparison of the derivatives of the PESs (solid lines) and the atomic forces (broken lines) of the CO dimer obtained from the VMC and LRDMC calculations, corresponding to Fig. 1. RE, VD, and Hyb refer to the Reynolds17 approximation, the variational-drift23 approximation, and hybrid forces.48 The vertical black broken line indicates the equilibrium bond length (req) obtained from the experiment.60 The CO dimer was depicted by VESTA.61 

FIG. 2.

Comparison of the derivatives of the PESs (solid lines) and the atomic forces (broken lines) of the CO dimer obtained from the VMC and LRDMC calculations, corresponding to Fig. 1. RE, VD, and Hyb refer to the Reynolds17 approximation, the variational-drift23 approximation, and hybrid forces.48 The vertical black broken line indicates the equilibrium bond length (req) obtained from the experiment.60 The CO dimer was depicted by VESTA.61 

Close modal

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.

TABLE I.

Left and right VMC forces of the Li2 dimer (d = 2.400 Å) obtained with and without the SWCT.

Li2Left 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) 
Li2Left 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) 
TABLE II.

Left and right VMC forces of the SiH dimer (d = 1.300 Å) obtained with and without the SWCT, where the left and right forces refer to the forces acting on the Si and H atoms, respectively.

SiHLeft 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) 
SiHLeft 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 req 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 (req) and harmonic frequencies (ω) estimated from the PESs and forces for the six dimers. The mean absolute errors (MAEs) of req 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 req and ω, respectively. These are much better than the corresponding VMC values, which are 0.021 50(92) Å and 27(10) cm−1 for req 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.

FIG. 3.

Deviations of the obtained bond lengths req from the experimentally obtained values.60 The values obtained from the VMC PESs and atomic forces are shown in the left-hand panel, while those obtained from the LRDMC PESs and atomic forces with the RE, VD, and hybrid forces are shown in the right-hand panel. The gray bounds represent acceptable discrepancy, ±0.01 Å. The corresponding numerical values are shown in Table III.

FIG. 3.

Deviations of the obtained bond lengths req from the experimentally obtained values.60 The values obtained from the VMC PESs and atomic forces are shown in the left-hand panel, while those obtained from the LRDMC PESs and atomic forces with the RE, VD, and hybrid forces are shown in the right-hand panel. The gray bounds represent acceptable discrepancy, ±0.01 Å. The corresponding numerical values are shown in Table III.

Close modal
FIG. 4.

Deviations of the obtained harmonic frequencies ω from the experimentally obtained values.60 The values obtained from the VMC PESs and atomic forces are shown in the left-hand panel, while those obtained from the LRDMC PESs and atomic forces with the RE, VD, and hybrid forces are shown in the right-hand panel. The gray bounds represent acceptable discrepancy, ±50 cm−1. The corresponding numerical values are shown in Table III.

FIG. 4.

Deviations of the obtained harmonic frequencies ω from the experimentally obtained values.60 The values obtained from the VMC PESs and atomic forces are shown in the left-hand panel, while those obtained from the LRDMC PESs and atomic forces with the RE, VD, and hybrid forces are shown in the right-hand panel. The gray bounds represent acceptable discrepancy, ±50 cm−1. The corresponding numerical values are shown in Table III.

Close modal
TABLE III.

Equilibrium bond distances req (Å) and harmonic frequencies ω (cm−1) calculated in the present work alongside experimental data.60 All calculations were carried out using the JDFT Ansatz. The subscripts E and F indicate values obtained from the PES and the atomic forces, respectively.

MoleculeMethodreq,E (Å)req,F (Å)ωE (cm−1)ωF (cm−1)
H2 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 ⋯ 
Li2 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 ⋯ 
MoleculeMethodreq,E (Å)req,F (Å)ωE (cm−1)ωF (cm−1)
H2 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 ⋯ 
Li2 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.

TABLE IV.

Mean absolute errors (MAEs) of req and ω.

ABMAE (req,Areq,B) (Å)MAE (ωAωB) (cm−1)
Expt. EVMC 0.021 50(92) 27(10) 
Expt. EDMC 0.009 13(49) 14.1(6.4) 
EVMC FVMC 0.004 90(93) 20(10) 
ELRDMC FVMC 0.015 73(51) 41.3(6.5) 
ELRDMC FRE-LRDMC 0.008 82(51) 23.5(6.6) 
ELRDMC FVD-LRDMC 0.011 0(11) 38(16) 
ELRDMC FHyb-LRDMC 0.007 83(57) 17.6(7.0) 
ABMAE (req,Areq,B) (Å)MAE (ωAωB) (cm−1)
Expt. EVMC 0.021 50(92) 27(10) 
Expt. EDMC 0.009 13(49) 14.1(6.4) 
EVMC FVMC 0.004 90(93) 20(10) 
ELRDMC FVMC 0.015 73(51) 41.3(6.5) 
ELRDMC FRE-LRDMC 0.008 82(51) 23.5(6.6) 
ELRDMC FVD-LRDMC 0.011 0(11) 38(16) 
ELRDMC FHyb-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 req 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 req 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 req 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 req and ω] slightly improve the corresponding RE forces [MAEs: 0.008 82(51) Å and 23.5(6.6) cm−1 for req 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, Li2, 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 (σF2/σ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 {nkn1α[ELxi+1ELxi]} 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 ∼42 = 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 {nkn1α[ELxi+1ELxi]}, 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 req 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.

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.23Table V shows the req and ω values of the Li2 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 req and ω values with the JDFT Ansatz for the Li2 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 req, ω, 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.

TABLE V.

Equilibrium bond distances req (Å) and harmonic frequencies ω (cm−1) of the Li2 dimer obtained with the JDFT, JSD, and JAGPs Ansätze. The subscripts E and F indicate values obtained from the PES and the atomic forces, respectively.

MethodAnsatzreq,E (Å)req,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 9c ⋯ 351.43c ⋯ 1.055 944(67)d 
MethodAnsatzreq,E (Å)req,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 9c ⋯ 351.43c ⋯ 1.055 944(67)d 
a

Binding energy at the obtained req,E without the zero-point energy.

b

These values were obtained by Orca62 with extrapolation using cc-pCVnZ basis sets (n = 3, 4).

c

These values are taken from Ref. 60.

d

This value is taken from Ref. 63.

FIG. 5.

Equilibrium bond lengths (left) and harmonic frequencies (right) of the Li2 dimer obtained with various Ansätze and methods. The deviations from the experimental values60 are plotted. The DMC forces here refer to those with the Reynolds approximation. The gray bounds represent acceptable discrepancies, ± 0.01 Å and ±50 cm−1.

FIG. 5.

Equilibrium bond lengths (left) and harmonic frequencies (right) of the Li2 dimer obtained with various Ansätze and methods. The deviations from the experimental values60 are plotted. The DMC forces here refer to those with the Reynolds approximation. The gray bounds represent acceptable discrepancies, ± 0.01 Å and ±50 cm−1.

Close modal

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—H2, Li2, N2, F2, P2, S2, Cl2, and Br2—and computed energy variances (σE2), force variances (σF2), and the ratios between the force and energy variances (σF2/σ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 Z2.65 and Z2.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 Cl2 (Z = 17, Zeff = 15) with the He-core ccECP pseudopotential64 at the VMC level. We obtained (σF/σE)21 with the SWCT, while (σF/σE)2317 without the SWCT, indicating that the SWCT is also effective for pseudopotential calculations.

FIG. 6.

Ratio of variances between force and energy (σF2/σE2) evaluated with the same run. The VMC and LRDMC forces were evaluated with/without the SWCT, where the LRDMC force refers to the force computed with the RE approximation. The VMC and LRDMC results were fitted with power-law scaling relations.

FIG. 6.

Ratio of variances between force and energy (σF2/σE2) evaluated with the same run. The VMC and LRDMC forces were evaluated with/without the SWCT, where the LRDMC force refers to the force computed with the RE approximation. The VMC and LRDMC results were fitted with power-law scaling relations.

Close modal
TABLE VI.

Scalings of σF2/σE2 with respect to the atomic number obtained by fitting the VMC and LRDMC forces of the homonuclear dimers (H2, Li2, N2, F2, P2, S2, Cl2, and Br2) with αZβ.

MethodSWCTαβ
VMC No 1.04 2.65 
Yes 1.88 −0.15 
LRDMC No 1.15 2.32 
Yes 1.97 −0.24 
MethodSWCTαβ
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.

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 (H2, LiH, Li2, 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 Li2 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.

See the supplementary material for PESs of the six dimers (H2, LiH, Li2, CO, MgO, and SiH).

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.

The authors have no conflicts to disclose.

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.

1.
W. M. C.
Foulkes
,
L.
Mitas
,
R. J.
Needs
, and
G.
Rajagopal
,
Rev. Mod. Phys.
73
,
33
(
2001
).
2.
C. J.
Umrigar
,
Phys. Rev. Lett.
71
,
408
(
1993
).
3.
M. L.
Stedman
,
W. M. C.
Foulkes
, and
M.
Nekovee
,
J. Chem. Phys.
109
,
2630
(
1998
).
4.
M.
Casula
,
C.
Filippi
, and
S.
Sorella
,
Phys. Rev. Lett.
95
,
100201
(
2005
).
5.
F.
Becca
and
S.
Sorella
,
Quantum Monte Carlo Approaches for Correlated Systems
(
Cambridge University Press
,
2017
).
6.
A.
Zen
,
J. G.
Brandenburg
,
J.
Klimeš
,
A.
Tkatchenko
,
D.
Alfè
, and
A.
Michaelides
,
Proc. Natl. Acad. Sci. U. S. A.
115
,
1724
(
2018
).
7.
K.
Hongo
,
M. A.
Watson
,
R. S.
Sánchez-Carrera
,
T.
Iitaka
, and
A.
Aspuru-Guzik
,
J. Phys. Chem. Lett.
1
,
1789
(
2010
).
8.
E.
Mostaani
,
N. D.
Drummond
, and
V. I.
Fal’ko
,
Phys. Rev. Lett.
115
,
115501
(
2015
).
9.
S.
Sorella
,
K.
Seki
,
O. O.
Brovko
,
T.
Shirakawa
,
S.
Miyakoshi
,
S.
Yunoki
, and
E.
Tosatti
,
Phys. Rev. Lett.
121
,
066402
(
2018
).
10.
T.
Frank
,
R.
Derian
,
K.
Tokár
,
L.
Mitas
,
J.
Fabian
, and
I.
Štich
,
Phys. Rev. X
9
,
011018
(
2019
).
11.
M.
Casula
and
S.
Sorella
,
Phys. Rev. B
88
,
155125
(
2013
).
12.
C.
Attaccalite
and
S.
Sorella
,
Phys. Rev. Lett.
100
,
114501
(
2008
).
13.
N. D.
Drummond
,
B.
Monserrat
,
J. H.
Lloyd-Williams
,
P. L.
Ríos
,
C. J.
Pickard
, and
R. J.
Needs
,
Nat. Commun.
6
,
7794
(
2015
).
14.
G.
Mazzola
,
R.
Helled
, and
S.
Sorella
,
Phys. Rev. Lett.
120
,
025701
(
2018
).
15.
C. J.
Umrigar
,
J.
Toulouse
,
C.
Filippi
,
S.
Sorella
, and
R. G.
Hennig
,
Phys. Rev. Lett.
98
,
110201
(
2007
).
16.
S.
Sorella
,
M.
Casula
, and
D.
Rocca
,
J. Chem. Phys.
127
,
014105
(
2007
).
17.
P. J.
Reynolds
,
R. N.
Barnett
,
B. L.
Hammond
,
R. M.
Grimes
, and
W. A.
Lester
, Jr.
,
Int. J. Quantum Chem.
29
,
589
(
1986
).
18.
A.
Badinski
and
R. J.
Needs
,
Phys. Rev. B
78
,
035134
(
2008
).
19.
R.
Assaraf
and
M.
Caffarel
,
J. Chem. Phys.
113
,
4028
(
2000
).
20.
R.
Assaraf
,
M.
Caffarel
, and
A. C.
Kollias
,
Phys. Rev. Lett.
106
,
150601
(
2011
).
21.
S.
Chiesa
,
D. M.
Ceperley
, and
S.
Zhang
,
Phys. Rev. Lett.
94
,
036404
(
2005
).
22.
C.
Filippi
and
C. J.
Umrigar
,
Phys. Rev. B
61
,
R16291
(
2000
).
23.
S.
Moroni
,
S.
Saccani
, and
C.
Filippi
,
J. Chem. Theory Comput.
10
,
4823
(
2014
).
24.
J.
van Rhijn
,
C.
Filippi
,
S.
De Palo
, and
S.
Moroni
,
J. Chem. Theory Comput.
18
,
118
123
(
2022
).
25.
C. J.
Umrigar
,
M. P.
Nightingale
, and
K. J.
Runge
,
J. Chem. Phys.
99
,
2865
(
1993
).
26.
K.
Nakano
,
R.
Maezono
, and
S.
Sorella
,
Phys. Rev. B
101
,
155106
(
2020
).
27.
D. M.
Ceperley
,
J. Stat. Phys.
43
,
815
(
1986
).
28.
B. L.
Hammond
,
P. J.
Reynolds
, and
W. A.
Lester
, Jr.
,
J. Chem. Phys.
87
,
1130
(
1987
).
29.
A.
Ma
,
N. D.
Drummond
,
M. D.
Towler
, and
R. J.
Needs
,
Phys. Rev. E
71
,
066704
(
2005
).
30.
J.
Tiihonen
,
R. C.
Clay
 III
, and
J. T.
Krogel
,
J. Chem. Phys.
154
,
204111
(
2021
).
31.
L.
Wagner
and
L.
Mitas
,
Chem. Phys. Lett.
370
,
412
(
2003
).
32.
D.
Alfè
,
M.
Alfredsson
,
J.
Brodholt
,
M. J.
Gillan
,
M. D.
Towler
, and
R. J.
Needs
,
Phys. Rev. B
72
,
014114
(
2005
).
33.
I.
Kylänpää
,
J.
Balachandran
,
P.
Ganesh
,
O.
Heinonen
,
P. R. C.
Kent
, and
J. T.
Krogel
,
Phys. Rev. Mater.
1
,
065408
(
2017
).
34.
L. K.
Wagner
and
L.
Mitas
,
J. Chem. Phys.
126
,
034105
(
2007
).
35.
J.
Koseki
,
R.
Maezono
,
M.
Tachikawa
,
M. D.
Towler
, and
R. J.
Needs
,
J. Chem. Phys.
129
,
085103
(
2008
).
36.
H.
Zheng
and
L. K.
Wagner
,
Phys. Rev. Lett.
114
,
176401
(
2015
).
37.
Y.
Luo
,
A.
Benali
,
L.
Shulenburger
,
J. T.
Krogel
,
O.
Heinonen
, and
P. R. C.
Kent
,
New J. Phys.
18
,
113049
(
2016
).
38.
J. R.
Trail
and
R. J.
Needs
,
J. Chem. Phys.
146
,
204107
(
2017
).
39.
H.
Shin
,
Y.
Luo
,
P.
Ganesh
,
J.
Balachandran
,
J. T.
Krogel
,
P. R. C.
Kent
,
A.
Benali
, and
O.
Heinonen
,
Phys. Rev. Mater.
1
,
073603
(
2017
).
40.
T.
Ichibha
,
Z.
Hou
,
K.
Hongo
, and
R.
Maezono
,
Sci. Rep.
7
,
2011
(
2017
).
41.
H.
Shin
,
A.
Benali
,
Y.
Luo
,
E.
Crabb
,
A.
Lopez-Bezanilla
,
L. E.
Ratcliff
,
A. M.
Jokisaari
, and
O.
Heinonen
,
Phys. Rev. Mater.
2
,
075001
(
2018
).
42.
S.
Song
,
M.-C.
Kim
,
E.
Sim
,
A.
Benali
,
O.
Heinonen
, and
K.
Burke
,
J. Chem. Theory Comput.
14
,
2304
(
2018
).
43.
T.
Ichibha
,
A.
Benali
,
K.
Hongo
, and
R.
Maezono
,
Phys. Rev. Mater.
3
,
125801
(
2019
).
44.
J.
Shee
,
B.
Rudshteyn
,
E. J.
Arthur
,
S.
Zhang
,
D. R.
Reichman
, and
R. A.
Friesner
,
J. Chem. Theory Comput.
15
,
2346
(
2019
).
45.
T.
Ichibha
,
A. L.
Dzubak
,
J. T.
Krogel
,
V. R.
Cooper
, and
F. A.
Reboredo
,
Phys. Rev. Mater.
5
,
064006
(
2021
).
46.
C. J.
Umrigar
,
Int. J. Quantum Chem.
36
,
217
(
1989
).
47.

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 RĤ.

48.
R.
Assaraf
and
M.
Caffarel
,
J. Chem. Phys.
119
,
10536
(
2003
).
49.
S.
Sorella
and
L.
Capriotti
,
J. Chem. Phys.
133
,
234111
(
2010
).
50.
C.
Filippi
,
R.
Assaraf
, and
S.
Moroni
,
J. Chem. Phys.
144
,
194105
(
2016
).
51.
O.
Valsson
and
C.
Filippi
,
J. Chem. Theory Comput.
6
,
1275
(
2010
).
52.
S.
Pathak
and
L. K.
Wagner
,
AIP Adv.
10
,
085213
(
2020
).
53.

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.

54.
L.
Hascoet
and
V.
Pascual
,
ACM Trans. Math. Software
39
,
1
(
2013
).
55.
K.
Nakano
,
C.
Attaccalite
,
M.
Barborini
,
L.
Capriotti
,
M.
Casula
,
E.
Coccia
,
M.
Dagrada
,
C.
Genovese
,
Y.
Luo
,
G.
Mazzola
,
A.
Zen
, and
S.
Sorella
,
J. Chem. Phys.
152
,
204121
(
2020
).
56.
M.
Casula
and
S.
Sorella
,
J. Chem. Phys.
119
,
6500
(
2003
).
57.
M.
Marchi
,
S.
Azadi
,
M.
Casula
, and
S.
Sorella
,
J. Chem. Phys.
131
,
154116
(
2009
).
58.
K.
Nakano
,
R.
Maezono
, and
S.
Sorella
,
J. Chem. Theory Comput.
15
,
4044
(
2019
).
59.
B. P.
Pritchard
,
D.
Altarawy
,
B.
Didier
,
T. D.
Gibson
, and
T. L.
Windus
,
J. Chem. Inf. Model.
59
,
4814
(
2019
).
60.
K.-P.
Huber
,
Molecular Spectra and Molecular Structure: IV. Constants of Diatomic Molecules
(
Springer Science & Business Media
,
2013
).
61.
K.
Momma
and
F.
Izumi
,
J. Appl. Crystallogr.
44
,
1272
(
2011
).
62.
F.
Neese
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1327
(
2018
).
63.
B.
Barakat
,
R.
Bacis
,
F.
Carrot
,
S.
Churassy
,
P.
Crozet
,
F.
Martin
, and
J.
Verges
,
Chem. Phys.
102
,
215
(
1986
).
64.
M. C.
Bennett
,
G.
Wang
,
A.
Annaberdiyev
,
C. A.
Melton
,
L.
Shulenburger
, and
L.
Mitas
,
J. Chem. Phys.
149
,
104108
(
2018
).

Supplementary Material