Accurate nuclear quantum statistics on machine-learned classical effective potentials

Molecular dynamics (MD) simulations provide insights into phenomena at atomic resolution and are an invaluable tool in chemistry, biophysics, and material science. Recent advances in force field,^1–3 algorithms,⁴ and hardware development⁵ allow for increasingly higher accuracy and a broader scope of application of molecular dynamics simulations.

Traditional MD simulations rely on the classical evolution of nuclear positions, neglecting an explicit treatment of so-called nuclear quantum effects (NQEs) such as zero-point fluctuations, tunneling phenomena, and isotope effects. This simplification has a historical justification as the parameters of classical empirical force field simulations are fit to reproduce (quantum) experimental data. Hence, NQEs are effectively included in the classical potential energy surface, at least partially, albeit with limited transferability. Meanwhile, bottom-up ab initio MD methods or machine learning interatomic potentials (MLIPs), based on the Born–Oppenheimer approximation,⁶ only account for the quantization of electrons. Hence, these simulations require an explicit treatment of quantum nuclear motion on the Born–Oppenheimer potential energy surface (PES) in theory. Particularly when accurate electronic structure methods are used to estimate the Born–Oppenheimer PES, the neglect of quantum nuclear motion may be the largest source of error for obtaining quantitative agreement with experiments.^7–10 Unfortunately, NQEs are typically neglected due to their significant computational overhead or because they are assumed to have small effects. However, even at ambient temperatures, NQEs significantly affect the properties of hydrogen-bounded systems, e.g., static and dynamic properties of water,^11–23 free energies of hydrogen bonding and proton transfer,^24–29 stability of molecular crystals,^30,31 strength of intermolecular interactions,^9,32 and enzyme activity.^19,33

The standard condensed-phase simulation technique for incorporating NQEs in equilibrium properties is Feynman’s imaginary time path integral (PI) formalism.³⁴ In this framework, a system’s quantum statistical mechanics is mapped onto the classical statistical mechanics of a ring polymer, which consists of multiple replicas of the system connected by harmonic interactions that depend on temperature and mass. The effective classical PI partition function of the ring polymer system can be sampled using MD or Monte Carlo techniques, yielding the PIMD³⁵ and PI Monte Carlo³⁶ methods. In theory, PI methods incorporate quantum statistics exactly with classical sampling, in the limit of the number of replicas going to infinity.³⁷ However, typically, the required number of replicas to converge structural properties, and thus the increased computational cost compared to classical MD, rises with increasing physical frequencies or lowering temperatures.³⁸ As a result, most molecular systems require 16–32 replicas at room temperature, which increases proportionally with lowering temperature, such as over 10⁴ replicas to obtain converged energies for molecular systems at cryogenic temperatures.³⁹ 

Over the past decades, several systematically improvable reduced-cost PI approaches have been proposed. These include high-order factorization of the Boltzmann operator,^40–45 perturbative expansion of the Boltzmann operator,^46,47 multiple time stepping,⁴⁸ ring polymer contraction,⁴⁹ and colored-noise thermostats.^50–52 Despite their notable success toward mainstreaming the use of PI methods,³⁸ these approaches do not yet present a “fix-all” solution. Issues still exist, including, but not limited to, lack of universal applicability of multiple time stepping and ring polymer contraction, tedious implementation of high-order path integral techniques, and poor ergodicity⁵³ or zero-point energy leakage⁵⁰ due to the coupling between physical and fictitious PI collective modes. Finally, these approaches do not eliminate the diverging cost of path integral simulations at low temperatures.

An appealing direction that avoids the computational overhead of PI methods is to avoid modeling the entire imaginary-time path. For instance, the Wigner–Kirkwood effective potential^54,55 expands the partition functions in powers of the reduced Planck’s constant and retains computationally tractable terms. Similarly, an effective quantum potential was suggested by Feynman and Hibbs^34,56 by integrating out all degrees of freedom except the centroid of the ring polymer, i.e., by coarse-graining the ring polymer to the centroid. Ever since, the idea of ring polymer coarse-graining has been explored for over 30 years^57–65 and has recently gained increasingly more attention thanks to the advent of MLIPs.^66–68 

