A “short blanket” dilemma for a state-of-the-art neural network potential for water: Reproducing experimental properties or the physics of the underlying many-body interactions?

Molecular mechanics (MM) force fields (FFs)^1,2 have been the workhorse in computational chemistry since the early days of Monte Carlo (MC)³ and molecular dynamics (MD) simulations.⁴ Continued progress in hardware technologies,⁵ accompanied by the development of more realistic representations of electrostatic interactions, has enabled not only molecular simulations of progressively larger systems, but also the use of more rigorous polarizable FFs^6–10 that go beyond the pairwise additive approximation adopted by conventional fixed-charge FFs.^11–14 

At the same time, the development of efficient algorithms for correlated electronic structure methods, such as coupled cluster theory,^15–17 has enabled routine calculations of interaction energies for molecular clusters with chemical accuracy.^18–20 This has led to the emergence of a new class of analytical potentials that quantitatively reproduce each individual term of the many-body expansion (MBE) of the energy²¹ calculated using correlated electronic structure methods. When applied to aqueous systems^22–35 and molecular fluids,^36–38 these many-body (MB) potentials exhibit unprecedented accuracy, enabling predictive simulations from the gas to the condensed phases.³⁹ Concurrently, machine learning (ML) approaches have gained popularity in computational molecular sciences mainly due to the rapid evolution of graphics processing unit (GPU) and tensor processing unit (TPU) architectures.⁴⁰ In particular, potentials represented by deep neural networks (DNNs) derived from electronic structure data are routinely used to model various molecular systems, from clusters to liquids and materials.^41–65 

State-of-the-art MB and DNN approaches use regression algorithms to construct data-driven representations of the multidimensional energy landscape of the system of interest. This process involves generating representative training sets of reference data calculated at an appropriate level of theory. While MB potentials require tailored parametric functions for each term of the MBE, DNN potentials are usually trained on the total energy and forces of the entire system. By applying embedding schemes to construct low-dimensional descriptors of molecular environments, DNN potentials can compute the gradients required for the propagation of the equations of motion in MD simulations more efficiently than MB potentials.⁶⁶ On the other hand, since MB potentials only use information about small clusters, the corresponding training data can be calculated at a higher level of theory than DNN potentials. As a matter of fact, MB potentials are usually trained on reference energies calculated at the coupled cluster level of theory, including single, double, and perturbative triple excitations, i.e., CCSD(T), which is currently referred to as the “gold standard” for chemical accuracy.⁶⁷ Furthermore, by construction, the functional form of MB potentials allows for accurately representing all physical contributions to the interaction energies, including both short- and long-range many-body effects.^68–70 

To account for long-range interactions, DNN potentials are often trained on condensed-phase configurations, which allows for modeling long-range effects either implicitly, by effectively encoding long-range contributions as short-range representations, or explicitly, by adding effective electrostatic terms.^{49,59,65,66,71,72} This implies that, due to the large number of molecules required to model condensed-phase systems, a lower level of theory than CCSD(T), usually density functional theory (DFT),⁷³ has to be used to retrieve the reference energies. In this context, it has recently been shown that the interplay between functional-driven and density-driven errors may impact the overall accuracy of DFT models and their transferability from gas-phase to condensed-phase systems.^74–79 By construction, these limitations also affect the ability of DNN potentials derived from DFT reference data to “extrapolate” to thermodynamic state points different from those used in the training process.⁸⁰ 

In this work, we investigate the possibility of integrating the best features of MB potentials (i.e., accuracy and transferability) and DNN potentials (i.e., speed and ease of use) into a computational framework that can enable large-scale MD simulations with chemical accuracy. To this end, we focus on the molecular modeling of water as a prototypical system, which has posed several challenges since the early days of MC and MD simulations^81,82 due to its rich phase behavior, characterized by several anomalies.⁸³ As a representative state-of-the-art MB potential, we selected MB-pol^30–32 due to its demonstrated ability to correctly predict the properties of water across the entire phase diagram,⁸⁴ including gas-phase clusters,^85–87 liquid water,⁸⁸ the vapor–liquid interface,^89–91 and ice.^92–95 MB-pol has also recently been used to predict structural and thermodynamic properties of supercooled water down to 200 K at 1 atm, which were found to be in excellent agreement with experimental data that are available above 225 K.⁹⁶ However, due to the associated computational cost, the MD simulations with MB-pol reported in Ref. 96 were limited in terms of both system size (up to 512 water molecules) and sampling time (up to ∼130 ns). The prospect of developing a fast DNN potential trained on MB-pol simulation data, which retains the same accuracy of MB-pol across the entire phase diagram, is thus particularly appealing. This will enable large-scale simulations of water as a function of temperature and pressure, which will provide further insights into water’s anomalous behavior and allow for full exploration of the so-called water’s “no man’s land” that has been proven difficult to probe experimentally.^97–102 To this end, we selected DeePMD^49,66,71 as a representative state-of-the-art framework for developing a DNN potential of water trained on the MB-pol simulations of Ref. 96. DeePMD-based DNN potentials have already been used in MD simulations of various molecular systems, including water,^103–106 ionic liquids,¹⁰⁷ and metals,^108–110 and enabled MD simulations with up to 10 × 10⁹ atoms.¹¹¹ 

