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

Fast, reliable, and accurate modeling of structural, physical, and chemical properties of water across all media - gas, liquid, interface, confined, and solid - and at different thermodynamic conditions is a long-standing challenge. This is largely because the quality of any theoretical prediction depends heavily on the underlying intermolecular potential energy surface (PES) utilized.^1–4 While in principle, ab initio molecular dynamics (AIMD) simulations carried out with correlated electronic structure methods can provide a correct description of water across all phases, these simulations are currently unfeasible due to the prohibitive associated computational cost. Despite recent progress in the development and implementation of simulation approaches based on the Möller-Plesset perturbation theory, density functional theory (DFT) still remains the most common approach used in ab initio simulations of water.^5–8 However, existing functionals have been shown to remain limited in the accuracy and predictive ability with which they can represent the properties of water in different phases.⁹ 

On the other hand, empirical water models based on molecular mechanics (commonly referred to as force fields) have been used extensively in computer simulations to investigate the properties of water under different thermodynamic conditions. In these models, the force field parameters are usually tuned to reproduce a (limited) set of experimental properties.^4,10–17 The most common empirical models are pairwise additive, assume the water molecules to be rigid, and use fixed point charges to describe the electrostatic interactions. Despite their simplicity, empirical models have had great success in reproducing, at least to some extent, structural, thermodynamic, and transport properties of liquid water over a broad range of temperatures and pressures.^4,10,12–17 Examples under this category include RWK,¹⁸ SPC,¹⁹ SPC/E,¹⁶ TIP4P,¹⁴ TIP4P-2005,¹² TIP4P-Ew,¹³ and TIP5P.¹⁵ The reader is referred to Refs. 4, 10, and 11 as well as the original studies for complete details of rigid water models. Flexible versions were also developed for some of these water models to investigate vibrational dynamics and energy transfer in the liquid phase.^20,21 To go beyond the pairwise approximation, several water models have been developed, which include either empirical three-body terms (e.g., E3B^22–24) or implicit many-body effects through classical polarization (e.g., BK3,²⁵ SWM4-DP,²⁶ SWM4-NDP,²⁷ SWM6,²⁸ COS,^25,29,30 TTMx-F,^31–33 AMOEBA,^34–37 GEM*,³⁸ and POLY2VS³⁹). The interested reader is referred to Ref. 1 for a systematic overview of recent water models.

While water models parameterized to reproduce experimental data have been instrumental in gaining qualitative insights into the behavior of liquid water, the lack of chemical accuracy results in limited predictive power across the entire phase diagram. Over the years, this has stimulated the development of analytical potential energy functions, originally introduced by Clementi and co-workers in the 1970s, which aim at representing the multidimensional potential energy surface (PES) associated with a system containing N water molecules through a rigorous representation of the many-body expansion (MBE) of the interaction energy,⁴⁰ 

where $V1B$ is the one-body (1B) contribution corresponding to the deformation energy and the $VnB$ are the n-body (nB) interaction energies defined as

Depending on the specific functional form adopted and the extent of the training sets used to fit the individual terms of the MBE, many-body potentials can approach the same level of accuracy as high-level correlated electronic structure calculations at a fraction of the computational cost.⁴¹ The most notable many-body potential energy functions are CC-pol,^2,3,42–45 WHBB,^46–50 HBB2-pol,^41,51,52 and MB-pol.^53–56 To date, MB-pol (and its precursor HBB2-pol) is the only many-body potential that has been consistently employed in molecular simulations, with explicit inclusion of nuclear quantum effects (NQE), which correctly reproduce the properties of water from the gas to the liquid and solid phases.^51–60 

The purpose of this study is to systematically assess the accuracy of MB-pol in predicting structural, thermodynamic, dynamical, and spectroscopic properties of water across all phases as well as to provide a metric by which these properties can be compared to experiment. The article is organized as follows. Section II describes the technical details of all calculations and simulations employed in this study. Section III reports the results on several physical properties of water in the three phases: gas, liquid, and solid. The overall performance of MB-pol is then discussed by assigning a score to each computed property. Section IV summarizes the results and provides an outlook of future applications of MB-pol.

