Transition path sampling: rearrangement of cage water hexamer

Chemical Physics 299 (2004) 123–129 www.elsevier.com/locate/chemphys

Transition path sampling: rearrangement of cage water hexamer Jin Yong Lee

*

Department of Chemistry, University of California at Berkeley, Berkeley, CA 94720, USA Department of Chemistry, Chonnam National University, 300 Yongbong-Dong, Bukgu, Gwangju 500-757, Republic of Korea Received 2 October 2003; accepted 10 December 2003

Abstract The hydrogen bond exchange reaction of cage water hexamer has been investigated. With an aid of transition path sampling method, a transition state for this reaction has been found. The transition state ensemble has been investigated as a function of energy, and the distribution of the structures of sampled transition states is consistent with that obtained by other method. From the ab initio calculations, it has been confirmed to be a transition state with one imaginary frequency along the reaction coordinate. The most important structural parameters for the transition are retardation (increase of ra and O


O distance) and C2 rotation of a water molecule. I obtained rate constants at temperatures ranging from 0 to 300 K considering quantum mechanical tunneling effect. At temperatures below 100 K, the tunneling effect is very important, while at high temperatures, the classical transition state theory approach gives reasonable data. Transition path sampling method is very useful to find transition states for many complicated chemical reactions in which we do not know about transition states or it is hard to find a transition state by quantum mechanical approaches. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Transition path sampling; Hydrogen exchange; Rate constant; Tunneling; Water hexamer

1. Introduction In the last decades there have been extensive studies on the water clusters ranging from dimer to over 108mer because of the fundamental role of water as a solvent. Water has a number of hydrogen bonds, thus understanding the bonding in water clusters is of great importance [1–3]. In particular, water hexamer has many isoenergetic structures with differences in the binding energies less than 0.5 kcal/mol, and different relative stabilities have been predicted at different levels of theory [4–7]. From high level calculations, Clary et al. [6] predicted the cage hexamer as the minimum energy form and they compared the experimental rotational constants by vibration–rotation–tunnelling (VRT) with those calculated for low-lying structures. Then, the VRT spectra agreed most closely with those predicted for cage

*

Tel.: +82-62-530-3384; fax: +82-62-530-3389. E-mail address: [email protected] (J.Y. Lee).

0301-0104/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2003.12.020

hexamer. On the other hand, modified IR predissociation experiments gave an observation of noncyclic water hexamer [8]. Recently, active studies have been focused on the rearrangement of water clusters both theoretically [9,10] and experimentally [11–16] ranging from dimer to hexamer. It is said that there are three different types of rearrangement, flipping, bifurcation, and the C2 tunneling in the order of rare event. For the rearrangement of water clusters, it is most desirable to answer the questions such as ‘‘What is the transition state?’’, ‘‘What is the reaction mechanism?’’, ‘‘What is the rate constants?’’, etc. To answer these questions, pathway calculations are highly desirable because it may be unreasonable to study the rearrangement from the transition state alone. Chandler and coworkers [17] developed a method to sample transition paths. For both stochastic [18] and deterministic [19] pathways, this method has been applied to study the dynamic processes and rate constants in 2-dim Lennard–Johns Disc rearrangements. This transition path sampling method has also been applied

124

J.Y. Lee / Chemical Physics 299 (2004) 123–129

to hydrogen bond breaking in water [20], and many other chemical processes [21–23]. Based on this transition path sampling method, I study the rearrangement of cage water hexamer, and try to obtain information on the transition state. To obtain the reasonable transition pathways for the reaction from A to B, all the possible pathways should be included with exact path probabilities. This may be realized by Monte Carlo sampling. Strictly speaking, the hydrogen bond rearrangement is not the rare event at room temperature. However, at low temperature like 100 K, this rearrangement can be considered as a rare event. One can often obtain some insight into mechanisms by determining the transition states [24–27]. However, it is problematic to search all the transition states for complex systems of current interest in chemical and biological processes, e.g., solvation dynamics [25], conformational transition [28], or transport processes in condensed matter systems [29] because there are too many saddle points in the potential energy surfaces. Usually, the transition states can be specified by searching for the saddle points on the potential energy hypersurfaces [30–33] in small systems. An alternative is to define transition state statistically by considering probabilities to relax into either one of stable states. Here, I study the transition state ensemble for the rearrangement of water hexamer. With this information, we find a transition state by quantum mechanical approach, and obtain the rate constant. The paper is organized as follows: Section 2 describes the details on the simulation, and Section 3 explains transition state ensemble. Section 4 presents results on the hydrogen bond rearrangement of cage hexamer. Section 1 concludes in Section 5.