This article is organized as follows: In Sec. II, we summarize the main features of the MB-pol potential (Sec. II A) and the DeePMD framework (Sec. II B). In Sec. III A, we first assess the ability of the DeePMD-based DNN potential to reproduce thermodynamic and structural properties of liquid water, calculated with MB-pol, from the boiling point down to deeply supercooled temperatures. We then use the DNN potential to characterize several vapor–liquid equilibrium properties, as well as many-body dependent properties of gas-phase clusters. In Sec. III B, we introduce two other DeePMD-based potentials—DNN(VLE10) and DNN(VLE20)—that are trained on an expanded training set that adds vapor–liquid configurations to the training set used to develop the DNN potential. The performance of both DNN(VLE10) and DNN(VLE20) potentials is assessed with the same structural and many-body-dependent properties used to assess the performance of the DNN potential. In Sec. III C, we introduce three DeePMD-based potentials— DNN(MB4), DNN(MB10), and DNN(MB20)—which are trained to incorporate low-order many-body interactions, and assess their performance with the same structural and many-body dependent properties used in the assessment of the DNN potential. Finally, in Sec. IV, we summarize our work and discuss possible future synergies between MB and DNN potentials.

Since the MB-pol potential of water has already been described in detail in the literature, we only overview here its salient features.^30–32 MB-pol was derived from the MBE that expresses the energy, E_N, of a system containing N (atomic or molecular) monomers as the sum of individual n-body energy contributions,

Here, ɛ^1B represents the distortion energy of an isolated monomer, such that ɛ^1B(i) = E(i) − E_eq(i), where E_eq(i) is the energy of the i-th monomer in its equilibrium geometry. The n-body energies, ɛ^nB, are defined recursively for 1 < n ≤ N by the expression

MB-pol approximates Eq. (2) as

The one-body term (ɛ^1B) is represented by the potential developed by Partridge and Schwenke.¹¹² The two-body term (ɛ^2B) describes four distinct contributions: permanent electrostatics, dispersion, 2B polarization, and 2B short-range interactions. The three-body term (ɛ^3B) describes two distinct contributions: 3B polarization and 3B short-range interactions. 2B and 3B short-range interactions are represented by the 2B and 3B permutationally invariant polynomial (PIPs)¹¹³ that were fitted in order for ɛ_2B and ɛ_3B to reproduce the 2B and 3B energies, respectively, calculated at the CCSD(T) level of theory in the complete basis set (CBS) limit.^30,31 2B and 3B polarization contributions are implicitly included in E_POL in Eq. (3), which represents classical many-body interactions at all orders through a polarization term. Further details of the MB-pol potential can be found in the original Refs. 30–32.

The DeePMD framework reads atomic positions and associated atom types as input features.⁶⁶ Neighbor information for each atom i is extracted from the input feature using a predefined cutoff radius (r_c) and stored as the coordinate difference of each ij atom pair into $Ri∈RNi×3$⁠, where N_i is the number of neighboring atoms. Each local feature is then mapped onto generalized coordinates $R̃i$⁠, as outlined in Ref. 71. A local embedding matrix, $Gi$⁠, is applied to each local feature $Ri$ in order to ensure rotation and permutation symmetry, while preserving translation symmetry. The resulting encoded feature matrix $Di∈RM1×M2$ takes the form

and is passed to a fully connected, feed-forward DNN, which maps it onto an “atomic energy” E_i.⁷¹ The total energy E is then calculated as the sum of all E_i, while the atomic forces F and virials Ξ are calculated from the derivative of the DNN with respect to the corresponding atomic positions. The DNN parameters are optimized by minimizing the loss function

where p_e, p_f, p_ξ are weighting factors, and ΔE, ΔF, and ΔΞ are the prediction errors for the reference energy, force, and virial values, respectively. The weighting factors p_e, p_f, and p_ξ are adjusted as the training progresses in order to improve the quality of the fit.

