We present a systematic analysis of state-of-the-art polarizable and flexible water models from a many-body perspective, with a specific focus on their ability to represent the Born–Oppenheimer potential energy surface of water from the gas to the liquid phase. Using coupled cluster data in the completed basis set limit as a reference, we examine the accuracy of the polarizable models in reproducing individual many-body contributions to interaction energies and harmonic frequencies of water clusters and compare their performance with that of MB-pol, an explicit many-body model that has been shown to correctly predict the properties of water across the entire phase diagram. Based on these comparisons, we use MB-pol as a reference to analyze the ability of the polarizable models to reproduce the energy landscape of liquid water under ambient conditions. We find that, while correctly reproducing the energetics of minimum-energy structures, the polarizable models examined in this study suffer from inadequate representations of many-body effects for distorted configurations. To investigate the role played by geometry-dependent representations of 1-body charge distributions in reproducing coupled cluster data for both interaction and many-body energies, we introduce a simplified version of MB-pol that adopts fixed atomic charges and demonstrate that the new model retains the same accuracy as the original MB-pol model. Based on the analyses presented in this study, we believe that future developments of both polarizable and explicit many-body models should continue in parallel and would benefit from synergistic efforts aimed at integrating the best aspects of the two theoretical/computational frameworks.

Since the early days of computer simulations, molecular models of water have provided the scientific community with a means of investigating aqueous environments at the microscopic level.1,2 Despite significant progress achieved in the last two decades, the accurate prediction of the molecular properties of water across the entire phase diagram remains an outstanding challenge, with the overall predictive power of any computer simulation heavily depending on the accuracy with which the underlying model is able to represent the interactions between water molecules.3–6 

Starting from the seminal work by Bernal and Fowler7 and taking advantage of increasing computational power, the field of computer simulations of water has witnessed the continued development of progressively more sophisticated models. Early models, such as ST28 and the subsequent TIPnP9–11 and SPC*12 families of models, have been instrumental in shedding light on the structural, thermodynamic, and dynamical properties of liquid water at the molecular level. These models adopt relatively simple functional forms, representing the water molecules as rigid bodies, and describe electrostatic interactions via Coulomb potentials between effective atomic point charges and repulsive and dispersion interactions via Lennard-Jones potentials.3–6 

Motivated by the recognition of the importance of many-body effects in determining the properties of the hydrogen-bond network, several models have been developed building upon the many-body expansion (MBE) of the total energy, EN, which is expressed as13,
(1)
where N is the number of water molecules and ri collectively denotes the coordinates of the oxygen (O) and hydrogen (H) atoms of the ith water molecule. Equation (1) shows that EN can be rigorously obtained by summing all contributions due to individual n-body energies, with the 1-body energy, ε1B, representing the deformation energy of a water monomer (which is set to zero for rigid models) and εnB corresponding to the n-body energies defined recursively as
(2)
From Eqs. (1) and (2), it follows that the interaction energy, Eint, between N water molecules can be expressed as
(3)
In polarizable models, all many-body energies, from ε3B to εNB in Eq. (1), are effectively accounted for by a single term representing classical many-body polarization.14 Building upon the pioneering work by Caldwell et at.,15 several polarizable water models, with different degrees of sophistication, have been developed and used in molecular dynamics (MD) and Monte Carlo (MC) simulations of aqueous systems.16–58 Among the most popular polarizable water models, which takes into account monomer flexibility, are the Thole type model (TTM)34–40 and AMOEBA45–50 families.

The TTM models are based on Thole’s polarization schemes59 and were primarily fitted to ab initio data for small water clusters.33–40 The development of the TTM models began with the parameterization of a polarizable model derived from ab initio data for the water dimer, which was then used to investigate the properties of water clusters and ice.33 A subsequent analysis of the water dimer potential energy surface (PES) demonstrated that models that were able to more closely reproduce the ab initio PES were also able to more accurately reproduce the experimental second virial coefficient.34 Starting from the polarizable model of Ref. 33, the TTM2-R model,36 a polarizable model of water with rigid monomers, was developed to reproduce the ab initio dimer PES of Ref. 34. TTM2-R was shown to reproduce both gas-phase properties and ab initio data calculated at the second-order Møller–Plesset (MP2) level of theory in the complete basis set (CBS) for small water clusters, as well as several bulk properties, including the diffusion coefficient and the radial distribution functions (RDFs) of liquid water, and the lattice constants of ice Ih.36 Building upon the TTM2-R parameterization, the TTM2-F model37 and its refined TTM2.1-F version38 were the first flexible water models to adopt a highly accurate 1-body term rigorously derived from high-level ab initio potential energy and dipole moment surfaces reported for an isolated water molecule reported in Ref. 60. While maintaining the same level of accuracy as TTM2-R in reproducing several properties of water,61–64 TTM2-F and TTM2.1-F were the first models able to correctly predict the widening of the HOH angle going from small gas-phase clusters to liquid water and ice.65 Subsequent developments led to two distinct TTM models, TTM3-F39 and TTM4-F,40 which were parameterized with specific focus on reproducing the dynamical properties and vibrational spectra of water. Both TTM3-F and TTM4-F were shown to provide better agreement with the experimental infrared spectra of liquid water and ice compared to TTM2-F, although some discrepancies with the experimental spectra were found in both line shapes and vibrational shifts, particularly in the OH stretching frequency region.39,40,66,67

The original AMOEBA model, AMOEBA03, was built upon a polarizable atomic multipole description of electrostatic interactions, including multipoles up to the quadrupole and representing polarization via self-consistently induced atomic dipoles within a modified Thole scheme.45 AMOEBA03 was shown to reproduce several experimental properties and ab initio data for both water clusters and liquid water.46 Subsequent developments and reparameterizations using the ForceBalance algorithm68 led to iAMOEBA48 and AMOEBA14.69 In contrast to AMOEBA03, iAMOEBA uses a direct approximation to polarization, with the induced dipoles being only determined from the permanent multipole electric fields.48 Besides being computationally less expensive, iAMOEBA was shown to provide better agreement with experiment than AMOEBA03 for several properties of water in the gas, liquid, and solid phases.48 It has recently been shown, however, that iAMOEBA is not able to correctly reproduce the sharp increase in the isothermal compressibility in the supercooled regime.70 Extensive comparisons with experimental data showed that AMOEBA14 correctly reproduces the temperature dependence of several structural, thermodynamic, and dynamical properties of water.69 Most recent developments of the AMOEBA models include AMOEBA+49 and AMOEBA+CF.50 Besides providing improved description of many-body polarization, repulsion, and dispersion compared to its predecessor models, AMOEBA+ incorporates an effective description of charge penetration by damping multipole–multipole interactions as well as an explicit pairwise charge transfer term fitted to reproduce the results of energy decompositions calculated at the symmetry-adapted perturbation theory (SAPT) level. AMOEBA+ was shown to perform well in both reproducing individual SAPT components and the properties of liquid water.49 AMOEBA+CF represents a further refinement of AMOEBA+, which was achieved by implementing a geometry-dependent charge flux (CF) scheme within each water molecule. This resulted in a noticeable improvement in the description of both the energetics of water clusters and the properties of liquid water.50 

Closely related to the AMOEBA models is MB-UCB, which was derived from a variational energy decomposition analysis.58 Similar to AMOEBA+, MB-UCB includes terms describing permanent electrostatics, polarization, repulsion, dispersion, and charge transfer, which were parameterized from electronic structure data for the water monomer and small water clusters. MB-UCB was shown to reproduce the relative energies of large datasets of water clusters as well as the structural, thermodynamic, and dynamical properties of liquid water.58 

It was recognized quite early that Eq. (1) also provides a rigorous theoretical framework for the development of explicit many-body models that aim to quantitatively reproduce each individual term of the MBE.71–74 With continued progress in both hardware technology and efficient algorithms for correlated electronic structure calculations, the development of explicit many-body models of water has witnessed a significant acceleration in the past decade. Since, on average, the sum of 2-body and 3-body energies contributes ∼95% to Eint, in these explicit many-body models, low-order n-body energies up to 3-body contributions (ε3B) are treated individually, while all higher-body energies are either approximated by an effective (n − 3)B-term or represented by classical polarization as in polarizable models. To date, the most notable explicit many-body models of water are CC-pol,75–78 which was the first explicit many-body water model with rigid monomers capable of reproducing both the properties of the water dimer and the radial distribution functions (RDFs) of liquid water at 300 K, WHBB,79–82 which was the first explicit many-body water model with flexible monomers, and HBB2-pol and MB-pol,83–87 which were the first explicit many-body water models with flexible monomers capable of correctly predicting the properties of water across all phases.88,89 In particular, MB-pol has been shown to quantitatively reproduce the vibration–rotation tunneling spectrum of the water dimer,85 the energetics, isomeric equilibria, tunneling splittings, and vibrational spectra of small water clusters,86,90–101 the structural, thermodynamic, and dynamical properties as well as the infrared and Raman spectra of liquid water,87,102–106 the sum-frequency generation spectra of the air/water interface,107–111 and the energetics and the infrared and Raman of various ice phases.112–114 More recently, molecular configurations extracted from MD simulations with MB-pol have been used in many-body perturbation theory calculations to model the x-ray absorption spectrum of liquid water115 and to determine the electron affinity of water, both in the bulk and at the air/water interface.116 MB-pol was also used as a reference for the development of an optimized exchange–correlation density functional for water.117 

Given the renewed interest in polarizable water models, we present here a systematic analysis of the TTM, AMOEBA, and MB-UCB models, with particular focus on the ability of these models to reproduce many-body effects in water systems from the gas to the condensed phase. Results obtained with the polarizable models are compared with the corresponding values obtained with MB-pol, and when possible, the accuracy of each model is assessed relative to reference values calculated at the coupled cluster level of theory, including single, double, and perturbative triple excitations, in the CBS limit, i.e., CCSD(T)/CBS, the current “gold standard” for chemical accuracy. To determine the effects of geometry-dependent representations of the 1-body charge distribution on interaction and many-body energies, we also introduce a simplified version of MB-pol (fq-MB-pol), which employs fixed point charges, and assess the accuracy of the new model relative to MB-pol and CCSD(T)/CBS. This study aims to investigate the accuracy of state-of-the-art polarizable and flexible water models through a many-body perspective and determine how possibly shortcomings inherent to a description of many-body effects entirely based on classical polarization may affect the predictions of the properties of water calculated from computer simulations.

This article is organized as follows: Sec. II provides a brief description of both polarizable and explicit many-body models analyzed in this study, while Sec. III provides computational details specific to the comparisons presented in Secs. IV A– IV D. Section IV presents comparisons for interaction and many-body energies calculated for small water clusters, associated harmonic frequencies, and interaction energies for liquid configurations. A general discussion along with an outlook is presented in Sec. V.

The TTM, AMOEBA, and MB-UCB models approximate Eq. (1) as
(4)
where EMBeff effectively accounts for all many-body interactions among N water molecules. As mentioned in the Introduction, the TTM models employ an accurate representation of ε1B(ri) derived from high-level ab initio calculations of the water monomer, which were then refined to reproduce the experimental line positions of the rovibrational spectrum of an isolated water molecule.60 The 1-body term in the AMOEBA and MB-UCB models is instead represented by empirically parameterized anharmonic bond-stretching and angle-bending potentials, along with a bond-angle cross term that was found necessary in the original developments of AMOEBA03 to improve the description of the splitting between symmetric and asymmetric OH stretching vibrations.45 
The TTM, AMOEBA, and MB-UCB models further decompose EMBeff into distinct physical contributions to the interaction energy according to
(5)
where Eelec, Epol, Erep, Edisp, and ECT are individual terms parameterized to represent permanent electrostatics, polarization, repulsion, dispersion, and charge transfer, respectively. Since the TTM, AMOEBA, and MB-UCB models adopt different functional forms to describe each contribution to EMBeff, we briefly summarize here the main features of all models. For specific details, see the original references.

1. TTM2.1-F

In TTM2.1-F, Eelec is described by Coulomb interactions between geometry-dependent point charges placed on the H atoms and the fictitious M-site site located along the bisector of the HOH angle.38 The geometry-dependent point charges were determined to reproduce the ab initio dipole moment surface of an isolated water molecule derived from high-level ab initio calculations in Ref. 60. Epol is represented according to a Thole’s damping scheme,59 with inducible point dipoles placed on the O and H atoms. Both Eelec and Epol in TTM2.1-F are fully damped to mimic charge penetration effects. TTM2.1-F does not include an explicit charge transfer term and combines Erep and Edisp in a single term representing pairwise additive van der Waals interactions between the oxygen atoms,
(6)
where Rij is the distance between the O atoms of molecules i and j, and D, α, A, B, and C are fitting parameters.

2. TTM3-F

Although adopting the same functional form as TTM2.1-F, the TTM3-F model was specifically parameterized to improve the description of vibrational spectra of water.39 For this purpose, the original geometry-dependent point charges of TTM2.1-F were empirically modified to mimic the increase in magnitude of the charges on the H atoms when the OH bonds elongate in liquid water. In contrast to TTM2.1-F, the TTM3-F model employs a single inducible point dipole placed on the M-site as well as a different expression for the van der Waals interactions,39 
(7)
where, as in TTM2.1-F, Rij is the distance between the O atoms of molecules i and j and ε, σ, and λ are fitting parameters.

3. TTM4-F