tained by the integration of Hamiltonian equations of motion. The statistics of the system can be described by an ensemble and the probability to find the system in a certain state x is given by the distribution qðxÞ. For example, the microcanonical (system is mechanically and thermally isolated with energy E) distribution function is qmc ðxÞ ¼ dðH ðxÞ  EÞ=Qmc ;

ð1Þ

where H is the hamiltonian and Qmc is the partition function which normalizes the distribution. Here, the system is limited to the microcanonical ensemble and the Hamiltonian equations of motion. I restrict the ensemble qðx0 Þ of all possible initial conditions to an ensemble fAB ðx0 Þ of those initial conditions located in region A, which are starting points of trajectories reaching region B within a maximum time s fAB ðx0 Þ  qðx0 ÞhA ðx0 ÞHB ðxs Þ=QAB ;

ð2Þ

where QAB is the normalizing factor Z QAB ¼ dx0 qðx0 ÞhA ðx0 ÞHB ðxs Þ:

ð3Þ

The function HB ðxs Þ is equal to 1 if the system visits B at least once in time s and vanishes otherwise:
1 if t exists in ½0; s satisfying hB ðxt Þ ¼ 1; HB ðxs Þ ¼ 0 otherwise; ð4Þ where hA;B ðxt Þ ¼




1 0

xt 2 A; B; xt 62 A; B:

ð5Þ

Each phase point for which fAB 6¼ 0 indicates a transition path connecting region A with region B in time s. 2.2. Sampling the path ensemble

2. Simulation details For the sampling of each path, I carry out constraint MD simulations with the time step of 1 fs and 800 time steps. The configurations and momenta out of every 10 time steps (10 fs) have been stored for the data analysis. In other words, I use 80 beads (which can be denoted as t0 ; t1 ; t2 ; . . . ; s) for each pathway and the time duration between adjacent beads is 10 fs. For the constraint of rigid water molecules, the rattle algorithm is used. I use TIP4P [34] potential which is often used in the study of water clusters. The system is restricted to be nontranslating and non-rotating, that is, the total linear and angular momentum are zero. 2.1. Deterministic transition pathways In classical systems specified by the phase coordinates x ¼ ðq; pÞ, given the initial state x0 at time t ¼ 0, the trajectories evolving in time tðxt ¼ xt ðx0 ÞÞ can be ob-

An ensemble of transition paths can be obtained by sampling with a Monte Carlo procedure based on the Eq. (2). In the Monte Carlo procedure, a new path xn0 is generated from an old path starting xo0 from an initial condition. Then this new path is accepted or rejected based on the detailed balance condition o!n o!n n!o n!o fAB ðX0o ÞPgen Pacc ¼ fAB ðX0n ÞPgen Pgen ;

ð6Þ

where Pgen and Pacc are the probabilities to generate and accept the path, respectively. The superscripts, o and n, denote the old and new path, respectively. This acceptance probability makes sure that the paths are sampled with the correct weight of fAB . Here, I generate a new path by the shooting method which has been successfully applied in samplings both with stochastic and deterministic dynamics. In the shooting algorithm, first, a phase configuration is randomly chosen from the phase configurations of the existing path (old path). And then the momentum is

J.Y. Lee / Chemical Physics 299 (2004) 123–129

