Variance reduction for forces and pressure in variational Monte Carlo

Abstract

We present simple and practical strategies to reduce the variance of Monte Carlo estimators. Our focus is on variational Monte Carlo calculations of atomic forces and pressure in electronic systems, although we show that the underlying ideas apply more broadly to other observables, like pair-correlation and angular-distribution functions, and other methods, including molecular dynamics. For Pulay-type contributions, we show that a minor modification based on the Metropolis acceptance ratio softens the power-law divergence of the variance to a logarithmic one, and that inexpensive regularizations can further suppress outliers at the price of a controlled small bias. For Hellmann-Feynman forces, we derive compact variance-reduced estimators for periodic systems that are straightforward to implement in standard Monte Carlo codes. The approach is illustrated for high-pressure metallic hydrogen with more than a hundred atoms described by neural quantum states, including an application to molecular dynamics driven by the improved forces.

Figures (6)

• Figure 1: Comparison of Pulay pressure ($P^{wf}$) estimators for a periodic system of $N=128$ hydrogen atoms. The estimators' names and associated variance-reduction factors are reported in the legend. The covariance estimator (denoted also as "Cov.") is given in \ref{['eq:pulay_covariance']}, the acceptance trick (denoted as "acc.") is given in \ref{['eq:acceptance_trick_estimator']}, and the (hard-cutoff) regularized estimator (denoted as "reg.") is given in \ref{['eq:hard_cutoff_regularized_pulay_pressure']}.
• Figure 2: Comparison between the pressure obtained from the direct estimator and from a finite-difference fit. The direct estimate is computed as the sum of the Virial and Pulay contributions, $P_{\mathrm{tot}}(\Omega_0) = P^{\mathrm{V}}(\Omega_0) + P^{wf}(\Omega_0)$ (indicated with a dotted black line), where $\Omega_0$ is the simulation-cell volume of the original (unscaled) nuclear configuration. The finite-difference estimate, denoted $P_{\mathrm{fit}}(\Omega_0)$ (indicated by a large red marker), is obtained by differentiating a polynomial fit to the VMC energy data and evaluating it at $\Omega_0$. In the inset, the original and scaled cells are shown as light green and light orange surfaces, respectively, with solid and dashed lines marking their boundaries. The scaled configuration is illustrated by displaced nuclear positions; the arrows indicate geometric displacements (not forces). The inset is not to scale and is intended only as a schematic, since the actual volume changes considered are too small to be visible.
• Figure 3: Comparison of the Hellmann-Feynman force estimators. The force component for the proton index 0 in the $x$-direction is shown without loss of generality. The estimators "IBP1" and "IBP2" are presented in \ref{['eq:ibp1_estimator', 'eq:ibp2_estimator']}, respectively, while "Previous work" (or "Prev." for short), refers to an estimator suggested in a recent paper 2024_qian_space_warp_forces, that is also described in \ref{['app:ac_pbc_alternative']}.
• Figure 4: Comparison of the Pulay force estimators. The covariance estimator is given in \ref{['eq:pulay_covariance']}, and the acceptance trick is given in \ref{['eq:acceptance_trick_estimator']}. The force component for the proton index 0 in the $x$-direction is shown without loss of generality.
• Figure 5: Molecular dynamics evolution of the hydrogen system, showing from left to right the electronic energy, the atomic forces, and the nuclear configurations. The molecular dynamics protocol consists of an iterative three-step procedure: (1) a short VMC optimization; (2) the evaluation of the forces associated with the optimized wave function and nuclear configuration; and (3) the propagation of the nuclear coordinates using the damped velocity Verlet integration scheme based on these forces, after which the cycle is restarted with the updated coordinates.