Directly building upon TTM2-F,37 the TTM4-F model was parameterized to reproduce the polarizability surface of the water monomer calculated at the MP2 level of theory.40 As a result, although TTM4-F employs the same functional form as TTM2-F to describe permanent electrostatics and polarization, it adopts a modified Thole’s damping scheme along with an expanded pairwise additive term representing van der Waals interactions,
(8)
where, as in TTM2.1-F and TTM3-F, Rij is the distance between the O atoms of molecules i and j and A2l are fitting parameters.

4. AMOEBA03

As mentioned above, the 1-body term of AMOEBA03, which includes anharmonic representations of both stretching and bending vibrations along with an additional Urey–Bradley term representing the coupling between stretching and bending vibrations, was parameterized to reproduce the gas-phase vibrational frequencies of the water monomer at the experimental geometry.45 AMOEBA03 employs a multipole expansion, with atomic monopole, dipole, and quadrupole moments placed on each atom, to represent Eelec in Eq. (5). These moments were derived from a distributed multipole analysis carried out on ab initio data for the water monomer. Epol in Eq. (5) is represented through the interactions of inducible point dipoles placed on the atomic sites. While AMOEBA03 adopts a Thole’s scheme to damp Epol at short range, no damping is applied to Eelec. As the TTM models, AMOEBA03 does not include an explicit charge transfer term but, different from the TTM models, includes van der Waals interactions between all pairs of atoms through a sum of pairwise additive buffered 14–7 potentials,
(9)
Here, i and j thus indicate O and H atoms on different water molecules within each dimer, εij is the potential well depth of the corresponding 14–7 potential, and ρij=Rij/Rij0, with Rij being the distance between atoms i and j and Rij0 being the corresponding minimum energy distance. The values of δ = 0.07 and γ = 0.12 were taken from Ref. 118. All fitting parameters in AMOEBA03 were optimized to reproduce ab initio data for gas-phase clusters as well as the temperature and pressure dependence of several properties of liquid water.

5. iAMOEBA

Although iAMOEBA adopts the same functional form as AMOEBA03, it employs a simplified expression for the polarization term in which the induced dipoles are only determined by the electric fields generated by permanent multipoles, and mutual polarization between the induced dipoles is neglected.48 Including only direct polarization implies that iAMOEBA only accounts for at most 3-body interactions where an induced dipole interacts with two other multipole moments. iAMOEBA was parameterized using the ForceBalance algorithm that allows for simultaneously fitting to experimental and ab initio training data.68 

6. AMOEBA14

AMOEBA1469 represents an updated version of the original AMOEBA03 model, which was obtained by using the ForceBalance algorithm68 to determine a refined set of fitting parameters.

7. AMOEBA+

Although it was built upon previous versions of the AMOEBA models, AMOEBA+ employs a different representation of EMBeff in which Thole’s damping is applied not only to Epol but also to the interactions between the permanent multipoles of Eelec.49 In addition, AMOEBA+ includes an explicit term describing charge transfer. Starting from energy decompositions derived from SAPT calculations, all fitting parameters in AMOEBA+ were optimized to reproduce both ab initio and experimental data using the ForceBalance algorithm.68 

8. AMOEBA+CF

Representing a further development of AMOEBA+, AMOEBA+CF adopts a more advanced description of permanent electrostatics which is achieved by implementing a geometry-dependent charge-flux (CF) scheme within each monomer.50 Similar to the TTM models, the charge flux in AMOEBA+CF explicitly depends on the local geometry of a water molecule and was found to improve the transferability of the model across different phases.

9. MB-UCB

Closely related to AMOEBA+ and AMOEBA+CF, the MB-UCB model adopts a similar functional form for Eelec, with fully damped interactions between permanent multipoles, although differently from the AMOEBA models, the multipole expansion is truncated at the dipole moment.58 Furthermore, MB-UCB implements anisotropic dipole polarizabilities to better reproduce anisotropic electric fields in nonhomogenous media. While MB-UCB adopts the same buffered 14–7 expression for EvdW as AMOEBA+ and AMOEBA+CF, differently from the AMOEBA+ models that account only for pairwise charge transfer, MB-UCB includes a many-body charge transfer term that is solved self-consistently.58 All parameters of the MB-UCB model were optimized on training data for the water monomer, dimer, trimer, tetramer, and pentamer clusters, with the ECT term being specifically parameterized to reproduce charge transfer energies derived from a variational absolutely localized molecular orbitals energy decomposition analysis (ALMO-EDA). Based on the analyses reported in Ref. 58, it was claimed that MB-UCB achieves comparable performance to MB-pol.

Both MB-pol and fq-MB-pol approximate Eq. (1) as
(10)
where ε1B represents the same ab initio PES of the water monomer of Ref. 60 adopted by the TTM models, ε2B and ε3B are explicit 2-body and 3-body terms, and ε>3B is an implicit term representing all many-body interactions involving more than three water molecules at a time.

1. MB-pol

The 2-body term of MB-pol includes three contributions,85,
(11)
Here, E2Bsr describes short-range 2-body interactions that are represented by a fourth-degree permutationally invariant polynomial (PIP)119 in variables that are functions of the distances between all pairs of the (six) sites of the MB-pol water monomer. As in TTM2.1-F and TTM4-F, Eelec in MB-pol is described by Coulomb interactions between geometry-dependent point charges that reproduce the ab initio dipole moment surface of a water molecule determined in Ref. 60. Edisp is a 2-body dispersion term expressed as
(12)
where the sum runs over all pairs of O and H atoms i and j on different water molecules within a dimer, Rij is the distance between atoms i and j, f(δij) is the Tang–Toennies damping function with fitting parameter δij,120 and C6,ij are the dispersion coefficients calculated from fits to the asymptotic ab initio reference energies of the water dimer as in CC-pol.75 
In Eq. (10),
(13)
describes 3-body short-range interactions that, as in the analogous 2-body term E3Bsr, are represented by a fourth-degree PIP119 in variables that are functions of the distances between all sites within a trimer.86 The coefficients of the 2-body and 3-body PIPs were optimized using Tikhonov regression (also known as ridge regression)121 to reproduce 2-body and 3-body energies calculated at the CCSD(T)/CBS level of theory. It was shown that the E2Bsr and E3Bsr terms of MB-pol quantitatively reproduce 2-body and 3-body quantum-mechanical effects that cannot be described by purely classical expressions (e.g., charge transfer and penetration and Pauli repulsion).122 Finally, Epol in Eq. (10) represents N-body polarization and is described by a modified version of the Thole-type model adopted by TTM4-F.85,86

2. fq-MB-pol

fq-MB-pol is introduced here to assess the importance of employing 1-body geometry-dependent atomic charges in reproducing CCSD(T)/CBS interaction and many-body energies. For this purpose, fq-MB-pol adopts the same functional form as MB-pol, with explicit representations of short-range 2-body and 3-body interactions in terms of PIPs that are combined with classical representations of permanent electrostatics, polarization, and dispersion. fq-MB-pol was developed from the same sets of CCSD(T)/CBS data for 1-body, 2-body, and 3-body energies used in the development of MB-pol.85,86 However, different from MB-pol that employs geometry-dependent charges to reproduce the ab initio dipole moment surface of Ref. 60, fq-MB-pol employs fixed multipole-derived charges123 calculated at the ωB97M-V/aug-cc-pVQZ level of theory124–126 for an isolated water molecule in the corresponding optimized geometry. The multipole-derived charges were then mapped into the M-site (qM) and the two H atoms (qH) according to the same scheme adopted by MB-pol and originally introduced in the TTM2.1-F model, resulting in qM = −1.151 790e and qH = 0.575 895e.

All TTM2.1-F, TTM3-F, TTM4-F, MB-pol, and fq-MB-pol calculations were carried out using in-house software based on DL_POLY2,127 while all AMOEBA03, AMOEBA14, iAMOEBA, AMOEBA+, and AMOEBA+CF calculations were carried out using Tinker.128 The 2-body and 3-body CCSD(T)/CBS reference energies are taken from Refs. 85 and 86. The reference structures for the water clusters, along with the associated interaction and many-body energies calculated at the CCSD(T)/CCSD(T)-F12 level of theory within the SAMBA scheme,129 are taken from Ref. 88. The reference harmonic frequencies calculated at the 2-body:many-body CCSD(T):MP2 level of theory are taken from Refs. 130 and 131, while the harmonic frequencies for all polarizable and explicit many-body models were calculated from the corresponding optimized cluster structures.

To assess the ability of the different models considered in this study to represent many-body effects in water from the gas to the condensed phase, Secs. IV A–IV D present a systematic analysis of interaction energies, many-body energy decompositions, and harmonic frequencies for small clusters (Fig. 1), as well as relative energies for configurations extracted from centroid molecular dynamics (CMD)132–136 simulations of liquid water carried out under ambient conditions with MB-pol in Ref. 105.

FIG. 1.

Structures of the (H2O)n clusters, with n = 1–6, used in the analysis of interaction and many-body energies as well as harmonic frequencies.

FIG. 1.

Structures of the (H2O)n clusters, with n = 1–6, used in the analysis of interaction and many-body energies as well as harmonic frequencies.

Close modal

The interaction energies [Eint in Eq. (3)] for the minimum-energy structures of the water dimer and trimer (isomers 1 in Fig. 1) calculated using both polarizable and explicit many-body models are listed in Table I. Also reported in Table I are the corresponding 2-body and 3-body energies contributing to the interaction energy of the water trimer. We note here that, by definition [Eq. (3)], the interaction energy for the water dimer also corresponds to the 2-body energy.

TABLE I.

Interaction energies (Eint) for the minimum-energy structures of the water dimer and trimer [isomers 1 of (H2O)2 and (H2O)3 in Fig. 1]. Also listed are the 2-body and 3-body energies contributing to Eint of the water trimer. All energies are in kcal/mol.

(H2O)2(H2O)3
ModelE2E3E2BE3B
CCSD(T) −5.06 −16.21 −13.74 −2.47 
TTM2.1-F −5.07 −15.82 −13.59 −2.22 
TTM3-F −5.30 −16.23 −15.01 −1.21 
TTM4-F −5.03 −14.31 −11.80 −2.51 
AMOEBA03 −5.21 −16.00 −13.18 −2.81 
iAMOEBA −5.22 −15.58 −13.91 −1.67 
AMOEBA14 −4.91 −16.23 −12.99 −3.25 
AMOEBA+ −5.20 −17.09 −14.29 −2.80 
AMOEBA+CF −5.00 −16.32 −13.75 −2.57 
MB-UCB −5.44 −16.40 −14.21 −2.20 
fq-MB-pol −5.07 −16.20 −13.74 −2.46 
MB-pol −5.05 −16.15 −13.69 −2.45 
(H2O)2(H2O)3
ModelE2E3E2BE3B
CCSD(T) −5.06 −16.21 −13.74 −2.47 
TTM2.1-F −5.07 −15.82 −13.59 −2.22 
TTM3-F −5.30 −16.23 −15.01 −1.21 
TTM4-F −5.03 −14.31 −11.80 −2.51 
AMOEBA03 −5.21 −16.00 −13.18 −2.81 
iAMOEBA −5.22 −15.58 −13.91 −1.67 
AMOEBA14 −4.91 −16.23 −12.99 −3.25 
AMOEBA+ −5.20 −17.09 −14.29 −2.80 
AMOEBA+CF −5.00 −16.32 −13.75 −2.57 
MB-UCB −5.44 −16.40 −14.21 −2.20 
fq-MB-pol −5.07 −16.20 −13.74 −2.46 
MB-pol −5.05 −16.15 −13.69 −2.45 

At the 2-body level, all polarizable models predict interaction energies for the water dimer that are close to the CCSD(T)/CBS reference value. Only the TTM3-F and MB-UCB models predict noticeably more attractive interactions, with 2-body energies of −5.30 kcal/mol and −5.44 kcal/mol, respectively, which are ∼5% lower than the reference CCSD(T)/CBS value. Although most polarizable models predict reasonably accurate interaction energies for the water trimer, large deviations from the CCSD(T)/CBS reference values are found for TTM4-F, which underestimates the interaction strength by 1.90 kcal/mol, and AMOEBA+, which instead overestimates the interaction strength by 0.88 kcal/mol. For both water dimer and trimer, MB-pol and fq-MB-pol provide quantitative agreement with CCSD(T)/CBS, with deviations being always within 0.06 kcal/mol of the reference values.

When the interaction energy of the water trimer is decomposed into the corresponding individual many-body contributions, it is found that TTM3-F overestimates 2-body effects and underestimates 3-body effects by more than 1.0 kcal/mol. Similar trend, although to a lesser extent, is found for AMOEBA14, which instead underestimates 2-body effects and overestimates 3-body effects by ∼0.8 kcal/mol. This analysis indicates that the apparent agreement with the CCSD(T)/CBS interaction energy for the water trimer exhibited by TTM3-F and AMOEBA14 results from nearly perfect error compensation between 2-body and 3-body effects.

On the other hand, TTM4-F significantly underestimates 2-body energies by ∼2.0 kcal/mol but is able to quantitatively reproduce 3-body energies. The opposite trend is predicted by iAMOEBA, which is able to correctly reproduce 2-body energies but underestimates 3-body contributions by ∼1.0 kcal/mol. The different performance exhibited by TTM4-F and iAMOEBA in reproducing 3-body energies is directly correlated with the different representations used by these two models to describe polarization effects, with TTM4-F being specifically parameterized to reproduce the polarizability of a water molecule and iAMOEBA neglecting mutual polarization between water molecules.

Table I shows that both MB-pol and fq-MB-pol quantitatively reproduce the CCSD(T)/CBS reference values for the 2-body and 3-body energies of the water trimer. This demonstrates that the ability of these models to accurately predict the interaction energies of both water dimer and trimer thus results from a quantitative account of both 2-body and 3-body effects.

