Deep neural network (DNN) potentials have recently gained popularity in computer simulations of a wide range of molecular systems, from liquids to materials. In this study, we explore the possibility of combining the computational efficiency of the DeePMD framework and the demonstrated accuracy of the MB-pol data-driven, many-body potential to train a DNN potential for large-scale simulations of water across its phase diagram. We find that the DNN potential is able to reliably reproduce the MB-pol results for liquid water, but provides a less accurate description of the vapor–liquid equilibrium properties. This shortcoming is traced back to the inability of the DNN potential to correctly represent many-body interactions. An attempt to explicitly include information about many-body effects results in a new DNN potential that exhibits the opposite performance, being able to correctly reproduce the MB-pol vapor–liquid equilibrium properties, but losing accuracy in the description of the liquid properties. These results suggest that DeePMD-based DNN potentials are not able to correctly “learn” and, consequently, represent many-body interactions, which implies that DNN potentials may have limited ability to predict the properties for state points that are not explicitly included in the training process. The computational efficiency of the DeePMD framework can still be exploited to train DNN potentials on data-driven many-body potentials, which can thus enable large-scale, “chemically accurate” simulations of various molecular systems, with the caveat that the target state points must have been adequately sampled by the reference data-driven many-body potential in order to guarantee a faithful representation of the associated properties.

Molecular mechanics (MM) force fields (FFs)1,2 have been the workhorse in computational chemistry since the early days of Monte Carlo (MC)3 and molecular dynamics (MD) simulations.4 Continued progress in hardware technologies,5 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 FFs6–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 energy21 calculated using correlated electronic structure methods. When applied to aqueous systems22–35 and molecular fluids,36–38 these many-body (MB) potentials exhibit unprecedented accuracy, enabling predictive simulations from the gas to the condensed phases.39 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.40 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.66 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.67 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),73 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.80 

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 simulations81,82 due to its rich phase behavior, characterized by several anomalies.83 As a representative state-of-the-art MB potential, we selected MB-pol30–32 due to its demonstrated ability to correctly predict the properties of water across the entire phase diagram,84 including gas-phase clusters,85–87 liquid water,88 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.96 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 DeePMD49,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,107 and metals,108–110 and enabled MD simulations with up to 10 × 109 atoms.111 

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, EN, of a system containing N (atomic or molecular) monomers as the sum of individual n-body energy contributions,

EN(1,,N)= i=1Nε1B(i)+i=1Nj>iNε2B(i,j)+i=1Nj>iNk>j>iNε3B(i,j,k)++εNB(1,,N).
(1)

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

εnB= En(1,,n)i=1Nε1B(i)i=1Ni<jNε2B(i,j)i=1Ni<jNi<j<kε3B(i,j,k)i<j<k<Nε(n1)B(i,j,k,,n1).
(2)

MB-pol approximates Eq. (2) as

EN(r1,,rN)=i=1Nε1B(i)+i>jNε2B(i,j)+i>j>kNε3B(i,j,k)+EPOL.
(3)

The one-body term (ɛ1B) is represented by the potential developed by Partridge and Schwenke.112 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)113 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 EPOL 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.66 Neighbor information for each atom i is extracted from the input feature using a predefined cutoff radius (rc) and stored as the coordinate difference of each ij atom pair into RiRNi×3, where Ni 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 DiRM1×M2 takes the form

Di=(Gi1)TR̃i(R̃i)TGi2
(4)

and is passed to a fully connected, feed-forward DNN, which maps it onto an “atomic energy” Ei.71 The total energy E is then calculated as the sum of all Ei, 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

L(pe,pf,pξ)=peNΔE2+pf3NΔF2+pξ9NΔΞ2,
(5)

where pe, pf, pξ are weighting factors, and ΔE, ΔF, and ΔΞ are the prediction errors for the reference energy, force, and virial values, respectively. The weighting factors pe, pf, 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),71 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.116 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.117 

FIG. 1.

Variation of the DNN training and validation RMSEs per atom, relative to the MB-pol values of the energy and force, as a function of the number of training steps. For visual clarity, we show values averaged over 200 training steps.

FIG. 1.

Variation of the DNN training and validation RMSEs per atom, relative to the MB-pol values of the energy and force, as a function of the number of training steps. For visual clarity, we show values averaged over 200 training steps.

Close modal

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 Å3. 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-kit114 plugin for Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).119 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, qtet,120 

qtet=138j=13k=j+14cos(ψjk)+13.
(6)

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 qtet = 1, the water molecules are in a perfect tetrahedral arrangement, while qtet = 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

FIG. 2.

Structures of the first eight low-lying energy isomers of the water hexamer used in the analysis of the interaction and many-body energies. The Cartesian coordinates of each isomer were taken from Ref. 84.

FIG. 2.

Structures of the first eight low-lying energy isomers of the water hexamer used in the analysis of the interaction and many-body energies. The Cartesian coordinates of each isomer were taken from Ref. 84.

Close modal

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.121 

FIG. 3.

Temperature dependence of the density (a) and isothermal compressibility (b), calculated from the present NPT simulations, carried out with the DNN potential (seed 2) at 1 atm (green), compared to the reference MB-pol values from Ref. 96 (blue). The associated shaded areas indicate 95% confidence intervals of the averages. Figure S1 displays the time dependence of the densities calculated along the NPT trajectories, carried out with the DNN potential, at each temperature.

FIG. 3.

Temperature dependence of the density (a) and isothermal compressibility (b), calculated from the present NPT simulations, carried out with the DNN potential (seed 2) at 1 atm (green), compared to the reference MB-pol values from Ref. 96 (blue). The associated shaded areas indicate 95% confidence intervals of the averages. Figure S1 displays the time dependence of the densities calculated along the NPT trajectories, carried out with the DNN potential, at each temperature.

Close modal

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.96 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 experiment99,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)].

FIG. 4.

(a) Oxygen–oxygen radial distribution functions and (b) tetrahedral order parameter distributions, calculated from the present NPT simulations, carried out with the DNN potential (seed 2) at 1 atm (green), compared to the reference MB-pol values from Ref. 96 (blue).

FIG. 4.

(a) Oxygen–oxygen radial distribution functions and (b) tetrahedral order parameter distributions, calculated from the present NPT simulations, carried out with the DNN potential (seed 2) at 1 atm (green), compared to the reference MB-pol values from Ref. 96 (blue).

Close modal

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 (Tc) and density (ρc) associated with DNN potential by fitting the vapor (ρv) and liquid (ρl) densities (Fig. 5), according to

ρl+ρv2=ρc+A(TcT),
(7)
ρlρv2=Δρ01T/Tcβ.
(8)