A number of coarse-grained effective potentials can be developed by identifying (a) a suitable mapping from the fine-grained ring polymer to the coarse-grained classical system and (b) the functional form of the interactions in the coarse-grained system. For instance, the so-called CG-PI theory, developed by Voth and co-workers,^69–71 proposes to map the imaginary time path to a two-replica coarse-grained system. This mapping gives access to both diagonal and nondiagonal elements of the thermal density matrix. However, the corresponding formalism is challenging to implement for realistic molecular systems in the regime of strong NQEs.⁷¹ Another common mapping involves coarse-graining of the ring polymer to the centroid.^65–68 Among others, we have recently shown that⁶⁶ the resulting effective potential can be used in the same manner as a classical force field. It is well known that coarse-graining the ring polymer to its centroid can still give rigorous access to equilibrium properties, in terms of the so-called centroid-constrained imaginary time correlation function.⁵⁷ However, estimating equilibrium properties, such as the average of a non-linear position-dependent operator, in this manner is usually complicated for generic systems because ensemble averages of the operators evaluated on the centroid coordinates need to be integrated in the high-dimensional configurational space with the centroid-constrained imaginary time correlation function.⁵⁷ On the other hand, the centroid-based description is most effective for an approximate treatment of dynamical properties of complex systems in terms of time correlation functions estimated on the centroid trajectory.⁵⁸ A prominent feature of our work is the path-integral coarse-grained simulation (PIGS) method, which uses many-body universal approximators⁷² to fit the coarse-grained potential of mean force. The PIGS approach utilizes machine-learned potentials, capable of approximating arbitrary complex functions, to learn the full many-body coarse-grained potentials as demonstrated recently on biomolecular simulations.^73–76 

Here, instead of dynamical properties, we focus on characterizing static position-dependent properties, including NQEs, and we develop an effective machine-learned potential to compute them accurately and at the cost of classical MD. Our approach builds upon the PIGS approach that uses multiscale coarse-graining⁷² and our work on the use of deep learning architectures for classical coarse-grained force-fields⁷⁷ but maps the ring polymer onto a single replica instead of the centroid. This coarse-grained mapping is distinct from earlier choices and allows quantum thermodynamic expectation values of position-dependent observables to be estimated rigorously as simple time averages akin to classical MD. We demonstrate that this approach yields accurate quantum nuclear statistics of bulk and isolated systems. The remaining sections of this paper are organized as follows: In Sec. II, we briefly summarize the imaginary-time PI theory and the application of thermodynamic coarse-graining to the problem. Section III discusses the numeric implementation of the approach for four model systems: one-dimensional Morse potential, single H₂O molecule, Zundel cation, and bulk water. The results are presented and discussed in Sec. IV. Section V provides concluding remarks.

Our goal is to reduce the computational cost of simulating ring polymer statistics by integrating additional degrees of freedom, such as (linear combinations of) the replicas of the system, in a way that retains the accuracy of Eq. (2). In the most extreme and computationally beneficial scenario, we can integrate the P − 1 replicas of the physical system to an effective potential felt by a single replica, reducing the dimensionality of the ring polymer to the physical system. In doing so, we reduce the cost of estimating thermodynamic properties to that of classical MD.

The choice of the mapping operator in Eq. (7) is distinct from earlier choices made in the context of coarse-graining of imaginary-time PIs. The most popular coarse-graining approach maps the ring polymer positions to its centroid. While this approach is useful for simulating dynamical properties, it does not give correct equilibrium averages of generic position-dependent operators when estimated as an ensemble average similar to Eq. (10). Similarly, the PICG theory introduced in Refs. 69–71 coarse-grains the ring polymer to two replicas, which are related to the positions at which the off-diagonal elements of the Boltzmann operator are estimated. This approach gives access to generic position- and momentum-dependent operators; however, it requires derivations of bespoke (complicated) operators. Meanwhile, our mapping reduces the ring polymer to the position at which the trace of the Boltzmann operator is estimated, with $Ueff(q)=−β−1⁡logqe−βĤq$⁠, modulo an additive constant, giving access to position-dependent operators as simple ensemble averages.

To develop the effective potentials that incorporate NQEs, we utilize the PIGS approach, outlined in Fig. 1. First, we generate training data via the PIMD simulation. The resulting ring polymer trajectories are then projected onto a coarse-grained space using a configurational and force map discussed above. Notice that the remaining replica j′ in Eq. (7) is selected arbitrarily, and for a ring polymer with P replicas, Eq. (7) defines P equivalent projections. Thus, each ring polymer configuration yields P data points for training.

