Recent work has shown that the many-body expansion of the interaction energy can be used to develop analytical representations of global potential energy surfaces (PESs) for water. In this study, the role of short- and long-range interactions at different orders is investigated by analyzing water potentials that treat the leading terms of the many-body expansion through implicit (i.e., TTM3-F and TTM4-F PESs) and explicit (i.e., WHBB and MB-pol PESs) representations. It is found that explicit short-range representations of 2-body and 3-body interactions along with a physically correct incorporation of short- and long-range contributions are necessary for an accurate representation of the water interactions from the gas to the condensed phase. Similarly, a complete many-body representation of the dipole moment surface is found to be crucial to reproducing the correct intensities of the infrared spectrum of liquid water.

It is known that the global potential energy surface (PES) of a system containing N interacting water molecules can be formally expressed in terms of the many-body expansion of the interaction energy as a sum over n-body terms with 1 ≤ nN,1 

V N ( x 1 , , x N ) = a V ( 1 B ) ( x a ) + a > b V ( 2 B ) ( x a , x b ) + a > b > c V ( 3 B ) ( x a , x b , x c ) + + V ( N B ) ( x 1 , , x N ) .
(1)

Here, xi collectively denotes the coordinates of all atoms in the i-th water molecule, V(1B) is the one-body potential describing the energy required to deform an individual water molecule from its equilibrium geometry, and all higher n-body interactions, V(nB), are defined recursively through

V ( n B ) ( x 1 , , x n ) = V n ( x 1 , , x n ) a V ( 1 B ) ( x a ) a > b V ( 2 B ) ( x a , x b ) a 1 > > a n 1 V ( ( n 1 ) B ) ( x a 1 , x a 2 , , x a n 1 ) .
(2)

The rapid convergence of Equation (1), demonstrated in several studies of water clusters,2–7 suggests that the PES associated with a system containing N water molecules can in principle be represented as a sum of low-order interactions that are amenable to accurate calculations using correlated electronic structure methods [e.g., coupled cluster with single, double, and perturbative triple excitations, CCSD(T), which currently represents the gold standard in quantum chemistry].

The theoretical modeling of many-body effects in water began in the late 1960s when self-consistent field (SCF) calculations were carried out on small water clusters.1,8–13 These studies concluded that nonadditive effects in water are generally non-negligible, with 3B contributions being as large as 10%–15% of the total pair interaction for ring structures. 4B effects were estimated to be, on average, ∼1% of the total pair interaction. The first attempt to derive potential energy functions for water from ab initio data was made by Clementi and co-workers who fitted dimer energies calculated at the Hartree-Fock (HF) level of theory to an analytical representation of the 2B interactions.14–16 Later, configuration interaction (CI) calculations were used to derive the pairwise additive MCY potential,17 to which 3B and 4B terms were subsequently added through a classical description of polarization effects.18 The first water potential including explicit 2B and 3B terms, derived respectively from 4th order Möller-Plesset (MP4) and HF calculations, along with a classical description of higher-body polarization interactions, was reported in Ref. 19. Following these pioneering efforts, the ASP potentials with rigid monomers were derived using intermolecular perturbation theory.20,21

In the 2000s, Xantheas and co-workers introduced the TTM (Thole-type model) water potentials, which for the first time made use of a highly accurate 1B term derived from ab initio calculations.22–25 The latest versions (TTM3-F26 and TTM4-F27) were shown to reproduce the properties of water clusters, liquid water, and ice reasonably well, although some inaccuracies were identified in the calculations of vibrational spectra.28–30Around the same time, 2B and 3B potentials with rigid water monomers were derived by Szalewicz and co-workers from symmetry-adapted perturbation theory (SAPT).31,32 These studies eventually led to the development of the rigid-monomer CC-pol potential, a 25-site model with explicit 2B and 3B terms fitted to CCSD(T)-corrected MP2 dimer energies and SAPT trimer energies, respectively, with higher-body terms being represented through classical polarization.33 CC-pol was shown to accurately reproduce the vibration-rotation tunneling spectrum of the water dimer and to predict the structure of liquid water in reasonable agreement with the experimental data. Within the CC-pol scheme, a refined 2B term with explicit dependence on the monomer flexibility has also been developed.34,35

More recently, the full-dimensional WHBB36 and MB-pol37–39 many-body potentials with flexible monomers were developed. WHBB includes explicit 2B and 3B terms fitted to CCSD(T) and MP2 data, respectively, with long-range many-body effects being represented by the same Thole-type model used in the TTM3-F potential. MB-pol was derived from fits to CCSD(T) energies calculated for both water dimers and trimers in the complete basis set (CBS) limit and includes many-body effects within a modified version of the polarization model employed by the TTM4-F potential. Although WHBB has successfully been applied to static calculations of several water properties, to date, MB-pol is the only many-body potential that has been consistently employed in quantum molecular dynamics (MD) simulations that accurately predicted the properties of water from the gas to the condensed phase.

While the calculations with many-body potentials reported in the literature demonstrate that Equation (1) can be effectively used to develop accurate molecular-level representations of the water interactions, the role of short- and long-range contributions at different orders has not been fully investigated. Here, we seek to address this issue by analyzing water potentials that treat the leading terms of the many-body expansion through implicit (e.g., TTM3-F and TTM4-F) and explicit (e.g., WHBB and MB-pol) representations, with a specific focus on how these terms are defined and incorporated as a function of the intermolecular separations. The article is organized as follows: The four potentials are briefly described in Section II, an analysis of the energetics as well as of the vibrational frequencies and intensities for water systems ranging from the gas-phase dimer to the liquid phase is reported in Section III, and the conclusions are given in Section IV.

The four potentials studied here were chosen because they treat the intramolecular distortions on equal footing, i.e., through the 1B PES developed by Partridge and Schwenke.40 In this way, the ability of the models to predict water properties is a direct reflection of the treatment of the intermolecular interactions. In addition, WHBB and MB-pol effectively share the same induction scheme employed by TTM3-F and TTM4-F, respectively, which thus allows for a systematic investigation of the effects of explicit low-order terms of the many-body expansion.

