Dispersion Interactions in Water Clusters - The Journal of Physical

Subscriber access provided by BALL STATE UNIV

Article

Dispersion Interactions in Water Clusters Emilie Brigitte Guidez, and Mark S. Gordon J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b11403 • Publication Date (Web): 05 Apr 2017 Downloaded from http://pubs.acs.org on April 6, 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.

The Journal of Physical Chemistry A 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 29

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 Journal of Physical Chemistry

Dispersion Interactions in Water Clusters Emilie B. Guidez and Mark S. Gordon* Department of chemistry, Iowa State University, Ames, IA 50011, United States *Email: [email protected] Phone: 515-294-0452

Abstract The importance of dispersion forces in water clusters is examined using the effective fragment potential (EFP) method. Since the original EFP1 water potential does not include dispersion, a dispersion correction to the EFP1 potential (EFP1-D) was derived and implemented. The addition of dispersion to the EFP1 potential yields improved geometries for water clusters that contain 2-6 molecules. The importance of the odd E7 contribution to the dispersion energy is investigated. The E7 dispersion term is repulsive for all of the water clusters studied here and can have a magnitude that is as large as half of the E6 value. The E7 term therefore contributes to larger intermolecular distances for the optimized geometries. Inclusion of many-body effects and/or higher order terms may be necessary to further improve dispersion energies and optimized geometries. Introduction The development of model potentials for the simulation of water remains an important research effort due to the omnipresence of water in the fields of biology, chemistry and engineering. Many water models have been developed in the past four decades.1 Molecular mechanics (MM) force fields are often parametrized to reproduce experimental bulk properties using classical Monte-Carlo or molecular dynamics simulations. Non-polarizable additive potentials represent an important class of force fields. In these methods, water molecules are 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

modeled with partial point charges. The most popular include TIP3P,2 TIP4P,2,3 TIP5P4,5, in which each rigid water molecule is assigned 3, 4 and 5 interaction sites respectively. Another example is the simple point charge model (SPC and its extended variant SPC/E) that does include polarization.6,7 Improvements on these force fields were developed by introducing internal molecular flexibility (in SPC/Fw8 and in TIP4P9 for example). The q-SPC/Fw10 and qTIP4P/F11 methods include nuclear quantum effects. “Ab initio” force fields are fitted to high level ab initio computations performed on molecular clusters. For instance, the pairwise nonpolarizable MCY12 potential is fitted to the water dimer energy using configuration interaction with single excitations (CIS). In non-polarizable force fields, polarization (many-body) effects are commonly introduced implicitly through fitting of the point charges, resulting in dipole moments that are overestimated. Three-body interactions may also be explicitly included. One example is the explicit three-body (E3B) method13-15 in which the three-body energy is fitted with parametrized exponential functions. One drawback of the above approaches is that the charge distribution of the molecules does not properly adjust to the molecular environment; this deficiency led to the development of polarizable force fields.16-18 Polarization can be modeled with, for instance, the induced point dipole (IPD) scheme, the Drude oscillator scheme or the fluctuating charge scheme. In the IPD scheme, polarizable point dipoles assigned to the molecules interact with the surrounding electric field and are computed in a self-consistent manner. Examples of force fields using this type of approach are the atomic multipole optimized energetics for biomolecular applications (AMOEBA),19-21 the Thole-type model (TTM)22-27 and the distributed point polarizable model (DPP).28,29 In the Drude oscillator scheme, the interaction between harmonic oscillators located at atomic centers is computed. Each oscillator is represented by a negatively charged, massless

2 ACS Paragon Plus Environment

Page 2 of 29

Page 3 of 29

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


particle attached with a spring to a heavy atomic core particle.30 Examples of force fields based on this scheme are SWM,31,32 COS,33,34 and BK3.35 In the fluctuating charge approach, the point charges that are used to model the molecules are not fixed: they are optimized by minimizing the total electrostatic energy. Examples include TIP4P-FQ and SPC-FQ.36 The polarizable force fields listed above are fitted to ab initio and/or experimental data. More sophisticated approaches combining polarization and explicit three-body interactions such as CC-pol,37-39 WHBB40 and HBB2-pol41 are also under active development. Some force fields are derived based on perturbation theory and compute the different contributions to the interaction energy. These methods have the advantage of providing some insight on the physical origin of the intermolecular interactions. Examples of such force fields include the anisotropic site potential (ASP),42,43 the sum of interaction between fragments ab initio (SIBFA)44 and the effective fragment potential (EFP) method.45 The EFP method was originally developed to model water solvation effects (EFP1).46,47 In the EFP1 method, water molecules are rigid and considered explicitly. The total interaction energy is divided into three contributions: E EFP1 = ECoul + E pol + E remainder

(1)

In eq.1, ECoul represents the Coulombic interaction energy computed using the multipole moments obtained with the Stone48,49 distributed multipole analysis (DMA) method. Charge penetration effects are included by using a Gaussian damping function.45 Epol is the polarization energy. Induced dipoles of the fragments in the field exerted by other fragments, or by a region that is represented by a quantum mechanical wave function, are iterated to self-consistency. Epol therefore includes many-body effects. Eremainder is a fitted remainder term. If the reference



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 29

interaction energy is restricted Hartree-Fock (RHF), Eremainder includes exchange repulsion and charge transfer. If the reference interaction energy is density functional theory, Eremainder includes short range electron correlation as well.50 The water potentials in EFP1 were generated using the DH(d,p) basis set.45 The EFP method has been generalized to model all molecules by relying on first principles derivations and avoiding empirical fitting. The general EFP method is called EFP2, or simply EFP. In EFP2, the intermolecular interaction energy is given by: E EFP 2 = ECoul + E pol + E exrep + E disp + Ect

(2)

In eq.2, ECoul and Epol are computed in a similar fashion to that in EFP1. Gaussian, exponential and overlap functions can be used for the screening of the Coulomb energy.51,52 Gaussian and exponential screenings are available for the polarization energy.51 Eexrep represents the exchange repulsion energy, computed from kinetic energy integrals and intermolecular localized molecular orbital (LMO) overlap integrals.53,54 Ect represents the charge transfer energy, obtained from the second order perturbative treatment of the intermolecular interactions on the antisymmetrized product of the SCF wavefunctions of two monomers.55 Edisp is the dispersion energy, computed using the distributed dynamic polarizabilities of the monomers and an overlap or Tang-Toennies damping function.51,56-58 In the EFP2 method, the R-6 and R-7 dispersion energy terms and their gradients have been derived and implemented.56-58 The R-6 dispersion energy term has also been implemented as a first-principles correction to density functional theory and Hartree-Fock energies.59 The optimization of a water trimer with EFP2 has revealed that the R-7 dispersion energy term yields longer intermolecular distances relative to the prediction using the R-6 term


Page 5 of 29

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