The DeePMD-based DNN potentials presented in this study were developed with the Deep Potential Smooth Edition (DeepPot-SE),⁷¹ following the procedure reported in Ref. 114, using 25, 50, and 100 neurons for the hidden embedding layers in the DeepPot-SE, while the submatrix of the embedding matrix uses 16 neurons. The distance cutoff was set to 6 Å, with a smoothing region of 0.5 Å. Each DNN potential is represented by a fully connected deep neural network with three layers of 240 neurons each.

Following Ref. 115, the training set for the DNN potential was constructed in an iterative fashion. Briefly, the training set comprises energies and forces of molecular configurations extracted from the MB-pol simulations of liquid water between 198 and 368 K, as reported in Ref. 96, as well as additional configurations extracted from simulations carried out with three successive iterations of the DNN potential. The final training set includes 94 770 configurations, each containing 256 molecules. All MB-pol reference data were computed using the MBX software package.¹¹⁶ Additional details about the training set are discussed in the supplementary material. To account for variations in training, validation, and testing errors, four distinct potentials (hereafter referred to as seed 1, seed 2, seed 3, and seed 4, respectively,) were trained using different random seeds. Only one of the four DNN potentials (seed 2) was then used in the MD simulations of liquid water and the vapor–liquid interface. Similarly, four distinct DeePMD-based potentials were trained on the expanded training sets containing vapor–liquid and cluster configurations, as described in Secs. III B and III C, respectively. Specific details about each of the four distinct DeePMD-based potentials trained on different training sets are reported in Table S3 of the supplementary material. Figure 1 shows the root-mean-square error (RMSE) curves of training and validation sets for the energies and forces per atom during the fitting process of the (seed 2) DNN potential. Overall, well-behaved learning curves are obtained for both quantities, with final errors of 0.01 kcal/mol and 1 kcal/mol Å on the energy and force validation errors, respectively. Similar errors have been reported for other state-of-the-art machine-learned potentials.¹¹⁷ 

We performed two sets of MD simulations to determine the ability of the DNN potential to reproduce both bulk and interfacial properties of liquid water calculated with MB-pol. The first set of simulations was carried out for a box containing 256 water molecules in the isothermal–isobaric (NPT) ensemble at 1 atm and in the temperature range between 198 and 368 K. The temperature was maintained using a global Nosé–Hoover thermostat chain of length 3, with a relaxation time of 0.05 ps, and the pressure was controlled by a global Nosé–Hoover barostat, with a relaxation time of 0.5 ps, which was thermostatted by a Nosé–Hoover thermostat chain of length 3. At each temperature, the last frame of the MB-pol trajectories reported in Ref. 96 was used as the initial configuration for the NPT simulations with the DNN potential. The second set of simulations was carried out in the canonical (NVT) ensemble between 400 and 575 K for a liquid slab of 512 water molecules in a box of dimensions 20 × 20 × 100 ų. The temperature was maintained using the same global Nosé–Hoover thermostat chain used in the NVT simulations. In both NPT and NVT simulations, the velocity Verlet algorithm was used to propagate the equations of motion with a time step of 0.5 fs, according to Ref. 118. All simulations were carried out using the DeePMD-kit¹¹⁴ plugin for Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).¹¹⁹ A complete set of input files is available on GitHub (see Data Availability).

Besides comparing the DNN and MB-pol radial distribution functions (RDFs), we also analyzed the ability of the DNN potential to describe the three-dimensional hydrogen-bond network by calculating the tetrahedral order parameter, q_tet,¹²⁰ 

Here, ψ_jk is the angle between the oxygen of the central water molecule and the oxygen atoms of two neighboring water molecules within a cutoff of 3.5 Å. When q_tet = 1, the water molecules are in a perfect tetrahedral arrangement, while q_tet = 0 represents the ideal gas limit.

In addition to the MD simulations for liquid water and the vapor–liquid interface, we also performed many-body decomposition analyses for two different sets of cluster structures. The first set consists of the first eight low-lying energy isomers of the water hexamer (Fig. 2), with geometries taken from Ref. 84. The hexamer occupies a special place in the description of many-body interactions in water because it is the smallest cluster, with low-lying isomers that display three-dimensional hydrogen-bonded arrangements similar to those found in liquid water and ice. The second set of clusters contains dimers and trimers extracted from the training sets used to fit the MB-pol 2B and 3B energy terms, respectively.^30,31