The dominant terms of the many-body water interactions are the 2B and 3B terms, which can be conveniently split into short- and long-range contributions, V ( 2 B , 3 B ) = V short ( 2 B , 3 B ) + V long ( 2 B , 3 B ) . At the two-body level, long-range interactions are dominated by electrostatics and dispersion, while exchange-repulsion and charge-transfer become increasingly important at short-range. 3B interactions in water, on the other hand, arise primarily from induced interactions at long-range and exchange-repulsion when the monomer electron densities overlap.41 Both TTM3-F and TTM4-F include permanent and induced electrostatics as well as dispersion and repulsion at the 2B level, while all higher order terms are represented through many-body induction.26,27

Although WHBB and MB-pol seek to exploit the range separation of the low-order water interactions, the two potentials are intrinsically different both in philosophy and by construction. In WHBB, the short-range 2B term is described by a 7th-degree permutationally invariant polynomial that smoothly transitions for an oxygen-oxygen separation between 6.5 and 7.5 Å into the long-range component described by the 2B TTM3-F potential. The 3B term only includes a short-range component, which is represented by either a 5th- (WHBB5) or a 6th-degree (WHBB6) permutationally invariant polynomial that dies off as the largest oxygen-oxygen separation between two molecules of the trimer approaches 8.0 Å. All higher-order interactions (n ≥ 4) are described by the polarization model employed in the TTM3-F potential. By construction, WHBB thus employs a strict separation at the 2B and 3B levels between short- and long-range interactions that are described in a completely independent way and are essentially disentangled from all higher-order contributions. It should be noted that a simplified version of WHBB including only 1B, 2B, and 3B interactions with shorter cutoffs was also used,42 which, however, cannot be implemented in a straightforward way in standard software for molecular simulations in periodic boundary conditions. For this reason, only calculations with the full WHBB5 implementation of Ref. 36 are reported in the following analysis.

On the other hand, MB-pol can be viewed as a classical polarizable model supplemented by short-range 2B and 3B terms that effectively represent quantum-mechanical interactions arising from the overlap of the monomer electron densities. Specifically, at all separations, V(2B) includes (damped) dispersion forces derived from ab initio computed asymptotic expansions of the dispersion energy along with electrostatic contributions due to the interactions between the molecular permanent and induced moments. At short-range, V(2B) is supplemented by a 4th-degree permutationally invariant polynomial that smoothly switches to zero as the oxygen-oxygen separation in the dimer approaches 6.5 Å. Similarly, the MB-pol 3B term, V(3B), includes a 3B polarization term at all separations, which is supplemented by a short-range 4th-degree permutationally invariant polynomial that effectively corrects for the deficiencies of a purely classical representation of the 3B interactions in regions where the electron densities of the three monomers overlap. This short-range 3B potential is smoothly switched off once the oxygen-oxygen separation between any water molecule and the other two water molecules of a trimer reaches a value of 4.5 Å. In MB-pol, all induced interactions are described through many-body polarization using a slightly modified version of TTM4-F.37 MB-pol thus contains many-body effects at all monomer separations as well as at all orders, in an explicit way up to the third order and in a mean-field fashion at all higher orders. It should be noted that alternative functions to the permutationally invariant polynomials (e.g., GAP43) have been suggested and successfully employed in modeling short-range many-body interactions.44 

The effects of the different representations of the many-body interactions employed by TTM3-F, TTM4-F, WHBB, and MB-pol are investigated through the analysis of the energetics of water systems ranging from the gas-phase dimer to the liquid phase. Figure 1(a) shows a comparison of the dimer binding energies calculated for the global minimum configuration of each PES. Although all potentials predict binding energies within ∼0.2 kcal/mol of the CCSD(T)/CBS reference value (independently of whether the dimer geometry is optimized within each PES or the CCSD(T)/CBS dimer geometry is used in the calculations),45 the analysis of the many-body water interactions reported in Ref. 46 shows that the TTM3-F and TTM4-F interaction energies deviate significantly from the CCSD(T) reference data when dimer configurations different from the corresponding equilibrium geometries are considered. Specifically, root mean square deviations (RMSDs) of 0.73 kcal/mol and 0.44 kcal/mol were obtained in Ref. 46 for TTM3-F and TTM4-F, respectively, from comparisons with ∼1400 CCSD(T)/aug-cc-pVTZ 2B energies corrected for the basis set superposition error. It was instead shown that both WHBB and MB-pol provide a highly accurate representation of the whole multidimensional dimer PES as demonstrated by RMSDs of 0.15 kcal/mol and 0.05 kcal/mol, respectively, obtained from comparisons with more than 40 000 CCSD(T)/CBS 2B energies.37 

FIG. 1.

(a) Binding energies for the global minimum configuration of the water dimer (in kcal/mol) obtained with the TTM3-F, TTM4-F, WHBB, and MB-pol potentials in comparison to the reference CCSD(T) value from Ref. 45 (filled symbols). Also shown as open symbols are the binding energies calculated using the reference CCSD(T) dimer geometry. (b) Comparison of the relative binding energies (ΔE) of the water hexamer isomers calculated with the TTM3-F, TTM4-F, WHBB, and MB-pol potentials using the ab initio geometries of Ref. 47. Also shown as a reference are the ab initio values from Ref. 47. (c) Mean absolute difference in the total energy between QMC and DFT with various exchange-correlation functionals and the TTM3-F, TTM4-F, WHBB, and MB-pol potentials calculated for configurations (in periodic boundary conditions) extracted from path-integral molecular dynamics simulations of water performed with the vdW-DF and vdW-DF2 functionals. The statistical errors in the mean absolute differences arising from the QMC calculations are less than 0.006 mHa/molecule. See Ref. 48 for specific details on the QMC and DFT calculations.

FIG. 1.

