Quantum statistics from classical simulations via generative Gibbs sampling

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.

💡 Research Summary

Accurately capturing nuclear quantum effects (NQEs) is essential for reliable molecular simulations, yet the standard approaches—path‑integral molecular dynamics (PIMD) and path‑integral Monte Carlo (PIMC)—are computationally expensive because the cost scales linearly with the number of beads used to discretize the imaginary‑time path. In this work the authors introduce GG‑PI (Generative‑Gibbs‑Path‑Integral), a framework that recovers quantum statistical distributions from inexpensive classical molecular dynamics (MD) data by learning a generative model of the single‑bead conditional density and then using Gibbs sampling to assemble the full ring‑polymer distribution.

The theoretical foundation rests on the observation that the Euclidean action of a ring polymer is local in bead index: for a given bead i only the on‑site potential τ V(x_i) and the harmonic couplings to its two neighbours appear. Consequently the exact conditional probability p_τ(x_i | x_{‑i}) depends only on the neighbour average y_i = (x_{i‑1}+x_{i+1})/2 and can be written as a Gaussian centered at y_i multiplied by the Boltzmann weight exp(−τ V(x_i)). Because the Gaussian term strongly localizes the bead, a modest‑capacity flow model is sufficient to approximate this conditional distribution.

Training data for the conditional density can be obtained in three ways: (i) a Bayesian construction that samples x_i from a classical MD run at the effective inverse temperature τ and draws y_i from the corresponding Gaussian; (ii) restrained MD where y_i is fixed and x_i is sampled under the potential; (iii) direct extraction of (x_i, y_i) pairs from existing PIMD or PIMC trajectories. The first two approaches require only cheap classical simulations, while the third leverages already‑available quantum data and enables temperature transfer without retraining.

A lightweight equivariant flow model (≈10⁵ parameters) is trained via flow‑matching to approximate p_τ(x_i | y_i). Once trained, multiple independent Markov chains are initialized and updated in parallel using odd‑even Gibbs sweeps: all odd beads are resampled simultaneously from the learned flow, followed by all even beads. Because the flow proposal is already very close to the true conditional, the authors omit the Metropolis‑within‑Gibbs correction in practice, though it could be added for exactness. A Rao‑Blackwellized estimator further reduces variance.

The method is benchmarked on three systems. For the Zundel cation (H₅O₂⁺) at 300 K, GG‑PI reproduces the quantum delocalization of the shared proton (radius of gyration ≈ 0.17 Å) and the O‑O distance fluctuations, matching PIMD while standard MD fails. In bulk liquid water (216 molecules, 300 K) GG‑PI accurately reproduces O‑O, O‑H, and H‑H radial distribution functions and the intramolecular H‑O‑H angle distribution, correcting the over‑structured results of classical MD. For para‑hydrogen modeled as distinguishable spherical particles, the authors demonstrate transfer along a fixed τ‑line: a single model trained on a 64‑molecule system at (T=100 K, P=8) is applied without retraining to a 172‑molecule system at (T=25 K, P=32). Kinetic, potential, and radius‑of‑gyration observables follow the PIMD reference across temperatures.

Performance is quantified by effective sample size per second (ESS/sec). GG‑PI achieves speed‑ups of ~50× for the Zundel ion, ~8.9× for liquid water, and ~1.6× for para‑hydrogen relative to PIMD. The overhead of data collection and model training is a one‑time cost that becomes negligible compared with the force‑evaluation bottleneck of PIMD, especially when expensive ab‑initio potentials are used. The authors also discuss extensions: the local imaginary‑time propagator is independent of boundary conditions, allowing straightforward generalization to open‑chain Path‑Integral Ground State (PIGS) simulations; incorporating exchange effects would enable treatment of indistinguishable fermions or bosons; and the Gibbs‑sampling framework can be applied to transition‑path sampling or diffusion‑bridge generation in classical contexts.

In summary, GG‑PI offers a conceptually simple yet powerful route to quantum statistics: learn the single‑bead conditional density from cheap classical data, then reconstruct the full quantum ensemble via Gibbs sampling. It eliminates the need for repeated force evaluations, transfers across temperatures without retraining, and delivers substantial computational savings while preserving the accuracy of full‑path‑integral methods. This approach opens new avenues for efficient quantum‑aware simulations in chemistry, materials science, and beyond.