As a first step in assessing the ability of the DNN potential to reproduce MB-pol, we analyze various properties of liquid water. In Fig. 3(a), we show the temperature dependence of the liquid density from 198 to 368 K. In general, the DNN potential reproduces the MB-pol results over the entire temperature range, predicting similar temperatures for the density maximum and the minimum. A more quantitative analysis indicates that the DNN potential underestimates the MB-pol density by $∼1%$ in the 220–290 K range, while it predicts a slightly denser liquid as the temperature approaches the boiling point. The DNN curve also displays a less negative slope for temperatures above ∼320 K, which suggests that it is relatively more difficult for the DNN potential to reproduce MB-pol, as the liquid properties become more gas-like. To put the present comparison between the MB-pol and DNN potentials in context, we note that the density of liquid water at 298 K, predicted by an analogous DeePMD-based DNN potential, trained on the SCAN functional, was found to be ∼5% lower than the corresponding value, calculated from the reference ab initio molecular dynamics (AIMD) simulations.¹²¹ 

The comparison between the DNN and MB-pol values for the isothermal compressibility, as a function of temperature, is shown in Fig. 3(b). Similar to the density, the DNN values are in agreement with the MB-pol reference data, reproducing the steep increase of the isothermal compressibility below 250 K and predicting a maximum at ∼230 K, which is ∼7 K higher than the temperature predicted by MB-pol.⁹⁶ As in the case of the liquid density, Fig. 3(b) also indicates that the ability of the DNN potential to reproduce MB-pol somewhat deteriorates as the temperature approaches the boiling point. In particular, the DNN potential predicts a nearly constant value of the isothermal compressibility above 300 K, with no indication of a distinct minimum that is instead found in both experiment^99,122 and MB-pol simulations.^84,96

A comparison between the structural properties of liquid water predicted by the DNN and MB-pol potentials between 198 and 368 K at 1 atm is shown in Fig. 4. The oxygen–oxygen radial distribution functions (RDFs) calculated from the NPT simulations, carried out with the two potentials [Fig. 4(a)], are nearly indistinguishable in the 238–368 K range. Small deviations are found at deeply supercooled temperatures, which become more apparent when analyzing the corresponding distributions of the tetrahedral order parameter in Fig. 4(b). These deviations appear to be consistent with the shift of the isothermal compressibility maximum to a slightly higher temperature predicted by the DNN potential [Fig. 3(b)].

Previous studies demonstrated that MB-pol correctly predicts structural, thermodynamic, and spectroscopic properties of the vapor–liquid interface, including the surface tension, vapor pressure, vapor and liquid densities,^84,91 as well as sum frequency generation spectra.^89,90 To assess the ability of the DNN potential to reproduce properties that are not directly related to the MB-pol liquid configurations used during the training process, in Fig. 5, we analyze the surface tension and liquid–vapor equilibrium densities as a function of temperature. These comparisons show that both surface tension and equilibrium densities predicted by the DNN potential deviate significantly from the corresponding MB-pol reference values as the temperature increases. Interestingly, while the liquid density predicted by the DNN potential decreases upon increasing temperature, in qualitative agreement with the expected physical behavior, the vapor density remains effectively constant over the entire temperature range. Following Ref. 91, we estimated the critical temperature (T_c) and density (ρ_c) associated with DNN potential by fitting the vapor (ρ_v) and liquid (ρ_l) densities (Fig. 5), according to

Here, A and Δρ₀ are system-specific parameters to be adjusted in the fitting, and β = 0.326 is the critical exponent of the three-dimensional Ising model.¹²³ The DNN potential predicts T_c = 857 ± 17 K and ρ_c = 0.302 ± 0.002 g cm⁻³, which are significantly different from the MB-pol values of T_c = 639 ± 14 K and ρ_c = 0.34 ± 0.03 g cm⁻³, respectively. As a reference, the experimental values are T_c = 647 K and ρ_c = 0.32 g cm⁻³, respectively.¹²⁴ 

In an attempt to rationalize the different performance of the DNN potential in reproducing bulk and interfacial properties calculated with MB-pol, we investigated the ability of the DNN potential to correctly describe many-body interactions. By construction, MB-pol quantitatively reproduces each term of the MBE [Eq. (1)] calculated at the CCSD(T)/CBS level.^30,31 In this context, we have shown that a correct representation of each individual n-body contribution to the interaction energies is required in order for a water model to be both accurate and transferable across different thermodynamic state points.^{8,75,80,125–127}

Following previous studies,^84,125 in Fig. 6, we present a many-body decomposition analysis of the interaction energies of the first eight low-lying energy isomers of the hexamer cluster (Fig. 2). As mentioned in the Introduction, among water clusters, the hexamer occupies a special place because it is the smallest cluster, with low-lying isomers that exhibit three-dimensional structures, reminiscent of hydrogen-bonding arrangements found in liquid water and ice. In addition, the large number of isomers with similar interaction energies makes the hexamer the ideal benchmarking system for determining the accuracy of water models.⁸ To provide a general perspective on DeePMD-based DNN potentials for water, in Fig. 6, we analyze the performance of four distinct DNN potentials trained on the same training set described in Sec. II B, but initialized using different random seeds, with seed 2 corresponding to the DNN potential used in Figs. 3–5.