changed by a small amount of dp by Gaussian distribution. Starting with the given space coordinates and modified momenta, a new path is generated by integration of Newtonian equations of motion backward and forward in time until the time t ¼ 0 and t ¼ s have reached, respectively. The new path is then accepted or rejected according to the Metropolis criterion [35]. The detailed description on the shooting algorithm is available in [17–19]. 2.3. Defining the stable states The reactant and product regions are defined by the hydrogen bond lengths which are directly related to the hydrogen bond rearrangement reaction (or hydrogen bond exchange). For the cage water hexamer, there are two isomers, X and Y. As shown in Fig. 1, the system is considered to be region A (reactant) if the hydrogen Ha is hydrogen bonded. On the other hand, it is considered to be region B (product) if the Hb is hydrogen bonded. Specifically, the system is said to be in region A if the hydrogen bond length ra is
and in region B if rb is less than 2.30 less than 2.30 A

is chosen A. The critical hydrogen bond length 2.3 A because most hydrogen bonds in water hexamer are
which has also been confirmed from the within 2.3 A, equilibrium molecular dynamics simulation of cage water hexamer using TIP4P potential. These two conditions are exclusive each other, but both conditions are used because those two sets, A and B, do not make a complete set. For example, a certain configuration may not belong to any of the two sets. The characteristic functions, hA ðxÞ and hB ðxÞ (x is a function of time), can then be written as
1 ra ; rb 6 2:30; ð7Þ hA;B ðxÞ ¼ 0 otherwise:

2.4. Rate constant The fluctuation–dissipation theorem relates the rate constant for a transition between stable states to microscopic correlation function [36]

D kðtÞ 

hA ðxð0ÞÞh_ B ðxðtÞÞ hhA ðxð0ÞÞi

125

E  kA!B expðt=trxn Þ;

ð8Þ

where hA and hB are the characteristic functions for the stable states. Whereafter, for convention, xð0Þ is neglected otherwise specified. This can be factorized as kðtÞ 

hhA h_ B ðxt Þi hhA hB ðxs Þi   mðtÞ  P ðsÞ: hhA hB ðxs Þi hhA i

ð9Þ

In the path sampling method, the rate constant can be obtained by the product of frequency factor, mðtÞ, and probability factor, P ðsÞ, which is the probability to find the system in product region at the end of time, s, starting from the reactant region. To calculate the frequency factor mðtÞ in the path ensemble, it is further factorized as mðtÞ ¼

hhA h_ B ðxt ÞH i hhA HB i hh_ B ðxt ÞiAB  ¼ ; hhA HB i hhA hB ðxt ÞHB i hhB ðxt ÞiAB

ð10Þ

where HB ¼ 0 only if hB ðtÞ ¼ 0 for all time t for a path. The notation hQðxt ÞiAB  hhA Qðxt ÞHB i=hhA HB i denotes the average of the quantity Q at time t over the sampled paths. Since mðtÞ contains the whole dependence of kðtÞ on time t, it has to reach a plateau in the time range tmol < t  trxn . This fact can be used to check whether the total simulation time s is long enough to sample all the relevant transition paths. For the calculation of the probability factor, the product region are divided into six windows (B1 –B6 ) as follows: 8 B1 : 1:50 6 rb 6 2:05; > > > > B > 2 : 2:05 6 rb 6 2:45; > < B3 : 2:45 6 rb 6 2:70; ð11Þ B4 : 2:70 6 rb 6 2:90; > > > > B : 2:90 6 rb 6 3:10; > > : 5 B6 : 3:10 6 rb 6 4:00: These windows (Bi s) are further divided into a number of thin shells and histogram the probability to find the end point in these shells. In addition, the adjacent windows should have overlap region so that the probability can be matched as a continuous function. After that, matching the histograms of each window results in the entire probability distribution to find the end point in

Fig. 1. Equilibrium and transition state of the rearrangement of cage water hexamer.

J.Y. Lee / Chemical Physics 299 (2004) 123–129

the shell characterized by the order parameter rb . The desired probability factor P ðsÞ is then easily found by R rb P ðrb Þdrb P ðsÞ ¼ R 01 ; ð12Þ P ðrb Þdrb 0 where rb denotes the constraint value to define the region B in the calculation of frequency factor (e.g., 2.30
in this study). Finally, the rate constant can be obA tained by the product of frequency and the probability factors.