(a) Binding energies for the global minimum configuration of the water dimer (in kcal/mol) obtained with the TTM3-F, TTM4-F, WHBB, and MB-pol potentials in comparison to the reference CCSD(T) value from Ref. 45 (filled symbols). Also shown as open symbols are the binding energies calculated using the reference CCSD(T) dimer geometry. (b) Comparison of the relative binding energies (ΔE) of the water hexamer isomers calculated with the TTM3-F, TTM4-F, WHBB, and MB-pol potentials using the ab initio geometries of Ref. 47. Also shown as a reference are the ab initio values from Ref. 47. (c) Mean absolute difference in the total energy between QMC and DFT with various exchange-correlation functionals and the TTM3-F, TTM4-F, WHBB, and MB-pol potentials calculated for configurations (in periodic boundary conditions) extracted from path-integral molecular dynamics simulations of water performed with the vdW-DF and vdW-DF2 functionals. The statistical errors in the mean absolute differences arising from the QMC calculations are less than 0.006 mHa/molecule. See Ref. 48 for specific details on the QMC and DFT calculations.

Close modal

Due to the different treatment of many-body effects, the relative accuracies of the TTM3-F, TTM4-F, WHBB, and MB-pol potentials are expected to differ as the number of molecules in the system of interest increases. This is shown in Figure 1(b), where the relative energies of low-lying hexamer isomers calculated with the four potentials are compared with the corresponding ab initio values.47 The hexamer isomers hold a special place in the experimental and theoretical studies of water because they are the smallest clusters in which the lowest energy configurations correspond to fully three-dimensional structures. For this reason, the hexamer serves as a prototypical model for the hydrogen-bond networks observed in condensed phases. While TTM3-F, WHBB, and MB-pol (qualitatively) reproduce the isomer energy ordering, large deviations are found between the TTM4-F predictions and the reference ab initio values. Similar results are also obtained for other small water clusters, including the tetramer and pentamer isomers, shown in Figures S1 and S2 of the supplemental material.49 Since it was shown that the Thole-type polarization scheme employed by the TTM4-F potential correctly captures higher-order water interactions,38,46 the low accuracy shown by TTM4-F in describing the relative energies of the hexamer isomers is possibly due to intrinsic deficiencies in the 2B term. In this context, the performance of the TTM3-F potential, which provides the best agreement with the reference ab initio values, is somewhat surprising given the large deviations already seen at both the 2B and 3B levels,46 and may suggest some fortuitous cancellation of errors. Based on this, it is interesting to note that the isomer energies predicted by WHBB, which employs the same electrostatic model as TTM3-F, deviate from the reference data by as much as 2.0 kcal/mol. The origin of these deviations may result from inaccuracies in the short-range 2B and 3B WHBB polynomials, incompatibility between the effective many-body representation encoded in the TTM3-F electrostatic model with the short-range WHBB polynomials, or both. On the other hand, MB-pol, which employs a slightly modified version of the TTM4-F electrostatic model, predicts isomer energies in good agreement with the reference data, suggesting that the MB-pol short-range 2B and 3B terms accurately correct for the deficiencies that are intrinsic to the purely classical representation of these interactions encoded in the TTM4-F electrostatic model.

To provide more quantitative insights into the different performance of the four potentials in describing the relative energies of the hexamer isomers, a many-body decomposition of the corresponding interaction energies (Table I) was carried out at the CCSD(T)-F12b level of theory50,51 with the VTZ-F12 basis set,52 which was shown to effectively provide the same accuracy as CCSD(T)/CBS calculations at a reduced computational cost.53,54 As has been well established,2–7 the interaction energies are dominated by the 2B and 3B terms, with 4B and higher terms making largest (but still relatively small) contributions for cyclic structures.55 It is interesting to note that, while the hexamer binding energies predicted by TTM3-F are in excellent agreement with the reference calculations (Figure 1(b)), the results of Table I demonstrate that this agreement is clearly achieved by a fortuitous cancellation of errors. For example, in the prism isomer, the TTM3-F 2B interaction is too attractive by ∼3.2 kcal/mol while the 3B interaction is too repulsive by ∼3.8 kcal/mol, resulting in a total interaction energy that differs from the reference value by merely 0.79 kcal/mol. In contrast, TTM4-F, which accurately represents the electrostatic interactions (as shown by 3B contributions that lie within 1 kcal/mol of the reference values for all low-lying hexamer isomers), is unable to correctly describe the 2B interactions, which have significant contributions from exchange/repulsion and charge transfer. Thus, without the benefit of fortuitous cancellation of errors (that can, to some extent, be encoded into models through empirical parametrization), potentials with accurate electrostatics but without a comparably accurate treatment of short-range quantum mechanical effects may be expected to have relatively poor performance. Having been explicitly derived from many-body expansions of the interaction energies obtained from correlated electronic structure data, both WHBB and MB-pol closely reproduce the CCSD(T)-F12 reference many-body terms. However, non-negligible discrepancies in the 3B and 4B terms result in WHBB being appreciably less accurate than MB-pol at reproducing the hexamer interaction energies.

TABLE I.

Many-body decomposition of the interaction energy of low-lying hexamer isomers. All energies are in kcal/mol. The CCSD(T)-F12b50,51 calculations were performed with the VTZ-F12 basis set52 and corrected for basis set superposition error through the site-site function counterpoise method56 using the Molpro quantum chemistry package.57 

(a) Prism (b) Cage
CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol
2B  −38.65  −41.83  −34.96  −38.59  −38.96  2B  −38.20  −42.01  −36.43  −37.99  −38.48 
3B  −8.78  −5.02  −8.41  −9.19  −8.74  3B  −9.04  −4.77  −8.92  −9.45  −8.93 
4B  −0.66  −0.42  −0.48  −0.42  −0.52  4B  −0.53  −0.26  −0.43  −0.26  −0.47 
5B  0.06  0.03  0.06  0.03  0.05  5B  0.01  0.01  0.03  0.01  0.03 
6B  0.00  0.00  0.00  0.00  0.00  6B  0.00  0.00  0.00  0.00  0.00 
Total  −48.03  −47.24  −43.79  −48.17  −48.17  Total  −47.77  −47.02  −45.74  −47.69  −47.85 
(c) Book-1  (d) Book-2 
  CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol    CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol 