As discussed in detail in previous studies,^41,53–56 the MB-pol potential is built upon the MBE of Eq. 1 and includes explicit terms that describe 1B, 2B, and 3B terms, along with classical N-body polarization to account for all higher-body contributions to the interaction energy. The polarization term is represented by a slightly modified version of the Thole-type model as introduced in the TTM4-F model.³³ The 1B term is represented by the spectroscopically accurate water monomer PES developed by Partridge and Schwenke.⁶¹ The 2B term is further divided into long-range and short-range interactions and is described using classical electrostatics, induction, and dispersion forces, which dominate the long-range part, supplemented by a set of multivariable polynomials in the short-range, to capture the more complex quantum mechanical effects arising from the overlap of monomer electron densities. Along the same lines, the 3B term is composed of 3B induction and a multi-dimensional function to accurately describe both long-range and short-range interactions. The polynomial functions used to describe the short-range 2B and 3B interactions are generated in such a way that they retain permutational invariance with respect to the hydrogen atoms within the same water molecule as well as to whole water molecules within all possible dimers and trimers. The permutationally invariant polynomials are trained to a large set of highly accurate correlated coupled-cluster energies via a supervised machine learning approach. Specific details about the functional form, training sets, and training algorithms can be found in the original studies.^{41,51,53–56} MB-pol is publicly available as an external plugin⁶² for the OpenMM toolkit for molecular simulations⁶³ and has recently been interfaced⁶⁴ to the i-PI wrapper.⁶⁵ As reported in Refs. 1, 57, and 66, MB-pol is ∼50 times computationally more expensive than flexible, point-charge models (e.g., q-TIP4P/F²⁰ and q-SPC/Fw²¹) and ∼3-7 times (depending on the induction scheme) more expensive than polarizable models (e.g., AMOEBA^34–37 and TTMx-F^31–33).

All electronic structure calculations of water clusters presented in Secs. III and IV were carried out using MOLPRO.⁶⁷ The reference interaction energies for (H₂O)_n clusters, with n = 4–6, optimized at the RI-MP2 and MP2 levels of theory in Refs. 68 and 69, were obtained using the MBE of the interaction energy⁴⁰ as described in Ref. 70, with individual MB contributions calculated with the CCSD(T)/CCSD(T)-F12b method in the complete basis set (CBS) limit.^71,72 The 2B interaction energies were computed at the CCSD(T) level of theory by extrapolating the values obtained with aug-cc-pVTZ and aug-cc-pVQZ basis sets supplemented with an additional set of (3s, 3p, 2d, 1f) midbond functions, with exponents equal to (0.9, 0.3, 0.1) for s and p orbitals, (0.6, 0.2) for d orbitals, and 0.3 for f orbitals, placed at the center of mass of each dimer configuration.^73–75 The following two-point formula^76,77 was used to extrapolate the interaction energies to the CBS limit:

with cardinal numbers X = 3 and 4, accordingly. The Hartree–Fock energy was not extrapolated separately since it was close to the CBS limit for either value of X. The 3B interaction energies were calculated at the CCSD(T) level of theory using the aug-cc-pVTZ basis set supplemented with the same set of midbond functions introduced above, which were placed at the center of mass of each trimer configuration. All higher (>3B) contributions were computed with the CCSD(T)-F12 method using the VTZ-F12 basis sets.^78–80 This method yields results close to the CBS values at lower computational cost than direct CCSD(T) calculations with large basis sets.^81,82 All 3B and higher-body energies were corrected for the basis set superposition error (BSSE) using the counterpoise method.⁸³ 

The MB-pol vibrational frequencies of the water clusters were calculated within the harmonic approximation from the diagonalization of the mass-weighted Hessian matrix. Each cluster structure was first energy minimized until the norm of the force vector reached a value smaller than 10⁻⁸ kcal mol⁻¹ Å⁻¹. The absence of any imaginary frequencies indicates that all structures reported in Section III A 2 correspond to either a local or a global minimum.

All molecular dynamics (MD) simulations of liquid water presented in Secs. III and IV were carried out at the classical level using in-house software based on the DL_POLY_2 simulation package,⁸⁴ which was modified to include the MB-pol potential.⁵⁵ Unless otherwise stated, the systems consisted of 256 water molecules placed in a cubic simulation box. The velocity-Verlet algorithm was used to integrate Newton’s equations of motion with a time step of 0.2 fs. All thermodynamic properties except the surface tension were calculated from simulations carried out in the isobaric-isothermal (NPT) ensemble. The temperature and pressure were kept constant using the Nosé-Hoover thermostat and barostat, respectively.^85–88 The short-range interactions were evaluated within a distance cutoff of 9 Å. Short-range electrostatic interactions were computed in real space using Coulomb’s law while the long range interactions were calculated in reciprocal space using the Ewald summation technique. A high precision was used (10⁻⁸) to generate the best Ewald vectors and the Ewald convergence parameter for MD simulations. The MD simulations were run in the NPT ensemble at atmospheric pressure (P = 1 atm) for twelve temperatures between 248 K and 368 K. The trajectory lengths at different temperatures are listed in Table XL of the supplementary material. During the production run, the positions and dipole moments were collected for analysis. To compute the isothermal compressibility, additional simulations were run for 4.5 ns at each temperature for the following pressures: −1.5, −1.0, −0.5, 0.5, 1.0, and 1.5 katm.