3. Transition state ensemble By the transition path ensemble, the properties of transition states which are not known can be investigated. Traditionally, the transition state is defined by the saddle point of the potential energy surface. More generally, it can be described by an ensemble of transition states with a distribution fTS ðxÞ. A configuration r can be defined as a transition state if Z dpHA ðxs ðr; pÞÞd½ðr; pÞ  E Z ¼ dpHB ðxs ðr; pÞÞd½ðr; pÞ  E; ð13Þ where the path length s is of the order of the transient time tmol . That is, a configuration r is a transition state if the probability to reach A equals to that to reach B starting from r with randomly distributed momentum. Using this criteria, it is possible to find at least one transition state for each transition pathway which connects A with B. To find the transition states, I calculate the probabilities, pA and pB , to reach A and B, respectively, at each configuration along the transition path. Practically, each configurations are stored for every time intervals. From each configuration, one shoots off a number of test trajectories with momenta reinitialized from a random isotropic distribution, and count how many shots reach AðnA Þ and BðnB Þ. Finally, the configuration is considered as a transition state if the nA is equal to nB at the configuration.

which satisfies the reactant and product region criteria after tens of trials. It has recently shown that how the path sampling method can be used for the calculation of rate constants and the study of transition state. Fig. 2 shows the characteristic functions (hhB xt ÞiAB ) as a function of time at E ¼ 41:0 kcal/mol. The linear region of the hhB ðxt ÞiAB ensures that the simulation time is enough to study the reaction of this system. From the first derivative of hhB ðxt ÞiAB , the frequency factor is obtained to be about 1.50  1012 s1 . The probability factor is obtained from the distribution of pðrb Þ as in Eq. (12). To obtain the probability factor, a million pathways have been sampled in the umbrella sampling calculations. At E ¼ 41:0 kcal/mol, this probability factor is 8.77  103 . Thus, the rate constant is 1.32  1010 s1 . Fig. 3 shows the free energy (which is defined by  lnðpðrb ÞÞ) as a function of rb at E ¼ 41:0 kcal/mol. The equilibrium energy at TIP4P potential is )47.2 kcal/ mol. Considering the TIP4P potential energy barrier of 2.88 kcal/mol, the total energy of )41.0 kcal/mol corresponds to the kinetic energy of 3.3 kcal/mol, which

1.0 0.8

AB

126

0.6 0.4 0.2 0.0 02

00

400

600

800

t(fs) Fig. 2. Characteristic function at E ¼ 41:0 kcal/mol.

20

4. Results and discussion To obtain an initial transition path, I assume that the hydrogen atoms may exchange simply by rotation along the C2 axis of the water monomer which involves the hydrogen exchange. First, I choose the configuration near the maximum potential during the C2 rotation. Then, starting from this configuration, the forward and backward movements are carried out by the integration of Newtonian equations of motion using randomly modified momenta until I obtain a good initial pathway

βF

15 10 5 0 1

2

3

4

5

rb Fig. 3. Free energy (bF ) distribution as a function of order parameter (rb ) at E ¼ 41:0 kcal/mol total energies.

J.Y. Lee / Chemical Physics 299 (2004) 123–129

further corresponds to temperature of about 1000 K roughly estimated from the fact that the average kinetic energy is equal to 3RT =2. For comparison, I obtained rate constants as a function of temperature with inclusion of tunneling transmission factor, and represented in Table 1 and Fig. 4. At very low temperatures, the rearrangement of water hexamer occurs very rarely and could not be observed from our simulations. Thus, the rate constants obtained by transition path sampling were not listed in Table 1. For the calculation of tunneling rate constant, I used the DOIT program [37], which have been applied successfully to many chemical reactions [38–41]. For the application of DOIT program, I obtained stable and transition states of the hydrogen exchange reaction of cage water hexamer with density functional theory (DFT) employing BeckeÕs three parameter hybrid method using the Lee–Yang–Parr correlation functional (B3LYP) [42] with 6-31G* basis sets by means of the Gaussian 98 suite of programs [43]. At temperatures

Table 1 Rate constants for the hydrogen exchange of cage water hexamer Temperature

kTST 2.90  10 1.31  103 2.69  106 1.34  108 1.46  109 7.34  109

50 100 150 200 250 300

k

j 7

11

4.61  104 1.18  105 6.41  106 2.09  108 1.97  109 9.27  109