Here, A and Δρ0 are system-specific parameters to be adjusted in the fitting, and β = 0.326 is the critical exponent of the three-dimensional Ising model.123 The DNN potential predicts Tc = 857 ± 17 K and ρc = 0.302 ± 0.002 g cm−3, which are significantly different from the MB-pol values of Tc = 639 ± 14 K and ρc = 0.34 ± 0.03 g cm−3, respectively. As a reference, the experimental values are Tc = 647 K and ρc = 0.32 g cm−3, respectively.124 

FIG. 5.

Surface tension (a) and vapor–liquid equilibrium densities (b), calculated from the present NVT simulations of a water slab, carried out with the DNN potential (seed 2, green), compared to the reference MB-pol values from Ref. 91 (blue).

FIG. 5.

Surface tension (a) and vapor–liquid equilibrium densities (b), calculated from the present NVT simulations of a water slab, carried out with the DNN potential (seed 2, green), compared to the reference MB-pol values from Ref. 91 (blue).

Close modal

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.8 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. 35.

FIG. 6.

Many-body decomposition analysis for the eight low-lying energy isomers of the water hexamer (Fig. 2), calculated with four distinct DNN potentials trained on the same MB-pol training data using four different seeds to initialize the fitting process. Panels (a)–(e) show the errors associated with the n-body energies (n = 2–6) calculated with the DNN potentials relative to the corresponding MB-pol values. Panel (f) shows the errors associated with the interaction energies, calculated with the DNN potentials, relative to the corresponding MB-pol values. The DNN potential with seed 2 is used in the comparisons shown in Figs. 35.

FIG. 6.

Many-body decomposition analysis for the eight low-lying energy isomers of the water hexamer (Fig. 2), calculated with four distinct DNN potentials trained on the same MB-pol training data using four different seeds to initialize the fitting process. Panels (a)–(e) show the errors associated with the n-body energies (n = 2–6) calculated with the DNN potentials relative to the corresponding MB-pol values. Panel (f) shows the errors associated with the interaction energies, calculated with the DNN potentials, relative to the corresponding MB-pol values. The DNN potential with seed 2 is used in the comparisons shown in Figs. 35.

Close modal

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,84 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.127 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 potentials65 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 Tc = 718 ± 4 K and ρc = 0.359 ± 0.003 g cm−3 and DNN(VLE20) predicting Tc = 674 ± 6 K and ρc = 0.347 ± 0.005 g cm−3. Adding vapor–liquid configurations to the training set was reported to enable simulations of “water along its liquid/vapor coexistence line with unprecedented precision.”53 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.128 

FIG. 7.

Surface tension (a) and vapor–liquid equilibrium densities (b) calculated from the present NVT simulations of a water slab, carried out with the DNN(VLE10) (seed 4, purple) and DNN(VLE20) (seed 4, pink), compared to the reference MB-pol values from Ref. 91 (blue).

FIG. 7.

Surface tension (a) and vapor–liquid equilibrium densities (b) calculated from the present NVT simulations of a water slab, carried out with the DNN(VLE10) (seed 4, purple) and DNN(VLE20) (seed 4, pink), compared to the reference MB-pol values from Ref. 91 (blue).

Close modal

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 qtet distributions shown for the DNN(VLE20) potential in Figs. S17 and S18 of the supplementary material.

FIG. 8.

Temperature dependence of the density (a) and isothermal compressibility (b), calculated from NPT simulations carried out with the DNN(VLE10) (seed 4, purple) and DNN(VLE20) (seed 4, pink) potentials at 1 atm. Also shown for reference are the corresponding MB-pol values from Ref. 96 (blue). The associated shaded areas indicate 95% confidence intervals of the averages.

FIG. 8.

Temperature dependence of the density (a) and isothermal compressibility (b), calculated from NPT simulations carried out with the DNN(VLE10) (seed 4, purple) and DNN(VLE20) (seed 4, pink) potentials at 1 atm. Also shown for reference are the corresponding MB-pol values from Ref. 96 (blue). The associated shaded areas indicate 95% confidence intervals of the averages.

Close modal

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.

FIG. 9.

Many-body decomposition analysis for the eight low-lying energy isomers of the water hexamer (Fig. 2), calculated with four distinct DNN(VLE20) potentials that were trained on the expanded MB-pol training set containing vapor–liquid configurations using four different seeds to initialize the fitting process. Panels (a)–(e) show the errors associated with the n-body energies (n = 2–6) calculated with the DNN(VLE20) potentials relative to the corresponding MB-pol values. Panel (f) shows the errors associated with the interaction energies calculated with the DNN(VLE20) potentials relative to the corresponding MB-pol values. The DNN(VLE20) potential with seed 4 is used in the comparisons shown in Figs. 7 and 8.

FIG. 9.

Many-body decomposition analysis for the eight low-lying energy isomers of the water hexamer (Fig. 2), calculated with four distinct DNN(VLE20) potentials that were trained on the expanded MB-pol training set containing vapor–liquid configurations using four different seeds to initialize the fitting process. Panels (a)–(e) show the errors associated with the n-body energies (n = 2–6) calculated with the DNN(VLE20) potentials relative to the corresponding MB-pol values. Panel (f) shows the errors associated with the interaction energies calculated with the DNN(VLE20) potentials relative to the corresponding MB-pol values. The DNN(VLE20) potential with seed 4 is used in the comparisons shown in Figs. 7 and 8.

Close modal

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,8 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.127 

FIG. 10.

Many-body decomposition analysis for the eight low-lying energy isomers of the water hexamer (Fig. 2), calculated with four distinct DNN(MB20) potentials that were trained on the expanded MB-pol training set containing cluster configurations using four different seeds to initialize the fitting process. Panels (a)–(e) show the errors associated with the n-body energies (n = 2–6) calculated with the DNN(MB20) potentials relative to the corresponding MB-pol values. Panel (f) shows the errors associated with the interaction energies calculated with the DNN(MB20) potentials relative to the corresponding MB-pol values. The DNN(MB20) potential with seed 4 is used in the comparisons shown in Figs. 11 and 13.

FIG. 10.

Many-body decomposition analysis for the eight low-lying energy isomers of the water hexamer (Fig. 2), calculated with four distinct DNN(MB20) potentials that were trained on the expanded MB-pol training set containing cluster configurations using four different seeds to initialize the fitting process. Panels (a)–(e) show the errors associated with the n-body energies (n = 2–6) calculated with the DNN(MB20) potentials relative to the corresponding MB-pol values. Panel (f) shows the errors associated with the interaction energies calculated with the DNN(MB20) potentials relative to the corresponding MB-pol values. The DNN(MB20) potential with seed 4 is used in the comparisons shown in Figs. 11 and 13.

