Accurate Modeling of Water Clusters with Density-Functional Theory

Subscriber access provided by Georgetown University | Lauinger and Blommer Libraries

Article

Accurate modeling of water clusters with densityfunctional theory using atom-centered potentials Jake D. Holmes, Alberto Otero-de-la-Roza, and Gino A. DiLabio J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00624 • Publication Date (Web): 11 Aug 2017 Downloaded from http://pubs.acs.org on August 14, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Accurate modeling of water clusters with density-functional theory using atom-centered potentials Jake D. Holmesa ID, Alberto Otero-de-la-Rozaa ID, and Gino A. DiLabioa,b* ID a

Department of Chemistry and b Faculty of Management, The University of British Columbia, 3247 University Way, Kelowna, British Columbia, Canada V1V 1V7 Abstract The ability of atom-centered potential (ACPs) to improve the modeling of water clusters using density-functional methods is explored. Water-specific ACPs were developed using accurate ab initio reference data to correct the deficiencies of the BHandHLYP density functional in the calculation of absolute and relative binding energies of water clusters. In conjunction with aug-cc-pVTZ basis sets and with or without dispersion corrections, it is possible to obtain absolute binding energies for water clusters containing up to 10 H2O molecules to within 0.44 kcal/mol, or 0.04 kcal/mol per water molecule. In contrast, dispersion-corrected BHandHLYP/aug-cc-pVTZ predicts binding energies with errors as large as 6 kcal/mol for (H2O)10 in absence of ACPs. Therefore, the ACPs improve predicted binding energies in these clusters by more than an order of magnitude. The conformers of (H2O)16 and (H2O)17 were used to validate the application of ACPs to larger clusters. ACP-based approaches are able to predict the binding energies in (H2O)16,17 within a range of 0.3 – 2.2 kcal/mol (or to within 1.3%) of recently revised ab initio wavefunction results. ACPs for basis sets smaller than aug-cc-pVTZ are also presented. However, the ability of the BHandHLYP/ACP approach to predict accurate binding energies deteriorates as the size of the basis sets decreases. Nevertheless, ACPs improve predicted binding energies by as much as a factor of 50 across the range of the basis sets studied. The BHandHLYP/aug-cc-pVTZ-ACP method was applied to (H2O)25 in order to identify the minimum-energy structure of a collection of proposed global minimum-energy structures. The BHandHLYP/aug-cc-pVTZ-ACP approach is an accurate and computationally affordable alternative to wavefunction theory methods for the prediction of the binding energies and energy ranking of water clusters.

ID

ORCID ID Jake D. Holmes: 0000-0003-1142-2535 A. Otero-de-la-Roza: 0000-0002-4866-5816 Gino A. DiLabio: 0000-0002-3778-3892

*

Author of correspondence. E-mail: [email protected]. Phone: +1-250-807-8617.

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction The importance of water in chemistry and biology has made it the subject of many computational studies. Recent works have focused on determining accurate binding energies (BE). for water clusters using ab initio wavefunction theory methods,1,2,3,4,5 and finding global minima of water clusters based on force field methods6,7 Monte Carlo methods,8,9 and ab initio approaches.4,10 First principles modeling of water is challenging because of the complex non-covalent interactions between molecules.11 For smaller clusters, the interplay between van der Waals and hydrogen bonding interactions can be accurately determined using currently available computational methods, such as second-order Møller-Plesset perturbation theory (MP2), provided that the electronic energies are extrapolated to the complete basis set (CBS) limit and higher-order electron correlation effects (e.g. from a CCSD(T) treatment) are incorporated12. A comparison of the water cluster BEs reported in references 1 and 2 to those obtained by Gordon’s group,12 for example, underscores the need for complete basis set extrapolation: The latter utilized CCSD(T) with augmented triple-zetaquality basis sets and, as a consequence, the predicted BEs are larger than those determined with CBS extrapolations by ~1 kcal/mol for the prism conformer of the water hexamer. Since large water clusters exhibit significant variations in structure within a sub-kcal/mol energy window, 1,5,10 the ability to determine very accurate binding energies is critical for global energy minima determination in water clusters. Although the ongoing development of computing technology and of ab initio wavefunction methods has extended the size of water clusters that can be simulated, the high computational cost associated with wavefunction approaches in general limits their applicability to fairly small clusters. Currently, the largest water cluster treated by ab initio techniques consists of 25 water molecules.10 A number of recent works have examined the applicability of density-functional theory (DFT) to determine the properties of water clusters. The computational cost of DFT methods scales formally as or with system size, much less than the (MP2) or (CCSD(T)) scaling associated with wavefunction methods. Consequently, DFT-based approaches have the potential to treat much larger water clusters than correlated wavefunction techniques. We note, however, that recent work by Valeev and Neese on domain based local pair natural orbital coupled cluster techniques have been shown to scale far more favorably than conventional CCSD(T) approaches and consequently the former approach could be applied to much larger clusters.13 Several recent studies have examined the ability of conventional DFT-based 2 ACS Paragon Plus Environment

Page 2 of 31

Page 3 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


approaches to predict the BEs of water clusters, with a particular focus on water hexamers. Multiple hexamer conformations lie within a small BE range, with several of the structures representing a cross-over from two-dimensional to three-dimensional hydrogen bonding networks. Dahlke et al. showed that a number of DFT-based methods are able to predict the geometries of six hexamers in very good agreement those obtained by MP2 with triple-zeta-quality basis sets. However, the mean errors in the BEs of five of the conformers relative to the “prism” (i.e. lowest energy) structure were shown to span a large range, from ca. 0.3 to 2.0 kcal/mol.14 Santra et al.15 studied four conformers of the water hexamer with the goal of assessing the performance of a large number of conventional DFT-based methods, and concluded that none of the tested approaches were capable of predicting the correct lowest-energy conformation. Ultimately, the authors demonstrated that the erroneous description of non-covalent interactions by the functionals they studied was responsible for the poor description of the hexamers. Santra et al. also showed that the inclusion of pair-wise dispersion corrections to the DFT energies resulted, in some cases, in a significant improvement in the predicted energy ordering of the water hexamer conformers. Several approaches to improve the accuracy of DFT for non-covalent interactions using dispersion corrections have been suggested, as discussed in recent reviews.16,17,18 Early efforts focused on the use of empirical, pair-wise dispersion parameters in conjunction with density functionals,19,20 an approach that was later popularized by Grimme’s group.21,22,23 An elegant and accurate pair-wise dispersion correction approach was developed by Becke and Johnson, whose method calculates dispersion coefficients from the dipole moment of the exchange hole.24,25,26 Zhao and Truhlar developed a series of “Minnesota” functionals, whose parameterizations cover noncovalent interactions.27,28,29 DiLabio has developed an atom-centered potential (ACP) approach to improve the ability of DFT-based methods in the treatment of noncovalently interacting systems (dispersion-correcting potentials, DCPs).30,31,32 DCPs are combined with a functional and basis set combination, and fit to minimize the errors in predicted non-covalent BEs. The DCP approach is similar to the DCACP technique introduced by von Lilienfeld and co-workers.33,34,35 Although all of these dispersion-corrected DFT methods were developed to improve the description of non-covalent interactions in general, they have shortcomings for particular systems, such as water clusters. For example, Santra et al. showed that the inclusion of empirical pair-wise dispersion corrections to certain functionals did not allow for the correct prediction of water hexamer conformer ordering.15 Rothlisberger’s group applied DCACPs in conjunction with the BLYP functional to four conformations of 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 31