The training data were generated from PIMD simulations using the i-PI code⁸² and a custom Python-based force provider implementing the one-dimensional Morse potential. Three simulations were performed in an NVT ensemble for 100, 300, and 600 K, using 64 replicas and a PILE-L thermostat with a time constant of 100 fs.⁸³ The integration was performed with a time step of 0.5 fs, and positions and forces acting on each bead were sampled every 20 steps. The first 2.5 ps of the simulation was discarded, and the positions and forces sampled in the subsequent 500 ps were used to train the model.

We use Eq. (11) to model the effective potential of the single water molecule, the Zundel cation, and bulk water. We represent U^ML with the MACE graph neural network architecture.⁸⁴ One of the hallmarks of the MACE architecture is the use of many-body features to represent atomic local environments, improving the model accuracy and data efficiency.

The models were trained with the AdamW optimizer, a learning rate of 0.001, and a weight decay of 0.001. The MACE neural network cutoff was set to 6.0 Å, with 15 radial basis functions, 6 polynomial, and 3 radial channels. The correlation order for Zundel cation and bulk water was set to 3, which corresponds to four-body features. For a single water molecule, three-body features were used. In all cases, 80% of the data available were used for training and the remaining 20% were used for validation. The training was performed for 60 epochs for the H₂O molecule and the Zundel cation and for 15 epochs for liquid H₂O. For each system, ten different models were trained, each using different parameter initialization and train–validation splitting.

The training data were generated from PIMD simulations using the i-PI code using the Partridge and Schwenke potential for a single water molecule,⁸⁵ the HBB potential⁸⁶ for the Zundel cation, and the MB-pol for bulk water.⁸⁷ The simulation parameters are summarized in Table S1. For each system, all the data generated during the production phase of PIMD simulation were used as a training dataset. Since all the replicas are equivalent, each frame of the trajectory yields P data points for model training. We use a force map optimized to minimize the average magnitude of the mapped forces, subject to the constraints of the valid force map. This force optimization was performed according to the procedure introduced in Ref. 77, with an L₂ regularization parameter set to 0.05.

For each of the trained models, we performed molecular dynamics simulations, using as a force field a combination of classical potential and a trained correction as defined in Eq. (11). The simulations were performed in the NVT ensemble, with a Langevin integrator and a friction thermostat with a time constant set to 100 ps⁻¹. The time step and periodic boundary conditions were the same as in the training simulations. All the frames were used for further analysis. For each of the systems, we performed additional simulations with only the classical potential as a control to compare the classical and quantum statistics.

First, we study a particle in a one-dimensional Morse potential, a model for a single O–H bond. The Schrödinger equation for this system can be solved analytically,⁸⁰ thus making Morse potential an ideal test system. It also exhibits significant nuclear quantum effects, as shown in Fig. 2(a), where the classical Morse potential [Eq. (17)] is compared with the corresponding quantum PMFs [Eq. (18)] at different temperatures. The classical approximation leads to the over-localization of the particle, and, as expected, the quantum effects are more prominent at low temperatures. However, for this system, even at a temperature of 600 K, there are significant deviations between the classical and quantum pictures.

The developed approach allows us to get a quantitative agreement with the reference, especially in the high-probability regions that are well-sampled in the training dataset [Figs. 2(b)–2(d)]. Further away from the minimum energy, the model’s prediction slightly deviates from the ground truth. This is to be expected, as these regions of the phase space are not as well represented in the training data as the high-probability regions.

We test the reweighting scheme proposed in Sec. II E by generating datasets corresponding to different temperatures (300 K) and mass of the particle that is either protium (corresponding to an O–H bond) or increased to tritium (corresponding to an O–T bond) from the 600 K O–H dataset previously used. The resulting reweighted datasets were used to train new models. As shown in Fig. 2(e), the models trained with the reweighted dataset are in excellent agreement with the analytical results for the corresponding temperature and isotope composition. Such data recycling allows for a further decrease in the computational cost of generating training data, although it becomes less effective for more complex systems. Nevertheless, we expect that this approach may be particularly useful for generating large datasets for training transferable models.