While the comparisons reported in Table I demonstrate that, although with different levels of accuracy, most polarizable models are able to represent 2-body and 3-body energies for the minimum-energy structures of the water dimer and trimer, they do not provide any information about 2-body and 3-body effects for distorted configurations that exist in liquid water at different temperatures and pressures. Therefore, to assess the ability of both polarizable and explicit many-body models to describe the multidimensional 2-body and 3-body energy landscapes, Figs. 2 and 3 show correlations between 2-body and 3-body energies calculated with each model and the corresponding CCSD(T)/CBS reference values for 1000 dimer and trimer configurations randomly extracted from the original MB-pol training sets.85,86 To guarantee fair comparisons between all models, the MB-pol and fq-MB-pol models shown in Figs. 2 and 3 were re-fitted to the reduced 2-body and 3-body training sets obtained after removing from the original training sets the same 1000 configurations used for the correlation analysis.

FIG. 2.

Correlation plots for the 2-body energies. Plotted on the x axes are the 2-body CCSD(T)/CBS reference energies from Ref. 85. On the y axes are the corresponding 2-body energies calculated with both polarizable and explicit many-body models. Also shown in the bottom right panel is the distribution of 2-body energies in the test set used for this analysis.

FIG. 2.

Correlation plots for the 2-body energies. Plotted on the x axes are the 2-body CCSD(T)/CBS reference energies from Ref. 85. On the y axes are the corresponding 2-body energies calculated with both polarizable and explicit many-body models. Also shown in the bottom right panel is the distribution of 2-body energies in the test set used for this analysis.

Close modal
FIG. 3.

Correlation plots for the 3-body energies. Plotted on the x axes are the 3-body CCSD(T)/CBS reference energies from Ref. 86. On the y axes are the corresponding 3-body energies calculated with both polarizable and explicit many-body models. Also shown in the bottom right panel is the distribution of 3-body energies in the test set used for this analysis.

FIG. 3.

Correlation plots for the 3-body energies. Plotted on the x axes are the 3-body CCSD(T)/CBS reference energies from Ref. 86. On the y axes are the corresponding 3-body energies calculated with both polarizable and explicit many-body models. Also shown in the bottom right panel is the distribution of 3-body energies in the test set used for this analysis.

Close modal

Figure 2 shows that all polarizable models exhibit significant deviations from the reference CCSD(T)/CBS 2-body energies. This is particularly evident for TTM4-F and AMOEBA+ that are associated with RMSDs of 2.28 kcal/mol and 4.38 kcal/mol, respectively. Overall, among the polarizable models, TTM3-F, AMOEBA03, iAMOEBA, AMOEBA+CF, and MB-UCB display similar RMSDs of ∼1.2 kcal/mol–1.3 kcal/mol, while TTM2.1-F and AMOEBA14 provide the closer agreement with the CCSD(T)/CBS reference energies, with RMSDs of ∼1.0 kcal/mol. In contrast, both MB-pol and fq-MB-pol exhibit the same accuracy as CCSD(T)/CBS over the entire range of 2-body energies, with RMSDs of ∼0.05 kcal/mol, which are more than one order of magnitude smaller than those associated with the polarizable models.

Importantly, both AMOEBA+CF and MB-UCB do not display the large (positive) deviations from the CCSD(T)/CBS 2-body energies that are instead associated with AMOEBA+. Considering the progressive development of the AMOEBA models described in Sec. II, these differences can be traced back to the fact that both AMOEBA+CF and MB-UCB include explicit terms parameterized to mimic charge transfer, which improves the description of short-range interactions in more distorted water dimers. In this regard, it should be noted that, although these parameterized charge transfer terms adopt relatively simple functionals and, therefore, are not sufficient to achieve CCSD(T)/CBS accuracy over a wide range of 2-body energies as shown in Fig. 2, they effectively mimic the same role played by the 2-body PIP adopted by MB-pol and fq-MB-pol. Using a multidimensional 2-body PIP to supplement a classical description of permanent electrostatics, polarization and dispersion were found to be necessary to guarantee an accurate representation of the global energy landscape of the water dimer over a wide range of molecular distortions.85,122

Figure 3 shows that, with the exception of TTM3-F, all polarizable models considered in this study are able to reproduce CCSD(T)/CBS 3-body energies relatively better than the corresponding 2-body energies. Specifically, all models are associated RMSDs smaller than ∼0.2 kcal/mol relative to the reference values. This is in agreement with previous observations that 3-body energies in water are primarily due to polarization, which can be properly accounted for using classical representations based on physically motivated parameterizations of molecular polarizabilities.83 TTM3-F not only displays the largest RMSDs (∼0.4 kcal/mol) but also provides a poor correlation with the reference CCSD(T)/CBS values, significantly underestimating attractive 3-body energies. Since Fig. 3 shows that the closely related TTM2.1-F and TTM4-F models perform significantly better than TTM3-F, the large deviations associated with TTM3-F can be attributed to the fact that this model employs a single polarizable site that may not be sufficient for a correct representation of the anisotropy of 3-body effects within a Thole’s damping scheme.

Figure 3 shows that all AMOEBA models, with the exception of iAMOEBA, tend to slightly overestimate the strength of 3-body interactions, which results in a slope greater than 1. Among the polarizable models, AMOEBA+CF overall provides the best correlation with the CCSD(T)/CBS reference data with an RMSD of ∼0.13 kcal/mol. As in the case of 2-body energies, Fig. 3 shows that both MB-pol and fq-MB-pol also provide quantitative agreement with the CCSD(T)/CBS 3-body energies, with RMSDs of ∼0.04 kcal/mol. Considering that the core representation of 3-body energies in both MB-pol and fq-MB-pol builds upon a slightly modified version of the polarization scheme adopted by TTM4-F, the different level of accuracy achieved by MB-pol and fq-MB-pol compared to that of polarizable models can be attributed to the ability of the 3-body PIP to correct for intrinsic deficiencies associated with purely classical representations of 3-body interactions based on many-body polarization. This is particularly relevant at short range where the overlap of the monomer electron densities is responsible for 3-body effects (e.g., 3-body exchange–repulsion, correlation, and dispersion, as well as more subtle effects such as coupling of exchange and correlation to the dispersion interaction) that do not have well-defined classical analogs.

The deviations from the reference many-body energies calculated for (H2O)n clusters, with n = 4–6, using the TTM, AMOEBA, MB-UCB, and MB-pol models are compared in Fig. 4. As shown in Eq. (3), the interaction energy for each cluster is defined as the difference between the total energy of the cluster and the energies of the individual water molecules kept at the same geometry as in the cluster. By removing the contributions associated with the monomer distortions, the interaction energies thus allow for direct comparisons of molecular interactions as predicted by the different models, without being “contaminated” by differences in the description of the 1-body energies specific to each model.

FIG. 4.

Deviations, ΔE, from the CCSD(T)/CCSD(T)-F12 reference values for individual many-body energies calculated for (H2O)n clusters, with n = 4–6, using both polarizable and explicit many-body models.

FIG. 4.

Deviations, ΔE, from the CCSD(T)/CCSD(T)-F12 reference values for individual many-body energies calculated for (H2O)n clusters, with n = 4–6, using both polarizable and explicit many-body models.

Close modal

Figure 4 shows that TTM2.1-F, TTM3-F, AMOEBA03, iAMOEBA, and AMOEBA14 benefit from significant error compensation between individual many-body terms. For example, the absolute 2-body deviations for the hexamer isomers are as large as 2.65 kcal/mol for TTM2.1-F, 5.31 kcal/mol for TTM3-F, 1.38 kcal/mol for AMOEBA03, 3.91 kcal/mol for iAMOEBA, and 3.81 kcal/mol for AMOEBA14. These models are correspondingly associated with large, but of opposite sign, 3-body deviations whose absolute values are up to 2.85 kcal/mol for TTM2.1-F, 5.12 kcal/mol for TTM3-F, 0.95 kcal/mol for AMOEBA03, 3.91 kcal/mol for iAMOEBA, and 2.31 kcal/mol for AMOEBA14. The observed error compensation in these models stems primarily from an insufficient representation of 2-body permanent electrostatics and van der Waals energies that need to be compensated by unphysically large contributions from many-body polarization, which is particularly evident at the 3-body level. Interestingly, the analysis of many-body energies shows that TTM3-F displays notably larger errors than its predecessor model TTM2.1-F. This can be partially attributed to the TTM3-F model employing a single polarizable site, which was a choice made during the parameterization as a compromise between accuracy and computational efficiency but may be insufficient to correctly reproduce anisotropic electric fields.39 In the AMOEBA03, iAMOEBA, and AMOEBA14 models, the aforementioned deficiencies in the representation of permanent electrostatics may be due to the adoption of undamped interactions between the multipoles. Since, by construction, iAMOEBA only includes up to 3-body contributions, the iAMOEBA errors for all εnB terms with n > 3 are equal, but opposite in sign to the corresponding CCSD(T)/CCSD(T)-F12 reference values. It should be noted that AMOEBA14 exhibits larger deviations than the original AMOEBA03 model, indicating that the representation of permanent electrostatics and polarization worsened during the parameterization of AMOEBA14 when the ForceBalance algorithm was used to simultaneously reproduce both ab initio and experimental data.

On the other hand, TTM4-F, AMOEBA+, AMOEBA+CF, and MB-UCB do not benefit from significant error compensation and exhibit, overall, relatively smaller deviations than their predecessor models. For example, the absolute deviations (AD) on the 2-body energies of the hexamer isomers are up to 4.0 kcal/mol for TTM4-F, 1.70 kcal/mol for AMOEBA+, 0.91 kcal/mol for AMOEBA+CF, and 2.06 kcal/mol for MB-UCB, while the corresponding absolute 3-body deviations are as large as 0.70 kcal/mol for TTM4-F, 0.90 kcal/mol for AMOEBA+, 0.30 kcal/mol for AMOEBA+CF, and 0.73 kcal/mol for MB-UCB. The improved description of many-body effects in these polarizable models can be attributed to the adoption of more advanced representations of both permanent electrostatics and polarization. Specifically, TTM4-F uses the same fully damped scheme for the electrostatics as TTM2.1-F and TTM3-F but employs a more flexible damping function along with distinct Thole parameters for intermolecular and intramolecular interactions. On the other hand, contrary to previous AMOEBA models, both AMOEBA+, AMOEBA+CF and MB-UCB adopt fully damped expressions for permanent electrostatics and include explicit charge penetration and charge transfer terms. Interestingly, the analysis of many-body energies for the water clusters indicates that the MB-UCB model predicts relatively larger 2-body errors for non-cyclic than cyclic isomers of the pentamer and hexamer clusters. This suggests that MB-UCB describes relatively more accurately water arrangements with cooperative hydrogen-bonding than those with anti-cooperative hydrogen-bonding, which is likely due to MB-UCB using anisotropic atomic polarizabilities that provide better representations of planar structures with concerted water arrangements.

Finally, independent of the cluster size and structure, MB-pol, which includes explicit 2-body and 3-body terms integrated with a classical representation of polarization similar to that adopted by the TTM4-F model, exhibits small deviations for each individual n-body energy. In particular, the 2-body and 3-body deviations for the hexamer isomers never exceed 0.27 kcal/mol and 0.14 kcal/mol, respectively, with the larger deviations (0.33 kcal/mol) being associated with 4-body energies. Importantly, fq-MB-pol, which adopts a simplified representation of permanent electrostatics, exhibits 2-body and 3-body deviations (∼0.15 kcal/mol) comparable to those associated with those predicted by the original MB-pol model.

Overall, the deviations on the interaction energies shown in Fig. 5 are found to increase with cluster size for all models, since the errors associated with each individual term of the MBE adds up. Importantly, the cyclic isomers of each clusters are associated, on average, with the largest deviations. Since these isomers are characterized by extended hydrogen-bonded chains of (nearly) equivalent water molecules, the difficulties in reproducing the CCSD(T)/CCSD(T)-F12 reference energies suggest that purely classical polarization terms may be insufficient for a quantitative description of water structures dominated by cooperative effects.

FIG. 5.

Comparisons between the interaction energies of the low-lying isomers of (H2O)n clusters, with n = 4–6 (see Fig. 1), obtained at the CCSD(T)/CCSD(T)-F12 level of theory within the SAMBA scheme in Ref. 88 and the corresponding values calculated with both polarizable and explicit many-body models.

FIG. 5.

Comparisons between the interaction energies of the low-lying isomers of (H2O)n clusters, with n = 4–6 (see Fig. 1), obtained at the CCSD(T)/CCSD(T)-F12 level of theory within the SAMBA scheme in Ref. 88 and the corresponding values calculated with both polarizable and explicit many-body models.