To calculate the surface tension of liquid water, a slab geometry was prepared from a fully equilibrated bulk simulation of 511 water molecules in a cubic box by expanding the z-axis of the box to 100 Å. After preparation of the initial slab geometry, equilibration trajectories of 500 ps were simulated in the isochoric-isothermal (NVT) ensemble by employing periodic boundary conditions and Ewald summation in all three dimensions. Following the equilibration, a 4 ns trajectory was generated at each temperature for analysis. The self-diffusion coefficient was calculated, at each temperature, by averaging over thirty independent 100 ps long trajectories carried out in the microcanonical (NVE) ensemble. The initial configurations for the NVE trajectories were obtained from thirty 10 ps long NVT trajectories started from configurations extracted at intervals of 50 ps from an equilibrated NPT trajectory.

For each ice phase, the system size was chosen such that the edges of the simulation box were always separated by more than 18 Å. The ice phases considered in this study are ice I_h, the proton disordered phase at ambient conditions, and several proton ordered phases (II, VIII, IX, XIII, XIV, and XV). The initial structures for the proton ordered phases were taken from Ref. 89, while for ice I_h, 13 independent configurations were generated by minimizing the net dipole moment following the algorithm proposed in Ref. 90. One additional configuration for ice I_h taken from Ref. 89 was also included in the calculations. All configurations satisfy the Bernal-Fowler rules.⁹¹ The densities of the different ice phases were calculated by averaging over 100 ps of NPT simulations.

All (H₂O)_n clusters, with n = 4–6, considered in the analysis of many-body contribution to the interaction energies are shown in Figure 1, while the errors associated with the individual terms are shown in Figure 2. As explained in Section II, the reference data were obtained with the CCSD(T)/CCSD(T)-F12b method in the CBS limit. Independent of the cluster size and geometry, MB-pol exhibits small errors, which are always within 0.5 kcal mol⁻¹ relative to the reference values, for all terms of the MBE of Eq. (1). The error in the total interaction energy increases with system size as the individual errors start to add up, most prevalently in ring-type configurations that consist of repeating dimer and trimer units. Due to extended hydrogen bonding and symmetry, the ring-type isomers also show larger higher-body contributions that can be non-negligible.⁵⁷ The comparison between the reference and MB-pol interaction energies for the tetramer, pentamer, and hexamer isomers shown in Figure 3 indicates that the relative interaction energy order for the different isomers of each cluster is retained by MB-pol, with a maximum deviation in the computed interaction energies of 0.84 kcal mol⁻¹. In this analysis, the interaction energy is defined as the cluster energy minus the energy of the individual water molecules kept at the same geometry as in the cluster. The comparison in Figure 3 thus directly probes the actual interaction between water molecules without being affected by differences in the CCSD(T)/CCSD(T)-F12b and MB-pol representations of the monomer distortion energies.

Scoring. To quantify the accuracy of MB-pol in describing many-body interactions, scores are assigned to each (H₂O)_n cluster based on the following three quantities: first, the maximum of the total unsigned error of the individual terms of the MBE, $ΔEunsigned,maxMB$⁠, across all isomers of an n-mer,

second, the maximum absolute error in the total interaction energy, $ΔEInt,max$⁠, across all isomers of an n-mer,

third, the relative energy difference between the interaction energy of the prism and cage hexamer isomers,

A small value for $ΔEInt,max$ with a larger error in $ΔEunsigned,maxMB$ would indicate that the MBE benefits from error compensation.

For the first two criteria, 100 points are assigned if the magnitude of the error is less than 1 kcal mol⁻¹, which is commonly defined as chemical accuracy. 10 points are deducted successively for every 0.5 kcal mol⁻¹ increase in the error. For the third criterion, the reference value calculated using the CCSD(T)/CCSD(T)-F12/CBS method is −0.29 kcal mol⁻¹. A score of zero is assigned for the wrong sign while a score of 100 is assigned if the relative energy is within 0.1 kcal mol⁻¹ from the reference value. 10 points are deducted for every additional 0.1 kcal mol⁻¹ difference. The MB-pol values for the three criteria and associated scores are reported in Table I, which demonstrates that MB-pol receives perfect scores in all three criteria for (H₂O)_n clusters with n = 4–6.

