Subscriber access provided by UNIVERSITY OF THE SUNSHINE COAST
Statistical Mechanics -
2
2
2
-
2
2
Anharmonic Densities of States for Vibrationally Excited I(HO), (HO), and I(HO) Xinyou Ma, Nan Yang, Mark A. Johnson, and William Louis Hase J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00300 • Publication Date (Web): 26 Jun 2018 Downloaded from http://pubs.acs.org on June 27, 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 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
Journal of Chemical Theory and Computation
Anharmonic Densities of States for Vibrationally Excited I-(H2O), (H2O)2, and I-(H2O)2 Xinyou Ma,a Nan Yang,b Mark A. Johnson,b and William L. Hase*,a a
Department of Chemistry and Biochemistry Texas Tech University Lubbock, Texas 79409, USA b
Sterling Chemistry Laboratory Department of Chemistry Yale University
New Haven, Connecticut 06520, USA
1 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 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 2 of 38
Abstract
Monte Carlo sampling calculations were performed to determine the anharmonic sum of states, Nanh(E), for I-(H2O), (H2O)2, and I-(H2O)2 versus internal energy up to their dissociation energies. The anharmonic density of states, ρanh(E), is found from the energy derivative of Nanh(E). Analytic potential energy functions are used for the calculations, consisting of TIP4P for H2O---H2O interactions and an accurate 2-body potential for I----H2O fit to quantum chemical calculations. The extensive Monte Carlo samplings are computationally demanding and the use of computationally efficient potentials was essential for the calculations. Particular emphasis is directed towards I-(H2O)2 and distributions of its structures versus internal energy are consistent with experimental studies of the temperature-dependent vibrational spectra. At their dissociation thresholds, the anharmonic to harmonic density of states ratio, ρanh(E)/ρh(E), is ~2, ~3, and ~260 for I-(H2O), (H2O)2, and I-(H2O)2, respectively. The large ratio for I-(H2O)2 results from the I(H2O)2 → I-(H2O) + H2O dissociation energy being more than two times larger than the (H2O)2 → 2H2O dissociation energy, giving rise to highly mobile H2O molecules near the I-(H2O)2 dissociation threshold. The work illustrates the importance of treating anharmonicity correctly in unimolecular rate constant calculations.
2 ACS Paragon Plus Environment
Page 3 of 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
Journal of Chemical Theory and Computation
I. Introduction Understanding structures and dynamics of hydrated anions, such as I-(H2O)n, from first principles is an important challenge for contemporary theoretical methods designed to predict the behavior of aqueous electrolytes.1-8 Even within the cluster regime, microsolvation has a dramatic effect on the reaction kinetics, with a recent example provided by the OH-(H2O)n + CH3Br and OH-(H2O)n + CH3I, SN2 reactions.9-11 Over the past twenty years, the experimental structures of small I-(H2O)1-4 structures have been determined by comparison of the linear vibrational spectra of the cold ions (isolated in the gas phase)12,13 with calculated geometries for minimum energy structures obtained with high level electronic structure calculations.1-8,14 Although the vibrational zero-point structures are established with this approach, simulation of realistic environments requires understanding large amplitude motions on an accurate potential energy surface (PES). In this context the I-(H2O)2 system, with the global minimum structure displayed in (Figure 1), has played an important role as vibrational spectra have been reported as a function temperature over the range T = 10-200 K.15 These spectra display a dramatic change at 125 K, which was interpreted in the context of the break-up of the intra-cluster water dimer into two essentially independent water molecules orbiting the large ion core.
This qualitative
suggestion raises the issue of how this ternary system behaves as a function of internal energy on a realistic PES. In this paper, we report results of a computational study exploring the accurate density of states, ρ(E), of I-(H2O)2 versus vibrational energy and the onset large amplitude motions. To compare with this study, ρ(E) is also calculated for I-(H2O) and (H2O)2. These ρ(E) are determined numerically using classical statistical mechanics, which is expected to be accurate for the “soft” forces and low vibrational frequencies for the intermolecular modes of these clusters. The calculations yield a picture of the populations of various configurations of I-(H2O)2 as a function E, and indicates a correction to the harmonic ρ(E) on the order of 103 near the I-(H2O)2 → I-(H2O) + (H2O)2 threshold. Harmonic RRKM theory is not applicable for this dissociation.
II. Classical Phase Space Volume and Density of States A. Monte Carlo Sampling The classical density of states ρ(E) is the energy derivative of the classical sum of states N(E), ρ(E) = dρ(E)/dE. N(E) is the phase space volume at E, PSV(E), divided by hs, where h is 3 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 38
Planck’s constant and s is the number of degrees of freedom. For a one-dimension system, such as A-B stretching, PSV(E) may be obtained by an analytic or numerical solution of the cyclic
integral ∮pdq. A straightforward extension of this approach is applicable for s separable degrees
of freedom, where PSV(E) is evaluated for each degree of freedom and their volumes convoluted. However, for a multi-dimensional, non-separable Hamiltonian, this approach is not tractable and the following multi-dimensional integral
=
⋯ d d
(1)
must be evaluated, where H is the Hamiltonian H(p, q).17 For the work reported here, this integral is solved by Monte Carlo sampling, for which the phase space volume is determined by random sampling configurations within the phase space (p, q); i.e.
=
⋯ d d
(2)
where Nsample is the total number sampled phase space configurations, NH≤E is the total count of
configurations within the hypersurface boundary, and the integration ⋯ d d is the volume
of a hyperrectangle, which contains the complete phase space manifold. For arbitrary chosen
boundaries of the hyperrectangle, the Monte Carlo sampling acceptance rate NH≤E/Nsample decreases for enlarged boundaries. Thus, the optimal hyperrectangle boundaries lie on the phase space boundary of each correspondent dimension. Additionally, if the phase space is symmetric for a dimension, it is usually sufficient to sample only one half of this dimension, and then a factor of 2 should be multiplied to compensate for the total phase space volume. For example, only one half of the H–H vibration phase space is necessary for sampling a Morse potential, and only a quarter is necessary for sampling a harmonic potential. Pseudo random numbers, following a uniform distribution are often used in Monte Carlo sampling for statistical accuracy.18-20 If the acceptance rate is extremely low, the sample size Nsample needs to be quite large, and, thus, the period of the pseudo random number generator
4 ACS Paragon Plus Environment
Page 5 of 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
Journal of Chemical Theory and Computation
(PRNG) must be much greater than the sample size Nsample to avoid repeatedly reusing pseudo random numbers. The Mersenne Twister (MT19937) pseudo random number generator proposed by Matsumoto and Nishumira21,22 is used here. It has a long period of 219937 – 1 and shows good uniformity in multi-dimensional sampling.18 Since the samplings are independent, Monte Carlo sampling may be implemented with a parallel algorithm for shorter wall-time execution and high parallel scalability. Here, message passing interface (MPI)23 and Nvidia’s Compute Unified Device Architecture (CUDA)24 with the Thrust library25 were implemented in the CPU and GPU codes, respectively. B. Accuracy of the Classical Density of States Given the large density of states ρ(E) for I-(H2O)2 the only practical way to calculate ρ(E) is by classical statistical mechanics, as described above. Nevertheless, it is important to identify the accuracy of the classical ρ(E) and this is illustrated in the Supporting Information for the normal mode harmonic oscillator models for I-(H2O)2, (H2O)2, and I-(H2O)2. The classical and quantum ρ(E) may be compared by considering the total energy Etotal which is a sum of the zero-point energy (ZPE) and the energy in excess of the ZPE.26,27 The classical ρ(E) has values beginning at Etotal = 0, while the quantum ρ(E) begins at the ZPE. As shown in the Supporting Information, the normal mode harmonic oscillator classical and quantum ρ(E) for I-(H2O) are in near exact agreement for all Etotal in excess of the ZPE. For (H2O)2 the two ρ(E) differ by a ratio of 2.7 at the ZPE and only 1.2 at its dissociation threshold. For I-(H2O)2 the two ρ(E) differ by 8 at the ZPE, but only 1.2 at its dissociation threshold. In addition, as shown in the Supporting Information, the normal mode harmonic oscillator classical and quantum partition functions with respect to the classical potential energy minimum are in excellent agreement for all three species; i.e. these respective partition functions are kBT/hν and exp[-hν/(2kBT)]/[1 – exp(-hν/kBT)]. III. Potential Energy Functions The potential energy function for I-(H2O)2 was written as a sum of two I----H2O potentials and a H2O---H2O potential. The Monte Carlo sampling described above is computationally demanding and, to facilitate the calculations, the potentials used are written as sums of two-body terms. The TIP4P model was used to represent the H2O---H2O potential.28 A number of 2-body2934
and many-body35-43 analytic potential energy functions have been developed for I-(H2O). For
the calculations reported here, a computationally efficient I-(H2O) analytic intermolecular 5 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 38
potential energy function was needed for the Monte Carlo samplings and was developed using DFT44, with the B97-1 functional45,46 and the ECP/d basis set.47 With H2O held rigid in its minimum energy geometry, i.e. r(O-H) = 0.964 Å and θ(HOH) = 104.4o, potential energy curves were calculated for I- approaching H2O in the five configurations shown in Figure 2. The dissociation enthalpy of I-(H2O) → I- + H2O predicted by B97-1/ECP/d is 10.59 kcal/mol, which agrees with the experimental value of 10.48 ± 0.14 kcal/mol.48,49 A sum of two-body analytic potentials was used to simultaneously fit the B97-1/ECP/d I---(H2O) potential energy curves in Figure 2. Each two-body potential consists of two generalized Lennard-Jones terms and the Buckingham exponential function, i.e.:
= +
+
!"
+
# + $%
(3)
where r is the interatomic distance, and A, B, C, D, c, d, m, and n are parameters determined by the fitting.50 To improve the fitting accuracy, a TIP4P type ghost atom (G) is introduced, which lies on the H-O-H angle bi-sector with a distance z from the oxygen atom. Thus, the I----(H2O) intermolecular potential consists of four two-body potential energy functions: two hydrogeniodine potentials (VHI), an oxygen-iodine potential (VOI), and a ghost atom-iodine potential (VGI), which respectively depend on the interatomic distances r(H1–I), r(H2–I), r(O–I) and r(G–I) as shown in Eq. (4).
&'() * = +
,
-./
(& ()& + *& *& + 0& 0&
(4)
A Particle Swarm Optimization (PSO) algorithm50 was used for fitting the I----(H2O) intermolecular potential energy curves, and the fitting results are shown in Figure 2. The fitted parameters are listed in Table 1.
IV. Anharmonic Density of States The total classical anharmonic sum and density of states for I-(H2O), (H2O)2, and I-(H2O)2 were calculated as described above. These quantities include all the equivalent potential energy minima for the species which are 2, 8, and 16 for I-(H2O), (H2O)2, and I-(H2O)2, respectively. 6 ACS Paragon Plus Environment
Page 7 of 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
Journal of Chemical Theory and Computation
Harmonic sums and densities of states were determined to compare with the anharmonic values. The harmonic quantities are the sum or density of states for one minima, multiplied by the number of equivalent minima. For all three species the ρanh(E)/ρh(E) approached 1 as E → 0. Symmetry numbers51,52 were not included in the analyses. A. I-(H2O) 1. Hamiltonian and sampling algorithm The water molecule is represented as a rigid asymmetric top rotor and the Hamiltonian for the three intermolecular modes of I-(H2O) is then
1=
, 6/,8 /, 4/, + + + + &'() * 23/ 23/ 5/, 298
(5)
8.:,;,