University of Southern California Dissertations and Theses

Computer simulation of chemical reactions in aqueous solutions and biological systems

(USC Thesis Other)

Computer simulation of chemical reactions in aqueous solutions and biological systems

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content COMPUTER SIMULATION OF CHEMICAL REACTIONS IN AQUEOUS SOLUTIONS AND BIOLOGICAL SYSTEMS by Jenn-Kang Hwang A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Chemistry— Chemical Physics) December 1988 Copyright 19-88 Jenn-Kang Hwang UMI Number: DP21970 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a com plete manuscript and there are missing pages, th ese will be noted. Also, if material had to be removed, a note will indicate the deletion. UMI DP21970 Published by ProQuest LLC (2014). Copyright in the Dissertation held by the Author. Dissertation Publishing Microform Edition © ProQuest LLC. All rights reserved. This work is protected against unauthorized copying under Title 17, United States Code ProQuest LLC. 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, M l 4 8 1 0 6 -1 3 4 6 UNIVERSITY OF SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIVERSITY PARK LOS ANGELES, CAUFORNIA 90089 under the direction of h..!.?. Dissertation Committee, and approved by all its members, has been presented to and accepted by The Graduate School, in partial fulfillm ent of re quirem ents for the degree of P K .P , C This dissertation, written by Jenn-Kang Hwancf D O C TO R OF PH ILOSOPH Y Dean of Graduate Studies Date Decem ber 2 0 , 1988 DISSERTATION COMMITTEE , A r - U h Chairperson To My Parents ACKNOWLEDGEMENTS I would like to express my gratitude to my wife Shan-Shan for her cheerful support through the years. Table of Contents Page DEDICATION ii ACKNOWLEDGEMENTS iii LIST OF TABLES vi LIST OF FIGURES vii ABSTRACT ix CHAPTER 1 General Formulation 1 1.1 Introduction 1 1.2 Semiclassical Green’s Function Method 2 1.3 Semiclassical Trajectory Approach 6 1.4 Application to Multidimensional System 14 1.5 The indirect Method 20 1.6 References 23 CHAPTER 2 Simulations of Electron Transfer in Water 26 2.1 Introduction 26 2.2 General Consideration 28 2.3 Computer Simulation of the Model System 29 2.4 Free Energy Perturbation Calculations 33 2.5 How Quadratic is the Free Energy Surface 35 2.6 The Deviations from Marcus’ Relationship 41 2.7 Discussion 42 2.8 References 45 CHAPTER 3 Energetics and Dynamics of SN2 Reaction 48 3.1 Introduction 48 3.2 Method 52 iv 3.3 Computational Details 3.4 Dynamics of SN2 Reaction 3.5 Free Energy Relationships 3.6 Dynamics Effect and the Solute-solvent Coupling 3.7 Results and Discussion 3.8 References CHAPTER 4 Simulations of Enzymatic Reaction 4.1 Introduction 4.2 Method 4.3 Free Energy Changes of Genetically Modified Subtilisin 4.4 Concluding Remark 4.5 References APPENDIX A List of Tables Table 2.1 Parameters used in calculations of benzene - water system 2.2 The Franck-Condon factors of quinone and biphenyl 3.1 Parameters used in Sn2 calculations List of Figures Figure Page 1.1 SGF and quantum mechanical rate constants 7 1.2 Vibronic transition of I2 8 1.3 ST and SGF rate constants in a one dimensional system 10 1.4 Quantum mechanical and ST rate constant for a two dimensional system 15 1.5 Classical trajectory on a two dimensional system 16 1.6 Nature of trajectories that contribute te kST 18 1.7 ST and SGF rate constants in a two dimensional system 19 1.8 Schematic illustrations of direct method and indirect method 21 2.1 Schematic description of electron transfer reaction 32 2.2 Snap shot of the configure of benzene-water system 34 2.3 Free energy profile of electron transfer in benzene-water system 36 2.4 Free energy profile of D“A and DA" states as a function of energy gap 37 2.5 Free energy profile of ET for AG0=0 kcal/mol and -10 kcal/mol 38 2.6 Calculated and Marcus’ free energy surfaces 39 2.7 The effect of Vibronic channels on the rate constant 40 3.1 Problems with calculations using gas phase charge distribution 50 3.2 Gas phase EVB potential of SN2 reaction 52 3.3 Polarization of solvent dipoles and energetics Sn2 reaction 59 3.4 Free energy profile as a function of X of SN2 reaction 64 3.5 Free energy profile as a function of energy gap of SN2 reaction 65 3.6 Free energy functional Agj and Ag2 as a function of energy gap 66 3.7 Relationship between activation energy and AGq. 75 3.8 Dependence of activation energy of SN2 reaction on solvation energy 78 3.9 The free energy surfaces of the solvent driven and the solute driven reaction 81 3.10 The solvent fluctuation and the activation barrier of S^2 reaction 82 3.11 Nature of the downhill trajectories for the SN2 reaction 85 vii 3.12 Autocorrelation of S^2 reaction 86 3.13 Time dependence of the solvent coordinate for the downhill trajectories 88 3.14 Time dependence of the solute coordinate for the downhill trajectories 89 3.15 Time dependence of the charge on the attacking Q ” atom 90 4.1 Schematic illustrations of the modelling of enzymatic reactions 100 4.2 Diagram showing the catalytic mechanism of subtilisin. 104 4.3 The snap shots of the active sites of native and mutant subtilisin 106 4.4 Free energy profile of the reaction of the native and Thr-155 subtilisms 111 4.5 Calculated free energy changes of the native and mutant subtilisins 113 viii ABSTRACT A general approach for simulations of chemical reactions in aqueous solution and in enzyme is developed and examined. The chemical reactions are described by resonance structures by means of the Empirical Valence Bond method, and reactions are treated as crossing between different resonance states. The probability of reaching the transition state is related to the activation free energy. To evaluate the activation free energy barrier, we use two powerful approaches to sample the phase space, (a) the direct approach (the semiclassical trajectory method) and (b) the indirect approach (the free energy perturbation method). In the semiclassical method the trajectories are propagated directly on one of the resonance states and the transition will occur whenever the surfaces cross. This approach can be applied to the case with small coupling between states, i.e., small activation barrier. In cases of large barriers, we have to switch to the free energy perturbation method, which forces the system to move along the reaction coordinate by means of the artificial mapping potentials. First chapter gives an introduction to the general formalism. The second chapter explores the electron transfer reactions in water by both the semiclassical trajectory method and the free energy perturbation method, and the relation between the Marcus’ theory and the result of our microscopic simulation is discussed. Chapter 3 describes the energetics and dynamics of the Sn2 reactions in aqueous solutions. The last chapter applies our Free Energy Perturbation method to the enzymatic reactions. The method is examined in simulations of the catalytic reactions in the native and in several artificially mutated subtilisin. ix Chapter 1 The General Formulation 1.1. Introduction The study of chemical reactions in solution has attracted significant theoretical interest in recent years1-26. This interest reflects in part the realization that many chemical and most biological reactions occur in the condensed phase, and that gas phase calculations cannot be expected to reliably describe condensed phase chemistry. Although great insight into the nature of solvent effects has been provided by early macroscopic approaches20’ 23-25, it seems now that progress of a more quantitative nature would require detailed, microscopic approaches. The dynamical role of the solvent should be described on a more detailed microscopic level. This is especially important in studies of charge transfer reactions, in which the energetics and the dynamics are determined by the strong interactions between the solute and the solvent. Despite the crucial role of the solvent in charge transfer reactions, it is not easy to include the solvent in microscopic quantum mechanical calculations of such reactions. The chemical reactions can be described as transitions between the reactant and the product states by the master equation I dNR / dt = - k+NR+k_Np (1.1) where and Np designate the concentrations of reactants and products respectively. The k+ (which will called k hereafter) can be written as27 k = kTSjF (1.2) ktst 1s the rate constant predicted by the transition state theory28 1 kTST = (kBT/h) exp(— Ag * p) (1.3) where P = l/k£T, kB is the Boltzmann constant and Ag* is the activation free energy. F in Eq. (1.2) is the transmission factor that accounts for the fraction of trajectories that will cross to the product state from the transition state. We will attempt to identify the factors that determine kTST and F and to relate them to the actual microscopic features of the reacting system. The quantum mechanical transition rate between electronic state a and electronic state b is given by the well-known expression k = C t m I ^ K < » n W amj J l%Eam-E bJ (1.4) where P(am) is the Boltzmann factor for the initial state and Vam bn is the matrix- element of the total Hamiltonian between am and bn zero-order states. To evaluate Eq. (2.4) one has to solve the Hamiltonian of the system, but for a multidimensional anharmonic system with thousands of degree of freedom, such a calculation is entirely beyond the limits of current computational ability. In order to overcome such difficulty the semiclassical approaches was developed, which can be very easily generalized to the multidimensional system. We will start by considering the Semiclassical Green’s Function (SFG) method2^, and from which we will derive the Semiclassical Trajectory (ST) method16ft,c. 1.2. Semiclassical Green’s Function Method The basic object is to evaluate the transition amplitude from state a to state b. To attack this problem we write the time dependent wave function of the system as y = Ca§a(r) exp {- W ) \ } + C ^ b{r) exp{-(ifh)^dt'Vb(t')} (2.1) where r describes the solvent coordinates and Va (r) is the potential of the a th state. Substituting this equation into the time dependent Schrodinger equation 2 ifKdy/dt) = //\|/ (2.2) we get Ca(x) = j ^ - W ) o exp {( m j ‘ o dt'AVab(r(t'))} Cb(t)dt Cb(x) = Jo T - ( ^ ) a exp {(i/fi) dt'AVba{v(t'))} Ca(t)dt (2.3) where AVba = Vb - VQ and C = <§a\H\§b>. Cb(t) is. the time dependent transition amplitude from state a to state b. With the initial conditions Cfl(0) = 1, Cf e (0) = 0, and for the time range when Cfl(t) ~ 1 one obtains Cb(x )= jW i/h )a e x p [W )d t'A V ba(t')]dt (2.4) When c is large we switch to the adiabatic representation166. In such case the expression for the transition amplitude is Cb{t) = J V r « y a (j)/a r > e x p [M )jtdt'AVba(t')']dt (2.5) We will first concentrate on the cases where diabatic approximation is valid and discuss the adiabatic case later. It should be noted that |C^(r)|2 is the possibility of moving from potential surface Va to potential surface V b along a yet unspecified path r(0- Since Cb can be identified as the propagator35 from xa in potential a to r^ in potential b, it can be reformulated by the Feynman’s path integral3 1 1 . The transition amplitude from ra(0) to rb(x) can be described by Feynman’s path integral K(rafb ’x) ~ <rb I exp[-(i/K)]Hx | ra> = J r*D r exp[{m ^L {r{t'))dt’} (2.6) ra ta where Dr designates integration over all possible paths which satisfy the boundary condition, r(0)=ra, r(x)=rb,and L is the classical Lagrangian, 3 L - T - V (2.7) where T is the kinetic energy and V is the potential energy of the system. In the stationary phase approximation, it can be shown35 that the paths r(t)’s that contribute the most to the Feynman’s path integral are the classical trajectories. The K(rb,ra,x) = - a m d X j K - i m - 1# S ^ r/x )) / dradrb]^ {K ih f [dA Vba{rft))ldt]1 ’ ^ } where Sc{rjlf)) is the classical action r}(t), r| • is the. number of focal points of the trajectory rft) and t* is the time when Va and Wb crosses. K(rbf a,x) is the probability amplitude of moving from ia to ib at a given, time x: Note that with the boundary condition r(0) = rQ and r(x) = rb, the trajectories, are not uniquelydefined.. It is. easier for us to work with the energy representation of K{r.bf a,% ). This can be done by carrying out the Fourier transformation of K(rb,ra,x), where the summation is over trajectories which begin from rQ and end at rb with the a little algebraic manipulation, the expression for G{rb,ra,E) can be rewritten as semiclassical approximation29 of K(rb,ra,x) is x exp{{ifh)[Ex-2 J ^ W a(r/r'))-2j^ (2.8) G(rb,raJE) = ~ W )\~ d x exp{m)Ex]K(rb,ra,x) (2.9) The semiclassical approximation29 of G(rb,raJE) is G(rb,ra,E) = - ^ a S y[ ^ ( ^ ) H ( ^ r 1/2][(-2^ ^ ) “1 325d ( ^ ) ) / 3 ^ ] 1/2 (2.10) given energy E, Pa = ^2m (E-Va) and r* is the crossing point of the potentials. With G(V a £ ) = - ^ o Y ll{ P ^ r y blirbr ln][l,-2Tamrli % f . r ^ ^ ra ^ a 4 x^p(/epF(e°)F(e°) (2.11) where I designate different types of trajectories and 9; specifies the phase of the particular type of trajectory. The function Fa is given by F(x(Q) = 'L Jaexp(ijaQ ° a) (2.12) where 6 ° = J (1 fr)Pa dq — 7 i: . Due to destructive interference Fa will become zero unless 0^ satisfies the Bohr-Sommerfeld quantization condition. The rate constant can be written as, k = jdE a(E )k(E ) (2.13) where k(E) is given by29 « £ ) “ < < „ ,/ I f„ (9 > f> (9“) I 2 < 2-I4> where G in a„ is the harmonic frequency of Va and Wuni is the semiclassical uniform approximation to the probability of surface crossing. Note that the last term is virtually zero if energy is not on the quantized level. Ww n i- is given by29 Wuni = { ^ !4/[2PaPb(dPa/d r - d P b/dr)]^2}Ai(-^) (2.15) where Pa = [2m(E— Va)]^2 and Az(— £,) is the Airy function which is defined as A /(-^ )= (l/3 W l{/1/3[(3 /2 )^ 3] + /_ 1/3[(3 /2 )^ 3] } (2.16) where ^=±(3/2)[f r\ p am d q A rHPbm d q b\2^ (2.17) rn r The ± signs depens on whether the slope of Va and Y b are of the same sign or 5 opposite sign. With a one-dimensional harmonic model system, we compare the dependence of SGF and quantum mechanical rate constants of on the reaction "exothermicity". It can be seen from Figure 1.1 that the agreement is very good. We also apply SGF to the one dimensional anharmonic system, X — *B vibronic transition in I2. This system is approximated by the Morse potentials. From Fig. 1.2 it can be seen that anharmonicity is correctly reproduced. The SGF approach is excellent in the one dimensional system but it becomes extremely complicated when applied to the systems with more than one dimension, which makes this approach impractical to the real systems. In the next section we will consider a more general approach which can be easily applied to the multidimensional system. 1.3. Semiclassical Trajectory Approach The SGF approach can be extended to several dimensions but at present it does not provide a practical way to treat quantum mechanical effects in large system. In order to relate the SGF to the more convenient approximations we have to compromise and to relax the requirement that the nuclear motion is quantized in the general anharmonic case. A hint for a possible approximation is obtained by considering the quantized SGF approach in the simple case where both Va and Wb are harmonic potentials of equal curvature. In this case one obtains W W - " ' L a * ) o (3.D Cb{t) can be related to rate constant, k(E)ST = (l/x)\Cb(x)\2 = (l/T )(a/ft)2|J J dt exp{m bat + d f «(r(r'»}|2 (3.2) where u = AVba— <AVba>Q and ^ ba = <AVba>Q. The relation between kST and kSGF can be realized by expanding the phase of the exponential function around its stationary points. This gives ^dt'AVba(t') = j^dt'AVba(n + 6 Rote constont (relative units) 0 . 108- 1000 2000 3 0 0 0 4 0 0 0 -A V ° 0 (cm-1) Figure 1.1 A comparison of the (dotted line) and quantum mechanical result (solid line). The calculations correspond to the zero temperature limit. The system is = (l/2)fco>?2 and Vb = (l/Z )*© ^ -*.)2 with to = 343cm-1 and X = 2.1. 7 I(w 0 )(R elative units) 0 .0 0 0 .2 0 0 .4 0 0 .6 0 0.B0 o o LU. i i! I 15200 16800 18400 co( cm-1) Figure 1.2 The calculated vibronic transition I(co0) for the 0 — *m,l — >m and 2 — >m transitions of X — >B transition in I2 at 300° K. 8 = + f ‘ * u W d t + (1 / 2 Jo (3.3) where v(t) = du(t)/dt. Substituting Eq. (2.23) to Eq. (3.2) will give | ■ # I k(E)ST « (l/r)(o/fc)2 v(r*) exp{ m'bat* + (i/K)JVd^«(fp)|2 (3.4) 1 The sum over j leads in the harmonic case to the same interference effects as in the quantized expression. kST will give destructive interference unless cofea=27tmco< 2 . The preexponetial factor of Eq. (3.4) can be shown to be the Landau-Zener transition probability31, T = 1 - expilno2 / h v)~ 2no2 / hv (3.5) A comparison of kST and kSGF in the harmonic case is shown in Fig. 1.3. We will examine in more detail the nature of the thermally average rate constant. We note that kST can be written as kST = lim J 1 /t)< |f ' d t d (f)l 2>0= ton J* V f <C b(0)C b(t)>0 = (a/H)2 j 0 ° dtexp(i<o'bat)<exp[(i/h)j^dt'u(t')]>0= j c ° d t exp(i<o’ bat)C(t) (3.6) which we use the factor that for ergodic systems the time average is identical to the ensemble average and the Wiener-Khintchine theorem which relates the correlation function to the corresponding power spectrum. C(t) can be expanded by the cumulant i expansion-’- 6 j C(r) = exp[ 7(0 ] j with 7(0 = X m w dt 1 J^1 dt2 JO-1 dth l <u(h ) ... u(tj)>c (3.7) where < • • • >c is the cumulant average. The second-order approximation for y(t) is 9 > 0 . 167- ° 0 . 1 12- < / j 6 0 120 180 240 - A V ° 0 ( c m 1) Figure 1 3 A comparison of kST (dotted line) and kSQF (solid line) for a one dimensional harmonic system (co= 34cm"*, X= 1.4) L 10 Y (!) » (//ft) df'cuCfj^-Cl/fc)2 dt1 J^1 dt2 dt2 <u(tx) u(t2)>c (3.8) where < u(t) >c = < u(t) >0 <u(tl)u(t2)>c - <u(tx)u(t2)>0 = <»(fi)>o <“(f2)>0 < 3-9) | It is easy to show that y(t) can written as, y (r ) = (-1/ft2 ) dtx J*i ^ r 2 < ^ (r1 )w (r 2 )> 0 = { -\[& )^ d tx df (t-t')<u(0)u(t')>0 and the rate constant is approximated by f~ dt exp[{ifh)<&Vba>Q t + 7(0] (3.11) J —o© This second order approximation appears to be more useful than what could have been assumed in view of the somewhat arbitrary neglect of the higher terms. This might reflect the fact that the fast oscillating functions expi— i/hj* udt') contributes to the rate constant only for small values of t-t*, so that the short time approximation is a reasonable approximation. Of particular interest is the nature of y(t) for harmonic systems of equal curvature. In this case we have v a = “ A v b = ( 1/2)2 ;* < ° /? r ^ 2+AVL " = W ba - <&v be>0 = - 2 / * < 3-12> with qj = (27ij¥l)^cos{(i3j+Q p 11 where r ^ - is the mean occupation number of the oscillator at a given temperature, | lij = l/[exp(fi C D y p )— 1] i j [ we obtain j <u(0)u(t)> = c ojxj(n.j+ l/2)cos(a>jt) I which gives upon substitution in Eq. (3.10) J y(t) = ^jXj(Jij+l/2)(cos(Djt-l) (3.13) H and the semiclassical rate constant kgT is given by Eq. (3.5) with y(t) of Eq. (3.13), kST = (G/h)2j°° dt exp{(0'bJ+^jy^{njfX/2)^a^^pr-l^.) (3.14) ...... • ■ this semiclassical rate constant has an interesting similarity to the exact quantum mechanical rate constant for the harmonic problem which is given by 12,13 I cq = dt exp{ Y y tOy1 " 1/2)(ccwcoy-1)]+ i'X (ij/2)sinoijt} (3.15) * -oo it is easy to show that in the high temperature limit (noting that n-~ kBT/tk£>j when temperature is large ) Eq. (3.14) is reduced to kST = (y/fc)2 J dt exp[i((Dba+afh)t-aJcBTt2 ’ /ftL ] (3.16) ' j where a is the reorganization energy, j a = (1 /2 jZ jtH O jtf By means of the well-known relation J°° dt expii&t-at2) = (n/a)1^ 2 exp(-(£p-/4a) We get the final expression which is identical to the quantum mechanical result k = (a/h)2(nh/ kBTa)112 exp[-(7mba + a)2 / 4kBTa\ (3.17) It is interesting to compare the expressions of kST and kq in the case of of the weak coupling, i.e. small y, in the one dimensional case. kST = (a/K)2 f °° dt exp{(co^r+ca0y2/2)+Y2(7H-l/2)(coi,C D -1)} (3.18) * — 00 Using the relation f °° dx exp(-ipx+ycosx) = 2n^J^__ooh(p-n)IJy) (3.19) » — C O where Ip is the modified Bessel functionV and using vx = coat, p = (0 0 ^ /2 + co^a)/cofl, and y = (n + 1/2)^, one obtains ksj. = (C T /& )2(27t/coa) exp[-(rn-1/2)^] * ' • : " x X v = - ~ W 7 * 1/2)72] 5(y2/2+co°a/coa+N) (3.20) for small y we use l/2)y2]«[(7i+-1/2)72/2] |A 7|/|W |! (3.21) I so that kST= (a/^)2(27t/cofl) «p[-(AH-l/2)y2] [(«+l/2)y2]lN|/lAT|!5(y2/2+co^co< 2 +A0 (3.22) This can be compared to the corresponding quantum mechanical result for the t | ! harmonical case 13 kQ = (c/ft)2(27t/<oa) expi-im -l/lrf] * L n = o exp(mx3am i N{ 1 ) ] V ¥ } 6 ( c o 6 a /c o a + A 0 (3.23) for small y one obtains kQ = Cct/^)2C 2 7 c /c o a) expl-ir*!!!)^ x Z “ =0 [(7H-l/2)1 /2y2/2]A /5(co^/coa+A0/iV! (3.24) The similarity between k jj and kg in the high temperature limit is shown in Figure 1.4. As seen from the figure the semiclassical rate constant reproduces the quantum mechanical results in the weak coupling - high temperature limit. The dependence on co^ is not identical, reflecting the fact that 0 — » m Franck-Condon factor of kST are smaller by 2m from the corresponding SGF and quantum mechanical Franck-Condon factors, this disagreement is not surprising since the derivation of Cb(t) of Eq. (2.11) did not address the question of the quantization of nuclear motion. The semiclassical approach, which might not be sufficiently :reliable for systems of few degrees of freedom, is expected to give reasonable, results33 for many-dimensional systems at the high temperature limit. In this limit the nuclear motion can be treated classically. Considering the fact that evaluation of kST by molecular dynamics method is probably the most practical way of simulating ET in solutions it is important to examine the validity of kST at the high temperature limit. 1.4. Application to Multidimensional System We will consider first the simple case of two-dimensional system. The SGF treatment of this system, Fig. 1.5, is formally simple. The exponential sum Ffl of the two dimensional action is accumulated every time the trajectories returns to its initial point. Next the trajectories are allowed to cross to V b evaluating the crossing probability at the crossing point. The rate is then given by k ^ < \W F ^ bf->E a (4.1) 14 0.029- 0.022 ° 0.015- o 0.007 -252 AVb 0 o (cm '1) Figure 1.4 Comparing kq (solid line) and (dotted line) for a two-dimensional harmonic system in the small displacement (Xj = 0.6, = 0.4) high temperature limit (T = 300 K, C D j = 34 cm-1 , C O 2 = 49 cm-1). 15 2 . 600- 1 .2 0 0 - - 0 . 200- -1 . 600- — , , 1 :.. -1 .600 - 0.200 1 200 2.600 Figure 1.5 Classical trajectories on a two-dimensional potential surface. To obtain the true quantized rate constant one has to propagate trajectories on both V a and V b and locates the exact quantized energies. 16 The rate constant is not zero only when both Ffl and are quantized. The surface crossing probability W can be approximated by W = W u n & 0 W um&2> (4.2) where % is defined in Eq. (2.19). Since to evaluate quantized actions in a two dimensional system is already very complicated and time consuming, the SGF approach in many-dimensional anharmonic systems becomes impractical in this aspect. But the alternative ST approach is very easy to extended to the many dimensional system. The ST rate constant for the two-dimensional case is evaluated by Eq. (3.2) exactly as in the one-dimensional case, running trajectories on the reactant surface V fl. The nature of the trajectories that contribute to kST is shown in Fig. 1.6. The calculation of kST for the zero temperature are. compared, to the corresponding SGF results in Fig. 1.7 We obtained the same interference "effects? for the: approximately harmonic case but the intensity pattern of the ST rate constant falls down faster than that of k£Qp. However, what we try to examine here is not the quantized rate constant for a single energy but the thermally averaged rate constant. This rate constant is given where < >0 denotes an average over the initial coordinates. In a many-dimensional system which can be described as a macrocanonical ensemble, the classical trajectories will span the phase space according to the Boltzmann distribution. Thus where < >0 means an ensemble average over the initial r and p. This average can be obtained in practical calculations by running trajectories with an average kinetic energy that corresponds to the assumed temperature of the system. by (4.3) we can write kST - <k{E)sT>0 (4.4) L 2 .6 0 0 - I .200- -0 200 - t .6 0 0 - — 1 ---------------- | 1 I ____ ~ ~ - 1 . 6 0 0 - 0 . 2 0 0 1 .2 0 0 2 . 6 0 0 0 .6 0 0 - 2 . 3 3 7 0 . 7 7 9 1 .5 5 8 Time (psec) Figure 1.6 Nature of the trajectories that contribute to The upper part of the figure (a) describes a typical trajectory r(t) that move on the two- dimensional harmonic system and cross from Vfl to Vfe . (b) The lower part of the figure describes the time dependent energy gap AK(r(r)>. 18 T ? 0.245- * - ‘c < U ~ 0 . 183- a C L > w o 0 . 122- + m C o u 4 1 £ 0.061- ot 60 120 180 240 -A V ° (cm-1) Figure 1.7 Comparing kST (dotted line) and kSGF (solid line) for the two dimensional harmonic system with c O j = 34cm_1, 0^ = 49 = 0.4 and X 2 = 1.4. The calculation correspond to the zero temperature case where the trajectories on V a move with zero point energy. 19 1.5. The Indirect Method Obtaining a rate expression which is related to free energy rather than potential energy is a key point of our semiclassical approach. Yet, Eq. (2.4) cannot be evaluated directly for cases with high activation barrier, since it will take extremely long computer time to reach points where the reactant and product surfaces intersect. Thus one has to evaluate the average number of time the system reaches the transition state by an indirect method. To do this we consider the statistical probability that the system will be at the region X * in phase space where the diabatic states intersect.The relevant probability factor PCX’4) is exp(-A g* (3), which relates the probabilty of reaching the transition state to the activation free energy. As argued above, one would like to evaluate the probability that an infinitively long trajectory on the reactant surfaces will reach the hypersurface where the reactant and product surfaces intersect. This probability can be evaluated by finding a sampling potential, em, that will force the trajectory to spend most of its time near the reaction: coordination, X. The free energy change from the sampling potential zm to the true potential; E 0 , along the reaction coordinate X is where <•*•>£ indicates an average over trajectories that are propagated on Em. m Using appropriate sampling potential the system will be forced to go from the reactant state to the transition state X = X * . The schematic illustration in figure 1.8 will help the readers understand these two approaches. The choice of the sampling potential and the reaction coordinate is arbitrary. Here we will define our sampling potential as exp(rAGm^ 0(,X)V) = < ex p {-(z0{X>-EmM X )> BJ (5.1) where £ ,• is the diabatic energy of the ith resonance structure and Xj is the mapping parameter which satisfies ^ Xi = 1. The choice of the X is a little bit tricky, 20 Energy Energy (a) Reaction Coordinate Reaction Coordinate Figure 1.8 Schematic description of (a) the direct method which propagates the trajectory directly on the true potential and (b) the indirect method which force the system on different points of reaction coordinate. The mapping potentials are represented by the dotted line and the trajectories are represented by curly lines. 21 especially in a multidimensional system which has thousands of degree of freedom. Here we will define our reaction coordinate as the difference of energies of the two relevant resonance diabatic states. It might appear strange to the readers that energy is used as the reaction coordinate. But in the latter chapters we will show that the energy gap reaction coordinate is extremely powerful for locating transition state and calculating activation free energy in a large system. 22 ------------------------------------------------------------------- 1.6. References 1. Warshel,, A. and Russell, S. T. (1984) Quart. Rev. Biophys. 17, 283. 2. a) Scrocco, E. and Tomasi, J. (1973) in Topics in Current Chemistry 42, 95. b) Tapia, O., Poulain, E. and Sussman, F. (1975) Chem. Phys. Lett. 33, 65. c) Rinaldi, D. and Rivail, J.-L., (1973) Theoret. Chim. Acta. 32, 57. I 3. a) Warshel, A. Qiem.Pys.LetL (1978) 55,. 454. - b) Warshel, A. (1979) J.Phys.Chem. 30, 285. 4. Warshel, A. and Weiss, R. M.(1980) I. Am. Chem. Soc. 102,6218. 5. Pross, A. and Shaik, S. S. (1983) Acc. Chem. Res. 16, 363. 6. a) Shaik, S. S. (1986) Can. J. Chem. 64,96. b) Shaik, S. S. (1987) J. Am. Chem. Soc. 106,1227 j c) Shaik, S. S., Hiberty, P. C., Lefour, J. and Ohabessian, G. (1987) j J. Am. Chem. Soc. 109, 363 7. Chandrasekhar, J. and Smith, S. F. and Jorgensen, W. L. (1985) J. Am. Chem. Soc. 107,155. 8. Tapia, O and Lluch, J. M. (1985) J. Chem. Phys. 83, 3971 9. Weiner, S., Singh, U. C. and Kollman, P. (1985) J. Am. Chem. Soc. 107, 2219. 10. Chiles, R. A. and Rossky, P. J. (1985) J. Am. Chem. Soc. 106, 6867. 11. Kong, Y. S. and Jhon, M. S. (1986) Theor. Chim. Acta. 70,123. 12. Creshmaschi, P., Gamba, A. and Simonetta, M. (1977) J. Chem. Soc., Perkin II, 162. 13. a) Adelman, S. A. (1983) Adv. Chem. Phys. 53, 61. b) Tully, J. C. (1980) J. Chem. Phys. 73, 1975. 14. van der Zwan, G. and Hynes, J. T. (1985) J. Phys. Chem. 89, 4181. 15. Calef, D. F. and Wolynes, P. G. (1983) Phys. Chem. 87, 3400. 16. a) Warshel, A. (1982), J. Phys. Chem. 86, 2218. b) Warshel, A. and Hwang, J.-K. (1986) J. Chem. Phys. 89,4938. v c) Hwang, J.-K. and Warshel A. (1987) J. Am. Chem. Soc. 109, 715. 17. Lee, J. and Robinson, G. W. (1985) J. Am. Chem. Soc. 107, 6153. 23 18. a) Warshel, A. (1984) Proc. Natl. Acad. Sci. U.S.A. 81,444 b) Warshel, A. (1984) Pontif. Acad. Sci. Script. Var. 55, 59. c) Warshel, A., Russell, S. T. and Sussman F.(1987) Israel J. Chem. 27, 217 19. a) Warshel, A. and Russell, S. T. (1986) J. Am. Chem. Soc. 108, 6595. b) Warshel, A. (1981) Biochemistry 20, 3167. c) Warshel, A. Sussman, F.(1986) Proc. Natl. Acad. Sci. U.S.A. 83, 3806. 20. Moore, J. W. and Pearson, R. G. Kinetics and Mechanism (Wiley, New York, 1981). 21. Bronsted, J. N. and Pederson, K. (1924) Z. Phys. Chem. 108, 185. 22. Hammond, G. S. (1955) J. Am. Chem. Soc. 77, 334. 23. a) Marcus, R. A. (1968) J. Phys. Chem. 72, 891. b) Marcus, R. A. (1964) Annu. Rev. Phys. Chem. 15, 155. 24. Hush, N. S. (1958) J. Chem. Phys. 28, 296i 25. Albery, J. W. and Kreevoy, M. M. (1978) Adv. Phys. Org. -Chem.. 16, 87. 26. Kruz, J. L. and Kruz, L. C. (1985) Isr. J. Chem. 26, 339. 27. a) Keck, J. C., (1966) Adv. Chem. Phys. 13, 85. b) Anderson, J. B. (1973) J. Chem. Phys. 5 8 ,4684; c) Bennett, C. H. in Algorithms fo r Chemical Computations, edited by Christofferson, R. E. (American Chemical Societ,. Washington, D.C., 1977). d) Rosenbert, R. O., Beme, B. J. and Chandler, D., (1980) Chem. Phys. Lett. 75, 162. e) Grimmelmann, E. K., Tully, J. C. and Helfand, E. (1981) J. Chem. Phys. 74,5300. 28. Glasston, S., Laidler, K. J. and Eyring, H. The Theory of Rate Processes (McGraw-Hill, New-York, 1941). 29. Warshel, A. and Hwang, J.-K. (1985) J. Chem. Phys. 82, 1756. 30. Feynman, R. P. and Hibbs, A. R. Quantum Mechanics and Path Integrals (McGraw-Hill, New York, 1965). 31. a)Landau, L. D. (1932) Phys. Z. Sowjetunion USSR 2 ,46. b)Zener, C. (1932) Proc. R. Soc. London Ser. A 17 562. 24 c)Stueckelberg, E. C. G. (1932) Helv. Phys. Acta. 5, 369. 32. Kubo, R. (1962) J. Phys. Soc. Jpn. 1 7 ,1100. 33. a)Warshel, A. (1976) Nature, 260, 679. b)Warshel, A. and Karplus, M. (1975) Chem. Phys. Lett. 32, 11. 34. Valleau, J. P. and Torrie, G. M. In Modern Theoretical Chemistry Berne, B. J. Ed,; Plenum: New York, Vol. 5, p 137. 35. a)Pechukas, P. (1969) Phys. Rev. 1 8 1 ,174. b)Miller, W. H. and George, T. F., (1972) J. Chem. Phys. 56, 5637. 25 Chapter 2 Simulations of Electron Transfer in Water 12.1. Introduction i j Electron transfer (ET) reactions play a major role in many chemical and biological processes. Extensive studies of this class of reaction1-12 have identified key factors and provide a qualitative and sometimes quantitative understanding of the relation between the rate constants and the nature of the reacting systems. The rate constant can be expressed as k ~ C 7 2(R) exp[-Ag5 4 (a J ,a s,AG0,R)/^71 (6.1) where kb is the Boltzmann constant and R is the distance between the donor and acceptor. The key factors in ET reactions are: (i) the coupling term between the donor and acceptor, a(R); (ii) the free energy difference between the products and reactants, AG0; (iii) the reorganization energy, a s, of the solvent around the donor and acceptor; (iv) the intramolecular reorganization energy, o ^ -, in the donor and acceptor. The search for the relationship between the parameters R, A G q , <x, and a s and the activation free energy for ET reactions, Ag*, has been the objective of many studies j of ET reactions. The basic theoretical guideline for this search was developed by | Marcus1 and Levich2 and later extended by others3-5. The key role of the solvent and | the importance of using free energy rather than potential energy was recognized by Marcus who gave (for structureless donor and acceptor, held at a fixed distance in a dielectric continuum) the relation Ag * = (a s+AG0)2/4 a s (6.2) The effect on Ag* by the molecular vibrations of the donor and acceptor was formally incorporated in quantum mechanical and semiclassical treatments3^’8. 26 In the past few years we have witnessed a very intense experimental effort in | ! studying ET reactions1-12. This impressive effort has started to push the field from j I the qualitative range to the quantitative direction. Model compounds9 have been built j | with fixed distances between donors and acceptors, thus isolating the effect of R on j j the rate constant. Changing solvents and A G q for fixed R started to give more j ! i quantitative information about the exact role of a and AG0. Ultra-fast time-resolved j I experiments started to give exciting results about the relation between the rate constant ; f ' I 1 and the dynamics of the surrounding solvents. The recent experiments appear to I | confirm, in a qualitative way, the available theories. However, on the detailed j | : quantitative level one still faces major open questions. One of the most fundamental j i t question is concerned with the validity of Eq. (6.2) on the microscopic level. This | ’ i ! : equation is based on the continuum description of the solvent, which is the basis of | j j most electrostatic theories1-5. Yet the solvent polarization by nearby charges might j I not follow the macroscopic assumptions. In fact, Eq. (6.2) is recovered if one j I describes the solvent in the harmonic approximation 14, and therefore cannot be exact j | for highly anharmonic systems. I > j I Experimental attempts to explore the range of validity of Marcus’ relationship j (Eq. (6.2) have moved a long way and gave encouraging results9-12. However, at | | present it appears that the available experimental approaches cannot provide an j | accurate estimate of the contribution of the solvent to Ag*. That is, the impressive ; j study of Miller and coworkers9 presented a fit of the calculated and observed j ; dependence of the rate constant on AG0, using three adjustable parameters (solvent! i • i S ■ reorganization energy, solute reorganization energy and an effective vibrational j frequency). Such a fit cannot be used to elucidate the exact nature of Ag*; in order to j f examine a theoretical expression one needs to evaluate the relevant parameters from j independent experiments or theory. Apparently, it is hard to determine quantitatively j the value of the solvent reorganization energy for different donors and acceptors27, j j Similarly it is not simple to obtain the solute reorganization energies and extracting j I ! such parameters by fitting the calculated and observed rate constants does not! constitute a unique examination of the given rate expression. t In view of the above difficulties, it seems useful to consider computer simulation approaches as alternative and complementary ways of examining the microscopic nature of Ag*. This is done in the present work which uses a semiclassical trajectory approach in exploring the validity of Eq. (6.2) on a microscopic level. Our study evaluates the actual activation free energy of ET reactions for a model system of two benzene molecules in water. It is found that Eq. (6.2) does provide a good description of the actual microscopic results. And the experimnetal deviations from Eq. (6.2) at the inverted region, j A G q I > ot, were reproduced by considering the solute vibronic channels. 2.2. General Consideration The process of ET in solution provides one of the most instructive tests for theories of condensed phases. This process is basically a radiationless process, which provides a useful test case for exploring the relation between quantum mechanical and semiclassical treatment. Furthermore, in "outer sphere" electron transfer reaction one deals with a simple case of a weak coupling between the reactants and the products states. In such case the solvent molecules respond to the electric field from the reacting molecules at the reactant state, as long as the system does not cross to the product state. Consider first a system of a donor (D) and acceptor (A) molecules surrounded by solvent molecules, we will take the "rigid solutes" case where there is no intramolecular motion within the solute molecules. The intramolecular vibronic transitions will be considered later. The wave functions of our system are approximated by the following expression, < t£ = (7.1) where a and (3 indicates various electronic states of the solvent molecules. Neglecting charge transfer interactions between the solvent and the solute molecules we can write the effective electronic Hamiltonian for the relevant states as H (R , r) a Haa( R , 0 a<R ,r) a (R ,r ) Hw <R,r> where r and R are the coordinates of the solvent molecules and the solute approximated by semiempirical potential functions that describe the energy of the given charge forms of the solute ( e.g., D+A+ and B++A) in the given solvent. These potential functions include the interactions between the solute molecules (D and A) and induced dipoles of the solvent. The off diagonal term a is approximated here by the corresponding term in the absence of the solvent a = «j)(D)<t>(A+)| H | < J > (D +)<))(A)>. The effect of the various excited states of the solvent molecules (that can be included in principle) is neglected. 2.3. Computer Simulation of the Model System One can consider an ET reactions as a crossing between two electronic surfaces. This surface crossing process is induced by fluctuation of the solvent dipoles, and the transition occurs at points where the two surfaces intersect. The rate constant for the transition from the reactant to the product state is determined by evaluating the probability C^; this is done by propagating trajectories r(r) of the solvent molecules by means of the classical Hamilton equations of the system and evaluating the energy gap Ae(r(t)) along the trajectories. The rate constant, k, is then calculated by using the following equation, where < • • ■ >o indicates an average over initial conditions which is automatically obtained by integrating a long time trajectory in a many dimensional system. It is shown in section 1.3 that this approach converges to the exact quantum mechanical result for the harmonic test case at the high temperature limit, and in the molecules, and are the zero order "diabatic" energies of the < j ) ^ and < j)^ . The are (8 .1) anharmonic systems at the statistical limit we obtain 8,14 k = (affi)2(7zft1/kbTa ) 1/2 exp(-A g*$) (8.2) 29 where kb is the Boltzmann constant , P = 1 /kbT and Ag* is the activation free energy associated with the probability of reaching the transition state. This important point J can be best understood by considering the fact that integration of Cb gives significant contributions only when the energy gap Ae(t) approaches zero. Thus we basically evaluate the number of times the system reaches points with zero energy gap weighted | by the average probability of surface crossing, which can be approximated by the Landau-Zener probability. The number of times of reaching the points in phase space with Ae = 0 is directly related to the probability (and the corresponding free energy) associated with this region. Eq. (8.1) provides a simple yet reliable way for computer simulation of electron transfer reactions in solutions. Here we choose as a model system two benzene molecules in water. This system can be used as a model for electron transfer between aromatic molecules by adjusting the oxidation and redox potentials of the donor and acceptor. The dynamics of the solute+solvent system is simulated by the surface constrained all-atom solvent (SCAAS) model22. This model represents the solute and the first few solvation shells by all atom model of the solute and water molecules. The boundary conditions are introduced by surrounding the system by a surface of water molecules which are subjected to constraint forces as well as the internal forces of the systems. The constraint forces are introduced in order to force the average polarization, radial I distances and the average fluctuations of the surface molecules to be close to those | expected from the corresponding infinite system. The specific calculations described j in this work were done by holding the two benzene molecules 5 & apart in an environment of 60 water molecules. A typical snapshot generated by the simulation is ! described in Fig. 2.1. The parameters used in the simulation are given in Table 2.1. j I The direct evaluation of Eq. (8.1) is particularly useful for studying fast ET j reactions, which can be simulated by investing a reasonable computer time. For 1 example we present in Fig. 2.2. a direct simulation of a case where the energy difference between the product and reactant states is taken as — 40 kcal/mol. This 30 Table 2.1 Parameters used in calculations0: Type Parameter Value Water charge Qo -0.8 Oh 0.4 van der Waals 629000,-625 AC0’ BC0 450000,-600 a A and B arc the coefficients of the r~12 and r-6 terms, in units of A1 2 kcal/mol and A6 kcal/mol, respectively. The van der Waals interactions involving the water hydrogens are taken as zero. Charges are given in units of electron charge. 31 i S w - j r > * Qmi ■ < Figure 2.1 Showing snapshot of one of the configurations generated by molecular dynamics. The model system consists of two benzene molecules separated by 5 X and 60 water molecules, which was simulated by SCAAS model. 32 ---------------------------------------------------------------------------------- ■ ------------------------------------------------------------------------------- j selection leads to a small activation barrier and the two surface intersect frequently during the simulation time. As shown from Fig. 2.2, the probability increases significantly only when the surfaces intersect, thus the rate constant is basically the j average of the contributions from the different crossing processes. j It should be pointed out that the actual value of a is not given in a quantitative | way by the macroscopic theory since the effective cavity radius is not defined uniquely30. For example, our calculated value of a is about 32 kcal/mol. This value should have been about 16 kcal/mol if the induced dipoles of the system were included in the simulation. On the other hand the macroscopic estimate1 with a reasonable cavity radius of 4 X for the benzene molecule is 8 kcal/mol. A very different result is obtained if the effective cavity radius is determined consistendy by following Ref. 29 and fitting the calculated solvation energy to the Bom ’s model. With solvation energy -7 0 kcal/mol, the calculated cavity radius is 2.31 X, and with the Bom ’s formula, a = 166(l/a— l/R), where a is the solute cavity radius and R the distance between the solutes, calculated reorganization energy is about 37 kcal/mol. ^.4. Free Energy Perturbation Method While the free energy perturbation method (see section 1.5) provides a means to evaluate free energy changes, it does not give a prescription for selecting an effective mapping potential. The selection of the proper potential is not obvious in cases of charge transfer reactions, where the reaction coordinate involves a concerted polarization of many solvent molecules. The charge transfer reaction can be described by two resonance,structures, *1 = [ < W c ^ 6 ] i i '¥2 = [C6H6 C6H6~] Our selection of a mapping potential is based on the idea that the solvent polarization can be mapped most effectively by changing the solute charges. This is done in a convenient way using a sampling potential of the form ■ 1 ■ 33- -1.52 \ -9.77- -18.01- -26.25- 3.37 2.53 1.68 0.84 Time(Psec) X B l 0.BI- c O & £ cT 0.38 04 1 o C o & 0.25 O T G c fl u E -* 0.13- 0.84 1.68 2.53 Time(Psec) 3.37 Figure 2.2 The relation between the reaction probability and the time dependent energy gap. The figure describes (top) the energy gap for our model system and (bottom) the corresponding probability. 34 t m = X z l + (l-X ) £3, with X € [0,1] where Ej and £2 3X6 diabatic energies of the two resonance structures defined above. By controlling the mapping parameter X from 1 to 0 the system will be forced to move from the reactant state, £m = £j + 0 £2, to the product state, Z m = 0 £ 1 + 62- There are two choices of reaction coordinates that we can use in our model system. The most simple and obvious choice is the mapping parameter X. This is done in Fig. 2 . 3 which corresponds to AG0 = 0. This figure describes a "charging" process in which e is changed as a function of X and consequently the charge will shift gradually from the first benzene molecule to the second benzene molecule. The free energy change associated with this process can be calculated by Eq. 5 . 1 . Thus, for example, the difference between the maximum and minimum of AG(A.) gives the difference between the solvation energy of the two benzene molecules with half charge each (held 5 X apart) and the solvation energy of benzenes molecules with one electron charge. The physical meaning of Eq. ( 5 . 1 ) can be clarified by an examination of Fig. 2 . 4 . The reaction coordinate is defined as, X = £2 — £2- Fig. 2 . 4 shows the free energy profile, Ag(X), as a function of the difference of diabatic energies of the two resonance states. With our choice of reaction coordinate it is very easy to identify the transition state and hence to calculate the activation free energy from the free energy profile. 2.5. How Quadratic is the Free Energy Surface The approach described above allows one to explore the validity of Eq. ( 6 . 2 ) on a microscopic level. This is done here by repeating the same calculations presented in Fig. 2 . 5 for many values of A G q . Some of these calculations are described in Fig. 2 . 6 . Collecting Ag* for different values of A G q gives the results presented in Fig. 2 . 7 . The actual calculations are compared in this figure to the prediction of Eq. ( 6 . 2 ) ,using the correct microscopic a 14, a = <A e> - AGq (9.1) 35 A g (K cal/m ol) 1 6 .0 0 12.00 Q .0 0 4 .0 0 0.88 0.66 0 .4 4 0.22 Figure 2.3 Free energy profile of the charging process for our model system. The figure describes the change of the free energy of the system as a function of the mapping parameter X. 36 6.00 - g -22.00 - r*L-38.00 - Q0 <3 c= 0 .8 ci 4-0.2e2 -54.00 - e=Q.lgj -4-0.9egl -46.00 -12.00 22.00 56.00 Ae(Kcal/mol) Figure 2.4 The free energy profiles of the D"A and DA- state as a function of the energy gap. The free energy is calculated by running the trajectories on the different mapping potentials which gradually force the system to change from the initial state to the final state (see the text for details). Different symbols stand for data generated different mapping potentials, for example, ® , + and A for Xj = 0.1, 0.8 and 1.0 where Xj + 5^ = 1. 37 Figure 2.5 7.00 Ag -32.00 - -48.00 -12.00 22.00 56.00 A e ( K c a l / m o l ) -4.00 o g -18.00 - c O M -32.00 <P -46.00 -46.00 -12.00 22.00 56.00 Ac ( K c a l / m o l ) Free energy profile of ET in our model system for the case when AGq — 0 kcal/mol and AG0 = — 10 kcal/mol. 38 21 November 1988 13.50 o 10.00 cd C J W 6.50 * Q O 3.00 26.20 53.40 80.60 107.60 -AG0 (Kcal/m ol) Figure 2.6 Comparing the calculated free energy surface with the one predicted by Marcus’ relationship. The calculation were obtained by repeating the calculations described in Figure 2.4 and 2.5 and obtained Ag* as a function of AGq. 39 L o g k (Sec 9 .0 0 8.00 7 .0 0 6.00 0 .5 0 1.00 1 .5 0 2.00 AG (eV) Figure 2.7 Showing the effect of intramolecular vibronic channels on the rate constant. The figure uses Eq. 10.2 with the calculated Franck Condon factors of table 2.2 for a biphenyl + quinone system. 40 The agreement between the free energy dependence obtained by the microscopic simulation and the one predicted by Eq. (6.2) is quite encouraging, considering the This finding does not mean, however, that the solvent can be described within the harmonic approximation or that the continuum approximation is rigorously valid. In | our opinion this result reflects the property of high dielectric solvents which have many equivalent ways for stabilizing charges. Such solvents can give almost the same solvation energies by either large anharmonic deformation of few solvent molecules near the solute, or small harmonic motions of many molecules far from the solute24. Since both type of free energy contributions are similar in magnitude it will require more studies and careful analysis to decide whether the actual contributions of the solvent are associated with harmonic or anharmonic polarization of the solvent. 2.6. The Deviations from Marcus’ Relationship Are Associated with Vibronic Channels. As indicated from the previous section one can approximate the solvent contribution to the activation barrier by a quadratic function. Yet experimental studies of the dependence of Ag* on AG0 seem to give very large deviations from Eq. (6.2) at the so called "inverted region" where | AG q | > a . The origin of this deviation can be best understood by applying our semiclassical approach to transitions between vibronic levels, rather than electronic levels. This gives: Here < tir is the normal mode frequency of the solutes and Smm' is the Franck - fact that the present calculations are subject to convergence error of about 2 kcal/mole. where the Ae is given by (10.1) 41 Condon factor for the vibronic transition between states. The effective activation free energy associated with Eq. (10.1) can be written as ASmm’ = [AG Q - ^ r?tor(m ’ r- m r)+ a ]2/4a.-R T ln(Smm,) (10.2) • o To check the "inverted region" of the Marcus formula, we evaluate here the free energy changes with oxidation of biphenyl and reduction of quinone, which are typical components of the molecules used in the recent experiments^. The corresponding Franck-Condon factors listed in Table 2.2 and a solvent of 0.8 eV give the result presented in Fig. 2.8. This calculation reproduces the observed experimental by considering the vibronic channels of the solute molecule. 2.7. Discussion This chapter explores the free energy relationship for ET reactions by using microscopic simulation approach. The free energy perturbation method gave converging results, indicating that it can be used in exploring the solvent contribution in series of donors and acceptors. The model system used here (the benzene + water system) gives free energy surface which follows to a good approximation the Marcus’ relationship. This type of quadratic relationship is associated with the contribution of the solvent to the activation barrier with the assumption that the geometry of solute remain relative rigid. But in realistic experimental systems the solute usually is flexible and undergoes significant geometry changes and the observed barrier at the inverted region is much smaller than that predicted by Eq. (6.2). The reason for this, as demonstrated in the previous section, is the availability of vibronic channels which eliminates the need to overcome the complete solvent barrier. The agreement between the simulated free energy curve and Eq. (6.2) might give the impression that polar solvents can be described reliably by the quadratic approximation or by the linear response theory. This agreement might be accidental in the following sense; The reorganization energy associated with polarizing the first solvation shells might be similar to that obtained by small displacement of many water molecules in distance solvation shells. Thus, while in reality the solvent polarization 42 Table 22 The Franck-Condon factors of reducing quinone3: Frequencies 0-1 0-2 0-3 0-4 0-5 1699 1.2252 0.8052 0.3758 0.1394 0.0436 1491 0.0390 0.011 1120 0.0783 0.0075 467 0.9284 0.4727 0.1739 0.0515 0.0130 The Franck-Condon factors of oxidizing biphenyl*2: .Frequencies 0-1 0-2 0-3 0-4 0-5 1660 0.4612 0.1275 0.0270 0.0048 1518 0.0365 0.0032 1369 0.0634 0.0057 1002 0.0906 0.0090 775 0.0427 0.0037 351 0.1114 0.0121 0.0011 a Frequencies in cm Franck-Condon factors in dimensionless units. The reported values were obtained by QCFF/PI method26. 43 might correspond to the anharmonic saturation limit, the observed results might be similar to those expected from harmonic solvents. If, however, the actual reorganization energy is due to the harmonic polarization of the solvent, then the microscopic basis of macroscopic dielectric theories might be more valid than previously suspected. 44 2.8. References 1. (a) Marcus, R. A. (1956) J. Chem. Phys. 24, 966. (b) Marcus, R. A. (1956) J. Chem. Phys. 24, 979. (c) Marcus, R. A. (1964) Annu. Rev. Phys. Chem. 15, 155. 2. Levich, V.G. (1966) Adv. Electrochem. Electrochem. Eng. 4, 249. 3. (a) Hopfield, J.J. (1974) Proc. Nad. Acad. Sci. USA 71, 3640. (b) Jortner, J. (1976) J. Chem. Phys. 64,4860. (c) Efiima, S and Bixon, M. (1976) J. Chem. Phys. 64, 3639. (d) Bixon, M. and Jortner, J. (1982) Faraday Discuss. Chem. Soc. 74,17. 4. Dogonadze, R. R. and Kuznetsov, A. M. (1967) Elektrokhimiya 3, 1324. ; (1967) Sov. Electrochem. Engl. Transl. 3, 1189. 5. (a) Ulstrup, J. Charge Transfer Processes in Condensed Media, (Springer-Verlag, Berlin, 1979). (b) Devault, D. (1980) Quart. Rev. Biophys. 13, 387. (c) Newton, M. D. and Sutin, N. (1984) Annu. Rev. Phys. Chem. 35, 437. 6. (a) Huppert, D., Kanety, H. and Kosower, E. M. (1982) Faraday Discuss. Chem. Soc. 74, 161. (b) Kosower, E. M. (1985) J. Am. Chem. Soc. 107, 1114. 7. Churg, A.K., Weiss, R. M., Warshel, A. and Takano, T., (1983) J. Phys. Chem. 87,1683. 8. Warshel, A. (1982) J. Phys. Chem. 8 6 , 2218. 9. (a) Miller, J. R., Beitz, J.V. and Huddleston, R.K. (1984) J. Am. Chem. Soc. 106, 5057. (b) Miller, J. R., Calcaterra, L. T. and Closs, G. L. (1984) J. Am. Chem. Soc. 106, 3047. 10. Yocom, K. M., Shelton, J. B., shelton, J. R., Schroeder, W. A., Worosila, G., Isied, S. S., Bordignon, E. and Gray, H. B. (1982) Proc. Nad. Acad. Sci. USA 79, 7052. 11. Wasielewski, M. R., Niemczyk, M. P., Svec, W. A. and Pewitt, E. B. (1985) J. Am. Chem. Soc. 107, 1080. 12. Isied, S.S., Worosila, G. and Atherton, S. J. (1982) J. Am. Chem. Soc. 104,7659. 45 | 13. Churg, A. K. and Warshel, A. in Structure and Motion: Membranes, Nucleic Acid and Proteins. Eds. Clementi, E., Corongiu, G., Sarnia, M. H. and Sarma, R.H. (Adenine Press, New York, 1985) p. 361. f 14. Warshel, A. and Hwang, J.-K. (1986) J. Chem. Phys. 84, 4938. 15. Kubo R. and Toyozawa, Y. (1955) Prog. Theor. Phys. 13, 160. 16. Lin, S. H. (1968) Theoret. Chem. Acta. 10, 301. 17. Kubo, R. (1962) J. Phys. Soc. Jpn. 17, 1100. 18. Mukamel, S. (1982) J. Chem. Phys. 7 7 ,173. j 19. Warshel, A. (1984) Proc. Natl. Acad. Sci. 8 1 ,444. 20. Valleau, J. P. and Torrie, G. M. in Modern Theoretical Chemistry, edited by Beme, B. J. (Plenum, New York), Vol. 5, p. 137 21. Warshel, A. and Sussman, F. (1986) Proc. Natl. Acad. Sci. 83, 3806. 22. Warshel, A. and King, G. (1985) Chem. Phys. Lett. 121, 124. 1 23. Warshel, A. (1980) Proc. Natl. Acad. Sci. U.S.A. 77, 3105. | 24. Experimently one finds that ions are stabilized in the gas phase by small cluster of water molecules as much as they are solvented in aqueous solutions.25 However in aqueous solutions the bulk solvent gives a significant part of the total solvation energy. 25. Kebarle, P. (1977) A. Rev. Phys. Chem. 28, 445. 26. (a) Warshel, A. in Modern Theoretical Chemistry, Segal, G. Ed., Vol. 7, Plenum Press, New York, 1977. (b) Warshel, A. and Levitt, M. QCPE 247, Quantum Chem. Programs Exchange, Indiana University (1974). 27. An interesting attempt28 to evaluate reorganization energy experimentally is based on the assumption that observed threshold energy for electron emission includes the energy of the unrelaxed solvent. However the threshold energy should be somewhat below the energy of the unrelaxed excited state due to the solvent Franck-Condon factors. Furthermore, such experiment are not available for donor and acceptor at a finite separation. 28. (a) Delahay, P. and Dziedzic, A. (1984) J. Chem. Phys. 1 1 ,5381. (b) Ibid. 11, 5793(1984) 29. Warshel, A. (1981) Isreal J. Chem. 21, 341. 30. Warshel, A. and Russell, S. T. (1984) Quart. Rev. Biophys. 17, 283. 46 31. Our independent estimate o f the solvation of a benzene ion is around -70 kcal/mole. This can be used to examine AG(X) by calculating the difference AAG = AG(l/2)-AG(0) where AG(0) is equal to the solvation of a benzene ion and AG(l/2) can be determined by the powerful relationship (p.335 of Ref.30), AGJo/(R) = AGJo/(oo)-332q1 q2/R, where A G s^ 00) is the solvation energy of two benzene ions each with a charge of 1/2. With AGsol(°°) = -35 kcal/mol one obtains AG(l/2) = -52 kcal/mol which gives AAG = 18 kcal/mol, in good agreement with our direct calculations. 47 Chapter 3 Energetics and Dynamics of reaction 3.1. Introduction Early attempts2 to represent the solvent by a macroscopic cavity could not be used for quantitative study, since the macroscopic solvation energy, AG5o/, depends on very strongly the solute cavity radius, a, which is an unknown quantity. Because of the current computational capability, quantum mechanical calculations on supermolecules12 (i.e. solute plus a limited number of solvent molecules) could only be applied to the relatively small solute molecule, and the many-dimensional solvent space could not be thoroughly samples. The first attempt to obtain a quantitative approximation for quantum chemical calculations of charge transfer reactions in solution was reported in calculations of the dissociation of a water molecule in aqueous solution3,4. This calculation introduced the effect of the solvent in the diagonal matrix elements of the SCF-MO Hamiltonian of the solute. The solvent was modeled by the Surface Constrained Soft Sphere Dipole (SCSSD) model. Later studies involved replacement of the SCF-MO method by the more reliable Empirical Valence Bond (EVB) method, which was implemented in both SCSSD calculations4 and in an all-atom water model16a. These calculations included studies of proton transfer and general acid catalysis reactions in aqueous solution. Similar approaches have been implemented in calculations on enzymatic reactionsL1^ 1^. Subsequent studies7 have introduced ab initio solute potential surfaces which interact, through the calculated gas phase charge distribution, with an all-atom water model. This approach was applied to studies of Sn 2 reactions in aqueous solution. Other works8-11 which are based on ab initio solute surfaces and empirical solvent potential functions have emerged recently, indicating that the field is entering a stage of rapid progress. 48 In view of the great interest in this exciting field and the emergence of several alternative approaches, we find it both useful and instructive to apply the EVB method to the S^2 class of reactions. This class has become a popular test case for theoretical simulation of charge transfer reactions in solution. This may reflect the fact that the Sn 2 reactions are among the simplest charge transfer processes, since the changes in | solvation energy are much smaller than those of other charge transfer reactions (such ! as proton transfer reactions). The combination of EVB method and free energy perturbation method will be used to study both the energetics and dynamics of reactions in solution. The concept of the linear free energy relationships will also be addressed. This concept is one of the cornerstones of physical organic chemistry, dating back to the Bronsted relationship21 and to t(ie more quantitative (yet macroscopic) formulations of Marcus23 and Hush24. But it is still very rare in the microscopic studies of the linear free energy relationships. In fact, only a very limited number o f microscopic studies1* * have attempted to explore this problem. The use of valence bond structures5 allows one to evaluate from first principles the free energy functionals used in Marcus’ and related theories24, and to test their validity. The SN2 systems offer a very extensive base of experimental information25 that can be used as a consistency constraint on the simulations. The combination of free energy calculations and quantum chemical treatments7,16’18 is a subject of significant interest. Our approach to this problem involves a rather elaborate treatment that combines the EVB method and a free energy perturbation method. This approach might appear to some readers as an unnecessary complication considering the more recent and seemingly simpler approaches7, that used gas phase ab-initio calculations to evaluate the solute charge distribution and then calculate the solvation free energy of the frozen solute charges in the given solvent. To realize the reason for our strategy and for many of the arguments made in this paper it is instructive to consider the heterolytic bond cleavage described in ! i Fig.3.1. The figure demonstrates why a completely rigorous gas phase calculation will produce a charge distribution that will lead to a very large error in calculating the Energy ( k c a l / m o l ) 100 - 80 GO l-EA 2 0 - This en erg y will be obtained ■*" by using g a s p b a se charges ob tain ed by reliable ab initio calcu latio n s O - - 2 0 - xv i C o n siste n t calcu la tio n s I w ill g iv e this e n e rg y AG - 6 0 - AG - 80 - - I O O - 10 2 0 3.0 4.0 50 r Figure 3.1 Demonstration of the major problems associated with using the gas phase charge distribution. This figure considers the heterolytic cleavage of the X-Y bond in polar solution. This bond breaking reaction can be described by use of the covalent state Y— X, Ef^s\ and the ionic state Y+X~, E8 2{s\ where superscripts g and s indicates the gas phase and solution phase. Note that in our scale the covalent state in gas phase and in solution state were overlapped with each other. The ground state obtained by mixing E | and E | describes the bond breaking reaction in the gas phase. Using the gas phase charge distribution in the calculations will give zero solvation energy and an energy D^y (typically 90 kcal/mol) for the bond breaking reaction in solution. On the other hand, if the solvent effect is included in the solvent Hamiltonian, i.e., by mixing and £^, one obtains a solvated ground state (Y+ X~) for the bond breaking reaction in solution, which will give about 40 kcal/mol. The corresponding error is around 60 kcal/mol in the case described here. 50 J dissociation energy in solutions. This crucial point, which is not so obvious in cases of SN2 reactions, emphasizes the need for a consistent coupling between the solvent charge distribution and the solute polarization. The coupling between the solvent charge distribution (i.e., the solvent polarization) and the solute charge distribution determines the energetics of the charge j transfer reactions, and thus strongly influences the dynamics of these reactions16. J Understanding the actual nature of this coupling could be crucial in determining the dynamics of a given reaction and the corresponding reaction pathway26. Although very significant progress has been made in experimental studies of charge transfer reactions27-31, it is difficult to "invert" the experimental information in a unique way to obtain the corresponding microscopic parameters. Here again it seems essential to augment the experimental information and the available phenomenological models through microscopic simulations of these systems. The combination of EVB method and the free energy perturbation method, has been very successfully applied to proton transfer and electron transfer reactions 16(see chapter 2), and to chemical reactions in proteins18(see chapter 4). The advantage of the EVB method is particularly apparent in studies of the coupling between the solute and solvent coordinates. This coupling is entirely different for diabatic and adiabatic processes. Approaches that explicitly include the interacting electronic states provide a convenient way of exploring the two limiting cases. The simpler diabatic limit, in . j which the mixing between the relevant charge transfer states is small, was explored recently in a formal way166. The adiabatic limit, however, requires much more extensive study, since in this case a simple separation of the solute and solvent fluctuations cannot be assumed. Simulating S^2 reactions might provide a powerful theoretical benchmark for various analytical models. 51 3.2. Method This section will describe the EVB procedure for the specific case of an Sn 2 reaction. X~ + CH3- Y X -C H 3 + Y~ (11.1) This reaction can described by the following resonance structures, \|t x = X~ CH3- Y (11.2) y 2 = X -C H 3 Y~ (11.3) The matrix elements of the above 2x2 Hamiltonian H will described by analytical functions, which were obtained by fitting both to the ab initio calculations and experimental results. The gas phase Hamiltonian is expressed as 4 = « ? != « ((.!> + v ® + v ® + i/2lm 4 1) ( i.™ - * 0) 2 + i/2 5 :m 4 1) ( e L 1) -e o )2 4 = H° 22= + V® + V ® + l/2 2 m 4 !)«’” > -*o)2 +1/2X m 4 2)(9 L 2)- 9o)2 + a<2> (11.4) h 1 2 = a exp[-\L(r - rQ )] in which the designate gas phase energies, b±, b2 and r are the X -C , C -Y and X -Y distances, M(bp is the Morse potential for the ]th bond, the bm are the C— H bond lengths, and the 0m are the X -C -H , H -C -H and H -C -Y angles defined by the covalent bonding arrangement for a given resonance structure. The nonbonded interactions are described by either Ae~^r or 6-12 van der Waals potential functions (see Table 3.1). The induced dipole term is C = - 1 « 2 v U ! $ I 2 (11.5) where is the vacuum field on the m th atom in state \\fit and ym is the atomic 52 polarizability of the mth atom. No induced dipole - induced dipole interactions are difference between \|/2 and V j/j with the fragments at infinite separation, which is simply the difference in heat of formation of X- + CH3Y and Y- + CH3X. The parameters used are given in Table 3.1. The key point about the functional form of and £2 is that the experimental properties of the reactants and products at infinite separation are reproduced exactly. The gas phase potential surface used in this paper for the Cl- CH3CI system is shown in Fig. 3.2. This surface presents the best fit obtained with the two state model to the ab initio calculations7 o f the barrier and the experimental estimate of the ion-dipole complexation energy. The formal treatment of the solute-solvent system can be done by writing the wavefunction of the system as a product of the solute and solvent wavefunctions: done for the inactive electrons of the solute)'. Such a treatment takes into account the effect of the excited solvent molecules by a perturbation treatment and results in a we can write the diagonal matrix elements of the solute-solvent system as the sum of the gas-phase matrix element for the solute and a term that represents the classical interaction between the solute and the solvent. The off-diagonal term was set to its its gas-phase value to a first approximation. That is, taken into account in the gas phase system. Finally, the a(2) parameter is the energy (11.6) where the < j > are the wavefunctions of the solvent molecules and designates that the m tfl solvent molecule is in it’s )th excited state. Neglecting charge transfer interaction between the solvent and solute we can treat the solvent molecules classically (as was classical interaction between the solute charges and the solvent induced dipoles1. Thus (11.7) H n = Hv n 53 Energy (kcal I mol) 3.6 _ _ 1.2 _ ■ 1.2 _ ■3.6 _ - 6.0 - 2.0 2.0 6.0 Figure 3.2 The gas phase EVB potential' surface for the Cl~ + CHjCl — » C/C/ / 3 + C T reaction, which was obtained by fitting to ab initio calculation6. 54 Table 3.1 Parameters for the potential functions'2 ^ Bonded function parameters6 atoms D (l/2)K f o b0 a C-Cl C-H 90.0 310.6 1.8 2 .0 1.1 Nonbonded function parameters'7 atoms e* r* A c Bc M - c c r c i e i - H c r 0.1459 0.1366 3.1500 5.5866 11297.2 120.45 3.6 55 Parameters of Water-Water and Water-solute interaction atom q 4 y8 H 0.328 6.14 1.650 0.10 O -0.656 648.5 23.17 1.23 H* 0.410 5.0 2.000 O* -0.820 793.3 25.02 F~ -1.000 11000.0 2.501 c r * -1.000 20000.0 2.501 Cl" -1.000 15000.0 2.501 i - -1.000 16000.0 2.501 56; Bending potential function parameters4 * atoms d / 2 )K0 eo H-C-H 36.0 109.5 H-C-Cl 30.0 109.5 EVB offdiagonal element EVB matrix6 A r0 H12 22.0 1.0 5.0 ^Parameters used for calculations without including induced dipole. \ o aEnergies are in kcal/mol, distances in A *The bonded energy is given either by Morse potential D(e~2a^ b~b0^-2era^ b'~ b0)) or by harmonic potential ( l/ 2)K^(b-bQ)2. ‘The nonbonded interaction is given either by e*[(r*/r)12— 2 (r*/r)6] or by A exp(— \ir)— B/r6. < # The bending function is given by 1/2K0(6 -6 O )2 e The EVB off diagonal element H 12 is given by A exp(-ji.(r-r0)). ■ ^ T h e solute-water van der Waals interaction is given by AjAy/r^-Bj-By/r6. Sy is the atomic polarizabilty. 57 where Vs is the interaction energy between the solvents, V g? is the electrostatic interaction between the solute charges Q in the given resonance structure and the solvent residual charges q, Vnb is the nonbond interaction between solute and solvent atoms, and Vind is the interaction of the solvent induced dipoles with the solute charges. ^ Q , = 332 'L ij Q fi/n j Vn„ = l/2 Z tf( V ; / 4 2 - B, V r//) <1L8) V i* = - ‘ 66 I j i j j /< 0 I % I 2 / W L t W j i ' 1 r% +Y,iQ fji / ^ where the eletrostatic field % includes the effect of both the solute and the solvent charges. The constant d accounts for the attenuation of the field due to induced dipole-induced dipole interactions1. The solvent-solvent interaction is modeled by charge-charge, charge-induded dipole electrostatic interactions and 6 -1 2 van der Waals interaction. The inclusion of the solvent induced dipoles allows us to use the gas phase charge distribution in assigning the residual charges for the solvent atoms. However, for the purpose of comparison with other treatments, we will also consider a parameter set without induced dipoles. All the solvent parameters are listed in Table 3.1. In treating the solvent molecules, we introduce spherical boundary conditions using the Surface Constrained All-Atom Solvent (SCAAS) model35. This model surrounds the reactive region with a sphere of solvent molecules. The surface of this solvent sphere is constrained by additional forces that simulate the effect of the bulk solvent. With analytical functions for Ej, and / / 12 we can now obtain an analytical surface for the solute-solvent system by diagonalizing H. This gives Eg = l/2[(El + 63) - [(El - £2)2 + 4H \2] ^ (11.9) The nature of this analytical surface can best be understood from Fig. 3.3. As seen from the figure, the reactant state is favored when the solvent molecules are polarized 58 Figure 3.3 Relationship between the polarization of the solvent dipoles and the energetics of the diabatic states, and adiabatic ground state, E^, for the X~ + CH^Y — » XCH3 + Y~ reaction. When the solvent dipoles are polarized towards X, state is stabilized (a), while the solvent dipoles point toward Y, X f t 2 state is stabilized (b). 59 Energy E nergy r » ■ < (a) S o lu te C oord. (R) (b) S o lu te C oord. (R ) 60 toward X, and the product state is favored when the solvent is polarized toward Y. The solvent will act to push down the energy of either state by dipole-charge interactions thus adding a new coordinate to the potential energy surface of the reacting molecule, the solvent coordinate. This coordinate can be conveniently represented as the difference of the interaction energy of state 1 and state 2 with the solvent. In the gas phase this coordinate always has a value of 0 and the potential energy surface of the system depends only on the nuclear coordinates of the solute. The multi-dimensional nuclear and electronic solvent coordinate can be projected onto this one dimensional energy gap which represents the extent to which the solvent is polarized toward one solute charge distribution or the other. It is instructive to note that this solvent polarization strongly influences the solute charge distribution at any given solute nuclear geometry. This is because the solvent will shift the energies of the two electronic states relative to their gas phase values and this causes the mixture of states 1 and 2 in the actual electronic wavefunction to deviate from the gas-phase mixture. Thus, at any given solute geometry, the ground state in Fig. 3.3a has more negative charge on X, and that in Fig. 3.3b has more negative charge on Y. This polarization effect (which is missing in most studies) becomes extremely important in charge separation reaction since it allows the solute to respond to the solvent polarization by changing its polarization and thus influencing the solvent motion. This is also the main reason that gas phase potential energy surfaces are not useful for solution simulations - they don’t have any information about the effect of the solvent fields on the solute charge distribution. 3.3. Calculation of Free Energy With analytical potential surface on hand, it is possible to evaluate the reaction rate by running molecular dynamics trajectories directly and the rate constant can be calculated by counting the number of times the system reaches the transition state. This, however, would require an extremely long simulation time for systems with relative high activation barriers. Fortunately, one can express the relevant rate constant by36,37 61 k = Fkj-sj = F(kBT/h) exp {-A g * 0} (12.1) where 0 = (kgT)-1 , and kg is the Boltzmann constant. is the rate constant given by the standard transition state theory36 and F is the transmission factor that contains all the true dynamical information. The factor exp{-A g*0} expresses the probability that a trajectory beginning in the reactant state will reach the transition state. The problem is to evaluate the activation free energy, Ag5 6 , for the quantum mechanical ground state surface E^. An effective method that allows one to tackle this problem is provided by a combination of the EVB method and a free energy perturbation technique. The details of this method are given in Chapter 2. j | Our task is to find the probability that the system will be at the transition state, j which of course is still unknown to us. This requires an effective way of forcing the system to spend a significant amount of time in selected regions o f configuration space. Unfortunately, none of our mapping potentials coincides with the actual ground state potential. Thus we take a two step strategy starting with a convenient mapping procedure and then using the relevant raw data to construct the ground state free energy surface. The mapping potential is defined as, E m = * ? E1 + - 2< ^ ^ 1C| H 12 I j X ™ + X£=l (12.2) The system is simulated in the mapping potential and the mapping vector j Xm = (X™,X£) is systematically changed from (1,0) to (0,1) in small increments. In this | way the system can be driven from the reactant state through the transition state and | « j into the product state. Consequently, the solvent is .forced to adjust its polarization to | the changing charge distribution of the solute and the free energy changes of the | system can be calculated by the free energy perturbation technique. J The H 12 term minimizes the difference between the mapping potential and the actual ground state potential and hence enhance the convergence of our calculations. The free energy associated with changing £j to Z ^ ^ calculated by the following formulas38,3^ 8 & K -> = -d /P )/n [< ^ (-(em - e m )p>J AGfrJm&GO* - + K> = (12.3) (12.4) where < • • • >m indicates an average overthe potential surface zm, and the free energy functional AG(A-) reflects both the electrostatic solvation effects associated with the changes of the solute charges, as well as the intramolecular effects associated with changing £j to Fig. 3.4 demonstrates the dependence of AG (A ,) on X for an exchange reaction where X = Y. The figure gives both the total AG(A.) and its electrostatic components. Obtaining AG(A.) with the perturbation procedure using the mapping potentials Z m is not sufficient for evaluating the activation free energy, Ag*, which requires konwledge of the position of the transition state. The quantity AG(A) is only the free energy change associated with the switching of mapping potential. However, one still needs to find the probability of reaching the transition state for trajectories that move on E^. We will define the reaction coordinate X in terms of the energy gap £2 “ el* X = A e = e 2 - E 1 The potential of mean force as a function of X is This expression gives the probability of the system at X = £ 2 - £ 1 on the ground state E . The evaluation of Ag(X) for the exchange reaction is shown in figure 3.5. o With the same formalism, we can also obtain the probability of the system at X on the mapping potential, z {. This probability defines the free energy functionals, Ag? - , exp{-Ag(X)$} = exp{-AG(Xmm < e x p { -(E g(X )-z m( X ) m > m (12.5) exp{-AgpC )$} = exp { - AG( Am)(3} <exp {-(£,(X )-£W (X)P} >m (12.6) Figure 3.6 shows Ag;(X) profile for CT+CH ^Cl — » ClCH-i + C\~ reaction. To 63 I I O • — o O — # o 24.0— 18.0 — 6. 0 - 0.4 0.6 0.8 0.2 Figure 3.4 Free energy profile as a function of the mapping parameters X. The figure shows both the total free energy AG(k) (solid line) and the electrostatic contribution AGel(X) (dotted line). The figure demonstrates that the largest contribution to AG(X) is due to electrostatic interactions. 64 Ag(kcal/mol) 20.00 - 1 4 . 0 0 B.OO 2.00 - 1 4 4 . 0 0 —4 8 . 0 0 4 8 . 0 0 1 4 4 . 0 0 A e^kcal/m ol) Figure 3.5 Actual free energy profile for the ground-state surface as a function of energy gap Ae. The calculations are done for the Ct~+C //3C/ — » ClCH-s + C r exchange reaction, with (•) and without (A) solvent-induced dipoles. 65 Ag(kcal/mol) 80.00 60.00 40.00 20.00 F ig u r e 3 .6 a m S ' X Imb- — i -----------------,— -150.00 — 50.00 1 ------------------- 1 — 50.00 150.00 A e (kcal/mol) 12 Free energy functionals, Ag1 and Ag2, as a function of energy gap Ae. The calculations are done with (•) and without (A) solvent-induced dipoles. 66 examine simple analytical approximations for the calculated dependence of Ag* on AG0 we consider the harmonic expansion of our diabatic surface. This expansion can be obtained from the power spectra of the time dependent ez (t) as was done in our dispersed polaron model (see chapter 1) and gives £l - ^ .(T O c o ^ + A ^ ) 2 + X / ^ V ^ + A ^ 2)2 « (m)a>R(R+ARp.)2 + (?U2)(£>q(Q+Aq /2)2 ( 1 2 . 7 ) £2 - S - ^ m ^ A ^ Z + X y ^ c ^ ^ A ^ ^ + A ^ « m )< * R{R - Ar 12? + W 2 ) c o e ( 2 - A g / 2 ) 2 + AV0 where the R ’s and Q’s are the dimensionless coordinates for the solute and the solvent respectively. The effective frequencies (D q and o n R are evaluated by to = f ° ° 0 L ) ' P(co') d (£> ' in which P(co) is the normalized power spectrum of the energy J0 gap, & 2 ~ e r The R is defined as R=(b 2 — bi)((£>R\x R/Ti)l/2 where is the effective mass of the system. The reaction coordinate Q is defined in terms of the electrostatic contribution to the energy gap, A e e / . 4 - z e l Thus we have Q = (2A£e//t o g ) 1^ 2. In this way the original multidimensional system is reduced to a simple two dimensional system. The reorganization energy a is given by: a = aR + clq~ (ft/2 )c o /?A ^ + ( 7 ^ 2 ) c o g A g ( 1 2 . 8 ) In the harmonic case we can evaluate the activation potential from eq. (11.9) and eq. (12.7), obtaining A V * = (AVq + a )2/4 a - 2HU (X*) ( 1 2 . 9 ) 67 3.4. Dynamics of SN2 Reaction A consistent description of the dynamics of charge transfer reactions in polar solvents is far from being a simple task. Major problems are associated with the dependence of the solvation energy on the solute reaction coordinate, as well as the necessity to describe the reaction in terms of several electronic surfaces. The simplest case appears to be associated with electron transfer reactions, in which H12 is small. In such a case one can consider the dynamics on a single diabatic surface and evaluate the surface crossing probability at the transition state region, where Ae approaches zero. This gives the rate constant in terms of the probability of reaching the transition state and a Landau-Zener type transmission factor16*. However, in the case of large of H 12 the situation is far more complicated and one must run trajectories on the adiabatic potential surface where the solute dipole moment varies along the reaction coordinate. In such cases one has to consider the mutual solute-solvent polarization, and the corresponding formal treatments are inherently complicated (see Appendics A). To explore dynamical effects, one would like to evaluate the transmission factor F of eq. (12.1) as a function of the solvent dielectric relaxation time. Numerical simulation of F can be obtained conveniently by propagating trajectories first on em^., which will constraint the system around the transition state, then switching the potential surface to E^ and continuing the trajectories. Monitoring the downhill trajectories which start at the transition state with momenta pointing to product side and end in the product state will give a direct estimate of F37. The transmission factor can also be evaluated in terms of the reactive flux autocorrelation function54’55’60. It is instructive to consider the dynamics in two dimensions in terms of one effective solute coordinate, R and one effective solvent coordinate, Q. Such effective coordinates are conveniently defined by the R and Q of Eq.(12.7) in terms of the corresponding solvent and solute contributions to An equivalent and more familiar definition of the solvent coordinate can be obtained in terms of the macroscopic reaction field, \ R, at the solute cavity, i.e., taking Q to be proportional to ^ R we obtain 68 Aeel = (13.1) where AE el is the electrostatic contribution to £2 " e l> M -i (X 2 are the dipole moments of the solute for the corresponding diabatic states and ^y is the projection of \ R on The solvent reaction coordinate Q is defined as Q = £,y/^, where £ is a proportional constant that will be specified later. Once the coordinates are selected we need a convenient description of the corresponding effective potential. As stated above one can describe the diabatic surfaces Ej and £2 in terms of shifted harmonic potentials. It is easy to show the ground state adiabatic potential can be approximated by, Eg(Q) = V( 0 « W 2 ) © Gfi2 - j ^ = (fi/2)co2 2 2-C <2M - (13.2) where the magnitude of ground state adiabatic dipole moment, p., is given by4® M - = Hi + OivaJZKl-cos 20) (13.3) where 0 = ( l/ 2)cof- 1((£2— £1)/2/ / 12) and M -m a x = \hr^V 3Tie parameter C , is related to the solvent reorganization energy a q, ? = - ( 2 o e /Ae t W ) (13.4) Eq.13.2 can be written as V{Q) = ( m ) a QQ2 - (2 a 2 /AG)(p/jimax)2 (13.5) With the potential of Eq.(13.5) one can write the equation of motion for the solvent coordinate as41: Q = -x qQ -<»q Q ~ (co2A)(2ae /A2 )[0 (i/8 2 )2 + p ] / ^ + A(r) (13.6) With the assumption that the friction term is a Dirac function, i L _ _ _ _ _ _ _ 69 Y(0=YQW ) (13.7) Similarly, the equation of motion of R can also be approximated as R = -^ Y rR ~ (aRm dE gIdR + AR(t) (13.8) The coupled equations (13.7) and (13.8) can be expressed in terms of the equivalent Smoluchowski equation and solved under simplified conditions (appendix A). It is also interesting to attempt to solve these equations numerically by Langevin dynamics. Here, however, we prefer to explore the nature of reactive trajectories in (Q,R) space using linear response theory. The rationale is based on the following intuitive argument. The rate constant can be written as65,66 ic(t) = <X(0)h(X(t))8(X(0) -X * )> exp(-Ag*$) ~ (<X+(t)>/[)exp{-A g *(3} (13.9) where the h{X{t)) is defined as h(X(t)) = l ifX (t)> X * h{X(f)) = 0 ifXit) < X * where / is a width of the transition state. The notation <X+(r)> designates the average • • of <X{t)h(X{t))>. The quantity / / <X+(f)> is the average time for a productive trajectory to pass the transition state, and the rate constant is ~*exp(— AG * p) (13.10) The time T is not much different than the time it takes a downhill trajectory to move from the transition state to the product state. T can be approximated by the following equation, —1 (d<Q+(t)>/dt)/AQ* (13.11) 70 The time derivative of <Q+> reflects the solvent dielectric relaxation time (see below), and therefore relates the rate constant to the solvent dynamical properties. Maybe more importantly, one can use <Q+(t)> to explore the mechanistic features of the reaction coordinate outside the transition state region. The relationship between <Q(t)> and <R(t)> can be used to define reactions as either "solvent driven" or "solute driven" reactions. With Eq.13.11 in mind we can evaluate <Q(t)> by propagating many downhill trajectories. Alternatively, it can also be evaluated by using the linear response theory41. This can be done by noting that the driving force in the equation of motion for the solvent coordinate can be described in terms of the zeroth-order solvent Hamiltonian and a perturbation Hamiltonian H’ given by: H \t) = (coe ^ )(2 a Q /^m a ;cAG)p e (0 ~ FextQ = Ap0(f) (13.12) According to linear response theory41 we have <Q{t)> = -A p j^ <2 ( 0 (13.13) Approximating the autocorrelation function of Q(t) by a single exponential, <Q(0)Q(t)> = B exp{-t/xG} -|< e (0 )2 (r)> = < 2 (0 )2 (r)> = (B/xQ)exp{-t/xQ} (13.14) Describing the time dependence of the solute dipole moment with a relaxation time x^, with the approximation that p = p 1(l-exp(-t/x), we obtain the recursive approximation, < Q (t> n = -(A g p p 1/xG)|^ x p {-(r -r')/t2 }(l-ex p {-r 7x [f-1)}dO = (<Q\°°)>/XQ)(xQ-x*exp(-t/x\?-V)-(xQ-x*)exp(-t/xQ}) (13.15) where (x*)-1 = (Xq-x^1), and x^”1^ is the relaxation time obtained by evaluating the 71 time-dependent solute coordinate using Eq. 13.8 with the <Q(t)>n_i obtained in the previous iteration of Eq. 13.15. In the special case of xL = = (13.16) This interesting behavior is different from the one obtained in related studies of fluorescence line shifts50’ 51. The difference reflects the fact that in realistic adiabatic systems the development of the dipole (i takes a finite time. The main point is that the dynamical effects can be extracted from the autocorrelation function (a.c.f.) of the time dependent energy gap. The solvent contribution to the a.c.f. is characterized on a macroscopic level by the longitudinal dielectric relaxation time. Our interest, however, is to explore the validity of the macroscopic estimates on a microscopic level. 3.5. Free Energy Surfaces and Free Energy Relationships To evaluate the free energy surface of the reaction it is necessary to sample the phase space of the solvent while the solute is constrained to a sequence of different positions along the reaction coordinate. The solute (nucleophile + a methyl halide) and a surrounding cluster of 60 water molecules constituted the system we studied. (A test study using 80 water molecules yielded similar results). The potential energy surface of the entire system was evaluated using the SCAAS model55 for water, and the EVB parameters of Table 3.1 for the solute. The SCAAS model applies surface constraint forces (to maintain the correct surface polarization and solvent density) and Brownian forces to the solvent molecules at or near the surface, which for this work were the outermost 34 molecules. The frictional and stochastic components of the Brownian force are adjusted so as to provide a "thermostat" for the system. For the present simulations the average temperature was maintained at ~ 300K. Trajectories of 4 psec were collected for each of a number of points along the solute reaction coordinate (after equilibration runs of ~ 6 psec) using a time step of 0.6 fsec. Each 4 psec trajectory required 8 hours of cpu time on a VAX 11/780. The free energy change associated with moving from state m to state m ’ was taken as the average of the forward and backward processes, [5G(A,m — » \ m>)+8G(km' — » Xm)]/2. The 72 convergence error of our simulations were assessed by integrating the free energy change of the forward and backward processes and was found to be less than 1 kcal/mole. As a first test case we consider the Ct~ + CH^Cl — > CICH3 + Ct~ exchange reaction. The gas phase potential surface, which was obtained by fitting the two state gas phase EVB calculations to the corresponding ab initio barrier7 and the experimental complexation energy52 is given in Figure 3.2. The free energy for this reaction in solution was obtained with our mapping procedure, and the results are described in Figs. 3.4, 3.5, and 3.6. Fig. 3.4 describes the free energy associated with changing the mapping potential £m of Eq.. 12.2 from to The electrostatic contribution to this free energy has been calculated by retaining only the electrostatic terms of £w (the solute-solvent interaction). The resulting free energy functional reflects the change in the free energy of the solute charges during the reaction, and is also given in Fig. 3.4. The increase in the electrostatic free energy for X = (1,0 ) — » X = ( 1/2 ,1/2) is due to the spreading of the solute charge over two centers, and can be estimated without any simulation by scaling the corresponding macroscopic expression. For X = (1,0) we basically have the solvation free energy of c r an isolated Cl~ ion, AGsol = -7 0 kcal/mol. By using the Bom equation, A Gsol = (166(22/fl3)(l-l/E ) (14.1) where Q is the charge of the solute, a is the cavity radius in X and £ is the dielectric constant of solvent. The effective cavity radius a is calculated to be 2.31 X. For X = (1/2 ,1/2 ) we have charges of -0.5 on each of the Cl atoms, separated by a distance R. The corresponding energy is given by1 AGsol(R) + 332 QiQ2/R * AG ^oo) + 332 Q & fzR = A G ^(~)} (14.2) where we use £ = 80. AGsol(°°) can be estimated with = Q2 = -0.5 by Eq. 14.1 as AGW<°°> = -166(0.52+0.52)/5 = 0.5A g£7 (14.3) This gives 73 AGsol(R) = AGsol(.°oy-Z'SlQ ^n = O SA G ^ -332/4R (14.4) If we define A G (^ ) = AGsol( R ^ ) -A G ^ = - 0 . 5 A - 332/4/?* (14.5) then with /?* = 4.6A one obtains AG(A,*) = 1 7 kcal/mol as compared to our calculated value o f 16.6 kcal obtained in our simulation. The mapping free energy AG(X) is converted by Eq. 12.6 to the actual ground state free energy Ag(Ae), which is shown in Fig. 3.6 Our calculated activation barrier Agfalc = 26.0 kcal/mol is very close to the observed values Ag*bs = 26.6 kcal/mol. The difference between the barrier in solution and in the gas phase is largely due to the corresponding change in solvation free energy. Another important factor in determining the free energy surface for this reaction is the influence of the solvent polarization on the solute charge distribution. This point will be discussed below. With a reasonable model we can explore several fundamental issues. In particular, it is interesting to examine the dependence of the activation free energy on the free energy of the reaction. The EVB method allows one to address this problem in a very convenient way. For instance, changing in eq. (11.4) changes in a parametric way the free energy difference, AG0, between the reactant and product states. Evaluating the Ag* that corresponds to the various AGq’s allows one to study the underlying free energy relationship on a microscopic level16,18. Such a study is described in Fig. 3.7, which presents the calculated microscopic dependence of Ag* on AGq for an hypothetical system where all the potential parameters except are set at their values for the Cl~ + CH^Cl — » CICH3 + Ct~ system. This corresponds to a variation of the electronegativity of the atom while keeping its ionic radius constant. Our computer experiments provide a convenient data base for examination of various phenomenological models. In particular, it eliminates the need for finding a common denominator for experiments with different X ’s and Y ’s, as the change of ligands involves changes in many parameters (e.g. H12), which makes a direct comparison to a 74 ■ * 3 4 8 - 1 3 6 - U > " 00 < 24 - 12 — 28.0 -4.0 -36.0 - 68.0 AG0 (kcal/mol) 0.8 — o O ■ q 0.6 — 00 < 0.4 - 0.2 _ 28.0 -4.0 -36.0 - 68.0 AG0 (kcal/mol) Figure 3.7 (a) Relationship between the activation free energy Ag* and the reaction free energy AGq . (b) The dependence of the "linear" correlation coefficient dAg*/dAGQ on AG0. 75 macroscopic model far from being simple. For example, the correspondence analysis would require a careful comparison to the relevant gas phase calculations in order to extract the solvent contribution. As seen from the figure, the dependence of Ag* on AG0 is quite monotonic, and can be described as AAg* = 0 AAG0 (14.6) where 0 changes from around 1 for very endothermic reactions to 0.5 for AGq = 0, and then to zero for very exothermic reactions. The change of 0 from around one for the case of a product-like transition state to around zero for reactant-like transition state has been predicted by many phenomenological approaches21-25. However, the main point in our approach is not in finding the universal 0 , but examining its microscopic nature. A formal starting point can be provided by Eq.12.9. Obviously this equation is useful for estimating Ag* only if the system is nearly harmonic and if Ag* is close to AV*. In this respect, it is important to mention that our Ag^ and Ag2 surfaces are not equal to £j and The Agj and Ag2 surfaces (Fig. 3.5) are uniquely defined in terms of the mapping procedure, and represent free energies. In the range of validity of the linear response theory, where Agj and Ag2 can be approximated by harmonic functions, one can exploit the fact that Agj = Ag2 at the hypersurface where & 2 = el> and write Ag* ~ (AG0 - a )2/4 a - 2 Hl2 (14.7) The correlation parameter 0 is obtained from Eq. 14.7 by: 0 = AAg*/AAG0 = (AG0-a )/2 a Eq. (14.7) without the perturbation 2H 12 is identical to Marcus’ equation for methyl transfer reactions23,25. The problem is, however, that Eq. 14.7 is not exact, nor is it 76 (14.8) based on a fundamental microscopic principle. In fact, the Ag curves will not be harmonic for systems that do not obey the linear response approximation. Furthermore, the effective parameter a is not reproduced correctly by macroscopic estimates. Even so, the solvation of ions by polar solvents can be approximated by the linear response theory, and the Ag curves can be fitted to a scaled quadratic relationship16c. Since the error range in the Sn 2 simulation is larger than in our previous study16c of electron transfer reactions (where H12 ~ 0)* it is hard to judge I conclusively from Fig. 3.5 how quadratic (or non-quadratic) the Ag curves are; longer simulations will be required to explore this question. Nevertheless, it appears from the present calculation that Eq. 14.3 gives a reasonable approximation provided that the key parameter a is obtained from microscopic calculations, and not from phenomenological macroscopic relationships where the cavity radius is not uniquely defined. In this respect, it is important to mention that a can be estimated quite reliably by the same approach that led to Eq. 14.3. That is, the AG(X,*) for our mapping potential is given by the intersection of the Ag curves. For harmonic Ag curves one obtains: A G(k*)sol = (ot-AG0)2/4 a (14.9) using this equation and Eq. 14.7 for exchange reactions (AG0 = 0) gives a = -2A Gsol(X~) - 332/R* (14.10) This powerful relationship overcomes the need to use the ill-defined cavity radius. The role of the solvent and its contribution to the overall reorganization energy can be examined by simulating several exchange reactions, i.e., X~ + CH^X — ¥ XCH3 + X~, where the ionic radius of the negative ion is changed in a parametric way. Such a study is described in Fig. 3.8. The figure shows the j correlation of the calculated Ag* with the calculated reorganization energy for three exchange reactions. In these three reactions the same intramolecular potential is used, while the ionic radius of X- is modified to reproduce the observed solvation free A g ^(kcal/mol) 30 Cl" 25 20 15 100 80 60 40 a (kcal/m ol) Figure 3.8 Dependence of the activation barrier Ag* for the symmetric exchange reaction X~ CH^X — > XCH3 X~ on the solvation energy of X~. In order to extract the pure solvent contribution we changed only the van der Waals solute-solvent parameters of F~ and I~, and kept the intermolecular solute potential unchanged. 78 energies of F", Cl“, and I". Employing the same intramolecular potential in all three cases guarantees that the entire change in Ag* is due to solvation effects. A comparison of Fig. 3.8 to direct experimental information requires one to consider the change in Ag* between the solution and gas phase reactions, and to take into account the change in a due to the change in the X- X distance. This interesting task is not, however, within the scope of the present study. Hence, we use this calculated correlation as an "experimental" microscopic check of the validity of Eq. 14.7. That is, for exchange reactions (AGq = 0) one obtains Ag* = a/4 - 2H12 (14.11) This gives, AAg*/Aa = 1/4 The validity of this correlation is examined in Fig. 3.8. 3.6. Dynamics Effect and the Solute-solvent Coupling A consistent treatment of chemical reactions in solution requires one to' take into account the effect of solvent polarizability. In order to examine the importance of this factor, we repeat the calculations with a model that does not include induced dipoles. Both models are parameterized to give the same solvation free energy for the Cl~ ion, by adjusting the C Cl- van der Waals radius. The results of the calculations are summarized in Figs. 3.5 and 3.6. It appears that the inclusion of the solvent induced dipoles reduces the reorganization energy by about 3 kcal/mol. This is due to the stabilization of the nonequilibrium excited diabatic states (e2 for mapping with Ae>0, and for Ae<0) by the solvent induced dipoles. The small induced dipole effect might seem to justify neglecting the solvent induced dipoles. This, however, would be permissible only for reactions that involve a charge transfer over a small distance. The same type of calculations for a hypothetical charge separation reaction, CICl — > Ct~Cl+ with 5X between C/’s, gave a reduction of ~ 30 kcal/mol as a result of including the solvent induced dipoles. Another point that requires special attention is the effect of the solvent on the solute polarization. In principle, one should not use gas phase calculations to obtain the charge distribution of the solute. The effect of the solvent polarization on the solute charges is enormous in cases of charge separation reactions (e.g. Fig. 3.1), and can also be significant for Sn 2 reactions. The demonstration of the importance of this effect is left to subsequent works on charge separation reactions. For the Sn 2 reaction this point is demonstrated in the next section by simulating dynamical effects with and without the solute polarization. The emphasis in our study of dynamical effects is on the nature of the solute- solvent coupling (which is missing in most treatments60). We start by considering the nature of the rare fluctuations that lead to a reactive events. In particular, it is important to determine whether the reactive fluctuations are driven by solvent fluctuations, by solute fluctuations, or by concerted solute-solvent fluctuations. This fundamental question26 is addressed in Fig. 3.9. In the "solvent driven" limit, the reaction does not occur until the solvent reaches a polarization that stabilizes either the charge distribution of the transition state or the product state, or the activation barrier along the solute coordinate becomes small. In the "solute driven" limit the solute changes its configuration (and charge distribution) to that of the transition state, while the solvent responds by moving toward its transition state polarization. If the system follows one of these limiting cases, it is possible to consider a simplified approximation for the reaction rate. For example, in the 'class of "solvent driven" processes, one can use the conceptual model given in Fig. 3.10, where the activation barrier along the solute coordinate fluctuates as a result of the solvent fluctuations. In this way one can obtain the activation free energy for the reaction by evaluating the probability of the solute reaching its transient transition state, R *(Q ), over the solvent fluctuations160. That is exp{-A g*pj « Jexp{-AE(R0-^R*,Q(t))p}dt (15.1) In fact, this approximation might be valid over a wider range16o,18a, since eq. (53) can be thought of formally as the result of a calculation in which the solute is 80 (a) Solvent Coord. (Q) o Solvent Coord. (Q) Figure 3.9 Hypothetical free energy surfaces for (a) the solvent driven and (b) the solvent driven cases. 81 Figure 3.10 Schematic description of the role of the solvent fluctuations in controlling the activation barrier for the solute coordinate. The reaction occurs only when a solvent fluctuation stabilizes the reactant state and reduces the barrier along the solute coordinate. 82 Tim e cF1 E n e r g y © constrained to stay at R0, and the probability of reaching R * is evaluated by an | ! umbrella sampling procedure. However, one would expect to obtain different transmission factors with Eq. 15.1 (and varying degrees of convergence) in the S s different limits of the solute-solvent dynamical coupling. j Regardless of the exact relationship between the nature of the reactive trajectories j and the validity of Eq. 15.1, one would like to understand the dynamics of the solute- ■ « solvent coupling from a mechanistic point of view. Furthermore, different limits may j t i j show different behaviors in terms of isotope effects and dependence on the dielectric \ ? i j relaxation of the solvent. A convenient way to explore the complicated nature of the j reactive trajectories is provided by propagating trajectories downhill from the j ; 1 [ transition state region. The time reversal of these trajectories generates the relevant set j I of reactive trajectories. The downhill trajectories are generated by using the mapping | ■ • f | potential with X = X* to prepare the system in the transition state region, and then j j I i initiating a relaxation process by changing the potential of the system from £m* to the ! I actual ground state potential E^. Fig. 3.11 describes typical downhill trajectories. As j j seen from the figure, the situation in our model reaction is not simple. Some ! * trajectories seem to be solvent driven, and some conceited or solute driven. In all \ j | cases, however, the solvent coordinate does play a central role. Some aspects of the \ ‘ solvent fluctuations will be explored below, but a more complete examination would I require running many more downhill trajectories as well as an evaluation of the actual j ■ : free energy along the solute and the solvent coordinates. i I I , ^ While the exact value of <Q(t)> requires propagating much more downhill j : ! j I trajectories, it is instructive to examine the dynamics of the average coordinate by the j a.c.f. of Q(t) and the linear response approximation of Eq. 13.13. To do this we j i evaluate the a.c.f of Q(t)(Fig. 3.12). The simulated autocorrelation time is about 0.23 | i I picoseconds in a good agreement with the estimate44 of 0.23 picoseconds for the i , j ] solvent longitudinal dielectric relaxation time % L (a simulation of the macroscopic j dielectric relaxation time of our model is required in order to examine the exact ; j i [ j relationship between % L and T q ) . With the simulated <Q(0)Q(t)> we evaluated <Q(t)> j using the following simplified version of Eq. 13.15. I 84 | R (Solute Coord.) 3 2 . 0 0 2 4 . 0 0 1 6 . 0 0 8.00 1 5 . 0 0 2 0 . 0 0 10.00 5 . 0 0 Q (S o lv e n t Coord..) Figure 3.11 Nature of typical downhill trajectories for the CI~ + CH^ClCICHj + Ct~ exchange reaction. The system was initialized by equilibrating the system at its transition state using e = 0.5 + 0.5 6 2- After the initial equilibration the potential was changed to the actual ground state potential, and the system was allowed to freely evolve in time. The solute coordinate R is defined as: R = 0.17(|±v)1 /,2 ( Z > 2-~^i)» where p . = M C/{Mc) in unit of amu, v is taken as 800 cm-1, approximately the stretching frequency of a C— Cl bond, The solvent coordinate Q is defined as: Q = ( 2 A w h e r e v = 100 cm-1 is the effective frequency obtained from the power € / W spectrum of AZel, and AZei = is the electrostatic contribution to the energy gap between the two states used in our calculations. .85 Correlation function 0.10 — 0.20 1.20 0.30 0.90 0.60 T im e ( p s e c ) Figure 3.12 Autocorrelation function (solid line) of the energy gap for the Ct~ + CH^Cl — » CICH^ + c r reaction, which was propagated at the transition state on the mapping potential 0.58! + 0 .582* The simulated correlation function is compared with the function expC-t/T^) (— ), where % L is the macroscopic dielectric relaxation time of the solvent (tl = 0.23 psec). 86 <Q(t» = ( 2 „ « /^ ) f '< iJ ( 0 ) e ( O x H ( '- < ') > c(, ^ (15.2) | where <M -(t)>ca/c is the average of calculated solute dipole moment for the j downhill trajectories. The calculated result is in good agreement with the simulated | time dependence of <Q(t)> in Fig.3.13. More systematic studies will be required to I | examine the validity of our approach. It is clear however from the analysis presented | in eq.(39) and (40) that the time dependence of <Q(t)> is affected by the fact that the i solute dipole is formed in a finite time. I j | Although the present study explores the nature of the reactive fluctuations only in J a preliminary way, it provides a clear indication of the importance of the solute j [ polarization effects. This point is described in Fig. 3.14, which compares typical ! | downhill trajectories, with and without consistent solute polarization. In the first J trajectory we use the correct ground state surface (Ep ,while in the second case we | evaluate the solute charge by using the gas phase EVB Hamiltonian (excluding the | solute-solvent interaction term). The solute charge distribution calculated in this way j (Fig. 3.15) is determined entirely by the solute geometry, and does not reflect the polarization of the solute by the surrounding solvent. Thus, the dynamics of the j corresponding trajectories reflects the interaction of the solvent with the gas phase | charge distribution of the solute. As seen from Fig. 3.14, the trajectories of the j polarized solute move toward the relaxed ground state faster than those of the unpolarized solute. The reason is quite obvious: when the polarization effect is ! | included, the magnitude of [i increases from zero to the fully polarized value faster | than when the polarization effect is not included. The larger p . is, the stronger the j relaxation force on the reaction field. The importance of this effect is also shown in > an instructive way in Fig. 3.13 which compares the results of the approach for the <(i> j calculated with and without the effect of the solute polarization. As seen from the j figure the solvent responds in a slower way in the inconsistent treatment. ( ! | f, ( Q(Solvent Coord.) 24.00 i a .00 12.00 6.00 - / 0.36 0.24 0.12 T im e (P s e c ) Figure 3.13 Time dependence of the solvent coordinate, Q, for the downhill trajectories. The (-----) lines represent the actual simulation results while the solid curves are the estimated <Q> from the model of eq. (15.2). The solid curve represents a model with consistent solute-solvent coupling while the dotted curve is a result of an inconsistent calculation. 88 R (Solute Coord.) 20.00 21.00 14.00 7.00 0.48 0.36 0.12 T im e (P s e c ) Figure 3.14 Time-dependent solute coordinate, R, in a typical downhill trajectory for simulations with (-) and without (•••) the solvent induced dipoles polarization effect. The figure demonstrates the importance of including polarization effects. 89 Charg — 0 . 6 0 - — 0 . 7 0 - 0 . 8 0 - 0 . 9 0 - 0 . 4 0 0.24 0.36 0.12 T im e (F s e c ) Figure 3.15 Charge on the attacking chlorine atom as a function of time for the downhill trajectories. The (-) and (•••) lines indicate simulations with and without solvent induced dipoles polarization. 90 | 3.7. Results and Discussion j Computer simulation approaches can provide detailed microscopic insight into * * 1 J chemical processes in the condensed phase. However, the exploitation of such j I approaches in a consistent and systematic way is far from being trivial. In particular, it \ t } | is not clear how to incorporate both the solute and the solvent in the quantum mechanical calculations of the reacting system. Calculations that explicitly include all j | degrees of freedom of both solvent and solute are still well beyond the capabilities of j j the current generation of supercomputers. On the other hand, the use of a rigorous j S solute gas phase Hamiltonian does not provide a consistent alternative to this problem, j : Thus, it is important to develop and examine more uniform approximations that j ' ■ i . capture the main coupling between the solute and the solvent. The EVB approach j described here and in our previous studies4,16,18,19 offers a practical and consistent j j way of incorporating the solvent in the solute Hamiltonian. The effectiveness of this ! | method is demonstrated here in a combined study of the energetics and dynamics of j i ! the Sn 2 class of charge transfer reactions. 1 • I A convenient combination of the EVB method and a free energy perturbation * | technique allows us to explore fundamental aspects of free energy relationships on a j ; microscopic level. The method is used to examine the validity of Marcus type ] \ relationships for adiabatic charge transfer reactions. It is found that the contribution of j j • \ the solvent to the activation barrier can be estimated qualitatively with such!' \ relationships. This is due to the fact that the relevant free energy functionals can be j ; i ] approximated by quadratic functions, which have the same apparent dependence on ' j the solute charges as the one obtained from macroscopic dielectric theories (provided j that the reorganization energy is not evaluated using an arbitrary cavity radius, but| ;j f [ either by microscopic calculations. This does not imply, however, that thej i macroscopic theory is rigorously valid. The previously used macroscopic models of; Marcus2- ^ and Hush24 evaluated free energies by using macroscopic models thatj | constrain the electrostatic free energy to change along the reaction coordinate. Here,I I on the other hand, we use the potential energy of the constraint Hamiltonian to drive? | the system and calculate the relevant free energy from the corresponding potential| I 91? j energies. The two approaches are clearly not identical; the present one is more i ; rigorous since it actually evaluates the relevant free energies from first principles. Using the microscopic approach should give one (when the convergence error is : reduced) the ability to explore the basic derivation of the macroscopic approaches, and i ; attack outstanding problems, such as the nature of the difference between the I Marcus^ and Hush^4 points of view. The simple incorporation of the solvent field in the solute Hamiltonian is one of the key advantages of the present treatment. I j Inclusion of the solute-solvent coupling, which is missing in most other approaches, is j ! more than a pedagogical exercise. While in the S^ 2 class of reactions the effect of the | solvent on the energetics and dynamics is rather small, it becomes a huge effect in ! ] charge separation reactions (e.g. proton transfer processes). | The present work explores some dynamical aspects of an Sn2 reaction, emphasizing the solute-solvent coupling. An examination of the adiabatic limit reveals a non-linear coupling, a consistent treatment of which might be essential for reliable simulation of charge transfer dynamics. It appears that the solvent response i function is determined by both the solvent dielectric relaxation time and the solute 'polarization time. Examining this point with the linear response approximation ! provides some interesting insight into the dynamics of the solvent. However, several major issues are left open. In particular, it would be interesting to find out (by running |many more downhill trajectories) whether or not the solvent reactive fluctuations, [which are determined by the actual many-dimensional potential surfaces, can be [reliably approximated by one-dimensional dynamics on the corresponding free energy I functionals. Such an assumption, which is implicit in most of the current formal .studies, has not been derived rigorously from first principles. ' l | In concluding this chapter we would like to point out that despite the significant progress being made in computer simulation approaches, one can still benefit greatly J [from the insight provided by phenomenological models. The results of the simulation [studies can be used (in conjunction with the available experimental information) to J examine and to "scale" these phenomenological models. In this respect, we find the [EVB method particularly useful in providing a convenient connection between the : l simulation studies and the relevant free energy functionals. 92 3.8. References 1. Warshel, A. and Russell, S.T. (1984) Quart. Rev. Biophys. 17, 283. 2. a) Scrocco, E. and Tomasi, J. (1973) in Topics in Current Chemistry 42, 95. b) Tapia, O., Poulain, E. and Sussman, F. (1975) Chem. Phys. Lett. 33,65 c) Rinaldi, D. and Rivail, J.-L., (1973) Theoret. Chim. Acta. 32, 57. 3. a) Warshel, A. Chem.Phys.Lett. (1978) 55,454. b) Warshel, A. (1979) J.Phys.Chem. 30, 285. 4. Warshel, A. and Weiss, R.M.(1980) JAm.Chem.Soc. 102, 6218. 5. Pross, A. and Shaik, S.S. (1983) Acc. Chem. Res. 16, 363. 6 . a) Shaik, S.S. (1986) Can. J. Chem. 64, 96. b) Shaik, S.S. (1987) J. Am. Chem. Soc. 106,1227 c) Shaik, S.S., Hiberty, P.C., Lefour, J. and Ohabessian, G. (1987) J. Am. Chem. Soc. 109, 363. 7. Chandrasekhar, J. and Smith, S.F. and Jorgensen, W.L. (1985) J. Am. Chem. Soc. 107,155. 8 . Tapia, O and Lluch, J.M. (1985) J. Chem. Phys. 83, 3971. 9. Weiner, S., Singh, U.C. and Kollman, P. (1985) J. Am. Chem. Soc. 107, 2219. 10. Chiles, R.A. and Rossky, PJ. (1985) J. Am. Chem. Soc. 106, 6867 11. Kong, Y.S. and Jhon, M.S. (1986) Theor. Chim. Acta. 7 0 ,123 12. Creshmaschi, P., Gamba, A. and Simonetta, M. (1977) J. Chem. Soc., Perkin II, 162. 13. a) Adelman, S.A. (1983) Adv. Chem. Phys. 53, 61. b) Tully, J.C. (1980) J. Chem. Phys. 73, 1975. 14. van der Zwan, G. and Hynes, J.T. (1985) J. Phys. Chem. 89, 4181. 15. Calef, D.F. and Wolynes, P.G. (1983) Phys. Chem. 87, 3400. 16. a) Warshel, A. (1982), J. Phys. Chem. 8 6 , 2218 b) Warshel, A. and Hwang, J.-K. (1986) J. Chem. Phys. 89,4938. c) Hwang, J.-K. and Warshel A. (1987) J. Am. Chem. Soc. 109, 715. 17. Lee, J. and Robinson, G. W. (1985) J. Am. Chem. Soc. 107, 6153. 93 18. a) Warshel, A. (1984) Proc. Natl. Acad. Sci. U.S.A. 8 1 ,444 b) Warshel, A. (1984) Pontif. Acad. Sci. Script. Var. 55, 59. c) Warshel, A., Russell, S.T. and Sussman F.(1987) Israel J. Chem. 27, 217 19. a) Warshel, A. and Russell, S.T. (1986) J. Am. Chem. Soc. 108, 6595 b) Warshel, A. (1981) Biochemistry 12.1, 3167 c) Warshel, A. Sussman, (1986) F. Froc. Natl. Acad. Sci. U.S.A. 83, 3806 20. Moore, J.W. and Pearson, R.G. Kinetics and Mechanism; Wiley: New York (1981) 21. Bronsted, J.N. and Pederson, K. (1924) Z. Phys. Chem. 1 0 8 ,185. 22. Hammond, G.S. (1955) J. Am. Chem. Soc. 77, 334. 23. a) Marcus, R.A. (1968) J. Phys. Chem. 72, 891. b) Marcus, R.A. (1964) Annu. Rev. Phys. Chem. 15, 155. 24. Hush, N.S. (1958) J. Chem. Phys. 28, 296. 25. Albery, J.W. and Kreevoy, M.M. (1978) Adv. Phys. Org. Chem. 16, 87. 26. Kurz, J.L. and Kurz J. C. (1985) Isr. J. Chem. 26, 339. 27. Kosower, E.M. and Huppert, D. (1986) Ann. Rev. Phys. Chem. 37, 127 28. Werst, D.W., Londo, W.F., Smith, J.L. and Barbara, P.F.,(1985) Chem. Phys. Lett. 118, 367 29. McGuire, M. and Mclendon G. (1986) J. Phys. Chem. 90, 2549. 30. Hicks, J., Vandersall, M., Babarogic, Z. and Eisenthal, K.B. (1985) Chem. Phys. Lett. 116, 18 31. Kenny-Wallace, G.A., Quitevis, E.L. and Templeton, G. Rev. of Chem. intermediates 6 , 197 (1985) 32. Goddard HI, W.A., Dunning Jr., T.H., Hunt, W.J., and Hay, P J. (1973) Acc. Chem. Res. 6 , 368. 33. Cooper, P.I., Gerratt, J., and Raimondi, M. (1986) Nature, 323, 699 34. Epiotis, N.D. Theory o f Organic Reaction Springer-Verlag, Berlin, (1978) 35. Warshel, A and King, G. (1985) Chem. Phys. Let. 121, 124 36. Glasstone, S., Laidler, K.J. and Erving H. (1941) The Theory of Rate Processes (McGraw-Hill, New York) 37. Anderson, J.B. (1973) J. Chem. Phys. 4 159 94 38. Zwanzig, R.W. (1954) J. Chem. Phys. 2 2 ,1420. 39. Valleau, J.P. and Tom e G.M. (1972) in Modern Theoretical Chemistry, ed. Beme, B. (Plennum, New York) Vol 5 pp. 137 40. Davydov, A.S. Quantum Mechanics, Pergamon, Oxford, 1976 41. Kubo R. (1966) Rep. Prog. Theor. Phys. 29, 255. 42. Nee, T.W. and Zwanzig, R.(1970) J. Chem. Phys. 52, 6353 43. a)Szabo, A., Schulten K. and Schulten Z.(1980) J. Chem. Phys. 72, 4350 b)Schulten K., Schulten Z. and Szabo, A. (1981) J. Chem. Phys. 74, 4426. 44. The longitudinal relaxation time, zL - 0.23 psec, of water is calculated by xL = t D(2eo o + l)/(2e0+ l) where Zx = 78.6 and e0 = 1.8 are the high- frequency and low-frequency dielectric constant and xD = 8psec is the Debye relaxation time. 45. Wang, M.c. and Uhlenbeck, G.E. (1945) Rev. Mod. Phys. 1 7 ,323 46. Kramers, H.A. (1940) Physica 7, 284 47. Ovchinnikova, M. Ya., (1986) Sov. J. Chem Phys. 4 1 48. Zusman, L.D. (1980) Chem. Phys. 49, 295 49. Feynman, R.P. and Hibbs, A.R. (1965) Quantum Mechanics and Path Integrals (McGraw-Hill, New York) 50. Mazurenko, Y.T. and Bakshiev, N.G., (1970) Opt. Spectrosc. 2 8 ,490 51.Bagchi, B., Oxtoby,D.W., and Fleming, G.R. (1984) Chem. Phys. 8 6 , 257 52. Dougherty, R.C. and Roberts, J.D. (1974) Org. Mass Spectrum. 8 , 77. 53. Warshel, A. (1977) Modern Theor. Chem; Segal, G., Ed. Plenum Press; New York, Vol. 7 54. Yamamoto, (1960) J. Chem. Phys. 33, 281. 55. Chandler, D. (1978) J. Chem. Phys. 6 8 , 2959. 56. Warshel, A. and Weiss, R. M. (1979) J. Am. Chem. Soc. 101, 6131. 57. Slater, J.C. (1963) Quantum Theory o f Molecules and Solids, vol 1, (McGraw-Hill, New York) 58. Coulson, C.A. and Danielsson, U. (1954) Ark. Fys. 25, 245. 59. Bergsma, J.P., Gertner, B.J. Wilson, K.R. and Hynes, J.T. (1986) J. Chem. Phys. 8 6 , 5625. 60. Montgomery, J.A., Chandler, D. and Beme, B. J. (1979) J. Chem Phys. 7 0 ,4056. 61. Bemardi, F. and Robb, M. A. (1984) J. Am. Chem. Phys. 106, 54. 62. Shaik, S. S. (1980) J. Am. Chem. Phys. 103, 3692. 63. Huheey, J.E. (1978) Inorganic Chemistry, 2nd ed., (Harper & Row, New York). 64. Benson, S., private communication. 65. Grimmelmann, E.K., Tully, J.C., and Helfand, E. (1981) J. Chem. Phys. 74, 5300. 6 6 . Bennett, C.H. (1977) in Algorithms fo r Chemical Computation, ed. by Christofferson, R.E. (ACS, Washington, D.C.) 96 Chapter 4 Simulations of Enzymatic Reaction 4.1. Introduction Enzymes are proteins specialized to catalyze chemical reactions in biological systems. They are among the most remarkable biomolecules known because of their extraordinary specificity and catalytic power, which are far greater than those of man- made catalysts. Much of the history of biochemistry is the history of enzyme research and understanding enzyme catalysis at the molecular level is one of the fundamental problems of biochemistry. Because of the the advance of genetic engineering, the manipulation of active site of enzyme by site-directed mutagenesis has become a reality in the past few years2’ 20-22. Changing systematically specific residues can identify those groups that are crucial for the reaction of an enzyme and determine how much they contribute to the overall catalysis. But it is still not clear in molecular terms how enzymes catalyze reactions with such efficiency, precision and specificity. It is well known that the enzymes catalyse chemical reactions by reducing the activation free energies, but the main question is how enzymes can perform such good jobs, and it is the main object of our study to pinpoint the key contributions to the catalytic free energy. The computer simulation offers a way to help us to elucidate the relationship between the structure and function of an enzyme at the molecular level. The changes of binding free energies and activation free energies of the mutant enzymes are calculated by molecular dynamics simulation using the combination of EVB method and free energy perturbation technique2-* . The results of our simulations are very close to the experiment data2’20. This demonstrates the power and utility of theoretical simulation methods in studies of the effects of site-directed mutagenesis on both enzyme binding and catalysis. 97 i 4.2. Methods \ j s Our basic philosophy in simulating the enzymatic reaction is to treat the reacting I species (which will be given more definite meaning later) as a hypothetical solute and j the surrounding water-enzyme environment as a supersolvent. To view the enzyme as j r I ; an optimal solvent for the charge distribution of the transition state can help us to j ; understand the mechanism and energetics of enzymatic reactions in a more intuitive ? I way. Enzyme is a better catalyst than water in biological system because that the j dipoles and charges of the surrounding enzyme-water system can solvate the transition j ; . i i state complex better the water molecules. This concept also offers us a clue that in j j order to correlate structure-function relationship of enzyme molecules we have to j ; evaluate the solvation energies of the reacting species by its microenvironment (the \ surrounding water-enzyme system). Our strategy to calculate the binding free energy j and activation free energy can be formulated as following: j The binding free energy can be evaluated by calculating the difference between j the solvation energies of the reacting species in the unbound state and the Michaelis \ t complex state, and the activation free energy can be obtained by calculating the j j i : difference between the solvation energies of the reacting species in the ground state \ i and the transition state. Of course one should not apply the above stated strategies j 1 - j direcdy to the simulations, or there will be serious convergence problem. In order to j carry out the practical simulations, one needs a convenient method to describe the ! j ! ■ chemical reactions happening in the enzyme-substrate system. We shall use EVB ; ■ i method14,17 which was applied to the simpler charge transfer reactions in the previous | : i | chapters. EVB describes the reaction states by the simple resonance structures and the i 2 reaction can be treated as crossing from one resonance state to another. In simulations ; [ i ! ! the free energy changes from one state to another can be calculated by the free energy j \ perturbation technique. The potential energy of each resonance structure includes the j ■ ! j j bonding energy, the electrostatic and nonbonded interaction energies, as well as the i : interactions between the charges of the given resonance structure and the surrounding j environment. The actual energy of the reacting system is obtained by mixing the! . s \ energies of the relevant resonance structures. Even though the EVB is described! before, the enzymatic reactions we consider now are much more complicated. To give the readers a more clear idea about EVB, we will give an example of how to describe reactions with EVB, consider an attack of an deprotonated serine RCH2- 0 - of a serine protease, where R designates the rest of the protein, on the C = O group of a polypeptide substrate to form a tetrahedral intermediate. RH2C -0~ + (C = 0)R'R" — » RH2C - 0 - C (0~)R'R" where R ’ and R ” are the rest parts of the substrate. This reaction is the rate determining step in catalysis of serine protease and the tetrahedral intermediate27-2^ is identified as transition state complex. The location of the anionic oxygen of the tetrahedral intermediate is called the "oxyanion hole"8,26. This reaction, which is basically a "nucleophilic attack", can be described by the two resonance structures: \\f2 = [0~ C =0] y 3 = [O -C -O - ] (16.1) The structures above can be treated as hypothetical solutes embedded in surrounding water-protein supersolvent (see figure 4.1). Each of these structures can be represented quite reliably by force field potentail that represents the bonding and nonbonding interactions of the relevant fragments in their given protein-solvent environment. Thus, for example, we can represent the energies of the two resonance structures by: where the AM are the changes of the Morse potentials of the specified bonds (relative to their minimum value). The 0 ’s terms represent the angle bending energies relative 99 Figure 4.1 In the treatment of nuceophilic attack reaction inside the enzyme, 0~ C =0 — » 0 - C - 0 ~ , the reaction region O C O is represented as a hypothetical solute embedded in the surrounding water-protein supersolvent. The interactions between the solute and the solvent are obtained by the force field calculations. With such conceptual aid,the activation free energy or binding energy can be obtained by calculating the difference between the solvation energies of the reacting species in the different states. The enzymatic reaction is simulated by use of the free energy perturbation method. 100 Enzyme SCAAS water o Enzyme Suh ^Reaction re^o SCAAS water to the unstrained equilibrium (120° and 109° for sp2 and sp^ hybridization, respectively), is the out of plane force constant for deforming a planar sp2 carbon, and V gg st- ® * respectively the nonbonded and the electrostatic interaction between the fragments of the indicated resonance structure, V is the force field of the protein- I solvent atoms not included in the reaction region and represents the electrostatic J and nonbonded interaction between the reacting atoms and the surrounding protein- | solvent atoms. The a® is the energy involved in the formation of the ith resonance | structure in solution at infinite separation of its resonance fragments, a ^ ’s are | adjustable parameters which should be so calibrated such that the calculated solvation energy in aqueous solution should reproduce the corresponding experimental values7* * . This simple yet powerful calibration guarantees that the calculated energies of the given fragments will be determined at the limit of large separation, only by their interaction with the given environment. Thus, when the environment changes from water to another solvent, or from one active site to another, the only change in the energy of each fragment is the change in The actual set of energy parameters and a detailed description of the potential functions is given in Ref. 17,19. With a convenient description of the potential energy of the system one can start to explore the activation free energies for the relevant catalytic reactions. The activation free energies represent the probability of reaching the transition state on the many dimensional potential surfaces E^. Exploring this potential by a direct simulation is quite hopeless. However, one can use the EVB potentials as a convenient way of constraining the system to be at different points between the reactant and product state and evaluate the activation barrier in an indirect way. As usually, the mapping potential is defined as where the Xm are mapping parameters and the E ,- are the EVB energies. Changing the Xm changes the charge distribution of the reacting substrate and forces the protein to polarize its dipoles in response to the new charge distribution. The free energy associated with changing the Xm is given by^’10 102 AG(Xq -» kN) - X m=0 ^m+l) 5G(Xm -» Xm+1) = -(l/p)/n[<exp(-(em +1 - em OP)>m ] (16.4) where (3 = 1/k^T, is the Boltzmann factor and < • • • >m indicates that the corresponding average is obtained by running trajectories over zm> . 4.3. Free Energy Changes of Genetically Modified Subtilisin The catalysis of serine protease has been the subject of more than 30 years of intense study. It is now well known that the different enzymes in this class have the same active center even though the rest parts of the their structures are quite different. The catalytic machinery of serine proteases has been approached by at least two genetically distinct families of enzymes: the Bacillus subtilisins and the mammalian and homologous bacterial serine proteases (for example, trypsin). These two families of serine proteases show remarkably similar mechanisms of catalysis; the invariant active site includes the triad of Asp, His and Ser, (see figure 4.2) which act in a general-base-catalysis reaction, involving nucleophilic attack by Ser on the amide substrate. These residues are positioned to faciliate nucieophilic attack by the serine hydroxylate on the carbonyl of the scissile peptide bond. However, despite the wealth of information about this reaction and its detailed description in many textbooks, it is still not entirely clear what the most important catalytic factors are. In order to understand the catalysis function of the catalytic triad, it is important to correlate the structure of the active site with their catalytic action. Such a correlation should be based on evaluation of the changes in interactions between the substrate and the enzyme during the reaction. Unfortunately it is hard to measure these interactions by experimental methods. Thus it is important to try to crack the problem by theoretical approaches. The general procedure described in the previous section was applied in studying the effect at the site directed mutagenesis on the catalytic reaction of one of the serine proteases, subtilisin. The first test case is taken as the effect of the A sn -1 5 5 ^ T h r ................ ............... ■ ‘ 103 Simulations of Enzymatic Reaction 21 November 1988 Asn-155 Asn-ISS © © H H -5 ©H -§(© >H » H Asp-32 er-221 Asp-32 ► 'er-221 His-64 Proton Transfer His-64 Asn-155 t >H Asp-32 er-221 His-64 Figure 4.2 Diagram showing the mechanism of the catalytic reaction of subtilisn. The scissile peptide bond of the substrate is oriented to faciliate nucleophilic attack by Ser- 221 hydroxylate on the carbonyl carbon. The ■Ser-221 hydroxyl proton is transfered to the imidazole of His— 64. Hydrogen bond is shown in dotted line. 104 mutation on the hydrolysis of amides by subtilisin 20. The rate determining step of this reaction is the formation of the tetrahedral intermediate, Im + R'-O H + X -{C = 0)R -+Im +-H + R '-0 ~ + X -(C = 0)R - » Im+-H + R'0(C0~)(R)X (17.1) where Im and R'— OH designate His— 64 and S er-221. This rate limiting step can be described as a transfer between three resonance structures: \|/1 = [I m H -0 C = 0 ] y 2 = [Im+— H 0~ C =0] (17.2) V3 = [Im+-H 0 - C - 0 - ] where Im, O -H and C - O are, respectively, imidazol ring of H is-64, hydroxyl group of Ser— 221 and the carbonyl group of the substrate. To evaluate the changes in activation free energy in this reaction we use a mapping potential of the form12’13’15 Vm = ^ le l + ^2% + ^3e 3 (17.3) where the e ’s are the force fields of the leading valence bond structures. The V s are the mapping parameters that change Vm in a continuous way from the reactant state to the product state. The actual EVB parameters and the protein force field are from Ref. 17. The solvent molecules are simulated by SCAAS16 where the solvent is confined to a sphere of 10 X around the active site. The coordinates of the enzyme were taken from Ref. 4 and the initial coordinates of the substrate were obtained by a general molecular modelling program SMART30. Typical configurations generated during the simulation, using a mapping potential that corresponds to the transition state region, are shown in figure 4.3. As seen from the figure the native enzyme stabilizes the oxyanion configuration y 3 by an hydrogen bond from Asn— 155. This electrostatic interaction is practically missing in the mutant 105 Figure 43 Typical configurations generated by the simulation for (a) the native, (b) the Thr— 155 mutant and (c) the Ala— 155 systems. The calculations were done using a mapping potential that corresponds to the transition state region. 106 107 , N ative Subtilisin H is-64 Ser-221 Oxyanion ion Asn-155 Mutant Subtilisin (Thr-155) His-64 Ser-221 Oxyanion ion Thr-155 Mutant Subtilisin (Ala-155) His-64 Ser-221 Oxyanion ion Ala-155 enzyme. The possible importance of this stabilization was recognized in the original crystallographic study4. However, this is the first time that computer simulation gave ! ! the quantitative information about the role of Asn-155 in the catalysis of subtilisin. j This problem can be quite complicated even in the simple case of a substitution of a group that forms a hydrogen bond with the charged oxyanion. The hydrogen bond can be replaced by a bound water molecules or by an hydrophobic group, and the resulting ! changes in stabilization could be quite different in the two cases. | To convert structure to catalytic free energy, we use the free energy mapping calculations described in figure 4.4. As seen from the figure we obtain a change in catalytic free energy AAGcat of 5 .1 + 0 .5kcal/mol in an excellent agreement with the experimental value (4.7 kcal/mol)20. It is important to note that the AAG*at obtained from figure 4.4 is almost exactly reproduced by considering only the electrostatic energy contribution to AAG * (the change of electrostatic energy at the transition state is 4.8 ±0.3 kcal/mol). It is also instructive to note that the catalytic difference between the mutant and native enzyme is manifested only at the transition state, where the carbonyl oxygen is significantly charged. Trying to evaluate this difference at other points along the reaction coordinate gives incorrect results. In fact the crossover of the ' free energy surface of the native and the mutant enzyme emphasizes the importance of a consistent evaluation of the position of the transition state. Obtaining the charge distribution of the transition state from gas phase calculations might lead to a substantial error. Our second test case is the Asn— 155 — $Leu m utation^. Here we examine a variant of our procedure; that is, since the position of the transition state and its charge distribution were already evaluated in the calculations reported above ,we can map the free energy by changing the protein rather than the reacting system. Now we use a mapping potential of the form = e jV jtf* ) + Q 2V2(k *) (17.4) where the Vj and V2 are the force fields of the native and mutant systems, 110 A G ( A ,) (kcal/mol) N a tiv e Thr-155 40 30 / • 20 - 10 (1, 0, 0) (0, 1, 0) (0, 0, 1) m apping param eter X [im H-O C = 0] [im+H O ' C=Ol [lm+ -H O'-C-O] Figure 4.4 Calculated free energy profile for the mapping process that takes the system from the ground state configuration to the oxyanion configuration. The calculations are done for the native (solid line) and the Thr— 155 mutant (dotted line) systems. The transition state parameter X*1 is found by searching for the X that leads to the most frequent intersection of and £3. Ill respectively, at the transition state (k * ) of the reacting system. This mapping potential gives the free energy associated with changing Asn-155 to Leu, where the reacting substrate is kept at its transition state. Similar calculations for the ground state complete the thermodynamic cycle needed to evaluate AAG*. The same calculations for the mutation without the substrate give the change in binding free energy , AAGbind. The basic shceme of our method can be understood by the following thermodynamical process, E + S A G j E' + S ■+ ES A G M + E S * A G N cat AG- + E'S E ’ S where E, S, ES and ES* designate the native enzyme, the substrate, the Michaelis complex and the transition complex, and E ’ denotes the mutant enzyme. The change M N in binding free energy for the mutant enzyme is AAGf e = AGb -AGb , and the change in catalytic free energy is AAG* = A G ^-A G ^. These free energy changes can be obtained from experimental measurements. The last three processes, which although are not physical meaningfully, are easy to simulated by means of free energy perturbation technique. Because free energy is a state function, it can be shown that AAGbind - AG2-AGj and AAG* = AG3-AG2. Here again we reproduce quantitatively the observed experimental results figure 4.5. It is interesting to note that also in this 112 Enzymes AAGbind A A G ^at exp. calc. exp. calc. ______ Thr 4.7 5.0 ...................... Ala -0.6 ■0.2 3.7 4.1 . . . . Leu -0. 0.5 3.3 3.7 Asn 0. 0. 0. 0. AGbind ES ES E+S Figure 4.5 Calculated and observed free energies changes for Asn-155 -^Thr, A sn -l55-±L eu and A sn-155— >Ala mutations. The calculated free energies are given in brackets near the corresponding observed values. 113 case we obtain basically the same results by only considering the electrostatic contribution to V0. Evaluating eq.(2) by considering only electrostatic terms converges much faster than calculations that evaluate the full V0. This finding supports the view that electrostatic energies are the key factor in enzyme catalysis 15 and explains why our early EVB calculations that used simpler solvent models14’15,17 gave reasonable estimate of catalytic energies. Finally we took as a test case Am— 155 — *Ala mutation, this calculation figure 4.5, which was completed before we knew the actual observed result, also reproduced the experimentally observed trend, A A G ^ = 3 .7 kcal/mol (Wells & Estell, personal communication). 4.4. Concluding Remark Our early examination of the current method gave an estimated error of about 5 kcal/mol. Using longer simulation time in a study of the catalytic reaction of trypsin ^ reproduced the observed change in AG*, but this could reflect an accidental success. The present study, however, has examined additional test cases and reproduced the observed change in AG* with error of 1 kcal/mol. This might mean that we are entering a stage of quantitative structure function correlation in macromolecules. This work indicate that the changes in catalytic free energies can be estimated by calculating the electrostatic energy (solvation energy) of the charges of the relevant resonance structures. This enforces the idea that enzymes should be viewed as ’supersolvents’ for the charge distribution of their transition states and that the evaluation of the solvation energy of the reacting substrate by the enzyme-water system is the key for understanding enzyme catalysis1 Recent studies of the oxyanion hole and related systems have been discussed in terms of transition-state stabilization2. This general concept7 does not resolve several key problems. Obviously, catalysis is determined by the energy difference between the energies of the ground state and the transition state. The question is, however, how is the energy reduced? It appears that catalysis can be accomplished in a very effective 114 way by a macrocyclelike environment with preoriented dipoles10. The resulting j storage of catalytic energy is not just in the stabilization of the transition-state charge j | distribution but in the fact that the dipoles are already partially oriented in the proper J direction. Thus the system does not have to invest a large free energy in polarizing its dipoles. Apparently, the catalytic free energy is stored in the unfavorable folding j S energy needed to prepolarize the protein dipoles. For example, in our calculations, the j free energy change of changing Asn— 155 to Leu in E + S state is only 5.6 kcal/mol, while the corresponding change for the two groups in water is 12.5 kcal/mol. This I means that the Asn dipoles are in an unfavorable environment in the folded native j enzyme. Binding the charged transition state converts the unfavorable folding energy j to a favorable dipole-charge energy. Our prepolarization concept is related to the j concept o f solvent "reorganization energy", which has been established as a key factor j for reactions in polar solvent5’ 4. Obtaining a small reorganization energy is not as i j simple as it might look. The enzyme has to stabilize charges and this means that it j | must provide a very polar environment and as large (or larger) solvation energy as j water does. A large solvation energy is usually associated with a large reorganization j energy; the same solvent dipoles that stabilize the charges tend to reorient when the | charge is absent in order to minimize their unfavorable > dipole-dipole interaction, j ) Preventing reorientation of the polar environment of the protein active site requires an j j investment o f a significant folding energy. | i Obtaining semiquantitative results for the free energies of enzyme-substrate * interaction may be viewed as a "miracle" by someone who is familiar with the j i dimensionality of the problem. That is, at present it is impossible to scan all the ; relevant minima of the system at any reasonable computer time. Nevertheless, the j study reported here and a more systematic study,19 give very stable and converging | results (with error range of about 2 kcal/mol) for the absolute values or electrostatic' energy in proteins.19 The reason for this remarkable result is probably associated with! jj the fact that we are calculating the energy associated with polarizing the system from! its given set of minima and many minima give similar results. Apparendy, fors systems with many degrees of freedom one expects similar effective response to aj developing charge for quite different local minima2^ . In this way one gets the correct change in free energy even if the system is not at the lowest minima. 116 4.5. References 1. Berendsen, H. J. C., Postma, J. P. M. and van Gusteren, W. in Molecular Dynamics and Protein Structure, Ed. J. Hermans, (Polycrystal Book Service, Western Spring, Illinois, 1985), pp. 43. 2. Bryan, P. ,Pantoliano, M. W., Quill, S. G., Hsiao, H-Y and Poulos,T. (1986), Proc. Natl. Acad. Sci. 83, 3743. 3. Chandrasekhar, J., Smith, S. F. and Jorgensen, W. L. (1984) J. Am. Chem. Soc. 109, 3049. 4. Hwang, J.-K. and Warshel, A. (1987) J. Am. Chem. Soc. 109, 715. 5. Marcus, R. A. (1964) Annu. Rev. Phys. Chem. 15, 155. 6 . Mezei, M., Mehrotra, P. K. and Beveridge, D. L. (1985), J. Am. Chem. Soc. 107,2239. 7. Pauling, L. (1946) Chem. Eng. News 263, 294. 8 . Robertas, J. D., Kraut, J., Richard, A. A., and Birktoft, J. J. (1972), Biochemistry, 11,4293. 9. Valleau, J. P. and Torrie, G. M. (1977) in Modern Theoretical Chemistry (Beme, B. J., Ed.) Vol. 5, pp 169-194, Plenum, New York. 10. Warshel, A. (1978), Proc. Nad. Acad. Sci. U.S.A. 75, 5250. 11. a) Warshel, A. (1981), Acc. Chem. Res. 14, 284. b) Warshel, A. (1981), Biochemistry, 20, 3167. 12. Warshel, A. (1982), J. Phys. Chem. 8 6 , 2218. 13. a) Warshel, A. (1984), Pontif. Accad. Sci. Script. Var. 5 5 ,59. b) Warshel, A. (1984), Proc. Nad. Acad. Sci. U.S.A. 81,444. 14. Warshel, A. and Weiss, R. M. (1980), J. Am. Chem. Soc. 102, 6218. 15. Warshel, A. and Russell, S. T. (1984) Q. Rev. Biophys. 17, 283. 16. Warshel, A. and King, G. (1985), Chem. Phys. Lett. 121, 124. 17. Warshel, A. and Russell, S. T., (1986), J.Am. Chem.Soc. 108, 6569. 18. Warshel, A. and Sussman, F. (1986), Proc. Nad. Acad. Sci. U.S.A. 83, 3806. 19. Warshel, A., Sussman, F. and King, G. (1986), Biochemistry, 25, 8368. 20. Wells, J. A., Cunningham, B. C., Graycar, T. P. and Estell, D. A. (1986), Proc. Trans. R. Soc.(London) A 317,415. 21. Wilkinson, A. J., Fersht, A. R., Blow, D. M., Carter, C. and Winter, G. (1984) Nature, 3 0 7 ,187. 117 22. Craik, C. S., Largman, C., Fletcher, T., Roczniak, S., Bar, P. J., Fletterick, R. and Rutter, W. J. (1985) Science, 228, 291. 23. Warshel, A. and Hwang, J.-K. (1986) J. Chem. Phys. 84,4938. 24. Warshel, A., Russell, S. T. and Weiss, R. M. (1982) in Biomimetric Chemistry and Transition State analogs, eds. Green, B. S., Ashani, V. and Chipman, D. (E;sevier, Amsterdam), pp. 267. 25. Hwang, J.-K. and Warshel A. (1987) Biochemistry 26, 2669. 26. Henderson, R. (1970) J. Mol. Biol. 54, 341. 27. Mathews, B. W., Sigler, P. B., Henderson, R. and Blow, D. M. (1967) Nature 214, 652. 28. Shotton, D. M. and Watson, H. C. (1970) Nature 225, 811. 29. Stroud, R. M., Kay, L. M. and Dickerson, R. E. (1974) J. Mol. Biol. 83, 185. 30. SMART is a general purpose molecular modelling software by J.-K. Hwang. 118 Appendix A The coupled equation of motion for the effective solute solvent coordinate depends of course on the definition of these coordinates. A convenient definition is provided by the Q and R of chapter 4. With this definition and the potential for Q in Eq. 13.5 we can write mQQ' = -j^ m QyQ(t-z)a(x)dx - dV/dQ' + A \t) (A l) where A ’(t) is a random force, and Q = (m^co^/h Approximating yq by a Dirac function y = y^S(t) gives Q = -y 2 < 2 - ( C OgmdV/dQ + A(t) (A2) where A = (cOg/h m g)1/ 2 A ’ and the mass m^ can be found by the equipartition relationship mQ = (A3) • * The friction term can be obtained from the a.c.f. of Q by1 yQ = (kBT/mQ)(j~<C? (0 )^ {t)>dt)~l = (kBTcoe /h ) (J ~ < £ (0)0 {t)>dt)~l (A4) and A(t) is related to yq by the fluctuation dissipation theorem1. We can also use the properties of Brownian Harmonic Oscillator4 to relate yq to the dielectric relaxation time of the solvent. That is, the autocorrelation time of the reaction field, % R y and therefore Q is given by the longitudinal relaxation time, xL, of the solvent. Although ________________________________________________________________________________________119 the microscopic nature of the first few solvation shells is likely to give somewhat different relaxation time, we can still use the approximation: X L ® = J “ < C (O X 2 (0 > * /(fi(0 )fi(0 )) = s in ^ t) = 2 2 2 where c O j = 00q - (Yq/4). The friction term can now be expressed as, V q = c o qtq Thus we have: Q = -to2 q * qQ ~ co q Q - 2(toQ/h)(iiaQ /}LrmxAQ)[(d\iJdQ)Q + |x] + A(f) (A6 ) Similarly, we can write the approximated equation of motion of R as R = -y J t-(( 0 R/h)dEg/dR + AR(t) (A7) where we assume that the main friction and random force contribution to R comes through the dependence of on Q. It should be instructive to compare the numerical results of eqs (A6) and (A7) to those of the actual many-dimensional simulations. It is also interesting to find approximated analytical solutions for these equations. For example, one can consider the Smoluchowski equivalent of eq. (A6 ) (with the 2 diffusion constant Dq = kgT/(m gC£>Q T q))- 3 r j ^ ^1/ |F (Q ,0 = (mGC O ^ )A (D (e)(^ + ^ ) )P (Q ,t) (A8) where P(Q,t) is the probability o f finding the system using the corresponding first passage approximation[3], which reduces to Kramer’s derivation5. In the high barrier case one obtains2,4,6 kQ(R) - 120 = (G > g /co 2)(27ttQ r 1 ejcp(-Ag,i(/?)3) (A9) where cO q is the curvature of Ag at the transition state and is the average diffusion constant. The problem is, however, that in contrast to the simple case of diabatic electron transfer, we deal here with a strong coupling of Q and R as well as with the complicated form of fi. A rather rough approximation can be obtained for the solvent- driven case by taking k as the value kg(R) at the ground state value of R. References 1. Kubo R. (1966) Rep. Prog. Theor. Phys. 29,255. 2. Calef, D.F. and Wolynes, P.G. (1983) Phys. Chem. 87, 3400. 3. a)Szabo, A., Schulten K. and Schulten Z.(1980) J. Chem. Phys. 72,4350 b)Schulten K., Schulten Z. and Szabo, A. (1981) J. Chem. Phys. 74, 4426. 4. Kramers, H.A. (1940) Physica 7, 284 5. Ovchinnikova, M. Ya., (1986) Sov. J. Chem Phys. 4 1 6 . Zusman, L.D. (1980) Chem. Phys. 49, 295 121,

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Computer modeling and free energy calculatons of biological processes

PDF

Computer simulation of polar solvents with the SCAAS model

PDF

Calculations of electrostatic interactions in proteins

PDF

Chemically modified electrodes and microelectrodes

PDF

Crystallographic stydies on Zn(2), Mn(2) and Pt(2) complexes of phosphonates

PDF

Computer simulation studies of charge transfer through biological and artificial membrane channels

PDF

Consistent evaluation of binding free energies and study of the role of electrostatic effects in the stabilization of protein complexes

PDF

A highly efficient and general method for the chemical synthesis of beta,gamma-fluoro- and difluoromethylene analogs of purine and pyrimidine 5'-ribo- and deoxyribonucleotides

PDF

A photophysical study of pyrene adsorbed onto silicas of variable surface area and porosity

PDF

Chemiluminescence of 9,10-dihydro-10-methylacridylic imidazolide and 9-alkylidene-10-methylacridans

PDF

Computer simulation of polymers and membranes

PDF

A study of the double sulphates of some rare-earth elements with sodium, potassium, ammonium and thallium

PDF

A general theory of auditor independence and empirical test of independence in fact.

PDF

Application of nonlinear spectroscopic techniques to the studies of molecular electronic states and reaction dynamics.

PDF

Chiroptical spectroscopy of nitrogenase

PDF

Carbocations, carboxonium dications and their role in superacid catalyzed reactions

Chemical reactions of C<sub>6</sub>-hydrocarbons on Au(111) surfaces.

PDF

Chemical reactions of C6-hydrocarbons on Au(111) surfaces.

PDF

Computer simulations of polymeric systems

PDF

13C NMR spectroscopic study of substituted diamondoid hydrocarbons and their carbocations (dications)

PDF

Computational investigation of the catalytic aspartic acids of HIV-1 protease

Asset Metadata

Creator Hwang, Jenn-Kang (author)

Core Title Computer simulation of chemical reactions in aqueous solutions and biological systems

Degree Doctor of Philosophy

Degree Program Chemistry (Chemical Physics)

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag chemistry, physical,OAI-PMH Harvest

Language English

Contributor Digitized by ProQuest (provenance)

Advisor Warshel, Arieh (committee chair), Cheng, Hsueh-Cheng (committee member), Singer, Lawrence (committee member)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c17-650716

Unique identifier UC11344165

Identifier DP21970.pdf (filename),usctheses-c17-650716 (legacy record id)

Legacy Identifier DP21970.pdf

Dmrecord 650716

Document Type Dissertation

Rights Hwang, Jenn-Kang

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA