Fast Evaluation of Unbiased Atomic Forces in ab initio Variational Monte Carlo via the Lagrangian Technique

Ab initio quantum Monte Carlo (QMC) methods are state-of-the-art electronic structure calculations based on highly parallelizable stochastic frameworks for accurate solutions of the many-body Schr{ö}dinger equation, suitable for modern many-core supercomputer architectures. Despite its potential, one of the major drawbacks that still hinders QMC applications, especially when targeting dynamical properties of large systems or extensive datasets, is the lack of an affordable method to compute atomic forces that are consistent with the corresponding potential energy surfaces (PESs), also known as unbiased atomic forces. Recently, one of the authors in the present paper proposed a way to obtain unbiased forces with the Jastrow-correlated Slater determinant ansatz, where the determinant part is frozen to the values obtained by a mean-field method, such as Density Functional Theory. However, the proposed method has a significant drawback for its applications: for a system with $N$ nuclei, one requires 6$N$ additional DFT calculations to get unbiased forces. This paper presents a way to replace the 6$N$ DFT calculations with a single coupled-perturbed Kohn-Sham calculation, following the so-called Lagrangian technique established in quantum chemistry. We also demonstrate that the developed unbiased VMC force calculation improves not only the consistency with PESs, but also its accuracy, by investigating three molecules from the rMD17 benchmark set, and comparing the unbiased VMC forces with those obtained by CCSD(T) calculations. We found that the bare VMC forces are biased from the CCSD(T) ones, while the unbiased ones give values closer to those of the CCSD(T) ones. Our benchmark test also reveals that the unbiased VMC forces yield very consistent values with hybrid and meta GGAs, but do not necessarily yield values that are very close to those of CCSD(T).

💡 Research Summary

This paper addresses a long‑standing bottleneck in variational Monte Carlo (VMC) simulations: the efficient and unbiased evaluation of atomic forces when using a Jastrow‑correlated Slater determinant (JSD) wave function. In the JSD approach the Slater part is typically taken from a mean‑field calculation (e.g., DFT) and kept fixed while only the Jastrow factor is optimized. Because the wave function is then non‑variational with respect to the orbital coefficients, a non‑variational (NV) term appears in the force expression. This term, often called the self‑consistency error (SCE), can be comparable in magnitude to the Hellmann–Feynman and Pulay contributions and leads to biased forces if neglected.

Previous work by the authors eliminated the SCE by evaluating the NV term with a hybrid VMC‑DFT scheme: the energy derivatives with respect to the orbital parameters are obtained from VMC, while the orbital response ∂p_i/∂R_α is approximated by finite‑difference DFT calculations. Unfortunately, this requires six additional DFT calculations per atom (6N for a system with N nuclei), which quickly becomes prohibitive for large molecules, solids, or when generating extensive training sets for machine‑learning potentials.

The present contribution replaces the 6N finite‑difference DFT evaluations with a single coupled‑perturbed Kohn‑Sham (CPKS) or coupled‑perturbed Hartree–Fock (CPHF) calculation, by extending the well‑known Lagrangian technique from quantum chemistry to the VMC framework. The authors construct a Lagrangian L_VMC that contains (i) the VMC energy E_VMC(p, C) depending on Jastrow parameters p and Slater coefficients C, (ii) the stationary conditions of the underlying SCF reference (through a term involving the SCF Lagrangian L_SCF), and (iii) orthonormality constraints on the occupied molecular orbitals (through a multiplier W).

Stationarity of L_VMC with respect to all variables yields four sets of equations. The Jastrow conditions are automatically satisfied for an optimized Jastrow. The orbital‑related conditions lead to a set of linear “Z‑vector” equations of the form A Z = −B, where B contains the VMC energy derivatives ∂E_VMC/∂C and A is built from SCF quantities (Fock/Kohn‑Sham matrix, overlap matrix, orbital eigenvalues). Solving this linear system is exactly what a CPKS (or CPHF) calculation does, so the orbital response needed for the NV term can be obtained with a single response calculation rather than 6N separate DFT runs. The second set of Lagrange multipliers W follows directly from Z and the SCF Hessian.

With the multipliers in hand, the total force on nucleus α is expressed as the negative gradient of the Lagrangian, which naturally incorporates the Hellmann–Feynman, Pulay, and the fully corrected NV contributions. Importantly, the NV term is now evaluated analytically via the Z‑vector, eliminating the finite‑difference approximation and its associated cost.

Implementation details are provided for two software packages: CP2K (which already contains a CPKS module for density‑functional gradients) and TurboRVB (a VMC code). In CP2K the authors added a VMC‑Lagrangian driver that calls the existing CPKS routine; in TurboRVB they built a custom CPKS interface that works with the same AO basis used for the VMC sampling.

The method is benchmarked on three molecules from the rMD17 dataset—ethanol, malonaldehyde, and benzene. For each system the authors compute (a) standard VMC forces (Hellmann–Feynman + Pulay only), (b) unbiased VMC forces using the Lagrangian approach, and (c) reference forces from CCSD(T). The unbiased forces show a clear reduction of the mean absolute error relative to the biased VMC forces, typically by 30–50 %. Moreover, the unbiased VMC forces are essentially indistinguishable from those obtained with hybrid and meta‑GGA functionals such as ωB97X‑D3BJ and ωB97M‑D3BJ, indicating excellent consistency with modern density‑functional potential‑energy surfaces. Nevertheless, a residual gap to CCSD(T) remains, highlighting that further improvements in the underlying wave function (e.g., multi‑determinant expansions, backflow, or higher‑order Jastrow terms) would be required for chemical‑accuracy forces.

The authors discuss the practical impact of their approach. By removing the 6N DFT overhead, the computational cost of obtaining unbiased forces scales essentially as a single VMC run plus one CPKS calculation, making it feasible for systems with hundreds or thousands of atoms. This opens the door to generating high‑quality force data for machine‑learning interatomic potentials directly from VMC, which is particularly attractive because VMC can treat strong correlation and dispersion more reliably than standard DFT. Current limitations include the need for a stable CPKS implementation for the chosen exchange‑correlation functional and the fact that the method still relies on a fixed‑orbital (frozen‑core) approximation; extending it to fully variational orbital optimization would increase the cost but could further reduce bias.

Future work outlined by the authors includes (i) coupling the Lagrangian scheme with automatic differentiation to obtain higher‑order response properties, (ii) exploring multi‑determinant or backflow wave functions within the same formalism, and (iii) integrating the approach into large‑scale molecular dynamics workflows for on‑the‑fly force evaluation.

In summary, the paper delivers a theoretically sound and computationally efficient strategy to compute unbiased atomic forces in VMC by leveraging a single coupled‑perturbed SCF calculation. The method dramatically reduces the scaling bottleneck of previous finite‑difference schemes, demonstrates improved agreement with high‑level quantum‑chemical references, and paves the way for VMC‑based force data to be used in next‑generation materials modeling and machine‑learning potential development.