Origin of Many-Body Vibrational Frequency Shifts in Water Clusters

Subscriber access provided by UNIVERSITY OF CALGARY

A: Spectroscopy, Molecular Structure, and Quantum Chemistry

Origin of Many-Body Vibrational Frequency Shifts in Water Clusters Joseph P Heindel, Elizabeth S Knodel, and Daniel P. Schofield J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b04380 • Publication Date (Web): 20 Jul 2018 Downloaded from http://pubs.acs.org on July 22, 2018

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 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 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.

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 37 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

Origin of Many-Body Vibrational Frequency Shifts in Water Clusters Joseph P. Heindel,†,‡ Elizabeth S. Knodel,†,‡ and Daniel P. Schofield∗,† Department of Chemistry and Biochemistry, Seattle Pacific University, Seattle, WA 98119 E-mail: [email protected]

∗

To whom correspondence should be addressed Department of Chemistry and Biochemistry, Seattle Pacific University, Seattle WA 98119 ‡ Contributed equally to this work †

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 We have demonstrated the application of many-body expansions to calculations of the anharmonic, local-mode, OH-stretching vibrational frequencies of water clusters. We focused on five low-lying isomers of the water hexamer and the DD*(20,1) isomer of (H2 O)21 . Our approach provides accurate OH-stretching vibrational frequencies when treating 1- and 2-body interactions with the CCSD(T)-F12 level of theory and the 3and 4-body interactions with the DF-MP2-F12 level. Additionally, we have investigated the physical origin of the large contribution that 2- and 3-body interactions make to the shifts of vibrational frequencies using symmetry-adapted perturbation theory (SAPT). We conclude that while 2-body vibrational frequency shifts can be correlated linearly with electrostatic energies, all strongly shifted 3-body interactions can be correlated to the induction energy with a single regression coefficient of approximately 70 cm−1 (kcal · mol−1 )−1 .

Introduction Recently, many-body expansions (MBEs) have been applied in computational chemistry both as a means of approximation and as an analysis tool in their own right. 1–15 The earliest studies were focused on the many-body interactions of various water clusters as a way of understanding the convergence of the many-body energy to the total interaction energy. 16–18 It soon became clear that many other frequently-studied properties could be couched in terms of a MBE. For example, schemes for finding the molecular dipole moment 2 and polarizability 2 as well as the contribution of basis-set superposition error (BSSE) to the binding energy in MBEs have been presented. 19 Following their demonstration of the accuracy of a 2-body:many-body approach to geometry optimizations, 13,14 Tschumper et al developed a framework for calculating the normal-mode vibrational frequencies using a MBE. 15 Their model requires a correct decomposition of the analytic Hessian into its many-body components and has proven to be very accurate for both geometry optimizations and frequency 2


Page 2 of 37

Page 3 of 37 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. 13–15 Bowman et al. have used a MBE to calculate the potential energy surface (PES) of arbitrary sized water clusters. This PES was constructed from 40,000 energy calculations on water dimers at the CCSD(T)/aug-cc-pVTZ level and 30,000 MP2/aug-cc-pVTZ energy calculations on water trimers. 20 Through the use of this PES and their previously-developed local-monomer approach to frequency calculations, 20 Bowman et al. have simulated spectra of 6 different water hexamers and some larger water clusters. 21,22 Their calculated spectra are well-aligned with experimental spectra in the fundamental OH-stretching region, especially when both the OH-stretching mode and the HOH-bending overtones are included in the vibrational Hamiltonian. 22,23 MBEs have also been used to parameterize force fields for molecular dynamics (MD) simulations. Paesani et al. parameterized a potential energy function (MB-Pol) for classical water simulations on the basis of many-body interactions in water. 9,10 This potential is built from some 40,000 CCSD(T)/CBS dimer calculations and over 12,000 CCSD(T)/augcc-pVTZ trimer calculations. In this model, long-range electrostatic and dispersion interactions are accounted for by surrounding point charges. The MB-Pol potential has been used in classical molecular dynamic simulations to accurately reproduce the second and third virial coefficients of water as a function of temperature. 10 A thorough review of the MB-Pol potential has been published recently. 1

Figure 1: From left to right, the cage and prism isomers of (H2 O)6 and the DD*(20,1) isomer of (H2 O)21 In this work, we use a MBE to study the vibrational properties of strongly interacting systems where both intermolecular interactions and vibrational anharmonicities will have 3



a strong effect on the final frequency, specifically the OH-stretching spectra for clusters of water molecules (such as those shown in Figure 1). We calculate the frequency shifts for a number of water clusters, where the frequency shift is defined as the difference between the OH-stretching frequency of the water cluster and that of an isolated gas-phase water molecule. We wish to avoid the diagonalization of Hessians required for normal mode analysis, as this becomes prohibitively expensive for large systems. Additionally, we wish to devise a scheme that could be extended to quickly and accurately compute the vibrational frequencies of heterogeneous systems without reparameterizing an entire PES. Lastly, we seek to use the many-body approach to gain insight into the interactions that explain the physical origin of the vibrational frequency shifts. In past work, the relationship between the geometric arrangement of a given water cluster and its resulting frequency shift has been investigated. 24,25 These studies illustrate the impact of the total cluster environment on the OH-stretching frequency, but do not explicitly consider the correlation of the 2- and 3-body components to the overall vibrational frequency shift. Recently, Ouyang and Bettens have demonstrated the effect of induction on the interaction energy in non-branching chains of hydrogen-bound water molecules. 26 We use a series of statistical models to demonstrate the relationship between the 2- and 3-body energetic interactions, calculated with symmetry adapted perturbation theory (SAPT), and the corresponding 2- and 3-body frequency shifts calculated from a MBE.

Theoretical Details The OH-stretching vibrational frequencies of five water hexamers and the DD*(20,1) isomer of (H2 O)21 were calculated with a MBE and a local-mode model Hamiltonian. The local mode model treats each stretching vibration as being decoupled from the other vibrational modes. This model has been successfully applied to many systems containing OH-bonds. 27–30

4


Page 4 of 37

Page 5 of 37 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


Many-body properties were calculated according to the scheme below:

∆x(1) a (2) ∆xa,j

= x(a) − xw =

N X j=1

(3) ∆xa,j,k

=

(1a)

[x(a, j) − x(a) − x(j)] (j 6= a)

N −1 X N X j=1 k>j

(1b)

[x(a, j, k) − x(a, j) − x(a, k) − x(j, k) + x(a) + x(j) + x(k)] (j 6= a; k 6= a) (1c)

(4) ∆xa,j,k,l

=

N −2 N −1 X N X X j=1 k>j l>k

[x(a, j, k, l) − x(a, j, k) − x(a, j, l) − x(a, k, l) − x(j, k, l) + x(a, j) + x(a, k) + x(a, l) + x(j, k) + x(j, l) + x(k, l) − x(a) − x(j) − x(k) − x(l)] (j 6= a; k 6= a; l 6= a)

(1d)

Where x is the property of interest (e.g. energy, frequency, or dipole moment) and each term, ∆x(n) , represents the n-body contribution to the total shift or change from the reference, xw . The vibrational frequency is determined by the shape, or curvature, of the potential energy as a function of the OH-bond length. Similarly, the vibrational intensity depends on the shape of the dipole moment function as a function of OH-bond length. As the OH-bond length of monomer a changes, all terms in Equations 1a-1d that do not involve index a contribute an identical energy and dipole at each value of the OH-bond length and therefore do not affect the vibrational potential energy or dipole moment curves. These terms can accordingly be ignored. This reduces the number of terms that need to be calculated. Notice also that the sums in Equations 1a-1d only run over the indices not including the monomer containing the bond of interest, a. We have used two alternate approaches to generate anharmonic local mode OH-stretching frequencies via a MBE. These are described as method I and II below:

5