Close modal
Motivated by a similar analysis carried out for density functional theory (DFT) models of water,137 we apply the same scoring scheme developed in Ref. 88 to quantify the accuracy of both polarizable and explicit many-body models in reproducing the CCSD(T)/CCSD(T)-F12 reference data. The resulting scores assigned to each model based on the ability to reproduce interaction and many-body energies for (H2O)n clusters with n = 4, 5, and 6 are reported in Table II. In this analysis, the maximum total unsigned error over all isomers of each (H2O)n cluster is defined as
(14)
The maximum absolute error in the interaction energy over all isomers of each (H2O)n cluster is defined as
(15)
The relative energy difference between isomer 1 and isomer 2 of each (H2O)n cluster is defined as
(16)
Following Ref. 88, the scores for the maximum total unsigned error and maximum absolute error are equal to 100 points if these errors are less than 1.0 kcal/mol, with additional ten points deducted for each increment of 0.5 kcal/mol in the error. Using the CCSD(T)/CCSD(T)-F12 reference values of −0.87 kcal/mol, −1.29 kcal/mol, and −0.29 kcal/mol for the tetramer, pentamer, and hexamer clusters, respectively,88 a score of 100 points is assigned if the relative energy difference between isomers 1 and 2 is within 0.1 kcal/mol of the reference value and ten points are deducted for every additional 0.1 kcal/mol thereafter. If a positive value is predicted for the relative energy difference, zero points are assigned to the model. For scores that are negative, the sign of the relative energy difference is correct, but the difference is greater than 1.0 kcal/mol from the reference, and the value falls outside the bounds of chemical accuracy.
TABLE II.

Errors and associated scores assessing the accuracy of both polarizable and explicit many-body models in describing interaction and many-body energies for (H2O)n clusters, with n = 4–6. ΔEunsigned,maxMB [Eq. (14)] and ΔEint,maxMB [Eq. (15)] are the maximum total unsigned error and maximum absolute total error in the many-body energy decomposition, respectively, while ΔEdiff [Eq. (16)] is the energy difference between the interaction energies of isomers 1 and 2 of each cluster. All energies are in kcal/mol.

ModelClusterΔEunsigned,maxMBScoreΔEint,maxMBScoreΔEdiffScore
 (H2O)4 1.78 80 1.30 90 −0.58 70 
TTM2.1-F (H2O)5 3.90 40 1.50 80 −1.55 70 
 (H2O)6 5.45 10 1.86 80 0.64 
 (H2O)4 5.87 0.67 100 −0.67 80 
TTM3-F (H2O)5 9.06 −70 0.51 100 −0.95 60 
 (H2O)6 11.18 −110 1.12 90 −0.21 100 
 (H2O)4 3.15 50 3.15 50 −0.65 70 
TTM4-F (H2O)5 3.67 40 3.67 40 −3.00 −80 
 (H2O)6 4.46 30 4.45 30 1.96 
 (H2O)4 1.58 80 0.49 100 −0.83 100 
AMOEBA03 (H2O)5 1.75 80 0.73 100 −1.29 100 
 (H2O)6 2.94 60 0.59 100 −0.65 60 
 (H2O)4 4.50 20 1.16 90 −0.95 100 
iAMOEBA (H2O)5 7.47 −30 1.38 90 −0.72 40 
 (H2O)6 9.66 −80 2.14 70 −0.75 50 
 (H2O)4 3.48 50 0.70 100 −0.57 60 
AMOEBA14 (H2O)5 4.79 20 0.82 100 −1.14 80 
 (H2O)6 7.05 −30 0.59 100 −0.03 70 
 (H2O)4 1.58 80 1.58 80 −1.56 30 
AMOEBA+ (H2O)5 2.43 70 2.40 70 −0.96 60 
 (H2O)6 2.32 70 2.27 70 −0.88 40 
 (H2O)4 0.47 100 0.43 100 −0.71 80 
AMOEBA+CF (H2O)5 0.96 100 0.60 100 −1.09 70 
 (H2O)6 1.72 80 0.53 100 −0.07 70 
 (H2O)4 1.77 80 0.83 100 −0.72 80 
MB-UCB (H2O)5 3.15 50 1.34 90 −0.91 60 
 (H2O)6 2.85 60 1.39 90 −0.59 70 
 (H2O)4 0.14 100 0.14 100 −0.79 100 
fq-MB-pol (H2O)5 0.31 100 0.31 100 −1.01 70 
 (H2O)6 0.62 100 0.62 100 −0.31 100 
 (H2O)4 0.26 100 0.26 100 −0.80 100 
MB-pol (H2O)5 0.51 100 0.51 100 −1.10 80 
 (H2O)6 0.84 100 0.84 100 −0.32 100 
ModelClusterΔEunsigned,maxMBScoreΔEint,maxMBScoreΔEdiffScore
 (H2O)4 1.78 80 1.30 90 −0.58 70 
TTM2.1-F (H2O)5 3.90 40 1.50 80 −1.55 70 
 (H2O)6 5.45 10 1.86 80 0.64 
 (H2O)4 5.87 0.67 100 −0.67 80 
TTM3-F (H2O)5 9.06 −70 0.51 100 −0.95 60 
 (H2O)6 11.18 −110 1.12 90 −0.21 100 
 (H2O)4 3.15 50 3.15 50 −0.65 70 
TTM4-F (H2O)5 3.67 40 3.67 40 −3.00 −80 
 (H2O)6 4.46 30 4.45 30 1.96 
 (H2O)4 1.58 80 0.49 100 −0.83 100 
AMOEBA03 (H2O)5 1.75 80 0.73 100 −1.29 100 
 (H2O)6 2.94 60 0.59 100 −0.65 60 
 (H2O)4 4.50 20 1.16 90 −0.95 100 
iAMOEBA (H2O)5 7.47 −30 1.38 90 −0.72 40 
 (H2O)6 9.66 −80 2.14 70 −0.75 50 
 (H2O)4 3.48 50 0.70 100 −0.57 60 
AMOEBA14 (H2O)5 4.79 20 0.82 100 −1.14 80 
 (H2O)6 7.05 −30 0.59 100 −0.03 70 
 (H2O)4 1.58 80 1.58 80 −1.56 30 
AMOEBA+ (H2O)5 2.43 70 2.40 70 −0.96 60 
 (H2O)6 2.32 70 2.27 70 −0.88 40 
 (H2O)4 0.47 100 0.43 100 −0.71 80 
AMOEBA+CF (H2O)5 0.96 100 0.60 100 −1.09 70 
 (H2O)6 1.72 80 0.53 100 −0.07 70 
 (H2O)4 1.77 80 0.83 100 −0.72 80 
MB-UCB (H2O)5 3.15 50 1.34 90 −0.91 60 
 (H2O)6 2.85 60 1.39 90 −0.59 70 
 (H2O)4 0.14 100 0.14 100 −0.79 100 
fq-MB-pol (H2O)5 0.31 100 0.31 100 −1.01 70 
 (H2O)6 0.62 100 0.62 100 −0.31 100 
 (H2O)4 0.26 100 0.26 100 −0.80 100 
MB-pol (H2O)5 0.51 100 0.51 100 −1.10 80 
 (H2O)6 0.84 100 0.84 100 −0.32 100 

For each cluster, a small value for ΔEint,maxMB accompanied by a large value for ΔEunsigned,maxMB indicates that the model benefits from error compensation between individual many-body contributions. Table II shows that TTM3-F displays the largest difference between the scores of the maximum unsigned error (−110) and maximum absolute error (90) for the hexamer cluster. This indicates that, while the total interaction energy is apparently accurate, it is the result of significant error compensation between individual terms of the MBE. Similar trends, although less pronounced, are found for iAMOEBA and AMOEBA14. On the other hand, among the polarizable models, AMOEBA+CF and MB-UCB earn high scores for both ΔEint,maxMB and ΔEunsigned,maxMB. It should be noted that MB-pol and fq-MB-pol score perfectly for each cluster size, indicating that these models do not rely on error compensation to correctly reproduce the corresponding CCSD(T)/CCSD(T)-F12 reference interaction energies.

In the case of the tetramer, all models predict the correct sign for ΔEdiff, and the associated errors fall within 1.0 kcal/mol of the reference value, with only AMOEBA+ having a notably smaller score of 30. For the pentamer, the TTM4-F score is −80 points, indicating that, while getting the ranking of the isomers correct, this model predicts an overly large energy difference. In the case of the hexamer, both the TTM2.1-F and TTM4-F models score zero points, meaning that they predict an incorrect order for isomers 1 (prism) and 2 (cage). AMOEBA03 earns a perfect score for the tetramer and pentamer clusters but fails at correctly reproducing the energetics of the hexamer isomers. MB-pol and fq-MB-pol, on the other hand, score perfectly for the tetramer and hexamer clusters but perform slightly less well for the pentamer cluster. Since MB-pol and fq-MB-pol are built on top of a classical representation of permanent electrostatics and polarization similar to that adopted by TTM4-F, their performance for the pentamer isomers indicates that the 2-body and 3-body PIPs have more difficulties compensating for the relatively large errors predicted by TTM4-F for these clusters.

While the comparisons reported in the previous sections allow for a quantitative assessment of the ability of the different models to represent interaction and many-body energies for various minimum-energy structures of small water clusters, they do not provide any information about the overall morphology of the underlying multidimensional PES. Specific insight into the PES curvature in the neighborhoods of the minimum-energy structure of each isomer can be obtained from the analysis of the corresponding harmonic frequencies. For this purpose, Fig. 6 shows the deviations associated with both polarizable and explicit many-body models relative to the reference 2-body:many-body CCSD(T):MP2 harmonic frequencies reported in Refs. 130 and 131.

FIG. 6.

Deviations from reference 2-body:many-body harmonic frequencies of (H2O)n clusters, with n = 1–6,130,131 calculated with both polarizable and explicit many-body models. The light red and blue backgrounds correspond to red and blue shifts relative to the reference values. Each harmonic mode, indicated by a stick line, is colored based on the associated frequency according to the color bar displayed at the top of the figure.

FIG. 6.

Deviations from reference 2-body:many-body harmonic frequencies of (H2O)n clusters, with n = 1–6,130,131 calculated with both polarizable and explicit many-body models. The light red and blue backgrounds correspond to red and blue shifts relative to the reference values. Each harmonic mode, indicated by a stick line, is colored based on the associated frequency according to the color bar displayed at the top of the figure.

Close modal

Among the polarizable models, all TTM models, which employ a 1-body term exhibiting spectroscopic accuracy,60 are able to quantitatively reproduce the harmonic frequencies of an isolated water molecule. However, differences in the representation of many-body effects discussed in the previous sections lead to large variations in the ability of these models to reproduce the harmonic frequencies of larger clusters. In particular, independent of the cluster size and structure, TTM2.1-F systematically overestimates the harmonic frequencies of both free and hydrogen-bonded OH stretches. TTM3-F displays the opposite behavior, predicting relatively large redshifts for all OH stretching vibrations, while TTM4-F exhibits intermediate behavior, overestimating the vibrational frequencies of the free OH stretches while underestimating those of the hydrogen-bonded OH stretches. These trends indicate that TTM2.1-F predicts hydrogen bonds that are too weak while both TTM3-F and TTM4-F overestimate the strength of the hydrogen bonds. All TTM models accurately reproduce the harmonic frequencies of the bending vibrations. Qualitative differences are found in the low-frequency region of the spectra, with TTM4-F, in general, overestimating (by up to ∼150 cm−1) the harmonic frequencies of vibrational modes between 500 cm−1 and 1000 cm−1 and both TTM2.1-F and TTM3-F systematically underestimating the harmonic frequencies of most of these modes. Since vibrational modes below 1000 cm−1 correspond to collective rearrangements of the hydrogen-bond network, these differences suggest that all TTM models, for different reasons, may have difficulties in correctly representing the hydrogen-bond dynamics in liquid water and its temperature dependence, which is in agreement with previous observations.64,67,138

Figure 6 shows that AMOEBA03, AMOEBA+, AMOEBA+CF, and MB-UCB tend to largely underestimate, up to 200 cm−1, the frequencies of the free OH stretching vibrations above 3700 cm−1, which is traced back to the inability of these models to correctly predict the harmonic frequencies of an isolated water molecule. On the other hand, while these models still underestimate by ∼100 cm−1 the frequencies of weakly hydrogen-bonded OH stretches (at ∼3500 cm−1), they significantly overestimate, by up to ∼150 cm−1–200 cm−1, the frequencies of strongly hydrogen-bonded OH stretches of the prism, cage, book-1, and cyclic chair isomers of the hexamer cluster. Considering that these models underestimate the harmonic frequencies of the symmetric and asymmetric OH stretches of an isolated water molecule, the relatively large blueshifts predicted for OH stretches at ∼3200 cm−1 imply that the effective “solvochromatic shifts” for these vibrational frequencies can be up to ∼400 cm−1, pointing to some deficiencies of these models in representing strong hydrogen bonds. AMOEBA03, AMOEBA+, AMOEBA+CF, and MB-UCB predict similar deviations within 100 cm−1 of the reference values for the bending vibrations (∼1700 cm−1), although larger variations are found in the ability of these models to reproduce collective vibrations (below 1000 cm−1). In particular, while AMOEBA03, AMOEBA+, and AMOEBA+CF still exhibit deviations within 100 cm−1 of the reference values, relatively larger deviations of up to ∼150 cm−1 are associated with MB-UCB in this frequency region. These relatively large deviations suggest that, as the TTM models, MB-UCB may not be able to correctly describe the dynamical rearrangements of the water hydrogen-bond network.

As shown in Fig. 6, independent of the cluster size, both MB-pol and fq-MB-pol accurately reproduce the harmonic frequencies of all vibrational modes, with deviations always being within 50 cm−1 of the reference values. Since MB-pol and fq-MB-pol adopt the same 1-body term as the TTM models, the higher accuracy displayed in reproducing the harmonic frequencies is a direct consequence of both MB-pol and fq-MB-pol being able to quantitatively reproduce each individual term of the MBE in Eq. (1), as discussed in Sec. IV B.