1.59  10 9.00  101 2.38 1.56 1.36 1.26

kTST , j, and k represent the rate constants obtained classical transition state theory, tunneling transmission factors, and the total rate constants, respectively. All units are in s1 .

12

-1

rate (sec )

8 4 0

y = -0.9822x + 13.048

-4 -8 0

5

10

15

20

1000/T(K) Fig. 4. Rate constants for the rearrangement of cage water hexamer. Solid line with cross represents the results obtained by classical transition state theory, and the dotted line with open circle by inclusion of quantum tunneling transmission factor.

127

below 100 K, the tunneling effect comes to be important, and around 50 K, the reaction is only controlled by tunneling. The reaction barriers for B3LYP/6-31G* and TIP4P potential energy surfaces are 5.03 and 2.88 kcal/ mol, respectively. The B3LYP/6-31G* method has been often used to study the mechanism of chemical reactions, that is, energetics of reactants, products, and transition states [44]. Since the purpose of this paper is to test the new path sampling method to search the transition state for the rearrangement with small barriers, the difference of those potential energy barriers was not further manipulated by other calculational methods. The rate constants obtained by path sampling with TIP4P should be compared with the classically measured limiting value corresponding to temperature of about 1000 K in Fig. 4 as discussed earlier. At temperature of 200 K which corresponds to the average kinetic energy of about 0.6 kcal/mol, an initial transition path connecting the reactant with the product could not be obtained for the rearrangement of water hexamer. For more accurate rate constants, the path sampling method should be implemented with first principle Hamiltonians. As a matter of fact, recently, path sampling method has been successfully implemented to the Car–Parrinello molecular dynamics (CPMD) based on ab initio potential energy surface [45]. The structures of a thousand transition states determined in this sampling have been investigated at different total energies, and the distribution of important structural parameters (which show prominent changes at transition state from the equilibrium) were shown in Fig. 5 along with equilibrium state distribution. This TS distribution is very useful to find an accurate transition state quantum mechanically. Indeed, I had a difficulty to find a transition state by ab initio method because the cage structure is often destroyed and changed to other water hexamer conformers. However, with an aid of TS information obtained by path sampling, I could obtain a transition state. The transition states distributions are sharper than equilibrium state distribution, which implies that the transition happens abruptly. The ra and rb are highly correlated as seen in Figs. 5(a) and (b). These two parameters change in a very symmetrical way. The most probable values of ra , rb , and r(O2


O6 ) of the
respectransition structures are 2.65, 2.65, and 3.05 A, tively. And the average values are 2.595, 2.711, and
respectively. These values are in good agree3.078 A,
of ra for the transition state ment with the value (2.64 A) obtained by Wales from the eigenvector following method [9]. These values can be compared with ab initio
respecresults, which are 2.458, 2.476, and 2.835 A, tively. From these structural parameters of the transition states, the reaction mechanism can be understood. When the hydrogen exchange reaction occurs, five water molecules can be considered to be fixed. At first, the water molecule labeled O6 backs off, and then it rotates

128

J.Y. Lee / Chemical Physics 299 (2004) 123–129 1.0

4.0 E=-41. 0

probability

0.8

Eq uil

3.0

0.6

0.4

2.0 0.2 0.0 1.5

2.0

(a)

2.5

3.0

3.5

ra

probability

Eq uil

0.6

0.2

2.5

(b)

3.0

3.5

4.0

rb 1.0 E=-4 1.0

probability

0.8

Equil

0.6 0.4 0.2

(c)

400

600

800

5. Conclusion

0.4

0.0 2.0

200

Fig. 6. Changes of average hydrogen bond lengths (solid line, ra ; dotted line, rb ; opened circles, r(O2


O6 )) along the trajectories at total energy of )41.0 kcal/mol.

E=-4 1.0

0.0 2.0

0

t (fs)

1.0 0.8

1.0

2.5

3.0 r (O2...O3)

3.5

4.0

Fig. 5. Structural parameters representing the transition states distributions at E ¼ 41:0 kcal/mol: (a) ra , (b) rb , and (c) r(O2


O6 ).