Close modal
FIG. 11.

Surface tension (a) and vapor–liquid equilibrium densities (b), calculated from NVT simulations of a water slab, carried out with the DNN(MB4) (seed 4), DNN(MB10) (seed 4), and DNN(MB20) (seed 4) potentials, compared to the reference MB-pol values from Ref. 91 (blue).

FIG. 11.

Surface tension (a) and vapor–liquid equilibrium densities (b), calculated from NVT simulations of a water slab, carried out with the DNN(MB4) (seed 4), DNN(MB10) (seed 4), and DNN(MB20) (seed 4) potentials, compared to the reference MB-pol values from Ref. 91 (blue).

Close modal
FIG. 12.

Vapor–liquid equilibrium densities calculated using the same four variants of each of the DNN, DNN(VLE10), DNN(VLE20), DNN(MB4), DNN(MB10), and DNN(MB20) potentials used in the analyses of Figs. 6, 9, and 10, and Figs. S8, S11, and S12, respectively. The densities are compared with the reference MB-pol values from Ref. 91 (blue dashed line) at 400 K (a), 500 K (b), and 575 K (c).

FIG. 12.

Vapor–liquid equilibrium densities calculated using the same four variants of each of the DNN, DNN(VLE10), DNN(VLE20), DNN(MB4), DNN(MB10), and DNN(MB20) potentials used in the analyses of Figs. 6, 9, and 10, and Figs. S8, S11, and S12, respectively. The densities are compared with the reference MB-pol values from Ref. 91 (blue dashed line) at 400 K (a), 500 K (b), and 575 K (c).

Close modal

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 Tc = 655 ± 2 K and ρc = 0.325 ± 0.002 g cm−3, Tc = 605 ± 10 K and ρc = 0.32 ± 0.01 g cm−3, and Tc = 660 ± 6 K and ρc = 0.338 ± 0.005 g cm−3, respectively, which are in better agreement with the MB-pol values (Tc = 639 ± 14 K and ρc = 0.34 ± 0.03 g cm−3) 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 qtet distributions calculated with DNN(MB20), which are shown in Fig. S18 of the supplementary material.

FIG. 13.

Temperature dependence of the density (a) and isothermal compressibility (b) calculated from NPT simulations carried out with the DNN(MB4) (seed 4), DNN(MB10) (seed 4), and DNN(MB20) (seed 4) potentials at 1 atm, compared to the reference MB-pol values from Ref. 96 (blue). The associated shaded areas indicate 95% confidence intervals of the averages. Figure S2 shows the density fluctuations along the NPT trajectories carried out with the DNN(MB20) potentials at each temperature.

FIG. 13.

Temperature dependence of the density (a) and isothermal compressibility (b) calculated from NPT simulations carried out with the DNN(MB4) (seed 4), DNN(MB10) (seed 4), and DNN(MB20) (seed 4) potentials at 1 atm, compared to the reference MB-pol values from Ref. 96 (blue). The associated shaded areas indicate 95% confidence intervals of the averages. Figure S2 shows the density fluctuations along the NPT trajectories carried out with the DNN(MB20) potentials at each temperature.

Close modal

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,129 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.

1.
S.
Lifson
and
A.
Warshel
, “
Consistent force field for calculations of conformations, vibrational spectra, and enthalpies of cycloalkane and n-alkane molecules
,”
J. Chem. Phys.
49
,
5116
5129
(
1968
).
2.
A.
Warshel
and
S.
Lifson
, “
Consistent force field calculations. II. Crystal structures, sublimation energies, molecular and lattice vibrations, molecular conformations, and enthalpies of alkanes
,”
J. Chem. Phys.
53
,
582
594
(
1970
).
3.
N.
Metropolis
,
A. W.
Rosenbluth
,
M. N.
Rosenbluth
,
A. H.
Teller
, and
E.
Teller
, “
Equation of state calculations by fast computing machines
,”
J. Chem. Phys.
21
,
1087
1092
(
1953
).
4.
B. J.
Alder
and
T. E.
Wainwright
, “
Phase transition for a hard sphere system
,”
J. Chem. Phys.
27
,
1208
1209
(
1957
).
5.
M. S.
Gordon
,
G.
Barca
,
S. S.
Leang
,
D.
Poole
,
A. P.
Rendell
,
J. L.
Galvez Vallejo
, and
B.
Westheimer
, “
Novel computer architectures and quantum chemistry
,”
J. Phys. Chem. A
124
,
4557
4582
(
2020
).
6.
J. W.
Ponder
,
C.
Wu
,
P.
Ren
,
V. S.
Pande
,
J. D.
Chodera
,
M. J.
Schnieders
,
I.
Haque
,
D. L.
Mobley
,
D. S.
Lambrecht
,
R. A.
DiStasio
, Jr.
,
M.
Head-Gordon
,
G. N. I.
Clark
,
M. E.
Johnson
, and
T.
Head-Gordon
, “
Current status of the AMOEBA polarizable force field
,”
J. Phys. Chem. B
114
,
2549
2564
(
2010
).
7.
K.
Vanommeslaeghe
and
A. D.
MacKerell
, Jr.
, “
CHARMM additive and polarizable force fields for biophysics and computer-aided drug design
,”
Biochim. Biophys. Acta, Gen. Subj.
1850
,
861
871
(
2015
).
8.
G. A.
Cisneros
,
K. T.
Wikfeldt
,
L.
Ojamäe
,
J.
Lu
,
Y.
Xu
,
H.
Torabifard
,
A. P.
Bartók
,
G.
Csányi
,
V.
Molinero
, and
F.
Paesani
, “
Modeling molecular interactions in water: From pairwise to many-body potential energy functions
,”
Chem. Rev.
116
,
7501
7528
(
2016
).
9.
Z.
Jing
,
C.
Liu
,
S. Y.
Cheng
,
R.
Qi
,
B. D.
Walker
,
J.-P.
Piquemal
, and
P.
Ren
, “
Polarizable force fields for biomolecular simulations: Recent advances and applications
,”
Annu. Rev. Biophys.
48
,
371
394
(
2019
).
10.
V. S.
Inakollu
,
D. P.
Geerke
,
C. N.
Rowley
, and
H.
Yu
, “
Polarisable force fields: What do they add in biomolecular simulations?
,”
Curr. Opin. Struct. Biol.
61
,
182
190
(
2020
).
11.
S. L.
Mayo
,
B. D.
Olafson
, and
W. A.
Goddard
, “
DREIDING: A generic force field for molecular simulations
,”
J. Phys. Chem.
94
,
8897
8909
(
1990
).
12.
A. K.
Rappé
,
C. J.
Casewit
,
K. S.
Colwell
,
W. A.
Goddard
 III