alone.57 In addition, it was shown that the R-7 term represents 20% of the magnitude of the dispersion energy for this system.57 Contrary to the EFP1 parameters which are internally stored, EFP2 parameters must be generated separately by the user, therefore requiring an extra step. In addition, while EFP2 is a more complete and physically accurate force field, EFP2 calculations are more expensive than EFP1 calculations. This is mostly due to the large computational cost of the EFP2 exchange repulsion and charge transfer terms, which are included in the parametrized remainder term in EFP1. The only interaction term that is completely missing in EFP1 is dispersion. An improved EFP1 force field which includes dispersion, called EFP1-D has therefore been developed. Just like EFP1, the EFP1-D parameters are internally stored, making EFP1-D just as easy to use as EFP1. The paper is organized as follows. First, the derivation and implementation of the new dispersion-corrected EFP1 method, EFP1-D, is described. Computational details are reported next. In the results section, the EFP1-D method and the EFP2 method are applied to the optimization of water clusters (H2O)n (n=2-6) and the importance of the R-7 term for these clusters is assessed. The final section concludes. Methods In the EFP1-D method, the dispersion energy, as computed in EFP2, is added to the EFP1 energy (eq.1). The total EFP1-D energy is thus given by: E EFP1− D = ECoul + E pol + Eremainder + Edisp

(3)

It should be emphasized that the remainder term in eq.3 is fitted to RHF,45 not DFT50 to avoid double counting of electron correlation. Since the importance of dispersion forces in water



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 6 of 29

clusters, and the ability of EFP1-D to capture these interactions accurately, are assessed in this paper, a brief description of the dispersion energy is given here. The total dispersion energy is given by: ∞

Edisp = ∑ En

(4)

n=6

In eq.4, E n is the nth order dispersion energy between two fragments, summed over all fragment pairs. In practice, this expression needs to be truncated. The E6 dispersion energy between two fragments A and B corresponds to the interaction between the induced-dipoles of the fragments and is expressed as:56

E6 = f 6 E6,0

(5)

In eq.5, f 6 is a damping function. The undamped dispersion energy E 6 , 0 can be derived from Rayleigh-Schrödinger perturbation theory.49 Within the isotropic approximation, it is given by:56 ∞

E6 ,0 = −

3h

π

∑∑

∫ α (iω )α (iω ) dω k

j

(6)

0

Rkj6

k∈A j∈B

In eq.6, Rkj represents the distance between the localized molecular orbital centroids k and j on j

fragments A and B respectively.

k

α ( iω ) and α ( iω ) represent the trace of the dynamic

polarizability tensors for LMOs k and j. Using the Gauss-Legendre quadrature formula, E 6 , 0 becomes:56 k

E6 ,0

( ) ( ) j

2ω 0 α iω i α iω i = − ∑ ∑ ∑Wi 2 π k∈A j∈B i=1 Rkj6 1− t 3h

12

(

i

)


(7)

Page 7 of 29

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


In eq.7, Wi and ti are the Gauss-Legendre weighting factor and abscissa, respectively. The damping function f6 is given by:57,59

−2ln S ) ( f ( k , j ) = 1− S ∑ 6

kj

2 kj

6

n/ 2

(8)

n!

n=0

In eq.8, Skj is the overlap integral between the LMOs k and j. The E 7 dispersion energy term corresponds to the interaction between induced dipoles and induced quadrupoles and has a form that is analogous to the E 6 term:

E7 = f 7 E7,0

(9)

The undamped dispersion energy E 7 , 0 is given by:58 ∞

x ,y ,z h kj E7,0 = − ∑ ∑ ∑ Tαβkj Tγσκ d ω α αγk iω Aβj ,σκ iω − α βκj iω Aαk ,γσ iω  ∫ 3π k∈A j∈B αβγσκ 0

( )

( )

( )

( )

(10)

Tαβkj and T γσk jκ are the second and third order electrostatic tensors given by:

Tαβkj = ∇α ∇ β

1 Rkj

kj Tγσκ = ∇γ ∇σ ∇κ

(11.a)

1 Rkj

(11.b)

Al represents the dynamic dipole-quadrupole polarizability tensor at LMO l. The damping

function f7 can be derived similarly to f6:57



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

(−2ln S ) f ( k , j ) = 1− S ∑ 7

7

n/ 2

kj

2 kj

n!

n=0

Page 8 of 29

(12)

Analytic gradients have been derived and implemented for both the E6 and E7 terms.56,57 In the original implementation, the E6 term was scaled by a factor of 4/3. This ad hoc correction was meant to empirically include higher order terms. This factor was removed for this study in order to provide a rigorous analysis of the effects of the exact E6 and E7 dispersion terms on water cluster geometries and energies. It should also be emphasized that the E7 term shows some angular dependence as this term is anisotropic. On the other hand, the isotropic approximation was used for E6 because the isotropic C6 coefficients can be compared directly with experimental data.56 In addition, this approach leads to considerable time savings as the full dynamic dipoledipole polarizability tensor does not need to be computed at each optimization step.56 Finally, it should be noted that the computation of the E7 term does not yield a significant increase in computational time compared to the E6 term alone. In EFP, the exchange-repulsion and charge transfer terms are by far the most computationally demanding.54,55 The most expensive components of the dispersion term are the overlap integrals used in the damping function. These integrals are computed for the exchange-repulsion term, stored and reused for dispersion. The cost of the E6,0 and E7,0 computations is very small in comparison. Computational details All calculations were performed using the GAMESS software package.60,61 The internally stored EFP1 geometry is used to model fragment water molecules: O-H bond lengths are 0.94 Å and the H-O-H bond angle is 106.7°. The EFP1-D water potentials were generated at the RHF/DH(d,p) level of theory as they were for EFP1.45 The recommended 8 ACS Paragon Plus Environment

Page 9 of 29

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


6-311G++(3df,2p) basis was used for EFP2. For EFP2 calculations, the localized orbitals of water were generated with the Boys procedure.62 An exponential screening was used for the Coulomb energy and Gaussian screening was used for polarization. Water clusters (H2O)n (n=26) optimized at the MP2/aug-cc-pVDZ level of theory by Temelso et al63 were used as starting structures for geometry optimizations. All of these water clusters were optimized with the EFP1, EFP1-D and EFP2 methods. Semi-numerical Hessians were computed for the optimized structures to confirm that they are minima on their respective potential energy surfaces. In order to evaluate the contribution of the E7 energy term to water cluster geometries, two separate optimizations were performed with both the EFP1-D and EFP2 methods. They are labeled EFP1D(E6+E7), EFP1-D(E6), and similarly for EFP2. The notation (E6) means that only the E6 term is included in the dispersion energy. The notation (E6+E7) means that both E6 and E7 dispersion energy contributions are included. Henceforth, the label EFP1-D will be used when referring to both EFP1-D(E6+E7) and EFP1-D(E6). Similarly, the label EFP2 will be used when referring to both EFP2(E6+E7) and EFP2(E6). No other contributions to the dispersion interaction energy were included in this work. Root mean squared deviations (RMSD) between EFP geometries and MP2 geometries were computed using the VMD software.64 RMSDs are given by:

∑ ( r ( EFP ) − r ( MP2)) N

RMSD =

i

2

i

i=1

N

(13)

In eq.13, ri(EFP) is the position of atom i of the EFP optimized structure (EFP refers to EFP1, EFP1-D(E6+E7), EFP1-D(E6), EFP2(E6+E7) or EFP2(E6)). ri(MP2) is the position of atom i of



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 29

the MP2 optimized structure.63 ri(EFP) and ri(MP2) are defined relative to the center of mass of the cluster. N is the number of atoms. In the following section, the percentage of the E7 contribution to the total dispersion energy in EFP1-D(E6+E7) and EFP2(E6+E7) calculations is computed as

E7 × 100 . Mean E6 + E7

absolute errors (MAEs) of EFP interaction energies are computed using the following formula: N

∑ E ( EFP) − E ( CCSD(T ) / CBS ) i

MAE =

i

i =1

N

(14)

In eq.14, Ei(EFP) and Ei(CCSD(T)/CBS) are the total intermolecular interaction energies computed with EFP and CCSD(T)/CBS63 for each water cluster i. EFP refers to EFP1, EFP1D(E6+E7), EFP1-D(E6), EFP2(E6+E7) or EFP2(E6)). N is the total number of water clusters investigated in this work (N=20).

