According to Feynman, the quantum partition function of a particle can be expressed in a path integral form, which, after discretization, is analogous to the classical partition function of a ring-polymer. For example, in the one-dimensional case (Tuckerman, 2010): Q = Tr [ e − β H ^ ] = ∮ D x ( τ ) × exp { − 1 ℏ ∫ 0 β ℏ d τ ⋅ [ m 2 x ˙ 2 ( τ ) + U ( x ( τ ) ) ] } ≈ ∫ d x 1 … d x P × exp { − 1 ℏ ∑ k = 1 P [ m P 2 β ℏ ( x k + 1 − x k ) 2 + β ℏ P U ( x k ) ] } | x P + 1 = x P ≈ Q P (1)

In here, Q P stands for the geometric part of the partition function of a classical cyclic ring polymer with P beads. U ( x ) is the PES employed in the simulation, computed using either force field models or ab initio approaches. Exact quantum Boltzmann distribution can be obtained by taking P → ∞ . In practice, we use limited number (typically P equals to 32 or 64 in room temperature) of slices to represent the continuous path of x ( τ ) . The conjugate momentum of the bead positions x k can be introduced in a variety of ways, leading to different flavors of path-integral-based methods such as RPMD, PA-CMD (partially adiabatic-CMD) (Hone et al., 2006), and CMD etc. In particular, CMD coarse grains the intra-ring vibrational modes and focuses on the physical motion of the geometric center of the ring-polymer (i.e., the “centroid”): x c = 1 ℏ β ∫ 0 ℏ β d τ   x ( τ ) ≈ 1 P ∑ k = 1 P x k (2)

When the motion of the centroid and the intra-ring vibrations are decoupled, the centroid can be considered as moving on an adiabatic PMF (denoted as U c ( x ) ), with the effective quantum centroid force defined as: F c ( x 0 ) = − d U c ( x 0 ) d x 0 = − ∮ D x ( τ ) δ ( x c − x 0 ) exp { − S / ℏ } [ − d U ( x ( τ ) ) d x ] ∮ D x ( τ ) δ ( x c − x 0 ) exp { − S / ℏ } ≈ 〈 F ( x k ) 〉 | x c = x 0 (3)

Here, the action S is simply: S = ∫ 0 β ℏ d τ ⋅ [ m 2 x ˙ 2 ( τ ) + U ( x ( τ ) ) ] (4)

In principle, one can conduct either MD or MC to sample the bead ( x k ) distributions while keeping the centroid positions ( x c ) fixed, then find the adiabatic centroid force by averaging the forces on the beads. But such scheme is computationally inefficient. Therefore, a common practice is to propagate all the motions simultaneously, but to set the masses of the intra-ring vibrations to be lighter, so these motions are decoupled from the physical centroid motion. However, such scheme requires a much smaller time step compared to regular classical MD simulations, thus causing extra computational cost.

Another simple alternative was proposed by Feynman and Hibbs (Feynman, 1998), when they proved that in a variational sense, the optimal effective classical potential K ( x ) that accounts for NQE is a Gaussian convolution of the underlying physical potential U ( x ) : K ( x ) = 6 m k T π ℏ 2 ∫ d x ′ ⋅ U ( x ′ ) exp [ − 6 m k T ( x − x ′ ) 2 ℏ 2 ] (5)

It is interesting to observe that the Gaussian kernel features a width of σ = ℏ 2 / 12 m k T , which happens to be the width of the ring-polymer bead distribution in the free particle limit (i.e., U ( x ) = 0 ). Therefore, FHC can be viewed as an approximation to CMD when the shape of U does not affect the shape of the ring-polymer significantly.

When the underlying potential is pairwise additive and spherically symmetric (i.e., U = ∑ i < j u ( r i j ) ), Eq. 5 leads to a quite simple correction to the original pairwise interactions. This correction, truncated at second order, can be written as: u F H ( r i j ) = u ( r i j ) + β ℏ 2 24 μ [ ∂ 2 u ( r i j ) ∂ r i j 2 + 2 r i j ∂ u ( r i j ) ∂ r i j ] (6)

Here, μ = 1 1 / m i + 1 / m j represents the reduced mass of the two interacting atoms. FHC provides an extremely simple solution to the NQE problem, avoiding the expensive explicit sampling of the ring polymer geometry. Therefore, FHC is widely used nowadays in complicated simulations such as the GCMC, which are often performed to compute the adsorption of gas molecules in porous materials. However, for a general many-body potential (e.g., ab initio PES), FHC does not give a simply analytical formula as given in Eq. 6, and cannot be easily computed without explicit MC sampling. Therefore, even though both CMD and FHC offer an effective quantum correction to the physical potential formally, neither of them can be conducted easily on a general many-body PES.