All four DNN potentials provide statistically equivalent training, validation, and testing errors (see Tables S2 and S3 of the supplementary material). Figure 6 shows that none of the four DNN potentials is capable of correctly reproducing individual n-body energies (n = 2–6) calculated with MB-pol, with significantly large errors, of the order of 10–20 kcal/mol, for two- and three-body contributions. Interestingly, these large errors compensate among the different n-body contributions in such a way that, when the n-body energies are added together, they result in interaction energies that are in relatively better agreement with the MB-pol values than the individual n-body contributions. By definition, the interaction energies are calculated as the difference between the energy of the cluster and the sum of the one-body energies of all six water molecules in the same distorted configurations as in the cluster. In this context, it should be noted that, besides MB-pol that reproduces the CCSD(T)/CBS reference energies of the hexamer isomers with chemical accuracy,⁸⁴ several modern, polarizable force fields predict the n-body and interaction energies of water clusters with significantly higher accuracy than the four DNN potentials examined here.¹²⁷ Direct comparisons between the n-body and interaction energies calculated with the four distinct DNN potentials, and the corresponding MB-pol reference values are shown in Fig. S3. Importantly, Fig. S3 shows that, besides displaying large errors, some of the DNN potentials (i.e., seed 1 and seed 4) also predict physically incorrect many-body contributions (e.g., positive three-body contributions), which indicates that, in their conventional implementation, DeePMD-based DNN potentials are not able to correctly disentangle the individual many-body contributions to the interaction energy of a given water system. Importantly, the inclusion of long-range effects through a classical electrostatic term does not improve the description of many-body energies, as shown in Figs. S4 and S5 of the supplementary material. It should be noted that this behavior is not specific to DeePMD-based DNN potentials, but appears to be common to other neural network potentials. For example, Figs. S6 and S7 of the supplementary material show that similar behavior is exhibited by Nequip-based potentials⁶⁵ trained on MB-pol. Interestingly, the Nequip-based potentials demonstrate superior accuracy in predicting the interaction energies of the water clusters, but also exhibit larger error compensation among different n-body energies, with errors on two- and three-body energies being as large as 20–30 kcal/mol.

In an attempt to improve the performance of the DNN potential on vapor–liquid equilibrium properties, we used active learning to incorporate vapor–liquid configurations extracted from simulations carried out with the DNN potential in the temperature range between 268 and 575 K. At the end of the active learning process, 2412 were added to the training set. The expanded training set was then used to train two potentials—DNN(VLE10) and DNN(VLE20)—with a 10% and 20% probability of selecting vapor–liquid configurations during training, respectively. Figure 7 shows that adding vapor–liquid configurations leads to more accurate predictions of both surface tension and vapor–liquid equilibrium densities. In particular, compared to the results obtained with the DNN potential, the surface tension predicted by both DNN(VLE10) and DNN(VLE20) shows the same temperature dependence as determined by MB-pol, although a systematic deviation from the reference values is still observed at all temperatures. Similarly, adding vapor–liquid configurations to the training set improves the ability of the DNN(VLE10) and DNN(VLE20) potentials to describe the equilibrium densities of both vapor and liquid phases. While both DNN(VLE10) and DNN(VLE20) potentials quantitatively reproduce the MB-pol vapor densities over the entire temperature range, the predicted liquid densities, however, increasingly deviate from the MB-pol reference values as the temperature increases. As a consequence, the critical point is still overestimated by both potentials, with DNN(VLE10) predicting T_c = 718 ± 4 K and ρ_c = 0.359 ± 0.003 g cm⁻³ and DNN(VLE20) predicting T_c = 674 ± 6 K and ρ_c = 0.347 ± 0.005 g cm⁻³. Adding vapor–liquid configurations to the training set was reported to enable simulations of “water along its liquid/vapor coexistence line with unprecedented precision.”⁵³ Inspection of Fig. 3 of Ref. 53, however, indicates that relatively large deviations, [similar to those found for DNN(VLE10) and DNN(VLE20) in Fig. 7], exist between the vapor and liquid densities predicted by the neural network potential used in the simulations and the corresponding reference RPBE-D3 values, which, when extrapolated, lead to very different estimates for the critical point.¹²⁸ 