, and
W. M.
Skiff
, “
UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations
,”
J. Am. Chem. Soc.
114
,
10024
10035
(
1992
).
13.
J.
Wang
,
R. M.
Wolf
,
J. W.
Caldwell
,
P. A.
Kollman
, and
D. A.
Case
, “
Development and testing of a general AMBER force field
,”
J. Comput. Chem.
25
,
1157
1174
(
2004
).
14.
K.
Vanommeslaeghe
,
E.
Hatcher
,
C.
Acharya
,
S.
Kundu
,
S.
Zhong
,
J.
Shim
,
E.
Darian
,
O.
Guvench
,
P.
Lopes
,
I.
Vorobyov
, and
A.
MacKerell
, Jr.
, “
CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields
,”
J. Comput. Chem.
31
,
671
690
(
2010
).
15.
T. B.
Adler
,
G.
Knizia
, and
H.-J.
Werner
, “
A simple and efficient CCSD(T)-F12 approximation
,”
J. Chem. Phys.
127
,
221106
(
2007
).
16.
G.
Knizia
,
T. B.
Adler
, and
H.-J.
Werner
, “
Simplified CCSD(T)-F12 methods: Theory and benchmarks
,”
J. Chem. Phys.
130
,
054104
(
2009
).
17.
J. G.
Hill
,
K. A.
Peterson
,
G.
Knizia
, and
H.-J.
Werner
, “
Extrapolating MP2 and CCSD explicitly correlated correlation energies to the complete basis set limit with first and second row correlation consistent basis sets
,”
J. Chem. Phys.
131
,
194105
(
2009
).
18.
U.
Góra
,
R.
Podeszwa
,
W.
Cencek
, and
K.
Szalewicz
, “
Interaction energies of large clusters from many-body expansion
,”
J. Comput. Phys.
135
,
224102
(
2011
).
19.
D.
Manna
,
M. K.
Kesharwani
,
N.
Sylvetsky
, and
J. M. L.
Martin
, “
Conventional and explicitly correlated ab initio benchmark study on water clusters: Revision of the BEGDB and WATER27 data sets
,”
J. Chem. Theory Comput.
13
,
3136
3152
(
2017
).
20.
J. P.
Heindel
,
K. M.
Herman
,
E.
Aprà
, and
S. S.
Xantheas
, “
Guest–host interactions in clathrate hydrates: Benchmark MP2 and CCSD(T)/CBS binding energies of CH4, CO2, and H2S in (H2O)20 cages
,”
J. Phys. Chem. Lett.
12
,
7574
7582
(
2021
).
21.
R. K.
Nesbet
, “
Atomic Bethe-Goldstone equations
,” in , edited by
R.
LeFebvre
and
C.
Moser
(
John Wiley & Sons
,
1969
), Vol. 14, pp.
1
34
.
22.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
van der Avoird
, “
Predictions of the properties of water from first principles
,”
Science
315
,
1249
1252
(
2007
).
23.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
van der Avoird
, “
Polarizable interaction potential for water from coupled cluster calculations. I. Analysis of dimer potential energy surface
,”
J. Chem. Phys.
128
,
094313
(
2008
).
24.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
van der Avoird
, “
Polarizable interaction potential for water from coupled cluster calculations. II. Applications to dimer spectra, virial coefficients, and simulations of liquid water
,”
J. Chem. Phys.
128
,
094314
(
2008
).
25.
X.
Huang
,
B. J.
Braams
, and
J. M.
Bowman
, “
Ab initio potential energy and dipole moment surfaces of (H2O)2
,”
J. Phys. Chem. A
110
,
445
451
(
2006
).
26.
Y.
Wang
,
B. C.
Shepler
,
B. J.
Braams
, and
J. M.
Bowman
, “
Full-dimensional, ab initio potential energy and dipole moment surfaces for water
,”
J. Chem. Phys.
131
,
054511
(
2009
).
27.
Y.
Wang
,
X.
Huang
,
B. C.
Shepler
,
B. J.
Braams
, and
J. M.
Bowman
, “
Flexible, ab initio potential, and dipole moment surfaces for water. I. Tests and applications for clusters up to the 22-mer
,”
J. Chem. Phys.
134
,
094509
(
2011
).
28.
Y.
Wang
and
J. M.
Bowman
, “
Ab initio potential and dipole moment surfaces for water. II. Local-monomer calculations of the infrared spectra of water clusters
,”
J. Chem. Phys.
134
,
154510
(
2011
).
29.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
, “
Toward a universal water model: First principles simulations from the dimer to the liquid phase
,”
J. Phys. Chem. Lett.
3
,
3765
3769
(
2012
).
30.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
, “
Development of a ‘first principles’ water potential with flexible monomers: Dimer potential energy surface, VRT spectrum, and second virial coefficient
,”
J. Chem. Theory Comput.
9
,
5395
5403
(
2013
).
31.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
, “
Development of a ‘first principles’ water potential with flexible monomers. II. Trimer potential energy surface, third virial coefficient, and small clusters
,”
J. Chem. Phys.
10
, 1599–1607 (
2014
).
32.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
, “
Development of a ‘first-principles’ water potential with flexible monomers. III. Liquid phase properties
,”
J. Chem. Theory Comput.
10
,
2906
2910
(
2014
).
33.
P.
Bajaj
,
A. W.
Götz
, and
F.
Paesani
, “
Toward chemical accuracy in the description of ion–water interactions through many-body representations. I. Halide–water dimer potential energy surfaces
,”
J. Chem. Theory Comput.
12
,
2698
2705
(
2016
).
34.
M.
Riera
,
N.
Mardirossian
,
P.
Bajaj
,
A. W.
Götz
, and
F.
Paesani
, “
Toward chemical accuracy in the description of ion–water interactions through many-body representations. Alkali-water dimer potential energy surfaces
,”
J. Chem. Phys.
147
,
161715
(
2017
).
35.
V. W. D.
Cruzeiro
,
E.
Lambros
,
M.
Riera
,
R.
Roy
,
F.
Paesani
, and
A. W.
Götz
, “
Highly accurate many-body potentials for simulations of N2O5 in water: Benchmarks, development, and validation
,”
J. Chem. Theory Comput.
17
,
3931
3945
(
2021
).
36.
M.
Riera
,
E. P.
Yeh
, and
F.
Paesani
, “
Data-driven many-body models for molecular fluids: CO2/H2O mixtures as a case study
,”
J. Chem. Theory Comput.
16
,
2246
2257
(
2020
).
37.
M.
Riera
,
A.
Hirales
,
R.
Ghosh
, and
F.
Paesani
, “
Data-driven many-body models with chemical accuracy for CH4/H2O mixtures
,”
J. Phys. Chem. A
124
,
11207
11221
(
2020
).
38.
S.
Yue
,
M.
Riera
,
R.
Ghosh
,
A. Z.
Panagiotopoulos
, and
F.
Paesani
, “
Transferability of data-driven, many-body models for CO2 simulations in the vapor and liquid phases
,”
J. Chem. Phys.
156
,
104503
(
2022
).
39.
F.
Paesani
, “
Getting the right answers for the right reasons: Toward predictive molecular simulations of water with many-body potential energy functions
,”
Acc. Chem. Res.
49
,
1844
1851
(
2016
).
40.
F.
Noé
,
A.
Tkatchenko
,
K.-R.
Müller
, and
C.
Clementi
, “
Machine learning for molecular simulation
,”
Annu. Rev. Phys. Chem.
71
,
361
390
(
2020
).
41.
J.
Behler
and
M.
Parrinello
, “
Generalized neural-network representation of high-dimensional potential-energy surfaces
,”
Phys. Rev. Lett.
98
,
146401
(
2007
).
42.
J.
Behler
, “
Neural network potential-energy surfaces in chemistry: A tool for large-scale simulations
,”
Phys. Chem. Chem. Phys.
13
,
17930
17955
(
2011
).
43.
J.
Behler
, “
Representing potential energy surfaces by high-dimensional neural network potentials
,”
J. Phys.: Condens. Matter
26
,
183001
(
2014
).
44.
K.
Hansen
,
F.
Biegler
,
R.
Ramakrishnan
,
W.
Pronobis
,
O. A.
Von Lilienfeld
,
K.-R.
Müller
, and
A.
Tkatchenko
, “
Machine learning predictions of molecular properties: Accurate many-body potentials and nonlocality in chemical space
,”
J. Phys. Chem. Lett.
6
,
2326
2331
(
2015
).
45.
J.
Behler
, “
First principles neural network potentials for reactive simulations of large molecular and condensed systems
,”
Angew. Chem., Int. Ed.
56
,
12828
12840
(
2017
).
46.
S.
Chmiela
,
A.
Tkatchenko
,
H. E.
Sauceda
,
I.
Poltavsky
,
K. T.
Schütt
, and
K.-R.
Müller
, “
Machine learning of accurate energy-conserving molecular force fields
,”
Sci. Adv.
3
,
e1603015
(
2017
).
47.
J. S.
Smith
,
O.
Isayev
, and
A. E.
Roitberg
, “
ANI-1: An extensible neural network potential with DFT accuracy at force field computational cost
,”
Chem. Sci.
8
,
3192
3203
(
2017
).
48.
J. S.
Smith
,
O.
Isayev
, and
A. E.
Roitberg
, “
ANI-1, a data set of 20 million calculated off-equilibrium conformations for organic molecules
,”
Sci. Data
4
,
170193
(
2017
).
49.
L.
Zhang
,
J.
Han
,
H.
Wang
,
R.
Car
, and
W.
E
, “
Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics
,”
Phys. Rev. Lett.
120
,
143001
(
2018
).
50.
B.
Cheng
,
E. A.
Engel
,
J.
Behler
,
C.
Dellago
, and
M.
Ceriotti
, “
Ab initio thermodynamics of liquid and solid water
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
1110
1115
(
2019
).
51.
C.
Schran
,
J.
Behler
, and
D.
Marx
, “
Automated fitting of neural network potentials at coupled cluster accuracy: Protonated water clusters as testing ground
,”
J. Chem. Theory Comput.
16
,
88
99
(
2019
).
52.
M.
Haghighatlari
and
J.
Hachmann
, “
Advances of machine learning in molecular modeling and simulation
,”
Curr. Opin. Chem. Eng.
23
,
51
57
(
2019
).
53.
O.
Wohlfahrt
,
C.
Dellago
, and
M.
Sega
, “
Ab initio structure and thermodynamics of the RPBE-D3 water/vapor interface by neural-network molecular dynamics
,”
J. Chem. Phys.
153
,
144710
(
2020
).
54.
C.
Schran
,
K.
Brezina
, and
O.
Marsalek
, “
Committee neural network potentials control generalization errors and enable active learning
,”
J. Chem. Phys.
153
,
104105
(
2020
).
55.
S.
Manzhos
and
T.
Carrington
, Jr.
, “
Neural network potential energy surfaces for small molecules and reactions
,”
Chem. Rev.
121
,
10187
10217
(
2020
).
56.
P.
Gkeka
,
G.
Stoltz
,
A.
Barati Farimani
,
Z.
Belkacemi
,
M.
Ceriotti
,
J. D.
Chodera
,
A. R.
Dinner
,
A. L.
Ferguson
,
J.-B.
Maillet
,
H.
Minoux
,
C.
Peter
,
F.
Pietrucci
,
A.
Silveira
,
A.
Tkatchenko
,
Z.
Trstanova
,
R.
Wiewiora
, and
T.
Lelièvre
, “
Machine learning force fields and coarse-grained variables in molecular dynamics: Application to materials and biological systems
,”
J. Chem. Theory Comput.
16
,
4757
4775
(
2020
).
57.
Z. L.
Glick
,
D. P.
Metcalf
,
A.
Koutsoukas
,
S. A.
Spronk
,
D. L.
Cheney
, and
C. D.
Sherrill
, “
AP-Net: An atomic-pairwise neural network for smooth and transferable interaction potentials
,”
J. Chem. Phys.
153
,
044112
(
2020
).
58.
C.
Schran
,
F. L.
Thiemann
,
P.
Rowe
,
E. A.
Müller
,
O.
Marsalek
, and
A.
Michaelides
, “
Machine learning potentials for complex aqueous systems made simple
,”
Proc. Natl. Acad. Sci. U. S. A.
118
,
e2110077118
(
2021
).
59.
T. W.
Ko
,
J. A.
Finkler
,
S.
Goedecker
, and
J.
Behler
, “
A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer
,”
Nat. Commun.
12
,
398
(
2021
).
60.
J.
Behler
, “
Four generations of high-dimensional neural network potentials
,”
Chem. Rev.
121
,
10037
10072
(
2021
).
61.
O. T.
Unke
,
S.
Chmiela
,
H. E.
Sauceda
,
M.
Gastegger
,
I.
Poltavsky
,
K. T.
Schütt
,
A.
Tkatchenko
, and
K.-R.
Müller
, “
Machine learning force fields
,”
Chem. Rev.
121
,
10142
10186
(
2021
).
62.
T.
Zubatiuk
and
O.
Isayev
, “
Development of multimodal machine learning potentials: Toward a physics-aware artificial intelligence
,”
Acc. Chem. Res.
54
,
1575
1585
(
2021
).
63.
V.
Zaverkin
,
D.
Holzmüller
,
R.
Schuldt
, and
J.
Kästner
, “
Predicting properties of periodic systems from cluster data: A case study of liquid water
,”
J. Chem. Phys.
156
,
114103
(
2022
).
64.
L. D.
Jacobson
,
J. M.
Stevenson
,
F.
Ramezanghorbani
,
D.
Ghoreishi
,
K.
Leswing
,
E. D.
Harder
, and
R.
Abel
, “
Transferable neural network potential energy surfaces for closed-shell organic molecules: Extension to ions
,”
J. Chem. Theory Comput.
18
,
2354
2366
(
2022
).
65.
S.
Batzner
,
A.
Musaelian
,
L.
Sun
,
M.
Geiger
,
J. P.
Mailoa
,
M.
Kornbluth
,
N.
Molinari
,
T. E.
Smidt
, and
B.
Kozinsky
, “
E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials
,”
Nat. Commun.
13
,
2453
(
2022
).
66.
J.
Han
,
L.
Zhang
, and
R.
Car
, “
Deep potential: A general representation of a many-body potential energy surface
,”
Commun. Comput. Phys.
23
,
629
639
(
2018
).
67.
J.
Rezac
and
P.
Hobza
, “
Benchmark calculations of interaction energies in noncovalent complexes and their applications
,”
Chem. Rev.
116
,
5038
5071
(
2016
).
68.
F.
Paesani
, “
Water: Many-body potential from first principles (from the gas to the liquid phase)
,” in
Handbook of Materials Modeling: Methods: Theory and Modeling
, edited by
W.
Andreoni
and
S.
Yip
(
Springer
,
2020
), pp.
635
660
.
69.
B. B.
Bizzarro
,
C. K.
Egan
, and
F.
Paesani
, “
Nature of halide–water interactions: Insights from many-body representations and density functional theory
,”
J. Chem. Theory Comput.
15
,
2983
2995
(
2019
).
70.
C. K.
Egan
,
B. B.
Bizzarro
,
M.
Riera
, and
F.
Paesani
, “
Nature of alkali ion–water interactions: Insights from many-body representations and density functional theory. II
,”
J. Chem. Theory Comput.
16
,
3055
3072
(
2020
).
71.
L.
Zhang
,
J.
Han
,
H.
Wang
,
W.
Saidi
,
R.
Car
, and
W.
E
, “
End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems
,”
Adv. Neural Inf. Process. Syst.
31
,
4436
4446
(
2018
).
72.
S.
Yue
,
M. C.
Muniz
,
M. F.
Calegari Andrade
,
L.
Zhang
,
R.
Car
, and
A. Z.
Panagiotopoulos
, “
When do short-range atomistic machine-learning models fall short?
,”
J. Chem. Phys.
154
,
034111
(
2021
).
73.
R. G.
Parr
, “
Density functional theory of atoms and molecules
,” in
Horizons of Quantum Chemistry
(
Springer
,
1980
), pp.
5
15
.
74.
S.
Vuckovic
,
S.
Song
,
J.
Kozlowski
,
E.
Sim
, and
K.
Burke
, “
Density functional analysis: The theory of density-corrected DFT
,”
J. Chem. Theory Comput.
15
,
6636
6646
(
2019
).
75.
S.
Dasgupta
,
E.
Lambros
,
J. P.
Perdew
, and
F.
Paesani
, “
Elevating density functional theory to chemical accuracy for water simulations through a density-corrected many-body formalism
,”
Nat. Commun.
12
,
6359
(
2021
).
76.
S.
Song
,
S.
Vuckovic
,
E.
Sim
, and
K.
Burke
, “
Density-corrected DFT explained: Questions and answers
,”
J. Chem. Theory Comput.
18
,
817
827
(
2022
).
77.
E.
Sim
,
S.
Song
,
S.
Vuckovic
, and
K.
Burke
, “
Improving results by improving densities: Density-corrected density functional theory
,”
J. Am. Chem. Soc.
144
,
6625
6639
(
2022
).
78.
E.
Palos
,
E.
Lambros
,
S.
Swee
,
J.
Hu
,
S.
Dasgupta
, and
F.
Paesani
, “
Assessing the interplay between functional-driven and density-driven errors in DFT models of water
,”
J. Chem. Theory Comput.
18
,
3410
3426
(
2022
).
79.
S.
Dasgupta
,
C.
Shahi
,
P.
Bhetwal
,
J. P.
Perdew
, and
F.
Paesani
, “
How good is the density-corrected scan functional for neutral and ionic aqueous systems, and what is so right about the Hartree–Fock density?
,”
J. Chem. Theory Comput.
18
,
4745
(
2022
).
80.
E.
Lambros
,
J.
Hu
, and
F.
Paesani
, “
Assessing the accuracy of the SCAN functional for water through a many-body analysis of the adiabatic connection formula
,”
J. Chem. Theory Comput.
17
,
3739
3749
(
2021
).
81.
J. A.
Barker
and
R. O.
Watts
, “
Structure of water; A Monte Carlo calculation
,”
Chem. Phys. Lett.
3
,
144
145
(
1969
).
82.
A.
Rahman
and
F. H.
Stillinger
, “
Molecular dynamics study of liquid water
,”
J. Chem. Phys.
55
,
3336
3359
(
1971
).
83.
P.
Gallo
,
K.
Amann-Winkel
,
C. A.
Angell
,
M. A.
Anisimov
,
F.
Caupin
,
C.
Chakravarty
,
E.
Lascaris
,
T.
Loerting
,
A. Z.
Panagiotopoulos
,
J.
Russo
,
J. A.
Sellberg
,
H. E.
Stanley
,
H.
Tanaka
,
C.
Vega
,
L.
Xu
, and
L. G. M.
Pettersson
, “
Water: A tale of two liquids
,”
Chem. Rev.
116
,
7463
7500
(
2016
).
84.
S. K.
Reddy
,
S. C.
Straight
,
P.
Bajaj
,
C.
Huy Pham
,
M.
Riera
,
D. R.
Moberg
,
M. A.
Morales
,
C.
Knight
,
A. W.
Götz
, and
F.
Paesani
, “
On the accuracy of the MB-pol many-body potential for water: Interaction energies, vibrational frequencies, and classical thermodynamic and dynamical properties from clusters to liquid water and ice
,”
J. Chem. Phys.
145
,
194504
(
2016
).
85.
J. O.
Richardson
,
C.
Pérez
,
S.
Lobsiger
,
A. A.
Reid
,
B.
Temelso
,
G. C.
Shields
,
Z.
Kisiel
,
D. J.
Wales
,
B. H.
Pate
, and
S. C.
Althorpe
, “
Concerted hydrogen-bond breaking by quantum tunneling in the water hexamer prism
,”
Science
351
,
1310
1313
(
2016
).
86.
W. T. S.
Cole
,
J. D.
Farrell
,
D. J.
Wales
, and
R. J.
Saykally
, “
Structure and torsional dynamics of the water octamer from THz laser spectroscopy near 215 μm
,”
Science
352
,
1194
1197
(
2016
).
87.
S. E.
Brown
,
A. W.
Götz
,
X.
Cheng
,
R. P.
Steele
,
V. A.
Mandelshtam
, and
F.
Paesani
, “
Monitoring water clusters ‘melt’ through vibrational spectroscopy
,”
J. Am. Chem. Soc.
139
,
7082
7088
(
2017
).
88.
S. K.
Reddy
,
D. R.
Moberg
,
S. C.
Straight
, and
F.
Paesani
, “
Temperature-dependent vibrational spectra and structure of liquid water from classical and quantum simulations with the MB-pol potential energy function
,”
J. Chem. Phys.
147
,
244504
(
2017
).
89.
G. R.
Medders
and
F.
Paesani
, “
Dissecting the molecular structure of the air/water interface from quantum simulations of the sum-frequency generation spectrum
,”
J. Am. Chem. Soc.
138
,
3912
3919
(
2016
).
90.
D. R.
Moberg
,
S. C.
Straight
, and
F.
Paesani
, “
Temperature dependence of the air/water interface revealed by polarization sensitive sum-frequency generation spectroscopy
,”
J. Phys. Chem. B
122
,
4356
4365
(
2018
).
91.
M. C.
Muniz
,
T. E.
Gartner
 III