In this work, we will use the embedded atom neural network (EANN) technique to learn the quantum centroid force correction (Cao and Voth, 1994a; Feynman et al., 2010): Δ F c ( x 0 ) = F c ( x 0 ) − F ( x 0 ) (7)

EANN, like many other popular ML methods, assumes that the total energy (in this work, the centroid PMF) can be written as a sum of atomic energies: U = ∑ a u a ( f → a ) (8) And each atomic energy u a is a function of its local environment, the information of which is encoded in an atomic feature vector f → a . Neural networks (NN) are then employed to fit the complicated many-body function u a ( f → a ) . Different ML methods differ in the ways of constructing f → a , but the basic mathematical structure showed in Eq. 8 is always retained. One limitation of this scheme is that f → a can only encode the local atomic environment within a certain range r c , and the computational cost of the ML model increases rapidly with increasing r c . Many deficiencies of the ML methods manifested in the molecular force field development problem root in this limitation. Therefore, the ML technique is in natural favor of extremely localized problems, which require smaller distance cutoffs, thus resulting in more efficient training processes and faster computational speeds.

In this work, we will show that in molecular systems, the quantum centroid force correction ( Δ F c ) happens to be an extremely localized target, so ML is a natural tool to learn this target. At first glance, it may not be very obvious why the machine learned quantum force correction is more efficient than the plain CMD. The training data collection could be painful: we need to fix the centroids at different geometries, and for each geometry, we have to run a thorough PIMD or PIMC simulation to fully converge the centroid force. This fully adiabatic procedure has been shown to be less efficient than the state-of-the-art CMD algorithm, which runs the ring-polymer dynamics on-the-fly using extended Lagrangian methods. However, in molecular systems, if the quantum force correction is localized enough, we can train the correction on small molecule clusters and apply it to bulk simulations. Such an approach brings a huge gain on computational efficiency, especially for potentials that are expensive to evaluate. Multipolar polarizable force fields (such as AMOEBA) and ab initio potentials are two typical examples: In the former case, bulk calculations involve the time-consuming multipolar Ewald summation, while cluster calculations with a simple cutoff scheme are much faster to run. In the latter case, the computational cost of regular DFT calculations scales as O ( N 3 ) , thus the bulk PIMD is also much harder to perform compared to small cluster sampling. Therefore, with the development of the new generation of accurate many-body force fields, it is becoming increasingly beneficial to train the NQE correction on small molecular clusters while using it in bulk simulations.

In this work, we use water as our test example, with both q-TIP4P/F (Habershon et al., 2009) and AMOEBA (Ren and Ponder, 2003) force fields. The two force fields were selected due to the following reasons: 1) Both of them are not too expensive to run, so we can obtain the rigorous bulk PIMD benchmark data at a reasonable cost; 2) The two force fields represent two distinct categories of PES. The q-TIP4P/F force field is one of the conventional force fields with simple point charges and Lennard-Jones interactions, which are completely spherical symmetric and pairwise-additive. Meanwhile, AMOEBA features multipole moments and explicitly polarizable atoms, thus it resembles a more general many-body PES. Therefore, it is interesting to compare the performances of the localized ML model on both cases, in order to examine the generality of the method.

For each force field, we first performed a bulk (1,000 water) NPT PIMD simulation at the corresponding condition (1 bar, 300 K or 100 K), and randomly drew 2000 cluster configurations from the bulk trajectory. Each cluster was consisted by one central water molecule and seven nearest molecules, so we guarantee that the first solvation shell of the central water was fully included. These water octamers were then used to generate the training data: with centroid positions fixed, we ran 1 ns PIMD simulation to sample the bead distribution, so the averaged centroid force is converged within 0.2 kJ/(mol Å). To examine the width of the ring-polymer, we also performed 1 ns PIMD simulations on a smaller bulk system (with 216 water) with a density of 0.997 g/ml at different temperatures. All the sampling simulations mentioned above were conducted using the OpenMM 7.4.0 program, with modifications made to enable the fixed-centroid PIMD simulations. All sampling simulations were performed using 0.5 fs timestep and a Langevin thermostat with a 1.0 ps⁻¹ friction constant. We used 32 beads at 300 K and 64 beads at 100 K, and the cutoff distance for the nonbonded interaction was set to 9 Å.