about the C2 axis. This implies that it is appropriate to call this hydrogen exchange reaction as C2 tunneling named by Saykally. The changes of ra , rb , and r(O2


O6 ) as reaction time are shown in Fig. 6 for one representative trajectory at E ¼ 41:0 kcal/mol. Transition time is the time duration in which the system does not belong to reactant either product when the transition occurs. The transition time for the hydrogen bond exchange reaction of cage water hexamer is about 100–200 fs as noted in Fig. 6.

I have investigated the rearrangement reaction of cage water hexamer. Using the transition path sampling method, I obtained the transition state ensemble for the hydrogen exchange reaction of cage water hexamer. With an aid of path sampling, I could find a transition state by ab initio calculations. The rate constants were obtained as a function of temperature considering tunneling transmission factor. At low temperatures below 80 K, the quantum mechanical tunneling factor is extremely important for the hydrogen exchange reaction of cage water hexamer. This transition path sampling method can be useful for the study of reaction kinetics for most chemical reactions ranging from unimolecular reaction to the protein folding and enzyme–substrate reaction dynamics.

Acknowledgements This work was supported by Department of Energy Grant (DE-FG03-99ER14987). I also wish to acknowledge the financial support of Korea Research Foundation made in the program year of 1997-0035. I am grateful to Prof. David Chandler and Prof. Christoph Dellago for the initiation of this work and detailed discussions. I also thank Prof. Kwang S. Kim for helpful discussions.

References [1] D.D. Nelson, G.T. Fraser, W. Klemperer, Science 238 (1987) 1670. [2] R.C. Cohen, R.J. Saykally, Annu. Rev. Phys. Chem. 42 (1991) 369. [3] J.M. Hutson, Annu. Rev. Phys. Chem. 41 (1990) 123.

J.Y. Lee / Chemical Physics 299 (2004) 123–129 [4] K.S. Kim, M. Dupuis, G.C. Lie, E. Clementi, Chem. Phys. Lett. 131 (1986) 451. [5] K. Kim, K.D. Jordan, T.S. Zwier, J. Am. Chem. Soc. 116 (1994) 11568. [6] K. Liu, M.G. Brown, C. Carter, R.J. Saykally, J.K. Gregory, D.C. Clary, Nature 381 (1996) 501. [7] J. Kim, K.S. Kim, J. Chem. Phys. 109 (1998) 5886. [8] R.N. Pribble, T.S. Zwier, Science 265 (1994) 75. [9] D.J. Wales, in: J. Bowman, Z. Bacic (Eds.), Advances in Molecular Vibrations and Collision Dynamics, vol. 3, JAI press, 1998, pp. 365–396. [10] T.R. Walsh, D.J. Wales, J. Chem. Soc. Faraday Trans. 92 (1996) 2505. [11] K. Liu, M.G. Brown, R.J. Saykally, J. Phys. Chem. 101 (1997) 8995. [12] N. Pugliano, R.J. Saykally, J. Chem. Phys. 96 (1992) 1832. [13] N. Pugliano, R.J. Saykally, Science 257 (1992) 1937. [14] M.R. Viant, J.D. Cruzan, D.D. Lucas, M.G. Brown, K. Liu, R.J. Saykally, J. Phys. Chem. A 101 (1997) 9032. [15] J.D. Cruzan, L.B. Braly, K. Liu, M.G. Brown, J.G. Loeser, R.J. Saykally, Science 271 (1996) 59. [16] K. Liu, M.G. Brown, J.D. Cruzan, R.J. Saykally, Science 271 (1996) 62. [17] C. Dellago, P.G. Bolhuis, F.S. Csajka, D. Chandler, J. Chem. Phys. 108 (1998) 1964. [18] C. Dellago, P.G. Bolhuis, D. Chandler, J. Chem. Phys. 108 (1998) 9236. [19] P.G. Bolhuis, C. Dellago, D. Chandler, Fraday Discuss. 110 (1998) 421. [20] P.L. Geissler, C. Dellago, D. Chandler, Phys. Chem. Chem. Phys. 1 (1999) 1317. [21] P.L. Geissler, C. Dellago, D. Chandler, J. Phys. Chem. B 103 (1999) 3706. [22] P.G. Bolhuis, C. Dellago, D. Chandler, Proc. Natl. Acad. Sci. USA 97 (2000) 5877. [23] (a) D. Laria, J. Rodriguez, C. Dellago, D. Chandler, J. Phys. Chem. A 105 (2001) 2646; (b) C. Dellago, P.G. Bolhuis, P.L. Geissler, Adv. Chem. Phys. 123 (2002) 1; (c) P.G. Bolhuis, D. Chandler, C. Dellago, P.L. Geissler, Rev. Phys. Chem. 53 (2002) 291; (d) T.J.H. Vlugt, C. Dellago, B. Smit, J. Chem. Phys. 113 (2000) 8791; (e) P.G. Bolhuis, C. Dellago, P.L. Geissler, D. Chandler, J. Phys. Condens. Matter 12 (2000) A147; (f) C. Dellago, P.G. Bolhuis, D. Chandler, J. Chem. Phys. 110 (1999) 6617. [24] P. H€ anggi, P. Talkner, M. Borkovec, Rev. Mod. Phys. 62 (1990) 251. [25] D. Chandler, J. Chem. Phys. 68 (1978) 2959. [26] J.C. Keck, Discuss. Faraday Soc. 33 (1978) 173. [27] J.B. Anderson, J. Chem. Phys. 58 (1973) 4684.