Table III reports the average deviations (AD), average absolute deviations (AAD), and maximum absolute deviations (MAD) from the reference 2-body:many-body CCSD(T):MP2 harmonic frequencies130,131 calculated for the same water clusters analyzed in Fig. 6. To quantify the accuracy of each model in reproducing the reference values, we follow Ref. 88 and assign a score of 100 if the MAD is below 20 cm−1 and deduct ten points for each additional increment of 20 cm−1, until zero points are reached. Since the TTM models use a highly accurate 1-body term, they score perfectly for the water monomer. However, the agreement with the reference harmonic frequencies deteriorates significantly for larger clusters, which can be attributed to deficiencies in the representation of the underlying interactions as discussed in Sec. IV B. On the other hand, the AMOEBA and MB-UCB models tend to suffer from large average absolute errors, with all models, except iAMOEBA, benefiting from error compensation as indicated by smaller average errors. Since these models already exhibit large deviations from the reference harmonic frequencies of the water monomer, the deviations observed for larger clusters result from inaccuracies in the description of both intramolecular and intermolecular potential energy surfaces. It should be noted that all polarizable models exhibit large maximum deviations on the order of 200 cm−1–400 cm−1, which, as shown in Fig. 6, are associated with OH stretching vibrations. Since these vibrations report directly on the local hydrogen-bonding environment,139 these large deviations indicate that the TTM, AMOEBA, and MB-UCB models may not be able to correctly describe the structural and dynamical properties of the hydrogen-bond network. In contrast, both MB-pol and fq-MB-pol perform uniformly well for all cluster sizes and structures, with the corresponding MAD never exceeding 50 cm−1, which is in agreement with previous studies of vibrational spectra of water in both gas and condensed phases.94,102,105–108,113,114

TABLE III.

Average deviation (AD), average absolute deviation (AAD), and maximum absolute deviation (MAD) in cm−1 for the harmonic frequencies of selected isomers of (H2O)n clusters with n = 1–6 calculated with both polarizable and explicit many-body models.

H2O(H2O)2(H2O)3(H2O)4(H2O)5(H2O)6(H2O)6(H2O)6(H2O)6
isomer 1isomer 1isomer 1isomer 1isomer 1isomer 1isomer 2isomer 3isomer 6
 AD −1.38 −5.39 −0.44 −1.01 1.59 2.25 2.12 1.75 3.80 
TTM2.1-F AAD 1.38 14.57 27.47 60.28 62.91 46.32 50.88 59.79 60.70 
 MAD 2.36 33.66 97.25 213.33 252.42 384.57 352.9 293.6 259.55 
 Score 100 90 60 −20 −90 −70 −40 −20 
 AD −1.38 −31.71 −56.61 −66.37 −65.18 −62.14 −64.37 −64.44 −62.25 
 AAD 1.38 31.71 56.62 66.37 65.22 64.64 64.37 64.44 62.26 
TTM3-F MAD 2.36 145.79 240.43 174.67 180.29 193.36 191.60 193.22 182.82 
 Score 100 30 −20 20 10 10 10 10 10 
 AD −1.38 7.75 7.75 10.41 1.82 7.29 9.32 9.13 −2.88 
 AAD 1.38 30.19 46.30 62.17 53.73 89.70 66.78 60.31 52.33 
TTM4-F MAD 2.36 86.99 101.60 148.92 128.08 449.83 327.32 192.82 124.14 
 Score 100 60 50 30 40 −120 −60 10 40 
 AD −140.08 −49.6 −29.7 −6.9 −4.0 −17.3 −11.1 −4.9 −5.9 
 AAD 140.08 76.0 54.3 47.3 48.8 46.0 46.3 49.4 46.5 
AMOEBA03 MAD 188.6 190.3 197.3 197.9 199.1 220.0 207.0 208.7 201.1 
 Score 10 10 10 10 10 −10 
 AD −141.95 −74.6 −73.0 −73.0 −74.4 −68.0 −72.1 −71.6 −72.1 
 AAD 141.95 84.6 80.1 80.4 82.6 76.8 78.6 79.1 82.0 
iAMOEBA MAD 189.7 269.2 264.0 215.8 200.5 239.5 239.8 234.3 200.8 
 Score 10 −30 −30 −10 −10 −10 
 AD −141.7 −73.4 −47.6 −39.2 −37.9 −44.4 −41.0 −38.5 −39.9 
 AAD 141.7 78.3 74.2 66.1 59.7 65.8 63.6 61.1 59.9 
AMOEBA14 MAD 189.8 198.8 215.6 195.9 196.7 229.2 209.7 198.2 195.4 
 Score 10 10 10 10 −10 10 10 
 AD −129.3 −60.2 −20.5 −9.7 −5.5 −23.3 −17.9 −9.0 −3.1 
 AAD 129.3 71.2 72.7 47.5 48.7 54.4 52.5 48.7 51.4 
AMOEBA+ MAD 188.0 191.4 194.7 189.0 186.2 187.0 194.3 186.4 184.3 
 Score 10 10 10 10 10 10 10 10 10 
 AD −141.61 −61.61 −28.82 −21.99 −22.86 −28.90 −23.17 −21.02 −26.31 
 AAD 141.61 77.29 78.01 46.00 47.10 48.36 47.03 50.04 47.06 
AMOEBA+CF MAD 189.94 193.98 190.30 182.44 183.37 196.42 184.40 183.95 184.36 
 Score 10 10 10 10 10 10 10 10 10 
 AD −143.54 −30.73 −4.59 27.54 26.73 −1.40 12.86 20.04 11.77 
 AAD 143.54 99.58 83.23 80.75 78.39 61.38 70.42 73.05 62.48 
MB-UCB MAD 192.23 195.65 203.18 199.47 199.15 202.58 216.20 210.44 206.40 
 Score 10 10 10 10 
 AD −1.39 0.45 0.56 −0.91 −3.64 0.35 −0.85 −3.31 −6.26 
 AAD 1.39 3.93 4.92 10.55 12.96 9.04 10.03 11.58 13.07 
fq-MB-pol MAD 2.37 12.95 19.81 35.92 38.20 28.58 30.00 33.05 45.06 
 Score 100 100 100 90 90 90 90 90 80 
 AD −1.38 −1.84 −1.67 −5.13 −6.00 −2.75 −3.28 −5.44 −7.29 
 AAD 1.38 4.36 4.67 11.58 15.61 7.78 8.88 12.65 16.52 
MB-pol MAD 2.36 11.97 16.53 31.40 46.22 24.10 23.02 38.55 49.37 
 Score 100 100 100 90 80 90 90 90 80 
H2O(H2O)2(H2O)3(H2O)4(H2O)5(H2O)6(H2O)6(H2O)6(H2O)6
isomer 1isomer 1isomer 1isomer 1isomer 1isomer 1isomer 2isomer 3isomer 6
 AD −1.38 −5.39 −0.44 −1.01 1.59 2.25 2.12 1.75 3.80 
TTM2.1-F AAD 1.38 14.57 27.47 60.28 62.91 46.32 50.88 59.79 60.70 
 MAD 2.36 33.66 97.25 213.33 252.42 384.57 352.9 293.6 259.55 
 Score 100 90 60 −20 −90 −70 −40 −20 
 AD −1.38 −31.71 −56.61 −66.37 −65.18 −62.14 −64.37 −64.44 −62.25 
 AAD 1.38 31.71 56.62 66.37 65.22 64.64 64.37 64.44 62.26 
TTM3-F MAD 2.36 145.79 240.43 174.67 180.29 193.36 191.60 193.22 182.82 
 Score 100 30 −20 20 10 10 10 10 10 
 AD −1.38 7.75 7.75 10.41 1.82 7.29 9.32 9.13 −2.88 
 AAD 1.38 30.19 46.30 62.17 53.73 89.70 66.78 60.31 52.33 
TTM4-F MAD 2.36 86.99 101.60 148.92 128.08 449.83 327.32 192.82 124.14 
 Score 100 60 50 30 40 −120 −60 10 40 
 AD −140.08 −49.6 −29.7 −6.9 −4.0 −17.3 −11.1 −4.9 −5.9 
 AAD 140.08 76.0 54.3 47.3 48.8 46.0 46.3 49.4 46.5 
AMOEBA03 MAD 188.6 190.3 197.3 197.9 199.1 220.0 207.0 208.7 201.1 
 Score 10 10 10 10 10 −10 
 AD −141.95 −74.6 −73.0 −73.0 −74.4 −68.0 −72.1 −71.6 −72.1 
 AAD 141.95 84.6 80.1 80.4 82.6 76.8 78.6 79.1 82.0 
iAMOEBA MAD 189.7 269.2 264.0 215.8 200.5 239.5 239.8 234.3 200.8 
 Score 10 −30 −30 −10 −10 −10 
 AD −141.7 −73.4 −47.6 −39.2 −37.9 −44.4 −41.0 −38.5 −39.9 
 AAD 141.7 78.3 74.2 66.1 59.7 65.8 63.6 61.1 59.9 
AMOEBA14 MAD 189.8 198.8 215.6 195.9 196.7 229.2 209.7 198.2 195.4 
 Score 10 10 10 10 −10 10 10 
 AD −129.3 −60.2 −20.5 −9.7 −5.5 −23.3 −17.9 −9.0 −3.1 
 AAD 129.3 71.2 72.7 47.5 48.7 54.4 52.5 48.7 51.4 
AMOEBA+ MAD 188.0 191.4 194.7 189.0 186.2 187.0 194.3 186.4 184.3 
 Score 10 10 10 10 10 10 10 10 10 
 AD −141.61 −61.61 −28.82 −21.99 −22.86 −28.90 −23.17 −21.02 −26.31 
 AAD 141.61 77.29 78.01 46.00 47.10 48.36 47.03 50.04 47.06 
AMOEBA+CF MAD 189.94 193.98 190.30 182.44 183.37 196.42 184.40 183.95 184.36 
 Score 10 10 10 10 10 10 10 10 10 
 AD −143.54 −30.73 −4.59 27.54 26.73 −1.40 12.86 20.04 11.77 
 AAD 143.54 99.58 83.23 80.75 78.39 61.38 70.42 73.05 62.48 
MB-UCB MAD 192.23 195.65 203.18 199.47 199.15 202.58 216.20 210.44 206.40 
 Score 10 10 10 10 
 AD −1.39 0.45 0.56 −0.91 −3.64 0.35 −0.85 −3.31 −6.26 
 AAD 1.39 3.93 4.92 10.55 12.96 9.04 10.03 11.58 13.07 
fq-MB-pol MAD 2.37 12.95 19.81 35.92 38.20 28.58 30.00 33.05 45.06 
 Score 100 100 100 90 90 90 90 90 80 
 AD −1.38 −1.84 −1.67 −5.13 −6.00 −2.75 −3.28 −5.44 −7.29 
 AAD 1.38 4.36 4.67 11.58 15.61 7.78 8.88 12.65 16.52 
MB-pol MAD 2.36 11.97 16.53 31.40 46.22 24.10 23.02 38.55 49.37 
 Score 100 100 100 90 80 90 90 90 80 

After determining the accuracy in reproducing interaction and many-body energies of small water clusters and associated harmonic frequencies, this section investigates the ability of the different models to reproduce the energy landscape of liquid water. Since CCSD(T)/CBS simulations for liquid water under periodic boundary conditions are currently unaffordable, in the following comparisons, MB-pol is used as a reference. Although this choice is arbitrary, it can be justified by considering that several studies have demonstrated that MB-pol correctly reproduces structural, thermodynamic, dynamical, and spectroscopic properties of liquid water across a wide range of temperatures.87,88,102,105–108 For this analysis, we use molecular configurations corresponding to 5 ps of a longer CMD trajectory, which includes nuclear quantum effects,132–136 carried out with MB-pol at 298.15 K for a system containing 256 water molecules.105 

Figure 7 shows that, among the polarizable models, TTM3-F and iAMOEBA more closely follows MB-pol, with average deviations per molecules of −0.15 kcal/mol and 0.31 kcal/mol, respectively, and maximum absolute deviations of 0.26 kcal/mol and 0.45 kcal/mol, respectively. These models are followed by AMOEBA14 that displays average and maximum absolute deviations of 0.48 kcal/mol and 0.57 kcal/mol, respectively. Interestingly, these three models perform rather well in this analysis despite displaying relatively poor agreement with the reference CCSD(T)/CCSD(T)-F12 many-body energies for small clusters. This suggests that, as already discussed in Sec. III B, the apparent better agreement with the MB-pol binding energies for liquid water likely results from error compensation between individual many-body contributions. Relatively larger average deviations per molecule are instead found for TTM2.1-F (−1.50 kcal/mol), TTM4-F (−1.98 kcal/mol), and AMOEBA+ (1.13 kcal/mol), while AMOEBA+CF and MB-UBC provides similar average deviations on the order of 0.6 kcal/mol.

FIG. 7.

Average (AD), average absolute (AAD), and maximum absolute (MAD) deviations of the binding energy per molecule relative to the corresponding MB-pol values along 5 ps of a longer MB-pol CMD trajectory carried out at 298.15 K in Ref. 105.

FIG. 7.

Average (AD), average absolute (AAD), and maximum absolute (MAD) deviations of the binding energy per molecule relative to the corresponding MB-pol values along 5 ps of a longer MB-pol CMD trajectory carried out at 298.15 K in Ref. 105.

Close modal

