Ab initio water pair potential with flexible monomers.

Subscriber access provided by MAHIDOL UNIVERSITY (UniNet)

Article

Ab Initio Water Pair Potential with Flexible Monomers Piotr Jankowski, Garold Murdachaew, Robert Bukowski, Omololu Akin-ojo, Claude Leforestier, and Krzysztof Szalewicz J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 17 Feb 2015 Downloaded from http://pubs.acs.org on February 17, 2015

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 89

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

Ab Initio Water Pair Potential with Flexible Monomers Piotr Jankowski,∗,†,¶ Garold Murdachaew,†,§ Robert Bukowski,†,k Omololu Akin-Ojo,†,⊥ Claude Leforestier,‡ and Krzysztof Szalewicz∗,† Department of Physics and Astronomy, University of Delaware, Newark, DE 19716, USA, and Institut Charles Gerhardt (CTMM)-UMR 5253, CC 15.01, Université Montpellier II-CNRS 34095 Montpellier, Cedex 05, France E-mail: [email protected]; [email protected]

∗

To whom correspondence should be addressed Department of Physics and Astronomy, University of Delaware, Newark, DE 19716, USA ‡ Institut Charles Gerhardt (CTMM)-UMR 5253, CC 15.01, Université Montpellier II-CNRS 34095 Montpellier, Cedex 05, France ¶ Faculty of Chemistry, Nicolaus Copernicus University in Torun, Gagarina 7, 87-100 Torun, Poland (permanent address) § Laboratory of Physical Chemistry, Department of Chemistry, University of Helsinki, P.O. Box 55 (A.I. Virtasen aukio 1), 00014 Helsinki, Finland (current address) k Bioinformatics Facility, Institute of Biotechnology, Cornell University, Ithaca, NY 14853 USA (current address) ⊥ Department of Theoretical Physics, African University of Science & Technology, Galadimawa, Abuja, Nigeria (current address) †

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

Abstract A potential energy surface for the water dimer with explicit dependence on monomer coordinates is presented. The surface was fitted to a set of previously published interaction energies computed on a grid of over a quarter million points in the 12-dimensional configurational space using symmetry-adapted perturbation theory and coupled-cluster methods. The present fit removes small errors in published fits and its accuracy is critically evaluated. The minimum and saddle-point structures of the potential surface were found to be very close to predictions from direct ab initio optimizations. The computed second virial coefficients agreed well with experimental values. The effects of monomer flexibility in the virial coefficients were found to be much smaller than the quantum effects.


Page 2 of 89

Page 3 of 89

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


1

Introduction

The importance and ubiquity of water 1 has inspired decades of computational studies 2 of this substance. Interest is broad enough that a “biography of water” for the layperson has been published. 3 Most of the investigations of water concern water clusters and liquid water with monomers in their ground electronic state. In the Born-Oppenheimer approximation, properties of such systems are determined by their potential energy surfaces, i.e., by their total electronic energies which are function of the relative positions of all nuclei in the system. If the total energies of monomers are subtracted, one gets an intermolecular potential energy surface, which we will refer to as “potential”. For clusters larger than the dimer, the potential can be separated into the pairwise-additive component, i.e., the sum of pair potentials, and the pairwise-nonadditive component. The potentials can be developed by fits to experimental data, and are then called “empirical” potentials, or by fits to ab initio interaction energies computed on a grid of cluster geometries. The ab initio interaction energies can be computed using the supermolecular approach at various levels of including the electron correlation effects or by using symmetry-adapted perturbation theory (SAPT). 4–8 An excellent approximation in developing intermolecular potentials is to assume rigid monomers. However, some properties, such as the red shifts of the bonded OH stretches, vibrational predissociation, and intramolecular vibrational redistribution cannot be described in this approximation. Although rigid-monomer potentials predict most other properties of water rather well, the inclusion of monomer flexibility should make these predictions still better. Reviews of all types of water potentials can be found in refs 2,9,10. In the present paper, we will restrict our attention to first-principles pair potentials and ignore the nonadditive effects. So far, these effects have been investigated at the three-body level and the number of published three-body potentials is rather small. 11–20 Such investigations, as well as direct ab initio calculations 21–23 have found that the nonadditive effects contribute on the order of 20% to the interaction energies in clusters and in condensed phases of water. The higher than three-body effects can be approximated by the polarization model. 24,25 It was recently 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

shown that whereas the polarization model recovers only about half of the total nonadditive interaction energy of liquid water at ambient conditions, 26 it predicts the four-body contributions to within a few percent. 19

1.1

Water pair potentials with rigid monomers

Since pair (or two-body) potentials recover the majority of interaction energies in all forms of water, their accuracy is critical for obtaining correct predictions of properties of water. These potentials can be compared directly to results of spectroscopic measurements performed on the water dimer, 27–32 but also to the second virial coefficients extracted from measurements of properties of bulk water. 33–35 Empirical pair potentials fitted to bulk properties of water, such as SPC, 36 SPC/E, 37 or TIP4P, 38 are biased by their attempt to effectively incorporate the many-body effects present in the liquid and therefore cannot reproduce well the properties of dimers, 10,39 although polarizable water potentials fitted to bulk properties of water 25,40 perform better in this respect. Some more recent polarizable empirical potentials were fitted simultaneously to theoretical data on small water clusters 41,42 and may perform still better, but one cannot expect any empirical potentials fitted to bulk properties to be competitive to modern ab initio pair potentials in predicting properties of water dimer. On the other hand, empirical potentials fitted exclusively to water dimer spectra 29,43–45 predict such properties very well. The earliest ab initio pair potential for water was the MCY potential of ref 46 developed in 1976 which was fitted to interaction energies computed at 66 grid points. MCY was succeeded in 1990 by the NCC 11 two-body potential fitted to 350 interaction energies. The ASP ab initio potentials, 47,48 the first version developed in 1992, used a low-order version of SAPT 5,49,50 called intermolecular perturbation theory (IMPT) and were fitted to 350 IMPT interaction energies. In 1997, two potentials for the water dimer were developed by Mas et al., 39 dubbed SAPT-pp and SAPT-ss, using a more complete version of SAPT and 1003 grid points. These potentials were improved in 2000, leading to a potential labeled 4 ACS Paragon Plus Environment

Page 4 of 89

Page 5 of 89

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


SAPT-5s, by using 2510 grid points and a site-site form of the fitting function, with five symmetry-distinct sites per monomer. The SAPT-5s potential was used to predict spectra of the water dimer, 51–53 the second virial coefficients, 54 and in simulations of liquid water. 13 In 2002, Burnham and Xantheas developed the polarizable rigid-monomer TTM2-R 55 potential which was obtained from an earlier empirical potential by a partial reparametrization to a set of 25 ab initio dimer energies. The consecutive improvements of the SAPT-5s potential were the SDFT-5s potential of ref 56 based on a density-functional version of SAPT, 57 the CCpol-5s potential of refs 58–61 based on coupled-cluster (CC) calculations including single, double, and noniterative triple excitations [CCSD(T)], and CCpol-8s of ref 62 fitted to the same data as CCpol-5s with a more flexible functional form. All these potential used the same set of grid geometries.

1.2

Water dimer potentials with flexible monomers

Since many published intermolecular potentials are of the site-site form, i.e., are represented as sum of functions depending on site-site distances with sites usually placed on atoms, a flexible-monomer potential can be trivially obtained from a rigid-monomer potential simply by allowing the sites to move with the atoms as the monomers are allowed to distort. 63,64 The off-atomic sites within a molecule can be moved in proportion to the motions of the atoms. We will refer to such “atom-following” potentials as containing implicit monomer-flexibility effects. Note that this method is of zero cost beyond the rigid-monomer calculations. Examples of atom-following empirical effective pair potentials are the potentials SPC/F, 65 TIP4F, 66 MCYL, 67 and NCCvib 68 (flexible-monomer versions of the SPC, 36 TIP4P, 38 MCY, 46 and NCC 11,69 potentials, respectively). Somewhat surprisingly, the accuracy of the atom-following approach was apparently never checked until the work of refs 63 and 64, where it was found to perform poorly for Ar– HF already at H-F separations rHF corresponding to the turning points for the ground rovibrational state of HF (v = 0). The dependence of this performance on the configuration 5 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

of the dimer was inconsistent, with reasonable predictions for some configurations, but with the errors of the atom-following potential sometimes larger than those of the rigid-monomer potential for some other configurations. Only for separations between Ar and HF larger than that of the van der Waals (vdW) minimum, the atom-following approach did recover the correct trends consistently. For the water dimer, the predictions of the atom-following form of the SAPT-5s potential, dubbed SAPT-5s-af, at the vdW minimum were reasonable for up to 10% changes in bond lengths, but, for the distortions in the HOH angles, SAPT-5saf sometimes predicted changes in the wrong direction, even for very small distortions. The overall conclusion was that the atom-following approach should not be expected to represent monomer-flexibility effects sufficiently well. For empirical potentials, there exists an option of using a flexible-monomer form in fitting experiments. Apparently, this has not been done in fitting to bulk water properties, but was in fitting to dimer spectra. Reference 70 described the development of such a potential, using an MCY-like functional form, denoted as VRT(MCY-5f). This intermolecular potential was supplemented by the intramonomer potential of Polyansky et al. 71 One possibility of improving the atom-following approximation is to include explicit monomer-flexibility effects in the asymptotic part of the potential. Since asymptotic interactions depend exclusively on monomer properties, such calculations are of significantly reduced dimensionality compared to calculations for the dimer (e.g., three-dimensional for the water monomer compared to 12-dimensional for the water dimer). The use of flexiblemonomer asymptotic terms was usually combined with the atom-following method for the exponentially-decaying terms of the potentials. Such approach was used to develop the flexible TTM2-F water model 72 (later improved in refs 73 and 74). The performance of this strategy was investigated for Ar-HF in ref 63 by comparisons to full-dimensional calculations. It was found to work well in predicting the unexpanded dispersion energies (i.e., including the charge-overlap effects), but less well for the induction energies. Overall, the improvements were insufficient for this method to be used in a nonempirical way.


Page 6 of 89

Page 7 of 89

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


Since the atom-following plus flexible-monomer asymptotics approach was found not sufficiently adequate, the “k-flex” method was proposed in ref 64. This method estimates exchange energies from their approximate proportionality to the overlap of monomer densities (see, e.g., refs 75–77) which are are inexpensive to calculate. The k-flex method continues to use flexible-monomer asymptotic information. This method was tested on Ar–HF, 64 and was found to perform much better than the method with atom-following plus flexible-monomer asymptotics. It recovered over 60% of the v = 0 → 1 and v = 0 → 2 red shifts obtained from the “exact” flexible-monomer potential of ref 78. This relatively good performance of the k-flex potential can be compared to that of the atom-following plus flexible-monomer asymptotics potential, which recovered only about 30% of these shifts (unpublished results from ref 64). An application of this method to the water dimer is a part of the present paper. Obtaining a flexible-monomer water pair potential without any of the approximations described above is a formidable task since such potential is 12-dimensional. If just 3 points were used per dimension, more than 500,000 grid points would have to be used. Thus, the first ab initio flexible-monomer water dimer potentials of this type appeared only in 2006. 79,80 One of these surfaces, denoted as SAPT-5s′ f was fitted to SAPT interaction energies obtained using the same level of theory and the same basis set as in ref 54, but calculated on a set of 239,928 flexible-monomer grid points, 81 with intramonomer coordinates sampling the region near the classical turning points of monomer’s ground rovibrational state. The asymptotic part of the surface was constructed from ab initio computed monomer properties dependent on intramonomer coordinates, whereas in the short-range part, the linear parameters were made dependent on intramonomer coordinates and fitted to ab initio interaction energies. In order to improve the description of the bottom region of the potential, a variant of this surface, denoted as SAPT-5s′ fIR, was fitted using additional grid points as well as Hessians computed near the potential minimum. This surface was used to calculate the second virial coefficient (only briefly discussed in ref 79, but described in detail in ref 81) and



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 spectra of water dimer, including the spectral shifts of intramonomer frequencies. The shifts agreed better with experiment than any earlier theoretical predictions. More recently, the spectral calculations were repeated 82 with a hybrid potential, denoted by CCpol-8sf, being a combination of the CCpol-8s and SAPT-5s′ fIR potential. CCpol-8sf predicts water dimer spectra in a very good agreement with experiment, on par with that produced by other recent water potentials. 83–85 The other surface, HBB0, 80 was fitted to CCSD(T) dimer energies on a grid of almost 20,000 points. The HBB0 surface was followed in 2008 by HBB1, 86 obtained using an additional 10,000 grid points. Further refinement of this surface was done in 2009 resulting in the HBB2 surface. 83 Moreover, to improve the potential for large distances, this surface was spliced at large separations with the TTM3-F surface 73 (this method is denoted as HBB2/H). The most recent member of the HBB family of surfaces is the WHBB potential, 17 which is based on the same ab initio data and has the same form as HBB2, but the fitted quantity was the dimer interaction energy rather than the dimer total energy. The sequence of the HBB surfaces have shown increasing accuracy in predictions of various spectroscopic and thermodynamic properties of water. A potential very similar to HBB2/H, dubbed HBB2pol, was proposed by Medders et al. 18 This potential is based on the same ab initio points as HBB2/H and differs in the way the asymptotic component is turned on. More recently, Babin et al. published the MB-pol surface. 85 The ab initio interaction energies for abut 40,000 grid points were obtained at the CCSD(T) level in aug-cc-pVXZ, X=T, Q, basis sets 87,88 supplemented by a set of [3s3p2d1f ] midbond functions, followed by an extrapolation to the CBS limit. The long-range part was built using the TTM models 74,89 (electrostatics and induction) and the SAPT model 59 (dispersion). The short-range form of the fit is similar to HBB fits. The surface includes 1153 linear and 16 nonlinear parameters.


Page 8 of 89

Page 9 of 89

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


1.3

Plan of the paper

The development of the SAPT-5s′ f and SAPT-5s′ fIR surfaces was only briefly discussed in ref 79. One of the goals of the present work is to provide a more detailed description of the methodology used in ref 79. It was found in ref 82 that the subroutine computing the potentials of ref 79 contains a minor symmetry-related error in one of the terms. The error was fixed, but the coefficients of the fit were not reoptimized. We have performed such reoptimizations in the present work using the original set of interaction energies from ref 79. We have also assembled the corresponding hybrid potentials, following the prescription of ref 82. The potentials were used to compute characteristic points on the dimer surface, harmonic frequency shifts of intramonomer vibrations, and the second virial coefficients. Section 2 briefly describes some approximate approaches that we explored (but were not used for the development of the final potential). A detailed account of this work is included in the Supporting Information. Section 3 describes the selection of monomer and dimer grid points and the embedding of flexible monomers in the dimer. Section 4 discusses how the ab initio SAPT flexible-monomer interaction energies were fitted to obtain the SAPT-5s′ f and SAPT-5s′ fIR potentials. This section contains also a discussion of symmetries in the flexiblemonomer water dimer and the choice of symmetry-adapted basis functions needed for fitting. The SAPT-5s′ f and SAPT-5s′ fIR fits, as well as the hybrid surfaces obtained by combining such fits with the CCpol-8s rigid-monomer potential, are compared to ab initio interaction energies in Sec. 5. This section also compares the harmonic shifts of intramonomer vibrations upon dimerization with the results of ref 82. The potential is further analyzed in Sec. 6 at the minimum and saddle points. Methods for calculations of the second virial coefficients with inclusion of monomer-flexibility effects are discussed in Sec. 7. The coefficients computed from our potentials are compared to the experimental results and the magnitude of monomerflexibility effects is quantified. Finally, Sec. 8 summarizes our work. Supporting Information contains further details concerning the embedding procedures, the selection of grid points, the analysis of the surfaces, and the calculations of the second virial coefficient. A part of 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 89

the research described here (calculations of ab initio energies at all grid points, development of the symmetry-adapted form of the fit, and derivation of methods for calculations of virial coefficients) has been previously published in refs 81 and 79.

2

Approximate methods of including monomer-flexibility effects

Due to the k D scaling of calculations for a D-dimensional surface with k points per dimension, SAPT calculations at the same level as was used to develop the rigid-monomer SAPT-5s potential would take a few CPU-years for the 12-dimensional flexible-monomer water dimer with just 3 points per dimension. Accurate supermolecular methods such as many-body perturbation theory at the fourth order (MBPT4) or CCSD(T) would require comparable resources. Whereas such large-scale calculations are possible with parallel machines—and we did eventually perform our calculations for a quarter million points—we had first explored methods of approximate incorporation of monomer-flexibility effects developed earlier in refs 63 and 64, i.e., a combination of atom-following, explicitly-flexible asymptotics, and the k-flex methods. It turned out that for the water dimer these methods do not produce sufficiently accurate potentials to be confronted with spectroscopic measurements. Still, such approaches can be useful for larger dimers with lower accuracy goals for which brute force ab initio calculations are not possible. Therefore, we present the results of such investigations in Sec. SI-II of the Supporting Information.

3 3.1

Selection of grid points Dimer grid

We have decided to construct the flexible-monomer potential as an extension of the rigidmonomer SAPT-5s potential, therefore, we have used the same 2510 dimer geometries, i.e.,


Page 11 of 89

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


sets of intermonomer coordinates. The dimer’s geometry is identified by the distance R between the monomers’ centers of mass (COMs) and the five Euler angles defining the relative orientations of monomers A and B: ωA = (αA = 0◦ , βA , γA ) and ωB = (αB , βB , γB ), see Fig. 1 of ref 39. The starting orientation (i.e., R and all Euler angles equal to zero) of each water molecule is such that all three atoms are located in the xz plane, the molecular COM is at the origin of the coordinate system, and the zˆ coordinate of O is positive. For a symmetric water molecule (i.e., equal OH bond lengths), the bisector of the HOH angle is ˆ For an asymmetrically distorted water molecule, there are several possible aligned along z. choices of the initial monomer orientation in the xz plane, as discussed below. After this orientation is chosen, Euler angle rotations (in the zyz convention, 90 see also Fig. 1 of ref 39) ωA and ωB are applied to monomers A and B, respectively. Finally, monomer B is ˆ We used (slightly rounded) CRC masses 34 of the translated by a distance R along +z. most common isotopes: m16 O = 15.994915 amu and m1 H = 1.007825 amu, while ref 54 used m16 O = 16 amu and m1 H = 1 amu, to determine the position of the COM. This resulted in a slight shift of monomer Cartesian site coordinates when monomers are placed as defined above. However, monomer’s geometry remains the same (including the relative positions of the off-atomic sites). Since SAPT-5s is a site-site fit, this mass change has therefore no effect on the actual potential.

3.2

Monomer grid

For each dimer grid point, we selected a number of associated flexible-monomer geometries. We restricted the range of variability of the intramonomer coordinates to that in the ground vibrational state of water and selected the monomer geometries which correspond to the classical turning points of the monomer intramolecular potential in this state. We used the symmetry coordinates, s1 , s2 , and s3 : 91 the symmetric stretch, the bend, and the asymmetric



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 89

stretch, respectively, defined as s1 =

√ ∆r1 + ∆r2 ∆r1 − ∆r2 √ √ , s2 = r1 r2 ∆θ1 , s3 = , 2 2

(1)

where r1 , r2 , and θ are bond lengths, and the bond angle, respectively, and ∆r1 = r1 − hrOH i0 , ∆r2 = r2 − hrOH i0 , ∆θ1 = θ1 − hθHOH i0 .

(2)

