J . Phys. Chem. 1993,97, 327Q-3277
3270
Calculation of the Rate Constant for the Reaction C1Framework of the Continuum Medium Model
+ CH3Cl -,ClCH3 + C1- in the
M. V, Basilevsky,' G. E. Chudinov, and D. V. Napolov Karpov Institute of Physical Chemistry, ul. Obukha 10, 1030644 Moscow K64, Russia Received: August 31 1992; In Final Form: December 14, 1992 I
+
A complete evaluation of the reaction rate constant for the S Nreaction ~ CI- CH3Cl- ClCH, + C1- in water is performed within the continuum solvent approximation. The sequence of calculational steps involves the quantum-chemical (in terms of the self-consistent reaction field theory) evaluation of the two-dimensional free energy surface (FES)in the subspace spanned by both solute and medium coordinates. The subsequent steps are the separation of inertial and noninertial polarization components, the formulation of the dynamical equation and its kinetic treatment in terms of the multidimensional Kramers-Grote-Hynes theory. The main numerical result, Le., the value found for the transmission factor x , is compared with the results of alternative earlier calculations. The charge-transfer reactions are classified according to the shape of the FES cross section along the medium coordinate near the saddle point. For single-well-shaped profiles the x value is expected to be close to unity. This is the case for the above reaction. For double-well-shaped profiles the reaction proceeds via an electron-transfer mechanism and deviations from the standard transition state theory may be significant.
1. Introduction Recently, the investigations of kinetics and reactivity in chemical reactions with charge transfer proceeding in polar solvents have been developing intensively. According to the traditional approach' the solvent affects the height of the reaction barrier by modifying the solvation energies of initial and transition states. However, a modern view of this problem as applied to systems with a strong interaction between a reacting solute and the medium implies that the medium effect can by no means be reduced to a simple modification of the adiabatic potential controlling a chemical conversion.* Rather, the medium variables must be treated on equal grounds with the variables of the reacting solute because the corresponding nuclear motions may become decisive in a critical region of the free energy surface ( F a ) , that is to say, in the vicinity of its saddle point. Several problems cannot be consistently solved at all if an explicit treatment of medium coordinates is ignored. Among them a very important one is the evaluation of the rate constant preexponential factor. In the present paper we explore the FES for the simplest nucleophilic substitution (SN2) reaction in water
C1-
+ CH,C1
-
C1,CH
+ C1-
the activation free energy AG*. The evaluation of partition functions is needed to obtain A . The calculation of transmission factor x is a much more sophisticated piece of work which has been performed by themolecular dynamical (MD) treatment.I6.'8 We calculated the FES of reaction 1 in terms of the selfconsistent reaction field (SCRF) method with the NDDO approximation (the PM3 versionz2) applied at the quantumchemical stage of the calculation. The nontrivial peculiarity of the present treatment is an explicit consideration of the FES cross sections along thecollective mediumcoordinate in thecontext of a quantum-chemical calculation. At the next step we used stochastic dynamical theory to estimate the value of the transmission factor x . Therefore, the approach elaborated here comprises a consistent continuum counterpart of the abovementioned microscopic calculations.14-'8 The relative merits and disadvantagesof the present continuum model are then discussed. 2. The Gas-Phase Potential Energy Surface We calculated the potential surface for the gas-phase reaction 1 using the semiempirical PM3 procedure. The following characteristics of the symmetrical transition state (TS)
(1)
(H20)
in the configurational space including a medium coordinate simultaneously with the natural solute variable, namely, its gasphase reaction coordinate. During the last decade the model reaction 1 has been the most popular object for theoreticalstudies of solvation effects on reactivity.3-21 Among numerous works of this group, the famous microscopic calculations using simulation m o d e l ~ ' ~are - I ~especially important. The calculation of the free energy profile combining ab initio quantum-chemical techniques with the Monte Carlo m e t h ~ d lresulted ~ . ~ ~ in the barrier height AG* in surprisingly nice agreement with the experimental activation energy. On the basis of this achievement, the rate constant is calculated as
m)
K = xA exp( - AG *
where R is the gas constant and T is the temperature. The prefactor A, as described in terms of the transition state theory (TST), deals only with the solute degrees of freedom because the contribution of medium variables is already incorporated into
(3) were found: (1) the formation energy of the electrostatic prereaction complex, AUI = -9.7 kcal/mol; (2) the height of the potential barrier relative to the reactant level, A&* = mol; (3) the lengths of C-C1 and C-H bonds, rc-cl =0*8 2.189 and rC-H = 1.081 A; (4)the charges on atoms C11 and Clz and on the CH3 group, pel, = p c ~ ?= 0.75 and ~ C H ,= 0.5. These results are rather close to those obtained by the AM1 method.s In particular, the quantity A&* is underestimated by 2-3 kcal/ mol as compared to the experimental value and by 4-7 kcal/mol in a comparison with the results of the best ab initio calculations including electron correlation (see Table I). The error arises due to the deficient PM3 parametrization for the C1 atom. In the calculation based on AM1,S this defect has been eliminated by a recalibration of one-center parameters for CI and C atoms. This allowed, as shown in Table I, a satisfactory barrier height to be gained but did not improve such an important TS characteristic as the length of the C-Cl bond. It was not the purpose of the present work to perfectly calculate the absolute barrier height. That is why we kept the PM3
0022-3654/93/2097-3270~04.00~0 0 1993 American Chemical Society
kcal
Rate Constant for C1-
+ CHjCl
-
CICH,+ C1-
The Journal of Physical Chemistry, Vol. 97, No. 13, 1993 3271
TABLE I: Several Parameters of tbe TS (Expression 3) in the cas Pbase method AUO’,kcal/mol rc ~ 1 A, refs STO-631G* account of the electron correlation experiment AMI improved AM1 PM3
3.6 5-8
2.383
14, 23
2.155 2.189
5, 14,23-25 5 5 present work
parametrization scheme untouched. Instead of that, at the stage of the solvation study presented in the forthcoming sections, we added 3-5 kcal/mol to the activation barrier height appearing in calculations, which is in accord with the data of Table I.
3. The Calculational Scheme for Equilibrium Solvation The treatment of the energetics of reaction 1 in solution was performed in terms of the previously elaborated versionZ6of the SCRF method. The medium was described by a double-valued dielectric permittivity function, 1
( r e Vi)
to
(rev.,)
c = c (r) =
u(r) = x A q i 6 ( r - ri)
(8)
i
23
3 f 1.8 0.5 3.1 0.8
orders found at the quantum-chemical stage of the calculation. At the numerical level,
where 6(r - ri) is the surface delta function. According to eq 8 the total surface S is divided into small pieces Psi centered at points ri. Then the solution of eq 5 provides elementary surface distributions ai and their charges Aqi as follows from eq 8. The distribution in eq 8 generates the potential CP to be inserted into the quantum-chemical calculation at the next step, (9)
Its matrix elements according to the NDDO approach are
(4)
where r is a position vector, V, and V, are the volumes inside and outside the cavity cut in the medium bulk, where the solute, Le., the system C1CHjCl, was accommodated. Under the limitation that thesolutecharge is entirely contained in thecavity, the potential field of the medium polarization is reduced to that created by the charge density a(r) distributed over the cavity boundary S. The quantity a(r) obeys the following integral equation:
+
By denoting as Fo,, the Fockian matrix elements in the absence of the medium we come out with the resulting matrix elements of the total Fockian: The parametrization of integrals (eq 6) was borrowed from the PM3 scheme. Theelectrostatic part of the free (Gibbs) solvationenergy Ace, was calculated as a difference of two quantities: the formation free energy AG for the compound interacting system ”the solute + the medium polarization field” and the enthalpy of formation AIIOof the isolated solute:
AG,I = AG - A H 0
The chemical substrate (the solute) has the charge distribution p(r). Operator a/an(r) represents the differentiation of argument r along the external normal to S. The whole calculation followed the earlier elaborated prescription.26-28 The cavity was constructed as a composition of interlocking spheres drawn around each of the solute atoms, the radius of each sphere being accepted to be 1.2 times the corresponding van der Waals radius. The NDDO approximation was used in the quantum-chemical calculations. When performing the electrostatic part of a calculation, the potential cp(r(p) is needed, which is created by the charge distribution p at point r. It is composed of elementary potentials generated either by basic electron distributions xrxv pv ( x , and xu are AOs of the solute) or by nuclear unit charges pA = A (A is the index of atomic nuclei). These quantities are defined as
(12) The evaluation of AG involves the estimation of the self-energy Er of the polarization field:
AG = AGO + E,,
+ E,,,,,, + E,
(13) Here AGOis the intramolecular solute component of AG from which the solute-medium interactions are eliminated. Nevertheless, the wave function and the geometry of the solute must be calculated in the presence of the field CP. The quantities E,, and E,,,,, are the interaction energies of the surface charges Aqi with the electrons and the nuclei of the solute:
The sum of the terms of eq 14 is always negative. It is utilized to compute the field self-energy:
In the subsequent study only the H F approximation was used. The corresponding potential, with the NDDO approximation invoked, is
4. The Free Energy Profile in Water
By using the methodology described in section 3 we calculated the FES cross section along the gas-phase reaction coordinate x in water ( € 0 = 80). Near the TS, x is an antisymmetrical mode: x = (l/2)1’2(~c-c,,- rc-c,,)
where Z A are the nuclear core charges and P,, are the bond
(16)
The barrier height so found was AG* = 20.4 kcal/mol. Comparing
Basilevsky et al.
3212 The Journal of Physical Chemistry, Vol. 97, No. 13, 1993 TABLE II: Variations of Charges on Chlorine Atoms and on the Methyl Group along tbe Reaction Coordinate X
pel,'"' PCl, p ~ "1l PCl~ p~ 14 ,'"I PC ti 1
0.2828
0.2121
0.1414
0.0707
0.0
-0.92 -0.95 -0.49 -0.44
-0.89 -0.93 -0.56 -0.52 0.45 0.44
-0.86 -0.89 -0.62 -0.60 0.48 0.50
-0.81 -0.84 -0.69 -0.69 0.50 0.53
-0.75 -0.71 -0.75 -0.17 0.50 0.55
0.4 I
0.39
TABLE 111: Force Constants' at the Saddle Point for Three FESs (Eqs 17, 29, and 30) the FES in the total polarization space BO approximation SC approximation
a' AG
(I
formula
A,,
Aoo
A,o
17
-85.8
122.7
-124.0
29 30
-131.9 -133.1
78.4 76.8
-78.4 -76.8
axaq etc. In kcal/(mol A').
A,, = -;ax2
The charges p , and pl'"' a r e calculated in solution and in the gas phase, respectively.
it with the experimental activation energy AG*exp= 26.5 kcal/ molI4J0we conclude that an acceptable agreement can be gained only after adding the correction of the error entered at the quantum-chemical stage of the calculation (addition of 3-5 kcal/ mol). Based on this reasoning we consider as satisfactory the calculation of the solvation contribution to the barrier height. The TS geometry is practically the same as that in the gas phase. The values of the charges on chlorine atoms and on the CH3 group varying along x around the TS are shown in Table I1 together with the same charges found in the gas phase. The observed discrepancy between the two sets of charges is indeed insignificant, which manifests a weak polarizability of the solute system (reaction 1) and explains the absence of solvent-induced geometry changes. The gas-phase values of the charges are very similar to those found in ab initio ~ a l c u l a t i o n s . ~The ~ ~ 'present ~ calculation neglected the energetic contributions coming from the change of the cavity shape in the course of the reaction. By noting that around the TS configuration the reaction coordinate x (eq 16) corresponds to the motion of the CH3 group with the two chlorine atoms fixed, we conclude that during such a motion the cavity change is negligible.
U kcal/mol -110 -1 14
-1 18
-1 22
-1 26
-28.28
5. Nonequilibrium Solvation and the FES in the Total Polarization Space
1-0.0
From a calculation of the free energy profile in water we have obtained for a number of values of x , (CY = 1,2,3, ...) sets of the corresponding elementary surface charges Aql (the set corresponding to x , will be denoted as (A&) and also the set of solute charge distributions denoted as p,(r). The charges ( A q J , and p,(r) are equilibrated to each other at a given x , value. Next we introduce the medium coordinate Q by associating each of its values Q = Q, with the corresponding set (Aq,),. It was found convenient to assign to Q, the numerical value of the corresponding coordinate x = x,, thus accepting Q, = x, as a definition of numerical values of Q. This opens the possibility of calculating the nonequilibrium solvation energy. To do this, let us consider the specified set (AqJp,the corresponding medium coordinate Q = Qp = xp, and also the solute charge distribution pcrp(r) as calculated for the solute configuration x , in the field created by surface charges (AqJa. This calculation of pad proceeds without iterational mutual adjusting of pep and (AqJ8. (In the case of a self-consistent procedure, one would come out with an equilibrium situation with mutually adjusted quantities pa, and (Aq,),,.)The nonequilibrium situation generates the free energy, which is calculated by a modification of equilibrium formula 13: W x , , Q . & = A G o ( 4 ) +E,,(&)
-14.14
+ E,cor,(~B) + Er(BP) (17)
Here the first three termsare obtained for thesolute configuration x, with a nonequilibrium wave function and the corresponding charge pad obtained in a H F calculation with x = x p and Q = Qp.
The self-energyEr(B/3)is calculated by using equilibrium formula 15 at x = x,, and Q = Q 6 = x~ With x , fixed, the energy profile (eq 17) passes through a minimum at the equilibrium point x,,, Q, = x(,.
i
4 x=10.81 l O - * i
0
14.14
28.28
Q A lo4 +r3.M
j+x=14.14
lO-'i
; I t ~ 1 7 . 0 7l O - ' i
lO-'i
Figure 1. Freeenergy profilesalong themediumcoordinate Qfordifferent values of solute coordinate x .
We applied the procedure described above to calculate the two-dimensional FES in the space x , Q and evaluated the matrix of force constants at the saddle point x = Q = 0. It is listed in Table 111. The main qualitative result of this calculation was the observation that all FES cross sections along Q appeared to be single-well-shaped curves. Double-well profiles along Q were not found even in the vicinity of the saddle, where the solute system (eq 3) was expected to have the maximum polarizability. The calculated cross sections for different x, are shown in Figure 1.
It should be stressed here that such a simple FES topography is a necessary condition for the above described introduction of coordinate Q to be legitimate. This definition of Q would become ambiguous, should the cross sections along Q be double-well shaped, because one would observe in this case three stationary sets (Aqi) for every value of x . 6. Separation of the Noninertial Polarization
The FES in the subspace of coordinates x , Q,as obtained by formula 17, needs a further refinement. A change of independent medium variables from total polarization P(r) to inertial polarization &(r) is necessary. It is w e l l - k n ~ w n ~that ~ - ~the ~ total polarization consists of inertial and noninertial P,(r)) ingredients. The indices "in" and "a"will now be extended onto other inertial and noninertial field variables. So, as a consequence of the linearity of electrostatic equations, similar decompositions exist for the surface charge density u(r) and its field @(r):
Rate Constant for C1-
+ CH3Cl
-
ClCH3+ C1-
The Journal of Physical Chemistry, Vol. 97, No. 13, 1993 3213 (23-26) one can express any desirable quantity in terms of the independent variable @in. In this way we can define a general nonequilibrium free energy functional as (x being the gas-phase solute reaction coordinate)
The true medium dynamical variable (PI u, or @) is always associated with nuclear solvent motions (first of all, with the orientational motions of H2O molecules), therefore it is identified with the inertial field component. The conventional method of the separation of the polarization components formulated by Pekar)' reads
where D is the dielectric displacement and 4, is the optical dielectric permittivity. These formulas work only within the uniform medium model. With a cavity present, the field D(r) loses its importance and becomes useless for the theory. As shown earlier,29the polarization separation procedure depends on the cavity shape in this case. One can express a formal relation between field @(r)and solute charge p(r) as @(r) = hp(r)
(20)
where k is a linear symmetrical operator. An actual realization of the operator relation (eq 20) needs a solution of integral eq 5 with thesubsequent application offonnula 9. Onecan introduce, by analogy, the similar operator K , making a change eo t, in formula 20. So, the noninertial reaction field is
-
a m ( r )= h m p ( r )
(21)
whereas the inertial field appears as a difference: @in(r) = (h-km)p(r)
(22)
This reasoning works only for an equilibrium situation when all reaction fields @, am, and a,, are equilibrated to one and the same charge p(r) found in a self-consistent fashion. Equation 22 fails to be valid for a nonequilibrium case. One may, however, introduce the new charge quantity pln(r) representing the selfcharge of the field @,, which is implied to be equilibrated to ai,,: aIn(r) = (k- h,)~,,,(r)
(23)
Here, the first two terms represent self-energies of fields @in and @m29 and the third term is the electrostatic interaction between the total field @ and the solute charge p as obtained according to the quantum-chemical calculation (eq 24). The latter quantity, AGO, is the solute self-energy calculated as a total energy corresponding to the HF procedure (eq 24) and depending on @in via a solution of this equation. To proceed further one needs to explicitly define the nonequilibrium charge density pmin eq 2 1. Such a definition cannot be consistent in general.35-39 There exist, however, two limiting cases when p, can be safely introduced. (a) One sets p- = p; that is to say, the noninertial field is adjusted to the solute charge distribution. Now becomes dependent on @in, and the system of equations (24 and 18) becomes nonlinear relative to @in;it should be solved by iterations. This case was called the selfconsistent (SC) a p p r ~ x i m a t i o n . (b) ~ ~ One denotes as pcq(x)the solute charge equilibrium at a given value of thecoordinate x and then sets p, = p&). Now 9, becomes independent of ai,,, and eq 24 becomes linear relative to ai,,. This is the so called BornOppenheimer (BO) appro~imation.3~ Note that such a definition of pm is legitimate only if the solute wave function is treated as a single Slater determinant, according to the HF scheme. If several determinants contribute to the wave function then in the BO limit no single charge distribution exists to which @, can be adjusted; rather it should be equilibrated to different charges paa, each corresponding to a given matrix element Ha, of the C I Hamiltonian matrix in terms of the basic determinant functions Do and Db. In other words the field amin the BO limit should be treated as an operator q ~ a n t i t y . ~ ~ . ~ ~ 7. The FES Calculation in BO and SC Limits
Based on the methodology developed in section 4, we now generalize the algorithm of the nonequilibrium free energy calculation, as given in section 5 , to a more complicated case involving the separation of the noninertial polarization. Let us begin with a simple BO limit. Consider two sets of mutually equilibrated fields:
Obviously, pln # p in a nonequilibrium case since p is defined in quite another way, as a solute charge obtained from the H F equation
bn= (Po-
= (24) are their energies, PObeing
where cc, are the solute MOs and 4, the gas-phase Fock operator. The external field is defined as a sum (eq 18), and what remains to be done is to decide at what charge p d r ) the field 0, is equilibrated. Given p,, one can extend eq 21 to am(r) = hmpm(r)
(25) Allchargescoincideat theequilibrium: p = pln= p-. In a general nonequilibrium situapon y e also i!trodu_ces two auxiliary fields, both adjusted to pln,@ = Kp,, and 0, = K,p,,,, and obtains from eqs 23 and 18
- -
a,,,= 3- 3, = @ - am
(26) The fields & and b, correspond toa fictitious equilibrium situation, as defined by a given field a,, and differ from nonequilibrium fields (q18) and am(eq 25). By using the set of equations
@in(2)
= (h- km)p2 (28)
Here the treatment of fields bo) and p I proceeds as in section 4, so that @ , ( I ) is adjusted to p l . The last quantity represents the self-charge of field @ i n ( l ) considered as a primary given variable. (The same is true for the quantities labeled by index 2.) The nonequilibrium energy is calculated as
+ @ - ( ' ) ) p d3r + A G o ( ~ l , @ i n ( 2 ) ) BO limit
JV,(einI2)
(29)
The meaning of separate contributions here is the same as in formula 27. The self-energies of fields @i,,(2) and @,(I) must be calculated with their own self-charges. (Note that after specifying the field and its self-charge the self-energy is calculated via
3214
Basilevsky et al.
The Journal of Physical Chemistry, Vol. 97, No. 13, 1993
7.0
7.0
4.0
4 0
.o
10
1
-2.0
-2 0
-5.0
-5.0
-8.0
-8 0
curve with a barrier having its maximum at x = 0,fo-1 is the inverse Pekar factor appearing above in eq 19, and Po(r,x) is the equilibrium polarization function at a given x. The equations of motion then are:4i.42 (3-fo)P = - -+ (Gaussian random force) 6P
Figure 2. Free energy surface (kcal/mol) in the vicinity of the saddle point x = 0, Q = 0. (The BO limit.)
equilibrium formula 15.) The charge p in the third term should be specially found in a H F calculation (eq 24) with the total field 9 = a,,(*) @ , ( I ) and without electrostatic iterations. Obviously, p = pI = p2 at equilibrium. The same quantum-chemical calculation yields the energy AGO,the fourth term of eq 29. Finally, in order to obtain the third term, one should calculate, as in the case of a self-energy calculation, the quantities E,, and E,,,,, defined by relations 14. Consider now the SC limit. We utilize again the techniques of equilibrium calculations, although changing the exact meaning of the respective fields, according to
+
x = x2: p2 = p(x2);
Ql,(2)
AG(x~,@~,,(~)) = -I/2Jv,@,n(2)p2d3r - I/2Jv,@mp d3r +
+ @,)p d3r + AGo(~,,91n(2)) SC limit
(30)
8. The Dynamical Equation and Its Parameters
The phenomenological equation of motion for our problem is formulated below in the framework of the uniform medium approximation, ignoring boundary effects of the cavity. Then theinertial polarization Pi,(r) 1P(r) (index'in" is canceled below) becomes a natural medium variable. Let us represent the energy functional in the form
+ '/2f0Jd3r
- 7)(...) d7
(33)
Po(r,x) = a(r)
+ xb(r)
(34) By calculating derivatives of functional 3 1 and performing standard transformationsg2the discrete medium coordinate Q is revealed as a functional of polarization: Q=
= Q[P]
Y = J ( P - a)b d3r; @ = Sb' d3r
Note that thequantities AGOin eqs 29 and 30 are different because eq 24 is solved with different fields 9, in these two cases. The FESs calculated by means of prescriptions 29 and 30, as found in our computations, differed quite insignificantly. Their cross sections along the medium coordinate Q are single-well curves qualitatively resembling those shown in Figure 1 for the total polarization calculation. However, their numerical distinctions from the FES (eq 17) reach several kcal/mol far away from the saddle point. At the saddle point (x = Q = 0) functions 17, 29, and 30 coincide, in accordance with the earlier proved theorem.29 Figure 2 displays the contour map of the BO limit FES. Finally, the values of force constants at the TS point for all three FESs are collected in Table 111.
AG[x,P] = U,(x)
3(...) = J-'Jr
We invoke a Taylor expansion of Po(r,x) at point x = 0:
= (h- hJp2
Now the field 9-is adjusted to thesolutecharge: CP, = k p . This charge, as before, comes from eq 24 which now becomes nonlinear because its matrix elements depend on e, the latter quantity depending onp. By this means fields 9,andp becomecomplicated nonlinear functionals of a given field = and their selfconsistent treatment is necessary. The formula for the nonequilibrium free energy is
Jvh(@ln(2)
whereJ= (4&)/(; - e,) is the basic dynamical operator and is an integral dielectric permittivity operator conventionally specified by the complex permittivity function e ( ~ ) . These ~ ~ definitions imply the following relations:
( P - Po)'
(35)
We also come out with the closed set of dynamical equations for x and Q:
Q-foMQ
= -foB(Q- X) + F ( t )
mx=f&(Q-x)--
dV, dx
(36)
Here F ( t ) is a Gaussian random force, the result of space integration with the weight b(r) of the Gaussian force entering the right-hand part of eq 32. Finally, regular forces on the right of eq 36 can be extracted by partial differentiation of the following potential function U(x,Q) = V,(X) + '/2f&(Q- x ) ~ (37) Equations 35-37 were derived earlier by B e r e z h k ~ v s k yfor ~ ~the Debye medium model. Most important for the following is the fact that the coordinate Qof eq 35 exactly coincides with our collective medium coordinate introduced in section 5. This is shown in Appendix 1. The immediate consequence is that, by invoking the uniform medium model with the simple energy function (eq 31), one can write in the explicit form of eq 37 the FES of reaction 1 as computed in section 7. We expect this representation to be sufficient around the saddle point x = Q = 0. So we derive from eq 37
fLl=&-t> 1 1 Here x is the gas-phase reaction coordinate, Ul(x) is a potential
where X = (dW,)/(dx2) a t x = 0. Now wecan evaluate the basic parameters of the phenomenological equation (eq 36) from the
+ CH,CI
Rate Constant for Cl-
-
CICH3+ C1-
The Journal of Physical Chemistry, Vol. 97, No. 13, 1993 3275
data listed in Table 111:
B = 3.0 kcal/(mol A*)
X = -210 kcal/(mol A') We also note a remarkable symmetry relation (d2U)/(aQ2) = -(aW)/(aQdx) following from eq 37. It appears as a consequence of the general form of the free energy function taken as eq 3 1. One can verify from the data of Table 111 that this symmetry is indeed obeyed with a high precision, which underlies our choice of formula 3 1 as a reasonable approximation. 9. Tbe Transmission Factor Let us write the rate constant of reaction 1 in the usual form kBT
'transl'
the equations of motion (eq 36) at the reactant minimum and at the TS. Among the TS frequencies there exists a single purely imaginary frequency having the positive imaginary part; it has the explicit form wg* = is2 thus defining the quantity s2 in eq 42. Now we assume that the reaction coordinate x is a normal mode of the solute and out of all the solute normal modes it is only the quantity x that is coupled to solvent variables. Moreover, we assume that this coupling operates only in the TS, being negligible for the initial state. (Note that in the initial state x must be different from the antisymmetrical expression (eq 16).) Under these assumptions, after some manipulation, we obtain the transmission factor as
'vib'
'rot'
K=xh--(0) z (0)z (0) 'trans1 rot vib ,
Here x is the transmission factor and Zbo)and Zi* are the partition functions of the reacting solute in its initial (isolated reagents) and transition states. The explicit inclusion of these partition functions is necessary because only the solvent coordinates contribute to the temperature dependence of the activation free energy AG* as calculated in section 4. Strictly speaking, the external medium field must be taken into account when the partition functions are calculated. In our treatment this complication was ignored, the partition functions being treated like those for the gas-phase case. The corresponding inaccuracy of the TST preexponential factor is expected to be relatively small. We can calculate the TST prefactor as A=---kBT
' '
'transl'
'vib'
'rot'
(0)
rot
'trans1
(0)
(0)
,
+ CH,Cl
k
[Cl-
+ CH,Cl]
k i
\ I1
medium+x
The product in the numerator counts only the isolated medium frequencies which are the roots of the free medium characteristic equation: f(w) = o (44) The frequencies in the denominator belong to the combined system "coordinate x coupled to medium coordinates at TS" and are the roots of the TS characteristic equation:
det F(w) = 0
(40)
vib
which is equal to A = 6.3 X los L/(mol s) (the PM3 calculation of the gas-phase frequencies and the moments of inertia) and A = 7.2 X lox L/(mol s) (the ab initio calculation with specially scaled f r e q ~ e n c i e s ~ ~One ) . may conclude that the value of A is indeed insensitive to details of a quantum-chemical calculation. The correction for zero-point temperature vibrations decreases the activation free energy in eq 39 by -0.9 kcal/mol. The main medium influence on the prefactor comes via the transmission factor x . We have calculated it by means of the two-dimensional version of the Kramers-Grote-Hynes (KGH) theory, which is also referred to as the generalized transitionstate theory (GTST).42.45-47 We postulate the following reaction mechanism: CI-
(43)
-
The mass corresponds to the antisymmetrical mode (eq 16):
m=
(41)
We separate here as an intermediate the state of reactants as a contact pair in the solvent cage (denoted as [...I) assuming this state to be at equilibrium with the state of isolated reactants, in the spirit of the conventional TS theory. The rate constant becomes (kI >> k 2 )
mCH,
(46)
+ 2mCl
The actual evaluation of expression 43 is described in Appendix
2. The so specified calculation was performed with the following two complex permittivity functions: e(w)
k2
products
mCH,mCI
e(u) = e,
= e,
- e, +1 - iwrD €0
Debye model
e0 + +1 - iwrD em
- e,)
(e,
+ w,'~,' (1 - iwr,)' + w,'r; 1 - iwr,
Saxton model4*
The following parameters values were used:48
k
K = s K z k- I
and it sufficesto find the transmission factor for the rate constant K 2 . According to the GTST42 K2 =
2)
E(- 7) exp( rIWY,O)
rI%
IJ2
= 3.84 X lO-I4 S; w, = 6.67 X 10l2 s-I
The results so obtained are
(42)
2* Here up(o) and up* are the characteristic frequencies of the combined system 'the reacting solute the solvent" related to the initial state (the reactant state in the solvent cage) and the TS, respectively. These frequencies are complex-valued quantities found as the roots of the characteristic equation extracted from
+
7,
x
= 0.81
Debye model
x
= 0.85
Saxton model
An important comment is that the characteristics of the state of the reacting particles in the solvent cage were eliminated from the calculation of Y . This is a consequence of the equilibrium hypothesis at the stage of formation of the contact reactant pair,
Basilevsky et al.
3276 The Journal of Physical Chemistry. Vol. 97, No. 13, 1993 according to mechanism 41. So, the rate constant is determined only by the properties of the TS and the isolated reagents, just as in the canonical TS theory.
10. Discussion In the present paper the equilibrium and nonequilibrium free energy changes were calculated for reaction 1 using a regular quantum-chemical methodology. This is the main distinction of our work from the earlier continuum treatment^^^.^^ where the shape of the free energy profile along the medium coordinate was postulated beforehand, without an underlying calculation. The explicit construction of the FES in the coordinate subspace including coordinates of both the reacting solute and the medium provided a background for the subsequent evaluation of the transmission factor in terms of two-dimensional KGH theory. An FES calculation of this type may be of general interest at the level of studying mechanisms of charge-transfer reactions proceeding in polar solvents. A classification scheme can be proposed according to which one may distinguish whether the FES cross section along the medium coordinate drawn in the vicinity of the TS is double-well or single-well shaped. In the first case one may state that the reaction proceeds via an electrontransfer mechanism, anticipating, as one of the consequences, nontrivial changes in the rate constant prefactor. In the second case, which has actually occurred in the present study, one concludes that the motion along the solute reaction coordinate is, as in the gas phase, the driving force of a chemical conversion. Then the medium influence on the reaction mechanism becomes insignificant. The collective coordinate Q introduced in secfion 5 is an appropriate medium variable as applied to FESs of the second (single-well-shaped respective to medium changes) type. After a linearization at the saddle point it becomes identical to the coordinate deduced earlier by Berezhk~vsky.~~ The linear approximation allows one to obtain a closed system of equations of motion, minimizing the number of dynamically essential variables. The coordinates (x,Q) are optimally fit to calculate the rate constant in the KGH approximation because it is this approximate theory that invokes a complete linearlization of all regular forces around FES stationary points. Away from the saddle point, where the linearization is hardly legitimate, our coordinateQdiffers from Berezhkovsky’scoordinate. Our version seems to be more convenient for practical purposes because the definition procedure (see section 5 ) , being universal, does not require an evaluationof specific integrals which appear in formula 35. We have to emphasize also that any attempt to circumvent the KGH approximation inevitably means that the linearization of regular forces is discarded. As soon as the barrier of reaction 1 is high enough, the multidimensionalKGH theory is an adequate tool for calculating the rate constant. Therefore, in terms of the accepted reaction model, the transmission factor value w H 0.8, found in section 9, looks like an accurate result. Its discrepancy from the value x = 0.55 found in MD calculations,’6-’*in our opinion, reflects the differences between the continuum (our case) and discrete (MDcalculations) models of the medium. Each of the models has its merits and disadvantages. The drawbacks of thecontinuum model are obvious: the representation of the water solvent as a dielectric continuum without an internal structure is a severe idealization. In particular, one may expect that the hydrodynamical type friction effects, completely ignored in our model, may be of real importance. Among the deficiencies of the discrete model the most serious one is, probably, that it cannot, for technical reasons, consistently separate the noninertial polarization effects. Taking into account only orientational motions of solvent molecules, the M D simulation thus requires that the relaxation of the electronic polarization, which must be instantaneous, actually follow the slow orientational relaxation. In order to
estimate the latter error in the frameworkof our continuum model we repeated the dynamical part of our calculations without separating the noninertial field. In doing so, we used the total polarization force constants generated by FES (eq 17) (see Table 1II)andset e, = 1 (insteadofthetruevaluee, = 1.8)indynamical equation 36. As a result, we obtained for the Debye medium model the new value x = 0.69 which is much closer to the MD result. Roughly speaking, the ratio 0.5k0.69 is the measure of the inaccuracy of our continuum model whereas the other ratio 0.69:0.81 reflects the deficiency of the MD calculations. Another attempt toestimate the transmission factor of reaction 1 was published by LeeandHynesin termsofacontinuummedium model.49 In this work the nonequilibrium free energy values were calculated by the phenomenological formula of the same type as eq 3 1. The dynamical treatment was essentially different from that adopted in the present paper: a different medium coordinate was introduced with a different (one-dimensional and nondissipative) resulting equationof motion. Theobtained transmission factor varied from 0.64 to 1.O depending on the parametrization of the dynamical equation. So, despite a significant difference in the methodology, the ultimate result of this calculation is close to that obtained in the present paper. That is why it is likely that expression 3 1 is fairly satisfactory as a phenomenological model of the free energy functional for reaction 1,although in its structure several circumstances are ignored, such as the boundary effects of the cavity and the nonlinearity of the medium response to the soluteelectricfield. Wecanalsoconcludethat the finalmagnitude of the transmission factor is quite insensitive to the details of the dynamical treatment. Probably, the most important conclusion must be that the true transmission factor for reaction 1 does not significantly differ from unity, so that the canonical TS theory is indeed a very good approximation. As we show in Appendix 3, this is equally true for any FES having a typical single-well shape along the medium coordinate. When calculating the rate coefficient for such reactions it suffices to proceed according to the old canonical prescription,’ only taking into account the modification of the gas-phase potential surface by the field of solvent molecules equilibrated to the solvated chemical substrate. Altogether, the rate constant prefactor as calculated according to the prescription of section 9 falls in the bottom of the range of the respective experimental values for simple S Nreactions ~ with halogen atoms.12
Acknowledgment. Theauthors aregrateful to Marshall Newton for helpful discussions. Appendix 1 Analytical Expressionof tbe Medium Coordinate in the Vicinity of the Saddle Point. Consider the quantity L(r,x) which depends both on the chemical coordinate x and space variable r. The typical examples are equilibrium field variables such as surface charge density qn(r,x),the corresponding reaction field @in(r,x) or the origin of these two fields, the equilibrium inertial polarization Po(r,x). The total (inertial + noninertial) quantities of this kind were used in section 5 for the definition of coordinate
Q. The Taylor expansion at point x = 0 gives
+
L(r,x) = L(r,O) xL’(r,O) (47) Repeating the reasoning of section 5 we define the coordinate Q as L(r,x=Q) - L(r,O) L’(r,O) The following identity is easily deduced from relations 47 and 48
Q=
Rate Constant for CI-
Q=
+ CH3CI
-
CICH3+ C1-
The Journal of Physical Chemistry, Vol. 97, No. 13, 1993 3277
$[L(r,x=Q) - L(r,O)]L'(r,O) d3r (49)
$ [L'(r,O)]' d3r
Note that no specification of quantity L is needed to obtain this result. By substituting here the equilibrium inertial polarization: L(r,x) Po(r,x) we come Out with 35. One bear in mind that relations 48 and 49 are legitimate only in the linear approximation (eq 47).
supplies us with the estimation 1 - x < which is in accord with the initial assumption (eq 52). We emphasize that the present derivation was essentially based on the KGH theory. So the corresponding conclusion is valid only when the anisotropy of time scales associated with the coordinates and is not significantand the free energy barrier (in the units of RT) is sufficiently large.50.51 These conditions Seem to be obeyed for the of the sNz reaction considered here.
Appendix 2
References and Notes
Derivation of Formula 43. To facilitate the calculation of the transmission factor according to formula 43 we separate in the characteristic polynomials 44 and 45 the terms having the leading power in w :
( I ) Glasstone, S.;Laidler, K. J.; Eyring, H. The fheoryojrateprocesses; McGraw-Hill: New York, 1941. (2) Hynes, J. T. In Theory of Chemical Reacfion Dynamics; Baer, M., Ed.; C R C Press Inc.: Boca Raton, FL, 1985; Vol. 4, p 171. (3) Morokuma, K. J. Am. Chem. Soc. 1982, 104, 3732. (4) Jaume, J.; Lluch, J. M.; Oliva, A.; Bertran, J. Chem. Phys. Left. 1984, 106, 232. (5) Gonzales-Lafont, A.; Truong, T . N.; Truhlar, D. G. J. Phys. Chem. 1991, 95, 4618. (6) Shaik, S. S. Acfa Chem. Scand. 1990, 44, 205. (7) Shaik, S. S. J . Am. Chem. SOC.1984, 106, 1227. (8) Burshtein, K. Y. J . Mol. Slruct. (THEOCHEM) 1987, 153, 209. (9) Kozaki, T.; Morishashi, K.; Kikuchi, 0. J . Am. Chem. SOC.1989, 111, 1457. (IO) Kong, Y. S.; Jhon, M. S. Theor. Chim. Acta 1986, 70, 123. (11) Tucker, S. C.;Truhlar, D. G. J. Am. Chem. SOC.1990, 112, 3338, 3347. ( I 2) German, E. D.; Kuznetsov, A. M. J. Chem. Soc., Faraday Trans. 2 1986.82, 1885. (13) Kuznetsov, A. M. J. Phys. Chem. 1990, 94, 8664. (14) Chandrasekhar, J.; Smith, S. F.; Jorgensen, J . W. J . Am. Chem. SOC.
det F(w) =
+ ...)
(A0O"+2
R(w)
where R ( w ) is a common denominator appearing because f l w ) is represented as a ratio of polynomials and n is the order of the polynomial appearing in its numerator. It easy to see that Ao/ao = -m& The following identity is true:
rI
me&uK+x
*
det F(w=O)
a0
*fi
So we infer that the calculation of the transmission factor x apmrding to eq43 reduces to an extraction of the decay frequency 0 from the characteristic equation (eq 45). Then, given function f l w ) , the ratio of frequencies is easily computed by eq 51. Appendix 3
TbeApproximate Estimationof the Transmission Factor When the FES Profile along the Medium CoordinateComprisesa SigleWell Curve. Let us confine our consideration to the case when the transmission factor is close to unity: 1-x