The ability of MB-pol in predicting accurate vibrational frequencies is assessed through the analysis of harmonic frequencies calculated for small water clusters, from the water, monomer to the hexamer. The comparisons are made with recently published benchmark data that are expected to closely approximate CCSD(T)/CBS values.^92,93

The water clusters included in this analysis are all structures with lowest energy for n = 2–6, along with the cage, book1, and cyclic-chair (also referred to as ring) isomers of the hexamer (see Figure 1). For n =3–5, the minimum energy structures correspond to cyclic structures in which each water molecule donates and accepts one hydrogen bond. In the water hexamer the lowest energy structure is the prism isomer, which is nearly isoenergetic with the cage isomer (see Section III A 1), while the book1 and cyclic-chair isomers are higher in energy. Prism and cage are the predominant isomers at very low temperatures while the cage and cyclic-chair isomers appear at higher temperatures.^52,94–99

Table II reports the average deviation (AD), average absolute deviation (AAD), and the maximum absolute deviation (MAD) between the reference and the MB-pol harmonic frequencies calculated for the water monomer and each of the eight water clusters. The complete list of harmonic frequencies can be found in the supplementary material. In all cases, MB-pol accurately reproduces the reference data, with AAD and MAD values never exceeding 17 cm⁻¹ and 50 cm⁻¹, respectively. The MB-pol deviations from the reference data increase with increasing cluster size, in particular, for cyclic structures, in which, as discussed in Section III A 1, small errors in the MBE can add up for repeating units due to the inherent symmetry. As a result, the deviations for the hexamer book1, cage, and prism isomers are smaller than for the cyclic-chair isomer. On average, the MB-pol harmonic frequencies are slightly redshifted as compared to the reference data. This redshift is mostly due to the low-frequency intermonomer and bending modes, while the OH stretch frequencies are slightly blueshifted (see the list of frequencies in the supplementary material). Both red and blue shifts associated, respectively, with the bending and stretching vibrations can be explained by considering that MB-pol does not allow water autoionization, and, consequently, underestimates the ionic character, and thus the strength, of the hydrogen bonds. Importantly, it has been shown that MB-pol provides a consistent representation of the vibrational structure of water independently of the system size, correctly predicting infrared,⁵⁶ Raman,⁵⁶ sum-frequency generation,⁵⁹ and two-dimensional infrared⁶⁰ spectra of liquid water at ambient conditions.

Scoring. To quantify the accuracy with which MB-pol reproduces the reference harmonic frequencies, a score for each cluster is assigned based on the corresponding MAD value. A score of 100 is assigned if the MAD is below the threshold value of 20 cm⁻¹, and 10 points are then deducted for each additional increment of 20 cm⁻¹. The MB-pol scores are reported in Table II.

The accuracy of MB-pol in reproducing the properties of liquid water is assessed through comparisons with the corresponding experimental data as a function of temperature. As discussed in detail in Section II, all MB-pol properties were calculated from classical MD simulations. The role played by NQE (e.g., zero-point energy and tunneling) will be discussed in a forthcoming publication. A summary of several thermodynamic properties of liquid water computed with MB-pol at atmospheric pressure (P = 1 atm) along with the corresponding experimental values is reported in the supplementary material (Tables XL and XLI).

The temperature dependence of the density of liquid water at 1 atm calculated from classical MD simulations with MB-pol is shown in Figure 4. The experimental data are taken from Refs. 100 and 101. At high temperature, the MB-pol results are in excellent agreement with the corresponding experimental values. As the temperature decreases, the difference between MB-pol and experiment increases nearly linearly. The maximum and average absolute deviations from the reference values are 0.042 and 0.013 g cm⁻³, respectively. The temperature of maximum density, obtained by calculating the analytical derivative of a fifth-order polynomial interpolating the simulation results, is located at 263 K, which is 14 K lower than experiment.^100,101 However, it should be noted that, based on the experimental data for D₂O and T₂O, the temperature of maximum density for hypothetical classical water should be ∼287 K.