To assess the ability of DNN(VLE10) and DNN(VLE20) to reproduce properties that do not directly depend on the coexistence between vapor and liquid phases, we examined the performance of both potentials on the same bulk and cluster properties used in Sec. III A to determine the accuracy of the DNN potential. Figure 8 shows the temperature dependence of the density and isothermal compressibility of liquid water predicted by DNN(VLE10) and DNN(VLE20). While both potentials are able to qualitatively reproduce the MB-pol trends, a comparison with the DNN results reported in Fig. 3 indicates that the inclusion of vapor–liquid configurations to the training set deteriorates the ability of the DeePMD-based potentials to reproduce the bulk properties. This is further confirmed by the analyses of the liquid density, RDFs, and q_tet distributions shown for the DNN(VLE20) potential in Figs. S17 and S18 of the supplementary material.

Finally, Fig. 9 reports the many-body decomposition analysis of the interaction energies of the hexamer isomers (Fig. 2) carried out with the DNN(VLE20) potential. The corresponding analysis carried out with DNN(VLE10) is reported in the supplementary material in Fig. S8. As for the DNN potential, we used four different seeds to develop four distinct DNN(VLE20) potentials that were trained on the expanded MB-pol training set containing vapor–liquid configurations. Seed 4 corresponds to the DNN(VLE20) potential used in the comparisons shown in Figs. 7 and 8. As in the case of DNN in Fig. 6, none of the DNN(VLE20) potentials is able to correctly reproduce the reference MB-pol many-body energies, with errors that are of the order of ∼10 kcal/mol for two-, three-, and four-body energies. Similar poor performance on the many-body decomposition analysis is exhibited by the DNN(VLE10) potential in Fig. S8 of the supplementary material. Analyses analogous to those shown in Fig. S3 for the DNN potential are reported in Figs. S9 and S10 for the DNN(VLE10) and DNN(VLE20) potentials, respectively, which lead to similar conclusions, i.e., both DNN(VLE10) and DNN(VLE20) predict physically incorrect many-body energies.

The analyses presented in this section demonstrate that, while the description of the vapor–liquid equilibrium properties can be improved by adding vapor–liquid configurations to the original DNN training set, this improvement is achieved at the cost of a less accurate representation of the bulk properties. Importantly, as in the case of the DNN potentials, the DNN(VLE10) and DNN(VLE20) potentials are unable to correctly capture the physics of many-body interactions in water.

Since the inability of a water model to correctly represent many-body contributions to the underlying molecular interactions appears to be correlated to the lack of transferability of the model across different thermodynamic state points,⁸ we investigated the possibility of “encoding” the many-body effects in the DNN potentials within the DeePMD framework. To this end, we supplemented the original training set used for developing the DNN potentials discussed in Sec. III A with a set of gas-phase clusters, including monomers, dimers, trimers, and tetramers, which provides direct information about the low-order and most important terms (i.e., one-body to four-body terms) of the MBE in Eq. (1). We then used the expanded training set to train three different DeePMD-based potentials, referred to as DNN(MB4), DNN(MB10), and DNN(MB20), with 4%, 10%, and 20% probability of selecting gas-phase cluster configurations during the training process, respectively. Specific details about the composition of the extended training set are provided in Sec. S1 of the supplementary material.

Figure 10 reports the same analysis reported in Fig. 6 for the DNN potential and shows the errors relative to the MB-pol reference values for n-body and interaction energies of the first eight isomers of the water hexamer (Fig. 2), calculated with four distinct versions of the DNN(MB20) potential that were trained on the same expanded training set, but initialized with four distinct seeds. Seed 4 corresponds to the DNN(MB20) potential used in the comparisons shown in Figs. 11 and 12. Analogous analyses carried out with the DNN(MB4), and DNN(MB10) potentials are reported in Figs. S11 and S12 of the supplementary material, respectively. The addition of the monomer, dimer, trimer, and tetramer configurations clearly allows the DNN(MB) potentials to become “aware” of the existence of distinct many-body contributions to the interactions energies, as demonstrated by the relatively smaller errors displayed by the DNN(MB20) two-body, and three-body energies compared to the corresponding errors associated with the DNN potentials in Fig. 6. Similar trends are observed in Figs. S13–S15, which show direct comparisons of individual many-body energies and interaction energies, calculated with the DNN(MB4), DNN(MB10), and DNN(MB20) potentials, respectively. Additional analyses of the error distributions associated with two-body and three-body energies calculated for the dimer and trimer configurations of the cluster training set are reported in Fig. S16 and demonstrate that the DNN(MB4), DNN(MB10), and DNN(MB20) potentials significantly improve on the DNN potential in their ability to represent many-body interactions in water. It should, however, be noted that the DNN(MB) potentials still rely on significant error compensation among individual n-body energy contributions, to minimize the error in the interaction energies of the hexamer isomers. The observed error compensation indicates that, similar to the DNN, DNN(VLE10), and DNN(VLE20) potentials, the DNN(MB) potentials are unable to “learn” that the interaction energy of an N-body system containing N water molecules is given by the sum of distinct n-body energy contributions (with n = 2 − N). To put things in perspective, the errors associated with the DNN(MB) predictions for each n-body energy contribution to the interaction energies of the hexamer isomers, in particular at the two-body and three-body levels, are still appreciably larger than those displayed by state-of-the-art, polarizable force fields for water.¹²⁷ 