,
M.
Riera
,
C.
Knight
,
S.
Yue
,
F.
Paesani
, and
A. Z.
Panagiotopoulos
, “
Vapor–liquid equilibrium of water with the MB-pol many-body potential
,”
J. Chem. Phys.
154
,
211103
(
2021
).
92.
C. H.
Pham
,
S. K.
Reddy
,
K.
Chen
,
C.
Knight
, and
F.
Paesani
, “
Many-body interactions in ice
,”
J. Chem. Theory Comput.
13
,
1778
1784
(
2017
).
93.
D. R.
Moberg
,
S. C.
Straight
,
C.
Knight
, and
F.
Paesani
, “
Molecular origin of the vibrational structure of ice Ih
,”
J. Phys. Chem. Lett.
8
,
2579
2583
(
2017
).
94.
D. R.
Moberg
,
P. J.
Sharp
, and
F.
Paesani
, “
Molecular-level interpretation of vibrational spectra of ordered ice phases
,”
J. Phys. Chem. B
122
,
10572
10581
(
2018
).
95.
D. R.
Moberg
,
D.
Becker
,
C. W.
Dierking
,
F.
Zurheide
,
B.
Bandow
,
U.
Buck
,
A.
Hudait
,
V.
Molinero
,
F.
Paesani
, and
T.
Zeuch
, “
The end of ice I
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
24413
24419
(
2019
).
96.
T. E.
Gartner
 III