The systematic deviation between MB-pol and experimental density values as the temperature decreases can be attributed, at least in part, to the neglect of NQE in the simulations, which, as expected, become increasingly important at lower temperature.¹⁰² In this context, at 298 K, the density predicted by classical MD simulations with MB-pol is 1.007 g cm⁻³ compared to the experimental value of 0.997 g cm⁻³. Using path-integral molecular dynamics (PIMD) simulations, which explicitly include NQE, it was shown that the density of water predicted by MB-pol at the same temperature decreases by 0.6%,⁵⁵ in excellent agreement with the experimental value.

The enthalpy of vaporization, $ΔHvap$⁠, is one of the properties usually included in the fitting procedures to parameterize empirical water potentials^12,13 and is defined as

where H, U, and V are the enthalpy, internal energy, and volume, respectively, and the subscripts denote that the molecules are in either the gas or liquid state. Since, at low pressure, the gas can be considered ideal, the contribution to U_gas due to interactions between water molecules can be neglected and $ΔHvap$ can then be rewritten (for 1 mol of water) as

where U_gas contains the average kinetic (i.e., 3/2 RT) and potential (i.e., intramolecular distortion) energies of the gas molecules. At each temperature, the average potential energy was calculated from a 2 ns long classical MD simulation of a single water molecule. It should be noted that, unlike point charge models, $ΔHvap$ calculated with MB-pol implicitly includes the self-polarization correction,^12,13 since MB-pol correctly describes the change in dipole moment from the gas to the condensed phase.

Figure 5 shows that $ΔHvap$ calculated from classical MD simulations with MB-pol is systematically larger than the corresponding experimental values, with an average absolute deviation of 0.41 kcal mol⁻¹. As for the density, the deviation from experiment increases as the temperature decreases, which can be attributed to the neglect of NQE in the simulations. At 298 K, $ΔHvap$ predicted by MB-pol is 10.93 kcal mol⁻¹, 0.42 kcal mol⁻¹ larger than the experimental value. Guillot and Guissani suggested that $ΔHvap$ for hypothetical “classical water” at 298.15 K should be 11.0 kcal mol⁻¹,¹⁰³ which is in excellent agreement with the value obtained from classical MD simulations with MB-pol.

The isothermal compressibility, $κT$⁠, is defined as

where V and P are the volume and the pressure of liquid water at a given temperature T, respectively. To solve Eq. (9), classical MD simulations with MB-pol were carried out at each temperature for seven pressure values (−1.5, −1.0, −0.5, 0.001, 0.5, 1.0, 1.5 katm), which were used to determine the corresponding average volumes. $κT$ was then calculated from the derivative of a third-order polynomial that was used to fit the volumes as a function of pressure at each temperature. As shown in Figure 6, the MB-pol values are in good agreement with the corresponding experimental data,¹⁰⁴ correctly predicting a minimum between 310 and 330 K. The average and maximum absolute deviations from experiment in the temperature range between the freezing and the boiling point are 0.5 × 10⁻⁶ and 1.8 × 10⁻⁶ atm⁻¹, respectively. As expected, the deviations from the experimental data become more pronounced at low temperature due to the neglect of NQE in classical MD simulations.

The thermal expansion coefficient, $αP$⁠, is defined as

where V, T, and $ρ$ are the volume, the temperature, and the density of liquid water at a given pressure P, respectively. From Eq. (10), it is evident that any water model that accurately predicts the temperature dependence of the density is also capable of correctly reproducing the variation of $αP$⁠. The MB-pol values of $αP$⁠, shown in Figure 7 along with the corresponding experimental data,¹⁰⁴ were determined numerically by solving Eq. (10) using the density values reported in Figure 4. Similar to $ρ$⁠, $αP$ derived from classical MD simulations with MB-pol is in close agreement with experiment above the freezing point, with average and maximum absolute deviations of 1.3 × 10⁻⁴ and 1.7 × 10⁻⁴ K⁻¹, respectively. At the classical level, MB-pol predicts $αP$ to be zero at 263 K, in line with the density maximum reported in Section III B 1 (experimentally, $αP=0$ at 277 K). Again,NQE appear to be important to quantitatively reproduce the variation $αP$ at low temperature, as already found in Secs. III B 1–III B 3 for other thermodynamic properties of liquid water.

The isobaric heat capacity, c_P, is defined as