2B  −35.76  −40.22  −35.21  −35.53  −35.80  2B  −35.86  −40.35  −35.61  −35.65  −35.93 
3B  −10.46  −5.65  −10.05  −9.99  −10.26  3B  −10.18  −5.28  −9.75  −9.71  −10.03 
4B  −1.08  −0.56  −0.89  −0.56  −0.92  4B  −1.00  −0.49  −0.81  −0.49  −0.85 
5B  −0.04  −0.02  −0.01  −0.02  −0.01  5B  −0.02  −0.01  0.00  −0.01  0.00 
6B  0.00  0.00  0.00  0.00  0.00  6B  0.00  0.00  0.00  0.00  0.00 
Total  −47.34  −46.46  −46.16  −46.10  −47.00  Total  −47.06  −46.14  −46.17  −45.86  −46.81 
(e) Bag  (f) Cyclic-chair 
  CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol    CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol 
2B  −35.01  −39.71  −34.99  −34.89  −35.24  2B  −32.42  −37.85  −32.43  −32.06  −32.44 
3B  −10.43  −5.42  −9.60  −9.86  −10.15  3B  −11.87  −6.66  −11.13  −10.71  −11.64 
4B  −1.16  −0.59  −0.86  −0.59  −0.90  4B  −1.78  −1.03  −1.43  −1.03  −1.45 
5B  −0.02  −0.02  0.00  −0.02  −0.01  5B  −0.19  −0.12  −0.10  −0.12  −0.10 
6B  0.01  0.00  0.00  0.00  0.00  6B  −0.01  −0.01  0.00  −0.01  0.00 
Total  −46.61  −45.74  −45.45  −45.36  −46.30  Total  −46.27  −45.66  −45.09  −43.92  −45.63 
(g) Cyclic-boat-1  (h) Cyclic-boat-2 
  CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol    CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol 
2B  −32.02  −37.54  −32.33  −31.79  −32.02  2B  −31.96  −37.55  −32.45  −31.66  −32.02 
3B  −11.44  −6.22  −10.90  −10.35  −11.30  3B  −11.43  −6.35  −10.93  −10.45  −11.29 
4B  −1.63  −0.89  −1.33  −0.89  −1.35  4B  −1.61  −0.90  −1.33  −0.90  −1.35 
5B  −0.16  −0.09  −0.09  −0.09  −0.09  5B  −0.16  −0.09  −0.09  −0.09  −0.09 
6B  −0.01  −0.01  0.00  −0.01  0.00  6B  −0.01  −0.01  0.00  −0.01  0.00 
Total  −45.26  −44.75  −44.65  −43.13  −44.76  Total  −45.16  −44.91  −44.81  −43.11  −44.76 
(a) Prism (b) Cage
CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol
2B  −38.65  −41.83  −34.96  −38.59  −38.96  2B  −38.20  −42.01  −36.43  −37.99  −38.48 
3B  −8.78  −5.02  −8.41  −9.19  −8.74  3B  −9.04  −4.77  −8.92  −9.45  −8.93 
4B  −0.66  −0.42  −0.48  −0.42  −0.52  4B  −0.53  −0.26  −0.43  −0.26  −0.47 
5B  0.06  0.03  0.06  0.03  0.05  5B  0.01  0.01  0.03  0.01  0.03 
6B  0.00  0.00  0.00  0.00  0.00  6B  0.00  0.00  0.00  0.00  0.00 
Total  −48.03  −47.24  −43.79  −48.17  −48.17  Total  −47.77  −47.02  −45.74  −47.69  −47.85 
(c) Book-1  (d) Book-2 
  CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol    CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol 
2B  −35.76  −40.22  −35.21  −35.53  −35.80  2B  −35.86  −40.35  −35.61  −35.65  −35.93 
3B  −10.46  −5.65  −10.05  −9.99  −10.26  3B  −10.18  −5.28  −9.75  −9.71  −10.03 
4B  −1.08  −0.56  −0.89  −0.56  −0.92  4B  −1.00  −0.49  −0.81  −0.49  −0.85 
5B  −0.04  −0.02  −0.01  −0.02  −0.01  5B  −0.02  −0.01  0.00  −0.01  0.00 
6B  0.00  0.00  0.00  0.00  0.00  6B  0.00  0.00  0.00  0.00  0.00 
Total  −47.34  −46.46  −46.16  −46.10  −47.00  Total  −47.06  −46.14  −46.17  −45.86  −46.81 
(e) Bag  (f) Cyclic-chair 
  CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol    CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol 
2B  −35.01  −39.71  −34.99  −34.89  −35.24  2B  −32.42  −37.85  −32.43  −32.06  −32.44 
3B  −10.43  −5.42  −9.60  −9.86  −10.15  3B  −11.87  −6.66  −11.13  −10.71  −11.64 
4B  −1.16  −0.59  −0.86  −0.59  −0.90  4B  −1.78  −1.03  −1.43  −1.03  −1.45 
5B  −0.02  −0.02  0.00  −0.02  −0.01  5B  −0.19  −0.12  −0.10  −0.12  −0.10 
6B  0.01  0.00  0.00  0.00  0.00  6B  −0.01  −0.01  0.00  −0.01  0.00 
Total  −46.61  −45.74  −45.45  −45.36  −46.30  Total  −46.27  −45.66  −45.09  −43.92  −45.63 
(g) Cyclic-boat-1  (h) Cyclic-boat-2 
  CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol    CCSD(T)-F12  TTM3-F  TTM4-F  WHBB5  MB-pol 
2B  −32.02  −37.54  −32.33  −31.79  −32.02  2B  −31.96  −37.55  −32.45  −31.66  −32.02 
3B  −11.44  −6.22  −10.90  −10.35  −11.30  3B  −11.43  −6.35  −10.93  −10.45  −11.29 
4B  −1.63  −0.89  −1.33  −0.89  −1.35  4B  −1.61  −0.90  −1.33  −0.90  −1.35 
5B  −0.16  −0.09  −0.09  −0.09  −0.09  5B  −0.16  −0.09  −0.09  −0.09  −0.09 
6B  −0.01  −0.01  0.00  −0.01  0.00  6B  −0.01  −0.01  0.00  −0.01  0.00 
Total  −45.26  −44.75  −44.65  −43.13  −44.76  Total  −45.16  −44.91  −44.81  −43.11  −44.76 

The differences between the four potentials become more evident when their accuracy is assessed against benchmark quantum Monte Carlo (QMC) interaction energies calculated for liquid water configurations extracted from path-integral molecular dynamics simulations of the vdW-DF and vdW-DF2 functionals.48 QMC has been shown to be a reliable benchmark in the study of small water clusters, producing relative energies with an accuracy comparable to that of CCSD(T).58,59 As a measure of accuracy, the mean absolute deviation (MAD) between the energies Ei obtained with each of the four potentials and the reference QMC energies EQMC was calculated as

MAD = 1 N c i = 1 N c E i E i QMC E E QMC .
(3)

Here, Nc is the number of water configurations and E E QMC = ( 1 / N c ) i = 1 N c ( E i E i QMC ) is the average energy difference for all configurations, effectively aligning the zero of energy with the reference QMC data. Figure 1(c) shows that the mean absolute deviation associated with WHBB is ∼4 times larger than the corresponding values obtained with both TTM3-F and TTM4-F, and ∼15 times larger than the MB-pol value. Interestingly, MB-pol also achieves better accuracy than all DFT models analyzed in Ref. 48. The TTM3-F and TTM4-F results indicate that effective representations of the many-body interactions through classical polarization models can provide a reasonable description of the energetics associated with the liquid phase (i.e., comparable with that provided by most of the DFT models currently used in water simulations), albeit through empirical parameterization. It should be noted that, as shown by the results reported in Table I and the many-body analysis of different DFT models carried out in Ref. 60, cancellation of errors between different terms of the many-body expansion of the interaction energies may also affect the energetics of bulk systems described by both DFT and classical polarizable models. On the other hand, the increased accuracy of MB-pol relative to TTM4-F indicates that short-range many-body interactions beyond a purely classical electrostatic representation are necessary for a correct, molecular-level description of the liquid phase. The significantly different accuracies associated with the WHBB/TTM3-F and MB-pol/TTM4-F potentials also suggest that the correct integration of explicit short-range and effective long-range many-body interactions is critical for ensuring the transferability of the potential across different phases.

The ability of TTM3-F, TTM4-F, WHBB, and MB-pol to represent the global multidimensional PESs of water systems with an increasing number of molecules is directly reflected in the accuracy with which the four potentials predict the associated vibrational frequencies. Harmonic frequencies for the water hexamer have been reported at the CCSD(T)/aug-cc-pVDZ level.62 More recently, ab initio reference harmonic frequencies for small water clusters have been obtained through two-body:many-body CCSD(T):MP2 calculations using the cc-pVQZ and aug-cc-pVQZ basis sets for H and O atoms, respectively, which enables a more complete treatment of the basis set through the hybrid two-body:many-body approach.61 As shown in Figure 2, the WHBB and MB-pol PESs result in harmonic frequencies for the water dimer which deviate by less than 10 cm−1 from the ab initio frequencies reported in Ref. 61. These results are consistent with the analysis of the dimer vibration-rotation tunneling spectra of WHBB and MB-pol, both of which exhibit nearly perfect agreement with the corresponding experimental data.37 In contrast, relatively large deviations are seen for vibrational frequencies of both the TTM3-F and TTM4-F potentials. The different performances of the four potentials in predicting the relative energies of the hexamer isomers clearly lead to different levels of agreement with the ab initio harmonic frequencies. Independently of the isomer, the MB-pol values consistently lie within 50 cm−1 of the reference values while the WHBB harmonic frequencies can deviate, in some cases, by more than 100 cm−1. Substantially larger deviations, up to ∼200 cm−1, are instead obtained with both TTM3-F and TTM4-F, reinforcing the notion that purely classical representations of the many-body water interactions are likely not sufficient to fully reproduce the complexity of the underlying multidimensional PESs. Interestingly, all potentials predict somewhat larger deviations for the cyclic-chair isomer of the water hexamer, supporting early observations that many-body effects in water are relatively more important for ring structures.13 

FIG. 2.

Deviations from the 2-body:many-body CCSD(T):MP2 harmonic vibrational frequencies61 of (H2O)n isomers with n = 2 and n = 6 calculated using TTM3-F, TTM4-F, WHBB5, and MB-pol.

FIG. 2.

Deviations from the 2-body:many-body CCSD(T):MP2 harmonic vibrational frequencies61 of (H2O)n isomers with n = 2 and n = 6 calculated using TTM3-F, TTM4-F, WHBB5, and MB-pol.

Close modal

While the differences in the PESs clearly affect the underlying vibrational structure, as shown in Figure 2, the infrared activity of those vibrations is ultimately dictated by the associated dipole moment surfaces (DMSs). Many-body representations of higher-order electric properties for molecular systems were characterized beginning in the early 1980s.64 The first analysis of many-body effects on the dipole moment of polar molecules was reported by Skwara et al.65 As noted in 1996 by Schwenke,66 molecular dipole moments can efficiently be represented in terms of geometry-dependent effective charges multiplied by their Cartesian positions. In line with these ideas, the first representation of the 2B dipole of water was reported as part of the WHBB suite.36 A refined version of this 1B + 2B DMS42 has been used to calculate the transition dipole moments which were then used to modulate the WHBB5 frequency distributions calculated within the local monomer approximation from configurations of liquid water extracted from MD simulations with the rigid E3B.67 Neglecting both 3B and higher-order contributions to the dipole moments and dynamical effects (e.g., motional narrowing), good agreement in the relative intensities was obtained with the measured IR spectrum of liquid water.42 Recently, a complete many-body representation of the DMS of water (MB-μ) was reported,68 consisting of the one-body DMS of Lodi et al.,69 an explicit two-body term, and a slightly modified version of TTM4-F polarization for long-range two-body and all higher-order many-body induced dipole moments. The particular functional form of MB-μ was derived from a systematic analysis of the many-body convergence of the electrostatic properties of water.46 