The next test case we considered is a single water molecule in vacuum. The quantum statistics of this three-body system can be exhaustively described in terms of the length of two O–H bonds and the H–O–H bond angle. The comparison of the probability distributions for the O–H distances and H–O–H bond angles is shown in Fig. 3, with additional details provided in Fig. S4. As was the case for the Morse potential, the classical distribution is significantly more localized than the quantum counterpart. The ML model consistently samples the correct quantum distribution: the results obtained from ten independently trained models produce similar results with low variance.

For all trained models, we run the MD simulation for 1 ns, which is four times longer than the PIMD simulation used to generate the training data. This allows us to obtain slowly converging properties by performing short PIMD simulations and using those to train a model that can be run for a much longer time.

We further investigated the impact of PIMD simulation length on model performance by training additional models using 2.5, 12.5, 25, and 125 ps of post-equilibration PIMD data, compared to the primary model trained with 250 ps. All models generated stable 1 ns trajectories. Notably, even with just 5% of the original dataset, the models showed on average good agreement with the reference PIMD data [Figs. S1 and S2]. As the dataset size increased, the Jensen–Shannon divergence decreased and eventually plateaued by 250 ps (Fig. S3).

Even at 600 K, the water molecule predominantly populates the ground vibrational state, with a negligible population of the excited states. Thus, as discussed in Sec. II E, the model trained at 600 K can be used to simulate the system at any lower temperatures. To demonstrate this point further, we report in Fig. S5 the results obtained with a model trained on PIMD data generated at 600 K and then simulated at lower temperatures, with forces rescaled according to Eq. (13).

A more challenging case is the Zundel cation, which consists of two water molecules sharing a hydrogen-bonded proton H*. This system has been extensively studied both theoretically^88,89 and experimentally.^90–94 The Zundel cation exhibits large conformational fluctuations, including anharmonic and floppy motion associated with the transfer of the hydrogen-bounded proton. Following Ref. 89, we characterize these fluctuations by means of the O–O distance d(O–O), the proton transfer coordinate δ = d(O₁H*) − d(O₂H*), the dihedral angle ϕ that measures the rotation of the H₃O⁺ fragments around the O–H*–O axis, and the d(O–plane) that measures the non-planarity of the H₃O⁺ fragment. The dihedral angle ϕ is defined by points X₁–O₁–O₂–X₂, where X₁ and X₂ represent a midpoint between the hydrogen atoms that are not involved in the hydrogen bond formation and bound to O₁ and O₂, correspondingly. The non-planarity of H₃O⁺ is characterized by a distance between an O atom and the plane, defined by the three closest H atoms (see the lower panel of Fig. 4).

At 300 K, the proton transfer coordinate and the rotational motion are almost classical, and only a small correction is needed to account for NQEs. The O–O distance and the O−plane distance require larger corrections, yet the NQEs are relatively small. At 100 K, however, significant quantum effects appear in all the collective variables, including changes in the equilibrium geometry. Nevertheless, our methods recover the accurate quantum distribution both at 300 K and at 100 K. Despite the diversity of conformational fluctuations, the model not only reproduces the target probability distributions but also leads to remarkably stable simulations: all the models were used to perform 1 ns simulation and did not exhibit any instabilities, often plaguing ML force fields.^95,96

Finally, we applied our approach to model liquid water to reproduce the NQEs manifesting in both bonded and nonbonded interactions. Since NQEs are highly local, the former are more affected than the latter. The quality of the machine-learned model is evaluated in Fig. 5 by computing radial distribution functions (RDFs). At 300 K, the first peak in the O–H [Fig. 5(b)] and H–H [Fig. 5(c)] RDFs predicted by classical MD (P = 1) are too narrow and significantly deviate from the correct quantum PIMD results. The deviations between the classical and quantum RDFs rapidly decrease with the distance.

Our neural-network-based potential is able to accurately predict the correct quantum RDFs, including the contributions that arise from both bonded and nonbonded interactions (Fig. 5), in contrast to another similar study.⁶⁸ There is also a remarkable agreement between independently trained models. The high accuracy and low uncertainty of our ML model prediction show that it reliably captures the subtle differences between classical and quantum RDFs.