the water hexamer36, and obtained reasonable relative BEs, but the absolute BEs were about 4 kcal/mol too high on average compared to the high-level values given in reference 1. The improved DCACPs of Karalti et al. give significantly improved results, with errors for the hexamer binding energies predicted to be within about 1.5 kcal/mol of the reference values.37 Likewise, the inclusion of dispersion corrections via the XDM scheme resulted in very good relative BEs for the water hexamer, but absolute errors for the prism ranged from ca. 0.6 to 4.4 kcal/mol, depending on the underlying functional.38 Truhlar’s group showed that the M06-2X functional predicts the relative BEs of five hexamer conformers with a mean unsigned error of 1.4 kcal/mol relative to the “prism” structure. In a more recent study, the Truhlar group found that M06-2X gives BEs for five conformers of the water 16-mer with errors ranging from 0.99 to 2.22 kcal/mol.39 However, it should be noted the reference CCSD(T)/aug-cc-pVTZ BEs used for the water 16-mer do not contain corrections for CBS extrapolations, which introduces some degree of uncertainty in their values.4 In fact, Miliordos and Xantheas recently downwardly revised the BEs in (H2O)16 and (H2O)17 by about 6 kcal/mol,5 thus indicating that M06-2X in fact performs rather poorly for the BEs in these conformers. As indicated in the recent thorough review by Gillan et al.,40 no DFT-based method can accurately predict the BEs in water clusters over all cluster size ranges. The approach we propose in this work uses DFT combined with ACPs, an approach that has been previously used to incorporate dispersion effects30,31,41 and also to mitigate the effects of basis set incompleteness.31,32,41,42, 43,44 ACPs also correct for the underlying deficiencies in certain DFT-based methods to predict molecular and cluster properties (see, for instance, the recent study by Torres and DiLabio44 on methane clusters). In this work, we describe the development of water-specific ACPs for accurately predicting BEs and structures of water clusters. We seek to develop a computational approach that has accuracy on par with CCSD(T)/CBS but at the computational cost of conventional DFT. We will outline how ACPs were fitted to highlevel ab initio data for a large number of (H2O)n for n=2,3,4 clusters to reproduce the individual two-, three-, and four-body interactions that occur in water clusters. The binding energies and structures of relatively stable water clusters containing up to six molecules are also included in the fit. The ACPs were developed for use with the BHandHLYP functional, and validated against (H2O)n clusters with n=7-10 and on (H2O)16 and (H2O)17 clusters, for which recently revised highly-correlated wavefunction BE data are available. As an example application, the newly developed ACPs were then applied to (H2O)25 clusters to ascertain the identity of the global minimum for this cluster size. 4 ACS Paragon Plus Environment

Page 5 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Theoretical and Computational Details Development of Atom-Centered Potentials. ACPs are atom-centered, one-electron potentials with the following form:

= ∑ + ∑ ∑ | |

where

= ∑

% !"# $

(1)

(2)

In Eqs. 1 and 2, A represents the atoms on which the potentials are centered (in this work oxygen and hydrogen), rA is the distance from the nucleus, l and L are angular momentum terms, and | are spherical harmonics. The coefficients and exponents & of the Gaussian function in Eq. 2 are empirical parameters that are determined through least-squares fitting. The forms of Eq. 1 and 2 are identical to that of standard effective core potentials,45 but no electrons are replaced.

The energy contribution due to the application of the potential, , associated with each

is: () , & = ∑+*+ ||*+

(3)