To investigate the effects of many-body dipole moments on the IR spectrum of liquid water, many-body molecular dynamics (MB-MD) simulations,68 within the centroid molecular dynamics (CMD) formalism, were carried out with the MB-pol PES in combination with the MB-μ DMS. As shown in Figure 3, while the overall shape of the IR spectrum is captured in the 1B + 2B representation of the DMS, 3B and higher-order contributions to the dipole moment significantly affect the IR intensities. 2B contributions to the DMS lower the bend intensity while they contribute to ∼70% of the intensity of the OH stretch band. Importantly, while 2B contributions are critical to appearance of the shoulder at ∼180 cm−1 corresponding to the hydrogen-bonding stretch, the 1B contributions are non-negligible. This indicates that some rotational motion is also involved in the hydrogen-bonding stretch. Nevertheless, the absolute intensities of the librational, bending, and stretching bands are only recovered when all many-body effects are included in the calculation of the dipole moment.30 These results support early observations based on molecular dynamics simulations with the SPC potential supplemented with a 3B dipole-induced dipole term showing that the calculated far infrared spectrum of liquid water was in better agreement with experiment relative to results obtained including only a 2B description of the dipole moment.71 

FIG. 3.

Decomposition of the IR spectrum obtained from CMD trajectories with the MB-pol PES in terms of the many-body components of the MB-μ DMS. 1B-Dip indicates that the one-body (gas-phase monomer) dipoles were used to calculate the dipole of the molecules sampled along the MB-pol CMD trajectories, from which the IR spectrum was calculated. (1B + 2B)-Dip indicates that short-ranged two-body dipoles were used in addition to the one-body dipoles. MB-Dip is the full MB-μ many-body dipole. The spectra were smoothed to facilitate the comparison between the line shapes obtained using the different approximations. The experimental IR spectrum was taken from Ref. 63.

FIG. 3.

Decomposition of the IR spectrum obtained from CMD trajectories with the MB-pol PES in terms of the many-body components of the MB-μ DMS. 1B-Dip indicates that the one-body (gas-phase monomer) dipoles were used to calculate the dipole of the molecules sampled along the MB-pol CMD trajectories, from which the IR spectrum was calculated. (1B + 2B)-Dip indicates that short-ranged two-body dipoles were used in addition to the one-body dipoles. MB-Dip is the full MB-μ many-body dipole. The spectra were smoothed to facilitate the comparison between the line shapes obtained using the different approximations. The experimental IR spectrum was taken from Ref. 63.

Close modal

While using a highly accurate PES is a prerequisite for a physically correct description of the water properties at the molecular level, an appropriate balance between accuracy and efficiency is often critical when deciding which potential to employ in computer simulations since the computational cost directly affects the ability to calculate statistically converged quantities. As shown in Table II, a performance analysis carried out on a single Intel Xeon E5-2640 processor for a system consisting of 256 water molecules in periodic boundary conditions indicates that MB-pol is able to achieve a high level of accuracy at a cost of 47 × that of the fixed point charge q-TIP4P/F model72 and ∼5.5 × that of the TTM4-F potential. WHBB, on the other hand, is ∼29 000 more expensive than q-TIP4P/F. A reference implementation of MB-pol is available as an independent plugin73 for the OpenMM toolkit for molecular simulations.74 

TABLE II.

Cost per molecular dynamics step for different PESs relative to q-TIP4P/F (a non-polarizable model with 3 point charges per molecule and Lennard-Jones interactions between oxygen atoms). The system examined contains 256 water molecules in periodic boundary conditions. Timings are presented for TTM3-F and TTM4-F, the underlying electrostatics model used by WHBB and MB-pol, respectively. The additional cost of WHBB and MB-pol beyond their baseline electrostatics represents the computational cost of the short-range 2B/3B polynomials. WHBB is implemented as described in Ref. 70 using the WHBB polynomial libraries provided by the authors. All timings were performed in a modified version of DL_POLY2 using a single core of a typical Intel Xeon E5-2640 based workstation.

Model Description Cost per MD step relative to q-TIP4P/F
q-TIP4P/F  Point charge  1.0 × 
TTM3-F  1 polarizable site/molecule  7.3 × 
TTM4-F  3 polarizable sites/molecule  8.5 × 
WHBB  TTM3-F electrostatics +  29 000 × 
  empirical 2B dispersion +   
  2B short-range 7th-degree polynomial +   
  3B short-range 5th-degree polynomial   
MB-pol  TTM4-F electrostatics +  47 × 
  ab initio 2B dispersion +   
  2B short-range 4th-degree polynomial +   
  3B short-range 4th-degree polynomial   
Model Description Cost per MD step relative to q-TIP4P/F
q-TIP4P/F  Point charge  1.0 × 
TTM3-F  1 polarizable site/molecule  7.3 × 
TTM4-F  3 polarizable sites/molecule  8.5 × 
WHBB  TTM3-F electrostatics +  29 000 × 
  empirical 2B dispersion +   
  2B short-range 7th-degree polynomial +   
  3B short-range 5th-degree polynomial   
MB-pol  TTM4-F electrostatics +  47 × 
  ab initio 2B dispersion +   
  2B short-range 4th-degree polynomial +   
  3B short-range 4th-degree polynomial   

The role of short- and long-range contributions to the total energy of water systems ranging from the gas-phase dimer to the liquid phase was investigated by considering four water potentials that treat the leading terms of the many-body expansion through implicit (i.e., TTM3-F and TTM4-F PESs) and explicit (i.e., WHBB and MB-pol PESs) representations. The analysis of the energetics, vibrational frequencies, and infrared intensity indicates that explicit short-range representations of 2B and 3B interactions along with a physically correct incorporation of short- and long-range contributions are necessary for an accurate representation of the water interactions, independently of the system size. These results thus suggest that atomistic water potentials built upon the many-body expansion of the interaction energy derived from “first principles” hold great promise to achieve the long-sought-after goal of describing the macroscopic properties of water across different phases from a rigorous microscopic viewpoint. A question that remains to be addressed is the extent to which these “first principles” many-body water potentials can be extended to the modeling of complex aqueous solutions.

