Quantum statistics from classical simulations via generative Gibbs sampling

TL;DR

The paper presents GG-PI, a data-driven framework that recovers equilibrium quantum statistics from classical MD by learning the single-bead conditional density of a ring-polymer and performing Gibbs sampling. By training an $E(3)$-equivariant flow to approximate $p_\tau(\mathbf x_i|\mathbf y_i)$ and updating beads in parallel, GG-PI achieves quantum-like statistics without per-step force evaluations and can transfer along the fixed imaginary-time step $\tau=\beta/P$ to other temperatures without retraining. Validation on Zundel ion, bulk water, and para-H$_2$ shows GG-PI reproduces key observables such as RDFs and radii of gyration while offering substantial speedups over traditional PIMD, especially for systems with expensive potentials. The approach generalizes to other problems with similar Markov structure and holds promise for broader applications in transition-path sampling and quantum simulations, with future work extending to indistinguishable particles and PIGS-style open-path formulations.

Abstract

Accurate simulation of nuclear quantum effects is essential for molecular modeling but expensive using path integral molecular dynamics (PIMD). We present GG-PI, a ring-polymer-based framework that combines generative modeling of the single-bead conditional density with Gibbs sampling to recover quantum statistics from classical simulation data. GG-PI uses inexpensive standard classical simulations or existing data for training and allows transfer across temperatures without retraining. On standard test systems, GG-PI significantly reduces wall clock time compared to PIMD. Our approach extends easily to a wide range of problems with similar Markov structure.

Figures (9)

• Figure 1: Zundel cation (H$_5$O$_2^+$) at $300~\mathrm{K}$. Joint distribution $P(\delta,R_{\mathrm{OO}})$ with $\delta=R(\mathrm{O_aH})-R(\mathrm{O_bH})$ and $R_{\mathrm{OO}}$. Panels: (a) standard MD, (b) GG–PI, and (c) reference PIMD. Path–integral simulations use $P=32$; identical contour levels are used across panels.
• Figure 2: Structural properties of bulk liquid water at 300 K, comparing standard MD (orange solid lines), our GG-PI (blue solid lines), and a reference PIMD simulation (green dashed lines). All path-integral simulations are performed with $P=32$. The four panels display: (a) the O-O radial distribution function (RDF), (b) the O-H RDF, (c) the H-H RDF, and (d) the intramolecular H-O-H bond angle distribution.
• Figure 3: Temperature dependence of three observables for para-H$_2$: (a) kinetic energy per molecule, (b) potential energy per molecule, and (c) radius of gyration. PIMD data are simulated at $25~\mathrm{K}$ ($P=32$), $50~\mathrm{K}$ ($P=16$), and $100~\mathrm{K}$ ($P=8$). GG-PI is shown as a purple solid line with circle symbols and PIMD as red square symbols; error bars indicate statistical uncertainties.
• Figure 4: Radial distribution functions for para-H$_2$ generated by GG-PI (blue solid lines) against reference PIMD simulation (green dashed lines) at three different temperatures with the same $\tau$: (a)$T$=25 K, $P=32$; (b)$T$=50 K, $P=16$; (c)$T$=100 K, $P=8$.
• Figure S1: Batch size vs. ESS/sec for GG-PI. The scanning results for (a) the Zundel cation, (b) water, and (c) para-H$_2$. The red stars indicate the optimal batch sizes that yield the highest sampling efficiency on an NVIDIA RTX 5090 GPU.