How good are polarizable and flexible models for water: Insights from a many-body perspective
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 at 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.
I. INTRODUCTION<!>A. Polarizable models<!>MB-UCB.<!>B. Explicit many-body models<!>III. COMPUTATIONAL DETAILS<!>IV. RESULTS<!>A. Analysis of 2-body and 3-body energies<!>B. Analysis of many-body energies is small water clusters<!>C. Harmonic frequencies of water clusters<!>D. Liquid water<!>V. CONCLUSION<!>VI. SUPPLEMENTARY MATERIAL
<p>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 of representing the interactions between water molecules. [3][4][5][6] Starting from the seminal work by Bernal and Fowler 7 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 ST2, 8 and the subsequent TIPnP [9][10][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.</p><p>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][4][5][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 manybody expansion (MBE) of the total energy, E N , which is expressed as 13</p><p>where N is the number of water molecules and r i collectively denotes the coordinates of the oxygen (O) and hydrogen (H) atoms of the i-th water molecule. Eq. 1 shows that E N 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</p><p>From Eqs. 1 and 2, it follows that the interaction energy, E int , between N water molecules can be expressed as</p><p>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 Kollman and coworkers, 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.  Among the most popular polarizable water models, which takes into account monomer flexibility, are the TTM [34][35][36][37][38][39][40] and AMOEBA [45][46][47][48][49][50] families.</p><p>The TTM models are based on Thole's polarization schemes 59 and were primarily fitted to ab initio data for small water clusters. [33][34][35][36][37][38][39][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)</p><p>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 I h . 36 Building upon the TTM2-R parameterization, the TTM2-F model 37 and its refined TTM2. 1-F version 38 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-F 39 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 lineshapes 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 represent-ing 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 reparametrizations using the ForceBalance algorithm, 68 led to iAMOEBA 48 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 of 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 manybody 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.</p><p>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 Similarly 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 data sets 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][72][73][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 last decade. Since, on average, the sum of 2-body and 3-body energies contributes ∼95% to E int , 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][76][77][78] which was the first explicit many-body water model with rigid monomers capable to reproduce both the properties of the water dimer and the radial distribution functions (RDFs) of liquid water at 300 K, WHBB, [79][80][81][82] which was the first explicit many-body water model with flexible monomers, and HBB2-pol and MB-pol, [83][84][85][86][87] which were the first explicit many-body water models with flexible monomers capable to correctly predict the properties of water across all phases. 88,89 In particular, MB-pol has been shown to quantitatively reproduce the vibrationrotation tunneling spectrum of the water dimer, 85 the energetics, isomeric equilibria, tunneling splittings, and vibrational spsctra 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][108][109][110][111] and the energetics and the infrared and Raman of various ice phases. [112][113][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 water 115 as well as 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 MBpol (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-ofthe-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.</p><p>The article is organized as follows: Section II provides a brief description of both polarizable and explicit many-body models analyzed in this study, while Section III provides computational details specific to the comparisons presented in the following sections. 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 Section V.</p><!><p>The TTM, AMOEBA, and MB-UCB models approximate Eq. 1 as</p><p>where E e f f MB 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 (r i ) 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 E e f f MB into distinct physical contributions to the interaction energy according to</p><p>where E elec , E pol , E rep , E disp , and E CT are individual terms parameterized to represent permanent electrostatics, polarization, repulsion, dispersion, and charge transfer, respectively. Since the TTM, AMOEBA, and MB-UCB models adopts different functional forms to describe each contri-bution to E e f f MB , we briefly summarize here the main features of all models. The interested reader is referred to the original references for specific details.</p><p>TTM2.1-F. In TTM2.1-F, E elec is described by Coulomb interactions between geometrydependent 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. E pol is represented according to a Thole's damping scheme, 59 with inducible point dipoles placed on the O and h atoms. Both E elec and E pol 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 E rep and E disp in a single term representing pairwise additive van der Waals interactions between the oxygen atoms,</p><p>where R i j is the distance between the O atoms of molecules i and j, and D, α, A, B, and C are fitting parameters.</p><p>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</p><p>where, as in TTM2.1-F, R i j is the distance between the O atoms of molecules i and j, and ε, σ , and λ are fitting parameters.</p><p>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,</p><p>where, as in TTM2.1-F and TTM3-F, R i j is the distance between the O atoms of molecules i and j, and A 2l are fitting parameters.</p><p>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 E elec in Eq. 5. These moments were derived from a distributed multipole analysis carried out on ab initio data for the water monomer.</p><p>E pol 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 E pol at short range, no damping is applied to E elec . As the TTM models, AMOEBA03 does not include an explicit charge transfer term, but, differently from the TTM models, includes van der Waals interactions between all pairs of atoms through a sum of pairwise additive buffered 14-7 potentials,</p><p>Here, i and j thus indicate O and H atoms on different water molecules within each dimer, ε i j is the potential well depth of the corresponding 14-7 potential, and ρ i j = R i j /R 0 i j , with R i j being the distance between atoms i and j, and R 0 i j 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.</p><p>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 AMOEBA14. AMOEBA14 69 represents an updated version of the original AMOEBA03 model which was obtained by using the ForceBalance algorithm 68 to determine a refined set of fitting parameters.</p><p>AMOEBA+. Although it was built upon previous versions of the AMOEBA models, AMOEBA+ employs a different representation of E e f f MB in which Thole's damping is applied not only to E pol but also to the interactions between the permanent multipoles of E elec . 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 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 Similarly 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.</p><!><p>Closely related to AMOEBA+ and AMOEBA+CF, the MB-UCB model adopts a similar functional form for E elec , 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 E vdW as AMOEBA+ and AMOEBA+CF, differently from the AMOEBA+ models that account only for pairwise charge transfer, MB-UCB includes a manybody 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 E CT 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.</p><!><p>Both MB-pol and fq-MB-pol approximate Eq. 1 as</p><p>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.</p><p>MB-pol. The 2-body term of MB-pol includes three contributions, 85</p><p>Here, E sr 2B describes short-range 2-body interactions that are represented by a 4 th -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, E elec 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 . E disp ) is a 2-body dispersion term expressed as</p><p>where the sum runs over all pairs of O and H atoms i and j on different water molecules within a dimer, R i j is the distance between atoms i and j, f (δ i j ) is the Tang-Toennies damping function with fitting parameter δ i j , 120 and C 6,i j 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,</p><p>describes 3-body short-range interactions that, as in the analogous 2-body term E sr 3B , are represented by a 4 th -degree PIP 119 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 E sr 2B and E sr 3B 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, E pol 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 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, differently from MB-pol that employs geometry-dependent charges to reproduce the ab initio dipole moment surface of Ref. 60, fq-MBpol employs fixed multipole-derived charges 123 calculated at the ωB97M-V/aug-cc-pVQZ level of theory, [124][125][126] for an isolated water molecule in the corresponding optimized geometry. The multipole-derived charges were then mapped into the M-site (q M ) and the two H atoms (q H ) according to the same scheme adopted by MB-pol and originally introduced in the TTM2.1-F model, resulting in q M = −1.151790e and q H = 0.575895e.</p><!><p>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 3body 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.</p><p>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.</p><!><p>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, the following sections 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][133][134][135][136] simulations of liquid water carried out at ambient conditions with MB-pol in Ref. 105.</p><!><p>The interaction energies (E int 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.</p><p>At the 2-body level, all polarizable models predict interaction energies for the water dimer 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.</p><p>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 TABLE I: Interaction energies (E int ) for the minimum-energy structures of the water dimer and trimer (isomers 1 of (H 2 O) 2 and (H 2 O) 3 in Fig. 1). Also listed are the 2-body and 3-body energies contributing to E int of the water trimer. All energies are in kcal/mol. 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.</p><p>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. Opposite trend is predicted by iAMOEBA that 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 to 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.</p><p>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.</p><p>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 manybody models to describe the multidimensional 2-body and 3-body energy landscapes, Figs.  energies, with RMSDs of ∼0.05 kcal/mol, which are more than one order of magnitude smaller than those associated with the polarizable models.</p><p>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 Section 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 functional forms 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 was 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 Fig. 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 2body 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.</p><p>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.</p><p>Fig. 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 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 exchangerepulsion, 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.</p><!><p>The deviations from the reference many-body energies calculated for (H 2 O) 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.</p><p>Fig. 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.</p><p>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 the ForceBalance algorithm was used to simultaneously reproduce both ab initio and experimental data.</p><p>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 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.</p><p>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 noncyclic 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 owed to MB-UCB using anisotropic atomic polarizabilities that provide better representations of planar structures with concerted water arrangements.</p><p>Finally, independently of the cluster size and structure, MB-pol, which includes explicit 2-body FIG. 5: Comparisons between the interaction energies of the low-lying isomers of (H 2 O) 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. 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.</p><p>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.</p><p>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 (H 2 O) 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 (H 2 O) n cluster is defined as:</p><p>The maximum absolute error in the interaction energy over all isomers of each (H 2 O) n cluster is defined as:</p><p>The relative energy difference between isomer 1 and isomer 2 of each (H 2 O) n cluster is defined as:</p><p>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 10 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, −1.29, 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 10 points are deducted for every additional 0.1 kcal/mol thereafter. If a positive value is predicted for the relative energy difference, 0 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.</p><p>For each cluster, a small value for ∆E MB int, max accompanied by a large value for ∆E MB unsigned, max</p><p>indicates that the model benefits from error compensation between individual many-body contri-butions. 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 ∆E MB int, max and ∆E MB unsigned, max . 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.</p><p>In the case of the tetramer, all models predict the correct sign for ∆E di f f , 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.</p><p>In the case of the hexamer, both the TTM2.1-F and TTM4-F models score 0 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.</p><!><p>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 insights 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 Ref. 130,131.</p><p>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, independently of the cluster size and structure, TTM2.1-F systematically overestimates the harmonic frequencies of both free and hydrogenbonded 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 hydrogenbonded 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 Fig. 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-200 cm −1 , the frequencies of strongly hydrogenbonded OH stretches of the prism, cage, book1, 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 point to deficiencies in the description of strong hydrogen bonds.</p><p>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 ).</p><p>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.</p><p>As shown in Fig. 6, independently 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  MAD 189.94 193.98 190.30 182.44 183.37 196.42 184.40 183.95   MAD 192.23 195.65 203.18 199.47 199.15 202.58 216.20 210.44  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 Section IVB.</p><p>Table III reports  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 Section IVB. 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-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][106][107][108]113,114</p><!><p>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 in 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][106][107][108] For this analysis, we use molecular configurations corresponding to 5 ps of a longer CMD trajectory, which includes nuclear quantum effects,  a system containing 256 water molecules. 105 Since the interaction energies for liquid water predicted by the different models can be on significantly different energy scales, the deviations, ∆E, from the MB-pol reference values are calculated relative to the corresponding average energies,</p><p>Fig. 7 shows that, among the polarizable models, AMOEBA14 more closely follows MB-pol, with a standard deviation of 7.0 kcal/mol, followed by AMOEBA+CF, with a standard deviation of 8.6 kcal/mol. TTM2.1-F, TTM3-F, AMOEBA03, iAMOEBA, and MB-UCB display similar standard deviations, between ∼9.5 kcal/mol and ∼11.0 kcal/mol. Relatively larger standard deviations are instead found for TTM4-F (16.7 kcal/mol) and AMOEBA+ (14.8 kcal/mol). Interestingly, iAMOEBA and AMOEBA14 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 shown for the interaction energies of small water clusters, the apparent better agreement with the MB-pol energies for liquid water is likely due to error compensation. It should be noted that the original TTM (TTM2.1-F) 38 and AMOEBA (AMOEBA03) 45 models, which made extensive use of ab initio data for water clusters during the parameterization process, display deviations from the MB-pol energies that are comparable to those associated with the most recent and sophisticated models such as AMOEBA+CF and MB-UCB.</p><p>The comparisons in Fig. 7 show that the use of fixed atomic charges does not prevent fq-MBpol 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.</p><p>To provide a more immediate assessment of the ability of each model to reproduce the energy landscape of liquid water, Table IV ranks all models based on the corresponding standard deviations shown in Fig. 7. In this analysis, a score of 100 points is given if the standard deviation from the MB-pol energies is within 0.59 kcal/mol (corresponding to k B T at 298.15 K, with k B being Boltzmann's constant) and 1 point is deducted for each additional increment of k B T. As already inferred from the analysis of Fig. 7, Table IV shows that most polarizable models, with</p><!><p>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.</p><p>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.</p><p>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.</p><p>On the other hand, our analysis indicates that, independently 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 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.</p><p>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,</p><p>while it has been shown that the same theoretical/computational framework adopted by MB-pol can be applied to more complex aqueous systems [140][141][142][143][144][145][146][147][148][149][150][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.</p><!><p>Tables with interaction energies and harmonic frequencies of water clusters calculated with both polarizable and many-body explicit models.</p>