Results and discussion I)

Water dimer.

The interaction energy of the water dimer (Figure 1) at the MP2/aug-cc-pVDZ geometry was computed with the EFP1, EFP1-D(E6+E7), EFP1-D(E6), EFP2(E6+E7) and EFP2(E6) methods. The results, shown in Table 1, are compared with the CCSD(T) energy extrapolated to the complete basis set limit (CCSD(T)/CBS), also computed at the MP2/aug-cc-pVDZ geometry.63 Table 1 also shows the energy decomposition analysis (EDA), performed by Su and Li at the CCSD(T)/aug-cc-pVQZ level of theory, with counterpoise correction, on the MP2/aug-cc-pVQZ


Page 11 of 29

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


optimized dimer.65 In the EDA calculation performed by Su and Li, the Coulomb, polarization, exchange and repulsion energies were obtained from the Hartree-Fock wavefunction while the dispersion term was computed as the difference between the CCSD(T) and Hartree-Fock interaction energies.65 It is therefore likely that the EDA approach overestimates the dispersion energy.

Figure 1. Water dimer optimized at the MP2/aug-cc-pVDZ level of theory.63

Table 1. EFP contributions to the water dimer energy at the MP2/aug-cc-pVDZ optimized geometry in kcal/mol. EFP2 Level of EFP1-D EFP1-D EFP2 EDAd EFP1 (E6) (E6+E7) theory (E6+E7) (E6) Coulomb -7.25 -7.25 -7.25 -8.14 -8.14 -8.41 Exchange 2.94a 2.94a 2.94a 5.59b 5.59b 7.16c Repulsion Polarization -0.69 -0.69 -0.69 -1.18 -1.18 -2.38 e Dispersion 0.00 -0.29 -0.53 -0.51 -1.01 -1.33 E6 0.00 -0.53 -0.53 -1.01 -1.01 N/A E7 0.00 0.24 N/A 0.50 N/A N/A Charge 0.00 0.00 0.00 -0.47 -0.47 N/A Transfer Total interaction -5.00 -5.29 -5.53 -4.71 -5.21 -4.95 energy a. Remainder term contains exchange repulsion and charge transfer. b. Exchange-repulsion from ref[53]. c. Sum of EDA exchange and repulsion energies from Ref[65] d. Obtained at the CCSD(T)/aug-cc-pVQZ //MP2/aug-cc-pVQZ level in ref[65] e. (E6+E7) for EFP1-D(E6+E7) and EFP2(E6+E7). E6 only for EFP1-D(E6) and EFP2(E6).



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 EFP1 method yields a total interaction energy of -5.00 kcal/mol, which is nearly identical to CCSD(T)/CBS value of -5.03 kcal/mol.63 Since the EFP1 remainder term is fitted to RHF, it does not account for dispersion. However, the remainder term does account for both exchange repulsion and charge transfer. The total EFP1 energy (which includes the Coulomb, polarization and “remainder” terms) and the sum of the Coulomb, exchange repulsion and polarization energies from the EDA both add up to the Hartree-Fock energy. However, these two quantities are not equal since they were computed with different basis sets (DH(d,p) for EFP1 and aug-cc-pVQZ for the EDA data). In fact, the total EFP1 energy is 1.37 eV larger than the sum of the Coulomb, exchange repulsion and polarization terms from the EDA. This difference is approximately equal to the EDA dispersion energy (-1.33 kcal/mol). The fact that the EFP1 energy is nearly identical to the CCSD(T)/CBS energy therefore probably results from a cancelation of errors. It is important to emphasize that the individual EFP1 and EDA energy contributions cannot strictly be compared, because of the fitted remainder term and since the energy decomposition schemes are not unique. The inclusion of dispersion interactions in EFP1-D(E6+E7) and EFP1-D(E6) results in a slight overestimation of the total interaction energy by 0.26 and 0.50 kcal/mol, respectively, in comparison to CCSD(T)/CBS. Both of these values remain well within accepted chemical accuracy (~1 kcal/mol). The E7 dispersion energy is positive, so this term is repulsive. The total dispersion energies obtained with the EFP1-D(E6+E7) and EFP1-D(E6) are 1.04 and 0.80 kcal/mol lower than the EDA value respectively. Since the EDA dispersion energy is defined as the difference between CCSD(T)/aug-cc-pVQZ and HF energies, thereby assuming that all correlation effects are due to dispersion,65 the EDA dispersion energy is expected to be too high. The total (absolute) EFP2(E6+E7) energy is too small by 0.32 kcal/mol relative to the


Page 12 of 29

Page 13 of 29

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