Once the averaged centroid forces were collected, their differences with the classical forces (i.e., Δ F c ) were used to train an EANN model. The maximum angular momentum in the EANN model was set to 2, and 11 Gaussian Type Orbitals (GTOs) were used in the radial dimension. We used two hidden layers in the network, with 20 neurons in each layer. All the parameters in EANN were trained using the hybrid extreme machine learning and Levenberg-Marquardt (ELM-LM) algorithm, with the convergence criterion for the loss function set to be 0.1.

The trained EANN model was then used in combination with the original force field to perform the machine learning CMD (ML-CMD) simulations. For validation, the simulation results were compared to the PIMD, RPMD, and the PA-CMD results, and all the production runs were conducted using the i-PI 2.2 program. In PA-CMD simulations, the lower limit of internal modes is set to 648 THz to ensure the adiabaticity between the centroid mode and the intra-ring modes, and the timestep was set to be 0.02 fs.

As we discussed in the theory section, the locality of the Δ F c is the key to the success of the ML models, and thus the key to the success of this work. Note that the quantum force correction is due to the smearing of the atom into a cyclic ring-polymer, so the correction would vanish if the size of the ring-polymer is small enough compared to the interaction range. Therefore, it is both interesting and necessary to rigorously quantify the distribution width of the ring-polymer at different conditions. Starting from Eq. 1, without losing generality, we assume that the centroid is located at x c = 0 , then expand and truncate the potential energy U around x c at the second order: Q ≈ ∫ D x ( τ ) exp { − 1 ℏ ∫ 0 β ℏ d τ ⋅ [ m 2 x ˙ 2 + 1 2 U ″ x 2 ] } (9)