We would like to thank Dr. Greg Schenter and Dr. Chris Mundy for several stimulating discussions on many-body effects in water, and Professor Rich Saykally for helpful discussions about the origin of the low frequency portion of the IR spectrum of liquid water. This research was supported by the National Science Foundation (Grant No. CHE-1453204) and used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (Grant No. ACI-1053575, allocation TG-CHE110009). AWG acknowledges support by the National Science Foundation (Grant No. CHE-1416571). GRM acknowledges the Department of Education for support through the GAANN fellowship program. MAM was supported through the Predictive Theory and Modeling for Materials and Chemical Science program by the Office of Basic Energy Sciences (BES), Department of Energy (DOE), and under the auspices of the US DOE by LLNL under Contract No. DE-AC52-07NA27344.

1.
D.
Hankins
,
J. W.
Moskowitz
, and
F. H.
Stillinger
,
J. Chem. Phys.
53
,
4544
(
1970
).
2.
S. S.
Xantheas
,
J. Chem. Phys.
100
,
7523
(
1994
).
3.
J. M.
Pedulla
,
F.
Vila
, and
K. D.
Jordan
,
J. Chem. Phys.
105
,
11091
(
1996
).
4.
M. P.
Hodges
,
A. J.
Stone
, and
S. S.
Xantheas
,
J. Phys. Chem. A
101
,
9163
(
1997
).
5.
J.
Cui
,
H.
Liu
, and
K. D.
Jordan
,
J. Phys. Chem. B
110
,
18872
(
2006
).
6.
U.
Góra
,
R.
Podeszwa
,
W.
Cencek
, and
K.
Szalewicz
,
J. Chem. Phys.
135
,
224102
(
2011
).
7.
R.
Khaliullin
,
E.
Cobar
,
R.
Lochan
,
A.
Bell
, and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
14
,
15328
(
2012
).
8.
K.
Morokuma
and
L.
Pedersen
,
J. Chem. Phys.
48
,
3275
(
1968
).
9.
P. A.
Kollman
and
L. C.
Allen
,
J. Chem. Phys.
51
,
3286
(
1969
).
10.
J.
Del Bene
and
J.
Pople
,
Chem. Phys. Lett.
4
,
426
(
1969
).
11.
J.
Del Bene
and
J. A.
Pople
,
J. Chem. Phys.
52
,
4858
(
1970
).
12.
D.
Hankins
,
J.
Moskowitz
, and
F.
Stillinger
,
Chem. Phys. Lett.
4
,
527
(
1970
).
13.
B. R.
Lentz
and
H. A.
Scheraga
,
J. Chem. Phys.
58
,
5296
(
1973
).
14.
H.
Popkie
,
H.
Kistenmacher
, and
E.
Clementi
,
J. Chem. Phys.
59
,
1325
(
1973
).
15.
H.
Kistenmacher
,
G. C.
Lie
,
H.
Popkie
, and
E.
Clementi
,
J. Chem. Phys.
61
,
546
(
1974
).
16.
G. C.
Lie
,
E.
Clementi
, and
M.
Yoshimine
,
J. Chem. Phys.
64
,
2314
(
1976
).
17.
O.
Matsuoka
,
E.
Clementi
, and
M.
Yoshimine
,
J. Chem. Phys.
64
,
1351
(
1976
).
18.
J.
Detrich
,
G.
Corongiu
, and
E.
Clementi
,
Chem. Phys. Lett.
112
,
426
(
1984
).
19.
U.
Niesar
,
G.
Corongiu
,
M.-J.
Huang
,
M.
Dupuis
, and
E.
Clementi
,
Int. J. Quantum Chem.
36
,
421
(
1989
).
20.
C.
Millot
and
A.
Stone
,
Mol. Phys.
77
,
439
(
1992
).
21.
C.
Millot
,
J.
Soetens
,
M.
Costa
,
M. P.
Hodges
, and
A.
Stone
,
J. Phys. Chem. A
102
,
754
(
1998
).
22.
C. J.
Burnham
and
S. S.
Xantheas
,
J. Chem. Phys.
116
,
1479
(
2002
).
23.
S. S.
Xantheas
,
C. J.
Burnham
, and
R. J.
Harrison
,
J. Chem. Phys.
116
,
1493
(
2002
).
24.
C. J.
Burnham
and
S. S.
Xantheas
,
J. Chem. Phys.
116
,
1500
(
2002
).
25.
C. J.
Burnham
and
S. S.
Xantheas
,
J. Chem. Phys.
116
,
5115
(
2002
).
26.
G. S.
Fanourgakis
and
S. S.
Xantheas
,
J. Chem. Phys.
128
,
074506
(
2008
).
27.
C. J.
Burnham
,
D. J.
Anick
,
P. K.
Mankoo
, and
G. F.
Reiter
,
J. Chem. Phys.
128
,
154519
(
2008
).
28.
S.
Habershon
,
G. S.
Fanourgakis
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
129
,
074501
(
2008
).
29.
F.
Paesani
,
S. S.
Xantheas
, and
G. A.
Voth
,
J. Phys. Chem. B
113
,
13118
(
2009
).
30.
G. R.
Medders
and
F.
Paesani
,
J. Chem. Phys.
142
,
212411
(
2015
).
31.
E.
Mas
,
R.
Bukowski
,
K.
Szalewicz
,
G.
Groenenboom
,
P.
Wormer
, and
A.
van der Avoird
,
J. Chem. Phys.
113
,
6687
(
2000
).
32.
E.
Mas
,
R.
Bukowski
, and
K.
Szalewicz
,
J. Chem. Phys.
118
,
4386
(
2003
).
33.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
van der Avoird
,
Science
315
,
1249
(
2007
).
34.
C.
Leforestier
,
K.
Szalewicz
, and
A.
van der Avoird
,
J. Chem. Phys.
137
,
014305
(
2012
).
35.
P.
Jankowski
,
G.
Murdachaew
,
R.
Bukowski
,
O.
Akin-Ojo
,
C.
Leforestier
, and
K.
Szalewicz
,
J. Phys. Chem. A
119
,
2940
(
2015
).
36.
Y.
Wang
and
J. M.
Bowman
,
J. Chem. Phys.
134
,
154510
(
2011
).
37.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
,
5395
(
2013
).
38.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
,
J. Chem. Theory Comput.
10
,
1599
(
2014
).
39.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
,
J. Chem. Theory Comput.
10
,
2906
(
2014
).
40.
H.
Partridge
and
D. W.
Schwenke
,
J. Chem. Phys.
106
,
4618
(
1997
).
41.
A. J.
Stone
,
Theory of Intermolecular Forces
(
Oxford University Press
,
1997
).
42.
H.
Liu
,
Y.
Wang
, and
J. M.
Bowman
,
J. Chem. Phys.
142
,
194502
(
2015
).
43.
A. P.
Bartók
and
G.
Csányi
,
Int. J. Quantum Chem.
115
,
1051
(
2015
).
44.
A. P.
Bartók
,
M. J.
Gillan
,
F. R.
Manby
, and
G.
Csányi
,
Phys. Rev. B
88
,
054104
(
2013
).
45.
G. S.
Tschumper
,
M. L.
Leininger
,
B. C.
Hoffman
,
E. F.
Valeev
,
H. F.
Schaefer
, and
M.
Quack
,
J. Chem. Phys.
116
,
690
(
2002
).
46.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
,
1103
(
2013
).
47.
D. M.
Bates
and
G. S.
Tschumper
,
J. Phys. Chem. A
113
,
3555
(
2009
).
48.
M. A.
Morales
,
J. R.
Gergely
,
J.
McMinis
,
J. M.
McMahon
,
J.
Kim
, and
D. M.
Ceperley
,
J. Chem. Theory Comput.
10
,
2355
(
2014
).
49.
See supplementary material at http://dx.doi.org/10.1063/1.4930194 for an analysis of the relative binding energies of the tetramer and pentamer isomers.
50.
G.
Knizia
,
T. B.
Adler
, and
H. J.
Werner
,
J. Chem. Phys.
130
,
054104
(
2009
).
51.
T. B.
Adler
,
G.
Knizia
, and
H. J.
Werner
,
J. Chem. Phys.
127
,
221106
(
2007
).
52.
K. A.
Peterson
,
T. B.
Adler
, and
H. J.
Werner
,
J. Chem. Phys.
128
,
084102
(
2008
).
53.
D. P.
Tew
,
W.
Klopper
,
C.
Neiss
, and
C.
Hattig
,
Phys. Chem. Chem. Phys.
9
,
1921
(
2007
).
54.
F. A.
Bischoff
,
S.
Wolfsegger
,
D. P.
Tew
, and
W.
Klopper
,
Mol. Phys.
107
,
963
(
2009
).
55.
S. S.
Xantheas
,
Chem. Phys.
258
,
225
(
2000
).
56.
B. H.
Wells
and
S.
Wilson
,
Chem. Phys. Lett.
101
,
429
(
1983
).
57.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
,
M.
Schütz
 et al, molpro, version 2012.1, a package ofab initio programs, 2012, see http://www.molpro.net.