where H and T are the enthalpy and the temperature of liquid water at a given pressure P, respectively. It should be noted that, unlike other thermodynamic properties, c_P is significantly affected by NQE even at room temperature.^105,106 Using Eq. (11), c_P was calculated as the temperature derivative of a fifth-order polynomial interpolating the values of H obtained at different temperatures from classical MD simulations carried out in the NPT ensemble. Consistent with previous studies, the classical MB-pol results shown in Figure 8 are larger than the corresponding experimental data^101,107 at all temperatures. Experimentally, it is known that the difference between the heat capacity of liquid H₂O and D₂O increases as the temperature decreases.¹⁰⁸ For example, the heat capacities of H₂O and D₂O at 250 K are 19.07 and 23.93 cal mol⁻¹ K⁻¹, respectively.^101,107,108 The corresponding value obtained from classical MD simulations with MB-pol is 29.25 cal mol⁻¹ K⁻¹, reinforcing the notion that explicit inclusion of NQE in simulations with MB-pol is required for quantitative, and physically correct, calculations of c_P at all temperatures. Levitt et al. estimated that 6 cal mol⁻¹ K⁻¹ must be subtracted from the classical value of the heat capacity at constant volume, c_V, to account for NQE at 298.15 K.¹⁰⁹ Assuming that the same correction can be applied to c_P, the quantum-corrected MB-pol result at 298 K becomes 21.85 cal mol⁻¹ K⁻¹, in closer agreement with the experimental value.

The surface tension, $γ$⁠, can be calculated from MD simulations as¹¹⁰ 

where the z-direction is defined along the normal vector to the surface, $α$ and $β$ refer to the liquid and the gas phase, respectively, and N and T denote the normal and tangential components of the pressure tensor, respectively. Considering the slab geometry defined in Section II, Eq. (12) reduces to

where P_xx, P_yy, and P_zz are the diagonal elements of the pressure tensor, L_z is the length of the simulation box in the z-direction, and the angular brackets denote an ensemble average. The comparison between experimental^111,112 and MB-pol values of the surface tension is shown in Figure 9. Since the MD simulations were performed using a single box with dimensions 26 Å ×  26 Å × 100 Å, it is not possible to correct the values of the surface tension for finite size effects.^113–117 However, for a simulation box with dimensions similar to those used in the present MD simulations with MB-pol, the correction for finite size effects has recently been estimated to be ∼5.7 $±$ 7.6 mJ/m².¹¹⁶ Overall, the MB-pol results are in good agreement with the experimental data over the entire temperature range considered in this study, correctly reproducing (within statistical uncertainty) the linear increase of $γ$ as the temperature decreases. At 298 K, the surface tension obtained from classical MD simulations with MB-pol is $68.4±4$ mJ m⁻² compared to the experimental value of 71.97 mJ m⁻².¹¹¹ It should be noted that, at each temperature, the agreement with the experimental data improves if the correction value estimated in Ref. 116 is added. Contrary to other thermodynamic properties, the differences between the experimental and MB-pol values of the surface tension remain effectively constant as a function of temperature, suggesting that $γ$ is less sensitive to NQE.

The static dielectric constant, $ϵ$⁠, is defined as

where M is the total dipole moment of the simulation box, k_B is Boltzmann’s constant, V and T are the volume and the temperature, respectively, and the angular brackets denote an ensemble average. The temperature dependence of the dielectric constant calculated from classical MD simulations with MB-pol is compared in Figure 10 with the corresponding experimental data.¹¹⁸ $ϵ$ calculated from the MB-pol simulations is in good agreement with experiment over the entire temperature range considered in this study. At 298 K, the value of $ϵ$ obtained from the MB-pol simulations is 68.4 which is ∼13% smaller than the experimental value of 78.5. These differences may be related, at least in part, to (small) differences in the liquid structure due to the neglect of NQE,⁵⁵ which, in turn, may affect the dipole moments of the water molecules.

The self-diffusion coefficient, D, of liquid water can be computed from the velocity autocorrelation function as

where v_i is the velocity of the center of mass of ith water molecule and the brackets indicate an ensemble average in the microcanonical (NVE) ensemble. The MB-pol results are compared in Figure 11 with the corresponding experimental values.^119–122 At 298 K, classical MD simulations with MB-pol predict D to be $0.23±0.02$ Ų ps⁻¹ in agreement with the experimental value of 0.229 Ų ps⁻¹. It should be mentioned that an incorrect value for the classical MB-pol diffusion coefficient at room temperature was previously reported by some of us in Ref. 55. The similarity between the present classical value and the corresponding quantum value obtained in Ref. 55 from centroid molecular dynamics simulations suggests the presence of competing NQE in the water diffusion as originally suggested in Ref. 20.

