J . Phys. Chem. 1988, 92, 3261-3267
3261
Stochastlc Molecular Dynamics Study of Cyclohexane Isomerization Robert A. Kuharski,+David Chandler,* John A. Montgomery, Jr,,t Faramarz Rabii, and Sherwin J. Singerg Department of Chemistry, University of California, Berkeley, Berkeley, California 94720 (Received: March 23, 1987: In Final Form: September 9, 1987)
Multidimensional BGK stochastic molecular dynamics calculations of reactive flux correlation functions have been performed with modified versions of the Rckett-Strauss potential and a variety of velocity randomization frequencies. The results provide an explanation of the experimental observations on the chair-boat interconversion rate as a function of solvent pressure or viscosity. Since the solvent contribution to the activation energy is positive and increases with increasing pressure, transition-state theory predicts that the rate will decrease with increasing pressure. The actual rate, however, increases. The stochastic molecular dynamics calculations show that the transmission coefficient increases with increasing pressure because additional coupling to the environment enhances the rate of stabilization of molecules after they have passed through the activated transition state. In other words, for the experimentally relevant regime, energy rearrangement between the molecule and the stochastic bath occurs with similar ease as the rearrangement among intramolecular modes. We conclude that cyclohexane isomerization in liquids is an activated process that does not coincide with the RRKM picture of unimolecular kinetics. The sensitivity of this conclusion to details in the intramolecular potential energy function is discussed.
I. Introduction This article presents the results of a stochastic molecular dynamics study of the chair-boat isomerization of cyclohexane. The results provide an explanation of the high-pressure liquid-phase measurements of interconversion rates by Jonas and co-workers.' They also bear on the issues of intramolecular vibrational energy transfer and suggest that coupling to a solvent environment alters significantly the course of this transfer from what might be anticipated in gas-phase unimolecular kinetics. The nature of intramolecular vibrational energy rearrangement or relaxation (IVR) is central to chemical dynamics of polyatomic molecules.2 If complete and swift, this relaxation leads to the RRKM theory of unimolecular kinetics,3 and the failure to complete the rearrangement during the course of an activated process can be viewed as a breakdown in the RRKM theory. In general, not every mode of a polyatomic molecule is relevant to the kinetics of a unimolecular reaction since the coupling between modes of highly disparate frequencies will generally be very weak if not negligible. As a result, one often regards as adjustable the number of relevant intramolecular modes. This adjustability is somewhat arbitrary, and it is possible that as the environment of a molecule changes, the number of relevant degrees of freedom may also change. Indeed, energy may.pass between the molecule and the environment as fast or faster than it can rearrange among the intramolecular degrees of freedom. If this possibility transpires, as noted by Borkovec, Straub, and Berne,4 then the effective dimensionality of a molecule may be lower when the molecule is in a liquid than when it is in a gas. Since this effective dimensionality or relevant number of modes may change as the state of the liquid changes, the two types of processes-intramolecular and intermolecular-can be competitive, and the assumption of full randomization of energy among a predetermined number of intramolecular modes can be an untenable assumption. The coupling between a polyatomic molecule and a condensed environment can therefore lead to a breakdown in RRKM interpretation of unimolecular reaction dynamics. Consider why this issue is of importance in the interpretation of the measured chair-boat interconversion rate for cyclohexane in different solvents: Jonas and co-workers' studied this interconversion rate and found that it increased when the pressure of the solvent was increased. In their interpretation, the observations verified the theoretical predictions5s6that the transmission coef'Permanent address: S-CUBED, P.O. Box 1620, La Jolla, CA 92037. * Author to whom correspondence should be addressed. 'Permanent address: United Technologies Research Center, East Hartford, CT 06108. 5 Permanent address: Department of Chemistry, Ohio State University, Columbus, O H 43210.
0022-3654/88/2092-3261$01.50/0
ficient for isomerization in a liquid would be an increasing function of pressure. The theoretical prediction was based initially on the results of a one-dimensional stochastic trajectory s t ~ d y .Sub~ sequent multidimensional calculations for the isomerization of +butane6,' also suggested that the transmission coefficient for isomerization in liquids could be an increasing function of pressure or viscosity. Nevertheless, the interpretation of the cyclohexane experiments was questioneds primarily because simple models exploiting the RRKM approximation seemed to imply that the transmission coefficient for most isomerizations should decrease with increasing liquid pressure. Indeed, it was s ~ g g e s t e d that *~~ the high dimensionality and vibrational energy transfer between modes of a polyatomic molecule would suppress the so-called inertial or Lindemann regimes of unimolecular kinetic^.^ Due to this suppression, for cyclohexane, one would observe in a liquid only the diffusive regime first elucidated by Kramers.Io In other words, when the coupling of a reaction coordinate to an intramolecular bath is sufficiently-strong-or, equivalently, when the dimensionality of the reaction surface near the transition state is sufficiently high-the introduction of further couplings to a fluctuating medium should only lower the transmission coefficient. The assumption of strong coupling between intramolecular modes is central to the validity of the RRKM picture of unimolecular kinetics. If this picture was correct for cyclohexane, it would mean that the experimentally observed increase in rate with increasing pressure was an equilibrium phenomenon that had little to do with intrinsic activated dynamics. In particular, it would mean that the transition state theory approximation to the rate constant was an increasing function of pressure, and increasing (1) (a) Hasha, D. L.; Eguchi, T.; Jonas, J. J. Chem. Phys. 1981, 75, 1571. (b) Hasha, D. L.; Eguchi, T.; Jonas, J. J . Am. Chem. SOC.1982, 104, 2990. (c) Hasha, D. L. Ph.D. Thesis, University of Illinois, 1981. (2) See, for example: Oref, I.; Rabinovitch, B. S. Acc. Chem. Res. 1979, 12, 166. McDonald, J. D. Annu. Rev. Phys. Chem. 1979, 30, 29. (3) Hinshelwccd, C. N. The Kinetics of Chemical Change; Oxford Press: Oxford, England, 1940. Robinson, P. J.; Holbrook, K. A. Unimolecular Reactions; Wiley-Intersciences: New York, 1972. (4) (a) Borkovec, M.; Berne, B. J. J . Chem. Phys. 1985, 82, 794. (b) Borkoveck, M.; Straub, J. E.; Berne, B. J. J . Chem. Phys. 1986,85, 146. (c) Straub, J. E.; Berne, B. J. J . Chem. Phys. 1986, 85, 2999. (5) Montgomery, J. A,, Jr.; Chandler, D.; Berne, B. J. J. Chem. Phys. 1979, 70, 4056. (6) Montgomery, J. A., Jr.; Holmgren, S. L.; Chandler, D. J . Chem. Phys. 1980, 73, 3688. (7) Rosenberg, R. 0.; Berne, B. J.; Chandler, D. Chem. Phys. Lett. 1980, 75, 162. (8) Hynes, J. T. J . Star. Phys. 1986, 42, 149. Nitzan, A. J. Chem. Phys. 1985, 82, 1614. Nitzan, A. J . Chem. Phys. 1986, 86, 2734. A closely related discussion though not specifically aimed at the cyclohexane experiments is presented by Shroeder, J.; Troe, J. Chem. Phys. Letr. 1985, 116, 453. (9) Zawadski, A. G.; Hynes, J. T. Chem. Phys. Lett. 1985, 113, 476. Zawadski, A. G.; Hynes, J. T., preprint 1985. (10) Kramers, H. A. Physica (Amsterdam) 1940, 7 , 284.
0 1988 American €hemica1 Society
3262 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
Kuharski et al.
where ri denotes the position of the ith carbon (the “center” of the CHI group), and the torsion and bond angles, t i and $i, respectively, are defined in Figure 1: (2.lb)
‘.V3
Figure 1. Coordinates for the cyclohexane molecule. The ri)s for 1 I i 5 n are the positions of the carbon atoms, f i is the torsion angle about the ith bond vector bi = ri+l - ci, the angle between bonds i-1 and i is J.r, and cos-] a2is the angle between the 3-4-5 plane and the 1-3-5 plane. a1and a3 (not explicitly pictured here) are similarly defined in terms of the projections of the 1-2-3 and 5-6-1 planes, respectively, onto the 1-3-5 plane.
so rapidly as to counteract the assumed decrease in the transmission coefficient. Recent calculations, however, have shown that the opposite is the case: the solvent contribution to the free energy of activation is positive and an increasing function of pressure.” Therefore, the experimental observations do indeed imply that, for cyclohexane, the transmission coefficient for isomerization in liquids is an increasing function of pressure or coupling between solvent and solute. The experiments therefore reveal a breakdown in both transition-state theory and the RRKM picture of unimolecular kinetics. What remains to be done, however, is to perform a quantitative study of the dynamics of this system to bring some form of theory into accord with experiment. Such a study employing a stochastic molecular dynamics model is described in this paper. Cyclohexane is an especially convenient system to examine because so much is known empirically about the relevant potentials of interaction. Pickett and Strauss proposed a useful model for the potential energy surface and an ingenious identification of the reaction coordinate for the chair-boat isomerization.I2 Other potential surface models are available, t00.l~ In the work presented herein, we have adopted the Pickett-Strauss coordinates and a modified version of their potential energy function. Reasonable interaction potentials between C H 2 groups and solvent molecules have been deve10ped.l~ We did not need this latter information in the current study since we consider only a stochastic caricature for the solvent. Knowledge of the intermolecular potentials was required, however, in the calculations of the solvent contribution to the free energy surface,“ and they will be important if and when one chooses to consider a dynamics simulation or calculation incorporating the molecular details of the solvent. Our modifications of the Pickett-Straws potential are presented in section 11, and the stochastic model of the environment is also described there. The reactive flux method for computing the transmission coefficient for this system is outlined in section 111. The results of our calculations are presented in section IV. Along with providing an explanation of the experimental results on cyclohexane, we compare the behaviors of reactive flux correlation functions for different models. These latter results should be useful when approximate theories for energy flow in nonlinear dynamical models are tested.
COS {,+I
=
X
b,+l)*(b,+l X bi+2)/h
(1 1) Singer, S. J.; Kuharski, R. A.; Chandler, D. J . Phys. Chem. 1986.90, 6015. (12) Pickett, H. M.; Strauss, H. L. J . Am. Chem. SOC.1970, 92,7281. A tutorial on the Pickett-Strauss conformational coordinates is given by Strauss, H. L. J . Chem. Educ. 1971, 48, 221. (13) Allinger, N. L. J . Am. Chem. SOC.1977, 99, 8127. Burkert, U.; Allinger, N. L. Moleculor Mechanics; ACS Monograph 177; American Chemical Society: Washington, DC, 1982. (14) See, for example, ref 1 1 and references cited therein.
bi+~Ilbi+~ X b1+21 (2.lc)
The equilibrium bond angle, q0,is 112.6O = 0 . 6 2 5 6 ~rad, and the constants h,f, w , and u have the values 125 kcal/(mol.rad2), -10 kcal/(mol.rad2), 5.044 kcal/mol, and 3.288 kcal/mol, respectively. Along with ignoring the protons, the PS potential leaves out C-C bond vibrations, and the height of the barrier between boat and chair conformers is slightly in error. We have modified these latter two features as described below, but first consider the coordinate characterizing the stable conformers. Pickett and Strauss observed that with their potential, the conformational energetics are dominated by two collective variables 0 and $ defined according to tan 0 = Y / X
(2.2a)
tan 4 = 3Il2(a3- cul)/(2cu2 - el - a3)
(2.2b)
with (2.3a)
and 013
x={
cos 6 - a2 cos (4
or a1 cos 4 - a* cos (4
+ + */pr, 4/gr)
(2.3b)
As depicted in Figure 1, cyl is the angle between the 1-3-5 and 1-2-3 planes; similarly a2and a3are the angles between the 1-3-5 and 3-4-5 planes and the 1-3-5 and 5-6-1 planes, respectively. The alternatives in (2.3) are equivalent, except at singular points where one choice is to be preferred over the other. Of the two primary angles, 0 is by far the most significant. The , and T, boat and two chair states are located near 0 = ~ / 2 0, respectively; the barriers are near 0 = 0* = 1 3 ~ / 4 5and T - 0*. We have utilized their coordinate 0 since the time evolution of step functions like ~ , [ e ( t ) l= 1,
o < e ( t ) < 0*
= 0, otherwise
(2.4)
provide a convenient measure the rates of transitions between stable conformers. Our choice of 0* = 52’ is based on the free energy calculation described below which gives a barrier maximum at that angle. The results of our dynamical calculations are independent of the convention we use-employing a transition state at the free energy barrier or at the nearby potential energy barrier will give the same results for the rate constants. In our modifications of the PS potential, we have added bond vibration potentials 6
11. Dynamical Model
The Pickett-Strauss (PS) potential for cyclohexane isI2
X
i= 1
with bj = (ri+l - ril, where r7 = rI, and 1
~ ( b =) - k ( b - bo)’ 2
+ Au(b)
Here, Au(b) is the deviation from the harmonic potential of a Morse potential with the same equilibrium bond length and harmonic force constant. The equilibrium bond length is bo = 1.54 A. The choice k = 662.4 kcal/(mol.A2) seems appropriate and that is the choice considering the frequency of C-C ~tretches,’~
Molecular Dynamics Study of Cyclohexane Isomerization
20.
1
The Journal of Physical Chemistry, Vol. 92, No. I I, 1988 3263 velocity is taken randomly from the Maxwell-Boltzmann distribution ~ M B ( U )a
I
I
e
77
Figure 2. Free energy (in units of k B T )of the cyclohexane molecule at the temperature 300 K as a function of the angle 8. It was computed from a Monte Carlo sampling of the Boltzmann factor for the potential energy V(r, ,...,r,) in (2.8).
we have used in most of our calculations. In some cases, however, we also performed calculations with an artificially softened C-C bond where kWft= 221.7 kcal/(mol.A2). Finally, the anharmonic constant in Morse potential, Au(b), was fit by considering the dissociation energy for the C-C bond.I6 This gave Au(b) = g(b - bo)3+ ..., where g = 629 kcal/(mol-A3). The results of our calculations performed with this potential were indistinguishable from those with Au(b) = 0. In other words, anharmonic terms in the C-C stretch potentials seem to be irrelevant for the activated unimolecular dynamics we consider herein. In a few instances, we have also studied the activated dynamics for the model in which the C-C bonds are rigid. We will see that bond vibrations do have a quantitative effect on the transmission coefficient and associated correlation functions. The other modification we have made concerns the barrier heights. Since the activation barrier in the PS potential is slightly high, we have added and use E = 2.0 kcal, y = 9/rad2, and 8* = 1 3 ~ / 4 5 When . Vba, and Vvibare added to Vps, the net potential energy has a barrier height in agreement with e ~ p e r i m e n t ,and ' ~ the behavior in the vicinity of the stable states is essentially indistinguishable from the V,. The total potential energy function we have adopted is therefore V(rl,...,r,J = VppS(rl,...,r,)
(2.10)
where m is 14 amu (the mass of a CH2 group). After assigning the new velocity and retaining all other phase space variables as unaltered from their values just before interruption, the Newtonian trajectory is started again and run until the next interruption. The times at which these instantaneous velocity alterations occur are identified at random according to a Poisson distribution with a = (average interruption frequency per carbon) (2.1 1 )
\I 0.
0
exp[-(l /2)pmu21
I
+ vvib(bl,***,bn)+ Vbar(@)
(2.8)
If vb,,(8) is omitted, we find some changes in the quantitative behavior for the transmission coefficient but nothing that changes our primary qualitative conclusions. The free energy F(0) is obtained from the Boltzmann integration exP[-PF(e)l = S d r , ... S d r , 6[e - 8(r,,...,r,)] exp[PV(rl ,...,r,)] ( 2 . 9 ) where (@ke)-'is the temperature T . The integration was computed by Monte Carlo umbrella sampling, and the free energy function obtained in this way is drawn in Figure 2. The calculation used the C-C force constant k. When ksoftor rigid C-C bonds are employed, the F(0)'s obtained are nearly identical with the one shown in Figure 2. To study interconversion rates, we examine Newtonian trajectories governed by this potential energy (2.8). Included in the trajectories is a stochastic element: After progressing from an initial phase space point for some period of time, the trajectory is interrupted and a new velocity is given t o the ith carbon. This (15) Infrared Handbook; Szymanski, H. A., Ed.; Plenum: New York, 1963. (16) A value of 83.2 kcal/mol for the C-C bond is given, for example, by Laidler, K. J.; Meiser, J. H. Physical Chemistry; Benjamin/Cummings: Menlo Park, CA, 1982. (17) With the specified parameters, the free energy barrier of our modified PS potential is 10.5 kcal/mol. Experimental estimates of the barrier height are all approximately 10 kcal/mol. See ref l b and also, for example Allerhand, A.; Chen, F.; Gutowsky, H. S.J . Chem. Phys. 1965, 42, 3040.
At a particular interruption point in time, the carbon for which the velocity is randomized is chosen at random and uniformly from among the six carbons in the molecule. This stochastic molecular dynamics scheme is often called the BGK model after the kinetic equation due to Bohm and Gross.I8 The solutions to that equation would give the same correlation functions as those determined by averaging over an ensemble of trajectories generated by the stochastic molecular d y n a m i c ~ . ~ J ~ The parameter a-l is the relaxation time for the velocity of a carbon atom in this stochastic model. In reality, there is not simply one time but a variety of time scales governing the solvent induced velocity relaxation of an atom in a solvated polyatomic molecule. Experience with the BGK model for normal alkanes, however, suggests that if we confine our attention to intramolecular dynamics, the role of a solvent is sufficiently characterized with just one parameter. Indeed, for n-butane as simulated in a Lennard-Jones fluid,7 the trans-gauche reactive flux correlation function is well reproduced by the BGK model6 with a = 0.2 X 1013 s-l. In general, since velocity relaxation is the result of fluctuating forces, and since viscosity is a measure of these forces, we anticipate u i= cv (2.12) where 7 is the solvent viscosity and c is a constant sensitive to the particular molecule under consideration but is relatively independent of that molecule's environment. Estimates based on Stokes law would give the momentum u therefore a randomization frequency relaxation time M / 2 ~ q and , u is the collision diameter and M is the of about 2 ~ 7 7 a / Mwhere total mass of the molecule. The constant c obtained from this Stokes law frequency can be much higher than the value we identify empirically. Of course, the relevant fluctuating forces for intramolecular relaxation are at much higher frequencies than the zero-frequency center of mass friction considered in Stokes law. This observation is in much the same spirit as those leading to models employing frequency-dependent friction. In our approach, however, the constant c is a parameter that we simply fit when comparing our stochastic molecular dynamics results with those of experiment. When we do fit this one adjustable parameter (in section IV), we obtain a value of c
= 1 x 10" s-1 cp-'
(2.13)
for the cyclohexane molecule. As in the n-butane this value for c is much less than the value that one might predict with Stokes law. The particular value given in (2.13) for cyclohexane is also less than that apparently appropriate for n-butane. We stress, however, that changes in c of one order of magnitude do not alter our primary conclusions. In summary, we use a BGK stochastic molecular dynamics model with modified versions of the Pickett-Strauss potential. While a variety of values for u have been examined and will be reported in section IV, those of primary importance for comparing with experiment are in the range 10" s-' 5 01 5 l O I 3 S-I. The randomization or interruption frequency, a,is much lower than would perhaps be predicted from a Stokes law estimate. It is (18) Bohm, D.; Gross, E. P. Phys. Reu. 1949, 75, 1864. Bhatnagar, P. L.; Gross, E. P.; Krook, M. Phys. Reu. 1954, 94, 5 1 1 . (19) Chandler, D. In Stochastic Molecular Dynamics, NRCC Proceedings No. 6; Ceperley, D., Tully, J. Eds.; USDOE and NSF: Washington, DC, 1979.
3264 The Journal of Physical Chemistry, Vol. 92, No. 1 I, 1988
Kuharski et al.
noted, however, that the efficiency of intermolecular couplings in altering intramolecular dynamics is both system specific and most likely far less than that for altering low-frequency diffusive motion of a molecule. 111. Method of Reactive Flux Simulations The dynamics of activated events are conveniently analyzed with the so-called "reactive flux" correlation function
k(t) =
(B(o) s[e(o) - e*]
~[e(t)])
(3.1)
-
= k(t)/k(o+),
0
4: \
n
4
W
4:
o* < q t )
0 and zero otherwise.
IV. Results Our results for the transmission coefficient of cyclohexane as a function of the randomization frequency, a,are plotted in Figure 4. They cover a wide range beginning in the inertial low-coupling (small a) regime and extending well-beyond the Kramers turn over into the diffusive regime (large a). The inertial or lowcoupling regime is most important for comparison with experiment, and so most of our data was collected there. An enlargement of that region is also depicted in the figure. The corresponding reactive flux correlation functions are illustrated in Figure 5 . Notice that, in general, the transient behavior is more complicated than a simple exponential relaxation. The behavior of the correlation function when a = 0, denoted by ko(t),can be used in a strong collision a p p r ~ x i m a t i o nto~ ~provide a theory for the k ( t ) at finite but small a. This approximation pictures dynamics in a fashion similar to the BGK model, but the interruption events thermally randomize both velocities and positions within a stable state. As a result, such an approximation does not describe the onset of spatial diffusion and is therefore 1. Nevertheless, it can incorrect at large a where it predicts K
-
( 2 3 ) Forst, W. Theory of Unimolecular Reaction; Academic: New York, 1982.
+ s ) i o ( s + a) + ai& + a)
(a
+ a/s)k,(s + a )
(4.1)
where i ( s ) denotes the Laplace transform of k ( t ) , and the approximate equality is virtually exact owing to the separation of time scales; Le., the activation barrier is high and deactivation is much more rapid than activation. The use of the approximate equality when performing the inverse Laplace transform introduces no significant error for times smaller than or as large as those required to obtain the plateau value of k ( t ) . An extension of this approximation employs different randomization frequencies for different stable states. The algebra involved in implementing this generalization is somewhat complicated but simplifies considerably for situations where multiple barrier recrossings are neglected. We use this simplification since, in view of ko(t)shown in Figure 6 , we see that it is accurate to better than 10%. With its incorporation, we find
k,,,(s) = (1
+ a+/s)R,c+)(s+ a+)+ (1 + a-/s)io’-’(s + ff-) (4.2)
where the subscript “gsc” denotes “generalized strong collision” approximation. The frequencies a+ and a- refer to the randomization rates in the boat and chair regions, respectively. The first term on the right-hand side of (4.2) is the gsc representation of the Laplace transform of ko(+)(t);similarly, the second term represents the Laplace transform of kJ-)(t). We have inverted the Laplace transforms numerically to compare the strong collision approximations with the BGK simulation results. Representative comparisons are shown in Figure 7 where we have taken a+ to be equal to the BGK a and a- = 3a+. We arrived at this choice by noting empirically that, with it, the first and second terms in (4.2) accurately reproduced the BGK k(+)(t) and k ( - ) ( t ) ,respectively, in the small a regime. No such close correspondence, however, could be found with any choice of a (24) Skinner, J. L.; Wolynes, P. G. J. Chem. Phys. 1978, 69, 2143. The specific formula we site, the first of (4.1), is discussed in ref 4b.
3266 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
Kuharski et al.
TABLE I: Transmission Coefficient, K, for Cyclohexane' a, 10"
K
(stochastic MD)
K
(W)*
0.01 0.408 f 0.030 0.407
0.02 0.432 f 0.027 0.441
0.05 0.1 0.2 0.487 f 0.543 f 0.616 f 0.021 0.019 0.008 0.496 0.548 0.610
"Trajectories performed with the modified PS potential. (2.8). a-k,(a-) / k ( O + ) .
s-l
atom-'
0.3 0.5 1.o 1.5 2.0 5 .O 10.0 0.664 f 0.722 f 0.750 f 0.746 f 0.728 f 0.555 f 0.369 f 0.008 0.009 0.013 0.006 0.013 0.022 0.023 0.667 0.728 0.827 0.876 0.905 0.952 0.970
bCeneralized strong co~~ision approximation, KgSC = a + l ; , + ( a + ) / k ( ~++ )
1.o
+-
0
W
3 \ 4
W
3
+-
3 \
*
0.5 0.0 0
40
80.
a= 0.05x10'~sec-'
0
v
2
l*O
n 3 W
3 2
1.o
0.5
0.0
0.0
t (PS4 Figure 7. Comparison of the generalized strong collision approximation
predictions (solid line) and the single-frequency strong collision approximation (dashed line) with the results of the stochastic molecular dynamics simulations (circles). These curves represent reactive flux correlation functions for cyclohexane at T = 225 K and two different values of the randomization frequency. when the simplified strong collision approximation, (4. l), was employed. Figure 7 illustrates the comparison with a set equal to the BGK a. Evidently, stabilization from stochastic collisions is 3 times more effective in the chair regions than in the boat regions. This observation could be understood if "collisions" with only one particular CH2 are effective in removing energy from the reaction coordinate when the molecule is in a boat configuration, but "collisions" with any one of three CH2 groups are operative when the molecule is in a chair conformer. Perhaps this is an idea worthy of future investigation. Numerical results for the transmission coefficient are given in Table I and compared with the strong collision theories. The failure of strong collision approximations to yield correct results in the large a regime can be avoided by employing a connecting formula similar to those proposed by Berne and co-workers4 and In particular, we write by Hynes and co-workers, K-'
i=
Kgscl+ a / W b
(4.3)
where Kgsc
+
= [a+ko'+'(a+)
CY&(-)(CY-)]
0.5
/k(o+)
(4.4)
is the transmission coefficient implied by (4.2), and abis a frequency associated with the barrier. We determine this barrier frequency by fitting (4.3) to match the value of K found at the highest randomization frequency considered in our simulations, a = lOI4 s-' atom-'. This empirical procedure gives wb = 5.9 X lOI3 s-', and the results of (4.3) graphed in Figure 4 use this value. With the aid of the connecting formula and the trajectory results for ko(t),we can quantify the qualitative discussion presented in the Introduction concerning RRKM or Lindemann-Hinshelwd behavior. In Figure 6 we see that ko(t) exhibits a short-time relaxation and a much longer time relaxation. The former persists for a few tenths of a picosecond. The latter corresponds to RRKM relaxation, and in the limit of low a,with either (4.1) or (4.2),
0.
2.
4.
Figure 8. Transmission coefficient for cyclohexane at T = 225 K as a
function of solvent viscosity. The lower panel expands the scale for the low-viscosity regime by employing a semilog plot. The experimental results (ref 1) are for three solvents: carbon disulfide (crosses),acetone (circles), and methylcyclohexane (squares). Trajectory results are represented with the solid line. this long-time relaxation will dowinate the rate process, giving the Lindemann-Hinshelwood formula for the rate constant. Comparison between gas-phase experimentsZ5and the strong collision model lets us estimate this long time as 300 ps. The short-time relaxation, however, is representative of non-RRKM behavior. With both the short- and long-time relaxations present, we find the turnover from inertial to diffusive behavior to occur at a = 10 ps-' as shown in Figure 4. Suppose, however, the non-RRKM short-time relaxation was actually absent. With the connecting formula we would then predict the turnover at a = 0.4 ps-I. If only a 20% contribution to ko(t) came from the short-time non-RRKM relaxation, we would predict a turnover at about 4 ps-l. It is essentially the neglect of this short-time non-RRKM relaxation that led to suppose that the turnover would not occur at randomization frequencies consistent with the behavior of cyclohexane in a liquid environment. The consequences of the non-RRKM behavior as we have found them here for cyclohexane were first hypothesized by Borkovec, Straub, and Note also that the non-RRKM relaxation in k , ( t ) is highly system dependent (see below). Therefore, the precise location of the turnover is sensitive to the details of the intramolecular potential energy function. Any correct quantitative predictions concerning the turnover must thus be based upon specific information concerning the potential function. These remarks presuppose that the BGK trajectory results do coincide with those of experiment. A comparison between our calculations and experiments are illustrated in Figure 8. The comparison was congtructed as follows: First, judging from the experimental results' for the rate constant, k , and our earlier (25) Ross, B. D.; True, N . S. J . Am. Chem. Sor. 1983,105, 1382,4871.
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3267
Molecular Dynamics Study of Cyclohexane Isomerization I
I
I
I
1
I
I
1.o
I
0
W
41 \ n
+n
0.0
0.6
0
W
4
a= o . ~ x I o ' ~ ~ ~ c - '
W
41
\
41
3
0.0
0.2
w
0.5
I
41
I t
U '
0.0
0.4
0.8
0.0
0.2
0.4
Figure 9. Reactive flux correlation functions for three different models of cyclohexane. The solid lines are for the model with potential energy (2.8). The dotted lines are obtained with rigid C C bonds (equations of motion are integrated with these holonomic constraints). The dashed lines are for when the soft C-C bond force constant, k,,, is used in place of k. All other aspects of the models remain the same.
calculations" of the pressure dependence for kTsT, we estimate that the experimental turnover for cyclohexane occurs at solvent viscosities around 100 cP. The trajectory results give a turnover a t CY = l O I 3 s-l. Therefore, the proportionality constant c = a / ~ is estimated to be 10" s-l cP-' as shown in (2.13). Next, recall that our earlier calculations" for the solvent liquid CS2 indicate that the solvation contribution to the free energy of activation is a linear function of pressure. Hence, for the three solvents considered in experiments (carbon disulfide, acetone, and methylcyclohexane), the pressure-dependent transition state theory rate constant is estimated as
~TST(P) = ~ T S T @ ) ~xP[-B(P - PI1
(4.5)
where p is the pressure at the state of median pressure for the given solvent, and B (essentially the activation volume per k B T ) is assumed to be the same for each solvent and set equal to the value computed by usll for CS2 solvent:
B = 0.025 kbar-'
(4.6)
Therefore, in placing the experimental points on Figure 8, for each of the three solvents one parameter, krsT@), is adjusted to agree with the trajectory results at the median points. The BGK simulation results are seen to agree qualitatively with those of experiment. Therefore, we can indeed conclude that the experimentally observed increase in rate with increasing solvent pressure is a manifestation of non-RRKM behavior. In view of this conclusion, it seems worthwhile to examine the isomerization with slightly different models of the intramolecular modes. The results for some of the alterations we have considered are shown in Figure 9. The a = 0 reactive flux correlation function shows a significant sensitivity to changing the C-C bonds from rigid to stiff to soft vibrating bonds. In the inertial regime, we have found that the effect of this sensitivity on k(t) for finite a can be predicted by using ko(t) and the generalized strong collision approximation. Beyond Kramers turnover, there are no significant differences in the k(t)'s with different C-C bonds. Incidentally, the structures observed in k ( t ) for a below the Kramers turnover (as illustrated in the two left panels of Figure 9) are real and not artifacts of poor sampling statistics. It is interesting to note the different effects the softening of C-C vibrations has on the forward and backward components to ko(t). These are shown in Figure 10. For the forward case, k&+)(t),the factor of 3 decrease in spring constants removes nearly entirely the short-time non-RRKM component to the reactive flux. However, this degree of softening is not sufficient to remove the
0.0
-0.2 0.0
0.4
0.8
Figure 10. Forward (upper) and backward (lower) contributions to k,(t) with different choices of C-C spring constant. The solid lines are for the model with potential energy (2.8). The dashed lines are for when the soft C-C bond force constant, kWft,is used in place of k .
short-time relaxation in k&-)(t). Thus, even in this case, the net non-RRKM behavior contributes 20-30% leading (as discussed earlier) to a turnover near CY = 4 ps-]. This frequency is still an order of magnitude greater than would be found if the non-RRKM were entirely absent from ko(t). Having examined the general trends, it may now be worthwhile to perform a full molecular dynamics simulation of cyclohexane solvated among several hundred CS2molecules. The success of a strong collision approximation in interpreting our stochastic molecular dynamics results suggests that some form of analytic theory should provide an explanation of the full molecular dynamics and also a quantitative explanation of the experiments. A simple and accurate theoretical model will be welcome as there are at least two important experimental implications of our findings. The first concerns nonstatistical specificity in the progress of a unimolecular reaction. While there are many to choose from, one interesting example is found in Carpenter's observation26that inertial nondiffusive trajectories will significantly alter the stereochemical outcome of unimolecular processes from what would be anticipated from transition-state theory. Our calculations show that inertial rather than diffusive barrier crossings are in fact the type of activated dynamics occurring with the cyclohexane mblecule. Second, as we pointed out several years ago: dynamical effects can dominate the pressure dependence of rate constants. This behavior will occur in the inertial Lindemann regime. (The corresponding dynamical effects contemplated by T r ~ for e the ~ ~ diffusive regime are possibly important, too, but probably smaller.) Hence, mechanistic interpretation of measured activation volumes-especially for those that are larger than the partial molar volume changes for completed reactionsZ8-must be reassessed in light of likely and possibly significant violations of transitionstate theory. Acknowledgment. In carrying out this research, we have benefited from conversations and advice from Albert Nichols. Several conversations with Bruce Berne and Michael Borkoveck were most enlightening, too. The work was supported financially by grants from the NSF and NIH. Registry No. Cyclohexane, 110-82-7. (26) Carpenter, B. K. J . Am. Chem. SOC.1985, 107, 5730. (27) Troe, J. High Pressure Chemistry; Kelm, H., Ed.; Reidel: Dordecht, Netherlands, 1978. (28) Asano, T.; le Noble, W. J. Chem. Reu. 1978, 78, 407.