Method I In the instance where ∆x is a change in the total vibrational frequency, the frequencies are calculated from a one-dimensional OH-stretching potential energy curve (PEC) for each n-mer (where n stands for di-, tri-, and tetra-) at the geometry they have within the overall water cluster (Eq. (1)). To find the energy at each point along our 1-D PEC, we fragment a cluster into all monomers, dimers, trimers, and tetramers that contain the O−H bond of interest. For example, for a given OH-stretch within a water hexamer, the cluster would be fragmented into one monomer, five dimers, ten trimers, and ten tetramers. In this work, we truncate the many-body expansion at fourth order and, as a result, the interactions from the five pentamers and one hexamer are not included. We obtain the 1-D PEC by displacing the hydrogen in the OH-bond of interest by 0.05 ˚ A around the equilibrium bond length, re , ranging from re − 0.4 ˚ A to re + 0.6 ˚ A. The OH-stretching frequencies of each n-mer are found by numerically solving the Schrodinger equation for each 1-D slice of the PEC using a finite element method as in previous work. 30 A PEC of this size provides converged energy and wavefunction solutions. The calculated frequencies are then combined as described in Eq. (1) to find the local mode OH-stretching frequency. This procedure is repeated for each of the O−H bonds and, in this way, a 1-D slice of the PES along each OH-stretching bond coordinate is constructed and each of the 12 OH-stretching frequencies in the water hexamer are determined.

Method II Alternatively, we have calculated the local-mode frequencies by constructing a single potential energy curve for each OH-bond using the many-body sum for energies (Eq. (1)). Rather than constructing individual potential energy curves for each n-mer, this alternative method uses a MBE to find the energy of the whole cluster as a function of the OH-bond length for a single OH-stretching coordinate. In this approach, the Schrodinger equation is numerically solved only once per local-mode as opposed to being solved for each of the monomers, dimers, 6


Page 6 of 37

Page 7 of 37 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


trimers, and tetramers, as described above for Method I. Another variation on each of the methods described above was to replace the numerical solution of the vibrational Schrodinger equation with an analytical solution through use of the Morse potential to describe the OH-stretching coordinate. The Morse potential has vibrational energy levels given by: E(v)/(hc) =

1 v+ 2

2 1 ω ˜− v+ ω ˜ x. 2

(2)

where ω ˜ is the OH-stretching frequency and ω ˜ x is the anharmonicity. In keeping with previous work, we find ω ˜ from the second derivative of the potential energy curve, while ω ˜x is found from the third and fourth derivatives of the potential energy. 31 These derivatives are found by fitting an 8th order polynomial to a calculated local mode potential energy curve consisting of nine points separated by 0.05 ˚ A from −0.2 to +0.2 ˚ A around the equilibrium bond length. The oscillator strengths of the OH-stretching vibrations were calculated using the local mode model and a many-body calculated dipole moment curve. The many-body dipole moment is converged at 2-body 2 and this dipole moment is used to compute the oscillator strength of the fundamental transition according to: 32

f10 = 4.702 × 10−7 [cm D−2 ]˜ ν10 h1|~µ|0i2

(3)

where ν˜10 is the fundamental OH-stretching frequency, µ ~ is the dipole moment function, and h1| and |0i are the first vibrational excited state and ground vibrational state wavefunctions respectively. In a manner similar to that used for the potential energy, we fit the dipole moment curve with a sixth-order, Taylor series polynomial which allows for evaluation of the transition dipole moment matrix element according to:

h1|~µ|0i =

∂~µ 1 ∂ 2µ ~ 1 ∂ 3µ ~ 2 h1|q|0i + h1|q |0i + h1|q 3 |0i · · · 2 3 ∂q 2 ∂q 6 ∂q 7


(4)


Page 8 of 37

where q is the vibrational displacement coordinate. The integrals h1|q n |0i are evaluated using either trapezoidal numerical integration in the instance where the wavefunctions are obtained from a numerical solution to the Schrodinger Equation, or analytically if Morse oscillator wavefunctions are used. 30 As will be detailed later, basis-set superposition error (BSSE) can have a large effect on the results obtained from fragmentation schemes. 19 All of our single point energies are counterpoise-corrected, unless otherwise noted, according to the following equations:

∗ EAB

AB = EAB + (EAA + EBB ) − (EAAB + EBAB )

∗ ABC EABC = EABC + (EAA + EBB + ECC ) − (EAABC + EBABC + ECABC )

(5a) (5b)

∗ ABCD D EABCD = EABCD + (EAA + EBB + ECC + ED ) ABCD − (EAABCD + EBABCD + ECABCD + ED )

(5c)

∗ Eq. (5a) gives the BSSE-corrected energy of dimer AB, denoted by EAB , by taking the

energy of dimer AB in the AB basis and finding the BSSE associated with each monomer by subtracting the energy of each in its own basis (i.e. EAA ) and subtracting the energy of that same monomer in the dimer basis (i.e. EAAB ). Eqs. (5b) and (5c) have an analogous interpretation. Note that this only counterpoise-corrects for BSSE due to the monomer basis as in the original Boys Bernardi scheme. 33 It has been shown that the BSSE is much larger for the 2-body term than for the 3- or 4-body terms. 34 As discussed later, it is more practical to minimize BSSE through the use of explicitly correlated ab initio methods and short OH-bond displacements than it is to use a many-body counter-poise correction. All water hexamer geometries were optimized with the DF-MP2 level of theory and the def2-TZVPP basis set 35 augmented with diffuse functions produced by an even-tempered extension of the last two exponents of each angular momentum type (ATZVPP). The geometry of (H2 O)21 was optimized at the DF-MP2/aug-cc-pVDZ level of theory. All subsequent

8


Page 9 of 37 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 for the many-body expansion were performed at both the CCSD(T)-F12/VDZF12 and DF-MP2-F12/VDZ-F12 levels of theory using MOLPRO. 36–40 The CCSD(T)-F12 calculations were performed using the CCSD(T)-F12a ansatz, and the triples were scaled as implemented in MOLPRO. 38 The dipole moments were also calculated at the CCSD(T)F12/VDZ-F12 level by finite differences using applied electric field strengths of ±0.005 atomic units. These methods were chosen since explicitly correlated coupled-cluster calculations have been shown to converge more quickly than other methods for a given basis set, and hence suffer less from BSSE. 41 The MBE we use gives explicit information about the contribution of each n-mer to the total frequency shift. To facilitate analysis of which interactions in each sub-cluster are most important, we have calculated interaction energies using 2- and 3-body symmetry-adapted perturbation theory (SAPT). SAPT has been reviewed extensively elsewhere, 42 so we only cover the pertinent details here. SAPT has the advantage of calculating the individual components of the interaction energy separately. Energetic contributions to the 2-body energy can be classified in terms of electrostatic, exchange, induction, and dispersion interactions. Likewise, 3-body contributions to the overall interaction energy are calculated in terms of induction, dispersion, and exchange. An explanation of the grouping of the individual terms is provided in the Supporting Information. All 2-body SAPT calculations were performed using the SAPT(2+3) approach as implemented in the Psi4 quantum chemistry package. 43–45 The ATZVPP basis set was used in all 2-body SAPT calculations. The 3-body SAPT calculations are performed using the DFT-SAPT approach with a PBE0 description of the monomers and an aug-cc-pVDZ basis set. To ensure accurate calculation of the dispersion energy, the PBE0 functional was asymptotically corrected with Fermi-Amaldi scheme and the experimental ionization potential of 0.4638 atomic units for the water monomer. 46 The SAPT code of Szalewicz et al 46 was used in conjunction with the DALTON program 47,48 for all 3-body SAPT calculations. We have statistically correlated the individual SAPT-calculated components of the in-

9



Page 10 of 37