Since it was shown that the calculation of the self-diffusion coefficient from MD simulations is particularly sensitive to the system size,¹²³ additional NVE simulations were carried out for a system containing 512 molecules. For this larger system, $D=0.24±0.02$ Ų ps^$−$1 at room temperature. The self-diffusion coefficient in the limit of an infinite system size can be calculated as^123,124

where k_B is Boltzmann’s constant, T is the temperature, η is the viscosity, and $ξ$ is a constant that depends on the shape of the simulation box (for cubic, $ξ=2.837 297$⁠). Plugging the values obtained for 256 and 512 water molecules into Eq. (16), and using the experimental value for η, $D(∞)=0.280±0.040$ Ų ps⁻¹ at 298 K.

Scoring. As in Section III A, scores are assigned to each property of liquid water calculated with MB-pol. Following the procedure described in Refs. 4, 17, and 9 the relative error (⁠$δx$⁠) associated with each property was used to define the corresponding score as

where $δx=|xsim−xexp|/|xexp|$ and x_tol is the threshold value. The score is 100 if the $δx$ value is between 0 and x_tol. The score is thus reduced by 10% for each successive increment by x_tol. A score of 100 indicates perfect agreement with experiment while a score of 0 indicates poor performance. The x_tol value for each property was taken from Ref. 4, except for $αP$ that was not included in that study. The scores were calculated at three different temperatures. The first temperature selected to assess the performance of MB-pol is 248 K since the maximum absolute deviations for all properties calculated from classical MD simulations with MB-pol are observed at this temperature. As mentioned above, since NQE are not included in the present MB-pol simulations, it is not surprising that the computed values deviate significantly from experiment, especially at low temperatures. The performance of MB-pol is also assessed at room temperature (298 K) and at higher temperature (360 K). Polynomial fits were used to obtain the values of all properties at 360 K. Table III shows the scores obtained from classical MD simulations with MB-pol for all properties discussed in Secs. III B 1–III B 7. The average scores of MB-pol are 71, 91, and 91 at 248, 298, and 360 K, respectively.

Figure 12 shows the oxygen-oxygen (OO) radial distribution function (RDF) calculated with MB-pol at several temperatures along with the corresponding experimental curves derived from the most recent X-ray diffraction measurements.^125,126 The temperature dependence of the oxygen-hydrogen and hydrogen-hydrogen RDFs is reported in the supplementary material. The MB-pol RDFs are in good agreement with the experimental data at all temperatures, although they slightly overestimate the height of the first peak. As discussed in previous studies, this difference, which increases with decreasing temperature, can be attributed to the neglect of NQE and lead to more structural order in the classical liquid than is experimentally observed. It has been shown that, while the position of the first peak moves to larger r_OO linearly with increasing temperature, the shift in the position of the second peak deviates from a linear dependence on r_OO above ∼320 K.¹²⁵ This trend is correctly captured by classical NPT simulations with MB-pol, which reproduce the experimental data nearly quantitatively (Figure 13). The numerical comparison between MB-pol and experimental OO RDFs is based on the position of the 1st (r₁) and 2nd (r₂) peaks, and the height of the 1st peak (g(r₁)) reported in Table IV. Altogether, these quantities are used to assess the accuracy of MB-pol at 268, 278, and 308 K, through comparisons with the corresponding experimental values.^125,126

Scoring. In addition to the numerical comparison, a score was assigned to each of the three quantities using Eq. (17) with tolerance values of 0.1%, 0.5%, and 5.0% respectively for the r₁, r₂, and g(r₁). The average MB-pol scores obtained for the three quantities at at 268, 278, and 308 K are listed in Table IV. The agreement between the calculated and experimental data is excellent at all temperatures, which provides further evidence for the accuracy and transferability of MB-pol.

The melting point (T_m) of hexagonal ice (I_h) at the pressure of 1 atm was calculated from classical MD simulations carried out in the NPT ensemble to characterize the ice/water coexistence as a function of temperature. This approach, originally proposed in Refs. 127–129, has been used to determine the melting point of several water models.^130–135 The water and ice I_h regions of the coexistence system, each containing 360 water molecules, were equilibrated independently at 300 and 100 K, respectively. The two regions were then combined into a single rectangular box of 720 water molecules (dimensions 22.57 × 23.46 × 44.24 ų) with the basal (0001) plane of ice I_h in contact with the liquid phase. Direct coexistence simulations were then carried out in the NPT ensemble at different temperatures, during which the enthalpy of the combined system was monitored. The melting point was determined as the temperature at which the enthalpy of the combined system remained constant throughout the entire trajectory. The classical melting point of MB-pol is found at 263.5 ± 1.5 K. The calculated T_m is ∼9 K lower than the experimental value and in good agreement with the classical MB-pol estimate of 263 K for the temperature of maximum density (see Section III B 1).

