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.
I. INTRODUCTION
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
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.
II. MOLECULAR MODELS
A. Polarizable models
1. TTM2.1-F
2. TTM3-F
3. TTM4-F
4. AMOEBA03
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
7. AMOEBA+
Although it was built upon previous versions of the AMOEBA models, AMOEBA+ employs a different representation of 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.
B. Explicit many-body models
1. MB-pol
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.
III. COMPUTATIONAL DETAILS
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.
IV. RESULTS
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.
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.
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.
A. Analysis of 2-body and 3-body energies
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.
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 . | ||
---|---|---|---|---|
Model . | E2 . | E3 . | E2B . | E3B . |
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 . | ||
---|---|---|---|---|
Model . | E2 . | E3 . | E2B . | E3B . |
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.
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.
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.
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.
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.
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.
B. Analysis of many-body energies in small water clusters
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.
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.
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.
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.
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.
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. [Eq. (14)] and [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.
Model . | Cluster . | . | Score . | . | Score . | ΔEdiff . | Score . |
---|---|---|---|---|---|---|---|
(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 | 0 | |
(H2O)4 | 5.87 | 0 | 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 | 0 | |
(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 |
Model . | Cluster . | . | Score . | . | Score . | ΔEdiff . | Score . |
---|---|---|---|---|---|---|---|
(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 | 0 | |
(H2O)4 | 5.87 | 0 | 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 | 0 | |
(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 accompanied by a large value for 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 and . 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.
C. Harmonic frequencies of water 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.
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.
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.
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
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 1 . | isomer 1 . | isomer 1 . | isomer 1 . | isomer 1 . | isomer 1 . | isomer 2 . | isomer 3 . | isomer 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 | 0 | −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 | 0 | 0 | 0 | |
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 | 0 | 0 | −10 | −10 | −10 | 0 | |
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 | 0 | 10 | 10 | −10 | 0 | 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 | 0 | 10 | 10 | 0 | 0 | 0 | 0 | |
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 1 . | isomer 1 . | isomer 1 . | isomer 1 . | isomer 1 . | isomer 1 . | isomer 2 . | isomer 3 . | isomer 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 | 0 | −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 | 0 | 0 | 0 | |
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 | 0 | 0 | −10 | −10 | −10 | 0 | |
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 | 0 | 10 | 10 | −10 | 0 | 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 | 0 | 10 | 10 | 0 | 0 | 0 | 0 | |
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 |
D. Liquid water
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.
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.
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.
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.
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.
Model . | AD (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 |
Model . | AD (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.
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) . | . |
---|---|---|
Model . | for the first hydration shell . | Score . |
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) . | . |
---|---|---|
Model . | for the first hydration shell . | Score . |
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 |
V. CONCLUSION
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.
SUPPLEMENTARY MATERIAL
See the supplementary material for tables with interaction energies and harmonic frequencies of water clusters calculated with both polarizable and many-body explicit models.
DATA AVAILABILITY
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.
ACKNOWLEDGMENTS
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).