Binding Energies from Diffusion Monte Carlo for the MB-pol H_2O and D_2O Dimer: A Comparison to Experimental Values
BBinding Energies from Diffusion Monte Carlo for the MB-pol H O and D ODimer: A Comparison to Experimental Values.
Joel D. Mallory and Vladimir A. Mandelshtam Department of Chemistry, University of California, Irvine, 1102 Natural Sciences II Irvine, California 92697,USA
The Diffusion Monte Carlo (DMC) method is applied to compute the ground state energies of the watermonomer and dimer and their D O isotopomers using MB-pol; the most recent and most accurate ab inito -based potential energy surface (PES). MB-pol has already demonstrated excellent agreement with high levelelectronic structure data, as well as agreement with some experimental, spectroscopic, and thermodynamicdata. Here, the DMC binding energies of (H O) and (D O) agree with the corresponding values obtainedfrom velocity map imaging within, respectively, 0.01 and 0.02 kcal/mol. This work adds two more valuabledata points that highlight the accuracy of the MB-pol PES.The development of a full-dimensional potential energysurface (PES) for a many-body system that extends toprogressively larger cluster sizes and, at the same time,is computationally feasible has been a longstanding issuein electronic structure theory . Many different empir-ical water PESs have been parametrized ranging fromfully coarse-grained to flexible atomistic models that in-clude polarization effects and charge transfer (see, e.g.,Refs. 2–9). Nevertheless, empirical PESs are generally in-adequate for capturing the behavior of water for a widearray of cluster sizes (i.e., from small clusters to the bulkliquid) and thus, have generated ample motivation for theconstruction of ab initio and ab inito -based surfaces withthe latter having a foundation in the many-body expan-sion of the interactions . Several notable PESs belong-ing to this family are DPP2 , CC-pol , the HBB0-2series for the water dimer , WHBB , and HBB2-pol . Along the same vein, the MB-pol PES hasmost recently emerged as an ab inito -based water surfacerigorously derived from the many-body expansion of theinteraction energy and expressed in terms of explicit one-,two-, and three-body contributions, with all higher-orderterms being represented by (classical) many-body induc-tion within a modified version of the polarization modeloriginally employed by the TTM4-F potential . Simi-larly to WHBB and HBB2-pol, the MB-pol one-, two-,and three-body terms were obtained from fits to largesets of CCSD(T) monomer, dimer, and trimer energiescalculated in the complete basis set limit.Paesani and co-workers (see Refs. 19–21) have demon-strated that MB-pol effectively reproduces high-level abinitio data for: (a) stationary point energies on the(H O) and (H O) PESs; (b) PESs for (H O) and(H O) interaction energies plotted versus O-O distanceand O-O-O angle, respectively; and (c) (H O) inter-action energies and (H O) − isomer energies. Like-wise, results from MB-pol were shown to exhibit goodagreement with experimental data for (H O) vibration-rotation tunneling splittings , as well as the structural,thermodynamic and dynamical properties of bulk waterat a fully quantum-mechanical level . Both infrared andRaman spectra calculated from centroid molecular dy-namics simulations with the MB-pol potential were found to be in good agreement with the corresponding experi-mental results. .In this paper, we use the diffusion Monte Carlo (DMC)method to compute the MB-pol (H O) and (D O) binding energies and compare them to the correspond-ing experimental values . Note here that Ref. 15reported a DMC result for the binding energy of theHBB2 water dimer, (H O) , that was later confirmedexperimentally .DMC uses a population of N W random walkers thatsample the configuration space bounded by a PES andcollectively represent the wavefunction of the many-bodysystem at projection time τ for each time step of length∆ τ . At sufficiently long τ the distribution of ran-dom walkers becomes stationary and the instantaneousenergy E ref ( τ ) fluctuates about its average value (cid:104) E ref (cid:105) .In the limit of τ → ∞ , ∆ τ →
0, and N W → ∞ , thisdistribution converges to the ground state wavefunctionwith (cid:104) E ref (cid:105) = E , the ground state energy. In a re-cent paper , we have undertaken a thorough analysisof the behavior and extent of the bias (systematic error)arising from the finite time step ∆ τ , as well as the biascaused by a finite random walker population N W for theq-TIP4P/F water monomer, dimer and hexamer. Thetime step bias in the estimate of E vanishes slowly, andfor the water monomer at ∆ τ = 10 . .
015 kcal/mol.However, this bias cancels nearly completely when theground state energy differences, such as the binding en-ergy, D := 2 E H O − E (H O) , (1)or the isotope shift, δD := D (D O) − D (H O) , (2)are computed. (Note that the same value of ∆ τ = 10 . .) The bias in N W also cancels for theisotope shift ( δD ), but for the binding energy ( D ) itdoes not cancel. This is because the DMC ground stateenergy estimate for the monomer converges much fasterwith respect to N W than that for the dimer. Thus, inorder to obtain an accurate estimate of the asymptoticvalue of D at N W → ∞ , one needs to perform a series of a r X i v : . [ phy s i c s . c h e m - ph ] S e p calculations using several different values of N W . It wasconcluded that 0 .
002 kcal/mol accuracy for the dimerbinding energy D could be achieved using N W ∼ . × .In this work, we compute the energies for the MB-polwater monomer and dimer adhering to the DMC protocolof Ref. 28 and performing studies for the time step with∆ τ = 5, 10, and 15 au and for the random walker popula-tion with N W = 1 . × , 3 . × , and 1 . × . All ofthe DMC simulations were run to a maximum projectiontime of 2 . × au. In the time step studies, the walkerpopulation was fixed at N W = 1 . × , while a timestep of ∆ τ = 10 . E (∆ τ ), in the ∆ τ → τ = 0. The latter is because the errorcaused by the split-operator approximation implementedin the DMC method is quadratic in ∆ τ . Consequently,since three different values of the time step ∆ τ have beenused, we consider the cubic interpolation: E (∆ τ ) ≈ A + B ∆ τ + C ∆ τ , (3)using three fitting parameters, with A giving the groundstate energy estimate.Here, we also comment that our walker number stud-ies employ larger values of N W than might seem to beneeded, as including more walkers in the population si-multaneously reduces the statistical uncertainty as wellas the systematic error .Fig. 1 and Fig. 2 show the bias in ∆ τ and N W forthe ground state energy estimates from DMC ( E ) asfunctions of the time step and the inverse walker popula-tion 1 /N W , respectively, as well as the associated bind-ing energies ( D ). The curves follow a pattern similarto that reported in Ref. 28 for the q-TIP4P/F dimerand monomer. A clear bias does exist for the groundstate energies ( E ) in both ∆ τ and N W , but for the re-ported range of time steps and walker numbers, the en-ergy change is in the second or third positions beyondthe decimal. Moreover, the binding energy bias in ∆ τ is virtually negligible due to error cancellations. This isbecause the largest contribution to the time step errorcomes from the intramolecular degrees of freedom, whichare essentially the same regardless of whether the sys-tem is comprised of a single water molecule or multiple,interacting water molecules. Because the intramolecularmodes are nearly invariant to changes in the number andspatial orientation of the water molecules in clusters, thebehavior and extent of the E bias curves in ∆ τ for theH O and D O monomer and dimer are very similar (seethe top four panels of Fig. 1), which leads to substan-tial elimination of the bias such that the binding energydisplays only a weak dependence on ∆ τ . On the otherhand, the bias in 1 /N W persists for the binding energyeven after the energy differences are taken. In this case,the E monomer bias in N W disappears much faster thanthat of the dimer (note the difference in scales between the monomer (top) and dimer (middle) panels of Fig.2), thereby causing imperfect cancellation of error andthe appearance of a strong residual bias in the bindingenergy.Table. I shows the DMC binding energies extrapo-lated to the N W → ∞ limit for three different PESs ex-amined in this work: q-TIP4P/F , TTM3/F , and MB-pol. The statistical errors for the binding energieswith the largest walker population used in this study of N W = 1 . × are all known to be on the order of10 − kcal/mol, but these numbers are not displayed inthe table because extrapolation of the error bar valueshas proven to be unreliable. Here, the same DMC proce-dure was used to calculate the ground state energies forthe TTM3/F H O and D O monomer and dimer as thosewith q-TIP4P/F and MB-pol. Additionally, we also pro-vide the (H O) TTM3/F and HBB2 binding energies asreported in Ref. 15. The D values for the q-TIP4P/Fdimers are ∼ . . This discrepancy is not surprising consid-ering that the q-TIP4P/F PES was not designed to rep-resent the microscopic properties of small water clustersaccurately. Our TTM3-F D value for (H O) agrees wellwith that reported by Bowman and co-workers with adifference of only 0 .
06 kcal/mol, but our results indicatethat the binding energy for this PES is still off from theexperimental values by 0 .
62 and 0 .
53 kcal/mol for (H O) and (D O) , respectively.Notably, the (H O) DMC binding energy, also fromthe ab initio -based two-body PES, HBB2 , (with valuesof ∆ τ and N W similar to those used in the present work)is highly accurate upon comparison to the experimentalresults. At the same time, we highlight that the MB-polPES likewise shows nearly perfect agreement with thevelocity map imaging D values for both MB-pol dimers,(H O) and (D O) . For the MB-pol PES, the bindingenergies differ from the experimental results by only 0 . O) (as is also the case for the HBB2dimer PES) and 0 .
02 kcal/mol for (D O) . Such an excel-lent correspondence adds two more valuable data pointswhich underscore a favorable assessment for the accuracyof the MB-pol PES.Therefore, efforts are currently underway to establishthe true ground state energy and wavefunction of theMB-pol water hexamer using DMC. Additionally, in ac-cordance with Ref. 28, the ground state energies for theisomers of the MB-pol hexamer will be computed. ACKNOWLEDGEMENTS
This work was supported by the National ScienceFoundation (NSF) Grant No. CHE-1152845. We thankSandra Brown for useful discussions, and Francesco Pae-sani and his co-workers for sharing the MB-pol code withus.
FIG. 1. DMC energies for the MB-pol water monomer and dimer as a function of the time step ∆ τ . The time steps used inthis study were ∆ τ = 5 .
0, 10 .
0, and 15 . N W = 1 . × , and the total projectiontime for all runs was 2 . × au. Data points were interpolated using a cubic polynomial fit with a vanishing derivative at∆ τ →
0. Top left: H O monomer energies. Top right: D O monomer energies. Middle left: H O dimer energies. Middle right:D O dimer energies. Bottom left: H O dimer binding energy. Bottom right: D O dimer binding energy. E ( k ca l/ m o l ) E ( k ca l/ m o l ) ∆τ D ( k ca l/ m o l ) ∆τ (H O) H O (D O) D O2*H O-(H O) O-(D O) TABLE I. DMC binding energy D for the H O and D O dimer with different PESs. The DMC binding energies computedin this work for a time step of ∆ τ = 10 . N W → ∞ limit. Binding energies marked with the “a”superscript were reported in ref. . The experimental binding energies (Expt.) are from refs. . The error bar magnitudesfor this work are on the order of 10 − . All energies are reported in kcal/mol. D Structure q-TIP4P/F TTM3/F MB-pol HBB2 Expt.(H O) . ± . a . ± . a . ± . O) . ± . G. Chalasinski and M. M. Szczesniak, Chem. Rev. , 4227(2000). L. X. Dang and B. M. Pettitt, J. Phys. Chem. , 3349 (1987). G. S. Fanourgakis and S. S. Xantheas, J. Phys. Chem. A ,4100 (2006). G. S. Fanourgakis and S. S. Xantheas, J. Chem. Phys. ,074506 (2008). C. Burnham, D. Anick, P. Mankoo, and G. Reiter, J. Chem.Phys. , 154519 (2008). V. Molinero and E. B. Moore, J. Phys. Chem. B , 4008 (2008). S. Habershon, T. E. Markland, and D. E. Manolopoulos, J.Chem. Phys. , 024501 (2009). C. Vega and J. L. Abascal, Phys. Chem. Chem. Phys. , 19663(2011). A. J. Lee and S. W. Rick, J. Chem. Phys. , 184507 (2011). J. E. Mayer and M. G. Mayer,
Statistical Mechanics ; John Wiley& Sons Inc. Hoboken, NJ, 1940.
FIG. 2. Same as Fig. 1 but as a function of reciprocal walker number 1 /N W . The walker populations used in this study were N W = 1 . × , 3 . × , and 1 . × . The time step was fixed at ∆ τ = 10 . E ( k ca l/ m o l ) E ( k ca l/ m o l ) × -4 × -4 × -4 W D ( k ca l/ m o l ) × -4 × -4 × -4 W (H O) H O (D O) D O2*H O-(H O) O-(D O) R. Kumar, F.-F. Wang, G. R. Jenness, and K. D. Jordan, J.Chem. Phys. , 014309 (2010). R. Bukowski, K. Szalewicz, G. C. Groenenboom, and A. Van derAvoird, Science , 1249 (2007). X. Huang, B. J. Braams, and J. M. Bowman, J. Phys. Chem. A , 445 (2006). X. Huang, B. J. Braams, J. M. Bowman, R. E. Kelly, J. Ten-nyson, G. C. Groenenboom, and A. van der Avoird, J. Chem.Phys. , 034312 (2008). A. Shank, Y. Wang, A. Kaledin, B. J. Braams, and J. M. Bow-man, J. Chem. Phys. , 144314 (2009). Y. Wang, X. Huang, B. C. Shepler, B. J. Braams, and J. M.Bowman, J. Chem. Phys. , 094509 (2011). V. Babin, G. R. Medders, and F. Paesani, J. Phys. Chem. Lett. , 3765 (2012). G. R. Medders, V. Babin, and F. Paesani, J. Chem. TheoryComput. , 1103 (2013). V. Babin, C. Leforestier, and F. Paesani, J. Chem. TheoryComp. , 5395 (2013). V. Babin, G. R. Medders, and F. Paesani, J. Chem. TheoryComp. , 1599 (2014). G. R. Medders, V. Babin, and F. Paesani, J. Chem. TheoryComp. , 2906 (2014). G. R. Medders and F. Paesani, J. Chem. Theory Comp. , 1145(2015). G. R. Medders and F. Paesani, J. Chem. Phys. , 212411(2015). B. E. Rocher-Casterline, L. C. Ch’ng, A. K. Mollner, andH. Reisler, J. Chem. Phys. , 211101 (2011). L. C. Chng, A. K. Samanta, G. Czako, J. M. Bowman, andH. Reisler, J. Am. Chem. Soc. , 15430 (2012). J. B. Anderson, J. Chem. Phys. , 1499 (1975). J. B. Anderson, J. Chem. Phys. , 4121 (1976). J. D. Mallory, S. E. Brown, and V. A. Mandelshtam, J. Phys.Chem. A , 6504 (2015). Y. Wang, V. Babin, J. M. Bowman, and F. Paesani, J. Am.Chem. Soc.134