Table V lists the lattice energies of different ice phases predicted by MB-pol in comparison with the corresponding experimental estimates. The MB-pol values were obtained from single point calculations carried out for the experimental crystal structures taken from Ref. 89. In all cases, the MB-pol results are in excellent agreement with the experimental values. The maximum deviation of ∼2% is found for ice VIII, a high density phase. A detailed analysis of the lattice energies in terms of many-body contributions will be reported in a forthcoming publication.

The densities of several ice phases calculated from classical NPT simulations with MB-pol at different temperatures and pressures are compared in Table V with the corresponding experimental data.^137–142 Excellent agreement is found between the theoretical predictions and the experimental values, with the largest relative error being less than 3% for ice VIII. It should be noted that the classical NPT simulations with MB-pol slightly overestimates the densities of all ice phases, which again indicates that explicit inclusion of NQE is necessary for more quantitative comparisons with the experimental data. In this context, it was demonstrated that quantum simulations are indeed strictly required to correctly compare the theoretical results with the experimental data at temperatures below 100 K.¹⁴³ Similar to liquid water,⁵⁵ NQE are expected to lower the density and, therefore, improve the agreement with experiment, especially for the ice phases (XIII, XIV, and XV) that are stable at lower temperature.

Scoring. Adopting the same scoring scheme used for the liquid properties (Eq. (17)), scores were assigned to both lattice energies and densities calculated with MB-pol for the different ice phases considered in this analysis, setting the tolerance value to 1% for both properties. For the melting point, it is set to 2.5% which is the same as used in Ref. 4. As shown in Table V, the average scores received by MB-pol are 95 and 89 for ice lattice energies and densities, respectively, which, combined with the analyses reported in Secs. III A–III C, demonstrate that MB-pol consistently achieves high accuracy in predicting the properties of water from clusters to the liquid phase and ice.

In this study, the accuracy of the MB-pol many-body potential is assessed from extensive analysis of the energetics as well as of spectroscopic, structural, thermodynamic, and dynamical properties of water from the gas to the condensed phase.

The analysis of gas-phase properties shows that both individual many-body contributions to the interaction energies and harmonic frequencies calculated for water clusters up to the hexamer with MB-pol are in excellent agreement with reference data obtained at the coupled cluster level of theory. The largest deviations are observed for cyclic clusters due to the accumulation of errors associated with repeating, symmetry-equivalent two-body and three-body units in the MB-pol representation of the underlying Born-Oppenheimer PES.

For the liquid phase, classical MD simulations carried out with MB-pol correctly reproduce the temperature dependence of all structural, thermodynamic, and dynamical properties analyzed in this study. The deviations from the experimental values increase as the temperature decreases. Since MB-pol was derived entirely from electronic structure data and thus represents the water Born-Oppenheimer PES, the differences between results obtained from classical MD simulations with MB-pol and experimental data are expected and confirm that nuclear quantum effects must be explicitly taken into account in the simulations with MB-pol for a quantitative (and rigorous) description of the molecular properties of liquid water. Finally, the melting point of hexagonal ice as well as both lattice energies and densities calculated with MB-pol for several ice phases are found in good agreement with the corresponding experimental data, which provides further evidence for the transferability of MB-pol.

Besides demonstrating the high and, in many respects, unprecedented accuracy with which MB-pol predicts the properties of water across different phases, this study also provides a series of rigorous tests that should be used to assess the ability of both empirical and ab initio water models “to get the right answers for the right reasons,”¹⁴⁴ which is the fundamental prerequisite for a physically correct understanding of the behavior of water at the molecular level.

See supplementary material for list of additional data related to the analysis of many-body interactions and vibrational frequencies of water clusters along with the corresponding molecular coordinates. Tables with numerical values for all water properties are shown in Figures 4–11.

This research was supported by the National Science Foundation through Grant No. CHE-1453204. This work was performed in part under the auspices of the US DOE by LLNL under Contract No. DE-AC52-07NA27344. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357, as well as the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (Grant No. ACI-1053575). S.S. acknowledges the University of California, San Diego for the financial support through the Frontiers of Innovation Scholars Program (FISP).

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.