The reference “hri0 ” monomer geometry, i.e., the geometry averaged over water monomer ground-state vibrational motion, (hrOH i0 = 0.9716257 ˚ A, hθHOH i0 = 104.69◦ ) was the same as used in refs 92 and 39. These values, derived from experiments, are somewhat different from those predicted by the PJT2 potential: 71 hrOH i0,PJT2 = 0.9754 ˚ A, hθHOH i0,PJT2 = 104.44◦ . Notice that the standard definition of si assumes the equilibrium coordinates, re , as the reference geometry. We will use the notations sA = (s1 , s2 , s3 ), rA = (r1 , r2 , θ1 ), and, for monomer B, sB = (s4 , s5 , s6 ), rB = (r3 , r4 , θ2 ). Also, s = (sA , sB ). To identify the turning points in the ground-state vibration, we have chosen the accurate empirical PJT2 71 potential. Since the publication of PJT2, a few even more accurate potentials have been developed, 93–95 but these improvements are insignificant for our purposes. All these potentials predict spectra of water monomer extremely well. 96 The PJT2 potential A, 104.500◦ ). The energies will be measured has the minimum geometry (re , θe ) = (0.95792 ˚ relative to this minimum. The zero-point energy (ZPE) is 4634.76 cm−1 (13.2515 kcal/mol), whereas energies of lowest excited states: the bend, symmetric stretch, and asymmetric stretch are 1594.68, 3657.13, and 3755.84 cm−1 (4.5529, 10.4563, and 10.7385 kcal/mol), respectively, above the ground state. The barrier to linearity is 31.35 kcal/mol. All these values were computed by us with complete inclusion of anharmonic effects and the excitations energies agree to within 0.02 cm−1 with calculations of ref 71. The classical turning points in the coordinates of the ground vibrational state were determined by identifying the energy contour on the PJT2 potential surface with the value of the ZPE energy. The 14


Page 13 of 89

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


selected distorted geometries all have energies within 1.3 cm−1 of ZPE. We included all 6 points with only one si distorted, as well as combinations of simultaneous distortions of s1 and s2 (with s3 = 0) and of s2 and s3 (with s1 = 0). For the s1 –s2 (s2 –s3 ) combination, we have arbitrarily chosen r1 = r2 (θ) and adjusted θ (s3 ). Table SI-I in the Supporting Information shows the geometries and energies of these 14 monomers, as well as of the re and hri0 monomers.

3.3

Embedding

Each of the intermolecular potentials developed in the present paper can be expressed entirely via distances between sites on different monomers. However, the most popular way to define geometry of a van der Waals complex is to use R, ωA , ωB , rA , and rB . To uniquely determine the Cartesian coordinates of the atoms, one has to specify how to “place” the molecule into “the frame of the dimer”. More precisely, one has to choose a monomerspecific coordinate system such that when its axes are parallel to the space-fixed system, the Euler angles describing this monomer are zero. This procedure is called an embedding and various possible choices for embedding and their implementation are discussed in details in Sec. SI-I of the Supporting Information. In the present work, we need embedding procedures at three stages. First, we have to produce the 12-dimensional grid for ab initio calculations. To combine the dimer and monomer grids, we have chosen the principal axis embedding (PAE). Note that with our form of the fitting function, the potential fitted to interaction energies computed on this grid does not depend on any embedding [except indirectly due to the choice of grid points (see discussion in Sec. 4)]. Therefore, when the potential was averaged over the monomer coordinates in the calculations of the virial coefficient, we could use a different embedding and the Eckart-axes embedding (EAE) 97–99 was employed. Finally, an embedding has to be chosen to define the hybrid interaction energy surface, see Sec. 4.4.



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

3.4

Page 14 of 89

Total grid

Our original intent was to calculate the interaction energies using only the k-flex method on the direct product of the dimer grid from ref 54 and the distortion grids for both monomers (monomers with indices 0–14 in Table SI-I of the Supporting Information. This gives 15 × 15 × 2510 = 564,750 grid points. When we decided to compute instead the SAPT interaction energies, due to computational demands we had to reduce this number in half. The reduction procedure, described in details in Sec. SI-III of the Supporting Information, led to the set of 239,928 grid points used in fitting of the SAPT-5s′ f potential. This number includes 2,660 grid points with fixed ωA and ωB and varying R computed to enable radial plots described in Sec. SI-V of the Supporting Information.

4 4.1

Fitting procedure Improvement of the rigid-monomer fit

We have modified the SAPT-5s potential by including in the fitting function a term modeling the induction energy via polarization of monomers. The modified potential, denoted as SAPT-5s′ , is in the form of the sum over sites (atomic and off-atomic) in both monomers and depends on intersite distances Rab : V

X

=

ab

g (Rab )e

αab −βab Rab

a∈A, b∈B

−

+

q a qb f1 (δ1ab , Rab ) Rab

+

1 B A2 α ¯ |F | + α ¯ A |F B |2 f6 (δ6d , R), 2

where g ab (Rab ) =

P3

n=0

Cnab ab fn (δn , Rab ) n Rab n=6,8,10 X

n ab ab aab n Rab (a0 = 1 kcal/mol), fn (δn , Rab ) = 1 −

(3)

Pn

i=0

ab R )i (δn ab i!

are Tang

and Toennies 100 damping functions, and the dipolar field F A of monomer A at COM of monomer B is FA =

ˆ R ˆ · µA ) − µA 3R( , R3


(4)

Page 15 of 89

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


with µA denoting the (permanent) dipole moment of monomer A. All the remaining symbols are parameters which will be discussed below. Only the last term in Eq. (3) is new, all the remaining terms are the same as in the SAPT-5s potential. We will use V to refer to a fit and Eint to refer to an ab initio computed interaction energy. Note that all interaction energies discussed in this paper will be the “vertical” interaction energies, i.e., the differences between the total dimer energies and the sum of monomer energies at the monomer configurations in the dimer. The first term in Eq. (3) models the exchange interactions and a part of overlap effects from other components; the second term models the damped long-range electrostatic energy with qx denoting partial charges; the third term represents the damped long-range dispersion energies, with Cnab denoting site-site van der Waals constants; the fourth term models the induction energy with α ¯ X denoting the static dipole-dipole isotropic polarizability of monomer X, α ¯ = (αxx + αyy + αzz )/3. The parameters qx , Cnab , and α ¯ X were obtained using ab d only the monomer information. The parameters αab , βab , aab n , δn , and δ6 were obtained in a

global fit to 2510 ab initio interaction energies from ref 54. The name of the potential, SAPT-5s, refers to its 5 symmetry-distinct site types. In addition to the sites on atoms, O, H1 , and H2 , there are also off-atomic site types dubbed D1, D2, and D3, all lying on the σv plane of the monomer, see Fig. 2 and Table I in ref 54. Due to our use of more accurate atomic masses, as discussed earlier, all z coordinates in this table are incremented by 0.000 901 563 6 bohr. Not all sites participate in all types of interactions, as indicated in Table I of ref 54. The COM dipole-dipole polarizabilities were computed at the coupled Hartree-Fock plus MBPT2 level of theory 101,102 using the monomer-centered part of the basis set applied in SAPT calculations. The polarization center was placed at the COM P of each monomer. The dipole moments were computed as µA = a∈A qa ra , where the ra are the site positions in monomer A relative to COM. Therefore, the distance used in Eqs. (3)

and (4) should be R. In practice, it was convenient to use ROO rather than R. The partial charges were taken from ref 54. The Cnab coefficients had to be refitted since these terms now



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 89

model almost exclusively the dispersion energy while the main part of the induction energy is modeled by the explicit induction term. We will refer to Cnab ’s as the dispersion coefficients (despite the fact that they also model a small amount of residual induction energy; vide infra in Section 4.2.3). We have used the COM induction and dispersion coefficients from ref 54 to compute the asymptotic induction and dispersion energies at the 2510 points from the regular grid, but with all R’s increased by 3 ˚ A. Then the site-site dispersion coefficients Cnab B A2 P P n − 12 α were found fitting these values by a∈A, b∈B n=6,8,10 Cnab /Rab ¯ |F | + α ¯ A |F B |2 and

keeping all parameters in the polarization term fixed. All the remaining fitting procedures were the same as in ref 54. As shown in Table 1, separating out the dominant portion of the asymptotic induction energy improved the long-range fit: the root-mean square error (rmse) is 65% of the original rmse, although in the whole range it is only marginally different from and not better than that of the original SAPT-5s fit. We will use the SAPT-5s′ potential as the “parent” rigid-

monomer potential for SAPT-5s′ f. However, when we compare to a rigid-monomer potential, it will be to the original rigid-monomer potential SAPT-5s, as published in refs 54 and 51.

4.2

SAPT-5s′ f flexible-monomer fit

Our flexible-monomer water dimer potential was obtained by directly fitting the quarter million ab initio points with a functional form explicitly dependent on intramonomer coordinates. This form was an extension of the rigid-monomer potential of Eq. (3), obtained by making all linear parameters, including the charges and the polarizability, polynomials in s. We have not expanded any of the nonlinear parameters since nonlinear least-square fits would have been a too difficult task given the large number of fitted data points. Thus, the general form of the fit is V

=

X

a∈A, b∈B

−

ab

g (Rab , s)e

αab −βab Rab

+

qa (sA )qb (sB ) f1 (δ1ab , Rab ) Rab

1 B α ¯ (sB )|F A (sA )|2 + α ¯ A (sA )|F B (sB )|2 f6 (δ6d , R). 2 16

ACS Paragon Plus Environment

C ab (sA , sB ) + fn (δnab , Rab ) n n Rab n=6,8,10 X

(5)

Page 17 of 89

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 zeroth-order term of each polynomial was kept constant, equal to the corresponding rigid-monomer potential value, i.e., the flexible-monomer potential becomes equal to the rigid-monomer SAPT-5s′ potential at s = 0. For deformed monomers, the off-atomic sites were kept at same distance from the molecular plane as in the rigid-monomer case, whereas their projections on the plane were on the HOH bisector and the distance between the projection and the O atom multiplied by (r1 + r2 ) cos(θ/2) . 2hri0 cos(hθi0 /2)

(6)

Since the distances Rab change as the monomers are deformed, the atom-following property is built into our fit as the zeroth-order description of the monomer flexibility effects. Symmetries of the fit and the choice of symmetry-adapted expansions in s are discussed in Sec. 4.5. There are too many parameters in the potential to list them in tables. These parameters are available within the FORTRAN potential subroutine in the Supporting Information. We have followed a four-step fitting procedure in developing the flexible-monomer potential, which is an extension of the rigid-monomer approach described in ref 54 and in Sec. 4.1. (1) Fit of qa (sA ) to ab initio computed COM multipole moments for 15 monomer geometries. (2) Fit of α ¯ A (sA ) to ab initio computed COM polarizabilities for 15 monomer geometries. (3) Fit of Cnab (sA , sB ) to the asymptotic expansion of the dispersion (and residual induction) energies obtained on a shifted flexible-monomer grid using ab initio COM van der Waals coefficients computed for all 225 combinations of monomer coordinates. (4) Fit of short-range site-site parameters to the SAPT ab initio interaction energies on the grid of points. Note that the asymptotic coefficients were computed at the level of theory consistent with the SAPT calculations. 39 4.2.1

Site charges

The first stage of our procedure was to fit the site charges qa , with a running only over the atomic and D1 sites, to reproduce the ab initio-computed tesseral components of the 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 89

multipole moments: Qlk , l = 0–8. As in ref 54, we included the monomer charge Q00 = 0 to impose charge neutrality. The number of nonzero Qlk components plus Q00 is 25 for symmetric monomers and 45 for the nonsymmetric ones. Our set of 14 flexible-monomer geometries includes 8 of the former and 6 of the latter type. Therefore, the total number of fitted data points is 470. The charge neutrality condition was imposed with a weighting factor of w = 105 , large enough so that the residual charge was of order of 10−14 atomic units for any monomer geometry. As in ref 54, the weighting factor for all other tesseral components was w = 10−3l−k . The partial charge on each of the involved sites was expanded as qa (s1 , s2 , s3 ) = qa(0) + qa(1) s1 + qa(2) s2 + qa(3) s3 + qa(4) s1 s2 + qa(5) s2 s3 + qa(6) s21 + qa(7) s22 + qa(8) s23 . (7) The term s1 s3 , describing the coupling between the two stretches, was neglected, consistent (i)

(i)

with our choice of the monomer grid. Due to symmetry, qD11 = qD12 for all i. For s3 = 0, the (i)

(i)

two H sites are equivalent, which implies qH1 = qH2 for all i except 3, 5, and 8. The planar symmetry also requires that the charges of O and D1 do not depend on the sign of s3 , i.e., (3)

(5)

qa (s3 ) = qa (−s3 ), for a = O, D11 , and D12 . This implies that qa = qa = 0 for these sites. For the pair of H sites, the charge on H1 at s3 should be the same as the charge on H2 at (i)

(i)

(0)

(2)

(3)

(5)

(7)

(8)

(0)

(2)

(3)

(5)

(7)

(8)

−s3 . Using the relationships between qH1 and qH2 derived so far, we can write qH1 (0, s2 , s3 ) = qH + qH1 s2 + qH1 s3 + qH1 s2 s3 + qH1 s22 + qH1 s23 , qH2 (0, s2 , −s3 ) = qH + qH1 s2 − qH2 s3 − qH2 s2 s3 + qH1 s22 + qH1 s23 . (3)

(3)

(5)

(5)

(8)

(8)

(8)

Therefore, qH2 = −qH1 , qH2 = −qH1 , and qH2 = −qH1 , which brings the total number of (i)

adjustable qa parameters to 20. As the rmse’s listed in Table 1 demonstrate, the resulting fit for the dipole moment is excellent. The l = 2 and l = 3 tesseral components were reproduced almost as well as in the rigid-monomer case in ref 54. The rmse increases with the rank of the moment, but the 18 ACS Paragon Plus Environment

Page 19 of 89

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


moments also become larger and thus the relative error remains small. Of course, higher moments contribute less than the lower ones, a fact that is reflected by the weighting factors used in the fit. A plot illustrating the quality of the fit for l = 1 can be found in Fig. SI-4 of the Supporting Information. 4.2.2

Dipole-dipole polarizability

In the next stage, we fitted the isotropic dipole-dipole polarizability α ¯ appearing in Eq. (3) to flexible-monomer ab initio polarizability values. The model was similar to that used for charges: α ¯ (s1 , s2 , s3 ) = α ¯ (0) + α ¯ (1) s1 + α ¯ (2) s2 + α ¯ (3) s3 + α ¯ (4) s1 s2 + α ¯ (5) s2 s3 + α ¯ (6) s21 + α ¯ (7) s22 + α ¯ (8) s23 . (9) Since due to our replacement of R by ROO in Eq. (3) α ¯ is placed on O atom, the symmetry arguments are the same as for the charged O site. Thus, α ¯ (3) = α ¯ (5) = 0, and this condition and the fact that α ¯ (0) was kept constant at the rigid-monomer ab initio value reduces the number of independent expansion coefficients in Eq. (9) from 9 to 6. The rmse given in Table 1 is less than 0.1% of the rigid-monomer value of α ¯ (0) , see also Fig. SI-4 of the Supporting Information. 4.2.3

Long-range dispersion coefficients

The flexible-monomer fit of the site-site asymptotic coefficients Cnab , with indices a and b running only over atomic sites, was based on ab initio flexible-monomer COM coeffi{Λ}

cients Cn , where Λ is the set of “quantum numbers”, Λ = {LA , KA , LB , KB , L}, labeling the dimer angular functions. 6 We added the induction plus dispersion coefficients, {Λ}

Cn

{Λ}

{Λ}

= Cn,ind + Cn,disp computed for all 225 products of monomer geometries and then

used these coefficients to generate a set of asymptotic induction plus dispersion energies on the original 12D grid with R shifted outwards by 3 ˚ A. This procedure gave a total of 564,750 asymptotic induction plus dispersion energies. This set was fitted by the expression 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

P P ab

n

Page 20 of 89

B n Cnab (sA , sB )/Rab − α ¯ (sB )|F A (sA )|2 + α ¯ A (sA )|F B (sB )|2 /2, with all parameters

in the polarization term constant. Each coefficient Cnab was assumed to be a polynomial in intramonomer variables sA , sB : Cnab (sA , sB ) =

X

Cnab,p

6 Y

spi i ,

(10)

i=1

p

where p = (p1 , p2 , . . . , p6 ) is a vector of integers representing the powers of internal coordinates and the summation runs over the chosen vectors. Among the powers of si included in the expansion are all linear terms, mixed quadratic terms such as s1 s4 , s2 s5 , s3 s6 , as well as various other terms, for example s23 , s26 , s23 s6 , and s3 s26 . To utilize the symmetries of the system, we wrote Eq. (10) in terms of symmetry-adapted functions, as described in Sec. 4.5. The total number of free parameters in Eq. (10) was 63. The fit was unweighted. A plot illustrating the quality of the long-range induction plus dispersion energies is given in Fig. SI-5 of the Supporting Information. The rmse of the fit on the flexible-monomer points is 0.0015 kcal/mol, see Table 1. 4.2.4

Fit of short-range parameters

The expansion of the site-site coefficients aab n was analogous to that used for the site-site asymptotic coefficients in Eq. (10): aab n (sA , sB )

=

X

aab,p n

p

6 Y

spi i ,

(11)

i=1

and the actual fit was performed in terms of symmetry-adapted functions discussed in Section 4.5. The number of vectors p in the expansion of each aab n was somewhat larger than in the case of Cnab and included also the following powers: s2i , s1 s2 , s4 s5 , s2 s3 , s5 s6 , s2 s23 , and s5 s26 . Note that the coefficient aab 0 , which was fixed at 1 in the rigid-monomer case, was ab,(0,0,0,0,0,0)

expanded according to Eq. (11), with an

= 1. This composition was established

after trying several choices, and selecting one that gave satisfactory performance in terms of 20 ACS Paragon Plus Environment

Page 21 of 89

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


rmse without introducing artifacts like an oscillatory behavior in the resultant fit. In fitting the 568 linear short-range parameters, the asymptotic parameters were held constant, as were the nonlinear parameters taken from the rigid-monomer SAPT-5s′ fit. The weighting factor was the same as in ref 54, 1/(Eint + E0 )2 , except that the parameter E0 was changed to 7.3 kcal/mol since the most negative energy in the set of 239,928 of interaction energies was about -6.4 kcal/mol. The linear least-squares problem turned out to be a difficult task due to its dimension of 239,928 × 568 and instabilities caused by the fact that the values of the fitted parameters were changing substantially with minor variations of the fitting procedures. To remove linear dependencies, the equations were solved using the singular value decomposition (SVD) method. 103 The list of the parameters used in all terms of the SAPT-5s′ f fit is given in Table 2. The symmetry adapted functions are defined in Section 4.5.

4.3

SAPT-5s′ fIR flexible-monomer fit

In addition to the SAPT-5s′ f surface described in Sec. 4.2, the authors of ref 79 developed also the surface denoted as SAPT-5s′ fIR. The reason for the creation of this surface was the finding that SAPT-5s′ f is not accurate enough near the minimum, as indicated by the discrepancies between the harmonic frequencies of the fit and those computed ab initio at the same level of theory as used in the development of SAPT-5s′ f. 79 To overcome this problem, an additional 2,896 points were computed 79 in the vicinity of the minimum and used in fitting with weights 100 (not 1000, as incorrectly stated in ref 79). Furthermore, the set of fitted data also included 171 second derivatives of the SAPT energies with respect to atomic Cartesian coordinates computed numerically at the minimum coordinates. For these data, the weights in the least-squares procedure were equal to 10. As in the case of SAPT-5s′ f, we recomputed SAPT-5s′ fIR after fixing the symmetry error discussed in Sec. 1.3. The effects of this error on the 2006 and 2012 potentials are analyzed in Sec. SI-VII of the Supporting Information. For all interaction energies below 20 kcal/mol, the discrepancies are smaller 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

than 1 kcal/mol, and only a very small fraction of points differs by more than 0.2 kcal/mol. The correction of the error in the code plus the reoptimization of the fit resulted in larger changes that the correction alone.

4.4

Hybrid potential energy surfaces

The CCpol-8s rigid-monomer surface developed in ref 62 performed very well in predictions of the spectra of water dimer, both due to the relatively high level of theory used in ab initio calculations and due to an accurate fit. Reference 82 proposed a hybrid surface which combines the high accuracy of CCpol-8s with the monomer-flexibility property of SAPT5s′ fIR. This surface, dubbed CCpol-8sf, is defined by the formula VCCpol-8sf (R, ωA , ωB , rA , rB ) = VCCpol-8s (R, ωA , ωB ) + [VSAPT-5s′ fIR (R, ωA , ωB , rA , rB ) − VSAPT-5s′ fIR (R, ωA , ωB , hrA i0 , hrB i0 )] , (12) so that the potential at a geometry with flexible monomers (i.e., in general distorted from their hri0 values) is obtained as the sum of the rigid-monomer CCpol-8s potential (which was developed with monomers at the hri0 geometry) plus a flexibility correction computed using the SAPT-5s′ fIR potential. Note that the definition (12) requires assumption of an embedding. A similar hybrid potential for the dimer of hydrogen molecules has been recently developed in ref 104. Reference 82 used the so-called Radau f = 1 embedding (see Sec. SI-I of the Supporting Information and refs 97 and 98). We included additionally the option of the Eckart embedding in the present version of the hybrid surface. The embeddings are programmed internally in our hybrid potential code, whereas calls to this code use exclusively the Cartesian coordinates of atoms. Thus, if the calling program works in terms of R, ωA , ωB , rA , and rB , these coordinates should be used to compute the Cartesian coordinates of atoms before calling the potential subroutine. The embedding used in the calling program can be


Page 22 of 89

Page 23 of 89

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


different from the one used in the definition of the hybrid potential. The procedure defined by Eq. (12) can be used to produce similar potentials based on any combination of rigid and flexible-monomer surfaces. For instance, we have obtained a hybrid potential from the CCpol-8s and SAPT-5s′ f potentials. To avoid confusion, we decided to rename the surface CCpol-8sf from ref 82 to CCpol-8sfIR, while the new surface, created from SAPT-5s′ f, will be called CCpol-8sf. In this way, we obtain a direct correlation between the names of the component surfaces and the hybrid ones. The fits obtained in the present work will be distinguished from literature fits denoted by the same acronyms by appending the year the potential was developed to the symbol. Thus, the potential used in ref 82 is now denoted as CCpol-8sfIR[2012] and the potentials obtained in the present work are distinguished by “[2014]”. The list of potentials developed and discussed in the present work is given in Table 3. All the results presented in the tables and figures were obtained with the [2014] version of the potentials, unless noted otherwise.

4.5

Symmetry-adapted flexible-monomer water dimer potential

In this section, we will consider the symmetry properties of the flexible-monomer water potential and introduce an analytic representation of this potential in terms of symmetryadapted functions, as used in the SAPT-5s′ f fit. This method was initially developed in ref 81. The simple cases of the terms that depend only on an isolated-monomer symmetry, i.e., the charges qa and the polarizability α ¯ , were already discussed in Sec. 4. The sum over sites a ∈ A, b ∈ B involving the exponential and dispersion contributions in Eq. (5) can be arranged into a sum involving types of atoms: Vexp/disp =

X

V IJ ,

(13)

I≤J

where the indices I and J go over the ordered list of the site types {O, D3, D2, H} for the exchange part and {O, H} for the dispersion part (the site D1 does not contribute to either 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 89

component) and are not specific for monomers A and B (i.e., there is only one term V OH and so on). Since the potential V must be invariant under the operations of the dimer symmetry group G16 , each of the terms V IJ must satisfy 1 A B (1 + PAB )(1 + P12 )(1 + P12 )(1 + E)V IJ = V IJ , 16

(14)

where PAB denotes the permutation of the molecules A and B (i.e., it interchanges all indices X labeled by A or belonging to A with those labeled by B or belonging to B in Eq. (5)), P12

stands for permutation of hydrogen atoms of monomer X, and E is the inversion operator. Each term V IJ can be written as a linear combination of terms of the form sk11 sk22 · · · sk66 w(β(sA , sB ); Rab ),

(15)

where w(β; R) = Rn exp(−βR) for the exponential component and w(β; R) = fn (β, R)/Rn for the dispersion component. In general, individual terms of the form (15) are not invariant under the symmetry operations of the G16 group. Thus, the problem is to find combinations of such terms which are symmetry-adapted according to Eq. (14). Such symmetry-adapted combinations, later on referred to as ψiIJ , are obtained by applying the G16 symmetrizer from Eq. (14) to terms of the form (15) and examining the structure of the resulting expressions. The index i in functions ψiIJ labels consecutive terms with different sets of powers km . Labels indicating the types (exponential or dispersion) of these functions and (inverse) powers n of site-site distances will be dropped since the symmetry properties do not depend on these features. However, one should remember that the actual formulas do contain summations over type- and n-specific symmetry-adapted functions. For instance, the dispersion term OO depending on the OO distance, Vdisp , can be written as

OO Vdisp =

XX n

OO,n OO,n Cdisp,i ψdisp,i ,

i


(16)

Page 25 of 89

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


OO,n are the fitting coefficients to be found. where Cdisp,i

The symmetry-adapted functions ψiIJ used in our fit are listed in Table 4. To derive these functions, it is useful to first determine basic symmetry properties of individual sitesite distances. First, we note that any interatomic distance is invariant under inversion, e.g., ERH1A OB = E |RH1A − ROB | = |−RH1A + ROB | = RH1A OB . This invariance also applies to any intramonomer interatomic distances and consequently to the si functions of Eq. (1) since such functions depend only on interatomic distances and on the angle θ which is given by the scalar products of interatomic vectors. In the SAPT-5s′ f fit, the nonlinear parameters β were assumed to be constants. However, as indicated by Eq. (15), we will consider here a more general case where these parameters are functions of (sA , sB ). To make final formulas simpler, we will constrain the function β to be expanded in terms of polynomials invariant under PAB and containing only even powers of s3 and s6 , e.g., the functions h+ i defined in captions of Table 4. With such constraints, the functions β(sA , sB ) are invariant under all symmetry operations and the symmetry-adapted functions ψiIJ have the same general form as the ones derived assuming constant nonlinear parameters. In some cases, the functions w will be explicitly symmetrized with respect to PAB and in such cases β(sA , sB ) do not need to be invariant under PAB . The distances involving off-atomic sites D2 and D3 can be written as functions of atomic coordinates. The position of a D3 site in molecule A can be expressed as (for simplicity, we drop the “A” subscript from atom labels) RD3 = RO + κ3

ROH1 ROH2 + ROH1 ROH2

,

(17)

where κ3 , proportional to the quantity from Eq. (6), depends symmetrically on the two OH vectors and on the HOH angle. It follows that ERD3 = −RD3 and thus any distance between D3 and an atom or the D3 site on another molecule is invariant under inversion. From Eq. (17), it is also clear that the position of D3 in a monomer (and thus, also the distance from this site to any other site on the other molecule) is invariant under permutation 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 89

of hydrogen atoms within this monomer. The positions of the two off-plane sites D2 in monomer A can be expressed as (here again we drop subscript “A” from all atom labels) RD21/2 = RO + κ2

ROH1 ROH2 + ROH1 ROH2

±σ

ROH1 × ROH2 , |ROH1 × ROH2 |

(18)

where κ2 is the same (up to a proportionality factor) as κ3 , σ is a constant, and + and − A refer to D21 and D22 , respectively. From this definition, it follows that P12 RD21 = RD22

and ERD21 = −RD22 . Therefore, if XB is an atom or a D3 site in monomer B, then A P12 RD21A XB = RD22A XB and ERD21A XB = RD22A XB . Note also that ERD21A D21B = RD22A D22B

and ERD21A D22B = RD22A D21B . A B The P12 and P12 operations change the sign of s3 and s6 , whereas all the remaining si ’s

are invariant under these operations. Since all si ’s are invariant under E, the action of the symmetrizer on any function of the form (15) can be written as A B (1 + PAB ) sk11 sk22 · · · sk66 (1 + (−1)k3 P12 )(1 + E)w(β; Rab ) . )(1 + (−1)k6 P12 4.5.1

(19)

O-O and D3-D3 terms

The term V OO depends on the distance ROA OB and (sA , sB ), all invariant under inversion, so that the operator E in Eq. (19) can be dropped. The distance ROA OB is also A B invariant upon operations PAB , P12 , and P12 . Therefore, the result of symmetrization of

sk11 sk22 · · · sk66 w(β OO ; ROA OB ) becomes (apart from an overall proportionality constant) (1 + (−1)k3 )(1 + (−1)k6 )w(β OO ; ROA OB )(1 + PAB )[sk11 sk22 · · · sk66 ]

(20)

so that the factor containing si can be expressed in terms of polynomials symmetric in PAB and containing only even powers of s3 and s6 , such as the functions h+ i which were actually used in the SAPT-5s′ f fit. A generalization to higher-order polynomials is straightforward.


Page 27 of 89

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 V D3D3 term can be obtained directly from V OO by replacing in the formulas all pairs of indices OO by the pairs D3D3. 4.5.2

O-D3 terms

For V OD3 , the symmetry properties are analogous to the OO case except that w is not symmetric under exchange of monomers. Therefore, general symmetry adapted functions are (1 + PAB )[sk11 sk22 · · · sk66 w(β OD3 ; ROA D3B )]

(21)

assuming even k3 and k6 . In practice, it is convenient to replace this set of functions by the − ones listed in Table 4, expressed through (anti)symmetrized polynomials h+ i and hi .

4.5.3

O-H and D3-H terms

Since the expression sk11 sk22 · · · sk66 w(β OH ; ROA H1B ) is invariant under inversion, the operator (1 + E) can be omitted from the G16 symmetrizer. In contrast to the cases considered above, B the symmetry operation P12 , replacing the label “1” of a hydrogen in monomer B with “2”

and vice versa, affects also intermolecular variables, so that we get B (1 + PAB )[sk11 sk22 · · · sk66 (1 + (−1)k6 P12 )w(β OH ; ROA H1B )],

(22)

where k3 must be even. The specific OH terms used in the fit (some of them expressed in terms of functions h± i ) are presented in Table 4. Notice that this is one of cases mentioned above when β does not have to be symmetric in PAB . To get the symmetry-adapted basis functions for the V D3H term, it is sufficient to replace label O with D3 in the expression above.



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

4.5.4

Page 28 of 89

H-H terms

Applying the G16 symmetrizer to sk11 sk22 · · · sk66 w(β HH ; RH1A H1B ) and noting that this term is invariant under inversion, we obtain A B (1 + PAB ) sk11 sk22 · · · sk66 [(1 + (−1)k3 P12 )(1 + (−1)k6 P12 )w(β HH ; RH1A H1B )] .

(23)

When k3 and k6 are of the same parity, the function in square brackets is symmetric in PAB and thus only symmetrized polynomials in s will contribute. Otherwise both symmetric and antisymmetric components will be present. Functions of this form that were used in the actual fit are listed in Table 4. 4.5.5

O-D2 and D3-D2 terms

Applying the G16 symmetrizer to an OD2 term sk11 sk22 · · · sk66 w(β OD2 ; ROA D21B ), we notice that since the effect of the inversion operator E on w is to change the index “1” on monomer B B into “2”, the function (1 + E)w(β OD2 ; ROA D21B ) = (1 + P12 )w(β OD2 ; ROA D21B ) is already B symmetric with respect to P12 . Thus, only even powers of s6 can contribute and the operator A (1 + E) can be dropped. Since w is invariant under P12 , also k3 has to be even. Thus, the

symmetry-adapted components are B (1 + PAB )[sk11 sk22 · · · sk66 (1 + P12 )w(β OD2 ; ROA D21B )],

(24)

or with the help of functions h± i , as in Table 4. Symmetry adapted functions representing the D3D2 terms are obtained by replacing O by D3 in the OD2 formulas.


Page 29 of 89

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


4.5.6

D2-D2 terms

Since the function f = (1 + E)w(β D2D2 ; RD21A D21B ) = w(β D2D2 ; RD21A D21B ) + w(β D2D2 ; RD22A D22B ), A B and similarly its counterpart obtained from w(β D2D2 ; RD21A D22B ), have the property P12 P12 f = A B f and P12 f = P12 f , k3 and k6 must be of the same parity in Eq. (19). On the other hand, A B any D2D2 primitive function acted upon by (1 ± P12 )(1 ± P12 ) is invariant under E. Thus,

the operator (1 + E) can be dropped in Eq. (19), leading to a symmetric function of the form A B (1 + PAB )[sk11 sk22 · · · sk66 (1 ± P12 )(1 ± P12 )w(β D2D2 ; RD21A D21B )]

(25)

A B The function (1 ± P12 )(1 ± P12 )w is symmetric in PAB , so that one can use the symmetrized

s polynomials, such as h+ , and drop PAB . In the actual fit, we neglected functions corresponding to minus signs in Eq. (25), as specified in Table 4. 4.5.7

D2-H terms

The action of the operator (1 + E) on w(β HD2 ; RD21A H1B ) produces a function symmetric A in P12 . Therefore, the product sk11 sk22 · · · sk66 in a generic D2H term can only contain even

powers of s3 . With this constraint, a general symmetric function can be written as B A (1 + PAB )[sk11 sk22 · · · sk66 (1 + (−1)k6 P12 )(1 + P12 )w(β HD2 ; RD21A H1B )],

A ]w. which is analogous to the OH case with function w replaced by [1 + P12


(26)


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

5

Quality of flexible-monomer potentials

5.1

Performance of fits on training data

The quality of a fit can be judged by computing its global rmse relative to the ab initio values. One can also examine the distribution of individual errors on scatter plots such as Fig. 1 displaying the values predicted by the SAPT-5s′ f[2014] fit versus ab initio energies. Panel (a) shows all points with energies less than 110 kcal/mol, while panel (b) shows points in the most important energy range [-10, +20] kcal/mol. In both panels, one sees the expected uniform increase of error with increase of energy (although for energies larger than about 20 kcal/mol, the error becomes constant). This behavior is due to the weighting factor used which is designed to emphasize the important low-energy region. The large discrepancies at highly repulsive energies are of little importance for the anticipated applications of SAPT5s′ f[2014]. To further analyze the accuracy of our fit, the rmse’s of the flexible-monomer potential can be compared to those of the rigid-monomer potential in all three ranges of interaction energy shown in Table 1. The ratios of the SAPT-5s′ f[2014] and SAPT-5s′ overall rmse’s is 7.8. This fairly large ratio is mainly due to the much larger number of points with very large energies—which always have large absolute errors—in SAPT-5s′ f[2014] compared to SAPT-5s′ . These relation can be seen in Table 1: if the 330 points with energies larger than 150 kcal/mol (irrelevant for all intended applications of our potential) are excluded, the rmse of SAPT-5s′ f[2014] drops to only 0.82 kcal/mol. Since there are no such points in the SAPT-5s′ case, its rmse remains the same and the ratio drops to 2.4. The analogous ratio is close to two for all the remaining ranges of energy. In general, one strives to make the errors of a fit smaller than the inherent errors of the ab initio interaction energies (due to incompleteness of basis sets and the truncation of the level of theory). In simpler cases (i.e., dimensionality of potential surface between 1 and 4), we were able to obtain fits with rmse’s 5–10 times smaller than the uncertainty of the fitted data. For the water dimer, our ab initio


Page 30 of 89

Page 31 of 89

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


calculations have an uncertainty of about 5–10% or 0.35 kcal/mol near the minimum, mainly due to the use of a relatively small basis set. The SAPT-5s′ fit has an rmse in the attractive region of 0.10 kcal/mol or about 2%, i.e., about 3 times smaller than the uncertainty of the ab initio calculations. Since the latter uncertainty is the same for flexible-monomer points, the 0.24 kcal/mol or 4% rmse of SAPT-5s′ f[2014] in the attractive region is still smaller than this uncertainty. Whereas such accuracy is satisfactory at the present time, in future work, as the accuracy of data points increases, the flexibility of the fitting function will have to be increased. The quality of the SAPT-5s′ f[2014] fit could be improved if the nonlinear parameters of the fit were expanded in powers of intramonomer coordinates. Table 1 shows that rmse’s of SAPT-5s′ fIR[2014] are from 22% to 26% larger than those of SAPT-5s′ f[2014], except for rmse’s on all points which are almost the same. In particular, for Eint < 0 the two uncertainties are 0.30 and 0.24 kcal/mol, or 25% different. This shows that the anchoring of the potential surface at the minimum reduced the accuracy in all other regions fairly significantly. In Fig. 2(a), we present a graphical comparison of SAPT5s′ fIR[2014] with the SAPT ab initio energies. The plot looks very similar to that shown if Fig. 1(b) despite the discussed differences in rmse’s. To examine the difference between the two potentials, Fig. 2(b) shows the correlation diagram of SAPT-5s′ f[2014] vs. SAPT5s′ fIR[2014] at the grid of points of the ab initio calculations. This plot demonstrates that the differences between the two surface are similar to the differences between a given surface and ab initio data. Clearly, the former difference may result in differences in the spectra of the dimer observed in ref 79, and will be seen in the differences in the calculated second virial coefficients in Sec. 7.4. A further discussion of accuracy of the SAPT-5s′ f[2014] and SAPT-5s′ fIR[2014] surfaces is given in Secs. SI-VIII and SI-XI of the Supporting Information. In Sec. SI-V A of the Supporting Information, we present an analysis of the quality of the potential energy surface through investigations of the radial scans for the three selected intermolecular spatial configurations with monomers deformed in various ways. In Sec. SI-V B of the Supporting Information, we analyze the performance of the potentials as functions



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 the internal coordinates.

5.2

Performance of fits on ab initio energies of Babin et al.

The supplementary material of ref 85 includes about 40 thousand interaction energies computed at the CCSD(T)/CBS level. These values provide an excellent test set for our potential. In fact, Babin et al. included in their work (also in supplementary materials) a graph showing the performance of the CCpol-8sfIR[2012] potential on these data. (Note that the CCpol-8sfIR[2012] interaction energies used in this comparison were incorrect. This error was corrected in an erratum, ref 105). Table 5 shows that the rmse of CCpol-8sfIR[2014] is 0.42 kcal/mol on all 40K points from ref 85. Thus, surprisingly, this rmse is much smaller than the 3.11 kcal/mol rmse of the SAPT-5s′ fIR potential on the training set used in fitting this potential. One reason for this behavior is that the latter set contains a significant number of very large interaction energies, in particular 330 energies are larger than 150 kcal/mol, while the largest interaction energy in the former set is below 150 kcal/mol. However, our potentials should not be compared to the full 40K set of interaction energies from ref 85 since these potentials were fitted to monomers’ distortions corresponding to the classical turning points in the ground rovibrational state of the water monomer, whereas the set of ref 85 contains many points beyond this range (some of these points include monomers deformed so significantly that their energies are more than 120 kcal/mol above equilibrium). In Fig. SI-16 of the Supporting Information, we present correlation plots illustrating the performance of the CCpol-8sfIR[2014] potential on the ab initio values from ref 85 for all reference geometries from ref 85, as well as results for the subsets of these geometries inside and outside of ground-state vibration. These plots show a very good correlation for the inside geometries. Table 5 shows that for these geometries the rmse drops to 0.19 kcal/mol. This rmse is smaller than it could have been expected taking into account the inherent uncertainty of the SAPT ab initio training data amounting to about 0.3 kcal/mol near the minimum 32 ACS Paragon Plus Environment

Page 32 of 89

Page 33 of 89

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 0.30 kcal/mol rmse of the SAPT-5s′ fIR potential on the training interaction energies below zero (whereas the test data are probably accurate to below 0.1 kcal/mol near the minimum). Clearly, the improved accuracy is due to the inclusion of the rigid-monomer CCpol-8s potential (rmse of 0.10 kcal/mol for all 2510 rigid-monomer grid points and 0.01 kcal/mol for points with interaction energies below zero) in the hybrid approach. If one views CCpol-8sfIR[2014] as CCpol-8s with a “correction” for the flexibility effects, if this correction is small and if the errors in the two terms in the square bracket in Eq. (12) cancel effectively, the hybrid potential can be more accurate than the underlying flexible-monomer potential. The rmse of 0.19 kcal/mol shows that our functional form is sufficiently elaborate relative to the uncertainties of the underlying training interaction energies. As mentioned earlier, fitting such data with accuracy much higher than their uncertainty is not advisable since it makes applications more time consuming.

5.3

Spectra of water dimer

In ref 82, the CCpol-8sfIR[2012] surface (called CCpol-8sf there) was used to calculate the vibration-rotation-tunneling spectra of the water dimer, including the intramonomer vibrational transitions. The results of those calculations were extensively compared with experimental spectra. We checked how the refit leading to the CCpol-8sfIR[2014] surface affected the theoretical spectra. Since the intermolecular frequencies are only mildly dependent on the intramonomer degrees of freedom, 82 such frequencies cannot be impacted by the refit. Thus, we have limited our test to the comparison of the frequency shifts of the intramonomer modes. Since the 12-dimensional fully anharmonic rovibrational calculations, such as those performed in ref 82, are very time consuming and at the same time the harmonic approximation works reasonably well for intramonomer vibrations, we have used the latter approximation. The results are presented in Table 6. For compatibility with ref 82, we used Radau f = 1 embedding, although we recommend now the Eckart embedding (results with the latter embedding are listed in the last column and are negligibly different 33 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

from Radau-embedding results). The original frequency shifts 82 differ by up to 0.4 cm−1 from those computed by us using the same potential. The reason for these discrepancies is that while we now use the standard approach of the diagonalization of the Hessian matrix of the potential, calculations of ref 82 arrived at the harmonic level treating it as a limiting case of the fully anharmonic solution, which amounts to a different numerical procedure. The effects of the refit are seen by comparisons of the third and fourth columns of Table 6: the largest differences appear for the modes O-Hb [D] of the donor and ss[A] of the acceptor and amount to 2.31 cm−1 and 1.05 cm−1 , respectively, both modes related to the symmetric stretch in monomer. Thus, there is a clear correlation with the symmetry error discussed in Sec. 1.3 which involved the coordinates s1 and s4 , representing such stretches. Although the error was corrected in ref 82, the fit was not reoptimized, so a partial effect of it remained. Since the observed discrepancies are within the experimental error bars for the frequency shifts 82 and we expect the discrepancies for the anharmonic frequency shifts to be very similar, the conclusions of ref 82 remain valid. The reasons for the discrepancies between theory and experiment for the frequency shifts have been discussed in ref 82 by Leforestier et al.. These authors gave three possible explanations. First, experimental data have unknown and possibly substantial uncertainties. Second, the adiabatic [6+6]D approximations used in calculations has also unknown uncertainties, although these uncertainties are expected to be small. There has been no new information available on these two points since ref 82 was published. Third possibility is, of course, inaccuracy of the potential energy surface. New light on this issue was recently shed by the work of Babin et al. 85 who evaluated the shifts using a different potential surface and the same [6+6]D model. They obtained only minor reductions of discrepancies with experiments compared to ref 82, which suggests that theoretical predictions cannot be changed significantly by further improvements of the water dimer potentials.


Page 34 of 89

Page 35 of 89

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

Characteristic points on energy surface

We have performed a search of the 12-dimensional potential energy surface to find characteristic points. The results and comparisons with literature are presented in Tables 7 and 8. The potential surface U is defined as U = UA + UB + V,

(27)

where UX is the monomer potential (PJT2 71 ) with the zero value at the equilibrium geometry of the molecule and V is the flexible-monomer intermolecular (vertical) interaction potential, defined by Eq. (5) in the case of SAPT-5s′ f/IR or by Eqs. (12) and (5) in the case of the hybrid CCpol-8sf/IR potentials. Thus, U is the energy required to bring two monomers initially at their equilibrium geometries from infinity to a given dimer geometry with deformed monomers. Equivalently, U is the total potential energy of the system measured with respect to the energy of two equilibrium-geometry monomers at infinite separation.

6.1

Minimum

The global minimum of U and its geometry for each of our potentials are given in Table 7 and compared to data from best literature calculations. 6.1.1

Geometry

Deformation of monomer geometry: For our potentials, the changes of intramonomer geometries upon dimerization relative to the re geometries can be seen comparing to the column SAPT-5s′ f[2014](re ). The most accurate calculation for the minimum configuration are those of Lane, 106 but, unfortunately, Lane does not give the best estimates of the isolated monomer geometry. Thus, we have to compare with the previous best benchmarks of Tschumper et al. from ref 107. These authors obtained the equilibrium geometry of an isolated water molecule with rOH = 0.9589 ˚ A and θ = 104.16◦ from CCSD(T) optimiza35 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

tions in a triple-zeta-quality basis set. The changes of his geometry upon dimerization are: -0.0008 ˚ A, 0.0064 ˚ A, 0.0008 ˚ A, 0.3◦ , and 0.4◦ , for the free OH of the donor, bonded OH, free OH’s of the acceptor, θ of the donor, and θ of the acceptor, respectively. As expected, the best agreement with these benchmarks is achieved by the CCpol-8sfIR[2014] potential, which gives the analogous deformations equal to -0.0004 ˚ A, 0.0039 ˚ A, 0.0011 ˚ A, 0.2◦ , and 0.1◦ . Taking into account the smallness of the deformations, the agreement is reasonable. In fact, since for some other quantities CCpol-8sfIR[2014] predictions are closer to Lane’s benchmarks than the Tschumper et al. predictions, it is not clear which set of numbers is more accurate. Note that the relative changes are very small, of the order of only 0.1% of bond lengths and angles. These small geometry changes result in the total energy U lowered by only 0.012 kcal/mol between SAPT-5s′ f[2014](re ) and SAPT-5s′ f[2014]: the vertical interaction energy V increases in magnitude by 0.024 kcal/mol, but this increase is compensated by 0.012 kcal/mol increase of the monomer energies in U . These energetic changes are of the order of only 0.5% of the interaction energy. The increase of monomer energies is larger, 0.031 kcal/mol, in calculations of ref 107, in accordance with the larger deformation of the bonded OH distance. Dimer geometry: For all potentials and calculations included in Table 7, the mutual orientations of monomers, i.e., the values of the angles α and β, are within experimental error bars, 108 so that the experiment cannot differentiate between current theoretical results. The angle α describing the nonlinearity of the hydrogen bond is within 0.3◦ of Lane’s value for all SAPT and CCpol potentials except for SAPT-5s which differs by 0.6◦ . Somewhat surprisingly, Tschumper’s et al. 107 result differs by 0.95◦ . The discrepancies are larger for the angle β describing the orientation of the acceptor plane relative to the O-O axis and the SAPT potentials give β from 3◦ to 6◦ larger than Lane’s value. For CCpol potentials, the largest discrepancy is much smaller, only 1.6◦ , and CCpol-8sfIR[2014] agrees with Lane’s value to all digits listed. Tschumper’s et al. result differs by 1.5◦ . It appears that the


Page 36 of 89

Page 37 of 89

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


inclusion of monomer flexibility has a minor effect on the values of the angles, in particular for SAPT-5s′ and SAPT-5s′ f[2014](re ) vs. SAPT-5s′ f[2014]. The performance for the ROO distance is similar to that for the angle β. All SAPT potentials give the value of about 2.96 ˚ A, 0.05 ˚ A larger than the benchmark. 106 This behavior is due to basis set deficiencies, as shown in ref 59. The CCpol-8s 62 potential, fitted to ab initio calculations in large basis sets, agrees with Lane’s value to 0.001 ˚ A. The inclusion of monomer-flexibility effects makes the agreement slightly worse, to within 0.003-0.005 ˚ A. However, even in the last case the relative error of 0.2% is very small. Tschumper’s et al. value of ROO agrees extremely well with Lane’s value, to 0.0003 ˚ A, which is surprising in view of larger discrepancies in the angular orientation. This is not due to the post-CCSD(T) effects included by Lane since his best estimate at the CCSD(T) level is 2.9095 ˚ A, not much different. We believe that this agreement is due to the cancellation of two errors: the BSSE included in the geometry optimizations of ref 107 tends to shorten the bonds, whereas the use of medium-size basis sets tends to make them longer. The results discussed above show that the CCpol-8sf/IR[2014] potentials make the minimum geometries of the dimer much closer to the benchmark values than is the geometry of SAPT-5s′ f[2014], supporting the hybrid method of constructing such potentials. One may also point out that the geometries from the hybrid potentials are of similar accuracy to those of the previous best benchmarks of Tschumper et al. 107 The parameters of the minimum obtained using the SAPT-5s′ rigid-monomer potential with hri0 monomer geometry, which is the same as SAPT-5s′ f[2014](hri0 ), can be compared to parameters produced by SAPT-5s′ f[2014](re ), which is equivalent to using equilibrium monomers. As already mentioned, the hri0 potentials provide a better rigid-monomer approximation to full dimensionality potentials in terms of reproducing the rovibrational spectra than the re potentials. The same feature does not need to be true for the intermonomer geometries at the global minima since, as discussed above, the intramonomer coordinates at such geometries are close to the re ones. Results in Table 7 show that the use of the hri0 monomer geometry instead of the re geometry does improve the values of ROO and 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

whereas the magnitude of the discrepancy in α is unchanged. However, the improvements are rather small. 6.1.2

Potential depth De

The potential depth De = −Umin is a single number characterizing the 12-dimensional total potential energy surface. This quantity should be distinguish from DeV = −Vmin which is different for each pair of monomer geometries and characterizes a given rigid-monomer 6-dimensional interaction potential. The difference between these two potential depths is equal to UA + UB at the global minimum and therefore is very small for all potentials included in Table 7. This table also shows an excellent agreement, to within 0.013 (0.005) kcal/mol, of Umin calculated from the CCpol-8sf[2014] (CCpol-8sfIR[2014]) with the best theoretical estimation from ref 106 of -5.0208 kcal/mol, obtained at the counterpoise-corrected CCSD(T)-F12b/CBS limit with corrections for higher order excitations, core correlation, relativistic effects, and diagonal Born-Oppenheimer effects. The sum of the post-CCSD(T) effects is only about -0.01 kcal/mol, and Lane’s result at CCSD(T) level amounting to 5.0103 kcal/mol agrees even better with our value: to within 0.002 (0.006) kcal/mol. Such very close agreement is better than one could have expected taking into account the uncertainty of SAPT-5s′ f, but apparently the high accuracy of CCpol-8s determines the overall accuracy. Clearly, the level of theory and the size of the basis set used in development of the CCpol-8s potential were sufficient. The results of Tschumper et al. 107 are 0.04 kcal/mol above Lane’s values. One interesting observation from Table 7 is that the vertical interaction energy Vmin from CCpol-8s agrees also quite well, to 0.08 (0.05) kcal/mol, with Umin (Vmin ) of Lane. This is perhaps surprising since the hri0 monomers used in CCpol-8s lie 2×0.2226 kcal/mol above the re monomers. Apparently, the rigid-monomer surface with monomers at their geometries from the global minimum optimized by Lane is to a high degree parallel to the CCpol-8s surface. Since the energies of these monomers are only 0.03 kcal/mol above the energies of


Page 38 of 89

Page 39 of 89

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


equilibrium monomers, this results in the observed agreement.

6.2

Saddle points

Since ref 106 does not provide information about other than minimum stationary points of the water dimer, we have used the results of Tschumper et al. 107 as the benchmarks. Table 8 lists the values of the interaction energies at stationary points relative to the potential minimum, i.e., the barrier heights. The CCpol-8sf[2014] and CCpol-8sfIR[2014] surfaces give very similar values of the barriers, with rmse for the former one 3 cm−1 smaller than for the latter. The largest discrepancy for CCpol-8sf[2014], observed for the stationary point number 10, is equal to 18 cm−1 , but this amounts to only 2% of the barrier height. The largest percentage error is for point 3 and amounts to 4% of the barrier height. The rmse of the CCpol-8sf[2014] fit is only 10 cm−1 , very small indeed taking into account that the estimated uncertainties of the benchmark barriers are between 5 and 23 cm−1 . It is also worth mentioning that Tschumper’s et al. barriers include the so-called Brueckner doubles contributions and relativistic effects. If both effects (amounting to 9 cm−1 or less) are subtracted from these barriers, the agreement with CCpol-8sf[2014] further improves and the rmse drops to 6 cm−1 . We believe this is an unprecedented agreement of predictions given by a potential function with single-point ab initio calculations. The barrier heights calculated with the rigid-monomers CCpol-8s surface are much further from the reference values, rmse equal to 47 cm−1 , than in the case of the flexible-monomers surfaces. Since both hybrid potentials CCpol-8sf[2014] and CCpol-8sfIR[2014] are constructed on the base of CCpol-8s, the observed significant improvements are due to the inclusion of the flexibility effects. The flexibility corrections for the barriers observed in our calculations are in the very good agreement with the estimates obtained in ref 62.



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

7

Second virial coefficient

Until very recently, the majority of published calculations of the second virial coefficients for molecules assumed rigid monomers. Most papers that used flexible-monomer potentials averaged these potentials before using them in calculations of virial coefficients. One approach was to average the potentials over quantum ground vibrational wave functions of monomers, see for example ref 109, and the first such calculations for water were performed in refs 81 and 79. The approximations made in deriving this approach are examined in Section 7.1 and were also recently discussed in ref 110. An alternative approach consisted in applying purely classical Boltzmann’s averaging over internal coordinates, 111 i.e., one performs integration of the Mayer function dependent on a flexible-monomer dimer potential with a weight equal to the product of the Boltzmann factors of intramonomer potentials UA and UB . A short derivation of this approach is presented in Section 7.1.1. This fully classical approach is commonly employed in the theory of flexible-monomer polymers (see, e.g., ref 112), where the intramolecular degrees of freedom are soft. Refson et al. 111 applied it to the case of water, which does not seem to be well justified due to the manifestly quantum nature of intramonomer vibrations in the temperatures involved. To our knowledge, no comparisons of the two approaches are available. After a flexible-monomer potentials has been averaged over the intramonomer coordinates, it can be used in the same way as a rigid-monomer potential in calculations of the classical second virial coefficient. However, for water it is necessary to go beyond this level and include quantum effects. To this end, one can use the perturbative approach, 39,113,114 the Takahashi-Imada (TI) approximation, 115–117 or finally calculate this effect virtually exactly (within the rigid-monomer or averaged-potential approximations) using the Path Integral Monte Carlo method (PIMC). 85,104,110,118,119 Very recently the problem of simultaneous accounting for both monomer-flexibility and quantum effects has been addressed in several papers, for water 85,117,119 and for other molecules (e.g., refs 104 and 110). The calculations presented here will use quantum-averaged potentials at the classical level and the TI method of including quantum effects applied with our rigid-monomer potential. 40 ACS Paragon Plus Environment

Page 40 of 89

Page 41 of 89

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


7.1

Second virial coefficient from flexible-monomer potentials

The general expression for the homogeneous second virial coefficient, true both in the classical and quantum statistical mechanics, reads 120,121 V B(T ) = − 2

2Q2 −1 , Q21

(28)

where Q1 and Q2 are the partition functions of the monomer and the dimer, respectively, and V is the volume. Let us denote the coordinates of the nuclei of the two (identical) monomers A and B in a space-fixed coordinate system by X A and X B , respectively, so that the potential energy of Eq. (27) can be written as: U (X A , X B ) = UA (X A ) + UB (X B ) + V (X A , X B ),

(29)

If we denote the kinetic energies of nuclei in A and B by TA and TB , respectively, a nuclear motion Hamiltonian for the whole complex reads H = TA + UA + TB + UB + V , while for the monomers we have HA = TA + UA and HB = TB + UB . These Hamiltonians fully determine the partition functions Q1 and Q2 . Note that we make here an assumption that N atoms in each monomer form a distinct molecule at all times, i.e., we ignore dissociation processes. 7.1.1

Purely classical approach

In classical statistical mechanics, one can factor out and perform analytically simple integrations over momenta of all atoms, giving a factor (2πmi /β)3/2 , where β = 1/kB T , for each atom of mass mi . These contributions are identical in Q2 and Q21 and therefore cancel in expression (28). What remains depends only on configurational integrals: R R dX A dX B e−β(UA +UB +V ) V V dX A dX B e−β(UA +UB ) fAB R R R Bcl (T ) = − −1 =− R 2 2 dX A e−βUA dX B e−βUB dX A e−βUA dX B e−βUB


(30)


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 42 of 89

where fAB = e−βV − 1. The factor 2 in front of Q2 in Eq. (28) cancels with the Gibbs permutational factor included in Q2 . This factor is just 1/2! for identical monomers as we assume that the monomers cannot separate into atoms. There can be additional Gibbs factors within monomers (for example, a factor 1/2! due the presence of two identical hydrogen atoms, but these factors cancel out between numerator and denominator in Eq. (28)). The integral in the numerator depends on 6N coordinates, but the potentials are independent of the dimer’s center of mass (COM) position and of overall rotations of the dimer, so one could immediately integrate over six coordinates. It is more convenient to proceed somewhat differently and compute first the COM’s of both monomers, express the position of the COM of monomer B relative to the position of the COM of monomer A (we will denote this vector by R), and define monomer-embedded coordinate systems for each monomer. The orientation of monomer’s A system with respect to the space-fixed system can be determined by three Euler’s angles ωA , and similarly for monomer B. The positions of atoms of each monomer in its coordinate system will be denoted by rA and rB (for the water dimer, these variables can be the same as defined in Sec. 3.2). The potentials of Eq. (29) can be written in terms of the new coordinates: UX = UX (rX ) and V = V (R, ωA , ωB , rA , rB ). We can now integrate over the center of mass of monomer A, which gives a factor V, and an analogous procedure for each of the monomers gives a factor V 2 in the denominator, so that the second virial coefficient can be written as R drA dωA drB dωB e−β(UA +UB ) dRfAB 1 R R R Bcl (T ) = − R 2 drA e−βUA dωA dωB drB e−βUB R drA drB e−β(UA +UB ) dRdωA dωB fAB 1 R R , = − 2(8π 2 )2 drA e−βUA drB e−βUB

(31)

where drX , dωX , and dR denote the appropriate volume elements, we have used the fact that R the monomer potentials do not depend on the orientation of monomers, and dωX = 8π 2 .

Thus, the classical second virial coefficient can be interpreted as a Boltzmann average of the Mayer factor containing the dimer interaction potential with flexible monomers (the fAB


Page 43 of 89

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


factor) over all possible configurations of each monomer. Note that there were no approximations made in deriving Eq. (31). An expression equivalent to Eq. (31) was apparently first derived by Refson et al. 111 and later used, for example, in refs 122 and 112. We have not used this formula in our work, but the one described in the next subsection that utilizes quantum averaging over monomer vibrations. 7.1.2

Quantum mechanical approach

Most calculations of molecular second virial coefficients use rigid-monomer potentials. Then the monomers are just rigid bodies and one can write the partition functions in Eq. (28) for the case of two interacting rigid bodies. The approximations leading to this form of B(T ) were recently discussed in detail in ref 110. However, for the case of averaged potentials hV i0 , the authors of this reference only stated that since hV i0 depends on the same set of coordinates as the rigid-monomer potential, it can replace the latter potential in the B(T ) formula for interacting rigid bodies. We will present a more convincing derivation of this formula below using quantum statistical mechanics to determine the partition functions in Eq. (28). One should mention that there is no direct way of averaging V in Eq. (28) which would give hV i0 since the quantum partition functions are traces of e−βH . Our derivation is based on representing this partition function as an expansion in eigenstates of a reduceddimensionality Hamiltonian and arguing that the use of hV i0 in this Hamiltonian provides the best set of eigenstates. The first step in deriving reduced-dimensionality formulas for B(T ) is a partial approximate separation of vibrational and rotational motions in monomers. The monomer partition functions separate rigorously into products of the translational and rovibrational parts: Q1 =

q1tr q1int

=V

2πmkB T h2

3/2

q1int

(32)

where m is the mass of monomer and h is the Planck constant. To approximately separate the vibrational and rotational motion, we partly neglect the coupling between the corresponding 43 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 44 of 89

degrees of freedom, so that the vibrational motion is assumed to be that for the zero angular momentum. We will denote the energy levels of such an oscillator by ǫi . However, while the rotational motion will be that of rigid monomer, for each vibrational state i, monomer’s geometry will be the geometry averaged over the ith vibrational wave function. Thus, we can write q1int as q1int ≈

X

e−βǫi q1rot,i .

(33)

i

The functions q1rot,i can be taken in the classical limit 121,123 since kB T at the temperatures of interest in investigations of B(T ) are high for most molecules compared to rotational excitation energies. Thus, for asymmetric tops we have: q1rot,i

29/2 π 7/2 ≈ (IA,i IB,i IC,i )1/2 σ

kB T h2

3/2

1 = σ

π Ai Bi Ci

1/2

kB T h

3/2

,

(34)

where Ai , Bi , and Ci are the average rotational constants of the molecule in the vibrational state i, defined in terms of principal moments of inertia IX as Ai = h/8π 2 IA,i , and σ = 2 for water is the so-called symmetry number, reflecting the fact that the two hydrogen atoms are indistinguishable. For the dimer, after separating the translational and internal motions, Q2 = V

2π2mkB T h2

3/2

q2int ,

(35)

the partition function q2int is defined in terms of the total rovibrational energies of the dimer, including the motions within monomers. Since the latter motions are much faster than the motions of monomers relative to each other, we can separate the two types of motion and write the energies as En;iA iB + ǫiA + ǫiB , where En;iA iB is the energy of dimer’s motion in coordinates R, ωA , and ωB only, while the monomers are in the vibrational states of energies


Page 45 of 89

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


ǫiX , same as in the monomer case discussed above. The internal partition function becomes: q2int =

1 X −β(ǫi +ǫi ) X −βEn;i i B A A B. e e 2! i ,i n

(36)

A B

We will assume that the sum over n runs over states which have no symmetry with respect to the exchange of two monomers, which results in the factor 1/2! The energies En;iA iB are not fully defined yet and several choices can be made. One choice is to assume that the dimer “sees” the monomers frozen in their averaged vibrational positions hrX iiX . Then En;iA iB are eigenenergies of the Hamiltonian with the (6-dimensional) interaction potential V = V (R, ωA , ωB , hrA iiA , hrB iiB ). However, it was shown in ref 78 that the energy spectrum of such a system agrees much better with the spectrum from full-dimensional calculations if the averaged potential hV iiA iB is used instead. With these assumptions, Eq. (28) becomes: # " P P h3 iA iB e−β(ǫiA +ǫiB ) n e−βEn;iA iB 1 B(T ) = − −V , P P 2 (2πµkB T )3/2 i e−βǫiA q1rot,iA i e−βǫiB q1rot,iB A B

(37)

where µ = m/2. Since the lowest excitation energy of water monomer is equal to kB T at T ≈ 2400, only a few first terms in the sums over iX need to be considered and only at the highest temperatures of interest for investigations of B(T ). Each such term can be computed relatively easily due to the reduced dimensionality. The method can be applied both with the rigid-monomer potential computed for appropriate averaged geometries V (R, ωA , ωB , hrA iiA , hrB iiB ) and with hV iiA iB . Thus, Eq. (37) offers a way for calculations of high-T virial coefficients which is an alternative to the methods developed in refs 104 and 110. One should emphasize that the high-T performance of Eq. (37) is achieved only with the noted dependence of the potential and of the rotational partition functions on iX . If this dependence is neglected, the vibrational partition functions cancel in the numerator and denominator and the result is the same as if the summations over iX were reduced to the first term only

" # P h3 n e−βEn;00 1 B(T ) = − −V , 2 (2πµkB T )3/2 q rot,0 2 1 45


(38)


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 46 of 89

This expression is of the same form as that for rigid-monomer potentials often used in literature. The expression (38) can by computed by direct summation over dimer energy levels or using the PIMC method. It can also be computed accounting for quantum effects at the level of the TI approximation. In fact, despite derivation from quantum mechanics, it can also be computed from the classical formula, replacing the sum over n in Eq. (38) by the phasespace integral of the Boltzmann factor with the rigid-monomer Hamiltonian containing the potential hV i0 . This Hamiltonian is the sum of the rotational Hamiltonians for monomers A and B, of the kinetic energy of a particle of mass µ (described by the vector R), and of hV i0 (R, ωA , ωB ). The integration of the first term over the momenta conjugated to ωA gives q1rot,0 /8π 2 , where the denominator is needed since the partition function of Eq. (34) does contain the factor 8π 2 resulting from position-space integration over ωA . The integration of the third term over the momenta conjugated to R gives (2πµkB T )3/2 /h3 . Thus, the virial coeffcient of Eq. (38) can now be written as Z Z 1 1 1 −βhV i0 Bcl (T ) = − −V =− dR dωA dωB e dR dωA dωB e−βhV i0 − 1 2 2 2 2 2 (8π ) 2(8π ) Z

1 (39) dR e−βhV i0 − 1 ω ω , = − A B 2 where we used V = (1/(8π 2 )2 )

R

dR dωA dωB 1. We added the subscript “cl” to indicate the

classical approximation, however, this is a different expression than Eq. (31) derived purely classically.

7.2

Final second virial coefficient

If a potential averaged over intramonomer coordinates is used, the calculation of the classical second virial coefficient reduces to the rigid-monomer case, see Eq. (39): 1 Bcl (T ) = − 2

Z

hf (T )iωA ωB dR,


(40)

Page 47 of 89

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


where T is temperature, < · >ωA ωB denotes averaging over monomer orientations, and f (T ) is the Mayer function written in terms of the vibrationally averaged potential hV i0 f (T ) = exp(−hV i0 (R, ωA , ωB )/kB T ) − 1,

(41)

with kB denoting the Boltzmann constant. The expression (40) includes monomer-flexibility effects in an averaged way, so that we can define a flexibility correction at the classical level, hV i

∆Bflex 0 , as hV i

hV i0

∆Bflex 0 = Bcl

V (hri0 )

− Bcl

,

(42)

where both terms on the right-hand side are given by Eq. (40), the first one with the Mayer function of Eq. (41) and the second one replacing hV i0 in this function by the rigid-monomer potential V (hri0 ). We will augment the virial symbols with acronyms of our surfaces. For hV i0

instance, Bcl

[CCpol-8sf] denotes the classical second virial coefficient calculated with the

(averaged) CCpol-8sf surface. Our final second virial coefficient, taking into account both monomer-flexibility and quantum effects, is defined as hV i

TI Bflex (T ) = BQM (T ) + ∆Bflex 0 (T ),

(43)

TI where BQM is the virial coefficient in rigid-monomer approximation with quantum effects

included in the TI approach, computed in the same way as in ref 60. The expression for TI is the same as the classical expression (40), but with the Mayer function defined in BQM

terms of an effective operator Veff of the form: 115,116 Veff

" # 2 2 X h ¯2 Fi2 X ταi = Vrigid + + , 24(kB T )2 i=1 M I α α

(44)

where Fi is the force on the molecule i exerted by the partner and ταi , α = 1, 2, 3, are the components of the torque on the molecule i along the principal axis α of this molecule.



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 48 of 89

The mass of the molecule and its principal moments of inertia are denoted by M and Iα , respectively. The forces and torques are computed from the potential Vrigid .

7.3

Vibrational averaging

A 12-dimensional water dimer potential will be reduced to the 6-dimensional averaged potential by integrating its product with with the squared ground state vibrational wave functions over intramonomer coordinates monomers: hV i0 (R, ωA , ωB ) = hV (R, ωA , ωB ; rA , rB )i0 =

Z

|χA (rA )χB (rB )|2 V (R, ωA , ωB ; rA , rB )drA drB . (45)

We used very accurate wave functions χX obtained using programs of Tennyson et al. 124,125 To our knowledge, it is the first time that the averaging defined by Eq. (45) is applied in calculations of water virial coefficients. However, this averaging has been extensively tested in literature on properties of smaller systems where comparisons to either full-dimensional results or to accurate experiments were possible 78,126,127,127–130 and always performed very well. To implement Eq. (45), we have to specify the embedding of monomers within the dimer. Although we have used principal axes embedding to generate grid points, the SAPT5s′ f[2014] and SAPT-5s′ fIR[2014] potentials are embedding independent. Therefore, we could use the more appropriate EAE in calculations of virial coefficients. As a consequence, we used EAE also for constructing the hybrid potentials CCpol-8sf[2014] and CCpol-8sfIR[2014], rather than the Radau (f = 1) embedding used in ref 82. The numerical details of our calculations of virial coefficients are given in Sec. SI-X of the Supporting Information. As described there, we have performed averages on-the-fly, i.e., for each numerical integration point in the six-dimensional integral of Eq. (40). This method is time consuming since one needs at least a few million points in such integration even using specialized quadratures such as the one developed by Conroy. 131 A much less expensive method would be to perform the average for a few thousand points, for example


Page 49 of 89

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


on the SAPT-5s grid of 2510 points, and then fit the averaged energies. We have not followed this way in order not to introduce additional uncertainties which would be resulting from this fitting. The averaging in Eq. (45) can be also performed using a truncated Taylor expansion (TE) of the potential in internal parameters s, see refs 126 and 132. Since the use of such expansion reduces the averaging of V to the averaging of the components of s and their products, it also dramatically reduces the costs. Our investigations of TE averaging in calculations of virial coefficients are presented in Sec. SI-XIII of the Supporting Information.

7.4

Results

The results of our calculations of the second virial coefficient are presented in Tables 9 and 10. In Table 9, we analyze the monomer-flexibility effect on Bcl (T ). The reference values are V (hri0 )

Bcl

V (re )

(T ), but we also list Bcl

(T ) to show to what extent the use of the potential with

monomers at re vs. at hri0 affects the quality of results. The rigid-monomer potentials are equally costly to obtain independently of the choice of monomer geometry. Note that the V (hri0 )

values of Bcl

(T ) are the same for both the CCpol-8sf[2014] and the CCpol-8sfIR[2014] V (re )

potentials. Table 9 shows that at the lowest temperature Bcl V (hri0 )

above Bcl

hV i0

. Thus, since Bcl

V (hri0 )

lies below Bcl

lies about 200 cm3 /mol

, the monomer-flexibility correction

computed relative to the former potential would have been much larger in magnitude for V (re )

low temperatures. The absolute difference between Bcl

V (hri0 )

and Bcl

gradually diminishes

hV i

with T and becomes small compared to ∆Bflex 0 above 1000 K, but the relative difference decreases only from 10% (9%) at T = 273.15 K to 4% (6%) at T = 1000 K for CCpol-8sf (CCpol-8sfIR). The monomer-flexibility contribution given by hV i0 computed from the CCpol-8sf[2014] potential, listed in Column 4, is negative and fairly small in magnitude for low T : it ranges from 1.0% at T = 273.15 K to 2.6% at 323.15 K. The relative importance of this correction continues to increase with T and becomes as large as 13% at 773.15 K. Above this tem49 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

perature, the relative correction becomes gradually larger and larger, but this is obviously related to the fact that B(T ) becomes small in magnitude to cross zero between 1400 K and 1500 K. The values of the correction remain negative and decrease smoothly with T . After the crossing of zero, the relative correction starts to decrease and at the highest T = 3000 K, it becomes equal to 8.3%. Thus, the monomer-flexibility correction behaves opposite to the quantum correction both in sign and in the dependence on T . The monomer-flexibility correction computed from the CCpol-8sfIR[2014] potential and listed in Column 6 is a few times larger in magnitude than that computed from the CCpol8sf[2014] potential. This finding was surprising for us since the performance of the two potentials on dimers spectra is not that much different. 79 Our explanation is that this behavior is due to the excessive emphasis on the minimum region in fitting the SAPT-5s′ fIR[2014] potential. As described earlier, this potential was fitted to additional points in the vdW minimum region taken with large weights and to ab initio computed second derivatives of interaction energy. This resulted in not so good description of the regions farther from the minimum which are important for calculations of virials, but not for calculations of dimer’s spectra, see the analysis of this issue in Secs. SI-VIII and SI-XI of the Supporting Information. Therefore, although we included the virial coefficients computed with the CCpol-8sfIR[2014] potential in Table 9, we do not recommend these values to be used. hV i

The relative importance of ∆Bflex 0 analyzed above is very different from that published by Donchev et al. 133 who found 20% to 40% effect, with the largest contribution at 300 K. A part of the discrepancy with our results may be due to using the re monomers as the reference point, 133 but the remaining differences are still large. For example, at 298.15 K our monomerV (re )

flexibility effect measured relative to Bcl

is only 12%. Since, as discussed below, our results

agree very well with experiment and with recent calculations by Babin et al. 119 (who have not, however, calculated the flexibility correction explicitly), we believe that the effects of monomer flexibility found by us must be close to true ones. In an independent calculation with the CCpol-8sfIR[2012] potential, one of us has found a 13% monomer-flexibility effect


Page 50 of 89

Page 51 of 89

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


at 300 K, 117 in agreement with our CCpol-8sfIR[2014] result in Table 9. It is worth to note that although the dependence of the CCpol-8sf/IR[2014] surface on the intramolecular degrees of freedom is inherited from SAPT-5s′ f/IR[2014], the flexibility hV i

correction ∆Bflex 0 calculated from these two surfaces is different due to the nonlinear depenhV i

dence of B on the potential via the Mayer function. The values of ∆Bflex 0 obtained from SAPT-5s′ f/IR[2014] in analogous way as the values in Table 9 are presented in Table SI-V of the Supporting Information. Table 10 and Fig. 3 compare our values of the second virial coefficient obtained from different potentials with the fit of B(T ) to multiple sets of experimental data, performed by Harvey and Lemmon. 35 These authors do not recommend the experimental fit to be used below 300 K, but we have include the values computed from this fit for all temperatures. The theoretical results for flexible monomers have been obtained as the sum defined by Eq. (42) hV i

TI calculated from the CCpol-8s and ∆Bflex 0 from the CCpol-8sf[2014] surface. This with BQM

sum will be denoted as Bflex [CCpol-8sf[2014]]. Column 2 contains the values of the rigidmonomer classical virial coefficient Bcl [CCpol-8s] computed using the same integration grid TI as for BQM , which enables a precise evaluation of the quantum effects. Panel (a) of Fig. 3

shows the total virial coefficients, whereas panels (b) and (c) show the differences from the experimental data of ref 35 for two different ranges of temperatures. In addition to the data from Table 10, Fig. 3 also presents the results of other recent theoretical calculations of the TI 117 second virial coefficient: B[HBB2-pol], 119 B[MB-pol], 85 and BQF .

The first observation from Table 10 and Fig. 3(a) is that theory and experiment achieved a remarkable degree of agreement for the second virial coefficient of water. This is particularly striking when comparing with analogous graphs published about 15 years ago, see, e.g., refs 48 and 54. The values of Bflex [CCpol-8sf[2014]] and the virtually identical values of B[HBB2-pol] lie overall closest to the experimental curve, with theory everywhere above experiment. At the lowest temperature where the fit to experimental data is valid, 323.15 K, the deviation between theory and experiment is about 1.7 times the experimental uncertainty.



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 discrepancy gets slightly larger (a factor of 2.7) at 373.15 K, and then gradually decreases, with theory becoming consistent with the experimental fit above 1800 K. However, for T > 1000 K, our present results probably provide the best estimate of the second virial coefficient for water. The reason is that the set of fitted data in Harvey and Lemmon work included for T > 700 K the values of B(T ) calculated from the SAPT-5s potential of ref 54. The Harvey-Lemmon B(T ) and SAPT-5s B(T ) differ by no more than 0.1 cm3 /mol above 873 K (see Table III of ref 56). Obviously, the virial coefficients computed from the present potential are superior to those computed from SAPT-5s. The strikingly good agreement between Bflex [CCpol-8sf[2014]] and B[HBB2-pol] is perhaps unexpected. The latter results have been obtained from a different potential energy surface, HBB2-pol, 18 using the PIMC method to calculate the second virial coefficient. This method does not use the three approximations assumed by us: the TI approximation for quantum effects, the neglect of quantum effects in the monomer-flexibility correction, and accounting for the flexibility effects via averaging of the potential. Apparently, all these effects are either small or cancel among themselves. The B[MB-pol] results 85,134 are included in only one of panels of Fig. 3 and not in Table 10 because of the narrow range and a different set of T ’s used by these authors (notice also that these values do not lie on a smooth curve). However, the differences between B[MB-pol], B[HBB2-pol], and Bflex [CCpol-8sf[2014]] are small compared to discrepancies of the virials from either potential relative to experiment. TI The values of BQF computed in ref 117 using the CCpol-8sfIR[2012] potential are in

excellent agreement with our values computed from the CCpol-8sfIR[2014] potential except for a 2.7% discrepancy at 323 K, the smallest temperature considered in ref 117. Since the former values were computed at a more advanced level of theory, this agreement confirms that the approximations made by us do not affect our predictions in a significant way, as already concluded from comparisons with results of ref 119. The approach used in ref 117 included monomer-flexibility effects at the adiabatic approximation level, i.e., the averaging


Page 52 of 89

Page 53 of 89

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


over intramonomer coordinates was performed using vibrational wave functions computed in the field of the interacting partner. In addition, the TI approximation was applied using the averaged potential, in this way partially including the coupling between quantum and monomer-flexibility effects. One might attribute the discrepancy between our and ref 117 results at the low T to this coupling, but this is probably not the case taking into account the uniformly good agreement with the results of ref 119.

8

Summary and discussion

The present work extends and supplements the development of a 12-dimensional potential energy surface for the water dimer initiated in ref 79 and then improved in ref 82. We used the original set of SAPT ab initio interaction energies, 79 larger than the sets used in more recent works. 17,85 The fits of ref 79, SAPT-5s′ f and SAPT-5s′ fIR, were reoptimized after fixing a minor symmetry error in the subroutine computing the potentials, but the changes due to the reoptimizations were very small. The two final hybrid potentials are denoted as CCpol-8sf[2014] and CCpol-8sfIR[2014]. We have performed a detailed investigation of the accuracy of the fits and of its components by analyzing scatter plots and rmse values on the training set and its subset, comparing to the ab initio interaction energies from ref 85, and analyzing selected radial cross-sections through the surface. This analysis shows that CCpol-8sf[2014] is accurate to about 0.2 kcal/mol for negative interaction energies in its region of validity in monomer’s distortions. This level of accuracy is consistent with the accuracy of ab initio training energies and with our attempt to keep the functional form of the fit relatively compact. With this form, we were not able fit the training data set extended by the additional points near the global minimum and the Hessian values equally well as the original data set, so that the rmse of SAPT-5s′ fIR[2014] is about 25% (15%) larger than that of SAPT-5s′ f[2014] on subsets of our training set (on Babin et al. 85 IN grid points). Therefore, the CCpol-8sf[2014] potential is our recommended one for all applications except for those probing the potential surface only near its minimum (e.g., dimer spectra and 53 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

low-energy characteristic points on the surface). For the latter applications, we recommend CCpol-8sfIR[2014]. The reoptimized potentials were applied to find the characteristic points on the 12dimensional water dimer surface. At the minimum of the surface, the results obtained with CCpol-8sf[2014] and CCpol-8sfIR[2014] differ by only 0.002 and 0.006 kcal/mol, respectively, from recent very accurate CCSD(T) calculations by Lane, 106 an agreement to within 0.1% or better. The geometries of the dimer’s minimum are also in excellent agreement with Lane’s results. 106 In fact, the potentials give the parameters of the minimum closer to Lane’s benchmarks than do the previous best benchmarks of Tschumper et al. 107 The only parameter that agrees less well is the elongation of the hydrogen-bonded OH bond relative to the re value: CCpol-8sfIR[2014] gives 71% of Lane’s value. However, this elongation is rather small, only 0.0055 ˚ A in Lane’s calculations. The rigid-monomer CCpol-8s potential characterises the minimum also very well, showing that this feature of the surface is not very sensitive to monomer-flexibility effects. The heights of the barriers defined by the nine saddle points and the global minimum obtained from the CCpol-8sf[2014] and CCpol-8sfIR[2014] surfaces are very close to the limit values computed ab initio by Tschumper et al. 107 For the slightly better performing CCpol8sf[2014] potential, the rmse is only 6 cm−1 relative to Tschumper’s et al. results at the same level of theory, well within error bars estimated by Tschumper et al. CCpol-8s, the rigid-monomer parent potential, gives a significantly larger rmse, amounting to 47 cm−1 , indicating larger effects of monomer-flexibility for the barriers than for the minimum. The developed surfaces were applied in calculations of the second virial coefficients. First, we investigated the monomer-flexibility effects. To account for these effects, we first averaged the flexible-monomer potential over the ground state rovibrational wave functions of the monomers, and then performed the 6-dimensional integral of the Mayer function of this averaged potential to obtain classical virial coefficients. The resulting values have been compared to the ones obtained from two rigid potentials, V (re ) and V (hri0 ), corresponding


Page 54 of 89

Page 55 of 89

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


to the equilibrium and the vibrationally-averaged geometries of the monomers, respectively. V (re )

If the flexibility effect is measured with respect to Bcl V (hri0 )

by Bcl

, a large part of this effect is recovered

. Since the two calculations are equally difficult, the use of the latter potential is V (hri0 )

advantageous and we defined the monomer-flexibility effect relative to Bcl

. The values

of ∆Bflex computed with CCpol-8sf[2014] and CCpol-8sfIR[2014] are quite different and we recommend the use of the former potential due to its better overall performance in fitting the training data set. The relative contributions of monomer-flexibility effects to our final values of the viral coefficients increase from about 1% at the lowest temperature to 9% at the highest, and are much smaller than those found in the work of Donchev et al. 133 The monomer-flexibility effects are an order of magnitude smaller than quantum effects at low T , but become larger in magnitude for T larger than about 600 K. Clearly, the former effects should be included to reach an agreement with experiment to within a few percent. To include quantum effects, the flexibility correction from the calculations described above was added to the coefficient computed with the rigid-monomer CCpol-8s potential at the Takahashi-Imada quantum level. 115 Our final results are in an excellent agreement (the two curves virtually overlap) with the results obtained from the HBB2-pol potential in ref 119 using a significantly different approach. Our virial coefficients were also compared to the analysis of experimental data by Harvey and Lemmon. 35 For the lowest temperatures where the experimental fit is valid, the largest deviations between theory and experiment are 2.7 times the experimental uncertainty and the deviations decrease with T , to reach near consistency at 973.15 K (the highest T for which experimental uncertainties were determined). For temperatures higher than 1000 K, our values may be actually more accurate than those of ref 35 since the experimental fit involved older, less accurate theoretical results of ref 54 in this range. Overall, monomer-flexibility effects were found to be only mildly important for determining the properties of water dimer except for calculations of the interaction-induced shifts of intramonomer vibrational transitions which require one to use a flexible-monomer poten-



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

tial. 82 For other properties, monomer-flexibility effects were found to be largest for virial coefficients at large temperatures, where these effects contribute about 10%. For low-T virials, intermolecular spectral transitions, and for the investigated characteristic points on the surface, the effects are generally amounting to only a couple percent. Nevertheless, at some point these effects have to be included to reach agreement with precise experiments. This relative smallness of monomer-flexibility effects is good news since, due to exponential scaling with the number of degrees of freedom, calculations of a flexible-monomer potential are orders of magnitude more difficult than calculations of a rigid-monomer potential at the same level of ab initio theory and are simply not practical for systems larger than the water dimer. Gora et al. 19 have recently shown that one can get better predictions of water cluster energetics with an accurate rigid-monomer force field than with less accurate flexible-monomer ones. All the investigations in the present work were limited to dimer properties. Our potentials can, however, be used in many-body applications, for example for water clusters and in simulations of liquid water. The explicit induction component of our potentials can be used to approximately model many-body polarization effects. The CCpol-8s potential was recently applied with such induction effects in simulations of liquid water, 26 and we expect CCpol-8sf[2014] to performed equally well or better. The next step would be to add to CCpol-8sf[2014] the complete three-body potential developed in ref 19.

Supporting Information Available Supporting Information contains a description of embedding procedures, presents results of calculations with atom-following and k-flex methods, explains our choices of grid points, analyzes performance of the fits on a range of diagrams, and gives numerical details of calculations of virial coefficients. Also, FORTRAN codes computing the potential are included. This material is available free of charge via the Internet at http://pubs.acs.org/.


Page 56 of 89

Page 57 of 89

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


Acknowledgement GM thanks David J. Nowak for invaluable computer assistance and advice and Alston J. Misquitta for helpful discussions. We thank Jonathan Tennyson for making available the PJT2 potential and the DVR3DJZ program and for advice on the use of this program. This work was supported by the NSF grant CHE-1152899.

References (1) See several papers in Water: A Comprehensive Treatise, edited by F. Franks (Plenum, New York, 1972), vol. 1. (2) Guillot, B. A Reappraisal of What We Have Learnt During Three Decades of Computer Simulations on Water. J. Mol. Liquids 2002, 101, 219–260. (3) Ball, P. Life’s Matrix: A Biography of Water ; Farrar: New York, 1999. (4) Szalewicz, K.; Jeziorski, B. Symmetry-Adapted Double-Perturbation Analysis of Intramolecular Correlation Effects in Weak Intermolecular Interactions. Mol. Phys. 1979, 38, 191–208. (5) Hayes, I. C.; Stone, A. J. Matrix Elements Between Determinantal Wavefunctions of Non-Orthogonal Orbitals. Mol. Phys. 1984, 53, 69–82. (6) Jeziorski, B.; Moszy´ nski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes. Chem. Rev. 1994, 94, 1887–1930. (7) Szalewicz, K.; Bukowski, R.; Jeziorski, B. In Theory and Applications of Computational Chemistry: The First 40 Years. A Volume of Technical and Historical Perspectives; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier, Amsterdam, 2005; Chapter 33, pp 919–962.



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) Szalewicz, K. Symmetry-adapted Perturbation Theory of Intermolecular Forces. Wiley Interdisc. Rev.–Comp. Mol. Sci. 2012, 2, 254–272. (9) Wallqvist, A.; Mountain, R. D. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; Wiley-VCH: New York, 1999; Vol. 13; pp 183–247. (10) Szalewicz, K.; Leforestier, C.; van der Avoird, A. Towards Complete Understanding of Water by First-Principle Computational Approach. Chem. Phys. Lett. 2009, 482, 1–14. (11) Niesar, U.; Corongiu, G.; Clementi, E.; Kneller, G. R.; Bhattacharya, D. K. MolecularDynamics Simulations of Liquid Water Using the NCC Ab Initio Potential. J. Phys. Chem. 1990, 94, 7949–7956. (12) Mas, E. M.; Bukowski, R.; Szalewicz, K. Ab Initio Three-body Interactions for Water. I. Potential and Structure of Water Trimer. J. Chem. Phys. 2003, 118, 4386–4403. (13) Mas, E. M.; Bukowski, R.; Szalewicz, K. Ab Initio Three-Body Interactions for Water. II. Effects on Structure and Energetic of Liquid. J. Chem. Phys. 2003, 118, 4404–4413. (14) 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. (15) Wang, Y.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Full-Dimensional, Ab Initio Potential Energy and Dipole Moment Surfaces for Water. J. Chem. Phys. 2009, 131, 054511–(1–8). (16) 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–(1–12). (17) Wang, Y.; Hunag, 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 22-mer. J. Chem. Phys. 2011, 134, 094509–(1–12). 58 ACS Paragon Plus Environment

Page 58 of 89

Page 59 of 89

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


(18) Medders, G. R.; Babin, V.; Paesani, F. A Critical Assessment of Two-Body and ThreeBody Interactions in Water. J. Chem. Theory Comput. 2013, 9, 1103–1114. (19) Gora, 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–(1:20). (20) Babin, V.; Medders, G. R.; Paesani, F. Development of a “First Principles” Water Potential with Flexible Monomers. II: Trimer Potential Energy Surface, Third Virial Coefficient, and Small Clusters. J. Chem. Theory Comput. 2014, 10, 1599–1607. (21) Anderson, J. A.; Crager, K.; Fedoroff, L.; Tschumper, G. S. Anchoring the Potential Energy Surface of the Cyclic Water Trimer. J. Chem. Phys. 2004, 121, 11023–11029. (22) 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. (23) Gora, U.; Podeszwa, R.; Cencek, W.; Szalewicz, K. Interaction Energies of Large Clusters from Many-Body Expansion: Water Clusters. J. Chem. Phys. 2011, 135, 224102–(1:19). (24) Stillinger, F. H.; David, C. W. Polarization Model for Water and Its Ionic Dissociation Products. J. Chem. Phys. 1978, 69, 1473–1494. (25) Lybrand, T. P.; Kollman, P. A. Water-water and Water-Ion Potential Functions Including Terms for Many Body Effects. J. Chem. Phys. 1985, 83, 2923–2933. (26) Akin-Ojo, O.; Szalewicz, K. How Well Can Polarization Models of Pairwise Nonadditive Forces Describe Liquid Water? J. Chem. Phys. 2013, 138, 024316–(1–13). (27) Zwart, E.; ter Meulen, J. J.; Meerts, W. L. The Submillimeter Rotation Tunneling Spectrum of (D2 O)2 . J. Mol. Spectrosc. 1990, 173, 115–121.



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

(28) Fraser, G. T. (H2 O)2 : Spectroscopy, Structure and Dynamics. Int. Rev. Phys. Chem. 1991, 10, 189–206. (29) Fellers, R. S.; Leforestier, C.; Braly, L. B.; Brown, M. G.; Saykally, R. J. Spectroscopic Determination of the Water Pair Potential. Science 1999, 284, 945–948. (30) Braly, L. B.; Liu, K.; Brown, M. G.; Keutsch, F. N.; Fellers, R. S.; Saykally, R. J. Terahertz Laser Spectroscopy of the Water Dimer Intermolecular Vibrations II: (H2 O)2 . J. Chem. Phys. 2000, 112, 10314–10326. (31) Keutsch, F. N.; Saykally, R. J. Water clusters: Untangling the Mysteries of the Liquid, One Molecule at a Time. Proc. Natl. Acad. Sci. USA 2001, 98, 10533–10540. (32) Keutsch, F. N.; Braly, L. B.; Brown, M. G.; Harker, H. A.; Petersen, P. B.; Leforestier, C.; Saykally, R. J. Water Dimer Hydrogen Bond Stretch, Donor Torsion Overtone, and “In-Plane Bend” Vibrations. J. Chem. Phys. 2003, 119, 8927–8937. (33) Kell, G. S.; McLaurin, G. E.; Whalley, E. PVT Properties of Water VII. Vapour Densities of Light and Heavy Water from 150 to 500

◦

C. Proc. R. Soc. London, Ser.

A 1989, 425, 49–71. (34) Lide, D. R., Ed. CRC Handbook of Chemistry and Physics; CRC, Boca Raton, New York, London, Tokyo, 1995; pp 6.30, 6.31. (35) Harvey, A. H.; Lemmon, E. W. Correlation for the Second Virial Coefficient of Water. J. Phys. Chem. Ref. Data 2004, 33, 369–374. (36) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; pp 331–342. (37) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269.


Page 60 of 89

Page 61 of 89

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


(38) 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. (39) Mas, E. M.; Szalewicz, K.; Bukowski, R.; Jeziorski, B. Pair Potential for Water from Symmetry-Adapted Perturbation Theory. J. Chem. Phys. 1997, 107, 4207–4218. (40) Cieplak, P.; Kollman, P.; Lybrand, T. A New Water Potential Including Polarization: Application to Gasphase, Liquid, and Crystal Properties of Water. J. Chem. Phys. 1990, 92, 6755–6760. (41) Wikfeldt, K. T.; Batista, E. R.; Vilazc, F. D.; Jonsson, H. A Transferable H2 O Interaction Potential Based on a Single Center Multipole Expansion: SCME. Phys. Chem. Chem. Phys. 2013, 15, 16542–16556. (42) Wang, L.-P.; Head-Gordon, T.; Ponder, J. W.; Ren, P.; Chodera, J. D.; Eastman, P. K.; Martinez, T. J.; Pande, V. S. Systematic Improvement of a Classical Molecular Model of Water. J. Phys. Chem. B 2013, 117, 9956–9972. (43) Fellers, R. S.; Braly, L. B.; Saykally, R. J.; Leforestier, C. Fully Coupled SixDimensional Calculations of the Water Dimer Vibration-Rotation-Tunneling States with Split Wigner Pseudospectral Approach. II. Improvements and Tests of Additional Potentials. J. Chem. Phys. 1999, 110, 6306–6318. (44) Goldman, N.; Fellers, R. S.; Brown, M. G.; Braly, L. B.; Keoshian, C. J.; Leforestier, C.; Saykally, R. J. Spectroscopic Determination of the Water Dimer Intermolecular Potential Energy Surface. J. Chem. Phys. 2002, 116, 10148–10163. (45) Keutsch, F. N.; Goldman, N.; Harker, H. A.; Leforestier, C.; Saykally, R. J. Complete Characterization of the Water Dimer Vibrational Ground State and Testing the VRT(ASP-W)III, SAPT-5st, and VRT(MCY-5f) Surfaces. Mol. Phys. 2003, 101, 3477–3492. 61 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

(46) Matsuoka, O.; Clementi, E.; Yoshimine, M. CI Study of Water Dimer Potential Surface. J. Chem. Phys. 1976, 64, 1351–1361. (47) Millot, C.; Stone, A. J. Towards an Accurate Intermolecular Potential for Water. Mol. Phys. 1992, 77, 439–462. (48) Millot, C.; Soetens, J. C.; Costa, M. T. C. M.; Hodges, M. P.; Stone, A. J. Revised Anisotropic Site Potentials for the Water Dimer, and Calculated Properties. J. Phys. Chem. A 1998, 102, 754–770. (49) Hayes, I. C.; Stone, A. J. An Intermolecular Perturbation Theory for the Region of Moderate Overlap. Mol. Phys. 1984, 53, 83–105. (50) Stone, A. J. The Theory of Intermolecular Forces; Clarendon Press, Oxford, 1996. (51) Groenenboom, G. C.; Mas, E. M.; Bukowski, R.; Szalewicz, K.; Wormer, P. E. S.; van der Avoird, A. The Pair and Three-Body Potential of Water. Phys. Rev. Lett. 2000, 84, 4072–4075. (52) Groenenboom, G. C.; Wormer, P. E. S.; van der Avoird, A.; Mas, E. M.; Bukowski, R.; Szalewicz, K. Water Pair Potential of Near Spectroscopic Accuracy: II. VibrationRotation-Tunneling Levels of the Water Dimer. J. Chem. Phys. 2000, 113, 6702–6715. (53) Smit, M. J.; Groenenboom, G. C.; Wormer, P. E. S.; van der Avoird, A.; Bukowski, R.; Szalewicz, K. Vibrations, Tunneling, and Transition Dipole Moments in the Water Dimer. J. Phys. Chem. A 2001, 105, 6212–6225. (54) Mas, E. M.; Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; Wormer, P. E. S.; van der Avoird, A. Water Pair Potential of Near Spectroscopic Accuracy. I. Analysis of Potential Surface and Virial Coefficients. J. Chem. Phys. 2000, 113, 6687–6701. (55) Burnham, C. J.; Xantheas, S. S. Development of Transferable Interaction Models for


Page 62 of 89

Page 63 of 89

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


Water. III. Reparametrization of an All-Atom Polarizable Rigid Model (TTM2-R) From First Principles. J. Chem. Phys. 2002, 116, 1500–1510. (56) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Interaction Potential for Water Dimer From Symmetry-Adapted Perturbation Theory Based on Density Functional Description of Monomers. J. Chem. Phys. 2006, 125, 044301–(1– 8). (57) Misquitta, A. J.; Podeszwa, R.; Jeziorski, B.; Szalewicz, K. Intermolecular Potentials Based on Symmetry-Adapted Perturbation Theory Including Dispersion Energies from Time-Dependent Density Functional Calculations. J. Chem. Phys. 2005, 123, 214103– (1:14). (58) 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. (59) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Polarizable Interaction Potential for Water from Coupled Cluster Calculations. I. Analysis of Dimer Potential Energy Surface. J. Chem. Phys. 2008, 128, 094313–(1:15). (60) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Polarizable Interaction Potential for Water from Coupled Cluster Calculations. II. Applications to Dimer Spectra, Virial Coefficients, and Simulations of Liquid Water. J. Chem. Phys. 2008, 128, 094314–(1:20). (61) van der Avoird, A.; Szalewicz, K. Water Trimer Torsional Spectrum from Accurate Ab Initio and Semi-Empirical Potentials. J. Chem. Phys. 2008, 128, 014302–(1:8). (62) Cencek, W.; Szalewicz, K.; Leforestier, C.; van Harrevelt, R.; van der Avoird, A. An Accurate Analytic Representation of the Water Pair Potential. Phys. Chem. Chem. Phys. 2008, 10, 4716–4731.



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

(63) Murdachaew, G.; Szalewicz, K. Intermolecular Potentials with Flexible Monomers. Faraday Discuss. 2001, 118, 121–142. (64) Murdachaew, G.; Bukowski, R.; Szalewicz, K. Efficient Generation of FlexibleMonomer Intermolecular Potential Energy Surfaces. Phys. Rev. Lett. 2002, 88, 123202–(1:4). (65) Toukan, K.; Rahman, A. Molecular-Dynamics Study of Atomic Motions in Water. Phys. Rev. B 1985, 31, 2643–2648. (66) Mahoney, M. W.; Jorgensen, W. L. Quantum, Intramolecular Flexibility, and Polarizability Effects on the Reproduction of the Density Anomaly of Liquid Water by Simple Potential Functions. J. Chem. Phys. 2001, 115, 10758–10768. (67) Lie, G. C.; Clementi, E. Molecular-Dynamics Simulation of Liquid Water with an Ab Initio Flexible Water-Water Interaction Potential. Phys. Rev. A 1986, 33, 2679–2693. (68) Corongiu, G. Molecular-Dynamics Simulation for Liquid Water Using a Polarizable and Flexible Potential. Int. J. Quant. Chem. 1992, 42, 1209–1235. (69) Niesar, U.; Corongiu, G.; Huang, M.-J.; Dupuis, M.; Clementi, E. Preliminary Observations on a New Water-Water Potential. Int. J. Quantum Chem. Symp. 1989, 23, 421–443. (70) Leforestier, C.; Gatti, F.; Fellers, R.; Saykally, R. Determination of a Flexible (12D) Water Dimer Potential via Direct Inversion of Spectroscopic Data. J. Chem. Phys. 2002, 117, 8710–8722. (71) Polyansky, O. L.; Jensen, P.; Tennyson, J. The Potential Energy Surface of H16 2 O. J. Chem. Phys. 1996, 105, 6490–6497. (72) 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 64 ACS Paragon Plus Environment

Page 64 of 89

Page 65 of 89

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


Dependent Charges Derived From an Ab Initio Monomer Dipole Moment Surface. J. Chem. Phys. 2002, 116, 5115–5124. (73) 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–(1–11). (74) Burnham, C. J.; Anick, D. J.; Mankoo, P. K.; Reiter, G. F. The Vibrational Proton Potential in Bulk Liquid and Ice. J. Chem. Phys. 2008, 128, 154519–(1–20). (75) Kita, S.; Noda, K.; Inouye, H. Repulsive Potentials for Cl− –R and Br− –R (R=He, Ne, and Ar) Derived from Beam Experiments. J. Chem. Phys. 1976, 64, 3446–3449. (76) Price, S. L.; Wheatley, R. J. An Overlap Model for Estimating the Anisotropy of Repulsion. Mol. Phys. 1990, 69, 507–533. (77) Hodges, M. P.; Wheatley, R. J. Application of the Overlap Model to Calculating Correlated Exchange Energies. Chem. Phys. Lett. 2000, 326, 263–268. (78) Jeziorska, M.; Jankowski, P.; Szalewicz, K.; Jeziorski, B. On the Optimal Choice of Monomer Geometry in Calculations of Intermolecular Potentials. Rovibrational Spectrum of Ar–HF Generated from Two- and Three-Dimensional SAPT Potentials. J. Chem. Phys. 2000, 113, 2957–2968. (79) Szalewicz, K.; Murdachaew, G.; Bukowski, R.; Akin-Ojo, O.; Leforestier, C. In Lecture Series on Computer and Computational Science: ICCMSE 2006 ; Maroulis, G., Simos, T., Eds.; Brill Academic Publishers, Leiden, 2006; Vol. 6; pp 482–491. (80) Huang, X.; Braams, B. J.; Bowman, J. M. Ab Initio Potential Energy and Dipole Moment Surfaces of (H2 O)2 . J. Phys. Chem. A 2006, 110, 445–451. (81) Murdachaew, G. Ph.D. Thesis, University of Delaware, 2005. 65 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

(82) Leforestier, C.; Szalewicz, K.; van der Avoird, A. Spectra of Water Dimer from a New Ab Initio Potential with Flexible Monomers. J. Chem. Phys. 2012, 137, 014305–(1:17). (83) Shank, A.; Wang, Y.; Kaledin, A.; Braams, B. J.; Bowman, J. M. Accurate Ab Initio and “Hybrid” Potential Energy Surfaces, Intramolecular Vibrational Energies, and Classical IR Spectrum of the Water Dimer. J. Chem. Phys. 2009, 130, 144314–(1– 11). (84) Leforestier, C. Infrared Shifts of the Water Dimer from the Fully Flexible Ab Initio HBB2 Potential. Phil. Trans. R. Soc. A 2012, 370, 2675–2690. (85) Babin, V.; Leforestier, C.; Paesani, F. Development of a “First Principles” Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient. J. Chem. Theory Comput. 2013, 9, 5395–5403. (86) Huang, X.; Braams, B. J.; Bowman, J. M.; Kelly, R. E. A.; Tennyson, J.; Groenenboom, G. C.; van der Avoird, A. New Ab Initio Potential Energy Surface and the Vibration-Rotation-Tunneling Levels of (H2 O)2 and (D2 O)2 . J. Chem. Phys. 2008, 128, 034312–(1–9). (87) Dunning, Jr., T. H. Gaussian-Basis Sets for Use in Correlated Molecular Calculations. 1. The Atoms Boron Through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007– 1023. (88) Kendall, R. A.; Dunning, Jr., T. H.; Harrison, R. J. Electron-Affinities of the 1st-Row Atoms Revisited — Systematic Basis-Sets and Wave-Functions. J. Chem. Phys. 1992, 96, 6796–6806. (89) Xantheas, S. S.; Burnham, C. J.; Harrison, R. J. Development of Transferable Interaction Models for Water. II. Accurate Energetics of the First Few Water Clusters from First Principles. J. Chem. Phys. 2002, 116, 1493–1499.


Page 66 of 89

Page 67 of 89

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


(90) Brink, D. M.; Satchler, G. R. Angular Momentum, 3rd ed.; Clarendon: Oxford, 1993. (91) Bernath, P. F. Spectra of Atoms and Molecules; Oxford: New York, 1995; p 228. (92) Mas, E. M.; Szalewicz, K. Effects of Monomer Geometry and Basis Set Saturation on Depth of Water Dimer Potential. J. Chem. Phys. 1996, 104, 7606–7614. (93) Partridge, H.; Schwenke, D. W. The Determination of an Accurate Isotope Dependent Potential Energy Surface for Water from Extensive Ab Initio Calculations and Experimental Data. J. Chem. Phys. 1997, 106, 4618–4639. (94) Polyansky, O. L.; Császár, A. G.; Shirin, S. V.; Zobov, N. F.; Barletta, P.; Tennyson, J.; Schwenke, D. W.; Knowles, P. J. High-Accuracy Ab Initio Rotation-Vibration Transitions for Water. Science 2003, 299, 539–542. (95) Shirin, S. V.; Polyansky, O. L.; Zobov, N. F.; Barletta, P.; Tennyson, J. Spectroscopically Determined Potential Energy Surface of H2 16 O up to 25000 cm−1 . J. Chem. Phys. 2003, 118, 2124–2129. (96) Coheur, P. F.; Bernath, P. F.; Carleer, M.; Colin, R.; Polyansky, O. L.; Zobov, N. F.; Shirin, S. V.; Barber, R. J.; Tennyson, J. A 3000 K Laboratory Emission Spectrum of Water. J. Chem. Phys. 2005, 122, 074307–(1–8). (97) Wei, H.; Carrington, T. Explicit Expressions for Triatomic Eckart Frames in Jacobi, Radau, and Bond Coordinates. J. Chem. Phys. 1997, 107, 2813–2818. (98) Wei, H.; Carrington, T. An Exact Eckart-Embedded Kinetic Energy Operator in Radau Coordinates for Triatomic Molecules. Chem. Phys. Lett. 1998, 287, 289–300. (99) Bunker, P. R.; Jensen, P. Molecular Symmetry and Spectroscopy, 2nd ed.; NRC: Ottawa, 1998.



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

(100) Tang, K. T.; Toennies, J. P. An Improved Simple-Model for the van der Waals Potential Based on Universal Damping Functions for the Dispersion Coefficients. J. Chem. Phys. 1984, 80, 3726–3741. (101) Wormer, P. E. S.; Hettema, H. Many-Body Perturbation Theory of FrequencyDependent Polarizabilities and van der Waals Coefficients: Application to H2 O–H2 O and Ar–NH3 . J. Chem. Phys. 1992, 97, 5592–5606. (102) Wormer, P. E. S.; Hettema, H. POLCOR Package, University of Nijmegen. 1992. (103) Golub, G. H.; Reinsch, C. Singular Value Decomposition and Least Squares Solutions. Num. Math. 1970, 14, 403–420. (104) Garberoglio, G.; Jankowski, P.; Szalewicz, K.; Harvey, A. H. Second Virial Coefficients of H2 and Its Isotopologues from Six-Dimensional Potential. J. Chem. Phys. 2012, 137, 154308–(1–11). (105) Babin, V.; Leforestier, C.; Paesani, F. Erratum: Development of a “First Principles” Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient. J. Chem. Theory Comput. 2014, 10, 3585–3585. (106) Lane, J. R. CCSDTQ Optimized Geometry of Water Dimer. J. Chem. Theory Comput. 2013, 9, 316–323. (107) Tschumper, G. S.; Leininger, M. L.; Hoffman, B. C.; Valeev, E. F.; Schaefer III, H. F.; Quack, M. Anchoring the Water Dimer Potential Energy Surface with Explicitly Correlated Computations and Focal Point Analyses. J. Chem. Phys. 2002, 116, 690–701. (108) Odutola, J. A.; Dyke, T. R. Partially Deuterated Water Dimers: Microwave Spectra and Structure. J. Chem. Phys. 1980, 72, 5062–5070. (109) Moszynski, R.; Korona, T.; Heijmen, T. G. A.; Wormer, P. E. S.; van der Avoird, A.;


Page 68 of 89

Page 69 of 89

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


Schramm, B. Second Virial Coefficients for Atom-Molecule Complexes from Ab Initio SAPT Potentials. Polish J. Chem. 1998, 72, 1479–1496. (110) Garberoglio, G.; Jankowski, P.; Szalewicz, K.; Harvey, A. H. Path-Integral Calculation of the Second Virial Coefficient Including Intramolecular Flexibility Effects. J. Chem. Phys. 2014, 141, 044119–(1–15). (111) Refson, K.; Lie, G. C.; Clementi, E. The Effect of Molecular Vibrations on Calculated Second Virial Coefficients. J. Chem. Phys. 1987, 87, 3634–3638. (112) MacDowell, L. G.; Vega, C. The Second Virial Coefficient of Hard Alkane Models. J. Chem. Phys. 1998, 109, 5670–5680. (113) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids Volume 1: Fundamentals; Clarendon Press, Oxford, 1984. (114) Wormer, P. E. S. Second Virial Coefficients of Asymmetric Top Molecules. J. Chem. Phys. 2005, 122, 184301–(1–7). (115) Takahashi, M.; Imada, M. Monte Carlo Calculation of Quantum Systems. II. Higher Order Correction. J. Phys. Soc. Jpn. 1984, 53, 3765–3769. (116) Schenter, G. K. The Development of Effective Classical Potentials and the Quantum Statistical Mechanical Second Virial Coefficient of Water. J. Chem. Phys. 2002, 117, 6573–6581. (117) Leforestier, C. Water Dimer Equilibrium Constant Calculation: A Quantum Formulation Including Metastable States. J. Chem. Phys. 2014, 140, 074106–(1–9). (118) Patkowski, K.; Cencek, W.; Jankowski, P.; Szalewicz, K.; Mehl, J. B.; Garberoglio, G.; Harvey, A. H. Potential Energy Surface for Interactions Between Two Hydrogen Molecules. J. Chem. Phys. 2008, 129, 094304–(1:19).



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

(119) 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. (120) McQuarrie, D. A. Statistical Mechanics; Harper: New York, 1976. (121) Hill, T. L. An Introduction to Statistical Thermodynamics; Dover, 1986. (122) Shaul, K. R. S.; Schultz, A. J.; Kofke, D. A. Mayer-Sampling Monte Carlo Calculations of Uniquely Flexible Contributions to Virial Coefficients. J. Chem. Phys. 2011, 135, 124101–(1–9). (123) Martin, J. M. L.; Franois, J. P.; Gijbels, R. On the Effect of Centrifugal Stretching on the Rotational Partition Function of an Asymmetric Top. J. Chem. Phys. 1991, 95, 8374–8389. (124) Tennyson, J.; Kostin, M. A.; Barletta, P.; Harris, G. J.; Polyansky, O. L.; Ramanlal, J.; Zobov, N. F. DVR3D: A Program Suite for the Calculation of Rotation Vibration Spectra of Triatomic Molecules. Comput. Phys. Commun. 2003, 163 . (125) DVR3D: A Program Suite for the Calculation of Rotation Vibration Spectra of Triatomic Tolecules by J. Tennyson, M. A. Kostin, P. Barletta, G. J. Harris, O. L. Polyansky, J. Ramanlal, and N. F. Zobov, CPC Program Library, http://cpc.cs.qub.ac.uk/summaries/ADTI. (126) Jankowski, P. Approximate Generation of Full-Dimensional Ab Initio van der Waals Surfaces for High-Resolution Spectroscopy. J. Chem. Phys. 2004, 121, 1655–1662. (127) Jankowski, P.; Szalewicz, K. Effects of Monomer Flexibility on Spectra of N2 –HF. Chem. Phys. Lett. 2008, 459, 60–64.


Page 70 of 89

Page 71 of 89

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


(128) Jankowski, P.; Szalewicz, K. A New Ab Initio Interaction Energy Surface and HighResolution Spectra of the H2 -CO van der Waals Complex. J. Chem. Phys. 2005, 123, 104301–(1–16). (129) Jankowski, P.; McKellar, A. R. W.; Szalewicz, K. Theory Untangles High-Resolution Infrared Spectrum of the orthoH2 -CO van der Waals Complex. Science 2012, 336, 1147–1150. (130) Jankowski, P.; Surin, L. A.; Potapov, A.; Schlemmer, S.; McKellar, A. R. W.; Szalewicz, K. A Comprehensive Experimental and Theoretical Study of H2 -CO Spectra. J. Chem. Phys. 2013, 138, 084307–(1:23). (131) Conroy, H. Molecular Schrödinger Equation. VIII. A New Method for the Evaluation of Multidimensional Integrals. J. Chem. Phys. 1967, 47, 5307–5318. (132) Jankowski, P. Exploring the New Three-Dimensional Ab Initio Interaction Energy Surface of the Ar-HF Complex: Rovibrational Calculations for Ar-HF and Ar-DF with Vibrationally Excited Diatoms. J. Chem. Phys. 2008, 128, 154311–(1–11). (133) Donchev, A. G.; Galkin, N. G.; Tarasov, V. I. Contribution of Quantum Molecular Flexibility to the Second Virial Coefficient of Water Vapor. Phys. Rev. Lett. 2006, 97, 220401–(1–4). (134) Babin, V.; Paesani, F. private communication.



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 72 of 89

Table 1: Comparisons of quality of fits. The unweighted root mean square errors (rmse) for various quantities and potentials are shown together with the number of data points contributing to a given rmse. Units of rmse of interaction energies are kcal/mol, whereas atomic units are used for tesseral multipole moments Qlk , l = 0–8 (e·bohrl ) and for isotropic dipole-dipole polarizabilities (bohr3 ). We have listed the numbers of points for individual tesseral components, but the training set always contained all points, including also those with l = 4-8.

All points SAPT-5s (ref 54) fit of charges to multipole moments — dipole tesseral components Q1k — quadrupole tesseral components Q2k — octupole tesseral components Q3k fit of induction plus dispersion asymptotics SAPT-5s

rmse

25 5.09 1 0.0000003 2 0.036 2 0.85 2510 0.0020 2510 0.37

SAPT-5s′ (only elements different from SAPT-5s) fit of induction plus dispersion asymptotics 2510 SAPT-5s′ 2510 SAPT-5s′ f fit of charges to multipole moments 470 20 — dipole tesseral components Q1k — quadrupole tesseral components Q2k 35 — octupole tesseral components Q3k 41 fit of isotropic dipole-dipole polarizability 14 fit of induction plus dispersion asymptotics 562,240 239,928 SAPT-5s′ f[2014] SAPT-5s′ fIR[2014] 239,928e

0.0013 0.39

Eint < 150 points rmse

Eint < 10 Eint < 0 number of points rmse points rmse parametersa 3a

2510

0.37

2083

0.18

1045

0.10

9 65b,c

2510

0.39

2083

0.18

1045

0.10

9d 66b,c

6.66 0.00086 0.025 0.76 0.0068 0.0015 3.06 239,598 3.11 239,598

20a

0.82 167,334 1.00 167,334

a

0.43 73,846 0.54 73,846

0.24 0.30

5 63d 568b,c 568b,c

Not including the position of site D1 which was optimized in ref 54 as a preliminary part of the fitting process. b As in footnote a for sites D2 and D3. c Not including the asymptotic parameters. d Not including the partial charges, the position of site D1, and the isotropic dipole-dipole polarizability. e The training set of SAPT-5s′ f[2014], additional interaction energies used to fit SAPT-5s′ fIR[2014] not included.


Page 73 of 89

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: List of the parameters defining the SAPT-5s′ f surface. The numbers of parameters refer to symmetry-distinct ones and show only the parameters actually optimized at a given level, i.e., the parameters from rigid-monomer fits kept fixed in flexible-monomer fits are not counted in the latter case. parameter

formula

qa, qb q a (sA ), q b (sB ) α ¯A, α ¯B A α ¯ (sA ), α ¯ B (sB ) Cnab , n = 6, 8, 10 Cnab (sA , sB )

(3) (5), (7), (8) (3) (5), (9) (3) (5), (10)

δ1ab δnab , n = 6, 8, 10 δ6d

(3) (3) (3)

aab n ab an (s)

(3) (5), (11) (5) (5)

αab βab

number type when optimized long range 3 linear rigid 20 linear flexible 1 linear rigid 6 linear flexible 9 linear rigid 63 linear flexible damping parameters 6 nonlinear rigid 9 nonlinear rigid 1 nonlinear rigid short range 30 linear rigid 568 linear flexible 10 nonlinear rigid 10 nonlinear rigid



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 3: The family of the interaction energy surfaces investigated in the present paper. The [2012] SAPT potentials differ from the [2006] ones by fixing a symmetry error in the potential subroutines without reoptimizing the parameters. Such reoptimization was performed for the [2014] potentials. This optimization is the reason for the [2014] potentials being different from [2012] ones. surface SAPT-5s′ SAPT-5s′ f[2006] SAPT-5s′ fIR[2006] SAPT-5s′ fIR[2012] SAPT-5s′ f[2014] SAPT-5s′ fIR[2014] CCpol-8sfIR[2012]a CCpol-8sf[2014] CCpol-8sfIR[2014] a

origin ref 79 ref 79 ref 79 ref 82 present present ref 82 present present

Denoted as CCpol-8sf in ref 82.


Page 74 of 89

Page 75 of 89

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 4: The list of the symmetry-adapted functions used to define the site-site components ± V IJ of Eq. (13). The functions h± i are defined as hi = (1 ± PAB )Si , i = 0, . . . , 10, where − − − 2 2 2 2 2 2 S = [1, s1 , s2 , s1 , s2 , s1 s2 , s1 s4 , s2 s5 , s3 , s2 s3 , s3 s6 ]. Note the h− 0 = h6 = h7 = h10 = 0. We drop the arguments of the β IJ (sA , sB ) functions for brevity. V IJ V OO

symmetry-adapted basis functions ψiIJ OO h+ ; ROA OB ) i w, i = 0, · · · , 10; w = w(β

V D3D3

like V OO with w = w(β D3D3 ; RD3A D3B )

V OD3

h+ i (1 + PAB )w, i = 0, · · · , 10; h− i (1 − PAB )w, i = 0, · · · , 10; w = w(β OD3 ; ROA D3B )

V OH

B h+ i (1 + PAB )(1 + P12 )w, i = 0, . . . , 7 − B hi (1 − PAB )(1 + P12 )w, i = 1, . . . , 5 i B (1 + PAB )[s6 (1 + (−1)i P12 )w], i = 1, 2 i i B (1 + PAB )[s5 s6 (1 + (−1) P12 )w], i = 1, 2 B (1 + PAB )[s23 s6 (1 − P12 )w] OH w = w(β ; ROA H1B )

V D3H

like V OH with w = w(β D3H ; RD3A H1B )

V HH

A B h+ i (1 + P12 )(1 + P12 )w, i = 0, . . . , 7 A B (1 + PAB )[s6 (1 + P12 )(1 − P12 )w] A B (1 + PAB )[s2 s6 (1 + P12 )(1 − P12 )w] A B (1 + PAB )[s3 s6 (1 − P12 )(1 − P12 )w] HH w = w(β ; RH1A H1B )

V OD2

B like V OD3 with w = (1 + P12 )w(β OD2 ; ROA D21B )

V D3D2

B like V OD3 with w = (1 + P12 )w(β D3D2 ; RD3A D21B )

V D2D2

A B like V OO with w = (1 + P12 )(1 + P12 )w(β D2D2 ; RD21A D21B )

V D2H

A like V OH with w = (1 + P12 )w(β D2H ; RD21A H1B )



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: Comparison of rmse values (in kcal/mol) for the interaction energies calculated with various versions of the CCpol-8sf potential on the training set of points from ref. 85, 105. These rmse’s were evaluated with respect to the interaction energies computed at the CCSD(T)/CBS level. 85,105 The values denoted by ALL have been obtained for all 42394 geometries. For 21714 of these geometries, the monomers were distorted within the classical turning points corresponding to the ground rovibrational state of the water molecule and the results obtained for this set are denoted by IN. If at least one of the molecules is more distorted, such dimer geometry belongs to the set denoted here by OUT. CCpol-8sfIR[2012] CCpol-8sfIR[2014] CCpol-8sf[2014] IN 0.1934 0.1933 0.1672 OUT 0.5721 0.5756 0.6100 ALL 0.4229 0.4251 0.4425


Page 76 of 89

Page 77 of 89

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 6: Comparison of harmonic frequency shifts of intramolecular modes in (H2 O)2 obtained with the CCpol-8sfIR[2012] and CCpol-8sfIR[2014] (in cm−1 ). Mode b[D] O-Hb [D] O-Hf [D] b[A] ss[A] as[A] a

CCpol-8sfIR[2012]a (Radau f = 1) +18.8 −48.3 −26.8 +4.2 −3.5 −9.3

CCpol-8sfIR[2012] (Radau f = 1) +18.50 −48.32 −26.38 +4.25 −3.51 −9.30

CCpol-8sfIR[2014] (Radau f = 1) +18.04 −50.63 −25.96 +4.05 −4.56 −9.56

Results from Table VII of ref 82.


CCpol-8sfIR[2014] (Eckart) +18.03 −50.63 −25.95 +4.05 −4.56 −9.56


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 78 of 89

Table 7: Global minimum on the water dimer potential. The [2014] versions of the SAPT5s′ f, CCpol-8sf, and CCpol-8sfIR surfaces were used and the PJT2 71 potential represented intramonomer parts. For rigid-monomer potentials, the intramonomer coordinates listed are fixed. For flexible-monomer potentials, Vmin is the vertical interaction energy at the geometry corresponding to Umin . All lengths are in ˚ A, angles in degrees, and energies in kcal/mol. The angle β is the angle formed by the OO vector and the bisector of the angle H3 O2 H4 , while α is the angle H2 O1 O2 (the hydrogen being donated is H2 ). α β ROO r O1 H 1 r O1 H 2 θA UA g r O2 H 3 r O2 H 4 θB UB g Vmin Umin g

SAPT-5sa SAPT-5s′ SAPT-5s′ f(re ) SAPT-5s′ f CCpol-8sb CCpol-8sf CCpol-8sfIR Tschumper et al.c Laned Expt./theory 6.3 5.9 5.5 5.6 6.0 5.7 5.5 4.74 5.686 2 ± 4e 127 128.3 129.1 128.9 121.88 122.5 123.5 124.92 123.458 123 ± 10e 2.955 2.9561 2.9623 2.9606 2.9105 2.9145 2.9118 2.9089 2.9092 2.91 ± 0.005f 0.9716 0.9716 0.9579 0.9596 0.9716 0.9596 0.9575 0.9581 0.9569 — 0.9716 0.9716 0.9579 0.9606 0.9716 0.9607 0.9618 0.9653 0.9641 — 104.69 104.69 104.50 104.77 104.69 104.72 104.71 104.45 104.854 — 0.0 0.0078 0.0076 0.0103 0.0299h 0.0271 — 0.9716 0.9716 0.9579 0.9583 0.9716 0.9584 0.9590 0.9597 0.9584 — 0.9716 0.9716 0.9579 0.9583 0.9716 0.9584 0.9590 0.9597 0.9584 — 104.69 104.69 104.50 105.00 104.69 105.05 104.67 104.58 104.945 — 0.0 0.0043 0.0053 0.0020 0.0016h 0.0037 — -4.861 -4.8484 -4.7450 -4.7693 -5.1040 -5.0211 -5.0285 −5.012i -5.0516i — -4.7450 -4.7571 -5.0081 -5.0163 −4.981j -5.0208 —

a

Calculated by Mas et al. in ref 54 with the rigid-monomer SAPT-5s surface. Calculated by Cencek et al. in ref 62 with the rigid-monomer CCpol-8s surface. c Tschumper et al., ref 107, using the CCSD(T) method and the TZ2P(d, f )+dif basis set. Unless otherwise stated, the zero of energy for this column is the energy of the monomer(s) at the monomer geometry of ref 107: bond length of 0.9589 ˚ A and angle of 104.16◦ . The values of minimum energies were computed at the geometries given using counterpoise-corrected CBS extrapolations from up to aug-cc-pV6Z bases. d Lane, ref 106, best estimate obtained at the counterpoise-corrected CCSD(T)-F12b/CBS limit with corrections for higher order excitations, core correlation, relativistic effects, and diagonal Born-Oppenheimer effects. The value of Umin at the CCSD(T) level is -5.0103 kcal/mol. e Odutola and Dyke, ref 108. The α angle quoted in ref 54 was incorrectly computed from the angles given in ref 108. f Mixed experimental/theoretical result as given in ref 54, based on experimental values from ref 108. g UX is the energy of monomer X at the geometry listed in the preceding rows computed from the PJT2 potential 71 unless otherwise stated. The zero value of UX is at PJT2 equilibrium value. h Calculated by us using Tschumper’s et al. 107 monomer geometries and the PJT2 potential. i Obtained by subtracting UA + UB calculated as in footnote g from Umin . j Counterpoise-corrected De from ref 107. b


Page 79 of 89

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 8: Energies (in cm−1 ) at stationary points relative to the potential minimum (barriers) calculated with the CCpol-8sf[2014] and CCpol-8sfIR[2014] potentials using the geometries optimized for individual potentials and compared to the barriers obtained by Tschumper et al. 107 Number 2 3 4 5 6 7 8 9 10 rmsed

CCpol-8sa 190.5 212.6 250.0 342.4 368.2 702.4 1293.6 701.0 1027.0 46.6

CCpol-8sfIR CCpol-8sf 175.6 185.2 198.6 206.1 245.4 244.8 322.0 333.2 343.4 357.1 c 646.2 647.2c 1238.5c 1253.5c 642.8 636.9 974.0 966.5 12.6 9.7

a

Tschumper et al.b 181 ± 5 198 ± 6 245 ± 15 333 ± 15 348 ± 18 634 ± 18 1249 ± 21 625 ± 16 948 ± 23

ref 62, rigid-monomer potential. ref 107, extrapolated nonrelativistic values with core effects included. c We were not able to converge our optimization procedure for these stationary points. Instead, the values calculated for the geometries optimized in ref 107 are given. d With respect to the Tschumper’s et al. results. b



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 9: Monomer flexibility effects in the classical second virial coefficient Bcl (T ) (in cm3 /mol) as functions of temperature (in K) for the CCpol-8sf[2014] and CCpol-8sfIR[2014] potentials. Definitions of the symbols used are given in Section 7.4. The flexibility correcV (hri ) hV i hV i tion to the virial coefficient is calculated from the formula ∆Bflex 0 = Bcl 0 − Bcl 0 . The calculations have been performed using the method of integration of Conroy. 131

T [K] 273.15 293.15 295.15 298.15 323.15 373.15 423.15 448.15 473.15 523.15 573.15 673.15 773.15 873.15 973.15 1000.00 1100.00 1200.00 1300.00 1400.00 1500.00 1600.00 1800.00 2000.00 2500.00 3000.00

CCpol-8s V (hri ) Bcl 0 -2236.422 -1494.340 -1440.404 -1364.614 -909.826 -484.254 -300.574 -245.915 -204.955 -148.522 -112.151 -69.113 -45.164 -30.247 -20.226 -18.086 -11.557 -6.690 -2.953 -0.017 2.336 4.252 7.152 9.211 12.318 13.940

CCpol-8sf[2014] hV i V (r ) ∆Bflex 0 Bcl e -2011.697 -21.394 -1359.576 -25.174 -1311.843 -25.234 -1244.676 -25.270 -838.962 -23.808 -453.480 -18.828 -284.269 -14.918 -233.412 -13.431 -195.108 -12.183 -142.016 -10.248 -107.565 -8.821 -66.516 -6.895 -43.513 -5.656 -29.120 -4.801 -19.417 -4.172 -17.341 -4.026 -11.000 -3.579 -6.265 -3.222 -2.624 -2.927 0.241 -2.687 2.538 -2.472 4.410 -2.304 7.248 -2.018 9.263 -1.795 12.310 -1.408 13.902 -1.158


CCpol-8sfIR[2014] hV i V (r ) ∆Bflex 0 Bcl e -2042.041 -249.812 -1377.710 -178.298 -1329.126 -173.008 -1260.775 -165.538 -848.258 -119.860 -457.117 -74.694 -285.841 -53.616 -234.445 -46.953 -195.770 -41.787 -142.223 -34.304 -107.526 -29.175 -66.253 -22.575 -43.172 -18.490 -28.754 -15.705 -19.048 -13.676 -16.975 -13.222 -10.643 -11.781 -5.920 -10.618 -2.294 -9.669 0.557 -8.887 2.841 -8.216 4.700 -7.656 7.513 -6.734 9.508 -5.999 12.514 -4.728 14.076 -3.908

Page 80 of 89

Page 81 of 89

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 10: Total second virial coefficient Bflex (in cm3 /mol) obtained by adding the monomer flexibility correction ∆Bflex computed with the CCpol-8sf[2014] potential to the rigidTI monomer BQM results computed with the CCpol-8s potential of ref 62. The values of Bcl calculated with CCpol-8s are also given to show the quantum effect on the second virial coefficient.

T [K] 273.15 293.15 295.15 298.15 323.15 373.15 423.15 448.15 473.15 523.15 573.15 673.15 773.15 873.15 973.15 1000.00 1100.00 1200.00 1300.00 1400.00 1500.00 1600.00 1800.00 2000.00 2500.00 3000.00

Bcl a -2257.580 -1506.552 -1452.009 -1375.381 -915.901 -486.650 -301.716 -246.741 -205.569 -148.880 -112.373 -69.207 -45.204 -30.260 -20.226 -18.083 -11.547 -6.678 -2.941 -0.002 2.351 4.268 7.168 9.224 12.320 13.922

hV i

TI a BQM ∆Bflex 0 b -1585.518 -21.394 -1115.594 -25.174 -1080.225 -25.234 -1030.191 -25.270 -720.153 -23.808 -408.522 -18.828 -263.561 -14.918 -218.588 -13.431 -184.195 -12.183 -135.683 -10.248 -103.610 -8.821 -64.712 -6.895 -42.552 -5.656 -28.547 -4.801 -19.040 -4.172 -17.000 -4.026 -10.750 -3.579 -6.070 -3.222 -2.461 -2.927 0.382 -2.687 2.666 -2.472 4.530 -2.304 7.357 -2.018 9.367 -1.795 12.399 -1.408 13.973 -1.158

Bflex c -1606.912 -1140.768 -1105.459 -1055.461 -743.961 -427.350 -278.479 -232.019 -196.378 -145.931 -112.431 -71.607 -48.208 -33.348 -23.212 -21.026 -14.329 -9.292 -5.388 -2.305 0.194 2.226 5.339 7.572 10.991 12.815

Bd -1638.029 -1157.746 -1121.907 -1069.044 -751.967 -427.835 -276.562 -230.330 -194.475 -144.169 -110.758 -70.350 -47.363 -32.566 -22.601 -20.458 -13.949

-2.101 0.262 2.284 5.358 7.484 10.875

a

Expt.e Expt. uncert.f -1916.87g -1307.84g -1262.93g -1199.69g -816.66 43.0 -451.62 9.0 -289.91 6.0 -240.65 5.0 -203.17 4.0 -150.49 1.5 -115.72 1.5 -73.61 1.0 -49.69 1.0 -34.62 1.0 -24.41 1.0 -22.22 -15.49 -10.43 -6.50 -3.37 -0.84 1.26 4.49 6.85 10.57 12.64

Rigid-monomer values computed using the rectangle method of integration. hV i ∆Bflex 0 (T ) computed with the CCpol-8sf[2014] potential using the method of integration of Conroy. 131 c Sum of columns 3 and 4. d Theoretical coefficients from refs 119 and 134 computed using the PIMC method and the HBB2-pol potential. e Fit to experimental data from ref 35. f Estimate of experimental uncertainty from ref 35. g Values beyond the range of validity of the experimental fit of ref 35. b



110

6000

100

(a)

90

5000

SAPT-5s’f[2014] [kcal/mol]

80 70

4000

60 50

3000

40 30

2000

20 10

1000

0 -10

0 -10

0

10 20 30 40 50 60 70 80 90 100 110 ab initio SAPT [kcal/mol]

20

1000 900

(b) SAPT-5s’f[2014] [kcal/mol]

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 82 of 89

800 700

10

600 500 400 0

300 200 100

-10

0 -10

0 10 ab initio SAPT [kcal/mol]

20

Figure 1: (a) Interaction energies predicted by SAPT-5s′ f[2014] versus the ab initio SAPT values. (b) Similar plot for the most important energy range. The contour lines enclose the areas where there were at least one (blue line) or five (red line) counts. The color shows the number of correlation points in each square of side 1.0 kcal/mol for (a) and 0.25 kcal/mol for (b). 82 ACS Paragon Plus Environment

Page 83 of 89

20

800

SAPT-5s’fIR[2014] [kcal/mol]

700

(a) 600 10 500 400 300 0 200 100

-10

0 -10


20

20

1000 900


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


800 700

10

600 500 400 0

300 200 100

-10

0 -10

0 10 SAPT-5s’fIR[2014] [kcal/mol]

20

Figure 2: (a) Interaction energies predicted by the SAPT-5s′ fIR[2014] surface versus the ab initio SAPT values. (b) The SAPT-5s′ f[2014] vs. SAPT-5s′ fIR[2014] energies on the same grid of the ab initio points. The plots are of the same type as in Fig. 1. 83 ACS Paragon Plus Environment


0 -200

(a) Experiment

-400

B(T) [cm3 mol-1]

Bcl[CCpol-8s] -600

BTI QM[CCpol-8s] -800

Bflex[CCpol-8sf] -1000

Bflex[CCpol-8sfIR] -1200

BTI QF -1400

B[HBB2-pol] -1600 -1800 -2000 400

600

800

1000

1200

1400

Temperature [K] 300

Experiment

B(T) - Bexp(T) [cm3 mol-1]

250

(b)

Bcl[CCpol-8s] TI [CCpol-8s] BQM

200

Bflex[CCpol-8sf] Bflex[CCpol-8sfIR]

150

TI BQF

100

B[HBB2-pol] B[MB-pol]

50

0

-50 250

300

350

400

450

500

550

600

650

700

750

800

2800

3000

Temperature [K] 10

Experiment

(c)

Bcl[CCpol-8s] 5


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 84 of 89

TI [CCpol-8s] BQM

0

-5

Bflex[CCpol-8sf] Bflex[CCpol-8sfIR] B[HBB2-pol]

-10

800

1000

1200

1400

1600

1800

2000

2200

2400

2600

Temperature [K]

Figure 3: Comparison of theoretical virial coefficients B(T ) with experimental values. (a) Virial coefficients; (b) and (c) Difference of theoretical predictions from the experimental values of Harvey and Lemmon 35 for low and high temperatures, respectively. The virial coefficients for CCpol-8s surface 62 calculated in this work. The coefficients with the flexibility effects included (all 2014 versions) calculated in the present work are defined as Bflex [PES] TI = BQM [CCpol-8s] + ∆Bflex [PES]. The values of B[HBB2-pol] are taken from refs 119 and TI 134, the values of B[MB-pol] from ref 85,134, and the values of BQF from ref 117. The lines only guide the eye. 84 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


Graphical TOC Entry 60 50

experiment

40

rigid classical rigid quantum

30

B(T)-Bexp(T) [cm3 mol-1]

Page 85 of 89

flexible classical 20

theory total

10 0 -10 -20 -30

Comparison of the effects of monomer flexibility to quantum

-40

corrections in the second virial coefficient, B(T), of water

-50 -60 350

400

450

500

550

600

650

700

750

800

850

900

950

1000

Temperature [K]



110

6000

100

(a)

90

5000

SAPT-5s’f[2014] [kcal/mol]

80 70

4000

60 50

3000

40 30

2000

20 10

1000

0 -10

0 -10

0

10 20 30 40 50 60 70 80 90 100 110 ab initio SAPT [kcal/mol]

20

1000 900


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 86 of 89

800 700

10

600 500 400 0

300 200 100

-10

0 -10



20

Page 87 of 89

20

800

SAPT-5s’fIR[2014] [kcal/mol]

700

(a) 600 10 500 400 300 0 200 100

-10

0 -10


20

20

1000 900


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


800 700

10

600 500 400 0

300 200 100

-10

0 -10

0 10 SAPT-5s’fIR[2014] [kcal/mol]


20


0 -200

(a) Experiment

-400

B(T) [cm3 mol-1]

Bcl[CCpol-8s] -600

BTI QM[CCpol-8s] -800

Bflex[CCpol-8sf] -1000

Bflex[CCpol-8sfIR] -1200

BTI QF -1400

B[HBB2-pol] -1600 -1800 -2000 400

600

800

1000

1200

1400

Temperature [K] 300

Experiment


250

(b)

Bcl[CCpol-8s] TI [CCpol-8s] BQM

200

Bflex[CCpol-8sf] Bflex[CCpol-8sfIR]

150

TI BQF

100

B[HBB2-pol] B[MB-pol]

50

0

-50 250

300

350

400

450

500

550

600

650

700

750

800

2800

3000

Temperature [K] 10

Experiment

(c)

Bcl[CCpol-8s] 5


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 88 of 89

TI [CCpol-8s] BQM

0

-5

Bflex[CCpol-8sf] Bflex[CCpol-8sfIR] B[HBB2-pol]

-10

800

1000

Plus2000 Environment 1200 ACS 1400 Paragon 1600 1800 2200 2400 Temperature [K]

2600

Page 89 of 89

60 50

experiment

40

rigid classical rigid quantum

30

B(T)-Bexp(T) [cm3 mol-1]

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


flexible classical 20

theory total

10 0 -10 -20 -30

Comparison of the effects of monomer flexibility to quantum

-40

corrections in the second virial coefficient, B(T), of water

-50 -60 350

400

450

500

550

600

650

700

750

Temperature [K] ACS Paragon Plus Environment

800

850

900

950

1000

Ab initio potential energy curves of the valence, Rydberg, and ion-pair states of iodine monochloride, ICl.

Ab initio RNA folding.

Flexible nuclear screening approximation to the two-electron spin-orbit coupling based on ab initio parameterization.

Analytical potentials from "ab initio" computations for the interaction between biomolecules. 1. Water with amino acids.

Flexible Force Field Parameterization through Fitting on the Ab Initio-Derived Elastic Tensor.

The benzoic acid-water complex: a potential atmospheric nucleation precursor studied using microwave spectroscopy and ab initio calculations.

Argon Interaction with Gold Surfaces: Ab Initio-Assisted Determination of Pair Ar-Au Potentials for Molecular Dynamics Simulations.

Structure and Raman spectra in cryolitic melts: simulations with an ab initio interaction potential.

Erratum: Development of a "First-Principles" Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient.

AIDA: ab initio domain assembly server.

Ab-initio SCF investigation ofγ-aminobutyric acid.

Ab initio non-relativistic spin dynamics.

Ab Initio Screening of CO2-philic Groups.

Ab initio H2O in realistic hydrophilic confinement.

Ab initio molecular dynamics simulation of liquid water by quantum Monte Carlo.

Relationship between structural order and water-like anomalies in metastable liquid silicon: Ab initio molecular dynamics.

Ab initio engineering of materials with stacked hexagonal tin frameworks.

Quantum fluctuations and isotope effects in ab initio descriptions of water.

Water-accelerated OH addition to sulfur dioxide SO₂: direct ab initio molecular dynamics (AIMD) study.

Photo-driven oxidation of water on α-Fe2O3 surfaces: an ab initio study.

Microsolvation of LiBO2 in water: anion photoelectron spectroscopy and ab initio calculations.

mineral interfaces from ab initio molecular dynamics simulations.

Characterization of adsorbed water in MIL-53(Al) by FTIR spectroscopy and ab-initio calculations.

Molecular Mechanical Molecular Dynamics Simulation.