teraction energy for a given dimer or trimer with the 2- and 3-body frequency shift of a particular OH-stretching local-mode. The data were correlated with a Bayesian, penalizedregression model using the Flexmix program in the R statistical software package. 49 This model allows for the possibility of distinct groupings of OH-stretches to be correlated with different statistical relationships to the SAPT-calculated interaction energies. For example, in a given dimer, up to four different OH-stretching frequencies will arise for a single set of SAPT-calculated energies. In this approach, the Flexmix algorithm sorts our calculated frequency shift-SAPT energy pairs into one of a number of groups. Each group is fitted with a linear regression model and the overall fit is simultaneously optimized in terms of both group population and regression coefficients. To account for the correlation between individual SAPT calculated energies and to prevent over-fitting, we used a penalized regression models that adds a penalty to the statistical model as the coefficient size increases. 50 This penalty leads to a preference for statistical models with fewer independent variables, but also biases the model towards fitting with only the variable/s of largest magnitude. In order to address this bias, the SAPT-calculated interaction energies were scaled prior to being fitted. The interaction energies were scaled by subtracting the respective mean from each component of the interaction energy and then dividing by the standard deviation. This allows each of the energetic components to be of approximately the same magnitude, and ensures that the penalized regression model is not biased toward fitting the variables of largest magnitude. In addition to the automatic group sorting performed by Flexmix, we used geometric criteria to pre-assign individual dimers and trimers into categories before fitting each category with a penalized regression model. Each OH-bond was classified as either hydrogen-bound or not hydrogen-bound using the definition proposed by Kumar et al: 24,51

N (r, ψ) = exp(−r/0.343˚ A)(7.1 − 0.050ψ + 0.00021ψ 2 )

10


(6)

Page 11 of 37

Here r is the distance between a hydrogen and the oxygen of another water molecule and ψ is the angle between this same OH vector and the normal vector of the plane bisecting the acceptor water. Through use of these N -values, we pre-sorted each OH-bond into one of four geometrically-defined groups (vide infra), and the associated OH-frequency was subsequently correlated to the SAPT energetic parameters using a penalized regression model.

Results and Discussion The calculated OH-stretching spectra for the five lowest-lying isomers of the water hexamer are shown in Figure 2. The OH-stretching frequencies were calculated using Method II and a Morse description of the vibrational potential. The spectra are in fairly good agreement with previously published spectra by Bowman et al. 20 Much of the discrepancy is likely due to the combined effects of ignoring coupling between the OH-stretch and other vibrational modes, and the use of different water hexamer optimized geometries (vide infra). bag book cage prism ring

Relative Intensity

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


3000

3200

3400

Frequency / cm

3600

3800

−1

Figure 2: Simulated OH-stretching vibrational spectra for each of the 5 water hexamers considered. The peaks were convoluted with a Lorentzian of full-width at half-maximum of 5 cm−1 . In Table 1 the n-body frequency shifts calculated with Method I for all twelve O−H stretches of the cage isomer of (H2 O)6 are listed. The predicted many-body frequencies that

11



result from adding the sum of these shifts to the calculated OH-stretching frequency of the water monomer are also listed in Table 1. The cage isomer is used here because it is one of the lowest energy hexamers that has a fully 3-dimensional hydrogen-bonding framework. It is therefore a better representation of larger clusters and liquid water than the ring or book hexamers where the hydrogen-bonding is largely 2-dimensional. The calculated frequencies for the bag, prism, ring, and book isomers of (H2 O)6 and the central water in the DD*(20,1) isomer of the (H2 O)21 cluster are given in the Supporting Information. One detail which is immediately striking about the vibrational frequency shifts in Table 1 is that, for strongly red-shifted vibrations, the contribution to the total shift from the 3-body interactions are often as large or larger than the 2-body shifts. Additionally, for the most strongly shifted OH-frequency, the 4-body frequency shift is close to 100 cm−1 . These significant 3- and 4-body frequency shifts are in contrast to the trends observed for dipole moment and binding energy. The 2-body term of the dipole moment and polarizability expansions have been observed to capture 95% of the total value for hydrogen-bonded clusters, and the 2-body term of the binding energy usually accounts for 85% of the total binding energy, with the 3-body term capturing nearly all of the remaining 15%. 2,3,14 Through the use of the MBE approach to the calculation of vibrational frequencies, we are able to analyze the effect of each n-mer on the calculated frequency shift. For the most strongly red-shifted peaks, all five dimers contribute to the shift in vibrational frequency, with a particularly large shift when the OH-bond of interest is acting as a hydrogen bond donor. If we extend this analysis to the other red-shifted OH-stretching transitions in both this cluster and the other hexamers we considered, this trend is also observed. In contrast to what we observed for the 2-body shift, only a subset of trimers contribute to the 3-body shift. In the hexamer structures considered, nearly all of the 3-body shift comes from four of the ten trimers, with the remaining six trimers exhibiting a near zero frequency shift. The four trimers that contribute to the calculated 3-body frequency shift all contain the hydrogen-bound dimer that produced the largest 2-body shift while the six

12


Page 12 of 37

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


13

2-Body -263.9 -215.4 -212.7 -182.3 -139.1 -105.2 -118.9 -94.8 -4.1 0.1 -3.4 2.7

CCSD(T)-F12 3-Body 4-Body -387.4 -96.6 -225.4 -37.4 -185.3 -22.4 -171.4 -37.5 -136.3 -25.1 -57.9 -4.4 -47.6 -4.3 -49.9 3.0 5.2 -0.4 4.5 -1.6 5.5 -1.2 5.8 -0.7 Many-bodya 2987 3257 3315 3339 3418 3546 3547 3571 3728 3731 3732 3739 1-Body 20.9 20.8 21.3 16.2 3.8 -1.4 2.4 -1.7 11.7 12.2 15.2 16.1

2-Body -343.2 -281.6 -272.5 -240.1 -189.2 -138.0 -149.8 -121.6 -4.9 -4.7 -8.3 1.9

DF-MP2-F12 3-Body 4-Body -389.9 -98.5 -225.5 -36.5 -187.5 -21.4 -173.4 -36.0 -138.4 -22.5 -60.5 -1.9 -51.1 -1.4 -52.6 5.1 2.2 1.5 0.4 1.0 2.3 1.3 3.1 1.1

Many-bodyb 2943 3231 3294 3320 3407 3552 3554 3583 3764 3763 3764 3776

The predicted many-body frequency is calculated from a sum of n-body frequency shifts, n=1-4, and a CCSD(T)-F12/VDZ-F12 calculated local mode frequency of 3724 cm−1 for the OH-stretching frequency of the water monomer. b The predicted many-body frequency is calculated from a sum of n-body frequency shifts, n=1-4, and a DF-MP2-F12/VDZ-F12 calculated local mode frequency of 3754 cm−1 for the OH-stretching frequency of the water monomer.

a

ν1 ν2 ν3 ν4 ν5 ν6 ν7 ν8 ν9 ν10 ν11 ν12

1-Body 11.5 11.4 11.8 -5.6 6.6 -11.0 -6.6 -11.2 3.1 4.5 6.9 7.7

Table 1: Calculated n-body frequency shifts (cm-1 ) and predicted many-body frequencies for the OH-stretching modes of the cage isomer of (H2 O)6

Page 13 of 37 The Journal of Physical Chemistry


that have a negligible effect on the 3-body frequency shift do not. Likewise, in the case of the 4-body interactions, nearly all of the 4-body frequency shift arises from those tetramers that contain the hydrogen-bound dimer that produces the largest frequency shift. The observations above suggest an approximate, truncated MBE for calculating the OHstretching frequency shifts is possible at a much lower computational cost than the full MBE. As noted, the majority of the vibrational frequency shift arises from those clusters which contain the OH-bond of interest acting as a donor in a hydrogen bound dimer. This suggests that the overall OH-stretching vibrational frequency could be approximated by considering only the trimers and tetramers containing that dimer. In order to test the efficacy of this procedure, the truncation scheme was used on the DD*(20,1) isomer of (H2 O)21 . The two O−H stretches studied are those of the water molecule located inside the dodecahedral cage of DD*(20,1). While this analysis was also done on the water hexamers (the results of which can be found in the Supporting Information), the (H2 O)21 structure, by virtue of its size, has a greater number of trimers and tetramers that will not be included in this truncation scheme. As a result it is likely to have a larger discrepancy between the truncated scheme and the full calculation than for the water hexamers. When considering the (H2 O)21 structure, only 19 of the 190 possible trimers contained the hydrogen bound dimer of interest, and 171 of the 1140 possible tetramers contain the same hydrogen-bound dimer. This leads to approximately a 7-fold decrease in computational expense. The results of three different truncation schemes are summarized in Table 2. First, only the 4-body terms are truncated. Secondly, 4-body terms are neglected entirely. Finally, 4-body terms are neglected entirely and 3-body terms are truncated. The non-truncated MBE frequencies out to the 4-body terms are included for the sake of comparison. It can be seen in Table 2 that the error accumulated from ignoring a very large number of terms in the expansion is quite small. The total 3-body frequency shift is on the order of 200 cm−1 , and the truncated decomposition reproduces about 94% of the total shift. The 4body case reproduces nearly the entire shift of ∼ 40 cm−1 , missing by only one or two cm−1 . 14


Page 14 of 37

Page 15 of 37 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: OH-stretching frequencies of the DD*(20,1) isomer of (H2 O)21 calculated using the full and truncated many-body expansions Expansion type OH-1 OH-2 Difference for OH-1 Difference for OH-2 a b

Full 4-body 3231 3271

Truncated 4-bodya 3233 3271 -2 0

Full 3-body 3274 3310

Truncated 3-bodyb 3265 3293 9 17

Difference calculated in comparison to the full 4-body reference. No tetramers are included. Difference calculated in comparison to the full 3-body reference.

This conforms to the intuition that higher-order interactions ought to be smaller in general, so whatever interactions we ignore are smaller in magnitude for the 4-body case. Since an error on the order of 10−20 cm−1 is likely less than the error due to our choice of vibrational model and the ab initio method itself, neither of these truncation schemes introduce significant uncertainty. The truncated frequencies tend to overestimate the vibrational frequency shift, a trend which was observed across all of the other water clusters as well. This tells us that the 3- and 4-body interactions which do not have our hydrogen-bound dimer are weakly destabilizing overall and shift the total calculated frequency towards the blue. To further decrease the computational cost it is possible to calculate the truncated 4body terms at a lower level of theory. For example, all 1- and 2-body terms could be calculated with the CCSD(T)-F12 method, all 3-body terms with the DF-MP2-F12 level, and a truncated set of 4-body terms also at the DF-MP2-F12 level of theory. As shown in Table 1 if one chooses to use the DF-MP2-F12 level of theory to calculate the frequencies, the 3- and 4-body frequency shifts differ by at most a few cm−1 from the CCSD(T)-F12 calculated values. This can save a significant amount of computer time with only a minimal loss of accuracy. In Table 3 we present a comparison of the OH-stretching frequencies of the prism isomer of (H2 O)6 calculated using Methods I and II, with either a numerical or Morse oscillator solution to the Schrodinger equation. The frequencies determined from a supramolecular local-mode analysis of the entire cluster are also presented. With both approaches the 15



16

b

a

Method 2966 3249 3354 3392 3537 3542 3580 3606 3636 3730 3731 3734

Method I 2978 3256 3358 3395 3538 3540 3579 3602 3636 3727 3729 3731

Morse Potential Method II Full Cluster 2956 2979 3253 3271 3356 3372 3394 3410 3541 3555 3542 3558 3579 3591 3602 3617 3636 3643 3727 3735 3729 3736 3731 3739

Full Cluster 3034 3282 3385 3424 3533 3539 3583 3600 3642 3724 3726 3728

b

3092 3256 3372 3442 3482 3521 3579 3588 3630 3697 3706 3728

LMM 20

Frequencies for Method I and Method II calculated at the CCSD(T)-F12/VDZ-F12 level from a n-body expansion of all terms up to 4-body Calculated at the CCSD(T)/aug-cc-pVDZ optimized geometry

Bond OH-1 OH-2 OH-3 OH-4 OH-5 OH-6 OH-7 OH-8 OH-9 OH-10 OH-11 OH-12

Numerical Solution I Method II Full Cluster 2936 2956 3247 3261 3354 3367 3393 3406 3541 3554 3545 3557 3581 3591 3607 3618 3637 3645 3730 3736 3731 3737 3734 3740

a

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

Table 3: OH-stretching frequencies of the prism isomer of the water hexamer.

The Journal of Physical Chemistry Page 16 of 37

Page 17 of 37 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


many-body frequencies calculated with terms up to fourth order are in good agreement with those obtained from the supramolecular approach. To further check the convergence of the many-body expansion we have calculated the 5- and 6-body contributions to the many-body OH-stretching frequencies of the bag isomer, the isomer for which the 4-body contributions are largest. These frequencies are presented in the supplementary information, and confirm that the 5-body terms for all but the most strongly redshifted modes are on the order of 1 cm−1 . In the case of the strongly redshifted modes, the 5-body shift is around 15-20 cm−1 to the red. For all OH-stretching modes, the 6-body shift is positive, likely a consequence of accounting for BSSE that has accumulated in the 3-, 4- and 5-body terms. As can be seen in Table 3 frequencies calculated using Method II are systematically more red-shifted for the strongest hydrogen bonds. Additionally, it is clear that the frequencies calculated using Method I are closer to the supramolecular values than those calculated from Method II. We attribute this to BSSE. Method II and the supramolecular approach each suffer from greater BSSE than Method I because the frequencies are calculated from only a single potential energy curve. This means that all frequencies inherit the 2-, 3-, and 4-body BSSE. On the other hand, using Method I, the largest frequency shifts, namely the 2-body and 3-body shifts, do not suffer from the BSSE present in 4-body interactions. The supramolecular approach has the advantage of incorporating 5- and 6-body interactions, but this is only a small effect. This result clarifies our motivation for using explicitly correlated methods, as they converge more quickly with basis set size and suffer less from BSSE even when using smaller basis sets. Additionally, it can be seen that the frequencies which come from treating each potential as a Morse potential shift the results from Method II much closer to the results from Method I than when the Schrodinger Equation is solved numerically. This is a consequence of the shorter OH-bond length range that is sampled when determining the frequencies from the Morse potential than the numerical potential. At long OH-bond lengths the BSSE is largest, and this is further evidence that Method II suffers more from BSSE than Method I.

17



These results make a good case for simply using the Morse oscillator as a description of the local modes of water clusters because it avoids BSSE as much or more than one can feasibly correct using typical counter-poise correction schemes. This is a significant timesaving approach since the Morse potential produces converged frequencies from a potential energy curve constructed from 9 single point energies compared to the 21 single point energies required in numerical solutions to the Schrodinger equation. To test the influence of geometry on the OH-stretching frequency, we have used the Morse oscillator method to calculate the OH-stretching frequency for the prism hexamer at the CCSD(T)/aug-cc-pVDZ optimized geometry calculated by Miliordos et al. 52 The most pronounced difference in frequency is for the most red-shifted mode, which changes from 2978 cm−1 (Table 3) when using the DF-MP2 geometry to 3034 cm−1 with the CCSD(T) geometry. The change for the other OH-stretches is on the order of 25 cm−1 or less. Our calculated frequencies are in reasonable agreement with those previously calculated using a local-monomer model (LMM). 20 The largest discrepancy is for the most strongly redshifted mode. The source of this discrepancy is likely two-fold: the LMM accounts for coupling between the modes of a given monomer while this is neglected in our vibrational model; and we use the CCSD(T)-F12 method to calculate the 4-body interactions whereas the LMM calculated values use the empirical TTM3-F forcefield to account for interactions beyond third-order. 20,21 As noted earlier, for OH-stretching transitions that are strongly red-shifted, both the 2and 3-body interactions have a large effect on the vibrational frequency shift. Additionally, the majority of this red-shift seems to result from the interactions in a small subset of the constituent trimers or tetramers. We have used both 2- and 3-body symmetry-adapted perturbation theory 42 to further study the relationship between energetic interactions and vibrational frequency shift for all dimers and trimers in five isomers of (H2 O)6 and those dimers and trimers containing the central water molecule of the DD*(20,1) isomer of (H2 O)21 . The 2- and 3-body interaction energies have been calculated using two different SAPT

18


Page 18 of 37

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


approaches. These interaction energies, and those determined from CCSD(T)-F12 calculations are tabulated for the prism isomer in the supplementary information. In agreement with previous work 53 we find that the 3-body interaction energy is well reproduced with modest basis set size, in this work aug-cc-pVDZ. Due to the significant contribution of the dispersion energy to the 2-body interaction energy, we find that the aug-cc-pVDZ basis set is of insufficient size to reproduce the CCSD(T)-F12 calculated interaction energies and therefore use the larger def2-ATZVPP basis set. Due to the larger basis set size required, we use the accurate and computationally efficient HF-based SAPT(2+3) approach 44 for 2-body calculations. The 2- and 3-body frequency shifts calculated using CCSD(T)-F12/VDZ-F12, Method I, and a Morse oscillator potential from all five water hexamers as well as the DD*(20,1) isomer of (H2 O)21 were fitted against a combination of SAPT energetic parameters using a Bayesian penalized regression model and the Flexmix algorithm as described in the theoretical methods section. For two body interactions the regression is performed with the calculated 2-body frequency shift as the dependent variable, and scaled SAPT-calculated electrostatics, exchange, induction and dispersion energies as independent variables. For three-body interactions the regression is performed with the calculated 3-body frequency shift as the dependent variable, and scaled SAPT calculated exchange, induction and dispersion as the independent variables. As part of the fitting routine, the Flexmix algorithm sorts the data into distinct groupings of OH-stretching frequency shift–SAPT-calculated interaction energy pairs and each group is fitted with different regression coefficients. For both 2- and 3-body fits, the geometric parameters r and ψ are used as concomitant variables that aid in Flexmix group assignment, but are not variables in the regression itself. Plots of the calculated versus the fitted frequency shifts for both 2- and 3-body are shown in Figure 3. A comparison to the calculated frequency shifts resulted in total root mean square error (RMSE) values of 5.9 cm−1 for the 2-body and 3.9 cm−1 for the 3-body vibrational frequency shifts.

19


0 −50 −100

Group Group Group Group

−150 −200 −200

−150

−100

−50

1 2 3 4

Fitted Frequency Shift / cm−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 48 49 50 51 52 53 54 55 56 57 58 59 60



0

Calculated Frequency Shift / cm−1

Page 20 of 37

100 50 Group Group Group Group Group Group

0 −50 −100 −150 −150

−100

(a)

−50

0

50

1 2 3 4 5 6 100

Calculated Frequency Shift / cm−1

(b)

Figure 3: Fitted vs. calculated many-body frequency shifts for (a) 2-body frequency shifts and (b) 3-body frequency shifts. The colors represent the data groupings as determined by the Flexmix algorithm. After obtaining converged regressions with the Flexmix algorithm, we analyzed the dimers and trimers within each Flexmix-determined group to identify geometric similarities that might provide a physical rationale to allow pre-assignment of dimers and trimers to particular groups. For both 2- and 3-body shifts, four dominant geometrical motifs emerged in the groupings where the OH-frequency shift and interaction energy could be correlated. Example dimers for three of these categories are illustrated in Figure 4. The groupings used to describe 2-body frequency shifts were: OH-bonds that are hydrogen bond donors (BD), OH-bonds where the oxygen is acting as a hydrogen bond acceptor (UA), the OH-bonds on hydrogen bond donor molecules that are not participating in the hydrogen bonding interaction (UD), and OH-bonds on unbound molecules (UU). For the ease of referencing these categories, we use a two letter key. The first letter can be either a ‘B’ or a ‘U’ indicating whether the OH of interest is bound or unbound. The second letter can be a ‘D’, ‘A’, or ‘U’ indicating that the molecule containing the OH-bond of interest acts as a hydrogen-bond donor, a hydrogen-bond acceptor, or is unbound to the other molecule in the dimer. The frequency shift-binding energy groupings are less distinct in the case of the 3-body interactions. The different geometric classes for 3-body interactions are shown in Figure 5. The most strongly red-shifted group contained trimers consisting of three water molecules in a hydrogen-bound chain (BDC). The strong red-shift of this chain system is in agreement 20


Page 21 of 37 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 the work of Ouyang and Bettens. 26 The most strongly blue-shifted group contained donor OH-bonds in double-donor and double-acceptor trimers (BD2). A third group included trimers where the frequency shift in the OH-bond of interest ranged from weakly red-shifted to weakly blue-shifted. This group consists of an OH-bond which acts as a hydrogen-bond donor in singly bound systems and has a 3-body OH-stretching vibrational frequency shift of magnitude less than ±30cm−1 . For the 3-body OH-vibrational frequency shifts, all non-bonded OH-bonds can all be classified as “non-shifted” frequencies (N). This includes OH-bonds on acceptor molecules, unbound OH-bonds on bound waters, and OHbonds on molecules with no hydrogen bonding interactions. Although these categories do have unique red- or blue- shifting patterns, all of their corresponding 2-body frequency shifts are less than 9 cm−1 . For this reason, they can be correlated to the SAPT-calculated interaction energies as a single group without any significant loss of accuracy.

Figure 4: Example water dimer structures that belong to different geometrical groupings. From left to right the OH-bond including the grey hydrogen are from the categories: BD (hydrogen-bound OH-bond on a donor molecule), UD (an unbound OH-bond on a donor molecule), and UA (an unbound OH-bond on an acceptor molecule). The UU category is not pictured, but corresponds to an OH-bond on a water molecule which is not hydrogen bound to another. In Figure 6 we plot the many-body frequency shift against SAPT calculated electrostatics (2-body) and induction (3-body). The data points are color coded to denote which geometric group they are assigned to. As can be seen, there is a distinct correlation between group identity and the relationship between interaction energy and OH-stretching vibrational frequency shift. The frequency shift and interaction energy data for each group were then fitted independently using penalized regression models and the results are shown in Figure 7. As

21



0

−100

BD UA UD UU

−200 −15

−10

−5


Figure 5: Examples of water trimers that belong to different geometrical groupings. From left to right, the OH-bonds including the grey hydrogen are from the categories BDC (hydrogenbound OH-bonds in a chain system), and BD2 (one double acceptor and one double donor system shown). Two classifications are not pictured here. The BD category is the same as that in Figure 4 with a third, unbound molecule. The N category consists of multiple geometric patterns all with near-zero frequency shifts.


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

Page 22 of 37

0

Electrostatics / kcal mol−1

100 50 0 −50

BDC BD2 BD N

−100 −150

−1.5

(a)

−1

−0.5

0

0.5

1

Induction / kcal mol−1

(b)

Figure 6: (a) 2-body OH-stretching frequency shift vs. SAPT-calculated electrostatics (b) 3-body OH-stretching frequency shift vs. SAPT-calculated Induction

22


Page 23 of 37

was the case when using Flexmix determined groupings, the regressions were performed using the calculating n-body frequency shifts as the dependent variable, and the SAPT-calculated many-body energies as the independent variables. The quality of these linear regressions are comparable to those obtained using the Flexmix algorithm. In Table 4, the RMSE for the Flexmix model and our pre-grouped fits are compared. It can be seen that the RMSE values for both approaches are comparable and on the order of approximately 5 cm−1 for the complete set of dimers and trimers. Within the geometrical structure groupings, the largest RMSE is observed for the most strongly shifted 2- and 3-body frequencies with maximum RMSE of around 15 cm−1 and 11 cm−1 respectively. The 2- and 3-body frequency shifts for these groups are approximately 150 cm−1 and 100 cm−1 respectively, indicating that our regression model results in frequency shifts that differ from the calculated value by at most

0 −50 −100

UU UA UD BD

−150 −200 −200

−150

−100

−50

Calculated Frequency Shift / cm


around 10%. Fitted Frequency Shift / cm−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 48 49 50 51 52 53 54 55 56 57 58 59 60


0 −1

100 50 0 −50

N BD2 BD BDC

−100 −150 −150

−100

−50

0

50

Calculated Frequency Shift / cm

(a)

100

−1

(b)

Figure 7: Fitted vs. calculated many-body frequency shifts for (a) 2-body frequency shifts and (b) 3-body frequency shifts. The regression was performed with a penalized regression model.

23



Page 24 of 37

Table 4: Root mean square errors for penalized regressions of many-body OHstretching frequency shifts vs SAPT-calculated energies 2-Body

3-Body

Categoriesa

Flexmix

Pre-Grouped

Categories

Flexmix

Pre-Grouped

Total

5.9

7.2

Total

3.9

6.3

Red (BD)

15.4

15.3

Red (BDC)

10.7

10.8

Blue (UD)

2.3

2.7

Blue (BD2)

4.5

9.9

Black (UU)

3.8

5.8

Black (N)

0.3

1.9

Green (UA)

2.3

4.9

Green (BD)

1.4

2.8

Light Blue

2.7

Grey

1.3

a

Note that the groups are listed by color with the corresponding group identity in parenthesis since the groups assigned by the Flexmix algorithm do not perfectly correspond with our pre-sorted group assignments. To allow for comparison, the Flexmix and pre-sorted groups groups that have the largest overlap in group population are aligned.

A challenge in determining which interaction energy component is responsible for the variation of OH-stretching frequencies is the strong degree of covariation in the individual SAPT-calculated energetic components. The correlation values between the individual SAPT components and the total energy are shown in Table 5. The Spearman correlation coefficient 54 allows for the degree of correlation to be identified even in instances where the correlation between variables is non-linear. 55 As evidenced by the coefficients, there is a strong degree of correlation between the energetic parameters, particularly for the 2-body interaction energy. In the case of the 3-body energy; there is a strong correlation between the induction energy and the total energy. This is not surprising given that for the systems we study, the induction energy comprises 111% of the 3-body energy. This value is over 100% since there are other energetic components which are of opposite sign to the induction, thus allowing the induction energy to be larger than the total energy on average. For comparison

24


Page 25 of 37 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


sake, the exchange and dispersion energy are 40% and 24% of the total energy respectively.

Table 5: Spearman correlation coefficients for the relationship between the total 2- and 3-body SAPT energy and individual SAPT calculated components.

2-body Total 3-body Total

Electrostatics 0.93 —

Exchange -0.80 0.60

Induction 0.81 0.996

Dispersion 0.80 -0.45

Although the calculated OH-stretching frequencies can be correlated when all SAPTcalculated energetic parameters are treated as independent variables, we carried out a linear regression of the calculated frequency shifts against the single components of the SAPT calculated interaction energy. In Table 6, RMSE values obtained from this approach are summarized. As anticipated, the RMSE values are larger than those in Table 4 due to the reduced number of independent variables. Furthermore, an additional restriction was placed on these models by constraining the regression to a zero intercept. This zero intercept was used in order to allow for a direct comparison of the regression coefficients obtained from fits to the data in each of the pre-sorted groups. Despite the additional constraints, the RMSE obtained from a fit of 2-body OH-stretching frequency shifts against 2-body electrostatics is comparable (within 3 cm−1 ) to the values obtained in Table 4 where all SAPT-calculated energies were considered as independent variables. This level of agreement is obtained for all four geometrical groupings considered. Similarly, the RMSE of the model wherein the 3-body frequency shifts were correlated to only the 3-body induction energy is within 7 cm−1 of the RMSE obtained when all interaction energy components were treated as independent variables.

25



Page 26 of 37

Table 6: Root mean square errors for linear regressions of SAPT calculated energies vs many-body OH-stretching frequency shifts Parameter

2-body

3-body

BD

UA

UU

UD

BDC

BD

N

BD2

All (scaled)

15.3

4.9

5.8

2.7

10.8

2.8

1.9

9.9

Electrostatics

17.8

3.5

5.8

5.3

—

—

—

—

Induction

24.3

4.1

6.5

5.8

14.4

7.7

1.9

16.2

Exchange

23.5

3.9

6.8

5.5

16.4

14.4

1.9

32.8

Dispersion

20.9

3.3

6.5

5.1

26.4

15.0

1.9

21.2

The smallest RMSE values for the 2-body shifts are, perhaps unsurprisingly, obtained from a fit against electrostatics. For 3-body shifts, the lowest RMSE is obtained when fitting 3-body shifts against the induction energy. The coefficients for the single-parameter fits are summarized in Table 7. Interestingly, the induction coefficient for three out of the four 3-body groups is approximately 70 cm−1 (kcal · mol−1 )−1 indicating a consistent relationship between OH-stretching frequency shift and the induction energy regardless of geometric arrangement. The only group with a different coefficient is the non-shifted group, which, as noted earlier, has almost negligible frequency shifts. The 3-body induction energy is the only component of the interaction energy which changed from a negative to a positive energetic interaction in conjunction with the shift from red- to blue-shift transitions. This allows for a single correlation coefficient to describe the relationship between OH-stretching frequency shift and 3-body induction energy for all trimers that exhibit frequency shifts of greater than 9 cm−1 as shown in Figure 8. This change of sign in induction indicates that the change from red- to blue-shifting in 3-body is primarily caused by the change in induction energy from stabilizing to destabilizing. In contrast, the 3-body exchange energy stays consistently negative and 3-body dispersion energy stays consistently positive regardless of the red or blue shift of the transition. As a result, the three-body exchange and dispersion energies 26


Page 27 of 37 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


are not correlated to the shift in vibrational frequency in the same manner and therefore we require a different regression coefficient for each group in the instance where they are used as the sole independent variable. Table 7: Linear regression coefficients (cm−1 (kcal mol−1 )−1 ) for single parameter SAPT energy–frequency shift models 2-body

Electrostatics

3-body

Induction

BD

13.7

BDC

72.0

UA

-0.724

BD

65.3

UU

3.87

BD2

71.6

UD

1.15

N

0.268

The lowest RMSE value when using a single independent variable to model the 2-body frequency shift is obtained when using electrostatic energy as the independent variable. However, in contrast to 3-body shifts, there is no single coefficient that can be used to fit the 2-body shifts of different groups because of the drastically different geometrical nature of these groups. In the 3-body fit, all three groups with significant frequency shifts can be classified as “hydrogen-bond donors” which contain the dimers found in the single 2-body BD group. The frequency shifts for OH-bonds where the oxygen acts as a hydrogen bond acceptor and for OH-bonds involving the non-bound hydrogens in a hydrogen-bound dimer exhibit different relationships to the energetic data as evidenced by the coefficients shown in Table 7. As noted earlier, in the case of 3-body frequency shifts, all non-bound trimers have very small 3-body frequency shifts, but 2-body frequency shifts can be sizable even for non-hydrogen bound dimers. Finally, it is important to note that the current pre-sorted model does not explicitly sub-categorize OH-bonds on the basis of the identity of their hydrogen bond geometry. In other words, distance and angle parameters are used in equation 6 to determine whether a