Compared to the simulations for the ML models of a single water molecule and of the Zundel cation, the simulations for bulk water exhibit a decreased stability. We note that some simulations of our ML models, trained on 35 ps of PIMD simulations, exhibit numerical instabilities. However, the ML models allow, on average, at least 38 ps of simulation time before the onset of any instability. Simulations can be reinitialized and restarted to improve statistics. The numerical instability of the model can be potentially cured by increasing the size of the training dataset. Alternatively, one can modify the learning problem by simplifying the ML correction, for example, by introducing physics-informed energy terms, as it is done in the learning of classical ML coarse-grained force-fields.^73–76 Given that the stability of the simulation strongly depends on the parameter α [see Eq. (11)], this approach may be especially promising.

We developed an approach to coarse-grain the ring-polymer dynamics describing NQEs, and we obtain results thermodynamically consistent with PIMD simulations. This approach allows us to accurately compute position-dependent NQEs at the cost of classical MD. We have demonstrated this method on four distinct systems where NQEs are prevalent—the Morse potential, a single water molecule, the Zundel cation, and bulk water—and reported results in excellent agreement with (much more expensive) reference PIMD simulations.

The proposed method builds an effective NQE correction potential by combining the rigorous framework of multiscale coarse-graining⁷² with the MACE architecture,⁸⁴ an expressive state-of-the-art graph neural network potential. We have shown that these effective potentials accurately reproduce the quantum statistics of nuclei, both near the classical limit and in the regime of strong NQE. Moreover, our method eliminates the need to simulate multiple replicas. Once parameterized, the model seamlessly integrates into existing classical force fields and does not require specialized algorithms other than ones already used for classical MD simulations. However, this approach is limited to reproducing quantum thermodynamics, so dynamical and momentum-dependent properties should be approximated using a previously proposed method.⁶⁶ 

For the results presented here, different sets of training data were used to train models for different molecular species at different temperatures. However, if the temperature associated with the training data is low enough to significantly populate only the quantum ground state, the trained model can be used to also simulate any lower temperature, as we have demonstrated in the case of the H₂O molecule. The proposed method is designed such that the forces originating from the coupling between the adjacent replicas have a non-zero contribution to the force projected onto the coarse-grained space. As the dependence of these forces on the atomic mass and temperature is known [Eq. (9)], the training dataset obtained for one isotope composition and temperature can be, in principle, reweighted to construct a dataset for different isotope compositions and temperatures. Thus, a single PIMD simulation can provide data for training several models, as we have shown in the case of the Morse potential.

The rapid proliferation of machine-learned force fields parameterized using accurate quantum chemical calculations to account for the electronic effects means that the missing NQEs are now becoming a dominant source of error. At the same time, NQEs are crucial to describe biomolecular processes accurately, for instance, in the case of proton transfer reactions. While the results presented here are model-specific, building upon the recent development of machine-learned force-fields (either coarse-grained^73–76 or atomistic^{81,84,97–100}), we believe that this study represents a significant step forward toward the design of transferable potentials that explicitly incorporate also NQEs but retain the computational cost of classical MD.

Section I of supplementary material details the Morse potential (Sec. I A) and ML model used for a particle in Morse potential (Sec. I B). Section II contains the simulation parameters used to generate the training data. Section III provides the derivation of coordinate rescaling. Section IV presents extended results for a water molecule.

We thank Gregory Voth and members of the Clementi’s group at FU for insightful discussions and comments on the manuscript. We acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) (Grant Nos. SFB/TRR 186, Project A12; SFB 1114, Projects B03, B08, and A04; and SFB 1078, Project C7), the National Science Foundation (Grant No. PHY-2019745), and the Einstein Foundation Berlin (Project No. 0420815101), and we also acknowledge the computing time provided on the supercomputer Lise at NHR@ZIB as part of the NHR infrastructure (Project No. beb00040), the HPC Service of FUB-IT, Freie Universität Berlin, the HPC Service of the Physics department, Freie Universität Berlin, and the Swiss National Supercomputing Centre (under Project Nos. s1209 and s1288). V.K. acknowledges the support from the Ernest Oppenheimer Early Career Fellowship and the Sydney Harvey Junior Research Fellowship, Churchill College, University of Cambridge.

The authors have no conflicts to disclose.

Iryna Zaporozhets: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Félix Musil: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Methodology (equal); Software (equal); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Venkat Kapil: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Cecilia Clementi: Conceptualization (lead); Funding acquisition (lead); Investigation (supporting); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal).

All the training data, code, and scripts that support the findings of this study are openly available in the online repository: cg_nuclear_quantum_statistics.