,
K. M.
Hunter
,
E.
Lambros
,
A.
Caruso
,
M.
Riera
,
G. R.
Medders
,
A. Z.
Panagiotopoulos
,
P. G.
Debenedetti
, and
F.
Paesani
, “
Anomalies and local structure of liquid water from boiling to the supercooled regime as predicted by the many-body MB-pol model
,”
J. Phys. Chem.
13
,
3652
3658
(
2022
).
97.
R. J.
Speedy
and
C. A.
Angell
, “
Isothermal compressibility of supercooled water and evidence for a thermodynamic singularity at −45 °C
,”
J. Chem. Phys.
65
,
851
858
(
1976
).
98.
C. A.
Angell
,
W. J.
Sichina
, and
M.
Oguni
, “
Heat capacity of water at extremes of supercooling and superheating
,”
J. Phys. Chem.
86
,
998
1002
(
1982
).
99.
K. H.
Kim
,
A.
Späh
,
H.
Pathak
,
F.
Perakis
,
D.
Mariedahl
,
K.
Amann-Winkel
,
J. A.
Sellberg
,
J. H.
Lee
,
S.
Kim
,
J.
Park
,
T.
Katayama
, and
A.
Nilsson
, “
Maxima in the thermodynamic response and correlation functions of deeply supercooled water
,”
Science
358
,
1589
1593
(
2017
).
100.
H.
Pathak
,
A.
Späh
,
N.
Esmaeildoost
,
J. A.
Sellberg
,
K. H.
Kim
,
F.
Perakis
,
K.
Amann-Winkel
,
M.
Ladd-Parada
,
J.
Koliyadu
,
T. J.
Lane
,
C.
Yang
,
H. T.
Lemke
,
A. R.
Oggenfuss
,
P. J. M.
Johnson
,
Y.
Deng
,
S.
Zerdane
,
R.
Mankowsky
,
P.
Beaud
, and
A.
Nilsson
, “
Enhancement and maximum in the isobaric specific-heat capacity measurements of deeply supercooled water using ultrafast calorimetry
,”
Proc. Natl. Acad. Sci. U. S. A.
118
,
e2018379118
(
2021
).
101.
K. H.
Kim
,
K.
Amann-Winkel
,
N.
Giovambattista
,
A.
Späh
,
F.
Perakis
,
H.
Pathak
,
M. L.
Parada
,
C.
Yang
,
D.
Mariedahl
,
T.
Eklund
,
T. J.
Lane
,
S.
You
,
S.
Jeong
,
M.
Weston
,
J. H.
Lee
,
I.
Eom
,
M.
Kim
,
J.
Park
,
S. H.
Chun
,
P. H.
Poole
, and
A.
Nilsson
, “
Experimental observation of the liquid-liquid transition in bulk supercooled water under pressure
,”
Science
370
,
978
982
(
2020
).
102.
L.
Kringle
,
W. A.
Thornley
,
B. D.
Kay
, and
G. A.
Kimmel
, “
Reversible structural transformations in supercooled liquid water from 135 to 245 K
,”
Science
369
,
1490
1492
(
2020
).
103.
T. E.
Gartner
 III