The comparisons in Fig.7 (and equivalently Table IV) that the use of fixed atomic charges does not prevent fq-MB-pol from quantitatively reproducing the original MB-pol trajectory. Considering that permanent electrostatics is purely pairwise additive in both MB-pol and fq-MB-pol, the agreement shown in Fig. 7 provides further evidence for the ability of the 2-body PIP adopted by these models to correct for inherent deficiencies in classical representations of molecular interactions at short range and demonstrates that (small) differences in the representation of the 1-body charge distributions do not affect the accuracy of explicit many-body models in describing the interactions between water molecules from small clusters to the liquid phase.

TABLE IV.

Average deviation (AD), average absolute deviation (AAD), and maximum absolute deviation (MAD) per molecule calculated with both polarizable and explicit many-body models relative to the MB-pol binding energies per molecule along 5 ps of a MB-pol CMD trajectory carried out at 298.15 K in Ref. 105.

ModelAD (kcal/mol)AAD (kcal/mol)MAD (kcal/mol)
TTM2.1-F −1.50 1.50 1.61 
TTM3-F −0.15 0.15 0.26 
TTM4-F −1.98 1.98 2.16 
AMOEBA03 0.80 0.80 0.94 
iAMOEBA 0.31 0.31 0.45 
AMOEBA14 0.48 0.48 0.57 
AMOEBA+ 1.13 1.13 1.28 
AMOEBA+CF 0.59 0.59 0.67 
MB-UCB 0.64 0.64 0.75 
fq-MB-pol −0.02 0.02 0.03 
ModelAD (kcal/mol)AAD (kcal/mol)MAD (kcal/mol)
TTM2.1-F −1.50 1.50 1.61 
TTM3-F −0.15 0.15 0.26 
TTM4-F −1.98 1.98 2.16 
AMOEBA03 0.80 0.80 0.94 
iAMOEBA 0.31 0.31 0.45 
AMOEBA14 0.48 0.48 0.57 
AMOEBA+ 1.13 1.13 1.28 
AMOEBA+CF 0.59 0.59 0.67 
MB-UCB 0.64 0.64 0.75 
fq-MB-pol −0.02 0.02 0.03 

To provide a more immediate assessment of the ability of each model to reproduce the energy landscape of liquid water, Table V ranks all models based on the corresponding standard deviations shown in Fig. 7. Considering that ∼4.5 molecules are within the first hydration shell as determined from the oxygen–oxygen radial distribution function (RDF) calculated with MB-pol at 298.15 K and 1 atm (which quantitatively agrees with the experimental RDF87), a score of 100 points is assigned if the average absolute deviation per molecule multiplied by 4.5 from the MB-pol binding energies is within 0.593 kcal/mol, which corresponds to kBT at 298.15 K, with kB being Boltzmann’s constant, and one point is then deducted for each additional increment of kBT. It should be noted that, while, rigorously, only free-energy differences should be compared with kBT, our comparisons between models are made for identical liquid configurations, which implies that the entropic contributions are the same for each model, and consequently, energy and free-energy differences are equivalent for these configurations. As already inferred from the analysis of Fig. 7, Table V shows that most polarizable models, with the exception of TTM4-F and AMOEBA+, provide similar performance relative to MB-pol, with scores between 81 and 88. This analysis shows that, already considering only the energetics of the first hydration shell in liquid water, this analysis the deviations of the TTM, AMOEBA, and MB-UCB models from MB-pol are, in all cases, larger than kBT, which suggests that MD simulations of liquid water carried out with these polarizable models would likely sample different regions of the underlying multidimensional energy landscape from those sampled with MB-pol and, consequently, would lead to different molecular structures.

TABLE V.

Average absolute deviations and corresponding scores calculated for the first hydration shell in liquid water using both polarizable and explicit many-body models relative to the MB-pol binding energies along 5 ps of a MB-pol CMD trajectory carried out at 298.15 K in Ref. 105. See the text for details.

Average absolute deviation (kcal/mol)
Modelfor the first hydration shellScore
TTM2.1-F 6.750 89 
TTM3-F 0.675 99 
TTM4-F 8.910 85 
AMOEBA03 3.600 94 
iAMOEBA 1.395 98 
AMOEBA14 2.160 96 
AMOEBA+ 5.085 91 
AMOEBA+CF 2.655 95 
MB-UCB 2.880 95 
fq-MB-pol 0.090 100 
Average absolute deviation (kcal/mol)
Modelfor the first hydration shellScore
TTM2.1-F 6.750 89 
TTM3-F 0.675 99 
TTM4-F 8.910 85 
AMOEBA03 3.600 94 
iAMOEBA 1.395 98 
AMOEBA14 2.160 96 
AMOEBA+ 5.085 91 
AMOEBA+CF 2.655 95 
MB-UCB 2.880 95 
fq-MB-pol 0.090 100 

In this study, we presented an assessment of state-of-the-art polarizable water models with flexible monomers based on a systematic analysis of interaction energies, many-body effects, and harmonic frequencies calculated for molecular systems of increasing size. For water clusters, comparisons were made with reference data calculated at the CCSD(T)/CBS level of theory, the current “gold standard” for chemical accuracy, while MB-pol was used as a reference for liquid water since it has been shown that this model accurately reproduces structural, thermodynamic, dynamical, and spectroscopic properties of water across the phase diagram.

Overall, our analysis indicates that significant progress has been made in the development of polarizable water models, although the path has not always been smooth and straight as demonstrated by early models (e.g., TTM2.1-F and AMOEBA03) often outperforming most recent models (e.g., TTM3-F, TTM4-F, iAMOEBA, and AMOEBA14). In particular, comparisons of interaction and many-body energies provide support for an early observation by Burnham and Xantheas that models able to more closely reproduce high-level ab initio data also appear to provide a more balanced description of the properties of water from the gas to the condensed phase.34 This suggests that some caution should be used in relying on the agreement with experimental data for bulk properties during the parameterization procedure since, given the high dimensionality of the parameter space, the apparent agreement with experiment may often be the result of fortuitous error compensation. In particular, several polarizable models considered in our analysis appear to benefit from significant error compensation between 2-body and 3-body energies, which consequently leads to incorrect representations of hydrogen-bond strengths depending on the system size and structure as demonstrated by relatively large deviations from the reference OH-stretch harmonic frequencies of water cluster, which are associated with all polarizable models.

The errors associated with 2-body and 3-body energies calculated with polarizable models for distorted dimer and trimer configurations indicate that it may be difficult to tune the relatively simple functional forms adopted by popular polarizable models (e.g., TTM and AMOEBA models) to quantitatively account for each individual many-body contribution and, at the same time, be able to correctly describe structural, thermodynamic, and dynamical properties of liquid water. Even the most advanced polarizable models (e.g., AMOEBA+CF and MB-UCB), which provide improved descriptions of the energetics of small water clusters by adopting relatively more sophisticated parameterizations of various physical contributions to the interaction energies, do not display significant improvement in the description of the harmonic frequencies of water cluster and the energy landscape of liquid water. This suggests that, while rigorous and systematic physics-based parameterizations certainly improves the description of water systems in the neighborhoods of minimum-energy structures, the relatively simple functional forms adopted by polarizable models may not allow for a quantitative representation of the entire multidimensional energy landscape.

On the other hand, our analysis indicates that, independent of the system size and specific water arrangements, explicit many-body models (e.g., MB-pol and fq-MB-pol) effectively provide CCSD(T)/CBS accuracy for both minimum-energy and distorted structures. As discussed in previous studies,89 this accuracy is owed to the flexibility of short-range 2-body and 3-body PIPs to correct for inherent deficiencies associated with classical representations of molecular interactions. Despite recent claims about certain polarizable models being comparable in performance to MB-pol,58 the present analysis demonstrates that, as it should be expected after considering the different theoretical frameworks these models are built upon and the nature of the approximations used in the parameterization process, many-body models provide significantly higher accuracy and predictive power than existing polarizable models. Importantly, the analyses presented in this study also demonstrated that, with the development of more advanced representations of low-order n-body interactions, foregoing geometry-dependent atomic charges appear to have negligible impact on reproducing the CCSD(T)/CBS reference values, as illustrated by the accuracy achieved by the fq-MB-pol model, which is able to accurately reproduce interaction and many-body energies and harmonic frequencies of small water clusters as well as the energetics of liquid water despite adopting fixed atomic charges.

Moving forward, we believe that future developments of both polarizable and many-body models should continue in parallel and exploit synergies instead of aiming to compete. In particular, many-body models should serve as a reference for side-by-side comparisons that would provide guidance for systematic parameterizations of new polarizable models. In this regard, Ref. 88 provides extensive reference data calculated with MB-pol for various properties of water across the entire phase diagram, which could be used to develop more balanced parameterizations of polarizable models capable of approaching “the right answers for the right reasons.” On the other hand, while it has been shown that the same theoretical/computational framework adopted by MB-pol can be applied to more complex aqueous systems140–151 and molecular fluids,152,153 future efforts should focus on the development of explicit many-body models for simulations of larger molecular systems in both gas and condensed phases, with the same accuracy exhibited by MB-pol for water.

See the supplementary material for tables with interaction energies and harmonic frequencies of water clusters calculated with both polarizable and many-body explicit models.

The coordinates of all optimized clusters used in the analyses of the interaction and many-body energies as well as harmonic frequencies, along with the coordinates of the CMD trajectory used in the analysis of the energetics of liquid water are available at https://github.com/paesanilab/polarizable_water_models. Any other data generated and analyzed for this study that are not included in this article and its supplementary material are available from the authors upon request.

We thank Colin Egan, Kartik Rallapalli, Kelly Hunter, and Jay Hu for useful discussion. We also thank Pengyu Ren, Jean-Philip Piquemal, and Chengwen Liu for guidance with the AMOEBA+ and AMOEBA+CF models and Akshaya Das and Theresa Head-Gordon for providing all the data for the MB-UCB model. This research was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Science, through Grant No. DE-SC0019490. All graphs were generated using Matplotlib.154 All calculations used resources of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation through Grant No. ACI-1053575, under Allocation No. TG-CHE110009, and the Triton Shared Computing Cluster (TSCC) at the San Diego Supercomputer Center (SDSC).