Having demonstrated that extending the training set by adding the monomer, dimer, trimer, and tetramer configurations allows the DNN(MB) potentials to recover a more balanced representation of many-body interactions, we next assess the ability of the DNN(MB4), DNN(MB10), and DNN(MB20) potentials to reproduce vapor–liquid equilibrium properties that were poorly predicted by the DNN potentials (Fig. 5). Figure 11 shows that all three DNN(MB) potentials more closely reproduce the MB-pol trends for the surface tension and the equilibrium densities of both vapor and liquid phases over the entire temperature range examined in this study compared to the DNN potential. The critical parameters predicted by the DNN(MB4), DNN(MB10), and DNN(MB20) potentials are T_c = 655 ± 2 K and ρ_c = 0.325 ± 0.002 g cm⁻³, T_c = 605 ± 10 K and ρ_c = 0.32 ± 0.01 g cm⁻³, and T_c = 660 ± 6 K and ρ_c = 0.338 ± 0.005 g cm⁻³, respectively, which are in better agreement with the MB-pol values (T_c = 639 ± 14 K and ρ_c = 0.34 ± 0.03 g cm⁻³) than with the results obtained not only with the DNN potential, but also with the DNN(VLE10) and DNN(VLE20) potentials. The structural differences at the vapor–liquid equilibrium between DNN and DNN(MB20) are further highlighted in Fig. S17, which shows the density profiles predicted by the two potentials at different temperatures. In particular, the interface structure predicted by DNN(MB20) is significantly different from that predicted by DNN and is in close agreement with the MB-pol results reported in Ref. 91.

Despite being able to provide more accurate estimates of the actual MB-pol critical point, Fig. 12 shows that the DNN(MB4), DNN(MB10), and DNN(MB20) potentials display a higher degree of variability in their predictions of the liquid density at high temperatures than the DNN(VLE10) and DNN(VLE20) potentials. This high variability at high temperatures, which is also displayed by the DNN potentials, can be traced back to the lack of explicit vapor–liquid configurations in the corresponding training sets, which, in turn, highlights the difficulty for DeePMD-based potentials to be transferable across different phases over a wide range of thermodynamic conditions.

The last question that remains to be addressed is whether the improved ability to represent many-body interactions and predict vapor–liquid equilibrium properties still allows the DNN(MB4), DNN(MB10), and DNN(MB20) potentials to accurately reproduce the liquid properties calculated with MB-pol. To this end, Fig. 13 shows comparisons between the temperature dependence of the density and isothermal compressibility calculated with the DNN(MB4), DNN(MB10), and DNN(MB20) potentials and the corresponding MB-pol reference values. The DNN(MB4), DNN(MB10), and DNN(MB20) potentials effectively predict indistinguishable (within statistical error) trends for both density and isothermal compressibility, which are in qualitative agreement with the MB-pol reference values. Specifically, they correctly predict a minimum at ∼300 K in the isothermal compressibility, which is, instead, absent in the DNN, DNN(VLE10), and DNN(VLE20) potentials. This suggests that, by being able to more accurately represent many-body interactions, the DNN(MB4), DNN(MB10), and DNN(MB20) potentials display a higher degree of transferability to the gas phase at ambient conditions. As in the case of the DNN(VLE10) and DNN(VLE20) potentials, the comparison between the results of Figs. 3 and 13, however, indicates that the addition of configurations different from bulk configurations (in this case, monomer, dimer, trimer, and tetramer configurations) to the training set overall deteriorates the ability of the DNN(MB4), DNN(MB10), and DNN(MB20) potentials to reproduce bulk properties calculated with MB-pol. Despite these differences, the liquid structure predicted by the DNN(MB4), DNN(MB10), and DNN(MB20) potentials is in close agreement with that of MB-pol, as demonstrated by the comparisons of the RDFs and q_tet distributions calculated with DNN(MB20), which are shown in Fig. S18 of the supplementary material.

