On the calculation of the dissociation rate constant of the water dimer by the
ab initio anharmonic RRKM theory
L. Yao
a,b,c,*, R.X. He
d, A.M. Mebel
e, S.H. Lin
b,c aDepartment of Physics, Dalian Maritime University, Dalian 116023, China
bInstitute of Atomic and Molecular Sciences, Academia Sinica, P.O. Box 23-166, Taipei 10764, Taiwan cDepartment of Applied Chemistry, National Chiao-Tung University, Hsin-chu, Taiwan
d
School of Chemistry and Chemical Engineering, Southwest University, Chongqing, China e
Department of Chemistry and Biochemistry, Florida International University, Miami, FL 33199, USA
a r t i c l e
i n f o
Article history:
Received 24 November 2008 In final form 27 January 2009 Available online 3 February 2009
a b s t r a c t
The dissociation rate constant of the water dimer was calculated by using the method proposed by Yao and Lin (YL method). The dividing surface method and RRKM theory are used to calculate the pseudo transition state and the rate constant, respectively. For the canonical case at temperature range of 243–1000 K and for the microcanonical system at total energy 1411–4000 cm1, the anharmonic rate constant, around 109–1011s1, is close to the experimental prediction. The present results indicate that the anharmonic effect should be included and YL method is suitable for calculating dissociation rate con-stants of such small flexible water cluster.
Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction
Due to the ubiquity of water clusters in nature, water cluster kinetics has been a focus of extensive research. The clustering of water molecules is relevant to a variety of physical and chemical processes, which are important in many industrial and environ-mental processes, including atmospheric aerosol formation, steam erosion of turbines, and atmospheric absorption of solar radiation. Long-lived water clusters composed of up to a few water molecules are known to exist in non-negligible amounts in water vapor and in the atmosphere[1].
Far-infrared vibration–rotation-tunneling spectroscopy tech-niques have been used to study dimers in experiments. Infrared (IR) spectroscopy has been also used to study larger water clusters. Lee and co-workers[1]predict an upper limit to the excited vibra-tional state lifetime, 1
l
s. The dissociation (evaporation) rate con-stant of water dimer was predicted to be around 109s1[1].The water dimer system has received a great amount of theoret-ical interest as it provides an important benchmark for further studies of the dynamics of water evaporation. The far-IR spectra of the water dimer, tetramer, and hexamer had been calculated by employing DFT-based MD simulations at two different temper-atures corresponding to the experimental measurements by Lee and co-workers[1]. Garrett and co-workers[2]reported the evap-oration rate constants for the water dimer to be around 1011s1of
at T = 243 K. They proposed a theoretical method called vapor-phase dynamical nucleation theory (DNT) and applied this ap-proach to evaporation of small water clusters and other systems [3,4]. In particular, Lee et al. considered the rate mechanism[5], whereas Ming et al.[6]used the free energy perturbation method to correct for the errors occurring when some standard models for the bulk water were applied to the water dimer, where the evapo-ration rate constants around 1010s1derived from the ab initio
po-tential energy surfaces.
Many chemists observed significant anharmonic effects in dis-sociation of molecular systems, especially of clusters and mole-cules with highly flexible transition states[7–17]. On the other hand, variational transition state theory (VTST) provides a power-ful computational tool for the studies of rates of chemical reac-tions. VTST has been widely used for the modeling of barrierless radical–radical recombination, dissociation processes and applied successfully to a variety of reactions both in the gaseous and con-densed phases[2–4,17,18].
The purpose of the present study is to compute accurate nharmonic rate constants for dissociation (evaporation) of the water dimer at the temperatures where earlier experimental mea-surements[1]were performed. To calculate the rate constants, we use the systematic method developed recently by Yao et al.[19] (below, the YL method) to examine the anharmonic effect on the molecular dissociation. In the calculations of the anharmonic effect in this Letter, the Morse oscillator has been chosen to fit the tial energy surfaces (PES) and to simulate the bonding. The poten-tial energy curve for the water dimer dissociation has also been calculated, and the variational transition state (VTS) has been
0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.01.074
*Corresponding author. Address: Department of Physics, Dalian Maritime University, Dalian 116023, China. Fax: +86 0411 84724335.
E-mail address:[email protected](L. Yao).
Contents lists available atScienceDirect
Chemical Physics Letters
located to obtain the dissociation rate constants. To our knowl-edge, this is the first anharmonic study of the temperature depen-dence of the evaporation rate constant for water clusters. 2. Computational methods
2.1. Ab initio calculations
The geometries of the reactants, products, various intermedi-ates and pseudo transition stintermedi-ates have been optimized by using the MP2 method with the 6-311++G**basis set. Vibrational
har-monic and anharhar-monic frequencies, calculated at the same level, are used for characterization of stationary points, zero-point-en-ergy (ZPE) corrections, and for calculations of reaction rate con-stants using transition state and RRKM theory. All stationary points have been positively identified as local minima or transition states. Single-point energies were recalculated by employing the coupled cluster CCSD(T)/CBS and CCSD(T)/aug-cc-pVTZ methods to achieve higher accuracy and reliability. In order to locate varia-tional transition states, we scanned the PES along the reaction coordinate corresponding to the O–O distance between two water monomers. The GAUSSIAN03 and MOLPRO2006 programs were
uti-lized for all ab initio calculations.
2.2. RRKM and variational transition state calculations
We used RRKM theory for computation of rate constants for unimolecular dissociation of the water dimer. Rate constant k(E) at an internal energy E for a unimolecular reaction A*?A–?P
(here A*and P are the reactant and products, respectively, and A–
denotes the activated complex) can be expressed as:
kðEÞ ¼
r
h W– ðE E– Þq
ðEÞ ; ð1Þwhere
r
is the symmetry factor (here we setr
= 1). h is Planck’s constant,q
(E) represents the density of states of the energized reac-tant molecules, and W–(E E–) is the total number of states for thetransition state (activated complex) A–with a barrier E–. E
repre-sents the total internal energy available to the system. We evaluate
q
(E) and W(E) by performing the inverse Laplace transformation of the partition function Q(b),where b ¼1kT. The steepest-descent
method has been used for carrying out the inverse Laplace transfor-mation up to the second-order approxitransfor-mation[19,20].
For a canonical system, according to the TS theory, the unimo-lecular rate constant is given by
kfðTÞ ¼ kT h QðTÞ– Q ðTÞ e E– kT: ð2Þ
For the case of Morse oscillator we have
Q ðbÞ ¼Y N i¼1 X niðmÞ ni¼0 ebEni ð3Þ
with a similar expression for Q–(b). Here n
i(m) represents the
max-imum value of quantum number ni. To consider the anharmonic
ef-fect on a unimolecular reaction, the anharmonic number and density of states for a system of Morse oscillators take a particularly simple form.
Because no distinct transition state exists on the PES for disso-ciation of water dimer H2O H2O (as for the case of a simple
bond-cleavage process), one can consider different positions for the tran-sition state along the reaction path and calculate rate constants corresponding to each of them. In the microcanonical VTST, the minimum in the microcanonical rate constant is found along the reaction path according to the following equation:
dkðEÞ
dq– ¼ 0 ð4Þ
in which q–is the reaction coordinate, such that a different
transi-tion state is found for each different energy. The individual micro-canonical rate constants are minimized at the point along the reaction path where the total number of states W–(E E–) has a
minimum value. Thus, the reaction bottleneck is located where the minimal W–(E E–) is found, that is, the sums of states for
tran-sition state candidates must be calculated along the reaction path. Each of these calculations requires the optimized structure, sin-gle-point energy, zero-point energy, and vibrational frequencies as functions of the reaction coordinate[21–24].
3. Results and discussion
The evaporation process is defined as the removal of a water molecule from a cluster. The reactant is a water dimer, while the product is two water monomers at infinite separation,
ðH2OÞ2! k
2ðH2OÞ: ð5Þ
The optimized structure of the water dimer at the MP2/6-311++G**
level is illustrated in Fig. 1, where the numbers show optimized bond lengths and bond angles.
Definitions of the energetic quantities begin with the PES or interaction potential that is a function of all internal coordinates of the water cluster. A series of energies at different distances be-tween two dissociating fragments corresponding to the length of a bond to be broken during dissociation is calculated, which is con-sidered to be the reaction coordinate. To obtain these energies, we performed partial geometry optimization with fixed values of the reaction coordinate and all other geometric parameters being opti-mized. No intrinsic barrier exists for dissociation of the water di-mer because the reverse process, condensation of a water molecule onto the cluster, is barrierless. We have used the follow-ing procedure for the VTST calculations. First, we scan the PES con-sidering the dimer separating to two water molecules with varying the bond O–O distances. We calculated 25 optimized geometries corresponding to the O–O distances from 2.9 to 10.0 Å with the step of 0.3 Å. To obtain most accurate position of the variational TS, additionally 40 optimized geometries were calculated at the O–O distances from 7.2 to 7.4 Å. The PES scan along the O–O dis-tance was carried out at the MP2/6-311++G**(Fig. 2a) and B3LYP/
6-311++G**(Fig. 2b) levels of theory. As expected,Fig. 2indicates
that no barrier exists for the reaction. The MP2 relative energies at each point along the reaction path inFig. 2a were scaled by a fac-tor obtained by comparing the water dimer dissociation energies at the CCSD(T)/CBS and MP2/6-311++G**levels and the results were
used to fit the Morse potential energy function. The reaction en-ergy calculated at the CCSD(T)/aug-cc-pVTZ and CCSD(T)/CBS lev-els is 3.19 and 2.33 kcal/mol, respectively, with ZPE, in close agreement with the other theoretical values for the evaporation of water dimer 2.42 ± 0.58 kcal/mol[5]and 3.00 kcal/mol[24]. In
1.948
0.966
0.959
0.961
0.961
103.6
1
2
3
4
5
6
Fig. 1. MP2/6-311++G**-optimized geometries (bond length in Å and angles in degrees) of the reactants water dimer.
this case, the TS, we called it pseudo TS, is located at the O–O dis-tance of 7.365 Å, where the energy barrier is 3.19 kcal/mol. We recomputed the single point energies at the CCSD(T)/CBS and CCSD(T)/aug-cc-pVTZ levels. Finally, the VTS with the minimal val-ues of the number of states were employed to calculate rate con-stants of the direct dissociation processes by using the above RRKM formalism, and the same minimum point was also taken into the calculation of the rate constant for the canonical and microcanonical case.
The geometric and energetic parameters of the reactant and variational TS are collected inTable 1including apparent activation energies, relative energies of the VTS with respect to the reactant. It is critical to compute the activation energy as accurate as possi-ble, because evaporation of the water dimer is an activated process and such activated chemical reactions have rates that depend exponentially on the height of an energy barrier. The anharmonic and harmonic evaporation rate constants calculated using VTST are presented inTable 2, at temperatures from 243 to 1000 K for the canonical system. Corresponding toTable 2, the evaporation rate constants for the water dimer are plotted inFig. 3. In the ear-lier theoretical works, the rate constants were predicted to be around 1011–1012s1by Garrett et al.[2–4]and 1010–1011s1by
Ming et al.[6], which is almost 100–1000 times higher than our
anharmonic rate constants (around 109s1at 243 K). We discuss
canonical evaporation rate constant in which the anharmonic ef-fect is included and compare with the DNT results of Garrett et al.[5].
The relationship between the total energy of a microcanonical system and the temperature of a canonical system is the following:
E ¼ @ln Q @b
: ð6Þ
The total energies corresponding to the canonical temperatures of 243, 273, 303, 333, 500, 600, 700, 800, 900, 1000 K are 1.287, 1.571, 1.867, 2.174, 4.035 (1411 cm1), 5.258 (1839 cm1), 6.551
(2291 cm1), 7.909 (2766 cm1), 9.33 (3263 cm1), 10.812
(3781 cm1) kcal/mol, respectively. The first four energies are all
lower than the calculated activation energy of 3.017 kcal/mol. Hence, we have to calculate the microcanonical rate constant at higher total energies of 1411–4000 cm1, corresponding to the
vibrational frequencies of the water dimer and high canonical tem-peratures above 500 K. Table 3shows the anharmonic and har-monic rate constants of evaporations water dimer at the total energies of 1411–4000 cm1and the data are also illustrated in
Fig. 4. From the results we can see that the rate constant increases with the increasing total energy and the anharmonic result is in the
2 3 4 5 6 7 8 9 10 11 -152.928 -152.926 -152.924 -152.922 -152.920 -152.918 -152.916 2 3 4 5 6 7 8 9 10 11 -152.560 -152.558 -152.556 -152.554 -152.552 -152.550 T
otal Energy (Har
tree)
B3L
YP/6-311++G**
distance between O-O (Angstrom)
T
otal Energy (Har
tree)
MP2/6-311++G**
Fig. 2. Changes of the relative MP2/6-311++G**
(a) and B3LYP/6-311++G**
(b) energies as function of coordinate distance between O1 and O4 in Å.
Table 1
The properties of water dimer and TS, frequencies in cm1
, bond length in Å, and angles in degrees.
E MP2/6-311++G**
B3LYP/6-311++G**
Minimum TS Minimum TS
Zero-point energy (Hartree) 0.047122 0.044165 0.046226 0.043218
Imaginary frequencies – 10.2873i – 9.7920i
O1–O4 distance 2.912 7.365 2.902 7.365
O1–H3 distance 0.966 0.960 0.970 0.962
Angle (H2–O1–H3) 103.6 103.2 105.1 104.9
Total energies (Hartree) 152.559525 152.550386 152.926350 152.917590
Barrier (kcal/mol) 3.19a
3.08b Barrier for MP2/6-311++G**
is 3.10 kcal/mol, and for using CCSD(T)/aug-cc-pvtz level the barrier is 3.07 kcal/mol. a
The scaling factor is 0.88 of MP2/6-311++G**
and CCSD(T)/CBS (energy barrier is 5.27 kcal/mol without zero point energy). b
range of 1011cm1which is closer to the predictions from the
lit-erature [6]. The harmonic rate constants increase sharply from 1012s1at 1411 cm1to 1014s1at 4000 cm1. Moreover, as seen
inTable 3, even though the anharmonic VTST evaporation rate con-stant is nearly insensitive to the available total energy in the given range, the harmonic rate constant is highly sensitive to the total energies. Much difference appears, in the anharmonic case, W = 169 at the total energy 4000 cm1, while W = 2039698 in the
harmonic case.Table 3shows the anharmonic and harmonic total number of states, density of states and rate constants of evapora-tions water dimer at the total energies of 1411–4000 cm1 for
microcanonical case (corresponding to 4.035–11.436 kcal/mol for canonical case). The dissociation rate constant with respect to the microcanonical total energies are also illustrated inFig. 4. From the results inTable 3, much difference appears that total number of states W and density of states
q
in the case is much lower than that of harmonic case. For example, at the total energy 4000 cm1, Wand
q
equal 169 and 20.97 separately for anharmonic case, while2039698 and 139.72 for harmonic case. The big difference is caused by the different model, harmonic and anharmonic poten-tial, which are used to simulate the vibrational bonds. For the dif-ferent models and the difdif-ferent vibrational states, we count the total number of states and density of states, which affect the disso-ciation rate constant. Both of the canonical and microcanonical rate constants increase with the increasing of total energy basically. The harmonic rate constants increase sharply from 2.99 1012s1 at 1411 cm1 to 4.38 1014s1 at 4000 cm1,
while the anharmonic result is in the range of 1.83 1011–
2.42 1011cm1, which are closer to the calculations from the
lit-erature[6]. Our anharmonic results are in reasonable agreement with the experimental predictions. For instance, the calculated va-lue is very close to the result measured by Lee et al.[1]for the hydrogen bond lifetime (reciprocal of the decay rate) of about 109s in water clusters.
Within both harmonic and anharmonic approaches, the canon-ical and microcanoncanon-ical calculations give very close results. The dividing surface VTST method and RRKM theory within the YL method can be effectively used in the studies of water clusters. In the liquid state, there exist many more ways for energy to be ap-plied to a bond, so one would expect the hydrogen bond lifetime to be shorter than that in clusters but still within two or three orders of magnitude. In contrast to the previous work, we have obtained the total number of states and density of states in order to calculate the microcanonical VTST. VTST evaporation comparison is made between the microcnonical VTST evaporation rate constants and canonical rate constants. The results show the behavior of anhar-monic rate constants VST is quite different to that of the haranhar-monic evaporation rate constants. For the flexible system, like water di-mer, the Morse potential is more suitable than the harmonic po-tential to simulate the vibrational bond. The anharmonic effect should be included to obtain the dissociation rate constant of such systems. The method used in this Letter can be extended to treat larger clusters of water or other molecules. The calculations then would become more complex than for the water dimer, not only due to increasing computation costs, but also because there are more contributions to the anharmonic effect from numerous vibra-200 300 400 500 600 700 800 900 1000 1100
10 12 14
Temperature (K)
rate constant log(k s
-1 ) Anharmonic
Harmonic
Fig. 3. Canonical rate constant versus temperature for the evaporation of water dimer.
Table 2
Rate constant of evaporation of water dimer for different temperatures for canonical system. Temperature
(K)
243.0 273.0 303.0 333.0 500.0 600.0 700.0 800.0 900.0 1000.0
Correspond energy (kcal/mol) 1.287 1.571 1.867 2.174 4.035 5.258 6.551 7.909 9.33 10.812
Correspond (cm1) 450 549 653 760 1411 1839 2291 2766 3263 3781
Anharmonic rate constant of canonical case (1/s) 4.75 109 8.42 109 1.31 1010 1.84 1010 4.65 1010 6.44 1010 8.06 1010 9.52 1010 1.09 1011 1.21 1011
Harmonic rate constant of canonical case (1/s) 2.84 1011 7.66 1011 1.72 1012 3.36 1012 3.42 1013 7.58 1013 1.35 1014 2.08 1014 2.93 1014 3.85 1014 Table 3
Rate constant of evaporation of water dimer for different total energies for microcanonical system. Total energy
(cm1 )
1411 1500 1839 2000 2291 2766 3000 3263 3781 4000
Correspond total energy (kcal/mol)
4.035 4.289 5.258 5.718 6.551 7.909 8.577 9.330 10.812 11.436
Anharmonic W(E) of TS 30 33 35 44 46 65 88 94 134 169
Anharmonicq(E) of reactant (1/cm1
)
4.68 4.79 5.71 6.74 7.20 9.53 12.20 12.77 17.10 20.97
Anharmonic rate constant (1/s) 1.95 1011 2.08 1011 1.83 1011 1.97 1011 1.94 1011 2.07 1011 2.16 1011 2.21 1011 2.35 1011 2.42 1011 Harmonic W(E) of TS 229 517 4944 10918 35282 152991 276825 503386 1393707 2039698 Harmonicq(E) of reactant
(1/cm1)
2.29 2.81 5.69 7.73 12.96 27.58 38.75 55.68 107.74 139.72
Harmonic rate constant (1/s) 2.99 1012 5.53 1012 2.60 1013 4.23 1013 8.16 1013 1.66 1014 2.14 1014 2.71 1014 3.88 1014 4.38 1014
tional frequencies. Another complication for bigger clusters is the possibility of existence of a large variety of different isomers and conformations. Nevertheless, the present study clearly demon-strates that the role of anharmonic effects in the kinetics of cluster decomposition is very important and should be taken into account. Acknowledgments
This work was supported by the National Natural Science Foun-dation of China (Grant No. 10704012), NSC (Taiwan), Academia Sinica and the Natural Science Foundation of Liaoning Province (Grant Nos. 2006Z064, 20082140).
References
[1] M.F. Vernon, D.J. Krajnovich, H.S. Kwok, J.M. Lisy, Y.R. Shen, Y.T. Lee, J. Chem. Phys. 77 (1982) 47.
[2] Shawn M. Kathmann, Gregory K. Schenter, Bruce C. Garrett, J. Chem. Phys. 111 (1999) 4688.
[3] Gregory K. Schenter, Shawn M. Kathmann, Bruce C. Garrett, Phys. Rev. Lett. 82 (1999) 3484.
[4] Gregory K. Schenter, Shawn M. Kathmann, Bruce C. Garrett, J. Chem. Phys. 116 (2002) 4275.
[5] Mal-Soon Lee, F. Baletto, D.G. Kanhere, S. Scandolo, J. Chem. Phys. 128 (2008) 214506;
S.M. Kathmann, B.J. Palmer, G.K. Schenter, B.C. Garrett, J. Chem. Phys. 128 (2008) 064306.
[6] Yi. Ming, Geeling Lai, Chinghang Tong, Robert H. Wood, Douglas J. Doren, J. Chem. Phys. 121 (2004) 773.
[7] R. Krems, S. Nordholm, Z. Phys. Chem. 214 (2000) 1467;
D. Shen, H.O. Pritchard, J. Chem. Soc. Faraday Trans. 92 (1996) 1297; P. Hobza, Z. Havlas, Chem. Rev. 100 (2000) 4253.
[8] J. C Tou, .S.H. Lin, J. Chem. Phys. 49 (1968) 4187; S.H. Lin, H. Eyring, J. Chem. Phys. 39 (1963) 1577; S.H. Lin, H. Eyring, J. Chem. Phys. 43 (1964) 2153.
[9] (a) S.A.C. McDowell, J. Mol. Struct.: Theochem. 770 (2006) 119; (b) L.B. Bhuiyan, W.L. Hase, J. Chem. Phys. 78 (1983) 5052. [10] G.H. Peslherbe, W.L. Hase, J. Chem. Phys. 105 (1996) 7432. [11] (a) S.E. Stein, B.S. Rabinovitch, J. Chem. Phys. 58 (1973) 2438;
(b) T. Beyer, D.F. Swinehart, Commun. Assoc. Comput. Mach. 16 (1973) 379; (c) Ian M. Mills, Quantum Chemistry, Theo. Chem., vol. 1, The Chemical Society, London, 1974. p. 110;
E.W. Schlag, R.A. Sandsmark, J. Chem. Phys. 37 (1962) 68. [12] W.L. Hase, Acc. Chem. Res. 31 (1998) 659;
Kihyung Song, W. Hase, J. Chem. Phys. 110 (1999) 6198.
[13] V.N. Bagratashvilli, V.S. Letokhov, A.A. Makarov, E.A. Ryabov, Laser Chem. 1 (1983) 211.
[14] S.S. Mitra, S.S. Bhattacharyya, J. Phys. B: At. Mol. Opt. Phys. 27 (1994) 1773. [15] J. Troe, Chem. Phys. 190 (1995) 381;
J. Troe, J. Phys. Chem. 83 (1979) 14; J. Troe, J. Chem. Phys. 79 (1983) 6017;
D. Romanini, K.K. Lehmann, J. Chem. Phys. 98 (1993) 6437. [16] M.R. Hoare, Th.W. Ruijgrok, J. Chem. Phys. 52 (1970) 113. [17] B.C. Garrett, D.G. Truhlar, J. Am. Chem. Soc. 101 (1979) 4534. [18] B.C. Garrett, D.G. Truhlar, J. Am. Chem. Soc. 101 (1979) 5207.
[19] Li. Yao, A.M. Mebel, H.F. Lu, H.J. Neusser, S.H. Lin, J. Phys. Chem. A 111 (2007) 6722;
Li. Yao, S.H. Lin, Mod. Phys. Lett. B 22 (2008) 3043; Li. Yao, S.H. Lin, Sci. China Ser. B 51 (2008) 1146.
[20] H.C. Chang, J.C. Jiang, C.W. Chuang, J.S. Lin, W.W. Lai, Y.C. Yang, S.H. Lin, Chem. Phys. Lett. 410 (2005) 42.
[21] H.C. Chang, J.C. Jiang, W.C. Tsai, G.C. Chen, C.Y. Chang, S.H. Lin, Chem. Phys. Lett. 432 (2006) 100.
[22] Hai-Chou Chang, Jyh-Chiang Jiang, Wei-Cheng Tsai, Guan-Ciao Chen, J. Phys. Chem. B 110 (2006) 3302.
[23] Chiachin Tsoo, Charles L. Brooks, J. Chem. Phys. 101 (1994) 6405.
[24] John C. Owicki, Lester L. Shipman, Harold A. Scheraga, J. Phys. Chem. 79 (1975) 1794. 1000 1500 2000 2500 3000 3500 4000 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0 Total energy (cm-1) rat e const a nt log(k s -1 ) Anharmonic Harmonic
Fig. 4. Microcanonical rate constant versus total energies for the evaporation of water dimer.