where *+ are the orbitals from the self-consistent solution of the Kohn-Sham equations. ACPs are developed by optimizing the coefficient and exponent of the Gaussian function in Eq. 2 for each atom and angular momentum channel in order to obtain an energy correction (Eq. 3) that yields the desired properties (in this case, BEs of water clusters) in agreement with high-level reference data. All angular momenta up to the maximum l for the corresponding basis set are used, including the local channels (, ). For instance, in combination with the aug-cc-pVTZ, we design ACPs that have local, s, p, and d channels for oxygen, and local, s, and p for hydrogen. Twenty exponents were selected for use in the ACP terms: 0.004, 0.0075, 0.01, 0.02, 0.04, 0.08 – 0.28 in 0.02 increments, and 0.40 – 1.00 in 0.2 increments. These exponents were chosen so that the resulting Gaussian functions cover a range of distances from the nuclei that is optimal to correct for non-covalent interaction energies in water clusters.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The coefficients in the ACPs were determined by a least-squares fitting procedure, which has been described in detail elsewhere.43 In short, the effect of each ACP term in Eq. 3 is calculated, and the coefficients are determined by a least-squares fit in which the ACP energy contribution is fit to a training set containing an extensive list of target properties (binding energies, dipoles, three- and four-body energies, etc.) The least-squares fitting procedure gives reliable ACPs only if the coefficients are small enough that the ACP does not excessively perturb the ground-state wavefunctions. This requires using some sort of constraint on the least-squares fit. In this work, we precomputed the maximum coefficient values for each atom and angular momentum channel such that, provided the ACP coefficients are below these values, the disagreement between the least-squares and the self-consistent ACP calculations is less than ca. 0.001 kcal/mol. For more details on the fitting procedure, the reader is directed to reference 43. Based on the good performance of its dispersion-corrected variants for water (see below), we chose to develop ACPs for the BHandHLYP functional. The fitting set comprises two kinds of systems. First, we use a fraction of the water clusters presented in the work of Temelso et al.,1 specifically those formed by (H2O)n with n=2-6. The second component of the fitting set is designed to ensure that the 2-, 3-, and 4-body interactions that occur in large clusters of water are well-reproduced.43 With this idea in mind, we took all clusters in the Temelso et al. set, i.e. (H2O)n with n=2-10 and the four (H2O)20 models proposed by Bryantsev et al.,3 and constructed the geometries of all dimers, trimers, and tetramers of water that can be formed by taking all combinations of two, three, and four molecules, respectively. The number of dimers, trimers, and tetramers resulting from this partitioning (1409, 5689, and 20720, respectively) is too large to be practical for ACP fitting. To reduce the number of data points, we consider the n-body energies (-( as a function of the square root of the sum of the squared intermolecular distances ./01 . Plots of -( vs. ./01 are shown in the Supporting Information (SI, Fig. S1). The -( vs. ./01 space was uniformly partitioned and the same number of points in each partition was chosen if available. This procedure resulted in a subset of fitting data containing 410 dimers, 580 trimers, and 58 tetramers that spans all types of interactions (both attractive and repulsive) and intermolecular distances. The relatively small number of tetramers (58) was chosen due to the small energy contribution of the 4-body energies to the total binding energy.46 The entire fitting set, which is provided in the SI, contains 1069 BE and corresponding geometry data points. The two-, three-, and four-body energy contributions to BEs associated with each 6 ACS Paragon Plus Environment

Page 6 of 31

Page 7 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of these fragments were calculated using RI-MP2 and a complete-basis-set (CBS) extrapolation using the 4-5 inverse polynomial formula47 with aug-cc-pVNZ, N=T,Q,5 (using the TURBOMOLE program48). Each contribution was modified by a correction for higher-level correlation by including a term for the difference between the CCSD(T)/augcc-pVDZ and MP2/aug-cc-pVDZ energies. In reference 1, this correction is described as 7789: 7789: 3456 . Temelso et al. demonstrated that the 3456 term contributes to water cluster BEs to varying extents depending on the size and specific structure of the cluster. 7789: For example, they showed that for three conformers of (H2O)4, the 3456 term ranges from about -0.08 to +0.22 kcal/mol. The plots shown in Figure 1 demonstrate similar 7789: variability in the 3456 terms calculated for the 2-, 3- and 4-body interaction energies and underscore the importance higher-order correlation effects in water clusters. For 7789: example, for the 2-body interaction terms, the 3456 term spans the range of ca -0.1 to +0.35 kcal/mol (Fig. 1A). Considering all n-body interactions (Fig. 1) suggests that higher-order correlation effects can contribute as much as 10% to the total BE. A



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 31

B

C

7789:

Figure 1. Plots of the 3456

term for the (A) 2-, (B) 3-, and (C) 4-body energies in the

cluster fragments used in our ACP fitting set, against the total CCSD(T)/CBS value of the 2-body, 3-body, 4-body BE contribution, respectively. 8 ACS Paragon Plus Environment

Page 9 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A second important conclusion that can be derived from Figure 1 is that the magnitudes 7789: and signs of the 3456 term do not depend strongly on the magnitude of the n-body BE contributions themselves. Therefore, in the absence of other information, the 7789: 3456 term should always be included in water cluster BE calculations. The Gaussian-09 program was used to compute the energy contributions associated with the ACP terms.49 Gaussian-09 was also used for subsequent validation calculations.

Methods and Basis Sets for which ACPs were developed. In principle, ACPs need to be developed for every atom in the system (in our case, oxygen and hydrogen) and for a particular functional and basis set combination. In previous work on ACPs designed to mitigate the effects of basis set incompleteness, we demonstrated that fairly good ACP transferability between density functionals can be achieved, although with some loss in the accuracy of predicted properties.43 However, this transferability arises from the nature of the basis-set incompleteness error and we do not expect that water-specific ACPs would show such a high degree of transferability because of the significant variability in the performance of different functionals for water (see Fig. 2) . In order to determine the best functional for our ACP, we calculated the BEs for the eight water hexamer conformations reported in reference 1 using a variety of functionals and the aug-cc-pVQZ basis sets, with and without dispersion corrections. The plot in Figure 2 summarizes the results of our survey, which shows the signed error in the BE of the prism conformer of the water hexamer against the mean absolute error (MAE) predicted for eight water hexamer conformers. The data in Fig. 2 supports the conclusions made in previous works that DFT-based methods generally perform poorly for water40. The figure also shows that BHandHLYP, which to the best of our knowledge has never been tested for water clusters, is the best performing functional among the ones we explored. We then selected a series of basis sets to use in conjunction with BHandHLYP in order to span the range from nearly complete to minimal basis sets. Our list includes aug-cc-pVTZ,50,51 pc-2-spd52 (the pc-2 basis set of Jensen53,54 without the oxygen ffunctions), 6-31+G(d,p),55 6-31G(d),55 and STO-3G.55 In addition, we considered three 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 31

different dispersion correction methods: no dispersion (i.e. the ACP accounts for dispersion, which is not a leading contribution to binding in water), Grimme’s D3(BJ) method56 (henceforth referred to as “D3”), and XDM.57,58,59 All combinations of basis set and dispersion model were considered. Our preliminary findings, summarized in the SI, indicate that the performance of BHandHLYP/basis-ACP degenerates fairly rapidly as the size of the basis sets decreases. Therefore, we focus our attention on the results obtained using the aug-cc-pVTZ basis sets, which are nearly complete for water.

Figure 2. Binding energy error for the prism hexamer (negative = overbinding) using several functional and dispersion correction combinations against mean absolute errors in the calculation of relative BEs for the eight hexamers in Reference 1 relative to the prism BE. All calculations use aug-cc-pVQZ.

Validation of ACPs. The quality of the proposed BHandhLYP/aug-cc-pVTZ ACPs was tested on the fitting set of 1069 data points. As noted above, this was done to ensure that the LS fitting process does not deviate significantly from the results obtained by using the resulting ACPs in a self-consistent calculation. In addition, the performance of our ACPs was assessed on the 17 water clusters comprising 7 to 10 monomers from reference 1. To the best of our knowledge, reference 1 contains the highest-level ab initio wavefunction BEs for water clusters in this size range. For clusters larger than 10 water 10 ACS Paragon Plus Environment

Page 11 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecules, high-level BE is scarce. As mentioned in the introduction, Xantheas et al. reported BEs for (H2O)16 and (H2O)17 using MP2 and CCSD(T) with the aug-cc-pVTZ basis sets.4 It is well known that correlation energies, and properties derived from them like binding energies, converge very slowly to the complete basis set limit.16,60 Therefore, the quality of the BE data presented in reference 4 is questionable. More recently, Xantheas et al. revised the BEs of (H2O)16 and (H2O)17 on the basis of the observation that averaging the counterpoise-61 and non-counterpoise-corrected energies obtained with incomplete basis sets provides BEs in much better agreement with the CBS BEs obtained for smaller clusters. Xantheas’s findings are consistent with our own in this connection,60 and those of others.62

Results and Discussion Fitting Results for the Generated ACPs Of all the basis sets studied in this work, aug-cc-pVTZ performs best, so our discussion will focus on BHandHLYP/aug-cc-pVTZ, with the remainder of the results largely relegated to the SI. The water-specific ACPs optimized for use with BHandHLYP/aug-cc-pVTZ (no dispersion correction) are presented in Table 1. A copy of this and the remaining ACPs in a format suitable for use with the Gaussian program, as well as a sample input file, are provided in the SI. The small magnitude of the ACP coefficients are consistent with the fact that BHandHLYP/aug-cc-pVTZ performs reasonably well on its own (see also Table 2). When the D3 or XDM dispersion corrections are applied, the ACP coefficients determined by the LS fitting tend to be smaller (see the SI). This is expected because dispersion corrections lead to somewhat improved descriptions of the bonding between water molecules, so a smaller energy correction from the ACPs is required. As the size of the basis sets are reduced, the energy corrections from the ACPs must account for both deficiencies in the ability of the underlying functional to describe water interactions and for basis set incompleteness error, so ACPs tend to have parameters that provide larger energy corrections (see the SI). The performance of BHandHLYP/aug-cc-pVTZ with and without ACPs and dispersion corrections is illustrated in Table 2, with additional results presented in the SI for smaller basis sets. In the absence of any corrections, BHandHLYP predicts the errors for the two-, three- and four-body interactions contributions to BEs with mean signed 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 31

errors (MSEs) of 0.0735, -0.0179, and 0.012 kcal/mol, respectively (positive is underbinding). These mean errors are deceptively small, and have the potential to contribute to large BE errors with increasing cluster size because of the factorial dependence on the number of the 2-, 3-, and 4-body terms that are generated as a function of n in (H2O)n. However, there are two factors that constrain the magnitude of BE errors. Firstly, the 3-body MSE has a sign that is opposite to those of the 2- and 4-body MSEs and this will result in some error cancellation. Secondly, the contribution of the errors in the 2-4 body interaction energy terms to the total BE decrease quickly when the distances separating two (or more) of the monomers increases. In the case of bare BHandHLYP/aug-cc-pVTZ, the BE errors in (H2O)n for n=3-6 are in the range 1.0-2.1 kcal/mol, which is lower than many other DFT-based methods.40 Table 1. Water-specific hydrogen and oxygen atom-centered potential (ACP) optimized for use with BHaHLYP/aug-cc-pVTZ without dispersion corrections.a Atom

Function Type

O

local

H

ζi

ci 0.08

-0.000 703 194

0.12

0.000 569 923

0.0075

-0.000 074 833

0.02

0.000 096 138

s

0.004

-0.000 250 272

p

0.8

0.039 662 879

0.12

-0.000 112 098

0.18

-0.000 241 892

0.22

0.001 433 737

0.004

0.000 019 100

s

0.004

-0.000 184 088

p

0.004

-0.000 206 693

local

a

The SI contains a copy of this and the rest the ACPs for all basis set/dispersion combinations in Gaussian format.

Inclusion of D3 or XDM dispersion corrections in BHandHLYP/aug-cc-pVTZ increases the MSE associated with the two-body terms and changes its sign, but affects the three- and four-body errors differently. D3 corrections leaves the three- and four12 ACS Paragon Plus Environment

Page 13 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


body errors effectively unchanged while the XDM correction reduces the three-body error while making the four-body error worse by about the same factor. Both D3 and XDM are pairwise dispersion corrections, so the changes to the three- and four-body terms arise from the dependence of the dispersion coefficients on the electronic environment, which seems to be much more pronounced for XDM than for D3. To some extent, the improvement or deterioration in the errors of the predicted BEs for (H2O)n=2reflect the errors in the 3- and 4-body terms. For the hexamers, for example, BHandHLYP/aug-cc-pVTZ with dispersion correction performs worse than the bare 6

functional.

Table 2. Mean signed error (MSE) in BHandHLYP/aug-cc-pVTZ predicted binding energies for the members of the fitting set (kcal/mol) relative to CCSD(T)/CBSa and improvement factorsb in the MSEs with the application of D3 and XDM dispersion corrections, with and without ACPs. Fitting (n)

Set BHandHLYP MSE

2-body (410) 3-body (580) 4-body (58) Dimer (1) Trimer (2) Tetramer (3) Pentamer (7) Hexamer (8)

Improvement Factor D3

XDM

ACP

D3/ACP

XDM/ACP

0.07352

-0.72

-0.78

3.36

4.66

5.75

-0.01791

1.00

1.19

1.53

1.67

2.15

0.01242

1.00

0.92

1.18

0.85

0.82

0.36

-3.44

-5.08

76.79

25.88

76.68

1.01

-2.03

-2.94

9.70

19.70

20.03

1.29

-1.00

-1.23

196.16

34.08

43.95

2.11

-1.26

-1.61

-935.82

149.30

177.75

1.88

-0.72

-0.83

-162.46

-84.62

-101.36

a

See text for method details. Defined as: MSEBHandHLYP/aug − cc − pVTZ TMSEmethod

b

The ACPs without dispersion corrections, improve the BEs of the two-, three- and four-body terms by 1.1-3.4 times. Only BHandHLYP/aug-cc-pVTZ-ACP improves all of these terms. All of the ACP-based methods reduce the BE errors (H2O)n=4-6 by double- or triple-digit factors. This is a demonstration of how ACPs improve BHandHLYP/aug-ccpVTZ, with and without dispersion corrections, in its ability to predict water cluster BEs. However, it is important to underscore the fact that all of the datasets listed in Table 2 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 31

were used to fit the ACPs. The results for smaller basis sets are not as promising as those obtained with aug-cc-pVTZ. Although ACPs generally improve the predicted BEs in all fitting datasets, the errors predicted by the underlying BHandHLYP/basis set combination when the basis set is significantly smaller than aug-cc-pVTZ can become too large for the ACPs to correct reliably (see the SI). It may be possible to find a DFT-based method and basis set combination for which the errors in 2-, 3- and 4-body terms neatly cancel. However, this kind of approach is, in general, not reliable since the contribution from each of these terms depends on the size and hydrogen bonding network of the particular water cluster under study. To ensure that the application of ACPs do not cause unusual changes in the structure of water, we applied them to the optimization of the water monomer. The results of the optimizations are presented in Table 3.

Table 3. Optimized geometries (R, Å; θ,⁰) and dipole moment (μ,Debye) for the water monomer using BHandHLYP/aug-cc-pVTZ with and without dispersion corrections and ACPs. Method

R

Θ

Expt.

0.9572

a

a

104.52

MP2/aug-cc-pV5Z BHandHLYPa BHandHLYP-D3 BHandHLYP-XDM BHandHLYP/ACP BHandHLYP-D3/ACP BHandHLYP-XDM/ACP

0.9582 0.9504 0.9504 0.9503 0.9612 0.9610 0.9612

104.40 105.77 105.77 105.78 104.72 105.17 105.09

μ 1.8546b 1.8629 1.8917 1.8917 1.8915 1.8572 1.8546 1.8525

a

Reference 63. bReference 64.

The results compiled in Table 3 demonstrate that the application of ACPs does not result in unusual geometric distortions or a significant change in dipole moment in the water monomer. This is to be expected based on the small magnitude of the energy corrections generated by the ACPs. In fact all quantities in Table 3 (O-H bond distance, angle, and dipole moment) are brought into better agreement with the experimental value upon application of the ACPs despite the fact that these quantities were not 14 ACS Paragon Plus Environment

Page 15 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


targeted in the fit. The BHandHLYP-ACP dipole moments are in even better agreement with the experimental value than that obtained from MP2/aug-cc-pV5Z. These findings lead to the expectation that ACPs can be used to compute water cluster structures in fairly good agreement with those obtained by nearly complete basis set wavefunction methods. However, more stringent tests, like geometry optimizations on a large water clusters, are required to make definitive conclusions in this regard. Performance of ACPs for (H2O)n, n=7-10 The performance of the ACPs is assessed by comparison to the CCSD(T)/CBSquality reference BEs from Temelso et al.1 for the (H2O)n, n=7-10 clusters. Table 4 compares BHandHLYP/aug-cc-pVTZ with and without dispersion corrections and ACPs to the reference data in terms of signed error. Uncorrected BHandHLYP/aug-cc-pVTZ underestimates the BEs by an amount that is roughly linear in the number of water monomers. For the 11 (H2O)7 clusters, the underbinding error spans a range of 1.7 – 3.7 kcal/mol, which is significant compared to the total BEs (around 55 kcal/mol). This suggests that the BE errors predicted by BHandHLYP/aug-cc-pVTZ are large enough that even the reliable prediction of the energy ordering of the heptamers is not possible. Indeed, seven of the 11 (H2O)7 structures are incorrectly ranked. The inclusion of D3 dispersion uniformly over-corrects the BEs, resulting in approximately 3 kcal/mol mean overbinding errors for the (H2O)n, n=7-10 clusters. The use of XDM has a similar effect on the BEs, with over-binding errors of about 3.4 kcal/mol on the same set. However, both XDM and D3 predict the correct ranking of the 11 (H2O)7 structures. The use of ACPs without dispersion corrections improves the calculated BEs; errors are reduced by a factor greater than 40. All (H2O)7 structures are correctly ranked with ACPs. However, BHandHLYP/aug-cc-pVTZ-ACP is unable to properly rank the (H2O)8 and (H2O)10 clusters, which are degenerate to within 0.01 kcal/mol. Given that the reference energies have some degree of uncertainty that is likely to be larger than 0.01 kcal/mol, the performance of BHandHLYP/aug-cc-pVTZ-ACP can be considered to be excellent: Over all, the method predicts mean absolute errors for clusters containing 710 molecules that average 0.009 kcal/mol per water monomer. Similar performance is seen for ACPs used with either of the dispersion corrections. Errors are slightly higher, and span a slightly larger range for D3-ACP than XDM-ACP. However, both methods correctly rank the (H2O)7 clusters. Nevertheless, mean absolute errors obtained with D3-ACP and XDM-ACP average 0.015 and 0.019 kcal/mol per water molecule. To our knowledge, this level of performance is better than 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 31

any DFT-based approach known to date.

Table 4. Reference binding energies (BE)a and signed errorsb using BHandHLYP/aug-ccpVTZ for (H2O)n, n=7-10. Also shown are the signed errors obtained with the inclusion of D3 or XDM dispersion corrections and/or ACPs. All values are in kcal/mol. Signed Error (H2O)nc

Reference

BHandHLYP

D3

XDM

ACP

D3/ACP

XDM/ACP

n=7 7-BI2

-53.33

-2.28

3.02

2.54

0.06

0.03

0.07

7-BI1

-53.40

-2.23

3.05

2.57

0.07

0.02

0.05

7-CA1

-55.68

-2.98

3.40

2.92

0.00

0.09

0.13

7-CA2

-54.88

-3.09

3.25

2.76

0.00

0.08

0.13

7-CH1

-54.38

-1.83

3.34

2.92

0.09

-0.02

0.03

7-CH2

-53.86

-2.05

3.22

2.81

0.03

-0.07

-0.02

7-CH3

-53.07

-1.97

3.22

2.80

0.06

0.07

0.12

7-HM1

-52.35

-1.73

2.99

2.59

0.07

-0.02

0.03

7-PR1

-57.37

-3.56

3.30

2.72

0.04

0.12

0.15

7-PR2

-57.10

-3.70

3.21

2.63

0.00

0.05

0.09

7-PR3

-56.90

-3.48

3.32

2.76

0.06

0.12

0.16

n=8 8-D2d

-72.54

-4.34

4.48

3.74

0.03

0.16

0.19

8-S4

-72.56

-4.36

4.47

3.73

0.01

0.13

0.17

n=9 9-D2dDD

-81.71

-4.31

5.29

4.54

0.14

0.27

0.30

9-S4DA

-81.46

-4.50

5.12

4.36

0.04

0.13

0.16

n=10 10-PP1

-92.89

-4.91

6.14

5.27

0.35

0.44

0.47

10-S4DA

-92.89

-4.93

6.10

5.23

0.30

0.40

0.42

-3.31

3.94

3.35

0.08

0.12

0.16

3.31

3.94

3.35

0.08

0.13

0.16

4.93

6.14

5.27

0.35

0.44

0.47

Mean Signed Error Mean Absolute Error Maximum Absolute Error a

b

Taken from reference 1. Signed error = Reference – DFT-based value. Negative means under-bound. cThe nomenclature used for specific water clusters is taken from 16 ACS Paragon Plus Environment

Page 17 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


reference 1.

Figure 3 reiterates some of the elements of the foregoing discussion on the performance of BHandHLYP/aug-cc-pVTZ with ACPs. The plot shows the mean absolute errors (MAE) as a function of cluster size for (H2O)n cluster with n=2-10, relative to the reference BEs.1 The errors obtained with non-ACP-corrected methods tend to increase roughly linearly as a function of cluster size, regardless of the choice of dispersion correction. The ACP-corrected methods present errors that are, in comparison, negligible. Although we expect the error in ACP-corrected BHandHLYP to scale linearly with the number of water molecules for sufficiently large clusters, the clusters in Figure 3 have errors so small that this is not apparent.

Figure 3. Mean absolute error in BHandHLYP/aug-cc-pVTZ binding energies for each of the (H2O)n clusters with n=2-10 clusters in Temelso et al. (reference 1), with and without dispersion corrections and/or ACPs. The inset shows a close-up of the ACP results.

Figure 4 provides similar information about BHandHLYP-D3 combined with a 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 31

variety of basis sets with varying degrees of incompleteness, and with and without the ACPs optimized for use with those basis sets. The plot demonstrates that the error in the calculated BEs increases rapidly with decreasing basis set size, and that for non-ACPcorrected methods, the BE errors are linear in the number of water molecules. ACPs are able to correct in all cases, with varying degrees of success. In the most extreme case, (H2O)10 using BHandHLYP-D3/6-31G(d), the error is reduced from ca. 54 over-binding to about 0.43 (under-binding) kcal/mol, resulting is a 126-fold improvement. However, all of the methods, including those that make use of ACPs, predict four of the 11 (H2O)7 clusters with incorrect ranking. The inability to discern the energy ranking of structures puts the general reliability of BHandHLYP/6-31G(d)-ACP into question.

Figure 4. Mean absolute error in BHandHLYP binding energies for each of the (H2O)n clusters with n=2-10 clusters presented in reference 1, with and without dispersion corrections and/or ACPs. The inset shows a close-up plot of the results obtained with the ACPs.

For the methods that use more complete basis sets, like BHandHLYP-D3/pc-2spd, Figure 4 shows that results resemble more closely those obtained using aug-ccpVTZ basis sets. However, no method predicts the energy ordering of the 11 (H2O)7 clusters correctly. The ACPs alone and with the D3 dispersion correction predict the 18 ACS Paragon Plus Environment

Page 19 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


order of two of the structures incorrectly, and the XDM-ACP combination predict the ordering of three structures incorrectly. BHandHLYP-D3/pc-2-spd with no ACPs performs worst, getting none of the structures with correct energy ordering. These results suggest that nearly complete basis sets are required to describe the subtle differences in the various interactions between the molecules in water clusters, and that the energy differences that a method is able to discern depends directly on the size of the basis sets.

Application of BHandHLYP/aug-cc-pVTZ-ACP to (H2O)n, n=16,17 In 2015, Miliordos and Xantheas described a computational protocol that, according to the authors, offers a relatively efficient ab initio approach to calculate the BEs in larger water clusters.5 They demonstrated that the MP2/CBS BEs for (H2O)n n=26,8 can be obtained to within 99.7 ± 0.7 % by averaging the counterpoise and noncounterpoise corrected BEs calculated using MP2/aug-cc-pVTZ. The same result was obtained from averaging the counterpoise and noncounterpoise corrected CCSD(T)/CBS BEs for (H2O)n n=2-4. Utilizing this averaging approach, Miliordos and Xantheas estimated the BEs for four structures of (H2O)16 and two structures of (H2O)17 clusters. Although some uncertainties associated with the assumption that the good performance of this averaging approach holds for increasing cluster sizes remain, the BEs obtained by Miliordos and Xantheas are 6-7 kcal/mol lower than those obtained using CCSD(T)/aug-cc-pVTZ, previously calculated by Yoo et al.4 These findings underscore the slow convergence of ab initio wavefunction BEs to the basis set limit, and illustrate the difficulty in reliably calculating high-quality BEs for moderately large water clusters. We applied the BHandHLYP/aug-cc-pVTZ-ACP approach to the (H2O)16 and (H2O)17 clusters described in Miliordos and Xantheas, using the geometries provided in the Supporting Information of reference 5. The results are summarized in Table 5. The BHandHLYP/aug-cc-pVTZ BE results continue the trend observed for smaller water clusters (see Table 4). Without ACPs or corrections for dispersion, the method predicts BEs that are under-bound by ca. 10-13.5 kcal/mol for the (H2O)16 and (H2O)17 clusters. On the other hand, and again in keeping with the results for smaller water clusters, D3and XDM-corrected BE for (H2O)16 and (H2O)17 are over-bound by 9.8-10.9 and 7.7-9.0 kcal/mol, respectively. These errors in the dispersion-corrected BEs are expected based the approximate linear scaling of BE error (see Figure 4). 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 31

The use of ACPs with BHandHLYP/aug-cc-pVTZ, with and without dispersion corrections, uniformly improves the predicted BEs, bringing most of the values to within the errors bars associated with the latest reference BEs. Given the uncertainties associated with the manner with which the reference BEs were obtained, the subkJ/mol errors associated with the ACP results for the (H2O)7-10 clusters, and the agreement between different ACP-corrected values, our DFT-based results may predict the BEs of the (H2O)16 and (H2O)17 clusters more accurately than any previous study using wavefunction methods.

Table 5. Reference binding energies (BE)a and signed errorsb in BHandHLYP/aug-cc-pVTZ BEs for (H2O)16 and (H2O)17. Also shown are the signed errors obtained with BHandHLYP with the inclusion of D3 or XDM dispersion corrections and/or ACPs. All values are in kcal/mol. Signed Errors (H2O)nb

Referencec

BHandHLYP

D3

XDM

ACP

D3/ACP

XDM/ACP

n=16 4444-a

-164.2 ± 1.1

-12.1

10.2

8.2

-1.6

-0.3

-0.6

Boat-a

-164.4 ± 1.6

-10.0

10.3

8.5

-1.5

-0.8

-1.0

Boat-b

-164.2 ± 1.6

-10.0

10.3

8.5

-1.4

-0.7

-0.9

Anti-boat

-164.6 ± 1.6

-10.5

9.8

8.1

-1.9

-1.2

-1.4

4444-b

-164.1 ± 1.6

-12.7

9.8

7.7

-2.2

-0.9

-1.2

n=17 Sphere

-175.7 ± 1.8

-13.5

10.5

8.4

-2.0

-1.5

-1.8

552ˊ5

-175.0 ± 1.8

-11.1

10.9

9.0

-1.4

-0.9

-1.0

a

Average of counterpoise and non-counterpoise corrected CCSD(T)/aug-cc-pVTZ binding energies from reference 5. bSigned error = Reference – DFT-based value. Negative means under-bound. cNomenclature as in reference 5.

The Minimum-Energy Structure of (H2O)25. As a final application of the BHandHLYP/aug-cc-pVTZ ACPs, we seek to resolve which of the structures for (H2O)25 presented in the literature represent the lowest energy structure. Several recent works have explored global optimization techniques coupled with force field6,7 and ab initio10 methods. Takeuchi predicted two structures for (H2O)25 based on the TIP3P and TIP4P force fields,65 while more recently Kazachenko 20 ACS Paragon Plus Environment

Page 21 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and Thakkar used TIP4P, TTM2.1-F,66 TIP4P-Ew,67 TIP4P/2005,68 and AMOEBA69,70 in conjunction with their search algorithm to predict the binding energies of global minima with each force field.7 In addition, Xantheas used MP2 with aug-cc-pVTZ and aug-ccpVDZ to extrapolate to the complete basis set limit to report the binding energy for seven low-energy structures for (H2O)25. Table 6. Binding energies of proposed global minimum structures for (H2O)25 predicted using various methods (Lit. BE), with those predicted by BHandHLYP-Disp/aug-cc-pVTZACP (“Disp” = D3, XDM). The lowest energy DFT-based clusters are indicated in bold font. All values are in kcal/mol. Method

Lit. BE

ACP

D3/ACP

XDM/ACP

-251.9

-250.8

-255.9

-254.6

Reference 6 TIP3P

a

TIP4P

a

-247.4 -251.9 Reference 7

AMOEBA

-268.23

-266.6

-269.4

-268.9

TIP4P/2005

-293.50

-254.2

-258.0

-256.9

TIP4P-Ew

-289.54

-252.6

-256.7

-255.5

TIP4P

-266.02

-251.9

-255.9

-254.6

TTM2.1F Structure

-266.2

-265.6

I

-277.48 -262.8 Reference 10: MP2 −281.38 -267.4

-269.1

-268.9

II

−281.74

-267.2

-268.8

-268.6

III

−281.03

-267.0

-268.7

-268.5

IV

−281.15

-267.1

-268.9

-268.7

V

−280.60

-267.0

-268.8

-268.6

VI

−280.57

-266.7

-268.4

-268.2

VII

−280.43

-266.4

-268.1

-267.9

a

Not provided.

We subjected each of the 15 distinct (H2O)25 structures previously reported to BHandHLYP/aug-cc-pVTZ-ACP calculations in order to rank them by BE and identify the global minimum. The results are summarized in Table 6. Using ACPs without dispersion correction, the BEs for each of the minimum structures proposed span a range of 20 kcal/mol. Another interesting feature of the reference BEs listed in the Table for the 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 31

MP2 derived BEs for the lowest seven structures of (H2O)25 is that they are ca. 15 kcal/mol higher in energy than those predicted using BHandHLYP. The latter point reinforces what is known about the ability to accurately extrapolate to the complete basis set limit using aug-cc-pVDZ energies. For example, using the extrapolation scheme employed in reference 10 with the RI-MP2/aug-cc-pVTZ and aug-cc-pVDZ BEs reported in reference 1 for (H2O)8 gives a BE that is lower (i.e. more binding) than the true CBS limit by about 0.36 kcal/mol per water molecule. Assuming that this error scales linearly with cluster size results in over-binding in (H2O)25 by 9 kcal/mol. Furthermore, the 7789: absence of higher-order correlation in the BEs (eg. through 3456 corrections, see Figure 1) could introduce even more over-binding error. The results of our survey of the proposed minimum-energy structures of (H2O)25 reveal that the structure proposed by Kazachenko and Thakkar using the AMOEBA force field has the lowest BE, if dispersion corrections are applied with the ACPs. Interestingly, the D3-ACP and XDM-ACP approaches predict the BEs for (H2O)25 AMOEBA and (H2O)25 I to be within 0.3 and 0.0 kcal/mol of each other, respectively. Note that the root mean square geometries of these two structures is 5 Å. With ACPs alone, the (H2O)25 I structure is lower than the (H2O)25 AMOEBA structure by 0.8 kcal/mol. We find that the BHandHLYP/aug-cc-pVTZ-ACP result is within 0.5 kcal/mol of the BE predicted using the AMOEBA force field, and that the inclusion of dispersion corrections results in slightly stronger binding. It is interesting to note that the non-ACP results also predict the AMOEBA structure to be the lowest in energy (results not shown), although there are much larger differences between the predicted BEs and those obtained with the force field. All methods also predict the same structures for the 2nd-(TTM2.1F) and 3rd-ranked (I) clusters. In general, the energy separation between energetically low-lying clusters of large molecules are expected to decrease with cluster size. As the clusters become large, it should be increasingly difficult to resolve the lowest energy structure by any method. The (H2O)25 results serve as a good example of this: The B3LYP-D3/aug-cc-pVTZ-ACP energy separation between (H2O)25 AMOEBA and (H2O)25 I of 0.3 kcal/mol represents just over 0.1% of the calculated binding energy. This energy difference is also about the 7789: magnitude of the 3456 term in water dimers and is smaller than the largest mean absolute error we predict using ACPs for the binding energy of (H2O)10. Therefore, Identification of the lowest-energy water structures for clusters beyond this size may be futile even for complete basis set ab initio methods. Conclusions 22 ACS Paragon Plus Environment

Page 23 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


This work explored the use of water-specific atom-centered potentials (ACPs) to correct the underlying deficiencies in the BHandHLYP density functional coupled with a variety of basis sets and two dispersion-correction schemes in order to accurately predict the binding energies (BEs) of water clusters. The ACPs were fit to 1069 highquality 2-, 3-, and 4-body interaction energies for dimers, trimers, and tetramers extracted from the clusters in the work of Temelso et al.1 and the four (H2O)20 clusters in Bryantsev et al.3 Also included in the fitting set were the binding energies for (H2O)n with n=2-6 from reference 1. It was demonstrated that ACPs used with BHandHLYP/aug-cc-pVTZ, with or without dispersion corrections, are capable for predicting the binding energies of water clusters containing up to 10 molecules to within 0.44 kcal/mol. ACPs were also developed for use with smaller basis sets, and these were also seen to perform well, although the errors in the binding energies are correspondingly higher. Our newly developed ACPs were applied to water clusters containing 16 and 17 molecules proposed by Yoo et al.4, which have been recently revised by Miliordos and Xantheas5, who proposed a downward correction to the binding energies of around 6 kcal/mol. We found that ACP-based approaches predict binding energies for (H2O)16,17 to within 2.2 kcal/mol of the revised, high-level ab initio BEs. Finally, we used BHandHLYP/aug-cc-pVTZ-ACP to survey a number of proposed minimum-energy structures for the (H2O)25 cluster, which were obtained using a variety of methods (force field and MP2) and conformer search strategies. Dispersion-corrected BHandHLYP/aug-cc-pVTZ-ACP methods predict one particular structure generated by Kazachenko and Thakkar using the AMOEBA force field to be the lowest-energy structure of those available in the literature, with the “I” structure being nearly isoenergetic. It is expected that the resolution of the lowest-energy structure for water clusters will become increasingly difficult with increasing cluster size, and that that ability to do so may not be possible using any currently available computational technique. The present work provides strong evidence that system-specific ACPs can be used to answer key questions in systems of broad interest and importance, such as water. The development of BHandHLYP/aug-cc-pVTZ-ACP, and its corresponding dispersion-corrected versions, provide a new tool for exploring the structures of water clusters ranging in size from two to several hundred molecules with unprecedented (ca. sub 0.1 kcal/mol per water molecule) accuracy. Our future efforts in this area will focus on refinements to the ACP approach that will allow for the use of smaller basis sets and 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 31

less computationally intensive DFT-based methods without accuracy loss.

Acknowledgements We thank Professors George Shields, Ajit Thakkar, and Gregory Tschumper for helpful discussions during the early stages of this work. GAD is grateful to the University of British Columbia and the National Sciences and Engineering Research Council (NSERC) for their financial support and to Compute Canada for a generous allocation of computing resources. JDH thanks NSERC for an undergraduate research award. Supporting Information Available Plots of the 2-, 3-, and 4-body interaction terms as a function of distance, fitting data used to generate ACPs, the results of ACPs applied to the systems described in this work in conjunction with BHandHLYP-Disp/basis (Disp = none, D3(BJ), XDM). This material is available free of charge from www.acs.org.

References

1. Temelso, B.; Archer, K. A.; Shields, G. C. Benchmark Structures and Binding Energies of Small Water Clusters with Anharmonicity Corrections. J. Phys. Chem. A 2011, 115, 12034-12046. 2. Bates, D. M.; Tschumper, G. S. CCSD(T) Complete Basis Set Limit Relative Energies for Low-Lying Water Hexamer Structures J. Phys. Chem. A 2009, 113, 3555-3559. 3. Bryantsev, V. S.; Diallo, M. S.; van Duin, A. C. T.; Goddard III, W. A. Evaluation of B3LYP, X3LYP, and M06-Class Density Functionals for Predicting the Binding Energies of Neutral, Protonated, and Deprotonated Water Clusters. J. Chem. Theory Comput. 2009, 5, 1016-1026. 4. Yoo, A; Aprà, E.; Zeng, X. C.; Xantheas, S. S. High-Level An Initio Electronic Structure Calculations of Water Clusters (H2O)16 and (H2O)17: A New Global Minimum for (H2O)16. J. Phys. Chem. Lett. 2010, 1, 3122-3127. 5. Miliordos, E.; Xantheas, S. S. An Accurate and Efficient Computational Protocol for Obtaining the Complete Basis Set Limits of the Binding Energies of Water Clusters at the MP2 and CCSD(T) levels of Theory: Application to (H2O)m, m=2-6, 8, 11, 16, and 17. J. 24 ACS Paragon Plus Environment

Page 25 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Chem. Phys. 2015, 142, 234303. 6. Takeuchi, H. Development of an Efficient Geometry Optimization Method for Water Clusters. J. Chem. Inf. Model. 2008, 48, 2226-2233. 7. Kazachenko, S.; Thakkar, A. J. Water Nanodroplets: Predictions of Five Model Potentials. J. Chem. Phys. 2013, 138, 194302. 8. Gillan, M. J.; Manby, F. R.; Towler, M. D.; Alfè. D. Assessing the accuracy of quantum Monte Carlo and density functional theory for energetics of small water clusters. J. Chem. Phys. 2012, 136, 244105. 9. Alfè, D.; Bartók, A. P.; Csányi, G. Gillan, M. J. Energy benchmarking with quantum Monte Carlo for water nano-droplets and bulk liquid water. J. Chem. Phys. 2013, 138, 221102. 10. Sahu, N.; Gadre, S. R.; Rakshit, A.; Bandyopadhyay, P.; Miliordos, E. Xantheas, S. S. Low Energy Isomers of (H2O)25 from a Hierarchical Method Based on Monte Carlo Temperature Basin Paving and Molecular Tailoring Approaches Benchmarked by MP2 Calculations. J. Chem. Phys. 2014, 141, 164304. 11. Howard, J. C.; Tschumper, G. S. Wavefunction Methods for the Accurate Characterization of Water Clusters. WIREs Comput. Mol. Sci. 2014, 4, 199–224. 12. Olson, R. M.; Bentz, J. L.; Kendall, R. A.; Schmidt, M. W.; Gordon, M. S. A Novel Approach to Parallel Coupled Cluster Calculations: Combining Distributed and Shared Memory Techniques for Modern Cluster Based Systems. J. Chem. Theory Comput. 2007, 3, 1312–1328. 13. Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory. J. Chem. Phys. 2016, 144, 024109. 14. Dahlke, E. E.; Olson, R. M.; Leverentz, H. R.; Truhlar, D. G. Assessment of the Accuracy of Density Functionals for Prediction of Relative Energies and Geometries of Low-Lying Isomers of Water Hexamers. J. Phys. Chem. A 2008, 112, 3976–3984. 15. Santra, B.; Michaelides, A.; Fuchs, M.; Tkatchenko, A.; Filippi, C.; Scheffler, M. On the accuracy of density-functional theory exchange-correlation functionals for H bonds in small water clusters. II. The water hexamer and van der Waals interactions. J. Chem. Phys. 2008, 129, 194111. 16. DiLabio, G. A.; Otero-de-la-Roza, A.; Noncovlaent Interactions in Density-Functional Theory. In Reviews in Computational Chemistry, Parrill, A. L.; Lipkowitz, K. B., Eds. John Wiley and Sons, Hoboken, NJ, 2016; Vol. 29, 1-97. 25 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 31

17. Hermann, J.; DiStasio Jr., R. A; Tkatchenko, A. First-Principles Models for van der Waals Interactions in Molecules and Materials: Concepts, Theory, and Applications. Chem. Rev. 2017, 117, 4714–4758. 18. Non-covalent Interactions in Quantum Chemistry and Physics: Theory and Applications. Otero-de-la-Roza, A.; DiLabio, G. A., Eds.; Elsevier: Cambridge, MA, 2017. 19. Wu, X.; Vargas, M. C.; Nayak, S.; Lotrich, V.; Scoles, G. Towards extending the applicability of density functional theory to weakly bound systems. J. Chem. Phys. 2001, 115, 8748-8757. 20. Wu, Q.; Yang, W. Empirical correction to density functional theory for van der Waals interactions. J. Chem. Phys. 2002, 116, 515-524. 21. Grimme, S. Accurate description of van der Waals complexes by density functional theory including empirical corrections. J. Comput. Chem. 2004, 25, 1463-1473. 22. Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a LongRange Dispersion Correction. J. Comput. Chem. 2006, 27, 1787-1799. 23. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu . J. Chem. Phys. 2010, 132, 154104. 24. Johnson, E. R.; Becke, A. D. A post-Hartree–Fock model of intermolecular interactions. J. Chem. Phys. 2005, 123, 024101. 25. Becke, A. D.; Arabi, A. A.; Kannemann, F. O. Nonempirical density-functional theory for van der Waals interactions. Can. J. Chem. 2010, 88, 1057-1062. 26. Otero-de-la-Roza, A.; Johnson, E. R. A benchmark for non-covalent interactions in solids. J. Chem. Phys. 2012, 137, 054103. 27. Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215-241. 28. Zhao, Y.; Truhlar, D. G. Exploring the Limit of Accuracy of the Global Hybrid Meta Density Functional for Main-Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2008, 4, 1849–1868. 29. Peverati, R.; Truhlar, D.G. Improving the Accuracy of Hybrid Meta-GGA Density Functionals by Range Separation. J. Phys. Chem. Lett. 2011, 2, 2810-2817. 30. Torres, E.; DiLabio, G. A. A (Nearly) Universally Applicable Method for Modeling Noncovalent Interactions Using B3LYP. J. Phys. Chem. Lett. 2012, 3, 1738–1744. 31. DiLabio, G. A.; Koleini, M. Dispersion-Correcting Potentials Can Significantly Improve 26 ACS Paragon Plus Environment

Page 27 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the Bond Dissociation Enthalpies and Noncovalent Binding Energies Predicted by Density-Functional Theory. J. Chem. Phys. 2014, 140, 18A542. 32. DiLabio, G. A.; Atom-centered potentials for noncovalent interactions and other applications. In Non-covalent Interactions in Quantum Chemistry and Physics: Theory and Applications. Otero-de-la-Roza, A.; DiLabio, G. A., Eds.; Elsevier: Cambridge, MA, 2017; 221-240. 33. von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U. Optimization of Effective Atom Centered Potentials for London Dispersion Forces in Density Functional Theory. Phys. Rev. Lett. 2004, 93, 153004. 34. von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U. Performance of optimized atomcentered potentials for weakly bonded systems using density functional theory. Phys. Rev. B 2005, 71, 195119. 35. Lin, I.-C.; Coutinho-Neto, M. D.; Felsenheimer, C.; von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U. Library of dispersion-corrected atom-centered potentials for generalized gradient approximation functionals: Elements H, C, N, O, He, Ne, Ar, and Kr. Phys. Rev. B 2007, 75, 205131. 36. Lin, I.-C.; Seitsonen, A. P.; Tavernelli, I.; Rothlisberger, U. Structure and Dynamics of Liquid Water from ab Initio Molecular Dynamics - Comparison of BLYP, PBE, and revPBE Density Functionals with and without van der Waals Corrections. J. Chem. Theory Comput. 2012, 8, 3902−3910. 37. Karalti, O.; Su, X.; Al-Saidi, W. A.; Jordan, K. D. Correcting density functionals for dispersion interactions using pseudopotentials. Chem. Phys. Lett. 2014, 591, 133-136. 38. Otero-de-la-Roza. A.; Johnson, E. R. Non-covalent interactions and thermochemistry using XDM-corrected hybrid and range-separated hybrid density functionals. J. Chem. Phys. 2013, 138, 204109. 39. Leverentz, H. R.; Qi,H. W.; Truhlar, D. G. Assessing the Accuracy of Density Functional and Semiempirical Wave Function Methods for Water Nanoparticles: Comparing Binding and Relative Energies of (H2O)16 and (H2O)17 to CCSD(T) Results. J. Chem. Theory Comput. 2013, 9, 995−1006. 40. Gillan, M. J.; Alfè, D.; Michaelides, A. Perspective: How Good is DFT for Water? J. Chem. Phys. 2016, 114, 130901. 41. DiLabio, G. A. Accurate Treatment of van der Waals Interactions Using Standard Density Functional Theory Methods with Effective Core-Type Potentials: Application to Carbon-Containing Dimers. Chem. Phys. Lett. 2008, 455, 348-353. 42. Van Santen, J. A.; DiLabio, G. A. Dispersion Corrections Improve the Accuracy of 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 31

Both Noncovalent and Covalent Interactions Energies Predicted by a Density-Functional Theory Approximation. J. Phys. Chem. A 2015, 119, 6703-6713. 43. Otero-de-la-Roza, A.; DiLabio, G. A. Transferable Atom-Centered Potentials for the Correction of Basis Set Incompleteness Errors in Density-Functional Theory. J. Chem. Theory Comput. 2017, 13, pp 3505–3524 44. Torres, E.; DiLabio, G. A. Density-Functional Theory with Dispersion-Correcting Potentials for Methane: Bridging the Efficiency and Accuracy Gap Between High-Level Wavefunction and Classical Molecular Mechanics Methods. J. Comput. Chem. Theory 2013, 9, 3342-3349. 45. Christiansen, P. A.; Lee, Y. S.; Pitzer, K. S. Improved ab initio effective core potentials for molecular calculations. J. Chem. Phys. 1979, 71, 4445-4450. 46. Otero-de-la-Roza, A.; DiLabio, G. A.; Johnson, E. Exchange-Correlation Effects for Noncovalent Interactions in Density Functional Theory. J. Chem. Theory Comput. 2016, 7, 3160-3175. 47. Klopper, W. Limiting Values for Møller–Plesset Second-Order Correlation Energies of Polyatomic Systems: A Benchmark Study on Ne, HF, H2O, N2, and He…He. J. Chem. Phys. 1995, 102, 6168-6179. 48. TURBOMOLE V6.5 2013, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com. 49. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, M. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision D.01; Gaussian, Inc., Wallingford CT, 2009. 50. Dunning Jr, T. H. Gaussian Basis Sets for use in Correlated Molecular Calculations. I. The Atoms Boron Through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007-1023. 51. Kendall, R. A.; Dunning Jr, T. H.; Harrison, R. J. Electron Affinities of the First-Row 28 ACS Paragon Plus Environment

Page 29 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796-6806. 52. Johnson, E. R.; Otero-de-la-Roza, A.; Dale, S.G.; DiLabio, G. A. Efficient Basis Sets for Non-Covalent Interactions in XDM-Corrected Density-Functional Theory. J. Chem. Phys. 2013, 139, 214109. 53. Jensen, F. Polarization Consistent Basis Sets: Principles. J. Chem. Phys. 2001, 115, 9113-9125. 54. Jensen, F. Polarization Consistent Basis Sets. II. Estimating the Kohn-Sham Basis Set Limit. J. Chem. Phys. 2002, 116, 7372-7379. 55. As implemented in reference 49. 56. Grimme, S.; Antony, J. Ehrlich, S.; Krieg, H. A Consistent and Accurate ab initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. 57. Johnson, E. R.; Becke, A. D. A Post-Hartree-Fock Model of Intermolecular Interactions. J. Chem. Phys. 2005, 123, 024101. 58. Becke, A. D.; Arabi, A. A.; Kannenmann, F. O. Non-Empirical Density Functional Theory for van der Waals Interactions. Can. J. Chem. 2010, 88, 1057-1062. 59. Otero-de-la-Roza, A.; Johnson, E. R. A Benchmark for Non-Covalent Interactions in Solids. J. Chem Phys. 2012, 137, 054103. 60. Mackie, I. D.; DiLabio, G. A. Approximations to complete basis set-extrapolated, highly correlated non-covalent interaction energies. J. Chem. Phys. 2011, 135, 134318. 61. Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553-566. 62. Burns, L. A.; Marshall, M. S.; Sherrill. C. D. Comparing Counterpoise-Corrected, Uncorrected, and Averaged Binding Energies for Benchmarking Noncovalent Interactions. J. Chem. Theory Comput. 2014, 10 , 49–57. 63. Hoy, A. R.; Mills, I. M.; Strey, G. Anharmonic Force Constant Calculations. Mol. Phys. 1972, 24, 1265-1290. 64. Landolt-Börnstein, Numerical Data and Functional Relationships in Science and Technology - New Series, II/14a, Springer-Verlag, Heidelberg; 1982. 65. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926–935. 66. Partridge, H.; Schwenke, D. W. The Determination of an Accurate Isotope Dependent 29 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 31

Potential Energy Surface for Water from Extensive ab initio Calculations and Experimental Data. J. Chem. Phys. 1997, 106, 4618-4639. 67. Horn H.W.; Swope, W.C.; Pitera, J.W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; HeadGordon T. Development of an Improved Four-Site Water Model for Biomolecular Simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665-9678. 68. Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. 69. Ren, P.; Ponder, J. W. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. J. Phys. Chem. B 2003, 107, 5933-5947. 70. Ponder, J. W.; Wu, C. ; Ren, P.; Pande, V. S.; Chodera, J. D.; Mobley, D. L.; Schnieders, M. J.; Haque, I.; Lambrecht, D. S.; DiStasio, Jr. R. A.; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; Head-Gordon, T. Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114, 2549-2564.


Page 31 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57


ACS Paragon Plus Environment

Accurate Modeling of Water Clusters with Density-Functional Theory

Recommend Documents