CCSD(T)/CBS value, whereas the EFP2(E6) energy is too large by 0.18 kcal/mol. The difference between the dispersion energies obtained with EFP1-D and EFP2 arises from the different basis sets used by the two methods. The E7 contribution represents 31% and 33%, respectively, of the total EFP1-D(E6+E7) and EFP2(E6+E7) dispersion energy, a significant contribution. Now, the impact of dispersion on the predicted geometries of the water dimer is considered. The water dimer was optimized with the five EFP methods described previously. The interaction energies at the optimized geometries are reported in Table 2. The oxygen-oxygen distance and the O(1)-H(2)-O(4) angle (See Figure 1) are reported in Table 3. Overall, the addition of dispersion to EFP1 yields geometries that are closer to the MP2/aug-cc-pVDZ geometries, compared to EFP1 alone. The positive E7 energy contribution yields intermolecular distances that are 0.04 Å longer for EFP1-D(E6+E7) and 0.05 Å longer for EFP2(E6+E7) than for EFP1-D(E6) and EFP2(E6), respectively. Accurate intermolecular orientations are important to properly model the dynamics of water. The O(1)-H(2)-O(4) angles obtained with EFP1, EFP1D(E6+E7) and EFP1-D(E6) are overestimated by more than 4 degrees in comparison with MP2/aug-cc-pVDZ. On the other hand, the angles obtained with the EFP2 methods are in excellent agreement with MP2. The geometry obtained with EFP2(E6+E7) is overall the closest to the MP2 structure. The EFP1-D(E6+E7) and EFP1-D(E6) interaction energies are larger (more negative) than the EFP1 interaction energies due to the overall attractive dispersion energy. The E7 dispersion energy represents nearly 40% and 50% of the magnitude of E6 for EFP1-D(E6+E7) and EFP2(E6+E7), respectively. Overall, the E7 dispersion energy contribution significantly affects the optimized dimer geometry and energy. Note also that the different geometries lead to differences in other energy terms, not just the dispersion interaction energy.



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

Table 2. Energies in kcal/mol at EFP optimized geometries EFP1-D EFP2 EFP1 EFP1-D(E6) EFP2(E6) (E6+E7) (E6+E7) Coulomb -6.65 -6.81 -7.73 -7.15 -8.48 Exchange 2.15 2.34 5.01 2.77 6.09 repulsion Polarization -0.60 -0.63 -1.08 -0.69 -1.26 E6 0.00 -0.46 -0.95 -0.51 -1.07 E7 0.00 0.18 0.45 0.00 0.00 Charge N/A N/A -0.43 N/A -0.50 Transfer Total interaction -5.10 -5.36 -4.72 -5.57 -5.22 energy

Table 3. Atom distances and angles in the water dimer optimized at different levels of theory. MP2/aug- EFP1 EFP1-D EFP2 EFP1-D(E6) EFP2(E6) cc-pVDZ (E6+E7) (E6+E7) O(1)-O(4) (Å) 2.92 2.99 2.97 2.95 2.93 2.90 O(1)-H(2)-O(4) 171.3 176.4 177.87 171.00 175.05 171.89 angle (in degrees)

II)

(H2O)n (n=3-6) clusters

Water trimers, tetramers, pentamers and hexamers were analyzed with EFP1, EFP1-D(E6+E7), EFP1-D(E6), EFP2(E6+E7) and EFP2(E6). The water clusters optimized at the MP2/aug-ccpVDZ level of theory63 are shown in Figure 2A-2T.


Page 14 of 29

Page 15 of 29

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


Figure 2. A) Up-up-down (UUD) trimer B) Up-up-up (UUU) trimer C) Cyclic Ci tetramer D) Pyramidal (Py) tetramer E) Cyclic S4 tetramer F) Pentamer cage A (CA_A) G) Pentamer cage B (CA_B) H) Pentamer cage C(CA_C) I) Cyclic pentamer (CYC) J) Pentamer fused ring A (FR_A) K) Pentamer fused ring B (FR_B) L) Pentamer fused ring C (FR_C) M) Hexamer prism (PR) N) Hexamer cage (CA) O) Hexamer BAG P) Hexamer cyclic chair (CC) Q) Hexamer book1 (BK1) R) Hexamer book 2 (BK2) S) Hexamer cyclic boat 1 (CB1) T) Hexamer cyclic boat 2 (CB2)

The EFP1, EFP1-D(E6+E7), EFP1-D(E6), EFP2(E6+E7) and EFP2(E6) relative energies of the water trimers, tetramers, pentamers and hexamers at the MP2/aug-cc-pVDZ optimized geometries are shown in Table 4. The mean absolute errors (MAEs) of the total interaction energies for the different EFP levels are also displayed. The total interaction energies used to compute the MAEs are listed in Table S1 of the supporting information. The EFP1-D(E6) total interaction energies are very similar to the CCSD(T)/CBS interaction energies, as shown by the 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 29

small mean absolute error of 0.24 kcal/mol. On the other hand, the largest MAEs are obtained for EFP1(4.52 kcal/mol). The E6 dispersion correction therefore overall significantly improves the total interaction energies for these systems. EFP1-D(E6+E7), EFP2(E6) and EFP2(E6+E7) have similar MAE values between 1.58 and 1.86 kcal/mol. Table 4. Relative water cluster energies and mean absolute error (MAE) of the total interaction energies in kcal/mol at the MP2/aug-cc-pVDZ optimized geometries. MP2/augEFP1-D EFP2 EFP1-D EFP2 EFP1 ccCCSD(T)/CBSa (E6+E7) (E6+E7) (E6) (E6) a pVDZ trimer_UUD 0.00 0.00 0.00 0.00 0.00 0.00 0.00 trimer_UUU 0.90 0.98 0.26 0.95 1.17 0.57 0.61 tetramer_S4 0.00 0.00 0.00 0.00 0.00 0.00 0.00 tetramer_Ci 1.00 1.05 0.67 1.05 0.71 0.82 0.85 tetramer_Py 3.68 2.91 3.22 3.62 4.83 3.72 3.55 pentamer_CYC 0.00 0.00 0.00 0.00 0.00 0.00 0.00 pentamer_FR_B 1.17 0.49 1.06 1.01 2.13 1.21 1.13 pentamer_CA_C 1.55 0.34 1.42 1.15 3.10 1.35 1.32 pentamer_CA_A 1.74 0.50 1.21 1.31 2.88 1.53 1.47 pentamer_CA_B 2.45 1.15 2.95 2.15 4.96 2.23 2.19 pentamer_FR_A 3.44 2.55 3.25 3.32 4.73 2.81 2.88 pentamer_FR_C 4.20 3.39 3.92 4.18 5.45 3.61 3.57 hexamer_PR 0.00 0.00 0.00 0.00 1.23 0.00 0.00 hexamer_CA 0.44 0.78 0.57 0.55 1.30 0.25 0.21 hexamer_BK1 0.66 1.87 0.22 1.25 0.14 0.75 0.63 hexamer_BK2 0.92 2.15 0.14 1.48 0.00 1.00 1.00 hexamer_CC 0.79 2.55 0.77 1.92 0.73 2.08 1.54 hexamer_BAG 1.75 2.93 1.38 2.26 1.25 1.47 1.55 hexamer_CB1 2.14 3.85 1.92 3.31 2.04 3.06 2.57 hexamer_CB2 2.04 3.71 2.19 3.23 2.43 3.18 2.63 b MAE 4.52 1.66 1.86 0.24 1.58 1.73 0 a. Ref63 b. Total EFP and reference CCSD(T)/CBS interaction energies used to compute the MAE are shown in the Supporting Information.

Trimers. The two trimers investigated are the up-up-down (UUD) isomer (Figure 2A) and the up-up-up (UUU) isomer (Figure 2B). Up-up-up refers to the non-bonded hydrogen atoms in each


Page 17 of 29

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


