J. Phys. Chem. 1982, 86, 2218-2224
2218
Dynamics of Reactions in Polar Solvents. Semiclassical Trajectory Studies of Electron-Transfer and Proton-Transfer Reactions Arleh Warshel Departmnt of Chemistry, University of Southern California,Los Angeles, California 90007 (Received:July 23, 1981)
A semiclassical trajectory approach for studies of the dynamics of reactions in polar solvent is developed. This approach focusses on the effect of fluctuations of the solvent dipoles on the rate of crossing between the solute resonance forms. The calculations of the surface crossing probability are based on a time-dependent quantum mechanical treatment which uses the classical time dependence of the energy gap (evaluatedalong the trajectory path of the solute-solvent system). The practical potential of the method is demonstrated in preliminary simulations of electron-transfer and proton-transfer reactions. The conceptual advantage of the method is demonstrated in exploring the molecular meaning of the activation free energy of electron-transfer reactions, and in a new derivation of the rate constant of electron-transfer reactions in the intermediate temperature range where the solvent motion is treated classically and the tunneling between the solute vibronic states is treated quantum mechanically. The applicability of the method for studies of bond-breakingreactions in polar solvents is demonstrated by calculating the rate of a proton-transfer reaction in the strong coupling limit. In this case the rate constant is expressed as a product of a factor that depends on the fluctuations of the solvent dipoles and a fador that depends on the energies of the solute resonance forms. This type of treatment is expected to be particularly useful in treating the dynamics of enzymatic reaction.
I. Introduction The dynamics of chemical processes in the condensed phase has long been a central point of interest for both experimental and theoretical ~tudies.l-~ Increased availability of high-speed computers has rendered microscopic simulation of such processes possible. Significant effort has been dedicated to exploring the dynamical effects of solvent4 (or protein) viscosity5 on the rate of activated events such as ground-state isomerizations. However, no attempt has been made to simulate the dynamics of real chemical reactions where the most important effect of the solvent is due to its polarity. In case of polar solvents, the fluctuations of the solvent dipoles are likely to play a much more important role than the solvent viscosity in affecting and controlling the reaction rate. The effect is expected to be of crucial importance in enzymatic reactions. The study of a reaction in a polar solvent requires a simple yet reliable potential surface that allows efficient simulation of the dynamics of such a reaction. Potential surfaces suitable for the study of reactions in solution and in enzymes have been obtained by the empirical valence bond (EVB) approach! In this paper we extend the EVB approach to semiclassical studies of the dynamics of reactions in polar solvents. The semiclassical treatment is based on evaluating the crossing probability between different electronic states quantum mechanically, as a function of the time-dependent energy gap between these states. The energy gap is evaluated classically along the multidimensional path of the trajectories of the solutesolvent system. This approach is demonstrated for two key processes: electron-transfer and proton-transfer reactions. (1)H. A. Kramers, Physica, 7 , 284 (1940). (2)S. Glasstone, K. J. Laidler, and H. Eyring, "The Theory of Rate Processes", McGraw-Hill, New York, 1941. (3)E. Eyring, J. Walter, and G. E. Kimball, "Quantum Chemistry", Wiley, New York, 1944. (4)B. J. Berne, J. L. Skinner, and P. G. Wolynes, J.Chem. Phys., 73, 4314 (1980). (5)J. A. McCammon and M. Karplus, Proc. Natl. Acad. Sci. U.S.A., 76,3585 (1979). (6)A. Warshel and R. M. Weiss, J.Am. Chem. Soc., 102,6218 (1980). 0022-365418212086-2218$01.25/0
Section I1 describes the evaluation of the EVB potential surface for reactions in polar solutions and derives the semiclassical trajectory approach in the two limiting cases of weak and strong coupling. This section also derives a rate expression that includes quantum mechanical tunneling through vibronic states of the solute. Section I11 presents preliminary applications of the method; this includes the following: (i) A trajectory study of the dynamics of an electron-transfer reaction is performed. (ii) The use of the semiclassical formalism in a derivation of the rate equation for electron-transfer reactions and in an analysis of the microscopic meaning of the activation free energy is explored. (iii) A preliminary study of a light-induced electron-transfer reaction is undertaken. This study demonstrates the potential use of the method in cases where the quantum yield can be determined by the reorientation time of the solvent dipoles. (iv) The application of the method in the strong coupling limit is demonstrated in a preliminary study of a proton-transfer reaction. This study simulates the dynamics of a real bond-breaking reaction in a polar solvent. It is found that the activation free energy of the reaction can be approximated by the s u m of the free energy needed to bring the solvent to its transition-state configuration, r f , and the activation energy of the solute at 6. This type of treatment is expected to be very useful in comparing the dynamics of enzymatic reactions to the corresponding reactions in aqueous solution. 11. Theoretical Approach 1. Potential Surfaces. Potential surfaces for dynamical studies of reactions in solutions are extremely hard to obtain reliably by standard quantum mechanical approaches. At present the most practical approach seems to be the empirical valence bond approach?,' This method represents the reacting system as a superposition of ionic and covalent resonance forms. The ffect of the solvent is introduced by adding the interaction with the solvent to the energies of the ionic resonance form. For example, for the bond-breaking reaction (7)A. Warshel, Biochemistry, 20,3167 (1981).
@ 1982 American Chemical Society
The Journal of Physical Chemistry, Vol. 86, No. 12, 1982 2219
Dynamics of Reactlons in Polar Solvents
AB
-
A+ + B-
(1)
--I
I
we consider the two resonance forms 41
= (A-B)
I#J2 = (A+ B-)
(2)
The effective Hamiltonian obtained for these is given (in a matrix representation) by W
z
/
Ibd2
where g and s denote gas phase and solution, respectively. R and r denote the solute and solvent coordinates and H g i j is the empirical gas-phase valence bond matrix element for the two states (see ref 6 for more details). The effect of the induced and permanent dipoles of the solvent is introduced through their interaction, V2)-,(R,r), with the ionic resonance forms of the solute V%R,r) = CQ,‘2’qj/lRi- rjl ij
+ CQjQj./lrj - rj.1 j>j’
vvdw
- f/zC7Ek2 k
(4)
where Qi are the solute charges in the ionic resonance form, qj are the solvent point charges, and Vvdw is the solutesolvent and solvent-solvent van der Waals interactions. y is the molecular polarizability of the solvent molecules and Ekis the electrostatic field at the center of the kth solvent molecule from the solute and the other solvent molecules. In the present work the solvent molecules are taken as water molecules with q of +0.25 and -0.5 for hydrogen and oxygen, respectively, and a y value of 1.4 A-3. The water-water van der Waals interactions are evaluated by use of a 6-12 potential function with eii* = 0.018 and 0.0519 kcal/mol and riio = 2.00 and 3.66 for H and 0, respectively. To save computer time we do not include explicitly in Ek contributions from the field due to the induced dipoles of the other solvent dipoles. These contributions are included implicitly by using an effective average dielectric constant of 1.36 = Ek/l.36) which was obtained previous19 in a self-consistent evaluation of the interaction between the solvent-induced dipoles. Our semiempirical approximation neglects electronic overlap (and therefore charge transfer) between the solvent and solute atoms. The off-diagonal matrix elements (H12) should include contributions of the solvent-induced dipoles which can be evaluated by the exciton formalism.27 These contributions are neglected in this work. 2. Semiclassical Trajectory Approach. The EVB effective Hamiltonian allows one to explore the dynamics of reactions in solutions. This can be done by the semiclassical trajectory approach*12 which was applied previously to photoisomerization reactions of medium-size m o l e c ~ l e s . ~We J ~ consider two limiting cases. a. The Weak Coupling Limit. The time-dependent wave function of the system can be written as
A
k12 = lim -lb2(At)l2 At-7 At
N
A lim At-7 At
= i
2
Wt) = Cbi(tMi(r,R) exp(-(i/h)s’ei[r(t?,R(t?l d t j (5)
where the 4i and ei are the diabatic valence bond wave functions and energies (i.e., H(r,R)dl = el(r,R)41 + (8) A. Warshel and M. Levitt, J. Mol. BioE., 103,227 (1976). (9)A. Warshel and M. Karplus, Chem. Phys. Lett., 32, 11 (1975). (10)A. Warshel, Nature (London),260,679 (1976). (11)J. C. Tully and R. M. Preston, J. Chem. Phys., 55, 562 (1971). (12)W. H.Miller and T. F. George, J. Chem. Phys., 56,5637(1972).
where At is the trajectory time, ti’ = (ti,l - ti)/2, and 7 is the time after which k12 reaches a constant value. Our study indicated that 16b2I2of eq 7 are given to a very good approximation (as long as HlZ2< Ih(&- &,)I by the Landau-Zener expression,20p21rLz.That is 16b612 N rLz(i) = 1 - e ~ p ( - H ~ ~ ~ / l h-( el’)il) e; (8) where Ad is the gradient of the Ae hypersurface at Ae = 0 and i is the velocity along this gradient. Equation 8 can be verified by an analytical continuation of eq 6 expanding
2220
The Journal of physical Chemistty, Vol. 86, No. 12, 1982
the JAthE around At(ti) = 0 (see ref 12 for a related treatment in the adiabatic approximation). This gives lSb',I2 = H122/1h(Ak)lwhich is identical with the expansion of eq 8 for H122< Ih(Ak)(where Ak = Ad?. Substituting eq 8 into eq 7 we obtain A k,, = lim -CrLz(?) At-r At i This type of expression will be used in section I11 in analyzing the multidimensional expression for electrontransfer reactions. b. The Strong Coupling Case. When H12is large the adiabatic basis set obtained by diagonalizing the Hamiltonian (i.e., H(r,R) $I(r,R) = EI(r,R) $I(r,R)) provides a better representation for the semiclassical studies. The corresponding time-dependent wave function is *(t) =
Substitution in the time-dependent Schrodinger equation an =
*I/dt)aI(t) exp((i/h)~'hE[r(t),R(t)l d t l (11) where hE = EII- EI. As in the diabatic case, the probability of crossing from state I to state I1 is determined by integrating aII along the path of the trajectory. 3. Quantum Mechanical Corrections. The evaluation of the rates of electron- and proton-transfer reactions may require significant tunneling corrections. These corrections can be obtained by separating the fast-moving, high-frequency solute modes (e.g., the intramolecular stretching vibrations) from the rest of the system, freezing these modes during the trajectory calculation, and treating them later quantum mechanically. In this case eq 6 can be extended to the treatment of transitions between vibrational states of different electronic states. For a typical case of a single effective vibrational quantum mode (e.g., the proton stretch in proton-transfer or C=C stretch in electron-transfer reactions) we obtain (*11(d
t
i)lm,2mt = -(i/h)H12Cmmte x p ( - ( i / h ) l (el
+ mhu, -
t2 -
mhv,) d t i (12) where Cmm,is the vibrational Franck-Condon factor for transition from the mth vibrational state of tl to the mth vibrational state of t2 and v, is the effective vibrational frequency of the given mode. This type of expression (which was used before in quantum mechanical studies of a radiationless transition^'^) is very useful for semiclassical studies since the integration of the energy gap term along the trajectory gives automatically the effective density of states (see section 111). 111. Results and Discussion 1. Outer-Shell Electron Transfer in Polar Solvents. a.
Semiclassical Treatment. Outer-shell electron-transfer reactions provide the simplest test case for our approach. In considering this type of reaction for a fixed and large distance between donor (D) and the acceptor (A) we can take a fixed H12(R)and consider two diabatic wave functions
41 = (D- A)
42 = (D A-)
(13)
(13)A. D.Brailaford and T. Y. Chang, J. Chem. Phys., 53,3108(1970).
Warshei
\
a!
/
/
i Reaction koordinate;
Figure 2. A schematic presentation of the potential surfaces t, and e2 for an electron-transfer reaction between a donor (D)and an acceptor (A). The reaction coordinate involves changes in the solvent coordinates (r) which are represented here by the reorientation of the dipoles in the lower part of the figure. In r i the solvent dipoles give maximum stabilization to D- and in r: they give maximum stabiliation to A-.
The diabatic energies (el and t2) for such a system are described schematically in Figure 2 as a function of the solvent coordinates. When the reaction coordinate can be reduced to a one-dimensional coordinate it is easy to derive (see below) a classical (high-temperature) expression for the rate constant of the electron-transfer reaction'"l6 k12 = (H122/h)(4TCYkbT)-1/2 eXp(-A&/kbT) (14) where k b is the Boltzmann constant and At* is given by q(d)- tl(ro)(see Figure 2). When t1 and t2 are harmonic functions of equal curvature one obtains At$ = (Ato + C Y ) ~ / ~ C Y (15) where Ac0 = t2(rz0) - El(rt)and CY is the solvent reorganization energy defined in Figure 2. Unfortunately eq 14 cannot be used in studying multidimensional cases since the use of energy differences other than free energy differences can be shown to violate the condition of microscopic reversibility. It is not entirely clear, however, whether one can simply replace AtO and CY by the corresponding free energy terms. Furthermore, the dynamics of electron-transfer reactions have not been explored on a microscopic level beyond the harmonic approximation. These types of problems can be explored with our semiclassical approach. As a test system we take two spheres of radii 2 A (representing D and A) held 4 A apart and embedded in 40 water molecules. The water molecules are surrounded by a surface of 15 fixed water molecules that simulate the bulk solvent in the spirit of the SCSSD model of ref 19 We take Ato = -30 kcal/mol and H12= 0.01 kcal/mol. Figure 1 describes the time evolution of the surfaces tl and e2 after initial equilibration at 300 K. The probability that the system will cross from surface 1 to 2 (the probability of electron transfer) is evaluated by eq 6 following the procedure described in section I1 where the rate is evaluated by eq 7. As much as the numerical evaluation of a specific rate constant is concerned, we consider the results of the present work only as a preliminary demonstration of the potential of a new approach. However, even without the numerical results, it seems that our model provides an (14)R. A. Marcus, Annu. Rev. Phys. Chem., 15, 155 (1964). (15)R. A. Marcus, J . Chem. Phys., 43,679 (1965). (16)V. G.Levich, Adu. Electrochem. Electrochem. Eng., 4,249 (1966). (17)J. J. Hopfield, h o c . Natl. Acad. Sci. U.S.A., 71, 3640 (1974). (18)A. Warshel, h o c . Natl. Acad. Sci. U.S.A., 77,3106 (1980). (19)A. Warshel, J. Phys. Chem., 83, 1640 (1979).
The Journal of Physlcal Chemistry, Vol. 86, No. 12, 1982 2221
Dynamics of Reactions in Polar Solvents
instructive way for analyzing (and thinking about) the multidimensional analogue of eq 14. This point is demonstrated below. We start by deriving eq 14 for the one-dimensional case (where r is represented by a single coordinate x ) . This is done by noting that the rate can be evaluated by eq 9 and reexpressing this equation as
k12 =
S t l mP(x*(t?)P(2)
dx dt’
(16)
where we assume formally that the rate is evaluated over a time much larger than the T of eq 9 and P(xf) dt’is the probability of reaching the crossing point xf (A&*) = 0) a t a time between t’ and t’ + dt’ and P(x) d i is the probability that the trajectory will have a velocity between 1 and 1 + &. Using P(x’(t 9)dt’ = P(x’)) dx which gives P(x’(t?) = P ( x * ( x ) ) iand the relation d / d x J x Y ( x ? dx’= Y(x)we obtain
k12 =
S P(x’(x))kP(k)rL&)dx
(17)
This expression is essentially the integral obtained by Levich16 which gives using the normalized expressions for P(x*)and P ( 3 ) the rate expression k12
= eXp(-(cl(d) - tl(.1°))/kbT1(H122/h)(4?r(ykbT)-1/2
(18)
where x o is the minimum energy coordinate of el. In order to treat the multidimensional case we divide the space of the system into classes, Si,which satisfy the equation, tl(Si)- t2(S’) = C’. The locus of points x i that minimize tl(Si)is the reaction coordinate. Thus, for example, d is the least energy point among the points S* that satisfy the relation tl(S*)- c2(S*)= 0. The relative probability of reaching xf (which is given in the one-dimensional case by exp(-(t,(xf) - tl(Xo))/kbT))is given in the multidimensional case by
P(xf)/P(xO)= exp(-ti(S*)/h,T) d S ’ / S exp(-ei(SO)/kbT) dSo (19) This can be rewritten as
Pbf)/P(XO) = exP(-Ag(x*)/kbT)
(20)
and the overall rate is given by k12
= (H12/h)(47&bT)-1/2 eXp(-&(d)/kbq
(21)
The key problem, of course, is the relationship between our Ag* and the thermodynamic parameters AGOand a. The difficult task of examining this relationship can be accomplished, in principle, by using trajectory calculations. That is, the trajectory can be used to obtain relative probabilities P(xi)/P(xi)by d / d where niis the number of times the value of At falls within the region Ci to Ci+ dC. The a’s give the relative free energy, Ag(xi) along the reaction coordinate by gI(X4) - gl(Xp) = -kbT In
(ff$/$)
(22)
where the subscript 1 indicates trajectories over el, and fit is the largest It among all the iii. The Ag of eq 22 can be used to evaluate the multidimensional analogue of eq 15. A preliminary study along this line is reported in the Appendix. This study indicates that Agf can be expressed as Agf = (AGO+ ~ ? ) ~ / 4 n
b. Tunneling Corrections. The current treatment represents the donor and acceptor as conducting spheres; however, it is quite simple to extend our approach to realistic systems evaluating the solute potential surfaces and Franck-Condon factors by the QCFF/PI method.23 In these cases one can represent the high-frequency modes with large Franck-Condon factors by one effective mode v,18 so that the rate can be evaluated by eq 7 and 12 giving
(23)
where AGO= G2 - G1 is the equilibrium free energy difference between state 2 and 1, while a is given by a = At(xY) + AGo.
A At i
k12 = CC lim -C16bim,2m12 exp(-mhvr/kbT) = mm’
CC(H12/ h )Cmmp2 exp(-mhvr/kbT) x m m‘
lim -CISt”exp(-(i/h)StAtmm. A d t ) dtI2 (24) At i t;
+
where Atmm,= (el hvrm - t2 - hvrm?. The direct evaluation of the tunneling correction is extremely simple in the above formulation by integrating eq 24 along the trajectory. Equation 24 provides a very simple way to obtain the rate expression in an approximate closed form. That is, any 16b1m,h,12in eq 24 is formally identical with the expression used in eq 7 where Hlz is replaced by H12Cmmt.andAt by Atmmt.Thus we can evaluate each of the terms in the s u m of eq 24 using the expression obtained in eq 21, where Agk,, = -kbT In (&,/&,) and iifnmlis the number of times that Acmm,becomes zero. This gives m m‘
exP(-&tm!/ kbT) e x p b h v r / kb T ) (25) By the same argument that Agf is given, at least approximately, by (AG, + ~ x ) ~ / 4we a obtain Agt,, = (AGO- hvr(m - m?
+
C Y ) ~ / ~ C= Y
-kbT In
(fikmr/&m,) (26)
with the analytical expressions for the Franck-Condon factors2s we obtain the expressions used in ref 18 and 24 for cases with lAGol > 1y (where the quantum corrections are crucialz4). 2. Light-Induced Electron Transfer. The present semiclassical approach is particularly useful for studies of light-induced electron-transfer reactions. These type of ultrafast processes can be simulated in a reasonable computer time. The detailed simulation can provide insight into interesting cases where the reorientation time of the solvent dipoles can determine the reaction yield (see below). To demonstrate the potential of our approach we study a model reaction of light-induced electron transfer in a system with a donor, D, composed of two atoms, ClC2(that simulates a double bond), and an acceptor sphere, A, held 4 8, from C2. The charge distribution of the relevant electronic states are indicated in the caption of Figure 3. The system is first equilibrated and then excited to the partially polar excited state (T*). The excited system is propagated on the a* state and the crossing probability to the charge-transfer state (CT) is evaluated semiclassi(20)L. Landau, Phys. Z. Sowjetunion, 2 , 46 (1932). (21)C. Zener, Proc. R. SOC.London, Ser. A, 137,696(1932). (22)J. P. Valleau and G. M. Torrie in ”Modern Theoretical Chemistry”, B. Beme, Ed., Vol. 5,Plenum Press, New York, 1977. (23)A. Warshel in “Modern Theoretical Chemistry”, G. Segal, Ed., Vol. 7,Plenum Press, New York, 1977. (24)A. Warshel and D. Schlosser,Proc. Natl. Acad. Sci. U.S.A.,78, 5564 (1981). (25)A. Warshel, Photochem. Photobiol., 30, 285 (1979). (26)A. Warshel and P. Dauber, J . Chem. Phys., 66, 5477 (1977). (27)A. Warshel, to be published.
The Journal of phvsical Chemistry, Vol. 86,No. 12, 1982
2222
. ........... . : ............................... .. .. . . . "(.,, D f A -
...
....
... "..
j I
\--
I
........... ....
---'-.
.-
Warshel
Solvent Coord I na tes
, 00
I
1
02
01 TIME (picoseconds1
Figure 3. Simulation of a llght-lnduced electron-transfer reaction in a system of a donor, D (represented by two carbon atoms), and a spherical acceptor, A. The trajectory of the solvent molecules Is propagated on the sound state (DA) potential surface and equUibrated at 300 K. After equilibration, at a time t = 0, the DA system is excited to the D"A state and the probability of the D'A D+A- transition is evaluated, integrating eq 6 along the trajectory. The assumed charge distributions are (-0.3, 0.3, O.O), (0.3, -0.3, O.O), and (0.7, 0.3, -1.0) for DA, D'A, and D'A-.
Flgure 4. A schematic description of potential surfaces for protontransfer reactions. The figure presents the diabatic potential surfaces t, and tP and the grwnd-state potential, ,EI, obtained from their mixing. r and R are the solvent and solute coordinates, respectively.
-
cally. For this model system we obtain the quantum yield, Y , for forming the charge-transfer state by
Y = ~+-.CT/(~~*.-.CT + k-)
where k- is the rate constant for the back-crossing to the ground state by radiationless transitions and fluorescence. The quantum yield of the model system in eq 27 was evaluated as a demonstration for the case with H12= 0.07 kcal/mol and k- = 10" s-l. With this Hlz, a very large fraction (-0.2) of the molecules cross during the first intersection of the ?r* and the CT states. The intersection time is determined by the solvent reorientation time (which is determined by the dielectric relaxation time). Thus, while the present calculation gave Y N 1,a solvent with slower reorientation time can give a much smaller quantum yield (e.g., for a reorientation time of 10-l' s the quantum yield will be -0.3). This preliminary calculation indicates that there should be a class of interesting cases where the dielectric relaxation of the solvent will determine the quantum yield of light-induced electron-transfer reactions. 3. Proton- Transfer Reactions. Proton-transfer reactions play a major role in biological systems including the action of hydrolytic enzymes7and light-induced proton pumps.25 For a typical proton-transfer reaction AH
+ B -AA-+
BH+
(28)
The most important resonance forms are
41 = (A-H B) 42 = (A- H-B') (29) The EVB matrix elements for these resonance forms can be determined by the EVB approach?' where (in contrast to ref 6) we do not include explicitly the third resonance and 42 as effective form 43 = (A-H+B) and consider states that already include the effect of d3 (in this representation 41is slightly polar). Typical potential surfaces t1 and t2 and the adiabatic potential surfaces EI and EII obtained from their mixing are described in Figure 4. The most important parameters in determining the shape of the surfaces are the hEoand the barrier height hEt. These two parameters are given approximately as7 AEo
uL
N
N
- Hi2
nl
15.25
15.50 TIME (picoseconds)
15 75
15 25
15 50 TIME (picoseconds)
1575
(30)
1600
(C)-
15.00
1600
Figure 5. A semiclassical trajectory calculation of the rate of a proton-transfer reaction in the strong coupling limit. The calculations are performed on a system of a formic acld molecule (A) and a water molecule (B) held 3 A apart in aqueous solution. Part (a) shows the fluctuations of the ground-state potential, ,E,(R) (for different values of the proton coordinate, R), along the pathway of the trajectory of the solvent molecules. Part (b) shows the variation of the proton coordinate as a function of the trajectory time. The figure indicates that the proton jumps over the barrier only when the solvent reaches the configurations r* that stabilize the charge distribution of the (A-HB+) resonance form.
a,,
where is approximated by the free energy of proton transfer at 300 K. MBH+is a Morse potential for the covalent BH+ bond and DBH+ is the corresponding dissociation energy.6 As seen from eq 30 the key factors are the pK,'s, which are strongly dependent on the solvation energy of A- and BH', and the distance between A and B which determine the barrier height; for short distance aEt is given approximately by hEowhile for large distance hEt increases as fast as the Morse potential increases. Most proton-transfer reactions occur at short distances between A and B, since is very large at large distances. At short distances (RA-B I4 A) H12is large and the reaction should be treated in the adiabatic representation by using eq 11. Our calculations indicated that for large H12the probability of crossing to state I1 is small and the system tends to stay in the ground state (state I). Thus it is possible to evaluate the reaction rate by propagating the trajectories of the solute solvent system on the adiabatic potential surface and counting the number of crossing n1-2 of R from R1 to R2 per unit time (see Figure 5). The rate will be given by 1
lim -Cnl,p t-+r t this can be expressed formally as kpT =
1.33[pKa(A)- pK,(BH+)]
+ (MBH+(RBH+) + &H+)
00
(27)
(31)
Dynamics of Reactions in Polar Solvents
kw = P(R*)VK
The Journal of Physical Chemistry, Vol. 86, No. 12, 1982 2223
(32)
where B is the average oscillation frequency of R, K is the average probability the trajectory that reaches R* continues to R!, and P(R*)is the relative probability of reaching R*. With P(R*)= exp(-Ag*/kbT) we obtain basically the familiar expression from transition-state theory. Unfortunately, such a brute force approach for evaluating Agt can be very time consuming for high barriers (although special procedures that involve propagation of the trajectory downhill from R* can simplify the problem). Furthermore, such an approach tends to mask the crucial role of the solvent in controlling the reaction dynamics. In order to gain more insight into the role of the solvent we take an alternative approach and treat the solute and solvent motions separately. Such a separation is justified by the fact that the solute vibrations are faster than the solvent fluctuations (see Figure 5) and that the energy barrier (EI(R*) - E@)) is almost the same for any solvent configuration that gives the same El(@) - (EI(Ry)), regardless of the particular value of the solute coordinate R. Treating R and r as uncoupled coordinates allows one to run a trajectory over E to evaluate the potential difference EI(R*)- EI(Rl) i forand each solvent configuration along the trajectory. The overall probability to reach the transition state is given in such a treatment by
I'
Figure 6. A preliminary examination of the validity of eq 23. The calculations were done for the system described in Fgue 1. The points near the minimum of the Ag, curve were evaluated with eq 22. The points on the Ag2 curve (which could not be evaluated by eq 22 with any reasonable trajectory time) were calculated with eq 47 and 48. The estimated error bars are indicated in the figure. They are relatively large since the calculations present only a preliminary study.
(e)
Appendix This section examines the validity of eq 23 by the semiclassical trajectory approach. Following section 111.1we divide the coordinate space to classes of coordinates Sithat satisfies the relation
P(R*) exp{-k(R*) - g(RP))/kJ'l = (exp{-(EI(R*,r)- EI(RP,r))/kbq)R,O (33)
e2(Si) - €'(Si) = A&
where the identity proposed by Valleau22is used to relate the Ag and AE terms, and the average ( ) R 1 ~is evaluated over the trajectory that moves on EI(Ry). In order to separate Ag* into solute and solvent contributions we use the fact that pI(R*)= (p(2)+ p('))/2 and pI(Ry) = p(') (where pI and p(') are the adiabatic and diabatic charge distributions of the solute). This gives EI(R*,r)- E,(RP,r) = AEAR*(r)),o~u,+ [V,,dpf(R*,r) - Vsol(~I(RP,r)l(34) where Vwl is the calculated solvation energy (obtained from eq 4)for the given p, and AE(R*),olu,is given by (&R*)Hiz(R*))- (d(Ry)- Hi,2/(tq(Ry) - 4(R!))). As long as H1&i:) is significantly smaller than E$(R!$ - er(Ry), the solute energy term AE,(R*) depends only weakly on the solvent coordinates and eq 33 can be approximated by kw
=~
=
~ ~ ~ ~ * ~ s o l ~ ~ ~ * ~ s o l l J ~ FK
exP{-&bl/kbT! eXPk&blute/kbT) (35)
where Agio]and Ag~ol, are the free energies obtained by 8 / k b q )R10 and (exp~-AEI(R*)solute/ (ex&( v w l b l ) - Vwl(p1) kbq)Rlo,respectively. Thus the probability of reaching Rt is expressed as the product of the probability that the solvent will reach the configurations r*(that stabilizes the charge distribution p ) times the probability that the solute will "jump" from R1to R* at the given r*. This point is best understood from Figure 5 which shows that the solute passes through R* only after the solvent reaches a configuration that stabilizes the (A-H-B+) resonance form. The present work indicates that eq 32 and 35 give similar order of magnitudes for kpT (i.e., log (kpT) = 11 f 1) but more careful studies are clearly needed (studies along this line are now underway in our laboratory). However, we already have a clear indication that the fluctuations of the solvent dipoles play in a key role in determining the overall activation barrier for proton transfer reactions. We believe that the present approach will be particularly useful in comparing the rates of related reactions in solutions and in enzymes.
(36)
where ell is a constant. The partition functions of the system are defined (for a given potential, e l ) as
where /3 = l/kbT and Qris the total partition function. The probability that the system will be between xi and xi dx is given by ( q ( x i ) / Q l )dx and the corresponding free energy is given by
+
k!r(x') - Gr = -kbT In ( q i ( x i ) / Q J
(38)
This free energy can be evaluated by the trajectory calculations with the relation gi(xi) - Gi = -kbT In [fii(x')/Ciii(x')] J
(39)
where fil(xi)is the number of times the value of Atzl falls between A& and Ati + b along the trajectory path. The gl's and g,s' (evaluated by eq 39 for trajectories that move over t1 and t2, respectively) are related to each other by gib') - gz(x') =
kibi)- GI)- (92(xj) - G2) + (GI - G2) (40)
and by the definition of x i we have
J
g2(xi)- gl(xi) =
A&
(41)
The above definitions allow one to examine the validity of eq 23 by evaluating g l ( x j ) - gl(xy)) and (g2(xj)- gl(xy)) (using eq 39 for different Atl2) as functions of ((G2- GI) a - A t ~ 1 ) / ( 2 ~ ) (see 1 / 2 Figure 6) where gl(xy) is the minimum of gl(x) and a is given by
+
= AQ~(X?)- (G2 - GI)
(42)
If the Ag's give quadratic functions with equal curvature then eq 23 is satisfied. However, direct evaluation of the Ag's might involve serious problems when Ag! = g , ( x f ) - gl(x!) is large and f i l ( d * ) / x i f i ( x iis) so small that only infinitely long trajectory can reach x*. This problem can be overcome by using the philosophy of the "umbrella
J. PhyS. Chem. 1982, 86, 2224-2231
2224
sampling" method.22 A reference state with an energy tl can be defined by (43) E l = V,l(P') + 4 where V,, denotes the solvation energy evaluated by eq 4 for a solute's charge distribution, pl, given by (44) P' = P' + 4 ( P 2 - P') where the parameter has a value between zero and one. The reference state is chosen so that trajectories over tl will pass the largest number of times through a given subspace Si (maximizing &(xi)). Thus, even if the probability to reach a given Si while moving on tl is very small we can use eq 41 and obtain (gl(d) - GI) = &'(xi) - GJ + (GI - GI) - Atfl (45) This gives g,(x') - gAx9 = (gi(x') - GJ - (gi(xP) - GI) + (GI - GI) - A& (46)
The free energy difference (GI - G,) is given by22
~XP~-P(GI - GI)) = (exp{-@(el-
(47)
where ( denotes an average over a trajectory that move on el. Our trajectory study indicates that &'(xi) - GI) and (g,(x?) - G,) are of similar order of magnitudes, thus we obtain the approximate expression gl(xi) - gl(xg)
E
GI - GI - A&
(48)
where p' is the solute charge distribution that maximizes the value of aI(xi)for a given xi. Equation 48 approximates &'(xi) by the difference in solvation energy of the charge distributions p' and p'. Since this free energy difference is a quadratic function of 01, at least in the dielectric continuum approximation, we conclude that eq 23 is a reasonable approximation. More numerical studies which involve direct calculations of eq 46 (Figure 6) are clearly needed for verification of eq 23.
Angular Momentum Decoupling Approximations. Current Status, Successes, and Difficulties D. J. Kauri. and D. E. Fltz Department of Chemnpby and Department of Mysics, Central Campus, Universlfy of Houston, Hwston, Texas 77004 (Received: July 27, 198 1)
The current status of angular momentum decoupling approximations is surveyed to identify those problems which are well treated and those which still present difficulties. Recent efforts to correct sudden approximations are briefly discussed.
I. Introduction During the past decade, there has been great progress in the development and testing of new approximations for treating atom-molecule collisions.' Some of the most successful and widely used approximations are the socalled "sudden" decoupling methods. These sudden approximations are based on the same type consideration employed in the Born-Oppenheimer separation of electronic and nuclear motions in the theory of molecular structure.2 Since the electrons move much faster than the nuclei, due to the wide disparity in the electronic and nuclear masses, one is able to solve for the electronic motion, in the field of the fixed nuclei, as a function of nuclear position. In the case of molecular collisions, it can sometimes be the case that the time scale for certain nuclear motions may be much slower than that of other nuclear motions, in particular, the time during which the potential is significant. Thus,one may solve the equations for the rapid collisional degree of freedom by holding the slower degrees of freedom fiied. The result is a scattering amplitude which depends parametrically on the slower (1) Reviews of a variety of methods can be found in "Atom-Molecule collieion Theory: A Guide for the Experimentalist",R. B.Bernetein, Ed., Plenum Press, New York, 1979. Also see "Dynamicsof Molecular Collisions",Parta A and B,W. H. Miller, Ed.,Plenum Press, New York, 1976. (2) M.Born and J. R. Oppenheimer, Ann. Phys., 84,457 (1927). 0022-3654/82/2086-222480 1.2510
degrees of freedom. Then one averages this amplitude appropriately over these slower degrees of freedom to obtain the physically meaningful scattering amplitude. Of course, the above is a highly simplified description of these methods (e.g., the averaging of the sudden amplitude can be viewed as a change of basis from the approximately diagonal coordinate representation basis to the internal state quantum number basis) and there is considerable subtlety in the exact details of how this procedure is carried The typical nuclear motions chosen as being slow (compared to the relative radial motion during the collision) are (a) the molecular rotation (yielding the energy sudden or ES approximation1t4),(b) the orbital rotation of the projectile relative to the molecular center of mass (yielding the centrifugal sudden or coupled states (CS) appr~ximation'.~), (c) a combination of (a) and (b) leading to the infinite order sudden (10s) approximation.lP6 (3) See the review and references therein, by D. J. Kouri in ref 1; also Kouri, and D. K. Hoffman, J. Chem. Phys., 74,2275 (1981); V.Khare, D. E. Fib, and D. J. Kouri, ibid., 73,2802 (1980), and references therein. (4) D. M. Chaw, Phys. Rev., 104, 838 (1956); S. I. Chu and A. Dalgarno, Roc. R. SOC.London, Ser. A , 342,191 (1975); D. Secrest, J . Chem. Phys., 62,710 (1975);V. Khare, ibid.,68,4631 (1978). See also references in the review by Kouri in ref 1. (5) P. McGuire and D. J. Kouri, J.Chem. Phys., 60, 2488 (1974);R. T.Pack, zbzd.,60,633 (1974). Again, see the review by Kouri in ref 1 and also ref 3. see V. Khare, D. J.
0 1982 American Chemical Society
