# Fileset

[PhysRevB.109.205151.pdf](https://mdr.nims.go.jp/filesets/b8177852-efa3-4871-97cf-40d46021d958/download)

## Creator

[Kousuke Nakano](https://orcid.org/0000-0001-7756-4355), Michele Casula, Giacomo Tenti

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Efficient calculation of unbiased atomic forces in <i>ab initio</i> variational Monte Carlo](https://mdr.nims.go.jp/datasets/9395e8ac-076a-4350-866b-050e1ebe3efc)

## Fulltext

Efficient calculation of unbiased atomic forces in ab initio variational Monte CarloPHYSICAL REVIEW B 109, 205151 (2024)Efficient calculation of unbiased atomic forces in ab initio variational Monte CarloKousuke Nakano ,1,* Michele Casula,2 and Giacomo Tenti 31Center for Basic Research on Materials, National Institute for Materials Science (NIMS), Tsukuba, Ibaraki 305-0047, Japan2Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC), Sorbonne Université, CNRS UMR 7590,IRD UMR 206, MNHN, 4 Place Jussieu, 75252 Paris, France3International School for Advanced Studies (SISSA), Via Bonomea 265, 34136 Trieste, Italy(Received 29 December 2023; revised 1 April 2024; accepted 7 May 2024; published 28 May 2024)Ab initio quantum Monte Carlo (QMC) is a state-of-the-art numerical approach for evaluating accurateexpectation values of many-body wave functions. However, one of the major drawbacks that still hinderswidespread QMC applications is the lack of an affordable scheme to compute unbiased atomic forces. In thisstudy, we propose an efficient method to obtain unbiased atomic forces and pressures in the variational MonteCarlo (VMC) framework with the Jastrow-correlated Slater determinant ansatz or the Jastrow antisymmetrizedgeminal power ansatz, exploiting the gauge-invariant and locality properties of their geminal representation. Wedemonstrate the effectiveness of our method for H2 and Cl2 molecules and for the cubic boron nitride crystal.Our framework has a better algorithmic scaling with the system size than the traditional finite-difference methodand, in practical applications, is as efficient as single-point VMC calculations. Thus, it paves the way to studydynamical properties of materials, such as phonons, and is beneficial for pursuing more reliable machine-learninginteratomic potentials based on unbiased VMC forces.DOI: 10.1103/PhysRevB.109.205151I. INTRODUCTIONAb initio quantum Monte Carlo (QMC) [1] is a state-of-the-art numerical approach for evaluating the expectation valuesof many-body wave functions. It usually provides extremelyaccurate energies. To date, QMC has been successfullyapplied to various materials for which other electronic struc-ture methods, such as the density-functional theory (DFT),lose predictive power. Examples are molecular crystals [2],two-dimensional materials [3–5], superconductors [6], andmaterials at extreme pressures [7–11]. Despite several suc-cessful applications done so far and the recent developmentof sophisticated QMC packages [12–16], this technique isnot as widely used as other established electronic structuremethods. If compared with DFT [17], one of the main QMCdrawbacks is the lack of an efficient and affordable scheme tocompute atomic forces consistent with the derivatives of thetotal energy with respect to atomic positions (a.k.a. unbiasedatomic forces). This problem is relevant in the constructionof machine-learning potentials (MLPs), which need largedatasets, where energy and forces are computed with themethod of choice. Recently, some QMC-driven MLPs havebeen reported [18–22], where the availability of unbiasedforces and pressures has been a major concern.*kousuke_1123@icloud.comPublished by the American Physical Society under the terms of theCreative Commons Attribution 4.0 International license. Furtherdistribution of this work must maintain attribution to the author(s)and the published article’s title, journal citation, and DOI.There are two main real-space QMC frameworks, the vari-ational Monte Carlo (VMC) and the fixed-node diffusionMonte Carlo (FN-DMC) methods [1]. In this study, we focuson VMC because the forces computation within the FN-DMCframework is much more difficult and it is still a highlydebated topic [23–30]. Let Rα be the atomic position of thenucleus α. The atomic force acting on α is defined as thenegative gradient of the energy with respect to Rα:Fα = − dEdRα= −〈∂∂RαEL〉(1a)− 2〈(EL − E )∂ log �T∂Rα〉(1b)−Np∑i=1∂E∂ pid pidRα, (1c)where �T is the variational wave function, 〈A〉 indicatesthe quantum average of the local operator A over theVMC sampling of |�T|2, EL is the so-called local en-ergy (EL ≡ Ĥ�T/�T), with E ≡ 〈EL〉, and {p1, . . . , pNp} isthe set of Np variational parameters included in the �Tansatz. Equations (1a)–(1c) are called the Hellmann–Feynman(HF), Pulay, and variational terms, respectively. One usu-ally ignores Eq. (1c) when evaluating atomic VMC forces,resulting inFVMCα = −〈∂∂RαEL〉− 2〈(EL − E )∂ log �T∂Rα〉. (2)The long-standing problem of obtaining a statistically mean-ingful FVMCα value with a finite variance and at the samecost as the VMC energy evaluation has been solved by the2469-9950/2024/109(20)/205151(7) 205151-1 Published by the American Physical Societyhttps://orcid.org/0000-0001-7756-4355https://orcid.org/0000-0002-0165-9056https://ror.org/026v1ze26https://ror.org/02en5vm52https://ror.org/004fze387https://crossmark.crossref.org/dialog/?doi=10.1103/PhysRevB.109.205151&domain=pdf&date_stamp=2024-05-28https://doi.org/10.1103/PhysRevB.109.205151https://creativecommons.org/licenses/by/4.0/NAKANO, CASULA, AND TENTI PHYSICAL REVIEW B 109, 205151 (2024)FIG. 1. Schematic picture of PESs as a function of the dimerbond length R. (a) The exact PES, not accessible in practice. (b) Thebest possible PES obtained in the VMC framework by minimizing allvariational parameters of �T. (c) The PES obtained with optimizedJastrow factor and Slater MOs yielded by DFT at each point R, whoseslope at R′ is exactly given by the VMC force F VMC supplementedby the variational term F c, as proposed in this study. (d) The PESobtained with frozen DFT orbitals computed at R′, whose slopecorresponds to F VMC without the additional term F c.zero-variance zero-bias principle [31] together with thespace-warp transformation [32] and reweighting techniques[24,31,33–36]. Hereafter, we will denote FVMCα as regularVMC force [Eq. (2)].Neglecting Eq. (1c) is justified only when the system is atits variational minimum for all parameters (i.e., ∂E/∂ pi = 0,∀i) or when the variational parameters, which are implicitlydependent on the atomic positions, accidentally or by con-struction become position independent (i.e., d pi/dRα = 0,∀i); otherwise, FVMCα can be biased. This bias is referred toas self-consistency error [36,37].In this paper, we propose a method to obtain unbiasedatomic forces and pressures that does not increase the com-putational complexity of the VMC energy calculation, bysupplementing the regular VMC force with a suitable varia-tional term, computed by exploiting the gauge-invariant andlocality properties of the antisymmetrized geminal power(AGP) ansatz [38]. For assessment, we demonstrate that thepotential energy surfaces (PESs) of the H2 and Cl2 molecules,and the equation of state (EOS) of the cubic boron nitride(cBN) are consistent with the forces and pressure obtained byour proposed method.II. ILLUSTRATING THE PROBLEMFor the sake of clarity, we present the case of the PES ofa dimer expressed as a function of the interatomic distance R,while the present discussion can be applied for any other sys-tem. Figure 1 shows a schematic picture of several PESs. LetE exact [Fig. 1(a)] be the exact PES of the dimer. E exact is theultimate goal of any electronic structure calculation, but it isunknown except for nodeless ground states. The best possiblePES within a given �T ansatz is E fullopt [Fig. 1(b)], yieldedby a VMC calculation with the fully optimized �T. This isachievable for rather small systems by optimization methodssuitable for noisy data [39,40], but becomes impractical forlarger ones. Therefore, a good compromise between accuracyand computational efficiency is EJSD [Fig. 1(c)], obtained bythe Jastrow correlated Slater determinant (JSD) ansatz, withone-body molecular orbitals (MOs) computed by DFT foreach interatomic distance R. The JSD is the most commonVMC ansatz: only the Jastrow factor is optimized at the VMClevel, while the DFT MOs are kept frozen in the Slater deter-minant (SD). However, in this case, the VMC force F VMC isnot consistent with the slope of EJSD, because the variationalparameters included in the SD are not at their VMC minima.Instead, F VMC corresponds to the slope of EbiasedJSD [Fig. 1(d)],where the DFT MOs obtained at R = R′ are used artificiallyfor all R, such that d pi/dR = 0, ∀i. In this work, we proposean efficient method to obtain atomic forces and pressures thatare unbiased, namely consistent with the slope of EJSD(c).As long as the Jastrow factor is at its variational minimum,the contribution to the bias comes only from the SD part. Thissuggests that a straightforward solution for correcting the biasis to compute the variational term Fcα ≡ −∑NSDpi=1∂E∂ pSDid pSDidRα,where pSDi are SD variational parameters. In the following,we introduce a method to evaluate these terms by combiningDFT with VMC gradients calculations.III. METHOD TO OBTAIN UNBIASED ATOMIC FORCESWe begin by introducing the AGP representation [38,41]of the SD ansatz made of MOs. The general AGPansatz for a system of Ne electrons is written as �AGP =Â[g(x1, x2)g(x3, x4) . . . g(xNe−1, xNe )], where  is the an-tisymmetrization operator and g is the so-called geminalfunction g(xl , xm) = f (rl , rm)(|↑↓〉 − |↓↑〉)/√2. The spatialpart f (r, r′) can be written in terms of MOs, such thatf (r, r′) = ∑Mk �k (r)�k (r′), where �k (x) is the kth MOexpressed as �k (r) = ∑Li ci,kψi(r), ψi(r) is the ith atomicorbital (AO), and ci,k are the AO coefficients obtained bya DFT calculation. If the Hilbert space is restricted to theoccupied states, i.e., M = Ne/2 for spin-unpolarized systems[42], the resultant AGP is equivalent to the SD ansatz. In otherwords, the SD ansatz can be treated as a special case of themore general AGP wave function. We assume that �T is realfor the sake of conciseness; thus, the variational parametersare also real. However, our method can be readily generalizedto complex �T [43]. In this work, the geminal function is con-structed from the MOs obtained from a DFT calculation, andthen converted to the AO representation, namely f (rl , rm) =∑L,Li, j λi, jψi(rl )ψ j (rm), with λi, j = ∑k ci,kc j,k . Thus, the vari-ational term needed for correcting the self-consistency errorreadsFcα = −L,L∑i, j∂E∂λi, jdλi, jdRα, (3)which is what one should compute to get unbiased atomicforces in the JSD ansatz, where λi, j are directly obtained by205151-2EFFICIENT CALCULATION OF UNBIASED ATOMIC … PHYSICAL REVIEW B 109, 205151 (2024)DFT calculations. As discussed later, the geminal representa-tion also allows one to compute unbiased forces and pressuresbeyond the JSD ansatz by optimizing a part of λi, j in the JAGPansatz at the VMC level.The first factor in the terms summed in Eq. (3), i.e.,∂E/∂λi, j , used for optimizing �T and often dubbed as gen-eralized force, can be efficiently computed by VMC [43]. Thesecond factor, the total derivative dλi, j/dRα , can be numer-ically evaluated using the finite-difference method (FDM),i.e., dλi, j/dRα ∼ (λRα+�Rαi, j − λRα−�Rαi, j )/2�Rα , or can be ob-tained by solving the coupled perturbed Hartree-Fock (CPHF)or Kohn-Sham (CPKS) equations [44], or the linear responseequations [45]. The second factor is N times more timeconsuming than the single-point DFT calculation, and thisis regardless of the number of variational parameters. In-deed, to compute the correction terms for the geometry R ≡{R1, . . . , RN }, one needs at least 3N-times HF/DFT calcula-tions, where 3 is the number of Cartesian components. In thisstudy, we employed the finite-difference approach becausethe gauge-invariant property of the AGP, inherited from itsclose relation with the reduced one-body density matrix [46],allows one to construct a robust workflow to compute thesecond factor. Indeed, thanks to the gauge-invariant propertyof the AGP, one does not suffer from (i) the global phase(or sign) indetermination of MOs, nor from (ii) their possibledegeneracy. As for (i), a global phase θ rotating the kth MO(�k → eiθ�k), which is reduced to a global sign eiθ = ±1 incase of a real �T, does not affect the total energy, but preventsthe calculation of orbital derivatives dci,k/dRα , based on finitedifferences. Indeed, the global phase (or sign) is sometimesinconsistent between DFT outcomes with different atomicdisplacements. Instead, the sign flip is not problematic in theAGP representation, because the relation λi, j = ∑k ci,kc j,kimplies that λi, j is invariant under a MO sign change [47].As for (ii), when two (or more) MOs are degenerate, cRα+�Rαi,kand cRα−�Rαi,k might have very different values due to thepresence of the other degenerate MOs. Nevertheless, it isstraightforward to show that a MOs degeneracy does notaffect the uniqueness of λi, j , making λi, j independent of thechoice of the particular DFT implementation for degenerateMOs. Thus, by exploiting the AO representation of the AGPwavefunction, one can always devise a well-defined method tocompute dλi, j/dRα , which will be superior to the calculationof dci,k/dRα .By combining the first and second factor in Eq. (3), thevariational term can be cast in a form suitable for a VMCestimate, as follows [48]:Fcα = −2〈EL(x)L,L∑i, j[(Oi, j (x) − Ōi, j )dλi, jdRα]〉, (4)where we made apparent the dependence of the local operatorson the total electronic coordinate x, sampled by VMC, todistinguish them from constant values. In Eq. (4), Oi, j (x) =∂ ln �T (x)/∂λi, j and Ōi, j ∼ 〈Oi, j (x)〉. We remark that Oi, j (x)can be efficiently computed in a VMC calculation using theadjoint algorithmic differentiation [34], and the divergences ofthe generalized forces can be cured by reweighting methods[33,49]. It is extremely important that the variational termis evaluated in a covariance form of random variables toreduce its fluctuations [39,50]. In addition, the expression inEq. (4) implies that if the variational wave function is an exacteigenstate of the Hamiltonian, Fcα vanishes regardless of theVMC sample, because the local energy coincides with the cor-responding eigenvalue E . Indeed, the zero-variance propertyholds in this expression, which is another way to recover theHellmann-Feynman theorem.IV. APPLICATIONS TO H2 AND Cl2 MOLECULESWe determine the interatomic force of the H2 and Cl2molecules, taken as first examples to assess the accuracyof our method. The ccECPs [51–54] accompanied with theuncontracted cc-pVDZ basis sets were employed for H2 andCl2 molecules. For Cl2, the He-core ccECP was employed.DFT MOs were prepared by PYSCF v2.0.1 [55,56] with theLDA-PZ exchange-correlation functional [57], and then con-verted to the TURBORVB wave function format [41] usingthe TURBOGENIUS package [58] via TREX-IO files [59]. Theinhomogeneous one-body, the two-body, and the three-bodyJastrow factors [41] were added to the SD with frozen DFTMOs and optimized using the linear method [39,40] imple-mented in TURBORVB [41]. The second factor, dλi, j/dRα ,was numerically evaluated using the displacements �R =±0.001 Å along the molecular bond direction.The simple H2 molecule highlights the importance of re-moving the self-consistency error in the forces calculation byadding the variational force term to the regular VMC expres-sion. In H2, the JSD ansatz with DFT MOs is, in principle,exact if the Jastrow factor is converged in the basis set [60].Indeed, the wave function is nodeless, so the difficulty offinding the optimal variational state can be fully transferredto the Jastrow factor determination. Thus, H2 allows one tostudy different situations, from a poor to a refined Jastrowfactor. In this study, we examined a small [1s] and a large[4s2p1d] basis set expansion, as a poor and refined Jastrowfactor, respectively. For the former, Fig. 2(a) shows that theDFT parameters are not optimal at the VMC level; thus, theself-consistency error is present. The equilibrium distanceobtained from the PES [0.7344(2) Å] and the one from regularVMC forces [0.7392(1) Å] are reported in Table I. Figure 2and Table I demonstrate that the self-consistency error ismitigated by the proposed force correction Fc, which givesa bond distance of 0.7341(1) Å, compatible with the onederived from the PES. Figure 2(b) shows that in the case ofa refined Jastrow factor, the self-consistency error is insteadnegligible, because the larger Jastrow expansion compensatesfor the DFT determinant, and all variational parameters areoptimal. Thus, the regular VMC force is already consistentwith the derivative of the PES, and the corresponding forcecorrection eventually vanishes, as reported in Table I. The H2example is illustrative of the capability of the variational termin Eq. (3) to correct the force bias due not only to the frozenDFT MOs, but also to an underconverged Jastrow factor.Figure 2(c) shows the PES of the Cl2 molecule, as yieldedby a [3s1p] Jastrow basis set. Table I reports the equilib-rium geometries obtained from the PES, regular VMC force,and corrected force. The figure and table show that the self-consistency error is more significant for Cl2 (Zeff = 15) than205151-3NAKANO, CASULA, AND TENTI PHYSICAL REVIEW B 109, 205151 (2024)FIG. 2. (a) and (b): H2 PESs (solid green curves), their numerical derivatives (dashed green curves), regular VMC forces (red diamonds),and corrected forces (purple squares) obtained with (a) small [1s] and (b) large [4s2p1d] Jastrow basis sets. The PESs and forces are computedfrom 0.30 to 2.00 Å with 18 equally spaced datapoints plus 5 additional datapoints (0.55 Å, 0.65 Å, 0.75 Å, 0.85Å, and 0.95 Å). The verticaldashed lines represent equilibrium bond lengths obtained by fitting the PES (forces) with a polynomial of 11th (10th) order. (c): Cl2 PES (solidgreen curve), its numerical derivative (dashed green curve), regular VMC forces (red diamonds), and corrected forces (purple squares). ThePES and forces are computed from 1.50 to 2.80 Å with 14 equally spaced datapoints. The vertical dashed lines represent equilibrium bondlengths obtained by fitting the PES (forces) with a polynomial of 6th (5th) order for energies (forces). In all panels, only the region in thevicinity of the equilibrium geometry is drawn. The plotted forces are Fx acting on the left atom of each dimer, where the x axis is aligned withthe direction of the molecular bond.H2 (Zeff = 1). This is consistent with the seminal work byTiihonen et al. [37], reporting that the self-consistency errorincreases with the effective nuclear charge. Figure 2 and Ta-ble I illustrate that the proposed force correction works alsofor heavier molecules.V. APPLICATION TO CUBIC BORON NITRIDENot only atomic forces, but also pressures can be correctedin solids using the same method, just by replacing dλi, j/dRαwith dλi, j/dV . To demonstrate it, we computed the cBNEOS. The ccECPs [51–54] with accompanying uncontractedcc-pVDZ basis sets were used for the cBN calculation. TheTABLE I. The equilibrium bond distances req (Å) of the H2 andCl2 molecules obtained from the PESs, the regular VMC force, andthe corrected force. The corresponding PESs are shown in Fig. 2.Dimers Source req (Å)H2 (Jas. [1s]) PES 0.7344(2)VMC force 0.7392(1)Corrected force 0.7341(1)Experiment 0.741aH2 (Jas.[4s2p1d]) PES 0.7418(3)VMC force 0.7408(6)Corrected force 0.7408(6)Experiment 0.741aCl2 (Jas.[3s1p]) PES 1.987(1)VMC force 1.9979(1)Corrected force 1.9864(1)Experiment 1.987 aaThese values are taken from Ref. [61].linear dependency of the basis sets is solved at the DFT levelby cutting basis set elements with exponents smaller than 0.20a.u. This is crucial to suppress the statistical errors on atomicforces and pressures for periodic systems in QMC calculations[62]. The 2×2×2 conventional supercell (256 valence elec-trons in the simulation cell) with k =  was employed. DFTMOs were prepared by the built-in DFT module implementedin TURBORVB [41] with the LDA-PZ exchange-correlationfunctional [57]. Then, the inhomogeneous one-body, the two-body, and the three-body Jastrow factors [41] were added tothe SD with frozen DFT molecular orbitals. [3s1p] Jastrowbasis sets were employed for B and N atoms. The Jastrowfactor was optimized using the linear method [39,40] im-plemented in TURBORVB [41] for each volume. The secondfactor, dλi, j/dV , was numerically evaluated using the built-inDFT module with volume variations �V = ±0.3%. Figure 3shows the cBN EOS, its volume derivative, the regular VMCpressures, and the corrected pressures. The obtained equilib-rium lattice parameters and volumes are reported in Table II.It is apparent that the self-consistency error in pressure is∼5 GPa, constant over the whole volume range. Our methodgives corrections that bring the estimated pressures very closeto the exact values for all volumes, as shown in Fig. 3. Thisresult illustrates the possibility to successfully correct not onlyatomic forces but also pressures in large systems with theexplicit evaluation of the variational pressure term.VI. DISCUSSIONWe first compare our method with the FDM, which is thetraditional way to obtain unbiased atomic forces in the VMCframework. The main drawback of the FDM is that it requiresat least 3N independent VMC runs to compute all 3N force205151-4EFFICIENT CALCULATION OF UNBIASED ATOMIC … PHYSICAL REVIEW B 109, 205151 (2024)FIG. 3. cBN EOS (solid green curve), its volume derivative(dashed green curve), the regular VMC pressure (red diamonds),and the corrected pressure (purple squares). The vertical dashedlines represent equilibrium volumes obtained by fitting the EOS andpressures with the Vinet forms [63].components, preventing its use in routine VMC calculations.Instead, our proposed method requires just a single VMC runto compute all 3N regular VMC forces, together with all Oi, jterms that appear in the expression for Fcα of Eq. (4). Thisis thanks to the algorithmic differentiation [34]. As we men-tioned before, the other terms in Eq. (4), namely dλi, j/dRα ,are computed by FDM using DFT as the driver, thus leadingto DFT calculations N times more time consuming than asingle DFT run. However, since the DFT cost is negligiblecompared to VMC, and it is mainly fast Fourier transformbound, with a favorable O(N2 log N ) scaling for a single run,the resulting algorithmic cost of our method is superior to theFDM evaluation of atomic VMC forces.Next, we discuss the scaling of the variance of the varia-tional term Fcα with respect to N . The variance of the localenergy EL scales with Ne [43], while the variance of thelogarithmic derivatives Oi, j is O(1) [34]. Thus, the variance ofFcα is bound by O(L2Ne ), where the factor of L2 comes fromthe double summation over the extended basis set elements[64]. The geminal representation needs the L2 summationinstead of the LNe summation of the SD representation forthe JSD ansatz. However, at variance with the SD represen-TABLE II. Equilibrium lattice parameters and volumes per atomobtained by fitting the EOS, and from the regular VMC pressure andthe corrected one. Zero-point energy and temperature effects are notincluded.Source Lattice (Å) Volume (Bohr3)EOS 3.5962(3) 39.232(9)VMC pressure 3.5800(1) 38.704(5)Corrected pressure 3.5943(2) 39.169(5)Experiment 3.594a 39.160aaThese values are taken from Ref. [61].tation, the geminal allows one to exploit the locality of theλi, j matrix. In other words, one can neglect |dλi, j/dRα| withsmall absolute values, obtained deterministically by DFT cal-culations. For instance, the percentage of elements such that|dλi, j/dRα|/ max |dλi, j/dRα| � 0.01% is 38.5%, 45.0%, and66.0% for 1 × 1 × 1 (8 atoms), 2 × 2 × 2 (64 atoms), and 3 ×3 × 3 (216 atoms) cBN supercells, respectively, demonstrat-ing that the larger a system becomes, the more terms can beneglected thanks to the locality. In this way, the summation inEq. (4) can be reduced from L2 to L terms, by lowering the sizescaling of the Fcα variance. Since L and Ne are proportional toN , the scaling of the variance of our method with respect to Nis bound by O(N2) in the N → ∞ limit, which is just N timeslarger than the variance of the regular VMC force calculationO(N ) [34]. However, in our H2, Cl2, and cBN calculations,we got the same error bars on the regular VMC forces andon the corrected ones with the same statistics. This pointsto a very small prefactor ε in the O(N2) variance term, suchthat in the total variance Var(Fα ) = Var(FVMCα ) + Var(Fcα ) ≈O(N ) + εO(N2), the O(N2) contribution can be neglected forany affordable N in VMC calculations.Finally, we emphasize the extensibility of the geminalrepresentation employed here, which allows one to readilygeneralize the method proposed in this work from the JSD tothe more general JAGP ansatz. A practically way to go beyondthe JSD ansatz for a large system is to optimize only a subsetof the variational λi, j parameters. The partially optimizedλi, j matrix will normally have a larger rank than the onecorresponding to the SD wave function, therefore includingAGP correlations. The subset of λi, j is chosen again based onthe AGP locality. Indeed, only the variational parameters λi, jcorresponding to atoms at a distance smaller than a reasonablecutoff can be optimized, while those with distance larger thanthe cutoff are kept fixed [43]. In this situation, only the fixedλi, j must enter in Eq. (3), thus correcting the force bias in theJAGP ansatz. In principle, our approach can also be extendedto more general antisymmetric wave functions, with the onlycaveat that, in order to compute the d pi/dRα derivatives,one has to consistently use the same auxiliary frameworkemployed to initialize the antisymmetric part of the VMCwave function.VII. CONCLUDING REMARKSIn this work, we analyzed the bias seriously affecting theregular VMC expression FVMCα . We then proposed a method toefficiently and robustly compute the missing contribution, i.e.,the variational term Fcα , to completely remove that bias for aJSD ansatz with DFT one-body orbitals, the most commonwave function in ab initio VMC calculations, usually thebest compromise between accuracy and computational cost.We demonstrated that the correction works very well for thesystems that have been tested here, namely the equilibriumgeometry of H2 and Cl2 molecules, and the EOS evaluationof the cubic boron nitride. Unbiased atomic forces within aJSD ansatz, which is in general much cheaper to optimize thanmore refined �Ts, will be particularly useful to generate VMCdatasets for MLPs construction, which would otherwise beaffected by the self-consistency error. Thus, our approach hasthe potential to open up new horizons for VMC applications,205151-5NAKANO, CASULA, AND TENTI PHYSICAL REVIEW B 109, 205151 (2024)also in the context of machine learning. Finally, the samescheme can be extended to more elaborated wave functions,once a suitable auxiliary method is used to generate theirantisymmetric part.ACKNOWLEDGMENTSK.N. is grateful for computational resources from theNumerical Materials Simulator at National Institute forMaterials Science (NIMS). The authors are grateful for com-putational resources of the supercomputer Fugaku providedby RIKEN through the HPCI System Research Projects(Project No. hp230030). K.N. acknowledges financial sup-port from Grant-in-Aid for Early Career Scientists (GrantNo. JP21K17752), from Grant-in-Aid for Scientific Re-search (Grant No. JP21K03400), and from MEXT LeadingInitiative for Excellent Young Researchers (Grant No. JP-MXS0320220025). K.N. is grateful for fruitful discussionwith Dr. Atushi Togo (NIMS) about the illustration of theself-consistency problem. The authors acknowledge valuablecomments from Dr. Jürg Hutter (UZH), Dr. Stefano Battaglia(UZH) and Dr. Emmanuel Giner (CNRS), Dr. Saverio Moroni(SISSA), and Dr. Claudia Filippi (UT). M.C. acknowledgesGENCI for granting access to French computational resourcesthrough the DARI Application No. 906493. The molecularand crystal structures were depicted using VESTA [65]. Thiswork is supported by the European Centre of Excellence inExascale Computing TREX - Targeting Real Chemical Accu-racy at the Exascale. This project has received funding fromthe European Unions Horizon 2020 Research and Innovationprogram, under Grant Agreement No. 952165. The ab initioQMC package used in this work, TURBORVB, is availablefrom its GitHub repository [66].[1] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev.Mod. Phys. 73, 33 (2001).[2] A. Zen, J. G. Brandenburg, J. Klimeš, A. Tkatchenko, D. Alfè,and A. Michaelides, Proc. Natl. Acad. Sci. USA 115, 1724(2018).[3] E. Mostaani, N. D. Drummond, and V. I. Fal’ko, Phys. Rev. Lett.115, 115501 (2015).[4] T. Frank, R. Derian, K. Tokár, L. Mitas, J. Fabian, and I. Štich,Phys. Rev. X 9, 011018 (2019).[5] Y. Nikaido, T. Ichibha, K. Hongo, F. A. Reboredo, K. H. Kumar,P. Mahadevan, R. Maezono, and K. Nakano, J. Phys. Chem. C126, 6000 (2022).[6] M. Casula and S. Sorella, Phys. Rev. B 88, 155125 (2013).[7] R. C. Clay, J. Mcminis, J. M. McMahon, C. Pierleoni, D. M.Ceperley, and M. A. Morales, Phys. Rev. B 89, 184106(2014).[8] R. C. Clay, M. Holzmann, D. M. Ceperley, and M. A. Morales,Phys. Rev. B 93, 035121 (2016).[9] 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).[10] G. Mazzola, R. Helled, and S. Sorella, Phys. Rev. Lett. 120,025701 (2018).[11] L. Monacelli, M. Casula, K. Nakano, S. Sorella, and F. Mauri,Nat. Phys. 19, 845 (2023).[12] L. K. Wagner, M. Bajdich, and L. Mitas, J. Comput. Phys. 228,3390 (2009).[13] A. Scemama, M. Caffarel, E. Oseret, and W. Jalby, J. Comput.Chem. 34, 938 (2013).[14] R. Needs, M. Towler, N. Drummond, P. Lopez Rios, and J.Trail, J. Chem. Phys. 152, 154106 (2020).[15] P. R. C. Kent, A. Annaberdiyev, A. Benali, M. C. Bennett,E. J. Landinez Borda, P. Doak, H. Hao, K. D. Jordan, J. T.Krogel, I. Kylänpää, J. Lee, Y. Luo, F. D. Malone, C. A. Melton,L. Mitas, M. A. Morales, E. Neuscamman, F. A. Reboredo,B. Rubenstein, K. Saritas et al., J. Chem. Phys. 152, 174105(2020).[16] W. A. Wheeler, S. Pathak, K. G. Kleiner, S. Yuan, J. N. B.Rodrigues, C. Lorsung, K. Krongchon, Y. Chang, Y. Zhou, B.Busemeyer, K. T. Williams, A. Muñoz, C. Y. Chow, and L. K.Wagner, J. Chem. Phys. 158, 114801 (2023).[17] R. M. Martin, Electronic structure: basic theory and practicalmethods (Cambridge University Press, Cambridge, 2004).[18] R. J. DiRisio, F. Lu, and A. B. McCoy, J. Phys. Chem. A 125,5849 (2021).[19] A. Tirelli, G. Tenti, K. Nakano, and S. Sorella, Phys. Rev. B106, L041105 (2022).[20] C. Huang and B. M. Rubenstein, J. Phys. Chem. A 127, 339(2023).[21] H. Niu, Y. Yang, S. Jensen, M. Holzmann, C. Pierleoni, andD. M. Ceperley, Phys. Rev. Lett. 130, 076102 (2023).[22] D. Ceperley, S. Jensen, Y. P. Yang, H. Niu, C. Pierleoni, and M.Holzmann, Electron. Struct. 6, 015011 (2024).[23] P. J. Reynolds Jr, R. Barnett, B. Hammond, R. Grimes, and W.Lester Jr, Int. J. Quantum Chem. 29, 589 (1986).[24] R. Assaraf and M. Caffarel, J. Chem. Phys. 113, 4028 (2000).[25] C. Filippi and C. J. Umrigar, Phys. Rev. B 61, R16291 (2000).[26] S. Chiesa, D. M. Ceperley, and S. Zhang, Phys. Rev. Lett. 94,036404 (2005).[27] A. Badinski and R. J. Needs, Phys. Rev. B 78, 035134 (2008).[28] R. Assaraf, M. Caffarel, and A. C. Kollias, Phys. Rev. Lett. 106,150601 (2011).[29] S. Moroni, S. Saccani, and C. Filippi, J. Chem. Theory Comput.10, 4823 (2014).[30] J. van Rhijn, C. Filippi, S. De Palo, and S. Moroni, J. Chem.Theory Comput. 18, 118 (2022).[31] R. Assaraf and M. Caffarel, J. Chem. Phys. 119, 10536 (2003).[32] C. J. Umrigar, Int. J. Quantum Chem. 36, 217 (1989).[33] C. Attaccalite and S. Sorella, Phys. Rev. Lett. 100, 114501(2008).[34] S. Sorella and L. Capriotti, J. Chem. Phys. 133, 234111 (2010).[35] C. Filippi, R. Assaraf, and S. Moroni, J. Chem. Phys. 144,194105 (2016).[36] K. Nakano, A. Raghav, and S. Sorella, J. Chem. Phys. 156,034101 (2022).[37] J. Tiihonen III, R. C. Clay III, and J. T. Krogel, J. Chem. Phys.154, 204111 (2021).[38] M. Casula and S. Sorella, J. Chem. Phys. 119, 6500 (2003).205151-6https://doi.org/10.1103/RevModPhys.73.33https://doi.org/10.1073/pnas.1715434115https://doi.org/10.1103/PhysRevLett.115.115501https://doi.org/10.1103/PhysRevX.9.011018https://doi.org/10.1021/acs.jpcc.1c10943https://doi.org/10.1103/PhysRevB.88.155125https://doi.org/10.1103/PhysRevB.89.184106https://doi.org/10.1103/PhysRevB.93.035121https://doi.org/10.1038/ncomms8794https://doi.org/10.1103/PhysRevLett.120.025701https://doi.org/10.1038/s41567-023-01960-5https://doi.org/10.1016/j.jcp.2009.01.017https://doi.org/10.1002/jcc.23216https://doi.org/10.1063/1.5144288https://doi.org/10.1063/5.0004860https://doi.org/10.1063/5.0139024https://doi.org/10.1021/acs.jpca.1c03709https://doi.org/10.1103/PhysRevB.106.L041105https://doi.org/10.1021/acs.jpca.2c05904https://doi.org/10.1103/PhysRevLett.130.076102https://doi.org/10.1088/2516-1075/ad2eb0https://doi.org/10.1002/qua.560290403https://doi.org/10.1063/1.1286598https://doi.org/10.1103/PhysRevB.61.R16291https://doi.org/10.1103/PhysRevLett.94.036404https://doi.org/10.1103/PhysRevB.78.035134https://doi.org/10.1103/PhysRevLett.106.150601https://doi.org/10.1021/ct500780rhttps://doi.org/10.1021/acs.jctc.1c00496https://doi.org/10.1063/1.1621615https://doi.org/10.1002/qua.560360826https://doi.org/10.1103/PhysRevLett.100.114501https://doi.org/10.1063/1.3516208https://doi.org/10.1063/1.4948778https://doi.org/10.1063/5.0076302https://doi.org/10.1063/5.0052266https://doi.org/10.1063/1.1604379EFFICIENT CALCULATION OF UNBIASED ATOMIC … PHYSICAL REVIEW B 109, 205151 (2024)[39] S. Sorella, Phys. Rev. B 71, 241103(R) (2005).[40] C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G.Hennig, Phys. Rev. Lett. 98, 110201 (2007).[41] 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).[42] The AGP ansatz has been generalized also to spin-polarizedsystems [38].[43] F. Becca and S. Sorella, Quantum Monte Carlo approaches forcorrelated systems (Cambridge University Press, Cambridge,2017).[44] F. Jensen, Introduction to computational chemistry (John wiley& sons, Hoboken, NJ, 2017).[45] J. Toulouse, Introduction to the calculation of molecular prop-erties by response theory, https://hal.science/hal-03934866/document (2018).[46] O. Goscinski, Int. J. Quantum Chem. 22, 591 (1982).[47] λi, j is invariant under a unitary transformation of MOs only ifλi, j = ∑k c∗i,kc j,k in complex cases.[48] Only its real part is taken in complex cases [43].[49] S. Pathak and L. K. Wagner, AIP Adv. 10, 085213 (2020).[50] C. J. Umrigar and C. Filippi, Phys. Rev. Lett. 94, 150201(2005).[51] M. C. Bennett, C. A. Melton, A. Annaberdiyev, G. Wang,L. Shulenburger, and L. Mitas, J. Chem. Phys. 147, 224106(2017).[52] M. C. Bennett, G. Wang, A. Annaberdiyev, C. A. Melton,L. Shulenburger, and L. Mitas, J. Chem. Phys. 149, 104108(2018).[53] A. Annaberdiyev, G. Wang, C. A. Melton, M. ChandlerBennett, L. Shulenburger, and L. Mitas, J. Chem. Phys. 149,134108 (2018).[54] G. Wang, A. Annaberdiyev, C. A. Melton, M. C. Bennett, L.Shulenburger, and L. Mitas, J. Chem. Phys. 151, 144110 (2019).[55] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo,Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma,S. Wouters, and G. K.-L. Chan, WIREs Comput. Mol. Sci. 8,e1340 (2018).[56] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt,N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen,Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval,S. Lehtola, Z. Li, J. Liu et al., J. Chem. Phys. 153, 024109(2020).[57] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).[58] K. Nakano, O. Kohulák, A. Raghav, M. Casula, and S. Sorella,J. Chem. Phys. 159, 224801 (2023).[59] E. Posenitskiy, V. G. Chilkuri, A. Ammar, M. Hapka, K. Pernal,R. Shinde, E. J. Landinez Borda, C. Filippi, K. Nakano, O.Kohulák, S. Sorella, P. de Oliveira Castro, W. Jalby, P. L. Ríos,A. Alavi, and A. Scemama, J. Chem. Phys. 158, 174801 (2023).[60] See Ref. [41] for the detail of the Jastrow factor implementationin the ab initio QMC package, TURBORVB, used in this study.[61] K.-P. Huber, Molecular spectra and molecular structure: IV.Constants of diatomic molecules (Springer Science & BusinessMedia, New York, NY, 2013).[62] K. Nakano, T. Morresi, M. Casula, R. Maezono, and S. Sorella,Phys. Rev. B 103, L121110 (2021).[63] P. Vinet, J. R. Smith, J. Ferrante, and J. H. Rose, Phys. Rev. B35, 1945 (1987).[64] In the variance estimation, covariance effects are neglected.Nevertheless, in the actual calculations, the full statistical errorwas estimated by the Jackknife method [43].[65] K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).[66] TURBORVB on GitHub, https://github.com/sissaschool/turborvb.205151-7https://doi.org/10.1103/PhysRevB.71.241103https://doi.org/10.1103/PhysRevLett.98.110201https://doi.org/10.1063/5.0005037https://hal.science/hal-03934866/documenthttps://doi.org/10.1002/qua.560220851https://doi.org/10.1063/5.0004008https://doi.org/10.1103/PhysRevLett.94.150201https://doi.org/10.1063/1.4995643https://doi.org/10.1063/1.5038135https://doi.org/10.1063/1.5040472https://doi.org/10.1063/1.5121006https://doi.org/10.1002/wcms.1340https://doi.org/10.1063/5.0006074https://doi.org/10.1103/PhysRevB.23.5048https://doi.org/10.1063/5.0179003https://doi.org/10.1063/5.0148161https://doi.org/10.1103/PhysRevB.103.L121110https://doi.org/10.1103/PhysRevB.35.1945https://doi.org/10.1107/S0021889811038970https://github.com/sissaschool/turborvb