Anomalous quantum and isotope effects in water clusters: Physical phenomenon, model artifact, or bad approximation?
AAnomalous quantum and isotope effects in water clusters: Physical phenomenon,model artifact, or bad approximation?
Sandra E. Brown and Vladimir A. Mandelshtam
Department of Chemistry, University of California, Irvine, CA 92697, USA
Free energy differences ∆ F := F − F prism are computed for several isomers of water hexamerrelative to the “prism” isomer using the self-consistent phonons method. We consider the isotopeeffect defined by the quantity δF D O := ∆ F D O − ∆ F H O , and the quantum effect, δF (cid:126) =0 :=∆ F (cid:126) =0 − ∆ F H O , and evaluate them using different flexible water models. While both δF D O and δF (cid:126) =0 are found to be rather small for all of the potentials, they are especially small for two of theempirical models, q-TIP4P/F and TTM3-F, compared to q-SPC/Fw and the two ab initio -basedmodels, WHBB and HBB2-pol. This qualitative difference in the properties of different water modelscannot be explained by one being “more accurate” than the other. We speculate as to whether theobserved anomalies are caused by the special properties of water systems, or are an artifact of eitherthe potential energy surface form/parametrization or the numerical approximation used. Changes in the properties of water upon isotopic sub-stitution provide clear evidence of the large influence ofnuclear quantum effects (e.g. zero-point motion, tunnel-ing, etc.) on the behavior of water. An oft-cited exam-ple is the shift in the melting point of water by 3.8 K indeuterated water, and 4.5 K in tritiated water [1]. Trendssuch as these have long suggested that nuclear quantumeffects destabilize the hydrogen bond network, “soften-ing” the structure of liquid water [2, 3]. However, withthe advent of improved water models, it has become in-creasingly clear that the situation is more complicatedthan may have been expected. Adding to the intrigueand complexity of the problem is the apparent competi-tion between quantum effects at play in hydrogen-bondedsystems [4], and in water in particular [5, 6], which canbe tipped one way or the other depending upon the tem-perature [4, 6]. Quantum fluctuations of the intramolec-ular bond stretching tend to strengthen hydrogen bondsas a result of larger monomer dipole moments, while in-termolecular quantum fluctuations tend to weaken them[4, 5]. Compared to the extensive body of work on quan-tum effects in bulk water systems, far fewer studies existexamining isotope effects in the water hexamer cluster,in spite of its recognition as a representative species forthree-dimensional hydrogen-bonded networks. The wa-ter hexamer is the smallest water cluster for which three-dimensional minimum energy structures are observed, asopposed to the two-dimensional ring structures favoredby smaller clusters, earning it the endearing distinctionof “smallest drop of water” and status as an importantbenchmark for understanding the structure and dynam-ics of water [7–10]. Aside from being an intrinsicallyinteresting physical problem, there is seemingly endlessmotivation for developing a complete understanding ofnuclear quantum effects in water, as isotope effects arerelevant even in the context of many large-scale biologicaland environmental processes [11, 12].A large number of water models have been proposed inthe past. However, partly due to numerical reasons (toexclude the fast intra-molecular degrees of freedom) mostof those were based on “rigid” water molecules. (The in- terested reader is referred to the review [13].) However,based on the previous studies [5, 6] we believe that aproper treatment of isotope effect must consider a flexi-ble water model. When working with empirical models awell-known concern is how well the model can reproduceproperties of water which were not used to parameterizeit. As far as we are aware, hardly any of the existingempirical models can be used to accurately predict anyspecific property of water, unless this very property wasactually used to design the model. At the same time,the current status of the electronic structure methodsdoes not appear to provide a sufficient level of accuracyfor an ab initio potential energy surface (PES) to havepredictive power. Moreover, even if we assume such amodel to be sufficiently accurate, the cost associated withthe evaluation of a high-quality ab initio
PES usuallymakes the inclusion of the quantum effects in a nucleardynamics simulation too expensive. In this context it isimportant to mention recent efforts in the constructionof accurate parametric fits of high-level ab initio waterPES’s, such as the WHBB [14] and HBB2-pol [15] po-tentials. In both WHBB and HBB2-pol, the monomercontribution is based on the same Partridge-Schwenkefit [16] used in the TTM3-F model. Most importantly,both WHBB and HBB2-pol include three-body terms(i.e., terms that explicitly couple coordinates of any threewater monomers), which are not present in the otherthree potentials. In addition, they are both permuta-tionally invariant in the sense that they allow for pro-ton exchange. While the HBB2-pol PES is substantiallyfaster and reportedly more accurate than WHBB [15], itis still much more expensive than most of the popular em-pirical models, including q-SPC/Fw [17], q-TIP4P/F [5],and TTM3-F [18]. Although it was not included in thepresent study, we note the recent development of MB-pol,a “first-principles” water potential which has been shownto be more accurate than HBB2-pol [19–21].In a recent paper [22] four of the above mentionedPES’s have been tested numerically using replica ex-change path integral molecular dynamics (RE-PIMD)simulations, specifically, by computing the free energy a r X i v : . [ phy s i c s . a t m - c l u s ] S e p differences for several isomers of water hexamer. Theconvergence of these RE-PIMD results is being ques-tioned elsewhere [23]. However, regardless of how accu-rate they are and regardless of how accurate any partic-ular water model is, our present focus is not so much onestablishing the absolute “truth” on whether or not thisor any other model correctly predicts that the cage iso-mer is energetically more favorable at low temperaturesthan the prism isomer, or vice versa. Rather, our pri-mary goal is to reveal possible generic properties of thesewater models which may be determined by the mannerin which they have been constructed and parametrized.In this study, in addition to the four PES’s consideredin Ref. 22, we also consider q-SPC/Fw [17], whose prin-cipal difference from q-TIP4P/F is the purely harmonictreatment of the OH-stretch. According to Habershon et al [5] this difference is responsible for the significantdifference in the magnitude of the isotope effect betweenthe two water models.All three empirical models, q-SPC/Fw, q-TIP4P/Fand TTM3-F, unlike the two ab initio -based PES’s, arenot permutationally invariant and include only two-bodyinteractions (i.e., interactions between not more than twowater monomers), except for the fact that TTM3-F ispolarizable. That is, by construction, the intramolec-ular forces in these models hold each hydrogen atomnearby its oxygen, while both WHBB and HBB2-pol al-low, in principle, for hydrogen exchange (albeit hardlycorrectly), thus leading to stronger couplings betweenthe intramolecular and intermolecular degrees of free-dom. With this said, it is not clear what the degree ofquasi-separability should be, and whether it is feasible toobtain a numerically practical parametrization that com-bines it correctly with the permutational invariance con-straint and a correct description of hydrogen exchange.The self-consistent phonons (SCP) method was firstproposed several decades ago as a means to include an-harmonic effects in the approximate treatment of the nu-clear dynamics of condensed phase systems [24, 25]. In-terest in SCP in the context of finite systems has emergedonly recently. In Ref. 26 it was used to compute the fun-damental frequencies of aromatic hydrocarbons, and inRefs. 27 and 28 for computing the ground states of verylarge Lennard-Jones clusters. Given a quantum many-body system localized in an energy minimum at thermalequilibrium, SCP maps it to a reference harmonic sys-tem by optimizing the free energy in the framework ofthe Gibbs-Bogolyubov variational principle. Recently wehave shown how to overcome the numerical bottleneck ofthe method, the accurate evaluation of multidimensionalGaussian averages of the potential and its derivatives[29]. This was done by implementing quasi-Monte Carlointegration which exhibits a superior quasi-linear scalingwith respect to the number of Monte Carlo points N MC ,compared to the √ N MC scaling of standard Monte Carlointegration. For the systems mentioned here, the incor-poration of quasi-Monte Carlo in SCP has been found toreduce the computational cost of the method by several orders of magnitude, expanding its applicability range to ab initio potentials.Though one can expect the accuracy of the SCPmethod to be significantly greater than that of the stan-dard harmonic approximation, the SCP method is itselfstill an approximation. The accuracy of SCP has beeninvestigated in Ref. [30] and Ref. [23] for the very caseof water hexamer. In Ref. [30] the calculation of theisomer energy differences at zero temperature (i.e. forthe isomer ground states), for which SCP turned out tobe accurate. In Ref. [23] we demonstrate that the clas-sical free energy differences (i.e., ∆ F (cid:126) =0 ) estimated bySCP agree well with those computed using a variant ofreversible scaling, an exact-in-principal method [31, 32].The quantities of interest are the free energies of var-ious isomers with respect to that of the prism isomer,ie., the differences ∆ F := F − F prism as a function oftemperature T . Fig. 1 shows the results for ∆ F for thecage and book isomers of classical (H O) and quantum(H O) and (D O) , for each of the five potentials spec-ified above. The calculations were carried out followingthe protocol described in Ref. 29, with the addition of arotational correction, described here in the appendix.While the absolute free energy for each isomer (notshown) was found to change dramatically and in a pre-dictable fashion, due to either isotopic substitution ofhydrogen in particular or to quantum effects in gen-eral, the effects are much smaller for the free energy dif-ferences. For all five potentials both the isotope shift δF D O := ∆ F D O − ∆ F H O and quantum shift δF (cid:126) =0 :=∆ F (cid:126) =0 − ∆ F H O are quite small compared to the ∆ F values, ie., neither isotopic substitution or even goingfrom classical water hexamer to quantum water hexamerwould change the energy ordering of the prism, cage, andbook isomers. A relatively small sensitivity of thermo-dynamic properties of water to isotopic substitution isthough a well established fact both experimentally andtheoretically (see, e.g., Refs. 5, 6). Arguably, the moststriking feature of these results is the much smaller quan-tum and isotope shifts for q-TIP4P/F and TTM3-F thanfor the other three potentials. The effect is actually sosmall that it is hardly possible to reproduce accuratelyusing a replica exchange path integral simulation due tostatistical errors [23].While a thermodynamic integration method, specif-ically that adapted to the path integral Monte Carloframework [33], would seem to be better suited, as it isdesigned to compute the isotope shift directly, very longsimulation times would still be needed in order to achievea sufficiently small statistical error. The SCP methodused in this work does not suffer from such convergenceproblems. Yet, before any reliable path integral simula-tions become available, which may or may not confirmthe present results, the manifestly approximate nature ofSCP requires some justification and discussion. In theabsence of any reliable path integral results, etc., to sup-port the (manifestly approximate) SCP results presentedhere, we consider the following two possibilities. First, in F - F ( p r i s m ) ( k ca l/ m o l ) H Oq-TIP4P/F 0 50 100 D OTTM3-F 0 50 100classicalWHBB0 50 100HBB2-pol 0 50 100q-SPC/Fw cage ca g e ca g eca g e book book book book book ca g e prism prism k B T FIG. 1: Free energy differences ∆ F := F − F prism of cage and book isomers of water hexamer with respect to the prism isomer,computed by SCP with rotational corrections ( cf. Eq. (A3)) for five different water potentials for both classical (H O) andquantum (H O) and (D O) . Note that the y-axis for q-SPC/Fw is shifted, as indicated by the dotted line corresponding tothe prism reference. spite of its approximate nature, the SCP estimates of thefree energy differences are actually accurate due to mas-sive cancellations of the systematic errors. Second, theobserved anomalies in the isotope effect are actually arti-facts of the inherent approximations of the SCP method.Since so far numerous numerical evidence is in favor ofthe SCP method (e.g., Refs. 23, 30), in the following dis-cussion we disregard the latter possibility (unless provedto be the case) as being not constructive.Supposing the SCP results to be correct, the questionis why the quantum and isotope effects for q-TIP4P/Fand TTM3-F are substantially smaller than for theother three potentials. According to the analysis ofHabershon et al. [5], the absence of anharmonic termsin the q-SPC/Fw description of the OH stretch couldaccount for a larger isotope effect than that seen in q-TIP4P/F. The argument then is that it is the an-harmonic terms in the OH stretch that are responsiblefor cancellations of the competing quantum effects inthe q-TIP4P/F water, resulting in a small isotope ef-fect. Since TTM3-F is also highly anharmonic, this sameargument can be used to explain the very small isotopeand quantum effects observed with this PES. However, inboth WHBB and HBB2-pol, the monomer contribution isbased on the same Partridge-Schwenke fit [16] used in theTTM3-F model, yet the larger isotope/quantum shifts forWHBB and HBB2-pol are comparable to those seen withq-SPC/Fw. As such, it appears that the relatively largeisotope effects in both ab initio potentials must be dueto their explicit inclusion of three-body terms. Finally,we note again the constraint of permutational invariance,present only in the parametrization of the ab initio po-tentials, which, in principle, allows for hydrogens to betransferred between oxygens. Neither q-TIP4P/F norTTM3-F are permutationally invariant, i.e., each hydro-gen remains assigned to its oxygen in these models. Still,the inclusion of this property does not guarantee thathydrogen exchange is accounted for accurately by eitherWHBB or HBB2-pol. In the context of this work, thekey implication of this property is the much greater cou-pling via intermolecular degrees of freedom in the ab ini-tio potentials. In other words, the quasi-separability ofthe intermolecular and intramolecular degrees of freedompresent in q-TIP4P/F and TTM3-F models is much lesspronounced in WHBB and HBB2-pol, leading to greaterquantum effects associated with much more flexible hy-drogen degrees of freedom. Acknowledgements
This work was supported by the National ScienceFoundation (NSF) Grant No. CHE-1152845. SEB waspartially supported by NSF Grant No. DMS-1101578.Volodymyr Babin and Francesco Paesani are acknowl-edged for discussing with us their results on water hex-amer and for sending us the source code for the HBB2-polPES.
Appendix A: Free energy within the quasi-harmonicapproximation with rotational correction.
Consider an N -atom cluster and its isomer correspond-ing to a relatively deep and stable potential energy mini-mum, i.e., we assume that it is separated from the rest ofthe configuration space by relatively large energy barri-ers. In the absence of an external field the translations ofthe center of mass can be separated so that we may con-sider the subspace R (3 N − that includes only the vibra-tional degrees of freedom and the rotations of the wholecluster. Because the potential energy U ( r ) is invariant tothese rotations, the energy minimum is a (3)-dimensionalmanifold.To further simplify the problem consider a harmonicapproximation or, more generally, a quasi-harmonicapproximation, represented by an effective (generallytemperature-dependent) harmonic Hamiltonian. To be concrete consider the quasi-harmonic approximation aris-ing within the SCP method [29], as used in this work.SCP considers the system in the so called “Eckart sub-space”, which is a reduced (3 N − H h ( T ) define the (generally temperature-dependent) ef-fective harmonic Hamiltonian, and ω k ( k = 1 , ..., N − F ( T ) ≈ (cid:104) U (cid:105) h + (cid:88) k (cid:34) k B T log (cid:18) (cid:126) ω k k B T (cid:19) − (cid:126) ω k (cid:126) ω k k B T (cid:35) , (A1)where (cid:104) U (cid:105) h defines the thermal average of the originalpotential over the reference harmonic system.The above approximation to the free energy completelyignores the rotational contribution, which may be impor-tant for small enough clusters. In order to include it wepropose the use of a rigid asymmetric top correction [34],which with the omission of terms that cancel when theenergy difference between isomers is considered, reducesto F rot ( T ) ≈ − k B T I I I ( k B T ) (cid:126) + k B T log σ , (A2)where I , I and I are the principal moments of inertia ofthe isomer evaluated at its minimum configuration, and σ is the order of the isomer point group, which is unityunless the isomer configuration has symmetries. Conse-quently, the rotational correction to the free energy dif-ference for two isomers A and B can be estimated using∆ F rot ≈ k B T (cid:20)
12 log I B I B I B I A I A I A + log σ A σ B (cid:21) . (A3)Note that for water hexamer the rotational contributionto the free energy difference between between differentisomers is of the order of ∼ . k B T . [1] D. R. Lide, CRC Handbook of Chemistry andPhysics, 85th Edition (Taylor & Francis, 2004),ISBN 9780849304859, URL http://books.google.com/books?id=WDll8hA006AC .[2] T. F. Miller and D. E. Manolopoulos, The Jour-nal of Chemical Physics , 154504 (2005), URL http://scitation.aip.org/content/aip/journal/jcp/123/15/10.1063/1.2074967 .[3] J. A. Morrone and R. Car, Phys. Rev. Lett. , 017801 (2008), URL http://link.aps.org/doi/10.1103/PhysRevLett.101.017801 .[4] X.-Z. Li, B. Walker, and A. Michaelides, Proceedingsof the National Academy of Sciences .[5] S. Habershon, T. E. Markland, and D. E. Manolopou-los, The Journal of Chemical Physics , 024501 (pages 11) (2009), URL http://link.aip.org/link/?JCP/131/024501/1 .[6] T. E. Markland and B. J. Berne, Proceedings ofthe National Academy of Sciences .[7] J. K. Gregory and D. C. Clary, The Jour-nal of Physical Chemistry , 18014 (1996),http://pubs.acs.org/doi/pdf/10.1021/jp9616019, URL http://pubs.acs.org/doi/abs/10.1021/jp9616019 .[8] S. S. Xantheas, C. J. Burnham, and R. J. Har-rison, The Journal of Chemical Physics , 1493(2002), URL http://scitation.aip.org/content/aip/journal/jcp/116/4/10.1063/1.1423941 .[9] C. P´erez, M. T. Muckle, D. P. Zaleski, N. A.Seifert, B. Temelso, G. C. Shields, Z. Kisiel,and B. H. Pate, Science .[10] R. J. Saykally and D. J. Wales, Science .[11] F. Paesani and G. A. Voth, The Journalof Physical Chemistry B , 5702 (2009),http://pubs.acs.org/doi/pdf/10.1021/jp810590c, URL http://pubs.acs.org/doi/abs/10.1021/jp810590c .[12] J. Hoefs,