- Email: [email protected]
Activation free energy for proton transfer in solution
Chemical Physics 180 (1994) 181-189 North-Holland
Activation free energy for proton transfer in solution Daniel Laria l, Giovanni
Ciccotti 2
Centre Europken de Calcul Atomique et Mokdaire,
Bbtiment 506, Universitk de Paris-Sud, 91405 Orsay, France
Mauro Ferrario Dipartimento di Fisica, Universitb di Messina, CP 50, 98166 S. Agata di Messina, Italy
and Raymond
Kapral
Chemical Physics Theory Group, Department of Chemistry Universityof Toronto, Toronto, Ontario, Canada M5S IA1
Received 18 May 1993
Path integral molecular dynamics methods are employed to compute the free energy for proton transfer reactions for strongly hydrogen bonded systems in a polar solvent. The free energy profile is calculated using several different techniques, including: integration of the mean force acting on the proton path with its centroid constrained at different values, the integral form of the free energy calculation in the constrained-reaction-coordinate-dynamics ensemble and direct simulation of the unconstrained dynamics. The results show that estimates of the free energy barrier obtained by harmonic extrapolation are likely to be in error. Both quantum and classical results for the free energy are obtained and compared with simulations using adiabatic quantum dynamics. Comparison of the quantum and classical results show that there are quantum corrections to the solvent contributions to the free energy.
1. Introduction
Proton transfer is a process that occurs in many chemical and biological reactions. The calculation of such transfer processes is difficult because quantum effects may be important and the reaction coordinate that is used to specify the progress of the reaction is often taken to involve the solvent degrees of freedom. In general, proton transfer involves the description of the dynamics of a quantum particle in a many-body solvent. In a recent paper [ 1] we presented a molecular dynamics calculation of the adiabatic transfer of a prol
On leave from Departamento Quimica de Reactores, Comision National de Energia Atomica, ( 1429) Buenos Aires, Argentina. * Permanent address: Dipartimento di Fisica, Universita “La Sapienza”, Piazzale Aldo Moro, 2, 00185 Rome, Italy. 0301-0104/94/$7.00
ton between two negatively charged ions immersed in a classical dipolar solvent. In this adiabatic limit the proton charge density can instantaneously adapt to the configuration of the classical particles; thus, the implementation of the dynamics is conceptually simple (but computationally demanding) and involves the solution of the Schrodinger equation for the proton, which is a function of the instantaneous solvent configuration, and Newton’s equations of motion for the solvent, which depend on the proton charge density. The proton transfer is governed by local solvent fluctuations in the vicinity of the proton-ion complex and the progress of the reaction can be monitored by choosing the solvent polarization, AE ( {R} 1, as a reaction coordinate [ 21. Consequently, the quantum problem reduces to a classical one with a complicated, many-body reaction coordinate. As a result of the preferential solvation of the complex with the proton charge localized near one of the two ions, the fluc-
@ 1994-Elsevier Science B.V. All rights reserved.
SSDI 0301-0104(93)00002-R
182
D. Lara et al. / Chemical Physics I80 (1994) 181-189
tuations of the solvent that induce the proton transfer are rare events. The asymmetry in the charge distribution of the complex leads to a larger dipole moment of the complex and to stronger solvation energy. The passage of the proton from reactant to product states requires considerable energy and involves substantial modification in the local structure of the solvent. The magnitude of the free energy barrier relative to the thermal energy determines the order of magnitude of the proton transfer rate. In our previous study, the free energy barrier was estimated from an extrapolation to the transition state of the logarithm of the probability density of the reaction coordinate computed near the reactant and product states. The free energy profile in the stable state regions was found to be approximately harmonic, in accord with Marcus’s theory [2] for charge transfer processes. In this article, we use Feynman’s [ 31 path integral formulation of quantum mechanics to treat the proton degrees of freedom, and choose the centroid as the reaction coordinate [4-61. Our calculations are carried out using constrained molecular dynamics methods. This study has several goals. We consider various constrained molecular dynamics methods for the quantum path integral calculation of the free energy. The calculations reported here using the centroid reaction coordinate are compared with our earlier results using AE as the reaction coordinate and adiabatic dynamics. We discuss features associated with the three-dimensional character of our reaction coordinate. Finally, we compare the quantum calculations with classical simulations on the same model. The present calculations are carried out in the same spirit as other recent calculations of the free energy of proton transfer; for example, we mention work by Lobaugh and Voth [7] for a full three-dimensional treatment of the centroid reaction coordinate for proton transfer in water, the analog of our calculations for a strongly hydrogen bonded system, and by Borgis et al. [8] for proton transfer in a dipolar solvent like the one considered here.
2. Activation free energy for proton transfer The methodology for the evaluation of the rate constants of simple classical activated processes is well
established and the computer simulation of such processes does not present a major challenge. Normally, the calculation proceeds in two steps: first, the free energy along the reaction coordinate is calculated. This free energy calculation is used to locate the bottleneck or transition state configuration from which the transition state theory estimate of the rate constant may be obtained [ 91. Second, corrections to transition state theory are obtained by evaluating the reactive flux correlation function k(t) [91. The transmission coefficient can be extracted from the plateau value of k (f ) provided the characteristic time scale for the reaction is much longer than that of any other microscopic relaxation times in the system. Only in this circumstance will a phenomenological description of the reactive process be meaningful. In this article we focus on the first of these steps, the calculation of the free energy for a quantum reaction coordinate. The calculations are considerably more complicated when quantum effects cannot be ignored. It is often convenient to cast the description of such quantum rate processes in terms of Feynman’s path integral formulation of quantum mechanics [ 31. Let 5 [r (t ) ] be the reaction coordinate, which is a functional of the coordinates, r (t 1, of the proton path at the imaginary time t, and 5 be its numerical value. The free energy II’ (5) is related to the probability density of < [r (t ) 1,
where (. . .) represents a canonical statistical average over all degrees of freedom of the system and /3-’ is Boltzmann’s constant times the temperature. Physically, I+‘(c) represents the reversible work required to move the reaction coordinate to the point <. If tt is the value of r at the transition state, then II’ s Wt can be used to compute the transition state theory approximation to the quantum rate constant. A discussion of the full rate constant and its breakup into transition state theory and correction terms can be found in refs. [ 5 ] and [ 6 1. In the case where quantum fluctuations can be neglected in all particles with the exception of the proton, the canonical average in eq. ( 1) can be expressed as
D. Lara et al. / Chemical Physics 180 (1994) 181-189
2.2. X
J.. J 9
wbl
s(ti[r(t)l
x exp(-S[rUL{R)l/h)
-0
,
(2)
where Vci({I?} ) represents the potential energy of all the classical particles with coordinates {R}, Q is the partition function of the system and D [ r (t ) ] represents the path integral. For a given configuration of the classical particles, the weight functional for the proton path is governed by the action
Simulation
+ J dt{~~pl~(t)12 &[r(f);{R}l},
method
Simulation techniques can be constructed to eMiciently sample fluctuations of the proton path. The proton path maps onto a ring polymer containing P interacting sites .with harmonic interactions between nearest neighbors [ lo]. This isomorphism can be exploited to construct a simulation algorithm. For a finite value of P, the problem reduces to the consideration of a classical system with an effective potential given by Kn = V,I({RI)
=
183
Pm, ’ + 2fph_2 c,=, b.1 -rl+l)* (41
(31
0
where mp is the mass of the proton, VP[r(t); {R} ] is the potential energy of the proton path in the field of the classical particles and h is Plan&s constant divided by 2~.
2. I. Model system The system studied is identical to that in our previous investigation [ 11. It consists of a proton and a pair of A- ions solvated by N = 342 dipolar molecules. The temperature of the system is 250 K and the length of the simulation box is 30.52 A. It is convenient to define a coordinate system with origin at the center of the simulation box. The two A- ions are fixed at Ri,z = (0,O) f 1.3 A. The solvent molecules are composed of two interaction sites with partial charges of opposite sign zu = f0.52e and a fixed bond length of 2 A. Site-site interactions arise from 6-12 LennardJones and Coulomb forces. The proton-ion potential was constructed to model strongly hydrogen bonded systems: the proton is confined by the potential to lie close to the ion-ion internuclear axis and there is a very small barrier in the potential (w 0.2 kcal/mol) along this axis (cf. fig. 1 of ref. [ 1] ) . Finally, the proton interacts with the solvent via simple Coulombic forces. The details of the potential energy parameters are given in ref. [ 11.
where r, are the coordinates of the polymer beads. The accuracy improves as P increases and the isomorphism becomes exact in the limit P --+ CO. The convergence in the quantities of interest determines the choice of P. We found no noticeable systematic changes for P > 20; this number has a magnitude similar to that used in previous path integral simulations of proton-like particles [ 11,12 1. Several methods exist for the calculation of free energy profiles for quantum activated processes. Lobaugh et al. [ 71 have used umbrella-sampling Monte Carlo [ 131 to treat proton transfer reactions in aqueous media. Here, we employ constrained molecular dynamics to compute the free energy of the proton transfer process. In particular, constanttemperature molecular dynamics trajectories were generated by the following Hamiltonian: P
H = c 1=1
N
fm*;f
+ c fM$: JCI
+ V, + Vk,
(5)
where m* is the fictitious mass of a polymer bead and Mj are the masses of the classical particles. The equations of motion were integrated using the Verlet algorithm [ 141 with a time step of 10pL4 s. Geomet-! rical constraints from intramolecular bonds and constraints acting on the reaction coordinate were treated with the SHAKE algorithm [ 151. Sampling from a canonical distribution was carried out using NosCHoover dynamics [ 161. We found it quite difficult to obtain proper thermalization of the proton-polymer
184
D. Lara et al. / Chemical Physics 180 (1994) 181-189
degrees of freedom even after long thermalization periods. Similar difftculties were pointed out in the muonium simulation described in ref. [ lo]. Although the total kinetic energy of the system was well controlled by the Nose-Hoover thermostat, we noticed a systematic drift towards higher temperatures in the kinetic energy associated with the degrees of freedom of the polymer. For this reason, a second thermostat acting exclusively on the proton-polymer degrees of freedom had to be included. Only in this case could we obtain reasonable thermal equilibrium between the proton and the solvent. The inclusion of a second thermostat introduces a net acceleration of the center of mass of the system; however, we verified that the center of mass oscillates and that the kinetic energy associated with it remained bounded below 2kBT during the entire simulation. It is worth emphazising that the dynamics generated by the Hamiltonian given in eq. (5 ) is, in reality, a fictitious one: it is simply a procedure for sampling configurations governed by the canonical distribution associated with V& . Consequently, the value of m* is arbitrary and can be chosen freely to explore efficiently the configuration space. We set m* = 40 amu so that the characteristic period of the intramolecular harmonic interactions corresponds roughly to 20 time steps. The choice of the reaction coordinate { deserves a brief comment. An especially convenient reaction coordinate is the centroid of the quantum path associated with the reactive quantum degree of freedom [ 46]. The centroid of the proton path, P,, is defined as d
ip = +
J
r(t)dt.
(6)
0
In our earlier study [I] we took the reaction coordinate to be the difference of solvent electrical potential energy, AE [ {I?} 1, between two points located along the interionic axis, a many-body function of the solvent coordinates. For adiabatic dynamics, we verified that the choice of AE as the reaction coordinate was dynamically equivalent to the position of the proton wave packet. In the present case we use the z component of the proton centroid, Z,, as the reaction coordinate. Considering the positions of the A- ions given earlier, configurations with negative and positive values of < = ZP correspond to reactant and product
states, respectively. By symmetry, the activated state lies at { = 0 = r+. The free energy profile, W(t), associated with the reaction coordinate was calculated using two independent methods. The first method involved the evaluation of the mean force F (l’) defined as F(c’)
= -9.
The free energy was obtained by simple integration of F ({’ ) calculated at selected values of <‘. The force F (5’) can be calculated conveniently using constrained molecular dynamics [ 171. By introducing the coordinate transformation {rl} ++ {ql, <}, where q1 is a set of generalized coordinates, F (t’) is given by [171
‘(“)
= (d(<(rt
- 5’))
-ly
P
(8) where lJ( is the Jacobian determinant (ar/a<(. The averages in eq. (8) can be calculated using the constrained-reaction-coordinate-dynamics (CRCD) ensemble [ 17 1:
_!%$!+ j+.!$
1) , f=f' . .
(9)
where (. . .)e=e, denotes a statistical average taken over the CRCD ensemble. Given the choice for the reaction coordinate, it is appropriate to define the coordinate transformation {rl } +-+{ t( } where
(10)
I=1
It is evident that r corresponds to the z component of r{ and the components of the remaining 4 are embodied in the set ql . The inverse of the transformation is straightforward to compute and it is easy to
D. Lara et al. / Chemical Physics 180 (1994) 181-189
185
show that the determinant ]J] is a constant. The final expression for the mean force takes the form
F(C) =
@),;,,=($$),.,.*
(11)
An alternative route for the calculation of free energy was proposed by Tobias and Brooks [ 181. Differences in free energy can be calculated’using the formula [ 18,191
-BOW
- W(b)1 = ln(ew(-PAV)t=t,, (12) -0.5 1 -1
where AI’ = KffLL{q(r)]l
- KdCl,{q(r))l.
Given the definition of the reaction coordinate, the value of AV is the difference in the potential energy between a configuration with the proton centroid constrained at r =
3. Results The free energy along the reaction coordinate is shown in fig. 1. This figure contains results obtained from the integration of the force, the direct estimate of the mean potential and from the histogram of a long ( 1.5 ns) unconstrained molecular dynamics run. The average force was calculated by constraining the z coordinate of the proton centroid at thirteen points spanning the range ZP = 0.0 A to ZP = 0.8 A and the free energy was obtained by integration of the average force acting on the proton-path using the trapezoidal rule. Very long runs, of the order of 2 ns, were required to equilibrate the force acting on the proton centroid. For the free energy calculated using the Tobias-Brooks scheme [ 181 we let G - {I = *St and performed a number of i constrained simulations, varying the value of C:by 26 in each simulation. In this way the potential of mean force could be computed over the range of interest within an additive constant. The agreement between the two simulation
I
I
I
-0.5
0
0.5
1
Z
Fig. 1. Free energy surfaces for the proton transfer reaction. The solid line corresponds to results obtained by integration of the mean force, the dashed line corresponds to the direct estimation of the mean potential and the open circles were obtained from an unconstrained MD simulation.
techniques is good and the activation energies differ by less than 5%. The average activation energy is estimated to be Wt = 1.36 kcal/ mol from these quantum simulations. The results obtained in the vicinity of activated state region exhibit larger fluctuations than for other values of the reaction coordinate. The open circles in fig. 1 were obtained by removing the constraint on the centroid and simply monitoring the value of the centroid in a long MD run. From these results the probability distribution of 5, could be estimated and the potential of mean force could then be obtained using eq. ( 1). Since the transfer process is activated, in such an unconstrained MD simulation the centroid will only take on values near the stable minima in the free energy and its shape in the barrier region cannot be determined. The figure shows that in the vicinities of the stable minima the results obtained from such unconstrained dynamics are in close accord with the two other constrained dynamics methods, providing additional confidence in the results. However, this method alone is not useful for the estimation of the barrier height or shape. One must rely on extrapolation from the regions of the minima and, as we have demonstrated, this procedure is subject to considerable uncertainty.
D. Lara et al. / Chemical Physics 180 (1994) 181-189
186
a
b
Fig 2. Proton polymer configurations with the centroid constrained in the product (a) and transition (b) states. The two ions between which the transfer occurs are depicted with arbitrary radius as are the polymer beads; solvent molecules are not shown but were included in the calculations used to produce the figure. Several features of the free energy profiles are worth noting. The curves are parabolic in only very narrow regions near the reactant and product states. Also, the free energy varies little within the range -0.2A < zp < 0.2 A. Since the free energy barrier is rather flat and the harmonic approximation in the stable well regions quickly breaks down, the free energy obtained by harmonic extrapolation from the reactant and product regions is likely to be an overestimate of the true value. Earlier [ 1 ] we estimated the free energy barrier to be Wt = 1.83 kcal/mol from adiabatic dynamics using AE as the reaction coordinate. We have computed AE in the centroid-constrained MD simulations and verified that the average value of AE tracks zp. Most of the free energy difference can be attributed to the fact that the barrier was determined by harmonic extrapolation from the behavior near the stable well regions, although there may be small corrections due to non-adiabaticity of the dynamics, especially in the vicinity of the transition state configuration. An impression of the quantum nature of the proton transfer process can be obtained from fig. 2 which shows two configurations of the polymer chain representing the proton with its centroid constrained in the
product (a) and transition (b) states. The stable positions of the centroid lie near zp x 4~0.65 A, which is the distance between the proton centroid in the stable regions and the transition state. From the figure one can see that the spatial extent of the proton is comparable to the distance over which it transfers: the proton cannot be treated as a point particle. It is possible to give a quantitative measure of the quantum dispersion of the proton as follows. The rms value of the displacement, R (t - t’ ), between two points on the proton path separated by a time increment 0 G t - t’ < /3h defined as R2(t - t’) = (Ir(t) - r(t’)l*),
(13)
provides an estimate of the spatial extent of the quantum particle. The characteristic size of the polymer is given by R = R(Ph/2). For a free particle,
(14) Rfree = flak
where ;1r =
(fi*Plm,) ‘I2 is the thermal wavelength associated with the proton. For a proton at T = 250 K, Rfree = 0.38 A. In fig. 3 R(t - t’) is plotted for a proton
D. Lara et al. / Chemical Physics 180 (1994) 181-189
.
.
I
.
.
‘._ ‘_ ‘.
“\
‘.
------mm____
L
0
3
0.3
’
..__
,U’ _’ _’
,
187
‘,,_
_-----
’
’
’
-..
Y cl
’
0.5
’
’
’
1
(t-t’)/ph Fig. 3. Root mean square displacement, R(t), as a function of imaginary time, t, for the proton. The solid line corresponds to the proton centroid constrained at the reactant or product states; the dashed line corresponds to the centroid constrained at the transition state; the dotted line corresponds to a free (non-interacting) proton. with its centroid contrained to lie at the stable and transition states. Both curves are similar in appearance with the proton somewhat more confined in the reactant-product states than in the transition state. In both cases the proton path has a size of roughly 0.3 A. The similar structure of both curves arises from the fact that there is almost no intrinsic barrier separating reactant and product configurations. For high barriers the proton will be more delocalized in the transition state region than in the stable-state regions. An examination of the behavior of R (t - t' ) as a funtion of t - t’ shows that the proton does not exhibit a clear dominance from the ground state; such ground state dominance is signalled by time independence of R (1 - t’ 1, except for t - t’ near 0 or /3h. For adiabatic dynamics to be valid an energy spacing large compared with ks T is required between the ground state and the manifold of excited states. In our model the energy spacing in the transition state region is comparable to /CBT and a partial breakdown of the adiabatic approximation is to be expected. The threedimensional character of the proton reaction coordinate plays an important role in the free energy of the
-0.3 4 3
I
0
0.3
X
Fig. 4. Scatter plot of the centroid position in the xy-plane with Z, = 0.
transfer process. Fig. 4 shows a scatter plot of the position of the centroid in the xy-plane with < = Z, = 0; i.e., with the centroid fixed at the transition state. Since only the z-coordinate of the centroid is fixed by the constraint the x and y values of its position are free to move. The scatter plot was obtained by recording the position of the centroid at equal time intervals from a long constrained MD simulation. It provides a measure of the centroid probability density in the xy-plane. While the dispersion does not appear to very large, x ~tO.2 A, it is sufficient for the centroid to explore regions of the potential comparable to and higher than the free energy barrier. The radial probability density is presented in fig. 5a and serves to illustrate that the average value of p, the radial distance from the line of centers joining the two ions, is non-zero. We have also carried out simulations of the transfer of a purely classical proton in the same solvent under the same conditions. The quantum and classical free energies are compared in fig. 6. It is clear from this figure that the quantum and classical free energy barriers differ by roughly 0.5kBT. The classical free energy profile is broader and the minima are displaced to larger separations. Like the quantum result, the classical free energy profile has broad relatively flat max-
188
D. Lara et al. / Chemical Physics 180 (1994) 181-189
0.00
0.05
0.10
0.15
0.20
0.25
0.30
035
P
Fig. 5. ( + ) Quantum radial probability density, p(p), for 5, = 0 (SOWp(p) dp = 1). (0) Classical radial probability density for the same constraint.
3.5 3 2.5 PW
2 1.5 1
imum and is harmonic only in the immediate vicinities of the minima. Other features of the classical calculation are similar to those of the quantum calculation; for example, the radial probability density for the classical proton with z-coordinate constrained at z = 0 is shown in fig. 5. Note the similarity with the centroid radial probability density: the classical proton is also able to explore the high energy regions of the potential normal to the line-of-centers. It is interesting to contrast some of the general features of our results with those of Lobaugh and Voth [7] on a different system. In their model system the intrinsic (bare) potential barrier is quite high. Most quantum effects arise from tunneling corrections which lower this bare barrier. The solvent contribution is not strongly affected by quantum corrections. The manner in which the quantum character of the proton affects the solvent contributions can be subtle since it depends on how the delocalized proton charge density interacts with the solvent degrees of freedom. In our strongly hydrogen bonded system there is virtually no barrier in the bare potential. The quantum effects come into play in the manner in which the proton charge density interacts with the surrounding classical solvent particles. We noted that the proton cannot be considered as a classical point particle since its charge density is delocalized over a distance that is comparable to the proton transfer distance. This delocalization is responsible for the lower and narrower barrier obtained in the quantum simulation. Consequently, our results show that the quantum character of the proton can indeed affect the solvent contribution to the free energy. We note that this correction lowers and broadens the free energy barrier, like tunneling corrections to the bare potential, but arises from a different mechanism.
0.5
4. Summary
0 -0.5
-1
-0.5
0
0.5
1
z
The results presented here demostrate the feasibility of using constrained molecular dynamics methods for the computation of free energy of quantum acti-
Fig. 6. Comparison
vated
tained by taking the average of the solid and dashed lines of fig. 1.
two methods we employed to obtain the free energy are in close accord and also compare quite well with our previous estimate for the free energy barrier. The path integral method is not limited to the adiabatic
of quantum (solid line) and classical (dashed line) free energies as a function of the reaction coordinate. The quantum result shown in this figure is ob-
transition
processes
in condensed
media.
The
D. Lara et al. / Chemical Physics 180 (1994) 181-189
regime, as were our earlier real-time dynamical calculations. It is also worth pointing out the important role that the three-dimensional character of the quantum reaction coordinate plays in the calculation. Motion normal to the line of centers joining the ions allows the system to access regions of configuration space of higher energy and, in fact, these trajectories carry most of the statistical weight. The system we have chosen to study presents a rather subtle interplay of solvent and quantum effects. Since the bare potential has a negligible barrier, it is solely the interaction of the proton charge density with the solvent molecules that is responsible for the activation free energy: small shifts in the solvent polarization are sufficient to induce localization of the proton, which in turn locks the solvent in a configuration that favors the localized state. The quantum effects arise from the dispersion of the proton charge density and the consequent interaction with the solvent that this dispersion produces. Our results show that because the proton density is dispersed on a length scale that is comparable to that over which the proton transfers, there are deviations from the purely classical calculation where the proton is a point particle. Since in our model there is a negligible intrinsic barrier, thus negligible tunneling corrections, all of these free energy changes can be attributed to solvent interactions. Such quantum effects on solvent contribution to the free energy will be more pronounced if the particle being transferred is more delocalized. This will be the case if a muon is substituted for the proton. Preliminary results indicate that there is a dramatic reduction and narrowing of the free energy barrier. Our calculations illustrate the delicate interplay between charge density dispersion and proton-solvent interactions that combine to determine the free energy barrier.
Acknowledgement This research was partially supported by the Italian CNR via the Cray Project on Statistical Mechanics and the P.F. “Sistemi Informatici e Calcolo Paral-
189
lelo”. D.L. is a recipient of a Postdoctoral Fellowship awarded by the Commission of the European Communities within the framework of the Scientific Cooperation Agreement CEE-Argentina. The research of R.K. was supported in part by a grant from the Natural Sciences and Engineering Research Council of Canada and a grant funded by the Network of Centres of Excellence program in association with the Natural Sciences and Engineering Research Council of Canada.
References [ 1 ] D. Laria, G. Ciccotti, M. Ferrario and R. Kapral, J. Chem. Phys. 97 (1992) 378. [2] R.A. Marcus and N. Sutin, Biochim. Biophys. Acta 811 (1985) 265. [ 31 R.P. Feynman, Statistical mechanics (Addison Wesley, Reading, 1972). [4] M. Gillan, J. Phys. C 20 (1987) 3621. [ 51 G.A. Voth, D. Chandler and W.H. Miller, J. Chem. Phys. 91 (1989) 7749. [6] D. Li and G. Voth, J. Phys. Chem. 95 (1991) 10425. [7] J. Lobaugh and G.A. Voth, Chem. Phys. Letters 198 (1992) 311. [8] D. Borgis, G. Tajus and H. Azzouz, J. Phys. Chem. 96 (1992) 3188. [9] T. Yamamoto, J. Chem. Phys. 33 (1960) 281; D. Chandler, J. Chem. Phys. 68 (1978) 2959. [lo] D. Chandler and P.G. Wolynes, J. Chem. Phys. 74 (1981) 4078. [ 111 M. Parrinello and A. Rahman, J. Chem. Phys. 80 (1984) 860. [ 121 B. de Raedt, M. Sprik and M.L. Klein, J. Chem. Phys. 80 (1984) 5719. [ 131 J.P. Valleau and G.M. Tonic, in: Statistical mechanics, part A, ed. B.J. Beme (Plenum Press, New York, 1977).
[ 141 L. Verlet, Phys.Rev. 159 (1967) 98. [ 151 J.P. Ryckaert, G. Ciccotti and H.J.C. Berendsen,
J. Comput. Phys. 23 (1977) 327. [16] S. Nose, Mol. Phys. 52 (1984) 255. [ 171 E.A. Carter, G. Ciccotti, J.T. Hynes and R. Kapral, Chem. Phys. Letters 156 (1989) 472. [ 181 D.J. Tobias and C.L. Brooks, Chem. Phys. Letters 142 (1987) 472. [ 191 E. Paci, G. Ciccotti, M. Ferrario and R. Kapral, Chem. Phys. Letters 176 (1991) 581.
Activation free energy for proton transfer in solution Daniel Laria l, Giovanni
Ciccotti 2
Centre Europken de Calcul Atomique et Mokdaire,
Bbtiment 506, Universitk de Paris-Sud, 91405 Orsay, France
Mauro Ferrario Dipartimento di Fisica, Universitb di Messina, CP 50, 98166 S. Agata di Messina, Italy
and Raymond
Kapral
Chemical Physics Theory Group, Department of Chemistry Universityof Toronto, Toronto, Ontario, Canada M5S IA1
Received 18 May 1993
Path integral molecular dynamics methods are employed to compute the free energy for proton transfer reactions for strongly hydrogen bonded systems in a polar solvent. The free energy profile is calculated using several different techniques, including: integration of the mean force acting on the proton path with its centroid constrained at different values, the integral form of the free energy calculation in the constrained-reaction-coordinate-dynamics ensemble and direct simulation of the unconstrained dynamics. The results show that estimates of the free energy barrier obtained by harmonic extrapolation are likely to be in error. Both quantum and classical results for the free energy are obtained and compared with simulations using adiabatic quantum dynamics. Comparison of the quantum and classical results show that there are quantum corrections to the solvent contributions to the free energy.
1. Introduction
Proton transfer is a process that occurs in many chemical and biological reactions. The calculation of such transfer processes is difficult because quantum effects may be important and the reaction coordinate that is used to specify the progress of the reaction is often taken to involve the solvent degrees of freedom. In general, proton transfer involves the description of the dynamics of a quantum particle in a many-body solvent. In a recent paper [ 1] we presented a molecular dynamics calculation of the adiabatic transfer of a prol
On leave from Departamento Quimica de Reactores, Comision National de Energia Atomica, ( 1429) Buenos Aires, Argentina. * Permanent address: Dipartimento di Fisica, Universita “La Sapienza”, Piazzale Aldo Moro, 2, 00185 Rome, Italy. 0301-0104/94/$7.00
ton between two negatively charged ions immersed in a classical dipolar solvent. In this adiabatic limit the proton charge density can instantaneously adapt to the configuration of the classical particles; thus, the implementation of the dynamics is conceptually simple (but computationally demanding) and involves the solution of the Schrodinger equation for the proton, which is a function of the instantaneous solvent configuration, and Newton’s equations of motion for the solvent, which depend on the proton charge density. The proton transfer is governed by local solvent fluctuations in the vicinity of the proton-ion complex and the progress of the reaction can be monitored by choosing the solvent polarization, AE ( {R} 1, as a reaction coordinate [ 21. Consequently, the quantum problem reduces to a classical one with a complicated, many-body reaction coordinate. As a result of the preferential solvation of the complex with the proton charge localized near one of the two ions, the fluc-
@ 1994-Elsevier Science B.V. All rights reserved.
SSDI 0301-0104(93)00002-R
182
D. Lara et al. / Chemical Physics I80 (1994) 181-189
tuations of the solvent that induce the proton transfer are rare events. The asymmetry in the charge distribution of the complex leads to a larger dipole moment of the complex and to stronger solvation energy. The passage of the proton from reactant to product states requires considerable energy and involves substantial modification in the local structure of the solvent. The magnitude of the free energy barrier relative to the thermal energy determines the order of magnitude of the proton transfer rate. In our previous study, the free energy barrier was estimated from an extrapolation to the transition state of the logarithm of the probability density of the reaction coordinate computed near the reactant and product states. The free energy profile in the stable state regions was found to be approximately harmonic, in accord with Marcus’s theory [2] for charge transfer processes. In this article, we use Feynman’s [ 31 path integral formulation of quantum mechanics to treat the proton degrees of freedom, and choose the centroid as the reaction coordinate [4-61. Our calculations are carried out using constrained molecular dynamics methods. This study has several goals. We consider various constrained molecular dynamics methods for the quantum path integral calculation of the free energy. The calculations reported here using the centroid reaction coordinate are compared with our earlier results using AE as the reaction coordinate and adiabatic dynamics. We discuss features associated with the three-dimensional character of our reaction coordinate. Finally, we compare the quantum calculations with classical simulations on the same model. The present calculations are carried out in the same spirit as other recent calculations of the free energy of proton transfer; for example, we mention work by Lobaugh and Voth [7] for a full three-dimensional treatment of the centroid reaction coordinate for proton transfer in water, the analog of our calculations for a strongly hydrogen bonded system, and by Borgis et al. [8] for proton transfer in a dipolar solvent like the one considered here.
2. Activation free energy for proton transfer The methodology for the evaluation of the rate constants of simple classical activated processes is well
established and the computer simulation of such processes does not present a major challenge. Normally, the calculation proceeds in two steps: first, the free energy along the reaction coordinate is calculated. This free energy calculation is used to locate the bottleneck or transition state configuration from which the transition state theory estimate of the rate constant may be obtained [ 91. Second, corrections to transition state theory are obtained by evaluating the reactive flux correlation function k(t) [91. The transmission coefficient can be extracted from the plateau value of k (f ) provided the characteristic time scale for the reaction is much longer than that of any other microscopic relaxation times in the system. Only in this circumstance will a phenomenological description of the reactive process be meaningful. In this article we focus on the first of these steps, the calculation of the free energy for a quantum reaction coordinate. The calculations are considerably more complicated when quantum effects cannot be ignored. It is often convenient to cast the description of such quantum rate processes in terms of Feynman’s path integral formulation of quantum mechanics [ 31. Let 5 [r (t ) ] be the reaction coordinate, which is a functional of the coordinates, r (t 1, of the proton path at the imaginary time t, and 5 be its numerical value. The free energy II’ (5) is related to the probability density of < [r (t ) 1,
where (. . .) represents a canonical statistical average over all degrees of freedom of the system and /3-’ is Boltzmann’s constant times the temperature. Physically, I+‘(c) represents the reversible work required to move the reaction coordinate to the point <. If tt is the value of r at the transition state, then II’ s Wt can be used to compute the transition state theory approximation to the quantum rate constant. A discussion of the full rate constant and its breakup into transition state theory and correction terms can be found in refs. [ 5 ] and [ 6 1. In the case where quantum fluctuations can be neglected in all particles with the exception of the proton, the canonical average in eq. ( 1) can be expressed as
D. Lara et al. / Chemical Physics 180 (1994) 181-189
2.2. X
J.. J 9
wbl
s(ti[r(t)l
x exp(-S[rUL{R)l/h)
-0
,
(2)
where Vci({I?} ) represents the potential energy of all the classical particles with coordinates {R}, Q is the partition function of the system and D [ r (t ) ] represents the path integral. For a given configuration of the classical particles, the weight functional for the proton path is governed by the action
Simulation
+ J dt{~~pl~(t)12 &[r(f);{R}l},
method
Simulation techniques can be constructed to eMiciently sample fluctuations of the proton path. The proton path maps onto a ring polymer containing P interacting sites .with harmonic interactions between nearest neighbors [ lo]. This isomorphism can be exploited to construct a simulation algorithm. For a finite value of P, the problem reduces to the consideration of a classical system with an effective potential given by Kn = V,I({RI)
=
183
Pm, ’ + 2fph_2 c,=, b.1 -rl+l)* (41
(31
0
where mp is the mass of the proton, VP[r(t); {R} ] is the potential energy of the proton path in the field of the classical particles and h is Plan&s constant divided by 2~.
2. I. Model system The system studied is identical to that in our previous investigation [ 11. It consists of a proton and a pair of A- ions solvated by N = 342 dipolar molecules. The temperature of the system is 250 K and the length of the simulation box is 30.52 A. It is convenient to define a coordinate system with origin at the center of the simulation box. The two A- ions are fixed at Ri,z = (0,O) f 1.3 A. The solvent molecules are composed of two interaction sites with partial charges of opposite sign zu = f0.52e and a fixed bond length of 2 A. Site-site interactions arise from 6-12 LennardJones and Coulomb forces. The proton-ion potential was constructed to model strongly hydrogen bonded systems: the proton is confined by the potential to lie close to the ion-ion internuclear axis and there is a very small barrier in the potential (w 0.2 kcal/mol) along this axis (cf. fig. 1 of ref. [ 1] ) . Finally, the proton interacts with the solvent via simple Coulombic forces. The details of the potential energy parameters are given in ref. [ 11.
where r, are the coordinates of the polymer beads. The accuracy improves as P increases and the isomorphism becomes exact in the limit P --+ CO. The convergence in the quantities of interest determines the choice of P. We found no noticeable systematic changes for P > 20; this number has a magnitude similar to that used in previous path integral simulations of proton-like particles [ 11,12 1. Several methods exist for the calculation of free energy profiles for quantum activated processes. Lobaugh et al. [ 71 have used umbrella-sampling Monte Carlo [ 131 to treat proton transfer reactions in aqueous media. Here, we employ constrained molecular dynamics to compute the free energy of the proton transfer process. In particular, constanttemperature molecular dynamics trajectories were generated by the following Hamiltonian: P
H = c 1=1
N
fm*;f
+ c fM$: JCI
+ V, + Vk,
(5)
where m* is the fictitious mass of a polymer bead and Mj are the masses of the classical particles. The equations of motion were integrated using the Verlet algorithm [ 141 with a time step of 10pL4 s. Geomet-! rical constraints from intramolecular bonds and constraints acting on the reaction coordinate were treated with the SHAKE algorithm [ 151. Sampling from a canonical distribution was carried out using NosCHoover dynamics [ 161. We found it quite difficult to obtain proper thermalization of the proton-polymer
184
D. Lara et al. / Chemical Physics 180 (1994) 181-189
degrees of freedom even after long thermalization periods. Similar difftculties were pointed out in the muonium simulation described in ref. [ lo]. Although the total kinetic energy of the system was well controlled by the Nose-Hoover thermostat, we noticed a systematic drift towards higher temperatures in the kinetic energy associated with the degrees of freedom of the polymer. For this reason, a second thermostat acting exclusively on the proton-polymer degrees of freedom had to be included. Only in this case could we obtain reasonable thermal equilibrium between the proton and the solvent. The inclusion of a second thermostat introduces a net acceleration of the center of mass of the system; however, we verified that the center of mass oscillates and that the kinetic energy associated with it remained bounded below 2kBT during the entire simulation. It is worth emphazising that the dynamics generated by the Hamiltonian given in eq. (5 ) is, in reality, a fictitious one: it is simply a procedure for sampling configurations governed by the canonical distribution associated with V& . Consequently, the value of m* is arbitrary and can be chosen freely to explore efficiently the configuration space. We set m* = 40 amu so that the characteristic period of the intramolecular harmonic interactions corresponds roughly to 20 time steps. The choice of the reaction coordinate { deserves a brief comment. An especially convenient reaction coordinate is the centroid of the quantum path associated with the reactive quantum degree of freedom [ 46]. The centroid of the proton path, P,, is defined as d
ip = +
J
r(t)dt.
(6)
0
In our earlier study [I] we took the reaction coordinate to be the difference of solvent electrical potential energy, AE [ {I?} 1, between two points located along the interionic axis, a many-body function of the solvent coordinates. For adiabatic dynamics, we verified that the choice of AE as the reaction coordinate was dynamically equivalent to the position of the proton wave packet. In the present case we use the z component of the proton centroid, Z,, as the reaction coordinate. Considering the positions of the A- ions given earlier, configurations with negative and positive values of < = ZP correspond to reactant and product
states, respectively. By symmetry, the activated state lies at { = 0 = r+. The free energy profile, W(t), associated with the reaction coordinate was calculated using two independent methods. The first method involved the evaluation of the mean force F (l’) defined as F(c’)
= -9.
The free energy was obtained by simple integration of F ({’ ) calculated at selected values of <‘. The force F (5’) can be calculated conveniently using constrained molecular dynamics [ 171. By introducing the coordinate transformation {rl} ++ {ql, <}, where q1 is a set of generalized coordinates, F (t’) is given by [171
‘(“)
= (d(<(rt
- 5’))
-ly
P
(8) where lJ( is the Jacobian determinant (ar/a<(. The averages in eq. (8) can be calculated using the constrained-reaction-coordinate-dynamics (CRCD) ensemble [ 17 1:
_!%$!+ j+.!$
1) , f=f' . .
(9)
where (. . .)e=e, denotes a statistical average taken over the CRCD ensemble. Given the choice for the reaction coordinate, it is appropriate to define the coordinate transformation {rl } +-+{ t( } where
(10)
I=1
It is evident that r corresponds to the z component of r{ and the components of the remaining 4 are embodied in the set ql . The inverse of the transformation is straightforward to compute and it is easy to
D. Lara et al. / Chemical Physics 180 (1994) 181-189
185
show that the determinant ]J] is a constant. The final expression for the mean force takes the form
F(C) =
@),;,,=($$),.,.*
(11)
An alternative route for the calculation of free energy was proposed by Tobias and Brooks [ 181. Differences in free energy can be calculated’using the formula [ 18,191
-BOW
- W(b)1 = ln(ew(-PAV)t=t,, (12) -0.5 1 -1
where AI’ = KffLL{q(r)]l
- KdCl,{q(r))l.
Given the definition of the reaction coordinate, the value of AV is the difference in the potential energy between a configuration with the proton centroid constrained at r =
3. Results The free energy along the reaction coordinate is shown in fig. 1. This figure contains results obtained from the integration of the force, the direct estimate of the mean potential and from the histogram of a long ( 1.5 ns) unconstrained molecular dynamics run. The average force was calculated by constraining the z coordinate of the proton centroid at thirteen points spanning the range ZP = 0.0 A to ZP = 0.8 A and the free energy was obtained by integration of the average force acting on the proton-path using the trapezoidal rule. Very long runs, of the order of 2 ns, were required to equilibrate the force acting on the proton centroid. For the free energy calculated using the Tobias-Brooks scheme [ 181 we let G - {I = *St and performed a number of i constrained simulations, varying the value of C:by 26 in each simulation. In this way the potential of mean force could be computed over the range of interest within an additive constant. The agreement between the two simulation
I
I
I
-0.5
0
0.5
1
Z
Fig. 1. Free energy surfaces for the proton transfer reaction. The solid line corresponds to results obtained by integration of the mean force, the dashed line corresponds to the direct estimation of the mean potential and the open circles were obtained from an unconstrained MD simulation.
techniques is good and the activation energies differ by less than 5%. The average activation energy is estimated to be Wt = 1.36 kcal/ mol from these quantum simulations. The results obtained in the vicinity of activated state region exhibit larger fluctuations than for other values of the reaction coordinate. The open circles in fig. 1 were obtained by removing the constraint on the centroid and simply monitoring the value of the centroid in a long MD run. From these results the probability distribution of 5, could be estimated and the potential of mean force could then be obtained using eq. ( 1). Since the transfer process is activated, in such an unconstrained MD simulation the centroid will only take on values near the stable minima in the free energy and its shape in the barrier region cannot be determined. The figure shows that in the vicinities of the stable minima the results obtained from such unconstrained dynamics are in close accord with the two other constrained dynamics methods, providing additional confidence in the results. However, this method alone is not useful for the estimation of the barrier height or shape. One must rely on extrapolation from the regions of the minima and, as we have demonstrated, this procedure is subject to considerable uncertainty.
D. Lara et al. / Chemical Physics 180 (1994) 181-189
186
a
b
Fig 2. Proton polymer configurations with the centroid constrained in the product (a) and transition (b) states. The two ions between which the transfer occurs are depicted with arbitrary radius as are the polymer beads; solvent molecules are not shown but were included in the calculations used to produce the figure. Several features of the free energy profiles are worth noting. The curves are parabolic in only very narrow regions near the reactant and product states. Also, the free energy varies little within the range -0.2A < zp < 0.2 A. Since the free energy barrier is rather flat and the harmonic approximation in the stable well regions quickly breaks down, the free energy obtained by harmonic extrapolation from the reactant and product regions is likely to be an overestimate of the true value. Earlier [ 1 ] we estimated the free energy barrier to be Wt = 1.83 kcal/mol from adiabatic dynamics using AE as the reaction coordinate. We have computed AE in the centroid-constrained MD simulations and verified that the average value of AE tracks zp. Most of the free energy difference can be attributed to the fact that the barrier was determined by harmonic extrapolation from the behavior near the stable well regions, although there may be small corrections due to non-adiabaticity of the dynamics, especially in the vicinity of the transition state configuration. An impression of the quantum nature of the proton transfer process can be obtained from fig. 2 which shows two configurations of the polymer chain representing the proton with its centroid constrained in the
product (a) and transition (b) states. The stable positions of the centroid lie near zp x 4~0.65 A, which is the distance between the proton centroid in the stable regions and the transition state. From the figure one can see that the spatial extent of the proton is comparable to the distance over which it transfers: the proton cannot be treated as a point particle. It is possible to give a quantitative measure of the quantum dispersion of the proton as follows. The rms value of the displacement, R (t - t’ ), between two points on the proton path separated by a time increment 0 G t - t’ < /3h defined as R2(t - t’) = (Ir(t) - r(t’)l*),
(13)
provides an estimate of the spatial extent of the quantum particle. The characteristic size of the polymer is given by R = R(Ph/2). For a free particle,
(14) Rfree = flak
where ;1r =
(fi*Plm,) ‘I2 is the thermal wavelength associated with the proton. For a proton at T = 250 K, Rfree = 0.38 A. In fig. 3 R(t - t’) is plotted for a proton
D. Lara et al. / Chemical Physics 180 (1994) 181-189
.
.
I
.
.
‘._ ‘_ ‘.
“\
‘.
------mm____
L
0
3
0.3
’
..__
,U’ _’ _’
,
187
‘,,_
_-----
’
’
’
-..
Y cl
’
0.5
’
’
’
1
(t-t’)/ph Fig. 3. Root mean square displacement, R(t), as a function of imaginary time, t, for the proton. The solid line corresponds to the proton centroid constrained at the reactant or product states; the dashed line corresponds to the centroid constrained at the transition state; the dotted line corresponds to a free (non-interacting) proton. with its centroid contrained to lie at the stable and transition states. Both curves are similar in appearance with the proton somewhat more confined in the reactant-product states than in the transition state. In both cases the proton path has a size of roughly 0.3 A. The similar structure of both curves arises from the fact that there is almost no intrinsic barrier separating reactant and product configurations. For high barriers the proton will be more delocalized in the transition state region than in the stable-state regions. An examination of the behavior of R (t - t' ) as a funtion of t - t’ shows that the proton does not exhibit a clear dominance from the ground state; such ground state dominance is signalled by time independence of R (1 - t’ 1, except for t - t’ near 0 or /3h. For adiabatic dynamics to be valid an energy spacing large compared with ks T is required between the ground state and the manifold of excited states. In our model the energy spacing in the transition state region is comparable to /CBT and a partial breakdown of the adiabatic approximation is to be expected. The threedimensional character of the proton reaction coordinate plays an important role in the free energy of the
-0.3 4 3
I
0
0.3
X
Fig. 4. Scatter plot of the centroid position in the xy-plane with Z, = 0.
transfer process. Fig. 4 shows a scatter plot of the position of the centroid in the xy-plane with < = Z, = 0; i.e., with the centroid fixed at the transition state. Since only the z-coordinate of the centroid is fixed by the constraint the x and y values of its position are free to move. The scatter plot was obtained by recording the position of the centroid at equal time intervals from a long constrained MD simulation. It provides a measure of the centroid probability density in the xy-plane. While the dispersion does not appear to very large, x ~tO.2 A, it is sufficient for the centroid to explore regions of the potential comparable to and higher than the free energy barrier. The radial probability density is presented in fig. 5a and serves to illustrate that the average value of p, the radial distance from the line of centers joining the two ions, is non-zero. We have also carried out simulations of the transfer of a purely classical proton in the same solvent under the same conditions. The quantum and classical free energies are compared in fig. 6. It is clear from this figure that the quantum and classical free energy barriers differ by roughly 0.5kBT. The classical free energy profile is broader and the minima are displaced to larger separations. Like the quantum result, the classical free energy profile has broad relatively flat max-
188
D. Lara et al. / Chemical Physics 180 (1994) 181-189
0.00
0.05
0.10
0.15
0.20
0.25
0.30
035
P
Fig. 5. ( + ) Quantum radial probability density, p(p), for 5, = 0 (SOWp(p) dp = 1). (0) Classical radial probability density for the same constraint.
3.5 3 2.5 PW
2 1.5 1
imum and is harmonic only in the immediate vicinities of the minima. Other features of the classical calculation are similar to those of the quantum calculation; for example, the radial probability density for the classical proton with z-coordinate constrained at z = 0 is shown in fig. 5. Note the similarity with the centroid radial probability density: the classical proton is also able to explore the high energy regions of the potential normal to the line-of-centers. It is interesting to contrast some of the general features of our results with those of Lobaugh and Voth [7] on a different system. In their model system the intrinsic (bare) potential barrier is quite high. Most quantum effects arise from tunneling corrections which lower this bare barrier. The solvent contribution is not strongly affected by quantum corrections. The manner in which the quantum character of the proton affects the solvent contributions can be subtle since it depends on how the delocalized proton charge density interacts with the solvent degrees of freedom. In our strongly hydrogen bonded system there is virtually no barrier in the bare potential. The quantum effects come into play in the manner in which the proton charge density interacts with the surrounding classical solvent particles. We noted that the proton cannot be considered as a classical point particle since its charge density is delocalized over a distance that is comparable to the proton transfer distance. This delocalization is responsible for the lower and narrower barrier obtained in the quantum simulation. Consequently, our results show that the quantum character of the proton can indeed affect the solvent contribution to the free energy. We note that this correction lowers and broadens the free energy barrier, like tunneling corrections to the bare potential, but arises from a different mechanism.
0.5
4. Summary
0 -0.5
-1
-0.5
0
0.5
1
z
The results presented here demostrate the feasibility of using constrained molecular dynamics methods for the computation of free energy of quantum acti-
Fig. 6. Comparison
vated
tained by taking the average of the solid and dashed lines of fig. 1.
two methods we employed to obtain the free energy are in close accord and also compare quite well with our previous estimate for the free energy barrier. The path integral method is not limited to the adiabatic
of quantum (solid line) and classical (dashed line) free energies as a function of the reaction coordinate. The quantum result shown in this figure is ob-
transition
processes
in condensed
media.
The
D. Lara et al. / Chemical Physics 180 (1994) 181-189
regime, as were our earlier real-time dynamical calculations. It is also worth pointing out the important role that the three-dimensional character of the quantum reaction coordinate plays in the calculation. Motion normal to the line of centers joining the ions allows the system to access regions of configuration space of higher energy and, in fact, these trajectories carry most of the statistical weight. The system we have chosen to study presents a rather subtle interplay of solvent and quantum effects. Since the bare potential has a negligible barrier, it is solely the interaction of the proton charge density with the solvent molecules that is responsible for the activation free energy: small shifts in the solvent polarization are sufficient to induce localization of the proton, which in turn locks the solvent in a configuration that favors the localized state. The quantum effects arise from the dispersion of the proton charge density and the consequent interaction with the solvent that this dispersion produces. Our results show that because the proton density is dispersed on a length scale that is comparable to that over which the proton transfers, there are deviations from the purely classical calculation where the proton is a point particle. Since in our model there is a negligible intrinsic barrier, thus negligible tunneling corrections, all of these free energy changes can be attributed to solvent interactions. Such quantum effects on solvent contribution to the free energy will be more pronounced if the particle being transferred is more delocalized. This will be the case if a muon is substituted for the proton. Preliminary results indicate that there is a dramatic reduction and narrowing of the free energy barrier. Our calculations illustrate the delicate interplay between charge density dispersion and proton-solvent interactions that combine to determine the free energy barrier.
Acknowledgement This research was partially supported by the Italian CNR via the Cray Project on Statistical Mechanics and the P.F. “Sistemi Informatici e Calcolo Paral-
189
lelo”. D.L. is a recipient of a Postdoctoral Fellowship awarded by the Commission of the European Communities within the framework of the Scientific Cooperation Agreement CEE-Argentina. The research of R.K. was supported in part by a grant from the Natural Sciences and Engineering Research Council of Canada and a grant funded by the Network of Centres of Excellence program in association with the Natural Sciences and Engineering Research Council of Canada.
References [ 1 ] D. Laria, G. Ciccotti, M. Ferrario and R. Kapral, J. Chem. Phys. 97 (1992) 378. [2] R.A. Marcus and N. Sutin, Biochim. Biophys. Acta 811 (1985) 265. [ 31 R.P. Feynman, Statistical mechanics (Addison Wesley, Reading, 1972). [4] M. Gillan, J. Phys. C 20 (1987) 3621. [ 51 G.A. Voth, D. Chandler and W.H. Miller, J. Chem. Phys. 91 (1989) 7749. [6] D. Li and G. Voth, J. Phys. Chem. 95 (1991) 10425. [7] J. Lobaugh and G.A. Voth, Chem. Phys. Letters 198 (1992) 311. [8] D. Borgis, G. Tajus and H. Azzouz, J. Phys. Chem. 96 (1992) 3188. [9] T. Yamamoto, J. Chem. Phys. 33 (1960) 281; D. Chandler, J. Chem. Phys. 68 (1978) 2959. [lo] D. Chandler and P.G. Wolynes, J. Chem. Phys. 74 (1981) 4078. [ 111 M. Parrinello and A. Rahman, J. Chem. Phys. 80 (1984) 860. [ 121 B. de Raedt, M. Sprik and M.L. Klein, J. Chem. Phys. 80 (1984) 5719. [ 131 J.P. Valleau and G.M. Tonic, in: Statistical mechanics, part A, ed. B.J. Beme (Plenum Press, New York, 1977).
[ 141 L. Verlet, Phys.Rev. 159 (1967) 98. [ 151 J.P. Ryckaert, G. Ciccotti and H.J.C. Berendsen,
J. Comput. Phys. 23 (1977) 327. [16] S. Nose, Mol. Phys. 52 (1984) 255. [ 171 E.A. Carter, G. Ciccotti, J.T. Hynes and R. Kapral, Chem. Phys. Letters 156 (1989) 472. [ 181 D.J. Tobias and C.L. Brooks, Chem. Phys. Letters 142 (1987) 472. [ 191 E. Paci, G. Ciccotti, M. Ferrario and R. Kapral, Chem. Phys. Letters 176 (1991) 581.