J . Phys. Chem. 1988, 92, 1075-1086 nonperfect crystal disorder and are restricted to an energy range of about 1 cm-I, are responsible for the long exponential decay observed in CARS experiment of ref 4. On the other hand, the faster decaying slope is characteristic of the destructive interference between all the individual components of the overall OPDS of the disordered crystal.
1075
Acknowledgment. We warmly acknowledge C. Caray for his help in crystal growth and R. Tierce for her technical assistance. We also thank Prof. R. M. Pick for stimulating discussions on this subject. Registry No. N 2 0 , 10024-97-2.
Dynamics of Reactlons Involving Diffuslve Multidimensional Barrier Crossingt Sangyoub Lee and Martin Karplus* Department of Chemistry, Harvard University, Cambridge, Massachusetts 02138 (Received: July 20, 1987)
An extension of the Kramers’ rate theory is presented for the case involving diffusive multidimensional barrier crossing. The theory is valid for a moderate to high reaction barrier and provides a clear description of the nature of the reactive mode. As an application of the theory, various aspects of the torsional dynamics of n-butane are investigated: they include the coupling between internal torsional motion and the overall translation and rotation of the molecule, the effect of hydrodynamic interaction among extended carbon atoms, and the effect of bond angle bending and bond stretching on the trans-gauche isomerization rate. The results are compared with stochastic dynamic simulation results.
I. Introduction Starting with the celebrated work by Kramers,’ the dynamics of activated rate processes in solution has been a subject of active investigations. Most studies use a model of a Brownian particle moving in a phase space under the effect of thermal noise and damping from the surrounding medium. The phase space can be divided into two regions, a reactant state and a product state in the language of chemical reactions, separated by a potential energy barrier. Either or both of the regions may be a bound state in a potential well or a continuum state. The objective is to calculate the rate of transition from the reactant state to the product state. Calculations based on such a model have been applied in many areas of chemistry and physics such as chemical reactions in condensed phases,* the decay of metastable state^,^ surface d e ~ o r p t i o nand , ~ diffusion of atoms or ions in solids5 or across membranes.6 Although activated processes occur in general on multidimensional potential energy surfaces, most studies have concentrated on one-dimensional systems. This is due in part to the computational complexities arising from the multidimensional formulation. Recent analytic works on one-dimensional barrier crossing include the kinetic theory approach of Skinner and Wolynes,’ the stable-state-picture formalism of Hynes et al.? the first passage time approach of Szabo et and the singular perturbation approach of Matkowsky and Schuss.l0 Generalization of these methods to multivariate systems is not trivial, and it is important for applications to realistic systems. Grote and Hynes” and van der Zwan and HynesI2 extended their formulation* to treat multidimensional barrier crossing dynamics in condensed phases for systems governed by the generalized Langevin equation. Brinkman,” Landauer and Swanson,I4 Langer,I5 and Blomberg,I6 have used the Fokker-Planck equation to derive a multidimensional version of Kramers’ rate theoryl for the high-friction regime. Of these, the formulations of Hynes et al. and Langer are of particular interest since they can be applied to the cases where hydrodynamic interactions should be taken into account; e.g., conformational transition rates in chain molecules are influenced significantly by the presence of hydrodynamic interactions.”J* The above-mentioned theories for multidimensional barrier crossing dynamics were derived with explicit or implicit assumptions that limit their validity to cases involving high barriers. ‘Supported in part by a grant from the National Science Foundation.
In this paper, we formulate a rate theory that is also valid for reactions involving relatively low barriers and explicitly accounts for the effects of hydrodynamic interactions on the reaction rate. Although we recognize that dynamical theory of multidimensional barrier crossing in the presence of long-time memory effects for the low to intermediate friction limit needs to be developed,I9 we will focus on the high-friction and Markovian limits in this paper. Such a description of dynamics is valid for reactions that involve the motion of large molecular units over a broad potential energy barrier in a viscous solvent. For example, many biochemical processes in proteins and lipids such as helix-coil transitionsZoand
(1) (2) (3) (4) (5) (6) (7)
Kramers, H. A. Physica (The Hague) 1940, 7, 284. Hynes, J. T. Annu. Reu. Phys. Chem. 1985, 36, 573. Hanggi, P. J . Stat. Phys. 1986, 42, 105. Jack, D. B.; Kreuzer, H. J. Phys. Reu. E 1982, 26, 6516. Dieterich, W.; Fulde, F.; Peschel, I. Adu. Phys. 1980, 29, 527. Rigos, A. A,; Calef, D. F.; Deutch, J. M. Eiophys. J . 1983, 43, 315. Skinner, J. L.; Wolynes, P. G. J . Chem. Phys. 1978, 69, 2143. (8) Northrup, S. H.; Hynes, J. T. J . Chem. Phys. 1980, 73, 2700. Grote, R. F.; Hynes, J. T. J . Chem. Phys. 1980, 73, 2715. (9) Szabo, A,; Schulten, K.; Schulten, 2.J . Chem. Phys. 1980,72,4350. Schulten, K.; Schulten, Z.; Szabo, A. J . Chem. Phys. 1981, 74, 4426. (10) Matkowsky, B. J.; Schuss, 2.S I A M J . Appl. Math. 1977, 33,365. Schuss, 2.SIAM Reu. 1980, 22, 119. (11) Grote, R. F.;Hynes, J. T. J . Chem. Phys. 1981, 74, 4465. (12) van der Zwan, G.; Hynes, J. T. J . Chem. Phys. 1982, 77, 1295. (13) Brinkman, H. C. Physica (Utrechr) 1956, 22, 149. (14) Landauer, R.; Swanson, J. A. Phys. Rev. 1961, 121, 1668. (15) Langer, J. S . Ann. Phys. (N.Y.) 1969, 54, 258. (16) Blomberg, C. Chem. Phys. 1979, 37, 219. (17) Pear, M. R.; McCammon, J. A. J . Chem. Phys. 1981, 74, 6922. (18) Ladanyi, B. M.; Hynes, J. T. J. Chem. Phys. 1982, 77, 4739. (19) Non-Markovian theories of activated rate processes were presented by Grote and Hynes.(ref 8 and 11) and more recently by Carmeli and Nitzan (Carmeli, B.; Nitzan, A. J. Chem. Phys. 1984,80, 3596 and their earlier works cited therein). The latter authors formulated a theory that is valid for the whole friction range but concentrated on one-dimensional systems. Reference 11 dealt with multidimensional systems but did not cover the small-friction limit and involved assumptions that limit the validity of the theory to highbarrier processes. It seems that a unified theory of non-Markovian activated rate processes in multivariate systems that is valid for the whole friction range needs to be developed. In this regard, the recent work of Straub et al. (Straub, J. E.; Borkovec, M.; Berne, B. J. J . Chem. Phys. 1986, 84, 1788; 1985,83, 3172) should be mentioned; it shows by a full stochastic simulation that the current theories of the non-Markovian one-dimensional problems may fail in certain friction ranges. See, also, Borkovec, M.; Berne, B. J. J . Chem. Phys. 1985, 82, 794; Zawadzki, A . G.; Hynes, J. T. Chem. Phys. Lett. 1985, lZ3, 476; Grote, R. F.; Hynes, J. T. J. Chem. Phys. 1982, 77, 3736; and Carmeli, B.; Nitzan, A. Chem. Phys. Lett. 1984, 106, 329.
0022-3654/88/2092-lO75$01.50/00 1988 American Chemical Society
1076 The Journal of Physical Chemistry, Vol. 92, No. 5, 1988 reactant j j barrier region
I!
II
region
!I product /i
region
Lee and Karplus semble of systems whose population on the potential energy surface. deviates from the equilibrium distribution at t = 0. We assume that after a transient period 7trany excess population in either the reactant or product region relaxes exponentially: N J t ) - W2 = SN,(O) exp(-A,,t)
bN,(t)
(2.1)
where N,(t) is the population of the reactant (a = R) or the is the corresponding product (a = P) region at time t , and equilibrium population. In the one-dimensional case, Northrup and Hynes22found that even for a barrier as low as 1.6 k B T( k B= Boltzmann constant, T = absolute temperature) there is a fairly well defined, although not complete, time scale separation that validates the above assumption: A r x ~ t r= 0.2. They also showed that T~~ is the time required for the barrier region population to reach its steady-state value and that the relaxation rate constant A,, is independent of the choice for the barrier region boundaries. In the multidimensional case, however, the height of the barrier alone cannot ensure the validity of (2.1). It is also required that the reaction coordinate should not be strongly coupled to modes that relax on time scales comparable to A,;]. We can then that A,, is independent of the choice for the barrier region boundaries by slightly generalizing the one-dimensional treatment of Northrup and Hynes. If the population of the barrier region remains constant for t > T ~ and , the law of detailed balancing holds for the rate constants kf and k, for the forward and reverse barrier crossing processes (Le., kf/k, = = Kq), (2.1) is equivalent to saying that the relaxation of population is governed by first-order kinetic^:^.^^
e
reaction coordinate
Figure 1. Stable-state picture of chemical reactions. The potential energy surface is divided into reactant and product regions separated by a barrier region. R* and P* denote activated states at the boundaries of the reactant and product regions on the reaction path, respectively.
enzyme-substrate reactions21properly fall into this category. Our theory is the extension of the one-dimensional formulation developed by Northrup and Hynes.22 We start with the coupled Langevin equations in a generalized coordinate system. When the coupled equations are written as a matrix equation, the couplings are seen to be due to off-diagonal components of the inertial, friction, and potential energy matrices. In the high-friction limit, the inertial terms in the Langevin equations may be neglected. Hence we can decouple the equations by applying a congruent transformation of the coordinate system that diagonalizes the friction matrix and the potential energy matrix simultaneously. By doing so, we also find a reactive mode that conforms to a variational criterion for the optimum reaction coordinate23in the neighborhood of the saddle point on the potential energy surface.. We then transform the resulting decoupled Langevin equations into an equivalent set of generalized Smoluchowski equations to obtain expressions for the forward and the reverse rate constants for the barrier crossing process. In section 11, we present a brief summary of the stable-statepicture formalism of Northrup and Hynes8-22and a derivation of the analytic expressions for the rate constants based on it. An application of the rate theory to conformational isomerization in n-butane is given in section 111. The effects of hydrodynamic interaction and flexibility of bonds and bond angles on the rate and the reactive mode are determined. Section IV states the conclusions of the present work. 11. Rate of Diffusive Multidimensional Barrier Crossing In this section we derive a steady-state rate expression for diffusive multidimensional barrier crossing. The basic framework corresponds to the one-dimensional formulation of Northrup and Hynes.22 According to their stable-state picture of chemical reactions,8 the full rate constants are expressed in terms of barrier rate constants and internal activation rate constants. We obtain the barrier rate constants by generalizing the Kramers' method] to multidimensional systems. Then generalizing the first passage time formalism of Szabo et a1.9 to multivariate systems, we obtain the internal activation rate constants. A . The Stable-State Picture (SSP) of Chemical Reactions. Although the SSP formalism8 is slightly more general than described here, the essential idea is illustrated in Figure l . The barrier crossing reaction is defined as the transition between the stable states R and P within the reactant and product regions. To observe the transition process, one usually starts with an en(20) McCammon, J. A.; Northrup, s. H.; Karplus, M.; Levy, R. M. Biopolymers 1980, 19, 2033. (21) Szabo, A.; Shoup, D.; Northrup, S. H.; McCammon, J. A. J . Chem. Phys. 1982, 77, 4484. (22) Northrup, S. H.; Hynes, J. T. J . Chem. Phys. 1978, 69, 5246. (23) Berkowitz, M.; Morgan, J. D.; McCammon, J. A.; Northrup, S. H. J . Chem. Phys. 1983, 79, 5563.
m/m dNR dt
_ _ -dNp - -kfNR + k,Np dt
with the relation, kf+ k, = L.It should be noted that, in contrast to A,, k f and k, depend on the definitions of the reactant and product regions as does the equilibrium constant Keq.2s In view of the above argument, the locations of barrier region boundaries can be chosen to simplify the dynamical problem for calculation of the relaxation rate constant Arx. If the forward and reverse rate constants are desired for some different definition of the reactant and product regions, they are readily obtained from the value of Kq that is evaluated or observed according to the new definition and A,, as calculated above. New dynamical considerations are not required. In the SSP the forward and reverse rate constants, kf and k,, are expressed in terms of barrier rate constants Kf and K, and internal activation rate constants K R and K ~ : kf = Kf[l + (Kf/KR) kr
+ (Kr/KP)l-'
= Kr[l + (Kf/KR)+ (Kr/KP)l-l
(2.3)
(2.4)
Here, Kf is a rate constant for the formation of P* from R maintained in internal equilibrium; K , has the corresponding ) a rate constant for the formation of R* (P*) definition. KR ( K ~ is starting from R (P) which initially has an internal equilibrium distribution. This condition for the calculation of KR and K~ implies that the stable states R and P in internal equilibrium are the reference initial and final states of the transition event.8,22This is not to say that development of internal nonequilibrium in R and P is neglected. Effects due to the internal nonequilibrium are reflected in the ratios Kf/KR and K,/Kp. In this regard, it is of interest to consider limiting forms of the above expressions for the full rate constants. When internal relaxation is rapid compared to the passage over the barrier (Le., Kf/KR, K,/Kp K f . Then, kf = ~ f ~ p / ( ~K,).p The factor K ~ / ( K ~ KJ measures the influence of the slow stabilization rate of P*. We note that for the reverse reaction the same factor measures the effect of depletion in the P* population due to the finite activation rate in the product region. E . Barrier Rate Constants. In this subsection we derive expressions for the barrier rate constants, K f and K,. We start with the coupled Langevin equations in Cartesian coordinates for a system with f degrees of freedom?
+
+
-
+
(2.5)
where i and j run over the particle coordinates (1 Ii, J IA, xi and m, are the Cartesian components of the position and the mass of the particle associated with index i , and dots over coordinates denote derivatives with respect to time. When hydrodynamic interaction is included, the friction tensor is a configurationdependent nondiagonal tensor. R” represents the randomly fluctuating force exerted by the surrounding medium and is described by a Gaussian distribution with mean and covariance: ( R N )=0
(2.6) (2.7)
From (2.5)-(2.7), we can derive the coupled Langevin equations in a generalized coordinate systemz7in which the potential energy function takes the simplest form and the potential energy profile along one of the generalized coordinates, say q,, is associated with the transition process from R to P. By this we do not mean that qr is the reaction coordinate. Only if there is no inertial or hydrodynamic coupling between qr and other coordinates is qr the reaction coordinate. Otherwise, the reaction coordinate does not follow qr strictly, as will be seen below. Here it must be noted that the overall translational and rotational degrees of freedom of the whole system need not and, in general, should not be separated out of the generalized coordinate system as was done As Helfand and SkolnickZ8noted in in Langer’s form~1ation.l~ their application of Langer’s theory to the calculation of polymer conformational transition rates, the advantage gained by including these unbounded modes is a great simplification of the programming of the calculations. However, the more important point is that the librational motion under investigation and the overall translational and rotational motions are in general coupled. In the freedraining limit (that is, when the hydrodynamic interactions between thef degrees of freedom of the system are neglected), only the overall rotational motion is coupled to the internal motions and the diffusive translational motion, referred to the center of diffusion, is not coupled to them.z9 However, in the presence of hydrodynamic interactions, all the degrees of freedom of the system are coupled. In the case of dihedral angle transitions in chain molecules, the motion along a dihedral angle is not a pure internal motion and is coupled to both translational and rotational motions even in the free-draining limit. Hence, the inclusion of overall translational and rotational degrees of freedom in the consideration of the torsional dynamics in chain molecules is essential rather than optional.30 If we include such unbounded modes, the barrier top is not the proper saddle point in the fullf-dimensional space and the saddle point description applies to a subspace spanned by the bounded modes. (26) Deutch, J. M.; Oppenheim, I. J . Chem. Phys. 1971,54, 3547. (27) Goldstein, H. Classical Mechanics; Addison-Wesley: Reading, MA, 1950. (28) Helfand, E.; Skolnick, J. J . Chem. Phys. 1982, 77, 5714. Skolnick, J.; Helfand. E. J . Chcm. Phys. 1980, 72, 5489. (29) Lee, S.; Karplus, M. J . Chem. Phys. 1984, 81, 6106. (30) Knauss, D. C.; Evans, G. T. J . Chem. Phys. 1980, 72, 1499.
The coupled Langevin equations in generalized coordinatesz7 are given by di (t ( daQj L ) - d Laqj = - E +aqrj Ry
( j = 1 , 2 ,...,j) (2.8)
where L = T - V(q), and the kinetic energy T is written as
3 is the Rayleigh dissipation function and is given by
The generalized random force on qj is given by (2.1 1)
with the relations (R,Q(t))= 0
(2.12)
(R?(t)Ry(t’))= 2 k ~ T ( 3 ~ ) , 6-( t’) t
(2.13)
Up to this point, we have not made any approximation except that we have assumed the validity of Langevin equations. Now we introduce several crucial approximations: (i) around the saddle point, we approximate the potential of mean force V(q) as V(q) = E*
+ (’/2)@P*q
(2.14)
(the saddle point is located at the origin of the generalized coordinate system) and (ii) we keep only the lowest nonvanishing contributions in the Taylor series expansions in writing the kinetic energy and the Rayleigh dissipation function, namely,
WO)
(2.15)
P ( q ) = 34(0)
(2.16)
Wq) =
Using matrix notation, we rewrite (2.8) in the form
P-ij
+ 3 4 4 + P * q = R4
(2.17)
We note that dynamic couplings among q-coordinates arise from the presence of off-diagonal components of Tq, 3q, and Vq. If we can find a coordinate system that diagonalizes the three matrices simultaneously, the Langevin equations are decoupled in that coordinate system. However, simultaneous diagonalization of three matrices is not possible in general. For a system made up of identical particles, the friction matrix 3 4 is proportional to the inertia matrix Ils in the free-draining limit, and thus the three matrices, 29, 3 4 , and P,can be diagonalized s i m ~ l t a n e o u s l y . ~ ~ In the diffusive regime, we can neglect the inertia term in (2.17) and write 34.q
+ Peq = R4
(high-friction limit)
(2.18)
We can then decouple the reduced Langevin equations in this limit by diagonalizing 3 4 and P simultaneously. The desired congruent transformation can be found by employing the same procedure that is used in the theory of small oscillations,z7~3z in which case F and P are diagonalized instead of 3 4 and P.First, applying a transformation by an orthogonal matrix B which diagonalizes the friction matrix
B3‘7.BB1.$+ B P - B B 1 . q = BR4
(2.19)
(31) Lee, S.; Karplus, M., unpublished results. (32) In ref 12, the diffusive normal modes were identified with the eigenvectors which diagonalize the nonsymmetric matrix 3-’.V, where 3 is the friction matrix and V is the second derivative matrix of potential energy around the saddle point. These eivenvectors are, in general, complex. Here, we find real diffusive normal modes by applying the same procedure used in the theory of small oscillations.
1078 The Journal of Physical Chemistry, Vol. 92, No. 5, 1988
we can write 3Y.y
+ VY*y = RY
(2.20)
Then, applying a transformation where y = B-lq and (W),,= by a diagonal matrix C whose components are given by (QIJ
= tl-1'261J
(2.21)
Lee and Karplus we assume that the distributions along all nonreactive modes are maintained in equilibrium. That is, expressing the potential of mean force [(2.30)] in the form f V(s) = Cvi(si) (2.31) i- 1
where
we get 2
+ P*z= R'
(2.22)
where z = C-l-y, VI = 2'-VY-C,and Rz = t?*RY. Finally, applying a transformation by an orthogonal matrix D which diagonalizes P,we get S
+ VS-s = Rs
(2.23)
where s = D-'.z and (VS), = X16,. Since VS is diagonal, (2.23) represents a set of decoupled Langevin equations for {s,]; that is, we have SI
+ X,s,
= R;
( i = 1, 2,
..., j )
(2.24)
From J2.12) and (2.13), we can show that the random force RS (=bC.d.Rq) satisfies the relations
( R X t ) )= 0
(2.25)
(RS(t)RS(t'))= 2k~T6,6(?- t')
(2.26)
The set of decoupled Langevin equations [(2.24)-(2.26)] implies that the probability of finding a particle at s on the barrier region at time t , P(s,t), can be written in the form P(s,t) = N f i P , ( s , , t )
(2.27)
1=1
Here N is the normalization factor for which the explicit expression is given in (2.39). P,(s,,t) should satisfy the generalized Smoluchowski equation since is Gaussian and purely randomd) that is c2.28) where Ji, the particle flux along the mode si, is given by
By solving (2.28) and (2.29) with appropriate boundary conditions, we can evaluate the barrier rate constants. We need first to find the saddle point reactive mode which plays a crucial'role in the evaluation. Since the friction tensor is a unit tensor in the s-coordinate system [Le., (39 = 6,], the reactive mode can be identified with the steepest liescent path in the s-coordinate system.23 The potential of mean force V [see (2.14)] as a function of is,] can be written in the form V(s) = E*
+ (y2)Fi.VS.s= E* + ( t / 2 ) ~ X , s I 2
(2.30)
1=1
where only one of the i l k , say XI, is negative.34 Hence, the correspnding diffusive normal mode sidefines the steepest descent path and is identified with the saddle point reactive mode. TO determine the probability density function P,(s,,t)along the nonreactive modes, (i = 2, 3, ...,A, (2.28) and (2.29) can be solved with the reflecting boundary conditions a t the barrier region boundary and some initial condition provided. The initial condition is specific to each problem under investigation. For simplicity, (33) Wang, M. C.; Uhlenbeck, G. E. Reu. Mod. Phys. 1945, 17, 323.
(34) The reaction path cannot go through a higher order saddle point where more than one of the X,'s are negative; see Pechukas, P. In Dynamics of Molecular Collisions, Part B; Miller, W. H.. Ed.; Plenum: New York, 1976; Chapter 6 .
As can be readily verified, these Pl's satisfy the Smoluchowski equations given by (2.28) and (2.29). In general, the equilibrium assumption for the distribution along all nonreactive modes may not be true if there are some modes that have relaxation times of the same order or longer than the inverse of the barrier crossing rate. In the case of the dihedral angle transition in chain molecules, there are modes such as overall twisting of the chain and end-to-end motion which have relaxation times longer than the dihedral angle transition time. In fact, in the presence of such long relaxation time modes, it may be difficult to define the transition rate unambiguously. To be rigorous, we must calculate the reaction rate separately for each configurational state in the subspace spanned by such modes and then take the statistically weighted average of the rates. Fortunately, the modes with such long decay times are longer wavelength modes (distortions spread over the distance comparable to the size of whole system) and therefore do not contribute strongly to the activation energy or other factors affecting the dihedral angle transition rate. The reaction rates obtained for different configurations of these nonreactive modes with long decay times are expected to be insensitive to the configuration. Then the assumption that all the nonreactive modes conform to equilibrium distributions can be practically justified. All these points are in fact related to the condition for the validity of (2.1). To calculate Kf and K ~ we , should solve (2.28) and (2.29) for P,(s,,r) with boundary conditions that conform to their definitions and an initial'conditipn provided.22 However, the problem can be put into a much easier form. : To calculate Kf, we imagine a steady state created from an equilibrium state in the following way: (i) all the particles at P, including P* are made invisible, (ii) on the way from R* to P*, only those particles that originate from R including R* are visible (this origin must be referred to the latest history of the particle since it may have been at P or P* for some time before it got to R or R*),(iii) particles are made invisible upon arriving at P*, and (iv) the invisible particles are made visible upon arriving at R * so that the reactant region maintains its equilibrium population that consists of the visible particles. Defining PI and J, in (2.28) and (2.29) by the population density and the flux of the visible particles, respectively, we can write
aP1
o=-=-at
d ~ , JI = const along sl as, -+
(2.35)
Rewriting (2.29) in the form'
and evaluating a line integralalong s1from sI(R*) to sl(P*),where R* and P* designate values at the boundaries of the reactant and the product regions on the reaction pathways (see Figure l ) , respectively [the method of evaluating sl(R*) and s,(P*) is described in the Appendix], we obtain
The Journal of Physical Chemistry, Vol. 92, No. 5, 1988 1079
Diffusive Multidimensional Barrier Crossing In the steady state under consideration, we have P l [ ~ l ( P * ) l= 0
(2.38)
we can get an expression for Kf that gives a more accurate numerical value than (2.47). We first rewrite Z* in the form
Here, ZR denotes the configuration integral over the reactant region given by
where qR ['(qp, q!, ..., q?)] are the generalized coordinates defined for the reactant region (see section IIC for details) and is the Jacobian associated with these coordinates [see (2.44)]. Substituting (2.38) and (2.39) into (2.37), we obtain
In the steady state we are considering, the flux across any surface s1 = constant [between s,(R*) and s,(P*)] is constant. Hence the flux across s, = 0 is equal to the reactive flux going into the product region. The barrier rate constant Kf is given by this flux divided by the population of the reactant region (= that is
w);l6
are the Jacobians
In (2.42) and (2.40), lJsl and ax,
ax,
where we multiplied both the numerator and the denominator by the same factor to recover a configuration integral over all barrier region modes in the numerator and then transformed the coordinate system back into the q-coordinate system. Only one of the vi's, say V,, can be negative for a simple saddle point.34 The final expression for the forward barrier rate constant Kf is now given by
where Z+ is a configuration integral over all degrees of freedom in the barrier region with the factor e-bE* separated out; that is z+ =
ax,
... -
lJql dql...dq,,dqr+l...dq,exp and (2.44) To the same degree of approximation as those made in (2.14) (2.16), we keep only the zero-order contribution to IJJ; that is, we write
Since q = B4.D.s [(2.18) - (2.26)] and Band D a r e orthogonal matrices, we have
Substituting (2.41) into (2.42), we have
where Z* is the configuration integral associated with the barrier region for all degrees of freedom orthogonal to the reactive mode and is given by
Z* =
s
lJSlods2...dsf exp
barrier region
[--E y2
Xis?
3
(2.48)
(2.47) together with (2.40), (2.46), and (2.48) gives the expression for the forward barrier rate constant K~ In practice, however, (2.47) may result in an inaccurate value of Kf especially when the number of degrees of freedom is large, since the relative error in the value of K~ is proportional to the sum of the relative errors in the values of the X i obtained from a numerical calculation. If we choose the q-coordinate system such that the potential energy matrix P in (2.14) is diagonal, that is (Ph, = VA, (2.49)
We note that, in the numerical evaluation of (2.51), only the relative error in X1 matters. Hence, (2.51) is expected to give a better value of Kf from the computational viewpoint than (2.47). In the same way, the reverse barrier rate constant K, is given by
1;:;
ds, exp( +PIX1lsl2/2)]-' (2.53)
where Z , is the configuration integral over the product region
C. Internal Activation Rate Constants. As introduced in section IIA, the internal activation rate constant KR ( K ~is) defined by the rate at which particles reach the reactant (product) boundary along the reaction pathways starting from R (P) in internal equilibrium. This rate is most easily calculated by using the first passage time (FPT)f ~ r m a l i s m .In ~ this subsection, we generalize the FPT approach to a multidimensional system. First we derive the expression for KR;the expression for K~ can be obtained in the same way. As in the evaluation of the barrier rate constant, we assume that (i) Around the minimum of the reactant potential well, the potential of mean force can be approximated as (2.55) V(qR) E R + '/2 ijR.WqR
(the position of the potential energy minimum is at the origin of the qR-coordinate system). (ii) W e keep only the lowest nonvanishing contribution in the Taylor series expansion around 0 in writing the Rayleigh dissipation function 3"R(qR)= P ( 0 ) (2.56)
Lee and Karplus
1080 The Journal of Physical Chemistry, Vol. 92, No. 5, 1988
Then neglecting the inertial term, the Langevin equations become 3R.qR
+
p . q R = RR
(2.57)
Applying a congruent transformation as has been done in section IIB, we find (we will hereafter omit the superscript R for designating variables associated with the reactant region in order to avoid notational complexities) SI
+ Xis, = R;
(2.58)
Here all the 1,’s are nonnegative. Although a single reactive mode cannot be selected out, in general, there are likely to be several diffusive normal modes which involve considerable movement along qr. These are the reactive-mode-like normal modes along which activation and deactivation will proceed. We denote these reactive-like modes by si and the nonreactive modes by sjn‘. The decoupled Langevin equations [(2.58)] leads to the result that the probability that a particle is at s at time t, given that it was at so [=(sIo, szo, ..., sfl)] at t = 0, has the form (2.59) i= 1
The Pi’s satisfy the Smoluchowski equation
(2.69)
Integrating both sides of this equation over s, we obtain [see (2.64)]
a
-\k(so,t)
at
J = E~+(s~~)\k(s,,t)
(2.70)
i=l
Finally, integration of this equation over all time gives the desired equation f -1 = CL+(si0)7(sO) (2.71) where we have used (2.65). The partial differential equation for by separation of variables:
7,
(2.71), can be solved
(2.72)
(2.61)
where Jj is the particle flux along the nonreactive mode si”’ [cf. (2.29)] and sY(R+) and sj”’(R-) are the values of s r at the reactant region boundary. The Pi’s for reactive-like modes are required to satisfy the absorbing boundary condition at R*
Pi(s;(R*),tls;o)= 0
This equation, which is the adjoint of the Smoluchowski equation [(2.60)], is called the (backward) Kolmogorov equation. The corresponding equation for the full probability P(s,tJso)is given by
i= 1
The Pis for nonreactive modes are required to satisfy the reflecting boundary conditiong
Jj(sy(R+),tls/?o’)= Jj(sjnr(R-),t)sj6) = 0
Pi(si,tlsio),regarded as a function of sio and t , satisfies the differential equation
(2.62)
The reason for the construction of this particular form of solution is that several reactive-like modes can make parallel contributions to the activation process in the reactant region. The magnitude ~ ) - ~ the activation rate along the mode si, and of T ~ ( s ~ measures the overall activation rate ~(s,,)-l is given by (2.72). Substituting (2.72) into (2.71), we obtain
and the reflecting boundary condition at the other reactant region boundary away from R*
J,(sf(R-),tls;O) = 0
(2.63)
The method used to evaluate the values of sf(R*) and s;(R-) is described in the Appendix. To calculate the first-order kinetic parameter K R , however, we need not have recourse to detailed time-dependent solutions of the Smoluchowski equations for the Pi‘s. The probability that a particle starting from so remains in the reactant region at time t is given by
where lJsl is the Jacobian associated with the reactant region s-coordinate system [cf. (2.43)]. The average residence time of the particle in the reactant region until it reaches R* for the first time (the mean first passage time from so to R * ) is given by 7(s0)
= &=dt t [ -:\k(s,,,r)]
= &-dr \k(so,t) (2.65)
-
where we have assumed that \k(so,t) decays faster than t-’ as t m . In fact, we presuppose the single-exponential decay of \k; i.e. W O J ) = exP[-t/T(so)l (2.66) which is equivalent to a first-order kinetic assumption since
J T -E(2.73) i=l7i
The second term in the summation on the left-hand side of (2.73) is much smaller than the first term especially deep inside the reactant region where the dominant contributions to T ] ’ S come from. The reason is that aT,/as,, has a very small value except in the vicinity of R*. This approximation becomes exact when there is only one s-mode that dominates the activation process, since in that case T , for the mode becomes nearly equal to T . We check the validity of the approximation numerically in section III (see Figure 4). Errors introduced by this approximation to the full rate constants, kfand k,, are found to be small. With neglect of the second term in (2.73) the equations for T , ( S , ~ )become
7i(sto) L+(Si0)Ti(Si(J= - SO) = -wi(so)-’
(2.74)
with f cwi=l
(2.75)
i= 1
The quantity Wi represents the relative contribution of the mode s, to the internal activation process and measures the projection
of the reaction Coordinate along the mode. The boundary conditions, (2.61)-(2.63), become: respectively (2.76)
We then identify the average of ~ ( 5 0 over ) the initial equilibrium distribution of the particle in the reactant region with the inverse of K R (see below). A differential equation for ‘(50) can be obtained in the following way. First, it can be shown35that the probability density function
(35) Stratonovich, R. L. Topics in the Theory of Random Noise; translated by Silverman, R. A,; Gordon and Breach: New York, 1963; Vol. 1, Chapter 4.
The Journal of Physical Chemistry, Vol. 92, No. 5, 1988 1081
Diffusive Multidimensional Barrier Crossing T ~ ( s ~ ( R *= ))0
(2.77) (2.78)
We note that Wi(so)is a much more slowly varying function of s,, and sio than the variation of T~ with sioespecially in the region where the main contribution to T ( s ~ )originates. Even in the vicinity of R*, Wiis expected to be quite smooth since both T~and T become small. Hence without much loss of accuracy we may assume constant values for the W:s. Then (2.74) together with the boundary conditions in (2.76)-(2.78) can be solved exactly. For reactive-like modes, we obtain Tj(Si0) = Wj-1 T p ( s i o )
(2.79)
where Tp(sio) =
pJ-YR"
du exp(PXiu2/2)
1"
of SSP reference initial and final states of the reaction to those for the real reaction event under i n ~ e s t i g a t i o nbecomes ~~ exact only if the internal nonequilibrium effects are negligibly small.s,22 This requires that the locations of the barrier region boundaries should be so chosen that a particle near the boundaries has a far greater chance of falling back into the potential wells than going over the barrier. The approximations made in section IIC [Le., the neglect of the second term on the left-hand side of (2.73) and the neglect of the dependence of Wi on so in (2.74)] can also be best satisfied by this rw irement on the locations of barrier region boundaries. We may rephrase this condition in the following way: in the vicinity of the reactant and product boundaries, the magnitude of an impulse due to the potential of mean force during a time interval At (the magnitude of At is comparable to the momentum relaxation time) is much larger than that due to the randomly fluctuating force. The magnitude of the impulse due to the potential of mean force along qr3' is given by
dw exp(-8Xiw2/2)
sXR-)
(2.80)
(2.87)
and Wj = 0 for nonreactive modes. Since we are interested in the average time to get to the R * zone starting from R in equilibrium, we average T ~ ( S #over ~ ) the equilibrium distribution. We have
where pr and 5, are effective mass and friction coefficients associated with the motion along qr. The magnitude of the impulse due to the random force along qr is given by
[2kBTPr]'/2 (2.88) Then, along qr there must be a point where the ratio of the former to the latter is a maximum; that is (2.89) (2.81) By substituting (2.79) and (2.80) into (2.81) and changing orders of integration we find Ti,q
=
(2.82)
6 - 1 Tg;
We define this point as the one where qrcrosses the barrier region boundary. If the variations in prand 5, along qr are much smaller than the variation in laV/aq,l, (2.89) reduces to the simple form
av = maximum Ian,/
where
(2.83) A reasonable estimate for the values ( Wi)is given by
The activation rate constant in the reactant region,
KR,
(2.90)
Since the error introduced by quadratic approximations of the potential energy may also be minimized by choosing the barrier region boundaries at the positions where the second derivative of the potential energy changes sign, (2.90) will give the optimum choice in many cases. 111. Conformational Isomerization Rate in n -Butane
is then (2.85)
where we have restored the superscript R to denote the quantity associated with the reactant region. The activation rate constant in the product region, K ~ is, given by exactly the same procedure and is expressed as (2.86) D . Optimum Locations of Barrier Region Boundaries. In calculating the barrier rate constants and the internal activation rate constants, we have approximated the potential energy as a quadratic function of coordinates around the saddle point and the minima of the reactant and product wells. The locations of the barrier region boundaries should be chosen such that the errors introduced by such approximations are minimized. However, there is another important aspect to choosing the proper locations of barrier region boundaries. The correspondence
The rate theory formulated in the previous section is applied here to evaluate the conformational isomerization rate in n-butane. Although its torsional dynamics is too simple to exhibit some interesting aspects manifested in the dynamics of longer chain molecules, it has been investigated by several authors with stochastic dynamics methods3840 and thus makes a good example to test the accuracy of the rate theory. In particular, the effects of hydrodynamic interaction and the flexibility of bonds and bond angles on the rate and the reactive mode are considered. A. Molecular Model and the Coordinate System. We employ the same molecular models as those in ref 38-40 to make possible a direct comparison with the earlier studies. We use the internal coordinate system as the q-coordinate system (see section 11). The first three coordinates (bl, 81, +1) define the Cartesian coordinates (36) The reference states for a real reaction event are in general stable states internally in a nonequilibrium steady state and depend on experimental specifics. See ref 8 and 22. (37) qr is the coordinate introduced in section IIB. (38) Helfand, E.; Wasserman, 2.R.; Weber, T. A. Macromolecules 1980, 13, 526. (39) Levy, R. M.; Karplus, M.; McCammon, J. A. Chem. Phys. Lett. 1979, 65, 4. (40) Montgomery, J. A,; Holmgren, S. L.; Chandler, D. J . Chem. Phys. 1980, 73, 3688. Montgomery, J. A,;Chandler, D.; Berne, B.J. J . Chem. Phys. 1979, 70,4056.
1082 The Journal of Physical Chemistry, Vol. 92, No. 5, 1988
i’
Lee and Karplus Here, it should be mentioned that the three minima of the potential energy function along the dihedral angle d4 are given in our convention by #4(gauche-) = 60’ #4(gauche+) = 300° (3.5)
#,(trans) = 180’
and that, precisely speaking, the q-coordinate system corresponds to the deviations of the internal coordinates from their values at the stationary positions:
(
=
:)iJ
2
(3.9)
can be calculated. The zero-order approximations to P and 9 are then given by N
P(0) =
\
(5) .M( 84
2)
0
with
(M), = misij
(3.10)
!
\ Figure 2. q-coordinate system used to investigate the torsional dynamics of n-butane.
In (3.1 l ) , the friction tensor [ is obtained by43
(x,, y,, zl) of the first carbon center with respect to the laboratory Cartesian coordinate system (called hereafter the xayoro-system):
(i;) (;;
(3.12)
bi sin Oi cos
=
sin
mi)
(3.1)
A new coordinate system is then introduced that is translated from the xggo-system by (x,, y,, zI) (i.e., it is located at the first carbon center) and rotated from the xoy,,zo-system by the Euler angles4’ (a = c#q, @ = 01, y = 0). We call this new Cartesian coordinate system the xlylzl-system. In the same manner, Cartesian coordinates of the ith carbon center and the associated Cartesian coordinate system are defined by ( b , Oi, &) with respect to the xi-lyi-lzi-l-system (see Figure 2). The six coordinates b,, O , , #,, 02, &, and & define the position and orientation of the molecule with a given conformation with respect to the xoy,go-system, and the others are the “internal” coordinates4* (in particular, 44is the dihedral angle, the librational motion of which is of interest and corresponds to the coordinate qr in the notation used in section 11). In terms of these internal coordinates, we can write the Cartesian coordinates T,, = (x,,yn,z,) of the nth carbon center with respect to the xGozo-system as rl = P1
where the 3 X 3 submatrices, Qij, are given by Qij
= F;’ 6,
7 + (1 - 6JA,
(3.13)
6, is the Kronecker delta, ti is the friction coefficient,of the ith vertex (extended carbon atom), and HiJ is the hydrodynamic interaction tensor between the ith and thejth vertices that is here approximated by the Oseen tensor44
(3.14)
or by the Rotne-Prager tensor4s -
r2 = A1-P2+ r l r3 = A1-AZ-P3 + r2
+
r4 = A,.A2-A3.P4 r3
(3.2)
where Pi and Ai are defined by Pi =
(
b, sin 0, cos 6, bi sin 0, sin bi)
bj
COS
Oi
cos B, cos 4, -sin di sin Bi cos mi cos Bi sin 6, cos bi sin 0, sin 6, -sin 0, 0 cos 0,
where r,, = r, - ri, rij = 131,uiis the hydrodynamic radius of the ith vertex, and qo is the solvent viscosity. In terms of the internal coordinates, the potential energy is assumed to be given by
(3.3)
(3.4)
(41) We define the Euler angles in the same way as Brink, D. M.; Satchler, G . R. Angular Momentum, 2nd ed.; Oxford University Press: London, 1968. (42) As mentioned in section IIB, these are not pure internal coordinates.
(43) Zwanzig, R. Adu. Chem. Phys. 1969, 15, 325. (44) Oseen, C. W. Hydrodynamik; Academisches Verlag: Leipzig, 1927. (45) Rotne, J.; Prager, S. J . Chem. Phys. 1969, 50, 4831.
The Journal of Physical Chemistry, Vol. 92, No. 5, 1988 1083
Diffusive Multidimensional Barrier Crossing
Ye(Bi)= (j/,)ke(cos Bi - cos
-
4")'
(3.18)
40.
0
-f! 30. P
= (j/~)lq(q$ - 4;)'
around the trans position
+ (!lz)k,(di- c$,B)~
.Ei
= E,
b
= E * - (/$c*(+~ 1 -
-20.
i? 10.
60.
120.
180.
240.
300.
360.
9, in degree
gauche-
Figure 3. Potential energy curve along the dihedral angle b4 which
I -0.6
0.0
kr
ki
ki
kr
z= trans t gauche+
(3.20)
In the temperature range of interest, the direct transition between two gauche states is almost impossible due to the high cis barrier (Ecis= 44.8 kJ/mol; see Figure 3). Barrier Rate Constants. The locations of barrier region boundaries for the transition of the dihedral angle 44 in n-butane were determined by using (2.90) (see Figure 3); the components of the inertia tensor and friction tensor along the dihedral angle are found to be constant. Following Skolnick and Helfand,26it is easy to show that for ~ ratios P / Z Rin (2.51) the transition of the dihedral angle 4 1 the and Z+/Zpin (2.53) are given by
defines the transition process. The dashed curves show the quadratic approximations of potential energy. $,(R-) = 156', $2 = 180', b4(R*) = 204'. $4* = 240', 44(P*) = 282O, $48' = 300', 44(P-) = 334'.
-1.0
around the saddle point
As can be seen from (3.16), no coupling among the coordinates is included in the potential energy. The potential parameters used in (3.17)-(3.19) are given in Table I. B. Isomerization Rate in the Free-Draining Limit. We are interested in the kinetics of the reaction
0.
0.
around the gauche positions
0.1
810
Figure 4. Plot to check the validity of the approximation made in section IIC [(2.73)]. The first term in the summation on the left-hand side of (2.73)is represented by the dashed curve and the second term by the solid curve. The dotted curve represents the values of exp(-PXisio2/2). See text.
TABLE I: Parameters of the System' system I b system 2< system 3d 300 300 425 1.0 x 105 1.0 x 105 1.0 x lo5 0.014 0.014 0.014 0.153 0.153 0.153 1.8 x 1 O l o 2.5 X lo9 6.6 X IO9 1.23 1.23 1.23 4.2 x 107 1.3 x 107 9.3 x io6
system 4e
298 1.5
x
0.014 0.153
m
1.23
io5
In the numerical evaluation of the configuration integrals we used the exact dihedral potential energy rather than its quadratic approximations [see (3.19)] to reduce the errors arising from the quadratic approximation. In Table 11, we list the barrier rate constants calculated from (2.51) and (2.53) for the four systems given in Table I, along with the representations of the reactive modes in the barrier region in terms of internal coordinates (see the Appendix). All the degrees of freedom are mixed into the reactive mode. Of these, the concerted movements of B3, B4, ?,, and in particular are crucial for the effectiveness of the reactive mode. Activation Rate Constants. The activation rate constants for the trans and gauche potential wells were calculated by applying the procedure described under section IIC. The activation rate constants calculated from (2.85) and (2.86) and the representations of the activation modes in terms of the internal coordinates (see Appendix) are given in Table 111. The activation rate constants are much greater than the barrier rate constants (at least by a factor of 10). Hence the conditions for the optimum locations
unit
K ns-l kg/mol nm ns-' rad J/kg
OThe potential parameters for the dihedral angle, tJ4, in all four systems are the same as those in Ryckaert, J.-P.; Bellemans, A. Chem. Fhys. Lett. 1975,30, 123. They are tabulated in ref 28 and 38 and will not be reproduced here. *This is the most realistic parameter set among the four systems (see ref 38).