,
L.
Zhang
,
P. M.
Piaggi
,
R.
Car
,
A. Z.
Panagiotopoulos
, and
P. G.
Debenedetti
, “
Signatures of a liquid–liquid transition in an ab initio deep neural network model for water
,”
Proc. Natl. Acad. Sci. U. S. A.
117
,
26040
26046
(
2020
).
104.
G. M.
Sommers
,
M. F. C.
Andrade
,
L.
Zhang
,
H.
Wang
, and
R.
Car
, “
Raman spectrum and polarizability of liquid water from deep neural networks
,”
Phys. Chem. Chem. Phys.
22
,
10592
10602
(
2020
).
105.
L.
Zhang
,
H.
Wang
,
R.
Car
, and
W.
E
, “
Phase diagram of a deep potential water model
,”
Phys. Rev. Lett.
126
,
236001
(
2021
).
106.
C.
Zhang
,
F.
Tang
,
M.
Chen
,
J.
Xu
,
L.
Zhang
,
D. Y.
Qiu
,
J. P.
Perdew
,
M. L.
Klein
, and
X.
Wu
, “
Modeling liquid water by climbing up Jacob’s ladder in density functional theory facilitated by using deep neural network potentials
,”
J. Phys. Chem. B
125
,
11444
11456
(
2021
).
107.
W.
Liang
,
G.
Lu
, and
J.
Yu
, “
Machine-learning-driven simulations on microstructure and thermophysical properties of MgCl2–KCl eutectic
,”
ACS Appl. Mater. Interfaces
13
,
4034
4042
(
2021
).
108.
T.
Wen
,
R.
Wang
,
L.
Zhu
,
L.
Zhang
,
H.
Wang
,
D. J.
Srolovitz
, and
Z.
Wu
, “
Specialising neural network potentials for accurate properties and application to the mechanical response of titanium
,”
npj Comput. Mater.
7
,
206
(
2021
).
109.
D.
Lu
,
H.
Wang
,
M.
Chen
,
L.
Lin
,
R.
Car
,
W.
E
,
W.
Jia
, and
L.
Zhang
, “
86 PFLOPS deep potential molecular dynamics simulation of 100 million atoms with ab initio accuracy
,”
Comput. Phys. Commun.
259
,
107624
(
2021
).
110.
H.
Niu
,
L.
Bonati
,
P. M.
Piaggi
, and
M.
Parrinello
, “
Ab initio phase diagram and nucleation of gallium
,”
Nat. Commun.
11
,
2654
(
2020
).
111.
Z.
Guo
,
D.
Lu
,
Y.
Yan
,
S.
Hu
,
R.
Liu
,
G.
Tan
,
N.
Sun
,
W.
Jiang
,
L.
Liu
,
Y.
Chen
 et al, “
