On the representation of many-body interactions in water
Gregory R. Medders, Andreas W. Gotz, Miguel A. Morales, Francesco Paesani
aa r X i v : . [ c ond - m a t . m t r l - s c i ] A ug On the representation of many-body interactions in water
Gregory R. Medders, Andreas W. G¨otz, Miguel A. Morales, Pushp Bajaj, and Francesco Paesani Department of Chemistry and Biochemistry, University of California, San Diego, La Jolla, California 92093,United States San Diego Supercomputer Center, University of California, San Diego, La Jolla, California 92093,United States Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, California 94550,U.S.A. (Dated: 3 October 2018)
Recent work has shown that the many-body expansion of the interaction energy can be used to developanalytical representations of global potential energy surfaces (PESs) for water. In this study, the role ofshort- and long-range interactions at different orders is investigated by analyzing water potentials that treatthe leading terms of the many-body expansion through implicit (i.e., TTM3-F and TTM4-F PESs) andexplicit (i.e., WHBB and MB-pol PESs) representations. It is found that explicit short-range representationsof 2-body and 3-body interactions along with a physically correct incorporation of short- and long-rangecontributions are necessary for an accurate representation of the water interactions from the gas to thecondensed phase. Similarly, a complete many-body representation of the dipole moment surface is found tobe crucial to reproducing the correct intensities of the infrared spectrum of liquid water.
I. INTRODUCTION
It is known that the global potential energy sur-face (PES) of a system containing N interacting wa-ter molecules can be formally expressed in terms of themany-body expansion of the interaction energy as a sumover n -body terms with 1 ≤ n ≤ N , V N ( x , . . . , x N ) = X a V (1B) ( x a ) + X a>b V (2B) ( x a , x b )+ X a>b>c V (3B) ( x a , x b , x c ) + . . . + V ( N B) ( x , . . . , x N ) . (1)Here, x i collectively denotes the coordinates of all atomsin the i -th water molecule, V (1B) is the one-body poten-tial describing the energy required to deform an individ-ual water molecule from its equilibrium geometry, and allhigher n -body interactions, V ( n B) , are defined recursivelythrough V ( n B) ( x , . . . , x n ) = V n ( x , . . . , x n ) − X a V (1B) ( x a ) − X a>b V (2B) ( x a , x b ) − . . . − X a > ··· >a n − V (( n − ( x a , x a , . . . , x a n − ) . (2)The rapid convergence of Equation (1), demonstrated inseveral studies of water clusters, suggests that the PESassociated with a system containing N water moleculescan in principle be represented as a sum of low-order in-teractions that are amenable to accurate calculations us-ing correlated electronic structure methods [e.g., coupledcluster with single, double, and perturbative triple exci-tations, CCSD(T), which currently represents the gold standard in quantum chemistry].The theoretical modeling of many-body effects in wa-ter began in the late 1960s when self-consistent field(SCF) calculations were carried out on small waterclusters. These studies concluded that nonadditiveeffects in water are generally nonnegligible, with 3B con-tributions being as large as 10 – 15% of the total pairinteraction for ring structures. 4B effects were estimatedto be, on average, ∼
1% of the total pair interaction.The first attempt to derive potential energy functionsfor water from ab initio data was made by Clementi andco-workers who fitted dimer energies calculated at theHartree-Fock (HF) level of theory to an analytical rep-resentation of the 2B interactions.
Later, configu-ration interaction (CI) calculations were used to derivethe pairwise additive MCY potential, to which 3B and4B terms were subsequently added through a classicaldescription of polarization effects. The first water po-tential including explicit 2B and 3B terms, derived re-spectively from 4 th order M¨oller-Plesset (MP4) and HFcalculations, along with a classical description of higher-body polarization interactions, was reported in Ref. 19.Following these pioneering efforts, the ASP potentialswith rigid monomers were derived using intermolecularperturbation theory. In the 2000s, Xantheas and co-workers introducedthe TTM (Thole-type model) water potentials, whichfor the first time made use of a highly accurate 1Bterm derived from ab initio calculations.
The lat-est versions (TTM3-F and TTM4-F ) were shown toreproduce the properties of water clusters, liquid wa-ter, and ice reasonably well, although some inaccura-cies were identified in the calculations of vibrationalspectra. Around the same time, 2B and 3B poten-tials with rigid water monomers were derived by Sza-lewicz and co-workers from symmetry-adapted pertur-bation theory (SAPT).
These studies eventually ledto the development of the rigid-monomer CC-pol poten-tial, a 25-site model with explicit 2B and 3B terms fittedto CCSD(T)-corrected MP2 dimer energies and SAPTtrimer energies, respectively, with higher-body terms be-ing represented through classical polarization. CC-polwas shown to accurately reproduce the vibration-rotationtunneling spectrum of the water dimer and to predict thestructure of liquid water in reasonable agreement withthe experimental data. Within the CC-pol scheme, a re-fined 2B term with explicit dependence on the monomerflexibility has also been developed.
More recently, the full-dimensional WHBB and MB-pol many-body potentials with flexible monomerswere developed. WHBB includes explicit 2B and 3Bterms fitted to CCSD(T) and MP2 data, respectively,with long-range many-body effects being represented bythe same Thole-type model used in the TTM3-F poten-tial. MB-pol was derived from fits to CCSD(T) energiescalculated for both water dimers and trimers in the com-plete basis set (CBS) limit and includes many-body ef-fects within a modified version of the polarization modelemployed by the TTM4-F potential. Although WHBBhas successfully been applied to static calculations of sev-eral water properties, to date MB-pol is the only many-body potential that has been consistently employed inquantum molecular dynamics (MD) simulations that ac-curately predicted the properties of water from the gasto the condensed phase.While the calculations with many-body potentials re-ported in the literature demonstrate that Equation (1)can be effectively used to develop accurate molecular-level representations of the water interactions, the roleof short- and long-range contributions at different ordershas not been fully investigated. Here, we seek to ad-dress this issue by analyzing water potentials that treatthe leading terms of the many-body expansion throughimplicit (e.g., TTM3-F and TTM4-F) and explicit (e.g.,WHBB and MB-pol) representations, with a specific fo-cus on how these terms are defined and incorporated asa function of the intermolecular separations. The articleis organized as follows: The four potentials are brieflydescribed in section 2, an analysis of the energetics aswell as of the vibrational frequencies and intensities forwater systems ranging from the gas-phase dimer to theliquid phase is reported in section 3, and the conclusionsare given in section 4. II. METHODS
The four potentials studied here were chosen becausethey treat the intramolecular distortions on equal foot-ing, i.e., through the 1B PES developed by Partridgeand Schwenke. In this way, the ability of the modelsto predict water properties is a direct reflection of thetreatment of the intermolecular interactions. In addition,WHBB and MB-pol effectively share the same inductionscheme 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 ex-pansion.The dominant terms of the many-body water inter-actions are the 2B and 3B terms, which can be con-veniently split into short- and long-range contributions, V (2B , = V (2B , + V (2B , . At the two-body level,long-range interactions are dominated by electrostaticsand dispersion, while exchange-repulsion and charge-transfer become increasingly important at short-range.3B interactions in water, on the other hand, arise primar-ily from induced interactions at long-range and exchange-repulsion when the monomer electron densities overlap. Both TTM3-F and TTM4-F include permanent and in-duced electrostatics as well as dispersion and repulsion atthe 2B level, while all higher order terms are representedthrough many-body induction.
Although WHBB and MB-pol seek to exploit the rangeseparation of the low-order water interactions, the twopotentials are intrinsically different both in philosophyand by construction. In WHBB, the short-range 2Bterm is described by a 7 th -degree permutationally invari-ant polynomial that smoothly transitions for an oxygen-oxygen separation between 6.5 and 7.5 ˚A into the long-range component described by the 2B TTM3-F poten-tial. The 3B term only includes a short-range component,which is represented by either a 5 th - (WHBB5) or a 6 th -degree (WHBB6) permutationally invariant polynomialthat dies off as the largest oxygen-oxygen separation be-tween two molecules of the trimer approaches 8.0 ˚A. Allhigher-order interactions ( n ≥
4) are described by thepolarization model employed in the TTM3-F potential.By construction, WHBB thus employs a strict separationat the 2B and 3B levels between short- and long-range in-teractions that are described in a completely independentway and are essentially disentangled from all higher-ordercontributions. It should be noted that a simplified ver-sion of WHBB including only 1B, 2B, and 3B interactionswith shorter cutoffs was also used, which, however, can-not be implemented in a straightforward way in standardsoftware for molecular simulations in periodic boundaryconditions. For this reason, only calculations with thefull WHBB5 implementation of Ref. 36 are reported inthe following analysis.On the other hand, MB-pol can be viewed as a classicalpolarizable model supplemented by short-range 2B and3B terms that effectively represent quantum-mechanicalinteractions arising from the overlap of the monomer elec-tron densities. Specifically, at all separations, V (2B) in-cludes (damped) dispersion forces derived from ab initio computed asymptotic expansions of the dispersion energyalong with electrostatic contributions due to the inter-actions between the molecular permanent and inducedmoments. At short-range, V (2B) is supplemented bya 4 th -degree permutationally invariant polynomial thatsmoothly switches to zero as the oxygen-oxygen separa-tion in the dimer approaches 6.5 ˚A. Similarly, the MB-pol 3B term, V (3B) , includes a 3B polarization term atall separations, which is supplemented by a short-range4 th -degree permutationally invariant polynomial that ef-fectively corrects for the deficiencies of a purely classicalrepresentation of the 3B interactions in regions wherethe electron densities of the three monomers overlap.This short-range 3B potential is smoothly switched offonce the oxygen-oxygen separation between any watermolecule and the other two water molecules of a trimerreaches a value of 4.5 ˚A. In MB-pol, all induced interac-tions are described through many-body polarization us-ing a slightly modified version of TTM4-F. MB-pol thuscontains many-body effects at all monomer separationsas well as at all orders, in an explicit way up to the thirdorder and in a mean-field fashion at all higher orders. Itshould be noted that alternative functions to the permu-tationally invariant polynomials (e.g., GAP ) have beensuggested and successfully employed in modeling short-range many-body interactions. III. RESULTS AND DISCUSSION
The effects of the different representations of themany-body interactions employed by TTM3-F, TTM4-F,WHBB, and MB-pol are investigated through the anal-ysis of the energetics of water systems ranging from thegas-phase dimer to the liquid phase. Figure 1a showsa comparison of the dimer binding energies calculatedfor the global minimum configuration of each PES. Al-though all potentials predict binding energies within ∼ the analysis of themany-body water interactions reported in Ref. 46 showsthat the TTM3-F and TTM4-F interaction energies devi-ate significantly from the CCSD(T) reference data whendimer configurations different from the correspondingequilibrium geometries are considered. Specifically, rootmean square deviations (RMSDs) of 0.73 kcal/mol and0.44 kcal/mol were obtained in Ref. 46 for TTM3-F andTTM4-F, respectively, from comparisons with ∼ 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 num-ber of molecules in the system of interest increases. Thisis shown in Figure 1b, where the relative energies of low-lying hexamer isomers calculated with the four potentialsare compared with the corresponding ab initio values. The hexamer isomers hold a special place in the exper-imental and theoretical studies of water because theyare the smallest clusters in which the lowest energy con- figurations correspond to fully three-dimensional struc-tures. For this reason, the hexamer serves as a proto-typical model for the hydrogen-bond networks observedin condensed phases. While TTM3-F, WHBB, and MB-pol (qualitatively) reproduce the isomer energy ordering,large deviations are found between the TTM4-F predic-tions and the reference ab initio values. Similar resultsare also obtained for other small water clusters, includ-ing the tetramer and pentamer isomers, shown in Fig-ures S1 and S2 of the supplemental material. . Sinceit was shown that the Thole-type polarization schemeemployed by the TTM4-F potential correctly captureshigher-order water interactions, the low accuracyshown by TTM4-F in describing the relative energies ofthe hexamer isomers is possibly due to intrinsic deficien-cies in the 2B term. In this context, the performanceof the TTM3-F potential, which provides the best agree-ment with the reference ab initio values, is somewhat sur-prising given the large deviations already seen at both the2B and 3B levels, and may suggest some fortuitous can-cellation of errors. Based on this, it is interesting to notethat the isomer energies predicted by WHBB, which em-ploys the same electrostatic model as TTM3-F, deviatefrom the reference data by as much as 2.0 kcal/mol. Theorigin of these deviations may result from inaccuraciesin the short-range 2B and 3B WHBB polynomials, in-compatibility between the effective many-body represen-tation encoded in the TTM3-F electrostatic model withthe short-range WHBB polynomials, or both. On theother hand, MB-pol, which employs a slightly modifiedversion of the TTM4-F electrostatic model, predicts iso-mer energies in good agreement with the reference data,suggesting that the MB-pol short-range 2B and 3B termsaccurately correct for the deficiencies that are intrinsic tothe purely classical representation of these interactionsencoded in the TTM4-F electrostatic model.To provide more quantitative insights into the differ-ent performance of the four potentials in describing therelative energies of the hexamer isomers, a many-bodydecomposition of the corresponding interaction energies(Table I) was carried out at the CCSD(T)-F12b/VTZ-F12 level of theory with the VTZ-F12 basis set, which was shown to effectively provide the same accu-racy as CCSD(T)/CBS calculations at a reduced compu-tational cost. As has been well established, the in-teraction energies are dominated by the 2B and 3B terms,with 4B and higher terms making largest (but still rel-atively small) contributions for cyclic structures. It isinteresting to note that, while the hexamer binding en-ergies predicted by TTM3-F are in excellent agreementwith the reference calculations (Figure 1b), the resultsof Table I demonstrate that this agreement is clearlyachieved by a fortuitous cancellation of errors. For ex-ample, in the prism isomer, the TTM3-F 2B interactionis too attractive by ∼ ∼ 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). Alsoshown as open symbols are the binding energies calculated using the reference CCSD(T) dimer geometry. b) Comparison ofthe relative binding energies (∆E) of the water hexamer isomers calculated with the TTM3-F, TTM4-F, WHBB, and MB-polpotentials using the ab initio geometries of Ref. 47. Also shown as a reference are the ab initio values from Ref. 47. c) Meanabsolute 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 frompath-integral molecular dynamics simulations of water performed with the vdW-DF and vdW-DF2 functionals. The statisticalerrors in the mean absolute differences arising from the QMC calculations are less than 0.006 mHa/molecule. See Ref. 48 forspecific details on the QMC and DFT calculations. which accurately represents the electrostatic interactions(as shown by 3B contributions that lie within 1 kcal/molof the reference values for all low-lying hexamer isomers),is unable to correctly describe the 2B interactions, whichhave significant contributions from exchange/repulsionand charge transfer. Thus, without the benefit of for-tuitous cancellation of errors (that can, to some extent,be encoded into models through empirical parametriza-tion), potentials with accurate electrostatics but withouta comparably accurate treatment of short-range quan-tum mechanical effects may be expected to have rela-tively poor performance. Having been explicitly derivedfrom many-body expansions of the interaction energiesobtained from correlated electronic structure data, bothWHBB and MB-pol closely reproduce the CCSD(T)-F12reference many-body terms. However, nonnegligible dis-crepancies in the 3B and 4B terms result in WHBB beingappreciably less accurate than MB-pol at reproducing thehexamer interaction energies.The differences between the four potentials becomemore evident when their accuracy is assessed againstbenchmark quantum Monte Carlo (QMC) interaction en-ergies calculated for liquid water configurations extractedfrom path-integral molecular dynamics simulations of thevdW-DF and vdW-DF2 functionals. QMC has beenshown to be a reliable benchmark in the study of smallwater clusters, producing relative energies with an accu-racy comparable to that of CCSD(T).
As a measureof accuracy, the mean absolute deviation (MAD) betweenthe energies E i obtained with each of the four potentialsand the reference QMC energies E QMC was calculated asMAD = 1 N c N c X i =1 (cid:12)(cid:12)(cid:12)(cid:16) E i − E QMC i (cid:17) − (cid:10) E − E QMC (cid:11)(cid:12)(cid:12)(cid:12) . (3) Here, N c is the number of water configurations and (cid:10) E − E QMC (cid:11) = (1 /N c ) P N c i =1 ( E i − E QMC i ) is the averageenergy difference for all configurations, effectively align-ing the zero of energy with the reference QMC data. Fig-ure 1c shows that the mean absolute deviation associatedwith WHBB is ∼ ∼
15 times larger than the MB-pol value. Interestingly,MB-pol also achieves better accuracy than all DFT mod-els analyzed in Ref. 48. The TTM3-F and TTM4-F re-sults indicate that effective representations of the many-body interactions through classical polarization modelscan provide a reasonable description of the energetics as-sociated with the liquid phase (i.e., comparable with thatprovided by most of the DFT models currently used inwater simulations), albeit through empirical parameteri-zation. It should be noted that, as shown by the resultsreported in Table I and the many-body analysis of dif-ferent DFT models carried out in Ref. 60, cancellationof errors between different terms of the many-body ex-pansion of the interaction energies may also affect theenergetics of bulk systems described by both DFT andclassical polarizable models. On the other hand, the in-creased accuracy of MB-pol relative to TTM4-F indicatesthat short-range many-body interactions beyond a purelyclassical electrostatic representation are necessary for acorrect, molecular-level description of the liquid phase.The significantly different accuracies associated with theWHBB/TTM3-F and MB-pol/TTM4-F potentials alsosuggest that the correct integration of explicit short-range and effective long-range many-body interactions iscritical for ensuring the transferability of the potentialacross different phases.The ability of TTM3-F, TTM4-F, WHBB, and MB-pol to represent the global multidimensional PESs of wa-
TABLE I. Many-body decomposition of the interaction energy of low-lying hexamer isomers. All energies are in kcal/mol. TheCCSD(T)-F12b calculations were performed with the VTZ-F12 basis set and corrected for basis set superposition errorthrough the site-site function counterpoise method using the Molpro quantum chemistry package. (a) Prism CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol2B − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . .
06 0 .
03 0 .
06 0 .
03 0 . .
00 0 .
00 0 .
00 0 .
00 0 . − . − . − . − . − . (b) Cage CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol2B − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . .
01 0 .
01 0 .
03 0 .
01 0 . .
00 0 .
00 0 .
00 0 .
00 0 . − . − . − . − . − . (c) Book-1 CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol2B − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . .
00 0 .
00 0 .
00 0 .
00 0 . − . − . − . − . − . (d) Book-2 CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol2B − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − .
01 0 . − .
01 0 . .
00 0 .
00 0 .
00 0 .
00 0 . − . − . − . − . − . (e) Bag CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol2B − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − .
02 0 . − . − . .
01 0 .
00 0 .
00 0 .
00 0 . − . − . − . − . − . (f) Cyclic-chair CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol2B − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − .
01 0 . − .
01 0 . − . − . − . − . − . (g) Cyclic-boat-1 CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol2B − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − .
01 0 . − .
01 0 . − . − . − . − . − . (h) Cyclic-boat-2 CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol2B − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − . − .
01 0 . − .
01 0 . − . − . − . − . − . ter systems with an increasing number of molecules isdirectly reflected in the accuracy with which the fourpotentials predict the associated vibrational frequencies.Harmonic frequencies for the water hexamer have beenreported at the CCSD(T)/aug-cc-pVDZ level. Morerecently, ab initio reference harmonic frequencies forsmall water clusters have been obtained through two-body:many-body CCSD(T):MP2 calculations using thecc-pVQZ and aug-cc-pVQZ basis sets for H and O atoms,respectively, which enables a more complete treatmentof the basis set through the hybrid two-body:many-bodyapproach. As shown in Figure 2, the WHBB and MB-pol PESs result in harmonic frequencies for the waterdimer that deviate by less than 10 cm − from the ab ini-tio frequencies reported in Ref. 61. These results are con-sistent with the analysis of the dimer vibration-rotationtunneling spectra of WHBB and MB-pol, both of whichexhibit nearly perfect agreement with the correspondingexperimental data. In contrast, relatively large devi-ations are seen for vibrational frequencies of both theTTM3-F and TTM4-F potentials. The different perfor-mance of the four potentials in predicting the relative energies of the hexamer isomers clearly leads to differentlevels of agreement with the ab initio harmonic frequen-cies. Independently of the isomer, the MB-pol valuesconsistently lie within 50 cm − of the reference valueswhile the WHBB harmonic frequencies can deviate, insome cases, by more than 100 cm − . Substantially largerdeviations, up to ∼
200 cm − , are instead obtained withboth TTM3-F and TTM4-F, reinforcing the notion thatpurely classical representations of the many-body waterinteractions are likely not sufficient to fully reproduce thecomplexity of the underlying multidimensional PESs. In-terestingly, all potentials predict somewhat larger devi-ations for the cyclic-chair isomer of the water hexamer,supporting early observations that many-body effects inwater are relatively more important for ring structures. While the differences in the PESs clearly affect theunderlying vibrational structure, as shown in Figure 2,the infrared activity of those vibrations is ultimately dic-tated by the associated dipole moment surfaces (DMSs).Many-body representations of higher-order electric prop-erties for molecular systems were characterized begin-ning in the early 1980s. The first analysis of many-
FIG. 2. Deviations from the 2-body:many-body CCSD(T):MP2 harmonic vibrational frequencies of (H O) n isomers with n= 2 and n = 6 calculated using TTM3-F, TTM4-F, WHBB5, and MB-pol.
500 1000 1500 2000 2500 3000 3500 4000Frequency (cm -1 )0e+001e+042e+04 I n f r a r ed I n t en s i t y FIG. 3. Decomposition of the IR spectrum obtained fromCMD trajectories with the MB-pol PES in terms of the many-body components of the MB- µ DMS. 1B-Dip indicates thatthe one-body (gas-phase monomer) dipoles were used to cal-culate the dipole of the molecules sampled along the MB-polCMD trajectories, from which the IR spectrum was calcu-lated. (1B+2B)-Dip indicates that short-ranged two-bodydipoles were used in addition to the one-body dipoles. MB-Dip is the full MB- µ many-body dipole. The spectra weresmoothed to facilitate the comparison between the line shapesobtained using the different approximations. The experimen-tal IR spectrum was taken from Ref. 63. body effects on the dipole moment of polar moleculeswas reported by Skwara et al. As noted in 1996 bySchwenke, molecular dipole moments can efficientlybe represented in terms of geometry-dependent effectivecharges multiplied by their Cartesian positions. In line with these ideas, the first representation of the 2B dipoleof water was reported as part of the WHBB suite. Arefined version of this 1B + 2B DMS has been usedto calculate the transition dipole moments which werethen used to modulate the WHBB5 frequency distribu-tions calculated within the local monomer approximationfrom configurations of liquid water extracted from MDsimulations with the rigid E3B. Neglecting both 3Band higher-order contributions to the dipole momentsand dynamical effects (e.g., motional narrowing), goodagreement in the relative intensities was obtained withthe measured IR spectrum of liquid water. Recently, acomplete many-body representation of the DMS of wa-ter (MB- µ ) was reported, consisting of the one-bodyDMS of Lodi et al. , an explicit two-body term, anda slightly modified version of TTM4-F polarization forlong-range two-body and all higher-order many-body in-duced dipole moments. The particular functional formof MB- µ was derived from a systematic analysis of themany-body convergence of the electrostatic properties ofwater. To investigate the effects of many-body dipole mo-ments on the IR spectrum of liquid water, many-bodymolecular dynamics (MB-MD) simulations, within thecentroid molecular dynamics (CMD) formalism, werecarried out with the MB-pol PES in combination withthe MB- µ DMS. As shown in Figure 3, while the overallshape of the IR spectrum is captured in the 1B + 2Brepresentation of the DMS, 3B and higher-order contri-butions to the dipole moment significantly affect the IR
TABLE II. Cost per molecular dynamics step for different PESs relative to q-TIP4P/F (a non-polarizable model with 3point charges per molecule and Lennard-Jones interactions between oxygen atoms). The system examined contains 256 watermolecules in periodic boundary conditions. Timings are presented for TTM3-F and TTM4-F, the underlying electrostaticsmodel used by WHBB and MB-pol, respectively. The additional cost of WHBB and MB-pol beyond their baseline electrostaticsrepresents the computational cost of the short-range 2B/3B polynomials. WHBB is implemented as described in Ref. 70 usingthe WHBB polynomial libraries provided by the authors. All timings were performed in a modified version of DL POLY2 usinga single core of a typical Intel Xeon E5-2640 based workstation.
Model Description Cost per MD step relative to q-TIP4P/Fq-TIP4P/F Point charge 1.0xTTM3-F 1 polarizable site/molecule 7.3xTTM4-F 3 polarizable sites/molecule 8.5xWHBB TTM3-F electrostatics + 29,000xempirical 2B dispersion +2B short-range 7 th -degree polynomial +3B short-range 5 th -degree polynomialMB-pol TTM4-F electrostatics + 47x ab initio
2B dispersion +2B short-range 4 th -degree polynomial +3B short-range 4 th -degree polynomialintensities. 2B contributions to the DMS lower the bendintensity while they contribute to ∼
70% of the intensityof the OH stretch band. Importantly, while 2B contribu-tions are critical to appearance of the shoulder at ∼ − corresponding to the hydrogen-bonding stretch, the1B contributions are non-negligible. This indicates thatsome rotational motion is also involved in the hydrogen-bonding stretch. Nevertheless, the absolute intensities ofthe librational, bending, and stretching bands are onlyrecovered when all many-body effects are included in thecalculation of the dipole moment. These results supportearly observations based on molecular dynamics simula-tions with the SPC potential supplemented with a 3Bdipole-induced dipole term showing that the calculatedfar infrared spectrum of liquid water was in better agree-ment with experiment relative to results obtained includ-ing only a 2B description of the dipole moment. While using a highly accurate PES is a prerequisitefor a physically correct description of the water proper-ties at the molecular level, an appropriate balance be-tween accuracy and efficiency is often critical when de-ciding which potential to employ in computer simulationssince the computational cost directly affects the abilityto calculate statistically converged quantities. As shownin Table II, a performance analysis carried out on a sin-gle Intel Xeon E5-2640 processor for a system consistingof 256 water molecules in periodic boundary conditionsindicates that MB-pol is able to achieve a high level ofaccuracy at a cost of 47x that of the fixed point chargeq-TIP4P/F model and ∼ ∼ for the OpenMM toolkit for molecular simulations. IV. CONCLUSION
The role of short- and long-range contributions to thetotal energy of water systems ranging from the gas-phasedimer to the liquid phase was investigated by consider-ing four water potentials that treat the leading terms ofthe many-body expansion through implicit (i.e., TTM3-Fand TTM4-F PESs) and explicit (i.e., WHBB and MB-pol PESs) representations. The analysis of the energet-ics, vibrational frequencies, and infrared intensity indi-cates that explicit short-range representations of 2B and3B interactions along with a physically correct incorpo-ration of short- and long-range contributions are neces-sary for an accurate representation of the water inter-actions, independently of the system size. These resultsthus suggest that atomistic water potentials built uponthe many-body expansion of the interaction energy de-rived from ”first principles” hold great promise to achievethe long-sought-after goal of describing the macroscopicproperties of water across different phases from a rig-orous microscopic viewpoint. A question that remainsto be addressed is the extent to which these ”first prin-ciples” many-body water potentials can be extended tothe modeling of complex aqueous solutions.
V. ACKNOWLEDGEMENTS
This research was supported by the National ScienceFoundation (grant number CHE-1453204) and used theExtreme Science and Engineering Discovery Environ-ment (XSEDE), which is supported by the National Sci-ence Foundation (grant number ACI-1053575, alloca-tion TG-CHE110009). AWG acknowledges support bythe National Science Foundation (grant CHE-1416571).GRM acknowledges the Department of Education forsupport through the GAANN fellowship program. MAMwas supported through the Predictive Theory and Mod-eling for Materials and Chemical Science program by theOffice of Basic Energy Sciences (BES), Department ofEnergy (DOE) and under the auspices of the US DOE byLLNL under Contract DE-AC52-07NA27344. We thankProf. Rich Saykally for helpful discussions about the ori-gin of the low frequency portion of the IR spectrum ofliquid water. D. Hankins, J. W. Moskowitz, and F. H. Stillinger, J. Chem.Phys. , 4544 (1970). S. S. Xantheas, J. Chem. Phys. , 7523 (1994). J. M. Pedulla, F. Vila, and K. D. Jordan, J. Chem. Phys. ,11091 (1996). M. P. Hodges, A. J. Stone, and S. S. Xantheas, J. Phys. Chem.A , 9163 (1997). J. Cui, H. Liu, and K. D. Jordan, J. Phys. Chem. B , 18872(2006). U. G´ora, R. Podeszwa, W. Cencek, and K. Szalewicz, J. Chem.Phys. , 224102 (2011). R. Khaliullin, E. Cobar, R. Lochan, A. Bell, and M. Head-Gordon, Phys. Chem. Chem. Phys. , 15328 (2012). K. Morokuma and L. Pedersen, J. Chem. Phys. , 3275 (1968). P. A. Kollman and L. C. Allen, J. Chem. Phys. , 3286 (1969). J. Del Bene and J. Pople, Chem. Phys. Lett. , 426 (1969). J. Del Bene and J. A. Pople, J. Chem. Phys. , 4858 (1970). D. Hankins, J. Moskowitz, and F. Stillinger, Chem. Phys. Lett. , 527 (1970). B. R. Lentz and H. A. Scheraga, J. Chem. Phys. , 5296 (1973). H. Popkie, H. Kistenmacher, and E. Clementi, J. Chem. Phys. , 1325 (1973). H. Kistenmacher, G. C. Lie, H. Popkie, and E. Clementi, J.Chem. Phys. , 546 (1974). G. C. Lie, E. Clementi, and M. Yoshimine, J. Chem. Phys. ,2314 (1976). O. Matsuoka, E. Clementi, and M. Yoshimine, J. Chem. Phys. , 1351 (1976). J. Detrich, G. Corongiu, and E. Clementi, Chem. Phys. Lett. , 426 (1984). U. Niesar, G. Corongiu, M.-J. Huang, M. Dupuis, andE. Clementi, Int. J. Quantum Chem. , 421 (1989). C. Millot and A. Stone, Mol. Phys. , 439 (1992). C. Millot, J. Soetens, M. Costa, M. P. Hodges, and A. Stone, J.Phys. Chem. A , 754 (1998). C. J. Burnham and S. S. Xantheas, J. Chem. Phys. , 1479(2002). S. S. Xantheas, C. J. Burnham, and R. J. Harrison, J. Chem.Phys. , 1493 (2002). C. J. Burnham and S. S. Xantheas, J. Chem. Phys. , 1500(2002). C. J. Burnham and S. S. Xantheas, J. Chem. Phys. , 5115(2002). G. S. Fanourgakis and S. S. Xantheas, J. Chem. Phys. ,074506 (2008). C. J. Burnham, D. J. Anick, P. K. Mankoo, and G. F. Reiter,J. Chem. Phys. , 154519 (2008). S. Habershon, G. S. Fanourgakis, and D. E. Manolopoulos, J.Chem. Phys. , 074501 (2008). F. Paesani, S. S. Xantheas, and G. A. Voth, J. Phys. Chem. B , 13118 (2009). G. R. Medders and F. Paesani, J. Chem. Phys. , 212411(2015). E. Mas, R. Bukowski, K. Szalewicz, G. Groenenboom,P. Wormer, and A. van der Avoird, J. Chem. Phys. , 6687(2000). E. Mas, R. Bukowski, and K. Szalewicz, J. Chem. Phys. ,4386 (2003). R. Bukowski, K. Szalewicz, G. C. Groenenboom, and A. van derAvoird, Science , 1249 (2007). C. Leforestier, K. Szalewicz, and A. van der Avoird, J. Chem.Phys. , 014305 (2012). P. Jankowski, G. Murdachaew, R. Bukowski, O. Akin-Ojo,C. Leforestier, and K. Szalewicz, J. Phys. Chem. A , 2940(2015). Y. Wang and J. M. Bowman, J. Chem. Phys. , 154510 (2011). V. Babin, C. Leforestier, and F. Paesani, J. Chem. Theory Com-put. , 5395 (2013). V. Babin, G. R. Medders, and F. Paesani, J. Chem. TheoryComput. , 1599 (2014). G. R. Medders, V. Babin, and F. Paesani, J. Chem. TheoryComput. , 2906 (2014). H. Partridge and D. W. Schwenke, J. Chem. Phys. , 4618(1997). A. J. Stone,
Theory of Intermolecular Forces (Oxford UniversityPress, 1997). H. Liu, Y. Wang, and J. M. Bowman, J. Chem. Phys. ,194502 (2015). A. P. Bart´ok and G. Cs´anyi, Int. J. Quantum Chem. , in press.DOI:10.1002/qua.24927 (2015). A. P. Bart´ok, M. J. Gillan, F. R. Manby, and G. Cs´anyi, Phys.Rev. B , 054104 (2013). G. S. Tschumper, M. L. Leininger, B. C. Hoffman, E. F. Valeev,H. F. Schaefer, and M. Quack, J. Chem. Phys. , 690 (2002). G. R. Medders, V. Babin, and F. Paesani, J. Chem. TheoryComput. , 1103 (2013). D. M. Bates and G. S. Tschumper, J. Phys. Chem. A , 3555(2009). M. A. Morales, J. R. Gergely, J. McMinis, J. M. McMahon,J. Kim, and D. M. Ceperley, J. Chem. Theory Comput. ,2355 (2014). “See supplemental material at [url will be inserted by aip] foran analysis of the relative binding energies of the tetramer andpentamer isomers.”. G. Knizia, T. B. Adler, and H. J. Werner, J. Chem. Phys. ,054104 (2009). T. B. Adler, G. Knizia, and H. J. Werner, J. Chem. Phys. ,221106 (2007). K. a. Peterson, T. B. Adler, and H. J. Werner, J. Chem. Phys. , 084102 (2008). D. P. Tew, W. Klopper, C. Neiss, and C. Hattig, Phys. Chem.Chem. Phys. , 1921 (2007). F. A. Bischoff, S. Wolfsegger, D. P. Tew, and W. Klopper, Mol.Phys. , 963 (2009). S. S. Xantheas, Chem. Phys. , 225 (2000). B. H. Wells and S. Wilson, Chem. Phys. Lett. , 429 (1983). H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Sch¨utz,and Others, “MOLPRO, version 2012.1, a package of ab initioprograms,” (2012). B. Santra, A. Michaelides, M. Fuchs, A. Tkatchenko, C. Filippi,and M. Scheffler, J. Chem. Phys. , 194111 (2008). M. J. Gillan, F. R. Manby, M. D. Towler, and D. Alf`e, J. Chem.Phys. , 244105 (2012). D. Alf`e, A. P. Bart´ok, G. Cs´anyi, and M. J. Gillan, J. Chem.Phys. , 014104 (2014). J. C. Howard and G. S. Tschumper, J. Chem. Theory Comput. , 2126 (2015). E. Miliordos, E. Apr`a, and S. S. Xantheas, J. Chem. Phys. ,114302 (2013). J. E. Bertie and Z. Lan, Appl. Spectrosc. , 1047 (1996). J. J. Perez, J. H. R. Clarke, and A. Hinchliffe, Chem. Phys. Lett. , 583 (1984). B. Skwara, W. Bartkowiak, A. Zawada, R. W. G´ora, andJ. Leszczynski, Chem. Phys. Lett. , 116 (2007). D. W. Schwenke, J. Phys. Chem. , 2867 (1996). C. J. Tainter, P. A. Pieniazek, Y. Lin, and J. L. Skinner, J.Chem. Phys. , 184501 (2011). G. R. Medders and F. Paesani, J. Chem. Theory Comput. ,1145 (2015). L. Lodi, J. Tennyson, and O. L. Polyansky, J. Chem. Phys. ,034113 (2011). Y. Wang, X. Huang, B. C. Shepler, B. J. Braams, and J. M.Bowman, J. Chem. Phys. , 094509 (2011). B. Guillot, J. Chem. Phys. , 1543 (1991). S. Habershon, T. E. Markland, and D. E. Manolopoulos, J.Chem. Phys. , 024501 (2009). “mbpol openmm plugin,” https://github.-com/paesanilab/mbpol openmm plugin, reference imple-mentation of the MB-pol water potential. Multicore and GPUimplementations for actual molecular dynamics production runswill become available in the near future. P. Eastman and V. Pande, Comput. Sci. Eng.12