Nuclear gradients from auxiliary-field quantum Monte Carlo and their application in geometry optimization and transition state search

TL;DR

This work develops a workflow to compute nuclear gradients in ph-AFQMC via reverse-mode automatic differentiation, achieving gradient costs near energy evaluation with $O(N^4)$ scaling per sample. It validates the gradients against finite difference and uses ML surrogates—most notably $\Delta$-KRR trained on UMA—trained on energies and forces to handle noisy AFQMC data. Applying these ML models to geometry optimization and NEB-based transition-state searches yields transition-state structures and barrier heights in close agreement with CCSD(T). The approach offers a scalable path toward AFQMC-based molecular dynamics and reaction-path simulations, with GPU acceleration and strategies to mitigate SR bias and memory overhead.

Abstract

In this article, we present a method for computing accurate and scalable nuclear forces within the phaseless auxiliary-field quantum Monte Carlo (AFQMC) framework. Our approach leverages automatic differentiation of the energy functional to obtain nuclear gradients at a computational cost comparable to that of energy evaluation. The accuracy of the method is validated against finite difference calculations, showing excellent agreement. We then explore several machine learning (ML) strategies for learning noisy AFQMC data. These ML potentials are subsequently used to perform geometry optimizations and nudged elastic band (NEB) calculations, successfully identifying the transition state of the formamide-formimidic acid tautomerization. The resulting transition state geometry and barrier heights are in close agreement with coupled-cluster reference values. This work paves the way for highly accurate geometry optimization, molecular dynamics, or reaction path calculations.

Figures (5)

• Figure 1: Difference between the nuclear gradients obtained from finite difference (FD) and reverse-mode automatic differentiation (AD) methods for the symmetric stretch of C–H bonds in methane, using different intervals for performing stochastic reconfiguration: (top) with orbital relaxation and (bottom) without orbital relaxation. The grey shaded region represents the stochastic error bars on the FD data.
• Figure 2: Mean absolute errors (MAEs) of (top) energies (in kcal/mol) and (bottom) forces (in kcal/mol/Å) obtained using different ML methods. The parameter $\Delta$ denotes the standard deviation of random Gaussian noise added to the training dataset, and the x-axis indicates the number of training samples used. The left panels show results for models trained on both energies and forces, and the right panels correspond to models trained exclusively on energies. Note that the training set sizes are scaled with $\Delta$ to approximate a fixed QMC computational budget across datasets (see main text).
• Figure 3: Formamide to formimidic acid hydrogen transfer reaction with a representative transition state shown above the arrow.
• Figure 4: Results of the NEB calculation for the formamide–formimidic acid tautomerization reaction. Differences in bond lengths (top) and bond angles (bottom) for all bonds in the transition state, as predicted by DFT (B3LYP) and AFQMC, relative to CCSD(T).
• Figure 5: Scaling of AFQMC energy and rev-AD force evaluations for a 1-D hydrogen chain in the STO-6G basis. Wall-clock timings versus system size are shown for energy-only (orange) and force (blue) calculations, comparing CPU (dashed) and GPU (solid) performance for identical sampling effort. The inset shows the ratio of the wall-clock time for force evaluations to that for energy-only evaluations.