molecule all pointing in the same direction outside the plane formed by the three oxygen atoms. Up-up-down means that one non-bonded hydrogen atom is pointing in the opposite direction. All of the EFP methods correctly predict the energy order of the water trimer isomers. The EFP1, EFP1-D(E6+E7), EFP1-D(E6) and EFP2(E6) methods predict similar relative trimer energies, with the UUU isomer nearly 1 kcal/mol higher in energy than the UUD isomer. The relative EFP2(E6+E7) energy is slightly lower (0.3 kcal/mol) but in the same direction. Tetramers. The energies of two cyclic tetramers (one with Ci and the other with S4 symmetry, shown in Figure 2C and 2D respectively) and the pyramidal tetramer are investigated (Figure 2E). All EFP calculations reproduce the correct order of the isomers with relative energies overall that are comparable to those obtained with CCSD(T)/CBS. Pentamers. Seven water pentamers are studied: the three cages CA_A, CA_B and CA_C, the cyclic pentamer (CYC), and the three fused rings FR_A, FR_B and FR_C. These isomers are shown in Figure 2F- 2L. All EFP calculations identify the cyclic isomer as the lowest energy structure. EFP1 and EFP1-D(E6) reproduce the CCSD(T)/CBS order of all of the pentamers with relative energies differing from the CCSD(T)/CBS values by at most ~0.6 kcal/mol. EFP1D(E6+E7) has two isomers out of order, but only by ~0.1 kcal/mol. The EFP1-D(E6+E7) relative energies differ by up to 1 kcal/mol from the CCSD(T)/CBS values. The EFP2(E6) method does not reproduce the CCSD(T)/CBS order of the isomers, and the relative energies are overestimated by up to 2.8 kcal/mol. Adding the E7 contribution to EFP2 mostly corrects this deficiency. Hexamers. Eight water hexamers are studied: the prism (PR), the cage (CA), the bag (BAG), the cyclic chair (CC), two books (BK1 and BK2) and two cyclic boats (CB1 and CB2). These isomers are displayed in Figure 2M-2T. The CCSD(T)/CBS energy range for the hexamers is 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

~2.6 kcal/mol from lowest to highest energy, with the four lowest energy isomers all within 1 kcal/mol of each other. The EFP1-D methods reproduce the CCSD(T)/CBS order except for the cyclic boat 1 and cyclic boat 2 isomers, which are in the reverse order. On the other hand, the book 2 and cyclic chair isomers are also reversed with EFP1. EFP1-D(E6+E7) and EFP1-D(E6) overestimate the relative energies compared to CCSD(T)/CBS by up to 1.4 and 0.7 kcal/mol respectively. Contrary to EFP2(E6+E7), EFP2(E6) does not identify the prism as the lowest energy hexamer; however, the overall energy order is in good agreement with that predicted by CCSD(T)/CBS and the specific differences are small. The EFP2(E6+E7) method also predicts an overall energy range that is similar to that obtained from CCSD(T)/CBS. Summary. Overall, adding dispersion to EFP1 does not significantly improve the water isomer relative energies but does improve total interaction energies as shown by the decreasing mean absolute errors. The lack of improvement of the relative energies might originate from error cancelation, as noted above for the dimer. The EFP1-D(E6) method yields relative energies that are closer to the CCSD(T)/CBS values than are the EFP1-D(E6+E7) values. The addition of the E7 term should, in principle, improve the dispersion energy. However, since the signs of the E6 and E7 contributions are opposite, it is likely that higher order dispersion terms and/or manybody contributions66-68 are needed to attain CCSD(T)/CBS accuracy. For EFP2, the E7 term clearly improves the relative energies of water clusters. Even so, the higher order and/or manybody terms may be important. All of the E7 dispersion contributions for all of the water clusters examined here are positive (repulsive), while the E6 values are negative (attractive). The magnitudes of the E6 and E7 dispersion energy terms for EFP1-D(E6+E7) and EFP2(E6+E7) are shown in Figure 3. The two methods yield different dispersion energies since the EFP1-D(E6+E7) parameters were generated 18 ACS Paragon Plus Environment

Page 18 of 29

Page 19 of 29

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


using the DH(d,p) basis set, whereas the EFP2(E6+E7) parameters were generated using the 6311G++(3df,2p) basis set. As illustrated in Figure 3, the 6-311G++(3df,2p) basis set used for EFP2 calculations yields E6 and E7 dispersion energies that are twice as large as those obtained with the DH(d,p) basis set. The E7 term contributes a non-negligible 18-33% of the total dispersion energy.

Figure 3. Magnitudes of the E6 and E7 contributions to the dispersion energy for water clusters optimized at the MP2/aug-cc-pVDZ level of theory.

The foregoing discussion assessed the performance of the EFP methods at the MP2 geometries. It is important to consider each method at its own optimized geometries. Therefore, all of the water clusters discussed above were optimized using the respective methods: EFP1, EFP1-D(E6+E7), EFP1-D(E6), EFP2(E6+E7) and EFP2(E6). The root mean square deviations (RMSDs) between the optimized geometries that were obtained using the different EFP methods



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 the MP2/aug-cc-pVDZ optimized molecules are shown in Figure 4 and listed in Table S2 of the supporting information. Figure 5 presents the average O-O distance (also defined as the intermolecular distance) between neighboring water molecules for all of the clusters optimized with the five EFP methods and MP2. The average O-O distances between neighboring molecules for all isomers are listed in Table S3 of the supporting information. Overall, RMSDs and intermolecular distances decrease in the order EFP1>EFP1-D(E6+E7)>EFP1-D(E6). Adding the positive E7 dispersion energy contribution results in larger intermolecular distances in comparison to E6 alone, as expected. The local minimum corresponding to the pyramidal tetramer could not be optimized with any of the EFP methods. Instead, the cyclic Ci isomer was obtained. The fused ring C pentamer optimized with EFP1 and EFP1-D(E6+E7) is very distorted, as shown by the large RMSDs. In fact, the nearly-planar structure formed by four of the oxygen atoms (O1, O4, O7 and O10) in the MP2/aug-cc-pVTZ optimized cluster is broken when this cluster is optimized with EFP1 or EFP1-D(E6+E7) (Figure S1 of the supporting information). For EFP2(E6), the neighboring O-O distances tend to be underestimated. Adding the repulsive E7 term therefore increases intermolecular distances, which then become larger than those obtained with MP2. Exceptions to this trend are the cage B and cage C pentamers, for which the average O-O distance obtained with EFP2(E6+E7) is shorter than the one obtained for EFP2(E6). In principle, including the E7 contribution should provide a more accurate description of dispersion forces and therefore more reliable geometries. It is possible that higher order terms and/or manybody contributions need to be included. Overall, the data in Figures 4 and 5 shows that the E7 term can significantly affect the geometry of water clusters.


Page 20 of 29

Page 21 of 29

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