In this study, we analyzed the performance and degree of transferability of a DeePMD-based DNN potential for water, trained on MB-pol reference configurations extracted from MD simulations of liquid water carried out from 198 to 368 K at 1 atm. We found that the DNN potential is able to reliably reproduce structural and thermodynamic properties of liquid water, as predicted by MB-pol, from the boiling point, down to deeply supercooled temperatures. However, while MB-pol exhibits remarkable accuracy from the gas to the condensed phase, the DNN potential does not share the same high level of transferability across phases. In particular, we found that the DNN potential is not able to accurately describe vapor–liquid equilibrium properties. More importantly, a many-body decomposition analysis of the interaction energies of the hexamer isomers indicates that the DNN potential is not able to correctly “learn” many-body interactions and effectively relies on error compensation among individual many-body energy contributions to reproduce the interaction energy of a given N-body system containing N water molecules.

To improve the performance of the DNN potential in vapor–liquid equilibrium properties, we expanded the initial DNN training set of bulk configurations, by adding configurations extracted from vapor–liquid equilibrium simulations carried out with MB-pol. While the new DNN(VLE10) and DNN(VLE20) potentials improve the description of the surface tension and equilibrium densities of both vapor and liquid phases, they predict less accurately bulk properties and are unable to correctly reproduce the individual many-body contributions to the interaction energies.

In an attempt to explicitly encode many-body interactions with the DNN potential, we also expanded the initial DNN training set by adding water monomer, dimer, trimer, and tetramer configurations, which provide direct information on the most important many-body contributions (i.e., one-body to four-body contributions) to the interaction energies in water systems. By improving the description of individual many-body contributions, the new DNN(MB4), DNN(MB10), and DNN(MB20) potentials are also able to reliably reproduce the vapor–liquid equilibrium properties predicted by MB-pol. We found, however, that all three potentials exhibit a high degree of variability in predicting the liquid density at high temperatures due to the lack of representative vapor–liquid configurations in their training sets, which limits their transferability over a wide range of thermodynamic conditions. Moreover, the improvement in the description of many-body interactions comes at the cost of a poorer representation of the liquid properties.

Although DeePMD-based potentials are intrinsically many-body in their functional form, our analyses show that they do not necessarily correctly represent the underlying many-body physics of the reference potentials. This suggests that some caution should be exercised when using DeePMD-based DNN potentials to predict thermodynamic properties for state points that are not explicitly and thoroughly included in the training sets. Although our study focuses on water, similar behavior is likely to be found in DeePMD-based DNN potentials for other molecular systems, including aqueous solutions, as well as molecular fluids and solids. In this context, we hope that our results can stimulate further developments of new training procedures and neural network architectures, capable of correctly capturing the physics of the many-body interactions in molecular systems.

With this caveat in mind, the computational efficiency provided by the DeePMD framework suggests that large-scale CCSD(T)-level MD simulations are possible by training DeePMD-based DNN potentials on data-driven many-body potentials derived from the MBE calculated at the CCSD(T) level of theory, such as MB-pol. However, for this to hold, the thermodynamic state points of interest in the DNN simulations must be adequately represented in the training sets generated using the reference data-driven many-body potentials. This suggests that a DNN potential trained on an extensive training set, including molecular configurations extracted from MB-pol simulations carried out over a wide range of thermodynamic conditions, is well suited for exploring the rich phase diagram of water,¹²⁹ particularly in the so-called “no man’s land” region at low temperature, which has been proven difficult to probe experimentally.^99,101,102

See the supplementary material for details about the DNN potentials and training procedure and about additional analyses of the performance of the different DNN potentials on bulk, vapor–liquid, and many-body properties of water.

We thank Maria Muniz, Athanassios Panagiotopoulos, and Vinicius Cruzeiro for stimulating discussions at the early stage of this research. This research was supported by the Air Force Office of Scientific Research, under Grant No. FA9550-20-1-0351, and used the computational resources of the Department of Defense High Performance Computing Modernization Program (HPCMP), as well as the Triton Shared Computing Cluster (TSCC) at the San Diego Supercomputer Center (SDSC). This work was supported in part by the UC Southern California Hub, with funding from the UC National Laboratories division of the University of California Office of the President.

The authors have no conflicts to disclose.

Y.Z., A.C., and S.L.B. contributed equally to this work.

Yaoguang Zhai: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Alessandro Caruso: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Sigbjørn Løland Bore: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Zhishang Luo: Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Francesco Paesani: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

Training sets, DNN models, MD input, analysis scripts, and key results are publicly available on Zenodo (https://doi.org/10.5281/zenodo.7577034). MD trajectories are available from the corresponding author upon request.