1.
J. A.
Barker
and
R. O.
Watts
, “
Structure of water: A Monte Carlo calculation
,”
Chem. Phys. Lett.
3
,
144
145
(
1969
).
2.
A.
Rahman
and
F. H.
Stillinger
, “
Molecular dynamics study of liquid water
,”
J. Chem. Phys.
55
,
3336
3359
(
1971
).
3.
B.
Guillot
, “
A reappraisal of what we have learnt during three decades of computer simulations on water
,”
J. Mol. Liq.
101
,
219
260
(
2002
).
4.
C.
Vega
and
J. L. F.
Abascal
, “
Simulating water with rigid non-polarizable models: A general perspective
,”
Phys. Chem. Chem. Phys.
13
,
19663
19688
(
2011
).
5.
I.
Shvab
and
R. J.
Sadus
, “
Atomistic water models: Aqueous thermodynamic properties from ambient to supercritical conditions
,”
Fluid Phase Equilib.
407
,
7
30
(
2016
).
6.
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
).
7.
J. D.
Bernal
and
R. H.
Fowler
, “
A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions
,”
J. Chem. Phys.
1
,
515
548
(
1933
).
8.
F. H.
Stillinger
and
A.
Rahman
, “
Improved simulation of liquid water by molecular dynamics
,”
J. Chem. Phys.
60
,
1545
1557
(
1974
).
9.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
, “
Comparison of simple potential functions for simulating liquid water
,”
J. Chem. Phys.
79
,
926
935
(
1983
).
10.
M. W.
Mahoney
and
W. L.
Jorgensen
, “
A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions
,”
J. Chem. Phys.
112
,
8910
8922
(
2000
).
11.
H. W.
Horn
,
W. C.
Swope
,
J. W.
Pitera
,
J. D.
Madura
,
T. J.
Dick
,
G. L.
Hura
, and
T.
Head-Gordon
, “
Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew
,”
J. Chem. Phys.
120
,
9665
9678
(
2004
).
12.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
, “
The missing term in effective pair potentials
,”
J. Phys. Chem.
91
,
6269
6271
(
1987
).
13.
D.
Hankins
,
J. W.
Moskowitz
, and
F. H.
Stillinger
, “
Water molecule interactions
,”
J. Chem. Phys.
53
,
4544
4554
(
1970
).
14.
A.
Stone
,
The Theory of Intermolecular Interactions
(
Clarendon
,
Oxford
,
1996
).
15.
J.
Caldwell
,
L. X.
Dang
, and
P. A.
Kollman
, “
Implementation of nonadditive intermolecular potentials by use of molecular dynamics: Development of a water–water potential and water–ion cluster interactions
,”
J. Am. Chem. Soc.
112
,
9144
9147
(
1990
).
16.
L. X.
Dang
and
T.-M.
Chang
, “
Molecular dynamics study of water clusters, liquid, and liquid–vapor interface of water with many-body potentials
,”
J. Chem. Phys.
106
,
8149
8159
(
1997
).
17.
C.
Millot
and
A. J.
Stone
, “
Towards an accurate intermolecular potential for water
,”
Mol. Phys.
77
,
439
462
(
1992
).
18.
C.
Millot
,
J.-C.
Soetens
,
M. T. C.
Martins Costa
,
M. P.
Hodges
, and
A. J.
Stone
, “
Revised anisotropic site potentials for the water dimer and calculated properties
,”
J. Phys. Chem. A
102
,
754
770
(
1998
).
19.
B.
Chen
,
J.
Xing
, and
J. I.
Siepmann
, “
Development of polarizable water force fields for phase equilibrium calculations
,”
J. Phys. Chem. B
104
,
2391
2401
(
2000
).
20.
K.-H.
Cho
,
K. T.
No
, and
H. A.
Scheraga
, “
A polarizable force field for water using an artificial neural network
,”
J. Mol. Struct.
641
,
77
91
(
2002
).
21.
H. A.
Stern
,
F.
Rittner
,
B. J.
Berne
, and
R. A.
Friesner
, “
Combined fluctuating charge and polarizable dipole models: Application to a five-site water potential function
,”
J. Chem. Phys.
115
,
2237
2251
(
2001
).
22.
L.
Ojamäe
,
I.
Shavitt
, and
S. J.
Singer
, “
Potential models for simulations of the solvated proton in water
,”
J. Chem. Phys.
109
,
5547
5564
(
1998
).
23.
S.
Iuchi
,
A.
Morita
, and
S.
Kato
, “
Molecular dynamics simulation with the charge response kernel: Vibrational spectra of liquid water and n-methylacetamide in aqueous solution
,”
J. Phys. Chem. B
106
,
3466
3476
(
2002
).
24.
T.
Ishiyama
and
A.
Morita
, “
Analysis of anisotropic local field in sum frequency generation spectroscopy with the charge response kernel water model
,”
J. Chem. Phys.
131
,
244714
(
2009
).
25.
G.
Lamoureux
,
A. D.
MacKerell
, Jr.
, and
B.
Roux
, “
A simple polarizable model of water based on classical Drude oscillators
,”
J. Chem. Phys.
119
,
5185
5197
(
2003
).
26.
W.
Yu
,
P. E. M.
Lopes
,
B.
Roux
, and
A. D.
MacKerell
, Jr.
, “
Six-site polarizable model of water based on the classical Drude oscillator
,”
J. Chem. Phys.
138
,
034508
(
2013
).
27.
H.
Yu
,
T.
Hansson
, and
W. F.
van Gunsteren
, “
Development of a simple, self-consistent polarizable model for liquid water
,”
J. Chem. Phys.
118
,
221
234
(
2003
).
28.
H.
Yu
and
W. F.
van Gunsteren
, “
Charge-on-spring polarizable water models revisited: From water clusters to liquid water to ice
,”
J. Chem. Phys.
121
,
9549
9564
(
2004
).
29.
A.-P. E.
Kunz
and
W. F.
van Gunsteren
, “
Development of a nonlinear classical polarization model for liquid water and aqueous solutions: COS/D
,”
J. Phys. Chem. A
113
,
11570
11579
(
2009
).
30.
P. T.
Kiss
and
A.
Baranyai
, “
A systematic development of a polarizable potential of water
,”
J. Chem. Phys.
138
,
204507
(
2013
).
31.
K. T.
Wikfeldt
,
E. R.
Batista
,
F. D.
Vila
, and
H.
Jónsson
, “
A transferable H2O interaction potential based on a single center multipole expansion: SCME
,”
Phys. Chem. Chem. Phys.
15
,
16542
16556
(
2013
).
32.
V. P.
Sokhan
,
A. P.
Jones
,
F. S.
Cipcigan
,
J.
Crain
, and
G. J.
Martyna
, “
Signature properties of water: Their molecular electronic origins
,”
Proc. Natl. Acad. Sci. U. S. A.
112
,
6341
6346
(
2015
).
33.
C. J.
Burnham
,
J.
Li
,
S. S.
Xantheas
, and
M.
Leslie
, “
The parametrization of a Thole-type all-atom polarizable water model from first principles and its application to the study of water clusters (n = 2–21) and the phonon spectrum of ice Ih
,”
J. Chem. Phys.
110
,
4566
4581
(
1999
).
34.
C. J.
Burnham
and
S. S.
Xantheas
, “
Development of transferable interaction models for water. I. Prominent features of the water dimer potential energy surface
,”
J. Chem. Phys.
116
,
1479
1492
(
2002
).
35.
S. S.
Xantheas
,
C. J.
Burnham
, and
R. J.
Harrison
, “
Development of transferable interaction models for water. II. Accurate energetics of the first few water clusters from first principles
,”
J. Chem. Phys.
116
,
1493
1499
(
2002
).
36.
C. J.
Burnham
and
S. S.
Xantheas
, “
Development of transferable interaction models for water. III. Reparametrization of an all-atom polarizable rigid model (TTM2-R) from first principles
,”
J. Chem. Phys.
116
,
1500
1510
(
2002
).
37.
C. J.
Burnham
and
S. S.
Xantheas
, “
Development of transferable interaction models for water. IV. A flexible, all-atom polarizable potential (TTM2-F) based on geometry dependent charges derived from an ab initio monomer dipole moment surface
,”
J. Chem. Phys.
116
,
5115
5124
(
2002
).
38.
G. S.
Fanourgakis
and
S. S.
Xantheas
, “
The flexible, polarizable, thole-type interaction potential for water (TTM2-F) revisited
,”
J. Phys. Chem. A
110
,
4100
4106
(
2006
).
39.
G. S.
Fanourgakis
and
S. S.
Xantheas
, “
Development of transferable interaction potentials for water. V. Extension of the flexible, polarizable, thole-type model potential (TTM3-F, v. 3.0) to describe the vibrational spectra of water clusters and liquid water
,”
J. Chem. Phys.
128
,
074506
(
2008
).
40.
C. J.
Burnham
,
D. J.
Anick
,
P. K.
Mankoo
, and
G. F.
Reiter
, “
The vibrational proton potential in bulk liquid water and ice
,”
J. Chem. Phys.
128
,
154519
(
2008
).
41.
A.
Defusco
,
D. P.
Schofield
, and
K. D.
Jordan
, “
Comparison of models with distributed polarizable sites for describing water clusters
,”
Mol. Phys.
105
,
2681
2696
(
2007
).
42.
R.
Kumar
,
F.-F.
Wang
,
G. R.
Jenness
, and
K. D.
Jordan
, “
A second generation distributed point polarizable water model
,”
J. Chem. Phys.
132
,
014309
(
2010
).
43.
P. K.
Mankoo
and
T.
Keyes
, “
POLIR: Polarizable, flexible, transferable water potential optimized for IR spectroscopy
,”
J. Chem. Phys.
129
,
034504
(
2008
).
44.
T.
Hasegawa
and
Y.
Tanimura
, “
A polarizable water model for intramolecular and intermolecular vibrational spectroscopies
,”
J. Phys. Chem. B
115
,
5545
5553
(
2011
).
45.
P.
Ren
and
J. W.
Ponder
, “
Polarizable atomic multipole water model for molecular mechanics simulation
,”
J. Phys. Chem. B
107
,
5933
5947
(
2003
).
46.
P.
Ren
and
J. W.
Ponder
, “
Temperature and pressure dependence of the AMOEBA water model
,”
J. Phys. Chem. B
108
,
13427
13437
(
2004
).
47.
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
).
48.
L.-P.
Wang
,
T.
Head-Gordon
,
J. W.
Ponder
,
P.
Ren
,
J. D.
Chodera
,
P. K.
Eastman
,
T. J.
Martinez
, and
V. S.
Pande
, “
Systematic improvement of a classical molecular model of water
,”
J. Phys. Chem. B
117
,
9956
9972
(
2013
).
49.
C.
Liu
,
J.-P.
Piquemal
, and
P.
Ren
, “
AMOEBA+ classical potential for modeling molecular interactions
,”
J. Chem. Theory Comput.
15
,
4122
4139
(
2019
).
50.
C.
Liu
,
J.-P.
Piquemal
, and
P.
Ren
, “
Implementation of geometry dependent charge flux into polarizable AMOEBA+ potential
,”
J. Phys. Chem. Lett.
11
,
419
426
(
2020
).
51.
G. A.
Cisneros
, “
Application of Gaussian electrostatic model (GEM) distributed multipoles in the AMOEBA force field
,”
J. Chem. Theory Comput.
8
,
5072
5080
(
2012
).
52.
R. E.
Duke
,
O. N.
Starovoytov
,
J.-P.
Piquemal
, and
G. A.
Cisneros
, “
GEM*: A molecular electronic density-based force field for molecular dynamics simulations
,”
J. Chem. Theory Comput.
10
,
1361
1365
(
2014
).
53.
H.
Torabifard
,
O. N.
Starovoytov
,
P.
Ren
, and
G. A.
Cisneros
, “
Development of an AMOEBA water model using GEM distributed multipoles
,”
Theor. Chem. Acc.
134
,
101
(
2015
).
54.
R. E.
Duke
and
G. A.
Cisneros
, “
Ewald-based methods for Gaussian integral evaluation: Application to a new parameterization of GEM
,”
J. Mol. Model.
25
,
307
(
2019
).
55.
N.
Gresh
, “
Inter- and intramolecular interactions. Inception and refinements of the SIBFA, molecular mechanics (SMM) procedure, a separable, polarizable methodology grounded on ab initio SCF/MP2 computations. Examples of applications to molecular recognition problems
,”
J. Chim. Phys. Phys.-Chim. Biol.
94
,
1365
1416
(
1997
).
56.
J.-P.
Piquemal
,
N.
Gresh
, and
C.
Giessner-Prettre
, “
Improved formulas for the calculation of the electrostatic contribution to the intermolecular interaction energy from multipolar expansion of the electronic distribution
,”
J. Phys. Chem. A
107
,
10353
10359
(
2003
).
57.
N.
Gresh
,
G. A.
Cisneros
,
T. A.
Darden
, and
J.-P.
Piquemal
, “
Anisotropic, polarizable molecular mechanics studies of inter- and intramolecular interactions and ligand–macromolecule complexes. A bottom-up strategy
,”
J. Chem. Theory Comput.
3
,
1960
1986
(
2007
).
58.
A. K.
Das
,
L.
Urban
,
I.
Leven
,
M.
Loipersberger
,
A.
Aldossary
,
M.
Head-Gordon
, and
T.
Head-Gordon
, “
Development of an advanced force field for water using variational energy decomposition analysis
,”
J. Chem. Theory Comput.
15
,
5001
5013
(
2019
).
59.
B. T.
Thole
, “
Molecular polarizabilities calculated with sa modified dipole interaction
,”
Chem. Phys.
59
,
341
350
(
1981
).
60.
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
).
61.
B.
Hartke
, “
Size-dependent transition from all-surface to interior-molecule structures in pure neutral water clusters
,”
Phys. Chem. Chem. Phys.
5
,
275
284
(
2003
).
62.
S. S.
Xantheas
and
E.
Aprà
, “
The binding energies of the D2d and S4 water octamer isomers: High-level electronic structure and empirical potential results
,”
J. Chem. Phys.
120
,
823
828
(
2004
).
63.
C. J.
Burnham
,
G. F.
Reiter
,
J.
Mayers
,
T.
Abdul-Redah
,
H.
Reichert
, and
H.
Dosch
, “
On the origin of the redshift of the OH stretch in ice Ih: Evidence from the momentum distribution of the protons and the infrared spectral density
,”
Phys. Chem. Chem. Phys.
8
,
3966
3977
(
2006
).
64.
F.
Paesani
,
S.
Iuchi
, and
G. A.
Voth
, “
Quantum effects in liquid water from an ab initio-based polarizable force field
,”
J. Chem. Phys.
127
,
074506
(
2007
).
65.
G. S.
Fanourgakis
and
S. S.
Xantheas
, “
The bend angle of water in ice Ih and liquid water: The significance of implementing the nonlinear monomer dipole moment surface in classical interaction potentials
,”
J. Chem. Phys.
124
,
174504
(
2006
).
66.
S.
Habershon
,
G. S.
Fanourgakis
, and
D. E.
Manolopoulos
, “
Comparison of path integral molecular dynamics methods for the infrared absorption spectrum of liquid water
,”
J. Chem. Phys.
129
,
074501
(
2008
).
67.
F.
Paesani
,
S. S.
Xantheas
, and
G. A.
Voth
, “
Infrared spectroscopy and hydrogen-bond dynamics of liquid water from centroid molecular dynamics with an ab initio-based force field
,”
J. Phys. Chem. B
113
,
13118
13130
(
2009
).
68.
L.-P.
Wang
,
J.
Chen
, and
T.
Van Voorhis
, “
Systematic parametrization of polarizable force fields from quantum chemistry data
,”
J. Chem. Theory Comput.
9
,
452
460
(
2012
).
69.
M. L.
Laury
,
L.-P.
Wang
,
V. S.
Pande
,
T.
Head-Gordon
, and
J. W.
Ponder
, “
Revised parameters for the amoeba polarizable atomic multipole water model
,”
J. Phys. Chem. B
119
,
9423
9437
(
2015
).
70.
A.
Späh
,
H.
Pathak
,
K. H.
Kim
,
F.
Perakis
,
D.
Mariedahl
,
K.
Amann-Winkel
,
J. A.
Sellberg
,
J. H.
Lee
,
S.
Kim
,
J.
Park
,
K. H.
Nam
,
T.
Katayama
, and
A.
Nilsson
, “
Apparent power-law behavior of water’s isothermal compressibility and correlation length upon supercooling
,”
Phys. Chem. Chem. Phys.
21
,
26
31
(
2019
).
71.
O.
Matsuoka
,
E.
Clementi
, and
M.
Yoshimine
, “
CI study of the water dimer potential surface
,”
J. Chem. Phys.
64
,
1351
1361
(
1976
).
72.
G. C.
Lie
and
E.
Clementi
, “
Molecular-dynamics simulation of liquid water with an ab initio flexible water–water interaction potential
,”
Phys. Rev. A
33
,
2679
(
1986
).
73.
M. W.
Evans
,
K.
Refson
,
K. N.
Swamy
,
G. C.
Lie
, and
E.
Clementi
, “
Molecular-dynamics simulation of liquid water with an ab initio flexible water–water interaction potential. II. The effect of internal vibrations on the time correlation functions
,”
Phys. Rev. A
36
,
3935
(
1987
).
74.
U.
Niesar
,
G.
Corongiu
,
E.
Clementi
,
G. R.
Kneller
, and
D. K.
Bhattacharya
, “
Molecular dynamics simulations of liquid water using the NCC ab initio potential
,”
J. Phys. Chem.
94
,
7949
7956
(
1990
).
75.
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
).
76.
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
).
77.
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
).
78.
K.
Szalewicz
,
C.
Leforestier
, and
A.
Van Der Avoird
, “
Towards the complete understanding of water by a first-principles computational approach
,”
Chem. Phys. Lett.
482
,
1
14
(
2009
).
79.
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
).
80.
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
).
81.
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
).
82.
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
).
83.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
, “
A critical assessment of two-body and three-body interactions in water
,”
J. Chem. Theory Comput.
9
,
1103
1114
(
2013
).
84.
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
).
85.
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
).
86.
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. Theory Comput.
10
,
1599
1607
(
2014
).
87.
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
).
88.
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
).
89.
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
).
90.
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
).
91.
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
).
92.
J. D.
Mallory
and
V. A.
Mandelshtam
, “
Diffusion Monte Carlo studies of MB-pol (H2O)2−6 and (D2O)2−6 clusters: Structures and binding energies
,”
J. Chem. Phys.
145
,
064308
(
2016
).
93.
P. E.
Videla
,
P. J.
Rossky
, and
D.
Laria
, “
Communication: Isotopic effects on tunneling motions in the water trimer
,”
J. Chem. Phys.
144
,
061101
(
2016
).
94.
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
).
95.
C. L.
Vaillant
,
D. J.
Wales
, and
S. C.
Althorpe
, “
Tunneling splittings from path-integral molecular dynamics using a Langevin thermostat
,”
J. Chem. Phys.
148
,
234102
(
2018
).
96.
M.
Schmidt
and
P.-N.
Roy
, “
Path integral molecular dynamic simulation of flexible molecular systems in their ground state: Application to the water dimer
,”
J. Chem. Phys.
148
,
124116
(
2018
).
97.
K. P.
Bishop
and
P.-N.
Roy
, “
Quantum mechanical free energy profiles with post-quantization restraints: Binding free energy of the water dimer over a broad range of temperatures
,”
J. Chem. Phys.
148
,
102303
(
2018
).
98.
P. E.
Videla
,
P. J.
Rossky
, and
D.
Laria
, “
Isotopic equilibria in aqueous clusters at low temperatures: Insights from the MB-pol many-body potential
,”
J. Chem. Phys.
148
,
084303
(
2018
).
99.
N. R.
Samala
and
N.
Agmon
, “
Temperature dependence of intramolecular vibrational bands in small water clusters
,”
J. Phys. Chem. B
123
,
9428
9442
(
2019
).
100.
N. R.
Samala
and
N.
Agmon
, “
Thermally induced hydrogen-bond rearrangements in small water clusters and the persistent water tetramer
,”
ACS Omega
4
,
22581
22590
(
2019
).
101.
M. T.
Cvitaš
and
J. O.
Richardson
, “
Quantum tunnelling pathways of the water pentamer
,”
Phys. Chem. Chem. Phys.
22
,
1035
1044
(
2020
).
102.
G. R.
Medders
and
F.
Paesani
, “
Infrared and Raman spectroscopy of liquid water through “first-principles” many-body molecular dynamics
,”
J. Chem. Theory Comput.
11
,
1145
1154
(
2015
).
103.
P.
Gasparotto
,
A. A.
Hassanali
, and
M.
Ceriotti
, “
Probing defects and correlations in the hydrogen-bond network of ab initio water
,”
J. Chem. Theory Comput.
12
,
1953
1964
(
2016
).
104.
S. C.
Straight
and
F.
Paesani
, “
Exploring electrostatic effects on the hydrogen bond network of liquid water through many-body molecular dynamics
,”
J. Phys. Chem. B
120
,
8539
8546
(
2016
).
105.
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
).
106.
K. M.
Hunter
,
F. A.
Shakib
, and
F.
Paesani
, “
Disentangling coupling effects in the infrared spectra of liquid water
,”
J. Phys. Chem. B
122
,
10754
10761
(
2018
).
107.
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
).
108.
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
).
109.
P. E.
Ohno
,
H.-f.
Wang
,
F.
Paesani
,
J. L.
Skinner
, and
F. M.
Geiger
, “
Second-order vibrational lineshapes from the air/water interface
,”
J. Phys. Chem. A
122
,
4457
4464
(
2018
).
110.
S.
Sengupta
,
D. R.
Moberg
,
F.
Paesani
, and
E.
Tyrode
, “
Neat water–vapor interface: Proton continuum and the nonresonant background
,”
J. Phys. Chem. Lett.
9
,
6744
6749
(
2018
).
111.
S.
Sun
,
F.
Tang
,
S.
Imoto
,
D. R.
Moberg
,
T.
Ohto
,
F.
Paesani
,
M.
Bonn
,
E. H.
Backus
, and
Y.
Nagata
, “
Orientational distribution of free OH groups of interfacial water is exponential
,”
Phys. Rev. Lett.
121
,
246101
(
2018
).
112.
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
).
113.
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
).
114.
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
).
115.
Z.
Sun
,
L.
Zheng
,
M.
Chen
,
M. L.
Klein
,
F.
Paesani
, and
X.
Wu
, “
Electron-hole theory of the effect of quantum nuclei on the x-ray absorption spectra of liquid water
,”
Phys. Rev. Lett.
121
,
137401
(
2018
).
116.
A. P.
Gaiduk
,
T. A.
Pham
,
M.
Govoni
,
F.
Paesani
, and
G.
Galli
, “
Electron affinity of liquid water
,”
Nat. Commun.
9
,
247
(
2018
).
117.
M.
Fritz
,
M.
Fernández-Serra
, and
J. M.
Soler
, “
Optimization of an exchange-correlation density functional for water
,”
J. Chem. Phys.
144
,
224101
(
2016
).
118.
T. A.
Halgren
, “
The representation of van der Waals (vdW) interactions in molecular mechanics force fields: Potential form, combination rules, and vdW parameters
,”
J. Am. Chem. Soc.
114
,
7827
7843
(
1992
).
119.
B. J.
Braams
and
J. M.
Bowman
, “
Permutationally invariant potential energy surfaces in high dimensionality
,”
Int. Rev. Phys. Chem.
28
,
577
606
(
2009
).
120.
K. T.
Tang
and
J. P.
Toennies
, “
An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients
,”
J. Chem. Phys.
80
,
3726
3741
(
1984
).
121.
T.
Hastie
,
R.
Tibshirani
, and
J.
Friedman
,
The Elements of Statistical Learning: Data Mining, Inference, and Prediction
(
Springer Science + Business Media
,
2009
).
122.
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
(
Springer
,
2020
), pp.
635
660
.
123.
A. C.
Simmonett
,
A. T. B.
Gilbert
, and
P. M. W.
Gill
, “
An optimal point-charge model for molecular electrostatic potentials
,”
Mol. Phys.
103
,
2789
2793
(
2005
).
124.
N.
Mardirossian
and
M.
Head-Gordon
, “
ωB97M-V: A combinatorially optimized, range-separated hybrid, meta-GGA density functional with VV10 nonlocal correlation
,”
J. Chem. Phys.
144
,
214110
(
2016
).
125.
R. A.
Kendall
,
T. H.
Dunning
, Jr.
, and
R. J.
Harrison
, “
Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions
,”
J. Chem. Phys.
96
,
6796
6806
(
1992
).
126.
D. E.
Woon
and
T. H.
Dunning
, Jr.
, “
Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon
,”
J. Chem. Phys.
98
,
1358
1371
(
1993
).
127.
W.
Smith
and
T. R.
Forester
, “
DL_POLY_2.0: A general-purpose parallel molecular dynamics simulation package
,”
J. Mol. Graphics
14
,
136
141
(
1996
).
128.
J. A.
Rackers
,
Z.
Wang
,
C.
Lu
,
M. L.
Laury
,
L.
Lagardère
,
M. J.
Schnieders
,
J.-P.
Piquemal
,
P.
Ren
, and
J. W.
Ponder
, “
Tinker 8: Software tools for molecular design
,”
J. Chem. Theory Comput.
14
,
5273
5289
(
2018
).
129.
U.
Góra
,
R.
Podeszwa
,
W.
Cencek
, and
K.
Szalewicz
, “
Interaction energies of large clusters from many-body expansion
,”
J. Chem. Phys.
135
,
224102
(
2011
).
130.
J. C.
Howard
,
J. L.
Gray
,
A. J.
Hardwick
,
L. T.
Nguyen
, and
G. S.
Tschumper
, “
Getting down to the fundamentals of hydrogen bonding: Anharmonic vibrational frequencies of (HF)2 and (H2O)2 from ab initio electronic structure computations
,”
J. Chem. Theory Comput.
10
,
5426
5435
(
2014
).
131.
J. C.
Howard
and
G. S.
Tschumper
, “
Benchmark structures and harmonic vibrational frequencies near the CCSD(T) complete basis set limit for small water clusters:(H2O)n=2,3,4,5,6
,”
J. Chem. Theory Comput.
11
,
2126
2136
(
2015
).
132.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. I. Equilibrium properties
,”
J. Chem. Phys.
100
,
5093
5105
(
1994
).
133.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. II. Dynamical properties
,”
J. Chem. Phys.
100
,
5106
5117
(
1994
).
134.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. III. Phase space formalism and analysis of centroid molecular dynamics
,”
J. Chem. Phys.
101
,
6157
6167
(
1994
).
135.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. IV. Algorithms for centroid molecular dynamics
,”
J. Chem. Phys.
101
,
6168
6183
(
1994
).
136.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. V. Quantum instantaneous normal mode theory of liquids
,”
J. Chem. Phys.
101
,
6184
6192
(
1994
).
137.
M. J.
Gillan
,
D.
Alfè
, and
A.
Michaelides
, “
Perspective: How good is DFT for water?
,”
J. Chem. Phys.
144
,
130901
(
2016
).
138.
F.
Paesani
,
S.
Yoo
,
H. J.
Bakker
, and
S. S.
Xantheas
, “
Nuclear quantum effects in the reorientation of water
,”
J. Phys. Chem. Lett.
1
,
2316
2321
(
2010
).
139.
H. J.
Bakker
and
J. L.
Skinner
, “
Vibrational spectroscopy as a probe of structure and dynamics in liquid water
,”
Chem. Rev.
110
,
1498
1517
(
2010
).
140.
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
).
141.
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
).
142.
P.
Bajaj
,
X.-G.
Wang
,
T.
CarringtonJr
, and
F.
Paesani
, “
Vibrational spectra of halide–water dimers: Insights on ion hydration from full-dimensional quantum calculations on many-body potential energy surfaces
,”
J. Chem. Phys.
148
,
102321
(
2017
).
143.
M.
Riera
,
S. E.
Brown
, and
F.
Paesani
, “
Isomeric equilibria, nuclear quantum effects, and vibrational spectra of M+(H2O)n=1−3 clusters, with M = Li, Na, K, Rb, and Cs, through many-body representations
,”
J. Phys. Chem. A
122
,
5811
5821
(
2018
).
144.
F.
Paesani
,
P.
Bajaj
, and
M.
Riera
, “
Chemical accuracy in modeling halide ion hydration from many-body representations
,”
Adv. Phys. X
4
,
1631212
(
2019
).
145.
D.
Zhuang
,
M.
Riera
,
G. K.
Schenter
,
J. L.
Fulton
, and
F.
Paesani
, “
Many-body effects determine the local hydration structure of Cs+ in solution
,”
J. Phys. Chem. Lett.
10
,
406
412
(
2019
).
146.
P.
Bajaj
,
M.
Riera
,
J. K.
Lin
,
Y. E.
Mendoza Montijo
,
J.
Gazca
, and
F.
Paesani
, “
Halide ion microhydration: Structure, energetics, and spectroscopy of small halide–water clusters
,”
J. Phys. Chem. A
123
,
2843
2852
(
2019
).
147.
P.
Bajaj
,
J. O.
Richardson
, and
F.
Paesani
, “
Ion-mediated hydrogen-bond rearrangement through tunnelling in the iodide–dihydrate complex
,”
Nat. Chem.
11
,
367
374
(
2019
).
148.
P.
Bajaj
,
D.
Zhuang
, and
F.
Paesani
, “
Specific ion effects on hydrogen-bond rearrangements in the halide–dihydrate complexes
,”
J. Phys. Chem. Lett.
10
,
2823
2828
(
2019
).
149.
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
).
150.
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
).
151.
Y.
Zhai
,
A.
Caruso
,
S.
Gao
, and
F.
Paesani
, “
Active learning of many-body configuration space: Application to the Cs+–water MB-nrg potential energy function as a case study
,”
J. Chem. Phys.
152
,
144103
(
2020
).
152.
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
).
153.
M.
Riera
,
A.
Hirales
,
R.
Ghosh
, and
F.
Paesani
, “
Data-driven many-body models for CH4/H2O mixtures
,” (unpublished) (
2020
).
154.
J. D.
Hunter
, “
Matplotlib: A 2D graphics environment
,”
Comput. Sci. Eng.
9
,
90
95
(
2007
).

Supplementary Material