Figure 4. RMSDs of the EFP optimized structures compared with the MP2/aug-cc-pVDZ optimized clusters. The pyramidal tetramer is not included since it could not be optimized with EFP.

Figure 5. Average O-O distance (Å) of the water clusters. Each bar represents an average between all isomers. The pyramidal tetramer is not included in the average since it could not be optimized with EFP.



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 relative energies calculated for the EFP optimized geometries are shown in Table 5, where they are compared with the MP2/aug-cc-pVDZ and CCSD(T)/CBS relative energies obtained using the MP2/aug-cc-pVDZ optimized structure. The mean absolute errors (MAEs) of the total interaction energies for the different EFP optimized molecules are also displayed (CCSD(T)/CBS is again used as the reference value). The total interaction energies used to compute the MAEs are listed in Table S4 of the supporting information. The MAE values obtained for EFP1-D(E6+E7) and EFP1-D(E6) are 0.80 and 0.66 kcal/mol respectively, much smaller than for EFP1 (2.87 kcal/mol). Therefore, adding dispersion to EFP1 yields again improved total interaction energies. The MAE obtained for EFP2(E6+E7) is about 1 kcal/mol lower than the one obtained for EFP2(E6). The energy difference between the UUD and UUU trimers is improved upon optimization for all EFP methods, getting closer to those obtained using CCSD(T)/CBS and MP2/aug-cc-pVDZ methods. The relative energies between the cyclic Ci and S4 tetramers are in better agreement with CCSD(T)/CBS upon optimization for all EFP methods. The pyramidal tetramer could not be optimized with any of the EFP methods, as explained previously. The relative EFP1 energies of the pentamers and hexamers that were obtained after geometry optimization are in excellent agreement with CCSD(T), except for the fused ring C pentamer, whose energy relative to the cyclic isomer is underestimated by over 2 kcal/mol. This is most likely due to the fact that this structure is highly distorted compared to the MP2 optimized geometry, as shown previously by the large RMSD. The relative ordering of the pentamers and hexamers computed with EFP2(E6+E7), EFP1-D(E6) and EFP1-D(E6+E7) remains mostly the same after optimization. One exception is again the fused ring C pentamer optimized with EFP1-D(E6+E7), whose relative energy to the cyclic isomer is underestimated due to the large distortions mentioned previously.


Page 22 of 29

Page 23 of 29

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


Table 5. Relative water cluster energies and mean absolute error (MAE) of the total interaction energies in kcal/mol at their optimized geometries. EFP1-D EFP2 EFP1-D EFP2 MP2/augEFP1 CCSD(T)/CBSa (E6+E7) (E6+E7) (E6) (E6) cc-pVDZ trimer_UUD 0.00 0.00 0.00 0.00 0.00 0.00 0.00 trimer_UUU 0.77 0.79 0.56 0.87 0.62 0.57 0.61 tetramer_S4 0.00 0.00 0.00 0.00 0.00 0.00 0.00 tetramer_Ci 0.86 0.86 0.79 0.96 0.87 0.82 0.85 tetramer_Py N/A N/A N/A N/A N/A 3.72 3.55 pentamer_CYC 0.00 0.00 0.00 0.00 0.00 0.00 0.00 pentamer_FR_B 0.90 0.46 0.59 0.60 1.65 1.21 1.13 pentamer_CA_C 1.22 0.36 1.18 0.65 2.49 1.35 1.32 pentamer_CA_A 1.51 0.58 0.83 0.89 2.51 1.53 1.47 pentamer_CA_B 1.90 0.95 2.22 1.42 2.35 2.23 2.19 pentamer_FR_A 2.78 2.29 2.16 2.45 3.46 2.81 2.88 pentamer_FR_C 1.49 0.60 2.28 3.06 3.53 3.61 3.57 hexamer_PR 0.00 0.00 0.00 0.00 0.86 0.00 0.00 hexamer_CA 0.52 0.80 0.67 0.68 1.11 0.25 0.21 hexamer_BK1 0.66 1.58 0.56 1.36 0.20 0.75 0.63 hexamer_BK2 0.98 1.92 0.50 1.66 0.00 1.00 1.00 hexamer_CC 1.27 2.61 1.15 2.60 0.41 2.08 1.54 hexamer_BAG 1.86 2.77 1.68 2.49 1.29 1.47 1.55 hexamer_CB1 2.33 3.50 2.34 3.76 1.78 3.06 2.57 hexamer_CB2 2.32 3.50 2.45 3.74 1.85 3.18 2.63 b MAE 2.87 0.66 1.34 0.80 2.38 1.73 0 63 a. Ref. b. Absolute EFP and reference CCSD(T)/CBS energies used to compute the MAE are shown in the Supporting Information. Conclusions Water clusters containing between two and six molecules were optimized using five versions of the effective fragment potential method: EFP1, EFP1-D(E6+E7), EFP1-D(E6), EFP2(E6+E7) and EFP2(E6). The first version, EFP1, includes electrostatic and polarization energies. It also includes an energy term that is obtained by fitting to restricted Hartree-Fock calculations. This term therefore includes exchange repulsion and charge transfer energies. The second version, EFP1-D(E6+E7), includes the E6 and E7 dispersion interactions between the fragments, which are computed from first principles. EFP1-D(E6) only computes the E6 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

dispersion energy term. EFP2(E6+E7) and EFP2(E6) include electrostatic, polarization, exchangerepulsion, dispersion and charge transfer energies, all computed from first principles. EFP2(E6+E7), includes the E6 and E7 dispersion terms between the fragments whereas EFP2(E6) only computes the E6 dispersion energy term. The E6 term is negative (attractive) whereas the E7 term is positive (repulsive) for all clusters. The magnitude of the E7 term on average represents half of the magnitude of E6, which is not negligible. The E7 dispersion energy term generally yields increased intermolecular distances. The intermolecular distances decrease in the order: EFP1>EFP1-D(E6+E7)>EFP1D(E6). The structures obtained with EFP1-D(E6) are closest to the MP2/aug-cc-pVDZ optimized structures. The EFP2(E6) method yields smaller intermolecular distances than the MP2/aug-ccpVDZ method. On the other hand, the EFP2(E6+E7) method yields intermolecular distances that are larger than those obtained with the MP2/aug-cc-pVDZ method, because of the repulsive nature of the E7 interaction energy. The EFP2(E6+E7) relative energies are overall closer to the CCSD(T)/CBS values than the EFP2(E6) relative energies. In general, the E7 term can strongly affect the optimized structures and energies of water clusters and should not be neglected. Since this term is anisotropic, the fact that it is consistently positive can be a result of the relative orientation of the water molecules that favors the formation of hydrogen bonds. As of now, the EFP2(E6+E7) method is the method of choice for modeling water as it is parameter-free and yields geometries that are quite similar to MP2. EFP1D(E6) is also a reasonable alternative as it presents a significant improvement over EFP1 geometries and yields total interaction energies that are closer to CCSD(T)/CBS, although relative energies are slightly overestimated (by about 1 kcal/mol). The parameters for this