58.
B.
Santra
,
A.
Michaelides
,
M.
Fuchs
,
A.
Tkatchenko
,
C.
Filippi
, and
M.
Scheffler
,
J. Chem. Phys.
129
,
194111
(
2008
).
59.
M. J.
Gillan
,
F. R.
Manby
,
M. D.
Towler
, and
D.
Alfè
,
J. Chem. Phys.
136
,
244105
(
2012
).
60.
D.
Alfè
,
A. P.
Bartók
,
G.
Csányi
, and
M. J.
Gillan
,
J. Chem. Phys.
141
,
014104
(
2014
).
61.
J. C.
Howard
and
G. S.
Tschumper
,
J. Chem. Theory Comput.
11
,
2126
(
2015
).
62.
E.
Miliordos
,
E.
Aprà
, and
S. S.
Xantheas
,
J. Chem. Phys.
139
,
114302
(
2013
).
63.
J. E.
Bertie
and
Z.
Lan
,
Appl. Spectrosc.
50
,
1047
(
1996
).
64.
J. J.
Perez
,
J. H. R.
Clarke
, and
A.
Hinchliffe
,
Chem. Phys. Lett.
104
,
583
(
1984
).
65.
B.
Skwara
,
W.
Bartkowiak
,
A.
Zawada
,
R. W.
Góra
, and
J.
Leszczynski
,
Chem. Phys. Lett.
436
,
116
(
2007
).
66.
D. W.
Schwenke
,
J. Phys. Chem.
100
,
2867
(
1996
).
67.
C. J.
Tainter
,
P. A.
Pieniazek
,
Y.
Lin
, and
J. L.
Skinner
,
J. Chem. Phys.
134
,
184501
(
2011
).
68.
G. R.
Medders
and
F.
Paesani
,
J. Chem. Theory Comput.
11
,
1145
(
2015
).
69.
L.
Lodi
,
J.
Tennyson
, and
O. L.
Polyansky
,
J. Chem. Phys.
135
,
034113
(
2011
).
70.
Y.
Wang
,
X.
Huang
,
B. C.
Shepler
,
B. J.
Braams
, and
J. M.
Bowman
,
J. Chem. Phys.
134
,
094509
(
2011
).
71.
B.
Guillot
,
J. Chem. Phys.
95
,
1543
(
1991
).
72.
S.
Habershon
,
T. E.
Markland
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
131
,
024501
(
2009
).
73.

“mbpol_openmm_plugin,” reference implementation of the MB-pol water potential. Multicore and GPU implementations for actual molecular dynamics production runs will become available in the near future, https://github.com/paesanilab/mbpol_openmm_plugin.

74.
P.
Eastman
and
V.
Pande
,
Comput. Sci. Eng.
12
,
34
(
2010
).

Supplementary Material