In here, U ″ = ∂ 2 U ∂ x 2 is the force constant of the potential at point x c along dimension x . Fourier transforming x (the n = 0 term is dropped since we assume that x c = 0 ): { x ( τ ) = ∑ n = − ∞ n ≠ 0 ∞ e i ω n τ x n ω n = 2 π n β ℏ   ( n = … , − 3 , − 2 , − 1,1,2,3 , … ) (10) and we can rewrite Eq. 9 as: Q = ∫ D x ( τ ) exp { − β 2 ∑ n [ m ω n 2 + U ″ ] | x n | 2 } (11)

Therefore, the distribution of x n is merely independent Gaussian distributions with the second moments: 〈 x n 2 〉 = 1 β [ m ω n 2 + U ″ ] = 1 m β [ ω n 2 + ω 2 ] (12) And in here, the intrinsic vibrational frequency of the atom is: ω = U ″ m .

Therefore, the total second moment of the distribution of x can also be written as a sum of all Matsubara frequencies: 〈 x 2 ( τ ) 〉 = ∑ n = − ∞ n ≠ 0 ∞ 〈 x n 2 〉 = 1 m β ∑ n = − ∞ n ≠ 0 ∞ 1 ω n 2 + ω 2 (13)

Define the dimensionless quantity t as: t = β ℏ ω 2 π (14)

Then the second moment of the ring-polymer bead distribution is: 〈 x 2 〉 = β ℏ 2 4 π 2 m ∑ n = − ∞ n ≠ 0 ∞ 1 n 2 + t 2 = β ℏ 2 2 π 2 m [ − 1 2 t 2 + π 2 t coth ( π t ) ] (15)

Apparently, the locality of the ring-polymer is determined by both the temperature and the curvature of the underlying potential. In the high temperature limit ( β → 0 ,   t → 0 ), Eq. 15 approaches the limit of ℏ 2 12 m k T , which is exactly the width used in the convolution kernel of FHC. When the temperature is high, the quantum effect is weak, so the ring-polymer is more similar to a classical particle with a very narrow width. In this case, the potential is smooth enough in the relevant neighborhood of the centroid, with the effects of its curvature being negligible. Therefore, the shape of the ring-polymer can be approximated well assuming the free-particle limit, and FHC can be viewed as an excellent approximation to CMD in the high-temperature region. However, in low temperature region, FHC fails miserably as it predicts a singular behavior with the distribution width going to infinitely large. Meanwhile, Eq. 15 in the low temperature limit ( β → ∞ ,   t → ∞ ) gives a asymptotic value of 〈 x 2 〉 = ℏ 2 m ω , which is exactly the distribution width of the ground state of a quantum harmonic oscillator (QHO). This is also not surprising, since the particle should only populate the vibrational ground state in the low temperature limit. It can be further proved that both the QHO and the FHC widths are rigorously the upper bounds of the actual distribution width. The most delocalized (thus the most problematic) dimension is the one that involves the lowest frequency mode of the lightest atom.

In Figure 1, using liquid water as example, we plot the widths of the H atoms given by the square root of Eq. 15, in conjunction with the values given by the FHC and the QHO formula, as well as the PIMD simulation results. The intrinsic frequency ω is set to be 670 cm⁻¹, which is approximately the frequency of the librational mode (i.e., the lowest frequency mode) of H in liquid water and ice (Tong et al., 2016). As it is shown, the analytical results agree with the simulation results in the entire temperature range, validating the utilization of Eq. 15 in the estimation of the bead distribution width. Numerically, bounded by both the QHO and the FHC limits, the average width of the H atoms is below 0.15 Å at all temperatures, showing that the PIMD bead distribution in a typical condensed molecular system is indeed extremely localized.

Using the results of the bead distribution width, we can roughly estimate the magnitude of the quantum force correction at a certain range. Let us consider a pair of particles interacting with Coulombic force (the force with the slowest decay rate in molecular systems), and evaluate how the quantum correction to the magnitude of the force decays with increasing distance. For convenience, without writing down the charges explicitly, we assume that the magnitude of the classical force between the two particles is: F ∝ 1 R 2 (16) And let us assume the three-dimensional ring-polymer bead shifts for particle i and particle j are x → i and x → j , respectively. Both x → i and x → j are with respect the centroids, thus we have: 〈 x → i 〉 = 〈 x → j 〉 = 0 .

Then, the averaged quantum force correction, truncated at second order, is: Δ F = 〈 1 | R → + x → i − x → j | 2 〉 − 1 R 2 = 〈 1 | R → + Δ x → | 2 〉 − 1 R 2 ≈ 〈 Δ x → ⋅ [ 4 ( R → R → ) − R 2 I R 6 ] ⋅ Δ x → 〉 = 1 R 4 ( 4 〈 ( n → ⋅ Δ x → ) ( n → ⋅ Δ x → ) 〉 − 〈 Δ x → ⋅ Δ x → 〉 ) (17)

Hence the relative correction is: Δ F F ≈ 1 R 2 ( 4 〈 ( n → ⋅ Δ x → ) ( n → ⋅ Δ x → ) 〉 − 〈 Δ x → ⋅ Δ x → 〉 ) (18)

Considering Δ x = n → ⋅ Δ x → is simply the one-dimensional projection of Δ x → along the direction of R → , and 〈 Δ x → ⋅ Δ x → 〉 ≥ 〈 Δ x 2 〉 , we have: Δ F F ≤ 1 R 2 ( 4 〈 Δ x 2 〉 − 〈 Δ x 2 〉 ) = 3 R 2 [ 〈 x i 2 〉 + 〈 x j 2 〉 − 2 〈 x i x j 〉 ] (19)

The largest force correction happens when the direction of R → (i.e., n → ) coincides with the lowest frequency mode (the librational mode) of hydrogen. Therefore, we can simply use the librational frequency and Eq. 15 to evaluate the upper bounds for 〈 x i 2 〉 and 〈 x j 2 〉 . According to Cauchy-Schwarz inequality, we further have: | 〈 x i x j 〉 | ≤ 〈 x i 2 〉 〈 x j 2 〉 = 〈 x 2 〉 . Therefore, we have: Δ F F ≤ 12 R 2 〈 x 2 〉 (20)

This inequality could be further improved, realizing that the correlation ( 〈 x i x j 〉 ) also decays as 1 / R 3 (see supporting information for more analysis to the position and force correlations). However, Eq. 20 is already good enough for the following numerical analysis.

In this work, we use a cutoff distance of 4 Å for our ML model, and at this distance, Eq. 20 gives a relative force correction as small as 2% at all temperatures. Furthermore, both analytical and numerical analysis indicate that the correlation between the force corrections decays fast with respect to interatomic distance. Therefore, it is strongly indicated that the quantum force correction is an extremely localized phenomenon and can be tackled using a very small localized ML model.

It is noted that previous studies (Rognoni et al., 2021) show the vibrational modes of liquid water is highly nonlocal. While there are certainly nonlocal NQE existing in aqueous systems, such result is not necessarily contradictory to our findings in this work. In here, we focus on the locality of the centroid force correction: just like localized force can lead to nonlocal vibrational modes, local force correction can certainly give rise to global NQE in bulk systems.

Using the water octamer PIMD centroid forces as training data, we train the EANN model with a small cutoff of 4 Å, and test its validity in both bulk and surface water systems. We examine the accuracy of the EANN model on the quantum force corrections for both Q-TIP4P/F and AMOEBA potentials, and the test results are shown in Figure 2. Excellent agreement between the fitted and the reference data are achieved in all cases, showing the great success of the methodology. This result is in consistent with our analysis in the last section, showing that the quantum force correction is virtually fully determined by the structure of the first solvation shell. Water octamer clusters that contain one fully solvated water molecule and seven surface molecules is all it needs to obtain a model that is transferrable to both bulk and surface environments. Such extrapolation capability does not rely on the pairwise additivity of the underlying potential, and is extremely important for PES that is difficult to evaluate (such as AMEOBA).

After examining the centroid force correction, we investigate the performance of the EANN model in the bulk ML-CMD simulations. We first investigate the centroid radial distribution functions (RDFs) for both O-O and O-H atom pairs, and the results are plotted in Figure 3. It is noted that the RDF reported here is not the regular RDF reported in literature, but the RDF computed by centroid positions. Although in principle, there is no essential technical difficulties to train the PMF of one bead (instead of the centroid) such that the real RDF can be obtained. Apparently, compared to pure classical MD, the machine learned quantum force correction reduces the error on the centroid RDFs significantly compared to the exact PIMD simulation. It is further observed that the NQE on the centroid RDFs are more significant in the AMOEBA simulations compared to the q-TIP4P/F simulations, and our correction is also more effective in AMOEBA. For both force fields, the quantum centroid RDF is generally less structured compared to the classical results, and such effect is accurately captured by the EANN model.

Besides the liquid structure, we further examine the NQEs on the system pressure (or, equivalently, the system density), and the results are shown in Figure 4. In this part, we performed the reference PIMD simulations using the more efficient OpenMM program, while the ML-CMD was conducted using i-PI for convenience. To rule out the potential numerical differences caused by different programs, we computed the classical baselines using both codes. The classical pressures computed by both OpenMM and i-PI are in excellent agreement, proving that the different integrators and thermostat implementations in the two codes do not affect our simulation results significantly. With the classical baselines established, it is clearly shown in Figure 4 that ML-CMD reproduces the exact PIMD results with a reasonable accuracy. Small overestimations can be observed in both the low-density and the high-density ends, potentially because that the EANN model was trained at 0.998 g/ml, so the model error is slightly larger when the density deviates from this density. In summary, it is shown that in both the centroid RDF and the pressure tests, the ML force correction trained by only small cluster data can reproduce the bulk PIMD results accurately, proving the capability of our methodology in the studies of thermodynamic properties.

As we discussed in the theoretical background section, the ML-CMD method is theoretically equivalent to adiabatic CMD, so it should also be able to capture the NQE in dynamic properties at CMD level of theory. We thus computed the self-diffusion constant of bulk water using classical MD, PA-CMD, and ML-CMD, the results of which are listed in Table 1. The effective mass used in PA-CMD is small enough such that the PA-CMD simulation we perform is a good approximation to the rigorous adiabatic CMD. Once again, the ML-CMD agrees with PA-CMD within the statistical uncertainty, validating our methodology.

To verify that ML-CMD indeed possesses significant advantages in computational efficiency, we compare the running speeds of classical MD, T-RPMD, PA-CMD, and ML-CMD in Table 2. AMOEBA force field is employed, and all simulations are run using i-PI, with one RTX 2080 graphic card, 8 CPU cores, and 80 GB of memory. The EANN model has not been implemented in CUDA platform, so this part of the calculation is performed on a separate CPU client, while all other calculations are run on GPU. Comparing ML-CMD with classical MD, we can see that even with a slow CPU implementation, the extra computational time due the EANN model is minor. Because of the extremely short cutoff distance (4 Å) we use, the EANN model is rather small and is very fast to evaluate. Meanwhile, both T-RPMD and PA-CMD are orders of magnitude slower than ML-CMD and classical MD, due to two major reasons. First, both simulations have to propagate the motions of 32 replicas of the system, which means at least 32 times slower due to more energy and force evaluations. Second, the timestep has to be smaller in both T-RPMD and CMD to maintain the stability of the simulations. We do acknowledge that a better integrator can be employed so a larger timestep is possible for at least T-RPMD, and techniques such as contraction can help to decrease the number of beads employed in force evaluation. All these techniques can alleviate, but not cure the heavy computational cost of path integral simulations. In comparison, ML-CMD is by far the simplest, and the most straightforward way to achieve the classical-like simulation speed, circumventing the complicated and highly specialized techniques used in PIMD.