Page 24 of 29

Page 25 of 29

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


method are all stored within GAMESS, making this method easier to use than EFP2, which requires the user to generate the parameters.

Supporting Information. Water cluster energies at their MP2/aug-cc-pVDZ optimized geometries, RMSDs of the EFP optimized clusters in comparison with the MP2/aug-cc-pVDZ optimized clusters, average oxygen-oxygen distances in MP2 and EFP optimized clusters, EFP optimized cluster energies, the FR-C pentamer optimized with EFP1-D and MP2 are supplied in the Supporting Information. Acknowledgement. This work was supported by a U.S. National Science Foundation Software Infrastructure (SI2) grant, ACI – 1450217. Some of the calculations that are reported here were performed on Cyence, a computer cluster obtained with funds from a National Science Foundation Major Research Instrumentation grant.

References (1) Cisneros, G. A.; Wikfeldt, K. T.; Ojamäe, L.; Lu, J.; Xu, Y.; Torabifard, H.; Bartók, A. P.; Csányi, G.; Molinero, V.; Paesani, F. Modeling Molecular Interactions in Water: From Pairwise to ManyBody Potential Energy Functions. Chem. Rev. 2016, 116, 7501-7528. (2) 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. (3) 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. (4) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910-8922. (5) Rick, S. W. A Reoptimization of the Five-Site Water Potential (TIP5P) for Use with Ewald Sums. J. Chem. Phys. 2004, 120, 6085-6093. (6) Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F.; Hermans, J. In Intermolecular Forces, Proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry Held in Jerusalem, Israel, April 13–16, 1981; Pullman, B., Ed.; Springer Netherlands: Dordrecht, 1981. (7) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269-6271. 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

(8) Wu, Y.; Tepper, H. L.; Voth, G. A. Flexible Simple Point-Charge Water Model with Improved Liquid-State Properties. J. Chem. Phys. 2006, 124, 024503. (9) Lawrence, C. P.; Skinner, J. L. Flexible TIP4P Model for Molecular Dynamics Simulation of Liquid Water. Chem. Phys. Lett. 2003, 372, 842-847. (10) Paesani, F.; Iuchi, S.; Voth, G. A. Quantum Effects in Liquid Water from an Ab InitioBased Polarizable Force Field. J. Chem. Phys. 2007, 127, 074506. (11) Habershon, S.; Markland, T. E.; Manolopoulos, D. E. Competing Quantum Effects in the Dynamics of a Flexible Water Model. J. Chem. Phys. 2009, 131, 024501. (12) Matsuoka, O.; Clementi, E.; Yoshimine, M. CI Study of the Water Dimer Potential Surface. J. Chem. Phys. 1976, 64, 1351-1361. (13) Tainter, C. J.; Pieniazek, P. A.; Lin, Y.-S.; Skinner, J. L. Robust Three-Body Water Simulation Model. J. Chem. Phys. 2011, 134, 184501. (14) Tainter, C. J.; Shi, L.; Skinner, J. L. Reparametrized E3B (Explicit Three-Body) Water Model Using the TIP4P/2005 Model as a Reference. J. Chem. Theory Comput. 2015, 11, 2268-2277. (15) Kumar, R.; Skinner, J. L. Water Simulation Model with Explicit Three-Molecule Interactions. J. Phys. Chem. B 2008, 112, 8311-8318. (16) Antila, H. S.; Salonen, E. In Biomolecular Simulations: Methods and Protocols; Monticelli, L., Salonen, E., Eds.; Humana Press: Totowa, NJ, 2013. (17) Halgren, T. A.; Damm, W. Polarizable Force Fields. Current Opinion in Structural Biology 2001, 11, 236-242. (18) Xu, P.; Wang, J.; Xu, Y.; Chu, H.; Liu, J.; Zhao, M.; Zhang, D.; Mao, Y.; Li, B.; Ding, Y. et al. In Advance in Structural Bioinformatics; Wei, D., Xu, Q., Zhao, T., Dai, H., Eds.; Springer Netherlands: Dordrecht, 2015. (19) Ren, P.; Ponder, J. W. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. J. Phys. Chem. B 2003, 107, 5933-5947. (20) Laury, M. L.; Wang, L.-P.; Pande, V. S.; Head-Gordon, T.; Ponder, J. W. Revised Parameters for the AMOEBA Polarizable Atomic Multipole Water Model. J. Phys. Chem. B 2015, 119, 9423-9437. (21) Ren, P.; Ponder, J. W. Temperature and Pressure Dependence of the AMOEBA Water Model. J. Phys. Chem. B 2004, 108, 13427-13437. (22) Burnham, C. J.; Li, J.; Xantheas, S. S.; Leslie, M. The Parametrization of a Thole-Type AllAtom Polarizable Water Model from First Principles and Its Application to the Study of Water Clusters (N=2–21) and the Phonon Spectrum of Ice Ih. J. Chem. Phys. 1999, 110, 4566-4581. (23) Fanourgakis, G. S.; Xantheas, S. S. The Flexible, Polarizable, Thole-Type Interaction Potential for Water (TTM2-F) Revisited. J. Phys. Chem. A 2006, 110, 4100-4106. (24) Fanourgakis, G. S.; Xantheas, S. S. Development of Transferable Interaction Potentials for Water. V. Extension of the Flexible, Polarizable, Thole-Type Model Potential (TTM3-F, V. 3.0) to Describe the Vibrational Spectra of Water Clusters and Liquid Water. J. Chem. Phys. 2008, 128, 074506. (25) Burnham, C. J.; Xantheas, S. S. Development of Transferable Interaction Models for Water. Reparametrization of an All-Atom Polarizable Rigid Model (TTM2–R) from First Principles. J. Chem. Phys. 2002, 116, 1500-1510. (26) Burnham, C. J.; Xantheas, S. S. Development of Transferable Interaction Models for Water. Iv. A Flexible, All-Atom Polarizable Potential (TTM2-F) Based on Geometry Dependent Charges Derived from an Ab Initio Monomer Dipole Moment Surface. J. Chem. Phys. 2002, 116, 5115-5124. (27) Burnham, C. J.; Anick, D. J.; Mankoo, P. K.; Reiter, G. F. The Vibrational Proton Potential in Bulk Liquid Water and Ice. J. Chem. Phys. 2008, 128, 154519. (28) Defusco, A.; Schofield, D. P.; Jordan, K. D. Comparison of Models with Distributed Polarizable Sites for Describing Water Clusters. Mol. Phys. 2007, 105, 2681-2696. 26 ACS Paragon Plus Environment

Page 26 of 29

Page 27 of 29

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