Extending the limit of molecular dynamics with ab initio accuracy to 10 billion atoms
,” in
Proceedings of the 27th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming
(
Association for Computing Machinery, New York,
,
2022
), pp.
205
218
.
112.
H.
Partridge
and
D. W.
Schwenke
, “
The determination of an accurate isotope dependent potential energy surface for water from extensive ab initio calculations and experimental data
,”
J. Chem. Phys.
106
,
4618
4639
(
1997
).
113.
B. J.
Braams
and
J. M.
Bowman
, “
Permutationally invariant potential energy surfaces in high dimensionality
,”
Int. Rev. Phys. Chem.
28
,
577
606
(
2009
).
114.
H.
Wang
,
L.
Zhang
,
J.
Han
, and
W.
E
, “
DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics
,”
Comput. Phys. Commun.
228
,
178
184
(
2018
).
115.
Y.
Zhang
,
H.
Wang
,
W.
Chen
,
J.
Zeng
,
L.
Zhang
,
H.
Wang
,
W.
E
, and
D. P.
Gen
, “
A concurrent learning platform for the generation of reliable deep learning based potential energy models
,”
Comput. Phys. Commun.
253
,
107206
(
2020
).
116.
Paesani Research Group
, MBX: A many-body energy and force calculator, available at https://github.com/paesanilab/MBX (
2022
).
117.
M.
Pinheiro
,
F.
Ge
,
N.
Ferré
,
P. O.
Dral
, and
M.
Barbatti
, “
Choosing the right molecular machine learning potential
,”
Chem. Sci.
12
,
14396
14413
(
2021
).
118.
W.
Shinoda
,
M.
Shiga
, and
M.
Mikami
, “
Rapid estimation of elastic constants by molecular dynamics simulation under constant stress
,”
Phys. Rev. B
69
,
134103
(
2004
).
119.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in ’t Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
, and
S. J.
Plimpton
, “
LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
120.
J. R.
Errington
and
P. G.
Debenedetti
, “
Relationship between structural order and the anomalies of liquid water
,”
Nature
409
,
318
321
(
2001
).
121.
P. M.
Piaggi
,
A. Z.
Panagiotopoulos
,
P. G.
Debenedetti
, and
R.
Car
, “
Phase equilibrium of water with hexagonal and cubic ice using the scan functional
,”
J. Chem. Theory Comput.
17
,
3065
3077
(
2021
).
122.
G. S.
Kell
, “
Isothermal compressibility of liquid water at 1 atm
,”
J. Chem. Eng. Data
15
,
119
122
(
1970
).
123.
J.
Zinn-Justin
, “
Precise determination of critical exponents and equation of state by field theory methods
,”
Phys. Rep.
344
,
159
178
(
2001
).
124.
P. J.
Linstrom
and
W. G.
Mallard
, “
The NIST chemistry WebBook: A chemical data resource on the internet
,”
J. Chem. Eng. Data
46
,
1059
1063
(
2001
).
125.
G. R.
Medders
,
A. W.
Götz
,
M. A.
Morales
,
P.
Bajaj
, and
F.
Paesani
, “
On the representation of many-body interactions in water
,”
J. Chem. Phys.
143
,
104102
(
2015
).
126.
M.
Riera
,
E.
Lambros
,
T. T.
Nguyen
,
A. W.
Götz
, and
F.
Paesani
, “
Low-order many-body interactions determine the local structure of liquid water
,”
Chem. Sci.
10
,
8211
8218
(
2019
).
127.
E.
Lambros
and
F.
Paesani
, “
How good are polarizable and flexible models for water: Insights from a many-body perspective
,”
J. Chem. Phys.
153
,
060901
(
2020
).
128.
P.
Schienbein
and
D.
Marx
, “
Liquid–vapor phase diagram of RPBE-D3 water: Electronic properties along the coexistence curve and in the supercritical phase
,”
J. Phys. Chem. B
122
,
3318
3329
(
2017
).
129.
S. L.
Bore
and
F.
Paesani
, “
Quantum phase diagram of water
,” chemRxiv, (
2023
).

Supplementary Material