129

[28] J.L. Cleland (Ed.), Protein Folding: In Vivo and In Vitro, American Chemical Society, Washington, DC, 1993. [29] S. Torquato, Appl. Mech. Rev. 44 (1991) 37. [30] J. McIver, A. Komornicki, J. Am. Chem. Soc. 94 (1972) 2625. [31] C.J. Cerjan, W.H. Miller, J. Chem. Phys. 75 (1980) 7800. [32] F. Jensen, J. Chem. Phys. 102 (1995) 6706. [33] K. M€ uller, Angew. Chem. Int. Ed. Engl. 19 (1980) 1. [34] W.L. Jorgensen, J. Chandraesekhar, J.W. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [35] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [36] D. Chandler, Introduction to Modern Statistical Mechanics, Oxford University Press, Inc, 1987. [37] Z. Smedarchina, A. Fernandez-Ramos, W. Siebrand, M.Z. Zgierski, DOIT 1.2, Software for Tunneling Reaction Dynamics, Steacie Institute for Molecular Sciences, National Research Council of Canada. Available from http://steacie.nrc-cnrc.gc.ca/ research_programs/theory_comp/doit_e.html. [38] A. Fernandez-Ramos, Z. Smedarchina, W. Siebrand, M.Z. Zgierski, J. Chem. Phys. 113 (2000) 9714; M. Fernandez-Ramos, Z. Zgierski, J. Phys. Chem. A 106 (2002) 10578. [39] A. Fernandez-Ramos, Z. Smedarchina, W. Siebrand, M.Z. Zgierski, M.A. Rios, J. Am. Chem. Soc. 121 (1999) 6280; A. Fernandez-Ramos, Z. Smedarchina, M.Z. Zgierski, J. Chem. Phys. 113 (2000) 2662. [40] Z. Smedarchina, M.Z. Zgierski, W. Siebrand, P. Kozlowski, J. Chem. Phys. 109 (1998) 1014. [41] A. Fernandez-Ramos, W. Siebrand, A. Fernandez-Ramos, L. Gorb, M.Z. Zgierski, J. Chem. Phys. 112 (2000) 566. [42] (a) C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785; (b) B. Miehlich, A. Savin, H. Stoll, H. Preuss, Chem. Phys. Lett. 157 (1989) 200; (c) A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [43] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery, R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, J. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P.M.W. Gill, B.G. Johnson, W. Chen, M.W. Wong, J.L. Andres, M. Head-Gordon, E.S. Replogle, and J.A. Pople, Gaussian 98 (Revision A), Gaussian, Inc., Pittsburgh, PA, 1999. [44] C.A. Tsipis, P.A. Karipidis, J. Am. Chem. Soc. 125 (2003) 2037. [45] P. Geissler, C. Dellago, D. Chandler, J. Huttler, M. Parrinello, Science 291 (2001) 2121.