(29) Kumar, R.; Wang, F.-F.; Jenness, G. R.; Jordan, K. D. A Second Generation Distributed Point Polarizable Water Model. J. Chem. Phys. 2010, 132, 014309. (30) Huang, J.; Lopes, P. E. M.; Roux, B.; MacKerell, A. D. Recent Advances in Polarizable Force Fields for Macromolecules: Microsecond Simulations of Proteins Using the Classical Drude Oscillator Model. J. Phys. Chem. Lett. 2014, 5, 3144-3150. (31) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell Jr, A. D. A Polarizable Model of Water for Molecular Dynamics Simulations of Biomolecules. Chem. Phys. Lett. 2006, 418, 245249. (32) Lamoureux, G.; MacKerell, A. D.; Roux, B. A Simple Polarizable Model of Water Based on Classical Drude Oscillators. J. Chem. Phys. 2003, 119, 5185-5197. (33) Yu, H.; Hansson, T.; van Gunsteren, W. F. Development of a Simple, Self-Consistent Polarizable Model for Liquid Water. J. Chem. Phys. 2003, 118, 221-234. (34) Yu, H.; van Gunsteren, W. F. Charge-on-Spring Polarizable Water Models Revisited: From Water Clusters to Liquid Water to Ice. J. Chem. Phys. 2004, 121, 9549-9564. (35) Kiss, P. T.; Baranyai, A. A Systematic Development of a Polarizable Potential of Water. J. Chem. Phys. 2013, 138, 204507. (36) Rick, S. W.; Stuart, S. J.; Berne, B. J. Dynamical Fluctuating Charge Force Fields: Application to Liquid Water. J. Chem. Phys. 1994, 101, 6141-6156. (37) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Predictions of the Properties of Water from First Principles. Science 2007, 315, 1249-1252. (38) Góra, U.; Cencek, W.; Podeszwa, R.; van der Avoird, A.; Szalewicz, K. Predictions for Water Clusters from a First-Principles Two- and Three-Body Force Field. J. Chem. Phys. 2014, 140, 194101. (39) Jankowski, P.; Murdachaew, G.; Bukowski, R.; Akin-Ojo, O.; Leforestier, C.; Szalewicz, K. Ab Initio Water Pair Potential with Flexible Monomers. J. Phys. Chem. A 2015, 119, 2940-2964. (40) Wang, Y.; Huang, X.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Flexible, Ab Initio Potential, and Dipole Moment Surfaces for Water. I. Tests and Applications for Clusters up to the 22Mer. J. Chem. Phys. 2011, 134, 094509. (41) Babin, V.; Medders, G. R.; Paesani, F. Toward a Universal Water Model: First Principles Simulations from the Dimer to the Liquid Phase. J. Phys. Chem. Lett. 2012, 3, 3765-3769. (42) Millot, C.; Stone, A. J. Towards an Accurate Intermolecular Potential for Water. Mol. Phys. 1992, 77, 439-462. (43) Hodges, M. P.; Stone, A. J.; Xantheas, S. S. Contribution of Many-Body Terms to the Energy for Small Water Clusters: A Comparison of Ab Initio Calculations and Accurate Model Potentials. J. Phys. Chem. A 1997, 101, 9163-9168. (44) Piquemal, J.-P.; Chevreau, H.; Gresh, N. Toward a Separate Reproduction of the Contributions to the Hartree−Fock and DFT Intermolecular Interaction Energies by Polarizable Molecular Mechanics with the SIBFA Potential. Journal of Chemical Theory and Computation 2007, 3, 824-837. (45) Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An Effective Fragment Method for Modeling Solvent Effects in Quantum Mechanical Calculations. J. Chem. Phys. 1996, 105, 1968-1986. (46) Day, P. N.; Pachter, R.; Gordon, M. S.; Merrill, G. N. A Study of Water Clusters Using the Effective Fragment Potential and Monte Carlo Simulated Annealing. J. Chem. Phys. 2000, 112, 20632073. (47) Chen, W.; Gordon, M. S. The Effective Fragment Model for Solvation: Internal Rotation in Formamide. J. Chem. Phys. 1996, 105, 11081-11090. (48) Stone, A. J. Distributed Multipole Analysis, or How to Describe a Molecular Charge Distribution. Chem. Phys. Lett. 1981, 83, 233-239. 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

(49)

Stone, A. The Theory of Intermolecular Forces; 2nd ed.; Oxford University Press: Oxford,

2013. (50) Adamovic, I.; Freitag, M. A.; Gordon, M. S. Density Functional Theory Based Effective Fragment Potential Method. J. Chem. Phys. 2003, 118, 6725-6732. (51) Slipchenko, L. V.; Gordon, M. S. Damping Functions in the Effective Fragment Potential Method. Mol. Phys. 2009, 107, 999-1016. (52) Slipchenko, L. V.; Gordon, M. S. Electrostatic Energy in the Effective Fragment Potential Method: Theory and Application to Benzene Dimer. J. Comput. Chem. 2007, 28, 276-291. (53) Jensen, J. H.; Gordon, M. S. An Approximate Formula for the Intermolecular Pauli Repulsion between Closed Shell Molecules. Mol. Phys. 1996, 89, 1313-1325. (54) Li, H.; Gordon, M. Gradients of the Exchange-Repulsion Energy in the General Effective Fragment Potential Method. Theor. Chem. Acc. 2006, 115, 385-390. (55) Li, H.; Gordon, M. S.; Jensen, J. H. Charge Transfer Interaction in the Effective Fragment Potential Method. J. Chem. Phys. 2006, 124, 214108. (56) Adamovic, I.; Gordon, M. S. Dynamic Polarizability, Dispersion Coefficient C6 and Dispersion Energy in the Effective Fragment Potential Method. Mol. Phys. 2005, 103, 379-387. (57) Guidez, E. B.; Xu, P.; Gordon, M. S. Derivation and Implementation of the Gradient of the R–7 Dispersion Interaction in the Effective Fragment Potential Method. J. Phys. Chem. A 2016, 120, 639-647. (58) Xu, P.; Zahariev, F.; Gordon, M. S. The R–7 Dispersion Interaction in the General Effective Fragment Potential Method. J. Chem. Theory Comput. 2014, 10, 1576-1587. (59) Guidez, E. B.; Gordon, M. S. Dispersion Correction Derived from First Principles for Density Functional Theory and Hartree–Fock Theory. J. Phys. Chem. A 2015, 119, 2161-2168. (60) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347-1363. (61) Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory: Gamess a Decade Later; Elsevier: Amsterdam, 2005. (62) Boys, S. F. Construction of Some Molecular Orbitals to Be Approximately Invariant for Changes from One Molecule to Another. Rev. Mod. Phys. 1960, 32, 296-299. (63) 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. (64) Humphrey, W.; Dalke, A.; Schulten, K. Vmd: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33-38. (65) Su, P.; Li, H. Energy Decomposition Analysis of Covalent Bonds and Intermolecular Interactions. J. Chem. Phys. 2009, 131, 014102.

TOC graphic:


Page 28 of 29

Page 29 of 29

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



Dispersion Interactions in Water Clusters - The Journal of Physical

Recommend Documents