27


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

Page 28 of 37

BDC BD2 BD

50 0 −50 −100 −150

−1.5

−1

−0.5

0

0.5

1


Figure 8: Linear regression fit for 3-body OH-stretching frequency shifts. The fitted line has a slope of 70 cm−1 (kcal mol−1 )−1 and is constrained to have a zero intercept. dimer is hydrogen-bound or not, but we do not make any distinction between the average bonding angles of 167° observed for the bag isomer and the average bonding angles of 153° in the prism isomer. Although we consider 380 OH-stretching frequencies calculated for 95 different water dimers, for strongly redshifted transitions our trial set contains more dimers with OH–O angles close to the idealized 180° than ‘bent’ dimers with non-linear hydrogen bonding geometries. As a result, the frequencies of the OH-stretching modes in the bag and prism isomers were predicted by pre-grouped, penalized regression models with an RMSE of 6.97 cm−1 and 14.65 cm−1 respectively. In agreement with Lenz et al. we observe that the OH-stretching frequency shifts of bent hydrogen bonds appear to be less than those of non-bent hydrogen-bonds. 24 This suggests that the dimers with hydrogen bonding interactions with OH–O angles far from 180° exhibit a different correlation between 2-body frequency shift and the SAPT-calculated energy components than do dimers with more idealized hydrogen bonded geometries. To correlate the 2-body frequency shifts and SAPT-calculated interaction energies in bent water dimers may require the use of more geometric groupings than were considered in this study. 28


Page 29 of 37 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


Conclusion We have presented two different variants of a many-body fragmentation scheme suitable for calculating the local-mode vibrational frequencies of strongly hydrogen-bound systems. Within this framework, we found that the 1- through 4-body interactions are needed to accurately describe the frequency shift for a hydrogen-bound OH-stretching vibration. Furthermore, we showed that the 3-body contribution to the vibrational frequency shift is often as large or larger than the 2-body contribution. This was shown to be caused by the large inductive interactions present when a third water molecule is placed near the dimer containing the hydrogen bond. The comparison of various linear regression models used to identify the dominant effects associated with 2- and 3-body frequency shifts predicts that the 2-body shift can be explained by the magnitude of the electrostatic interactions, while the 3-body shifts can be well correlated to the induction energy. We demonstrate that 3-body chain systems have large stabilizing inductive energies and correspondingly large negative frequency shifts, while 3-body double-donor and double-acceptor systems were shown to have destabilizing inductive energies and correspondingly positive frequency shifts. The discovery that a given vibrational frequency shift depends on only a very small number of the total trimers or tetramers allows a further approximation which greatly reduces the computational cost of employing a many-body aproach to calculating vibrational frequency shifts. Thus, we predict that this fragmentation scheme could be used to calculate the vibrational frequencies of molecules on the fly in fairly large systems.

Supporting Information Description Cartesian coordinates for all optimized structures. 2-body energies calculated with CCSD(T)F12 and SAPT(2+3), and 3-body energies calculated with CCSD(T)-F12 and DFT-SAPT for the prism isomer of (H2 O)6 . Calculated n-body frequency shifts for the bag, book, prism, 29



ring, and DD*(20,1) clusters. Table describing groupings of SAPT 2- and 3-body terms. Tabulated comparison of full and truncated many-body expansions for bag, book, cage, prism, and ring clusters.

Acknowledgements The authors acknowledge financial support from the M. J. Murdock Charitable Trust.

References (1) Reddy, S. K.; Straight, S. C.; Bajaj, P.; Huy Pham, C.; Riera, M., Moberg, D. R.; Morales, M. A.; Knight, C.; Gotz, A. W.; Paesani, F. On the Accuracy of the MB-pol Many-body Potential for Water: Interaction Energies, Vibrational Frequencies, and Classical Thermodynamic and Dynamical Properties from Clusters to Liquid Water and Ice. J. Chem. Phys. 2016, 145, 194504. (2) Medders, G. R.; Paesani, F. Many-body Convergence of the Electrostatic Properties of Water. J. Chem. Theory Comput. 2013, 9, 4844-4852. (3) Cui, J.; Liu, H.; Jordan, K. D. Theoretical Characterization of the (H2 O)21 Cluster: Application of an n-body Decomposition Procedure. J. Phys. Chem. B. 2006, 110, 18872-18878. (4) Bates, D. M.; Smith, J. R.; Janowski, T.; Tschumper, G. S. Development of a 3body: Many-body Integrated Fragmentation Method for Weakly Bound Clusters and Application to Water Clusters (H2 O)n=3,10,16,17 . J. Chem. Phys. 2011, 135, 044123. (5) Dahlke, E. E.; Truhlar, D. G. Electrostatically Embedded Many-body Expansion for Large Systems, with Applications to Water Clusters. J. Chem Theory Comput. 2007, 3, 46-53.

30


Page 30 of 37

Page 31 of 37 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) Chen, W.; Gordon, M. S. Energy Decomposition Analyses for Many-body Interaction and Applications to Water Complexes. J. Phys. Chem. 1996, 100, 14316-14328. (7) Li, W.; Li, S.; Jiang, Y. Generalized Energy-based Fragmentation Approach for Computing the Ground-state Energies and Properties of Large Molecules. J. Phys. Chem. A. 2007, 111, 2193-2199. (8) 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. (9) 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. (10) 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. (11) Mayhall, N. J.; Raghavachari, K. Many-overlapping-body (MOB) Expansion: A Generalized Many Body Expansion for Nondisjoint Monomers in Molecular Fragmentation Calculations of Covalent molecules. J. Chem. Theory Comput. 2012, 8, 2669-2675. (12) Mayhall, N. J.; Raghavachari, K. Molecules-in-molecules: An Extrapolated Fragmentbased Approach for Accurate Calculations on Large Molecules and Materials. J. Chem. Theory Comput. 2011, 7, 1336-1343. (13) Elsohly, A. M.; Shaw, C. L.; Guice, M. E.; Smith, B. D.; Tschumper, G. S. Analytic Gradients for the Multicentred Integrated QM: QM Method for Weakly Bound Clusters: Efficient and Accurate 2-body: Many-body Geometry Optimizations. Molec. Phys.. 2007, 105, 2777-2782. 31



(14) Bates, D. M.; Smith, J. R.; Tschumper, G. S. Efficient and Accurate Methods for the Geometry Optimization of Water Clusters: Application of Analytic Gradients for the Two-body: Many-body QM: QM Fragmentation Method to (H2 O)n , n= 3-10. J. Chem. Theory Comput. 2011, 7, 2753-2760. (15) Howard, J. C.; Tschumper, G. S. N-body: Many-body QM: QM Vibrational Frequencies: Application to Small Hydrogen-bonded Clusters. J. Chem. Phys. 2013, 139, 184113. (16) Pedulla, J. M.; Kim, K.; Jordan, K. D. Theoretical Study of the n-body Interaction Energies of the Ring, Cage and Prism Forms of (H2 O)6 . Chem. Phys. Lett. 1998, 291, 78-84. (17) Hodges, M. P.; Stone, A. J.; Xantheas, S. S. Contribution of Many-body Terms to the Energy for Small Water Clusters: A Comparison of Ab Initio Calculations and Accurate Model Potentials. J. Phys. Chem. A. 1997, 101, 9163-9168. (18) Xantheas, S. S. Ab Initio Studies of Cyclic Water Clusters (H2 O)n , n= 1–6. II. Analysis of Many-body Interactions. J. Chem. Phys. 1994, 100, 7523-7534. (19) Ouyang, J. F.; Bettens, R. P. Many-Body Basis Set Superposition Effect. J. Chem. Theory Comput. 2015, 11, 5132-5143. (20) Wang, Y.; Bowman, J. M. Ab Initio Potential and Dipole Moment Surfaces for Water. II. Local-monomer Calculations of the Infrared Spectra of Water Clusters. J. Chem. Phys. 2011, 134, 154510. (21) Wang, Y.; Huang, X.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Flexible, Ab Initio Potential, and Dipole Moment Surfaces for Water. I. Tests and Applications for Clusters up to the 22-mer. J. Chem. Phys. 2011, 134, 094509.

32


Page 32 of 37

Page 33 of 37 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


(22) Wang, Y.; Bowman, J.M. IR Spectra of the Water Hexamer: Theory, with Inclusion of the Monomer Bend Overtone, and Experiment are in Agreement. J. Phys. Chem. Lett. 2013, 4, 1104-1108. (23) Steinbach, C.; Andersson, P.; Melzer, M.; Kazimirski, J. K.; Buck, U.; Buch, V. Detection of the Book Isomer from the OH-stretch Spectroscopy of Size Selected Water Hexamers. J. Phys. Chem. Chem. Phys. 2004, 6, 3320-3324. (24) Auer, B.; Kumar, R.; Schmidt J. R.; Skinner, J. L. Hydrogen Bonding and Raman, IR, and 2D-IR Spectroscopy of Dilute HOD in liquid D2 O. PNAS. 2007, 104, 14215-14220. (25) Lenz, A.; Ojamäe, L. Theoretical IR Spectra for Water Clusters (H2 O)n (n=6-22,28,30) and Identification of Spectral Contributions from Different H-bond Conformations in Gaseous and Liquid Water. J. Phys. Chem. 2006, 110, 13388-13393 (26) Ouyang, J. F.; Bettens, R. P. When are Many-Body Effects Significant? J. Chem. Theory and Comp. 2016, 12, 5860-5867. (27) Kjaergaard, H. G.; Howard, D. L.; Schofield, D. P.; Robinson, T. W. OH- and CHStretching Overtone Spectra of Catechol. J. Phys. Chem. 2002, 106, 258266. (28) Schofield, D. P.; Kjaergaard, H. G. Calculated OH-stretching and HOH-bending Vibrational Transitions in the Water Dimer. J. Phys. Chem. Chem. Phys. 2003, 5, 3100-3105. (29) Ishiuchi, S.; Robinson, T. W.; Miller, B. J.; Kjaergaard, H. G. Vibrational Overtone Spectroscopy of Phenol and Its Deuterated Isotopomers. J. Phys. Chem. 2006, 110, 7345-7354. (30) Schofield, D. P.; Lane, J. R.; Kjaergaard, H. G. Hydrogen Bonded OH-stretching Vibration in the Water Dimer. J. Phys. Chem. A. 2007, 111, 567-572. (31) Schofield, D. P.; Kjaergaard, H. G.; Matthews, J.; Sinha, A. The OH-stretching and OOH-bending Overtone Spectrum of HOONO. J. Chem. Phys. 2005, 123, 134318. 33



(32) Atkins, P. W.; Friedman, R. S. Molecular Quantum Mechanics, 3rd edn.; Oxford University Press: Oxford, 1997. (33) Boys, S.F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Molec. Phys. 1970, 19, 553-566. (34) Mierzwicki, K.; Latajka, Z. Basis Set Superposition Error in N-body Clusters. Chem. Phys. Lett. 2003, 380, 654-664. (35) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297-3305. (36) Peterson, K. A.; Adler, T. B., Werner, H.-J. Systemically Convergent Basis Sets for Explicitly Correlated Wavefunctions: The Atoms H, He, B-Ne, and Al-Ar. J. Chem. Phys. 2008, 128, 084102. (37) Werner, H.-J.; Knowles, P.J.; Knozia, G.; Manby, F. R.; Sch¨ utz, M. Molpro: A GeneralPurpose Quantum Chemistry Program Package. WIREs Comput. Mol. Sci. 2012, 2, 242-253. (38) Werner, H.-J.; Knowles, P.J.; Knozia, G.; Manby, F. R.; Sch¨ utz, M.; et al. MOLPRO, version 2015.1, a Package of ab initio Programs. 2015. (39) Adler, T. B.; Knizia, G.; Werner, H.-J. A Simple and Efficient CCSD(T)-F12 Approximation. J. Chem. Phys. 2007, 127, 221106. (40) Knizia, G; Adler, T. B.; Werner, H.-J. Simplified CCSD(T)-F12 methods: Theory and Benchmarks. J. Chem. Phys. 2009, 130, 054104. (41) McMahon, J. D.; Lane, J. R. Explicit Correlation and Basis Set Superposition Error: The Structure and Energy of Carbon Dioxide Dimer. J. Chem. Phys. 2011, 135, 154309. 34


Page 34 of 37

Page 35 of 37 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


(42) Szalewicz, K. Symmetry-adapted Perturbation Theory of Intermolecular Forces. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 254-272. (43) Hohenstein, E. G.; Sherrill, C. D. Wavefunction Methods for Noncovalent Interactions. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 304-326. (44) Hohenstein, E. G.; Sherrill, C. D. Density Fitting of Intramonomer Correlation Effects in Symmetry-adapted Perturbation Theory. J. Chem. Phys. 2010, 133, 014101. (45) Hohenstein, E. G.; Sherrill, C. D. Efficient Evaluation of Triple Excitations in Symmetry-adapted Perturbation Theory via Second-order Møller–Plesset Perturbation Theory Natural Orbitals. J. Chem. Phys. 2010, 133, 104107. (46) Podeszwa, R.; Szalewicz, K. Three-body Symmetry-adapted Perturbation Theory Based on Kohn-Sham Description of the Monomers. J. Chem. Phys. 2007, 126, 194101 (47) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriana, S.; Dahle, P. et al. The Dalton Quantum Chemistry Program System. WIREs Comput. Mol. Sci. 2014, 4, 269-284. (48) Dalton, a molecular electronic structure program, Release Dalton2.0 2005 (49) Gr¨ un, B.; Leisch, F. FlexMix Version 2: Finite Mixtures with Concomitant Variables and Varying and Constant Parameters. J. Stat. Softw. 2008, 28, 1-34. (50) Mortier, F.; Ouédraogo; D.-Y.; Claeys, F.; Tadesse, M. G.; Cornu, G.; Baya, F.; Benedet, F.; Freycon, V.; Gourlet-Fleury, S; Picard, N. Mixture of Inhomogeneous Matrix Models for Species-rich Ecosystems. Environmetrics 2015, 26, 39-51 (51) Kumar, R.; Schmidt, J. R.; Skinner, J. L. Hydrogen Bonding Definitions and Dynamics in Liquid Water. J. Chem. Phys. 2007, 126, 204107. (52) Miliordos, E.; Aprà, E.; Xantheas, S. S. Optimal Geometries and Harmonic Vibrational Frequencies of the Global Minima of Water Clusters (H2 O)n , n= 2-6, and Several 35



Hexamer Local Minima at the CCSD (T) Level of Theory. J. Chem. Phys. 2013, 139, 114302. (53) 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. (54) Savicky, P. (2014). pspearman: Spearman’s Rank Correlation Test. R package version 0.3-0. (55) Artusi, R.; Verderio, P.; Marubini E. Bravais-Pearson and Spearman Correlation Coefficients: Meaning, Test of Hypothesis and Confidence Interval. The Int. J. Bio. Markers. 2002, 17, 148-151.

36


Page 36 of 37

Page 37 of 37

TOC Graphic Fitted Frequency Shift / cm−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 48 49 50 51 52 53 54 55 56 57 58 59 60


100

BDC BD2 BD

50

BD

0 BDC

−50 −100 −150

BD2

−1.5

−1

−0.5

0

0.5

1


37


Origin of Many-Body Vibrational Frequency Shifts in Water Clusters

Recommend Documents