J. Phys. Chem. 1993,97, 13469-13476
13469
Influence of Solvent on Intramolecular Proton Transfer in Hydrogen Malonate. Molecular Dynamics Simulation Study of Tunneling by Density Matrix Evolution and Nonequilibrium Solvation Janez Mavri,? Herman J. C. Berendsen,’ and Wilfred F. van Gunsteren* BIOSON Research Institute, Department of Biophysical Chemistry, the University of Groningen, Nijenborgh 4, 9747 AG Groningen, the Netherlands Received: November 3, 1992; In Final Form: May 7, I993
A density matrix evolution (DME) method (Berendsen, H. J. C.; Mavri, J. J. Phys. Chem. the preceding paper in this issue) in combination with classical molecular dynamics simulation was applied to calculate the rate of proton tunneling in the intramolecular double-well hydrogen bond of hydrogen malonate (HM) in aqueous solution. The one-dimensional time-dependent Schrdinger equation was solved for an ensemble of 50 000 configurations to obtain a time course of the splitting and coupling matrix elements between the product and reactant states. The rate of tunneling is in agreement with the approximate analytical solutions given by Borgis and Hynes (J. Chem. Phys. 1991, 94, 3619). When the in vacuo contribution to the free energy difference between the symmetric and asymmetric form of H M is omitted, the free energy difference has the meaning of the reversible work necessary for such a solvent reorganization that zero energy splitting between product and reactant states is achieved. The method of biased sampling is used to enhance the sampling of such configurations where tunneling probability is large. After the correction for biased sampling, the DME method results in practically the same rate of tunneling as obtained by an unbiased calculation. The rates of tunneling are compared to the classical proton-transfer rate in the nonadiabatic limit.
1. Introduction Proton transfer in hydrogen-bonded systems is an essential step in many important processes, including both biochemical processes as enzymatic reactions and transmission of signals in biologicalsystems and electrochemical processes. It is important that such processes are understood and can be studied on a first principles basis. Calculation of a rate constant for any reaction in the condensed phase is a nontrivial task. Accurate ab-initio calculations are usually limited to a few atoms directly involved in the chemical reaction. Ab-initio calculations provide information about transition states and activation energies for a static, truncated system. In addition, they form a valuable source for potentials for intramolecular and intermolecular degrees of freedom. Since the complexity of the system of interest usually does not allow analytical solutions, one must proceed with sampling of the conformationalspace by computer simulationsin order to calculate free energy differences associated with the change of reaction coordinate. On can then apply transition-state theory (TST) and calculate the classical rate constant directly from free energy differences. A nice example is the study of a simple S Nreaction ~ in aqueous solution.’ The procedure can be extended to correct the calculated rate constant for the trajectory recrossings at the transition state that are not incorporated in TST2and to take into account nonequilibrium sol~ation.3,~ Proton transfer in hydrogen-bonded systems can be. treated in a similar way, but there are at least two points that need to be considered first: the predicted activation energies and the nature of the proton motion. In ab-initio calculations, predicted proton potentials are extremely sensitive to the quality of the electronic wave function and are often not conclusive about the barrier height. Even for hydrogen-bonded systems of moderate size, the Hartree-Fock limit is not achieved with modern computers. The
* To whom correspondence should be addressed. t Permanent address: Institute of Chemistry, POB 30,61115, Ljubljana, Slovenia. Present address: ETH,Physical Chemistry, CH-8092,Zllrich, Switzerland.
*
0022-365419312097-13469$04.00/0
other point is the quantum nature of the proton. The assumption that atoms obey the laws of classical mechanics is valid until the thermal wavelength X A=
h
(1) (keTm)”2 is of the same order as the typical length characterizing the interatomic potential.2 For a proton at room temperature, the thermal wavelength is about 0.3 A, indicating that protons moving in a hydrogen bond will exhibit quantum behavior. In the context of a study on the proton-transfer rates, tunneling is the most interesting quantum phenomenon. In this study, we address the intramolecular proton-transfer rate in hydrogen malonate (HM) in aqueous solution. For HM, high-quality ab-initio calculations were done and some experimental data are available. Stability of the six-membered ring, formed by intramolecular hydrogen bonding in HM, was proposed to be the reason for the anomalous ratio of the first and second dissociationconstantsof malonic acid.5 Although HM does indeed appear in this form in some salts in the solid state (ref 6 and references therein), the Raman spectroscopicdata of Maury and Bardet did not confirm the existence of intramolecular hydrogen bonding in aqueous s ~ l u t i o n . However, ~ differences in 13C chemical shifts of the methylene carbon of malonic acid with progressive neutralization were interpreted in terms of intramolecular H-bondingin HM.8 In solution,oneexpectsan equilibrium between the planar E-conformer, which is characterized by an intramolecular hydrogen bond, and the open Z-conformer. We restrict our considerationsto the E-conformer and do not consider the general problem of proton transfer including reaction paths that involve solvent molecules. Experimental data give evidence of a symmetric double-well proton potential in solid guanidinium HM.6 All calculations on HM at a high level of ab-initio theory9 predict the proton potential to be of symmetric double-well character, although they are not conclusive about the precise height of the barrier. Inclusion of the solvent reaction field intoab-initio calculations by the method of Miertui and Tomasilo predicts an increase of the potential 0 1993 American Chemical Society
Mavri et al.
13470 The Journal of Physical Chemistry, Vol. 97, No. 51, 1993 , .
TABLE I: Parameters for the Nonbonding Interactions Used in the MD SimUl.ti0~8of Hydrogen Malonate and Water Molecules with Atomic Cbargea Given in Atomic Units atom 4" qb 4' (C12)1/2 (C6)1/2
05
01
x
C A J J
Figure 1. Structure and atom numbering of hydrogen malonate (HM).
Methylene group C3 was treated in the simulationsas a united atom. The energy profiles, shown in parts a and b of Figure 2, for the approach of a single water molecule to the atoms 0 5 (a) and 0 7 (b), respectively, were calculated with water molecules oriented as shown in the picture. The hydrogen bonds were linear, and water molecules were lying in the plane of the HM molecule. In case (a), the water molecule was approaching oxygen atom 0 5 perpendicular to the line connecting the atoms01 and 0 5 . In case (b), the water moleculewas approachingatom 0 7 parallel to the line 01-05. energy barrier. In addition, a HF/4-31G calculation with inclusion of the solvent reaction field resulted in a comparable or even more favorable energy of the Z-conformer depending on the level of geometry optimization. Hydration of H M is interesting also in the context of surfaceactive substances, since the alkylated malonic acid has found wide use as a micelle-forming substance." Using a simulation approach, we investigate the influence of hydration on the rate of the intramolecular proton-transfer process. We use both classical transition-state theory and the quantum mechanical density matrix evolution (DME) approach introduced in the preceding article in this issue.12 The organization of this article is as follows: Section 2 describes the interaction potentials used in the simulations. Section 3 treats the proton-transfer rate by transition-state theory, followed by the calculation of the rate of tunneling using DME in section 4. Section 5 is devoted to the classical nonadiabatic proton transfer. Section 6 gives concluding remarks. 2. Interaction Potentials
Structure and numbering of H M are shown in Figure 1. For the simulations, the following interaction potentials were used: The first term represents the water-water interaction energy that was given in terms of the SPC (Simple Point Charge) model.13 The second term corresponds to the solute-water interaction and was modeled in terms of a 12-6-1 potential. The last term represents the intramolecular interaction energy. The methylene group C 3 was treated as a united atom. Three sets of atomic charges were used in the simulations (Table I). The first set of symmetric charge corresponds to values chosen to be close to the present GROMOS force field.I4 The other two sets of charges were obtained by van Duijnen (personal communications) using themethod preserving thedipole moment15 a t the ab-initio HF/DZP level of theory. A planar geometry was used in the ab-initio calculations, and two hydrogen atoms were added to the methylene carbon atom, treated in the simulations as a united atom, with standard bond lengths and angles. The geometry of the heavy atoms was kept in the C2" symmetry with the values of internal coordinates corresponding to the equilibrium values of bond lengths and angles given in Table 11. The hydrogen bond was treated as linear, which indeed
01 c2 c3 c4 05 06 07 H8
-0.6 0.4 0.0 0.4 -0.6 -0.5 -0.5 0.4 -0.82 0.41
-0.5808 0.5666 -0.0578 0.5436 -0.7245 -0.5942 -0.6382 0.4854 -0.82 0.41
-0,6570 0.5526 -0.0530 0.5526 -0.6570 -0.6173 -0.6173 0.4969 -0.82 0.41
900 898 2906 898 900 900 900 0.003610 793.3 0
23.25 23.65 46.63 23.65 23.25 23.25 23.25 0.0 25.01 0.0
OWf HWf a Symmetric charges chosen to be close to the original GROMOS force field. During proton transfer,values of this set did not changewith the 01-H8 distance. Chargespreservingthe dipolemoment of the HF/ DZP wave function for the asymmetrichydrogen bond; Le., H8 is at 0 1 . Charges preserving dipole moment of the HF/DZP wave function for the symmetric hydrogen bond, Le., Cb symmetry,or H8 is in the middle between 01 and 05. In (AI2 kcal mol-l)I/z. In (A6 kcal mol-l)I/z. /Values for the SPC water model. 0 Value for the 01-H8 and 05-H8 interactions in order to reproduce the in vacuo proton potential energy profile. In the OW-H8 pair interaction, this value is 0.01 1 44 (AI2 kcal in order to securea proper hydronium ion equilibriumgeometry. mol-1)1/2 TABLE Ik Bonding Parameters* Used in the Simulations bonds
bo
01x2 c4-05 C2-06 c4-07 C243 c344 01-05 01-H8
1.30 1.30 1.23 1.23 1.52 1.52 2.50 b
kb 1000 1000 1000 1000 900 900 1000 lo00 ~
angles
eo
ke
01-H8-05 C342-01 c344-05 0142-06 0544-07 C342-06 c344-07 C24344
180 116 116 125 125 119 119 118
50
~
impr dihedrals C3-06-0142 c3-05-07-c4
120 120 120 120 120 120 110
€0
kc
0 0
40 40
a The values of the harmonic force constants are in kcal mol-' A-2 for valence bonds and in kcal mol-' rad-2 for valence and improper dihedral angles. Bond 01-H8 was chosen as a reaction coordinate, and its equilibrium bond length was changed in the course of the TI simulation from bo = 0.9 A (state A) to bo = 1.25 A (state B).
is close to the fully optimized geometries of HM.9 The ab-initio calculations were performed for two positions of the proton. In the first the proton H8 was positioned 1 A from the proton donor 0 1 , and in the second it was positioned in the middle of the line connecting the atoms 0 1 and 0 5 . The former geometry is refered to as asymmetric and the latter as symmetric. The LennardJones parameters were taken from the standard GROMOS force field and were the same for all sets of charges. In order to assess the quality of the intermolecular potential functions, interaction energies of H M with a singlewater molecule were compared for a few configurations with ab-initio energies at the HF/DZP level of theory. No basis set superposition error (BSSE)correction was taken into account. For a detailed description of the geometries, see Figure 1. The resulting two curves are shown in Figure 2. The agreement is fairly good: at least the minima appear a t about the same hydrogen bond length. One should bear in mind that the ab-initio interaction energies
The Journal of Physical Chemistry, Vol. 97, No. 51, 1993 13471
Molecular Dynamics Simulation of Tunneling
0.00.l
' e
-4.00-
-6.00'
W
-8.00-10.00-
-12.00' 2.20
'
I
'
2.40
I
2.60
'
'
I
I
3.00
2.80
'
I
3.20
'
'
3.40
0
dOW-05) A CFiDZP
A
E W
0
9\
-2.00 O 7
nonequilibriumsolvation study the 01-05 distance was restrained with a realistic force constant of only 250 kcal mol-' A-2 in order to take fluctuations in the proton potential energy function associated with this internal coordinate into account. The resulting proton potential energy function was characterized by a classical energy barrier of 2.52 kcal mol-' and equilibrium OH distance of 0.98 A. The proton potential energy function in this form is parametrically dependent on the distance between atoms 0 1 and 05. By increasing the distance to 2.60 A, the barrier increases to 4.00 kcal mol-', while by decreasing the distance to 2.40 A the barrier drops to 1.21 kcal mol-'. When the asymmetric charges were used (in the free energy perturbation and in the nonequilibrium solvation), the proton was excluded from the intramalonate atom pair list, and the in vacuo proton potential was later added to the free energy differences. As a reaction coordinate, the distance was defined as 01-H8. Molecular dynamics simulations were performed by a vectorized version of the GROMOS program14 and ab-initio calculations by a modified HONDO program (ref !6 and references therein). 3. Thermodynamic Integration by the Method of Slow Growth
-lL,W
'
220
I
2.40
2.60
2.00
3.00
3.20
3.40
0
dOW-07) A
A
tF/DZP
0
EW
Figure 2. Comparisons of the interaction energies when water is approaching (a, top) 0 5 and (b, bottom) 0 7 in hydrogen malonate. For geometries,see Figure 1. Abinitiointeractionenergies(A)were evaluated on the HF/DZP level without BSSE correction, while in the empirical calculations ( 0 ) the charges preserving the dipole moment in the CI symmetry of hydrogen malonate was applied.
are overestimated due to the BSSE error, and thus the agreement between empirical and ab-initio potentials is even closer. Comparison of the interaction energiescalculated by the HM potential corresponding to the GROMOS symmetric charges (not shown) with the interactionenergiescalculated by abinitiocharges reveals that the former are at minimum energy about 1 kcal mol-' less favorable. Stouten et a1.l' applied the original GROMOS force field in a study of conformational equilibria of carboxylic acids in aqueous solution and obtained satisfactory agreement with experiment. We have used the ab-initio charges in the force field for the MD simulations from which quantities for the DME calculations were extracted. Equilibrium values and corresponding force constants for the bonding parameters are collected in Table 11. The equilibrium values correspond to the values obtained from the HF/6-3 1G** geometry-optimized structureg as the average between both carboxylic groups when the proton was at one side. Force constants were taken from the standard GROMOS force field. In deriving the intramolecular potential for HM, special care was taken to reproduce the abinitiocharges predicted by in vacuo proton potential. Since the GROMOS package does not provide the possibilityof including a specialfunctionalform of the potential functions, we modeled the intramolecular proton potential by a 12-1 interaction between the proton and atoms 0 1 and 0 5 using theGROMOScharges collected in Table 1. The latter two atoms were excluded from mutual nonbonded interaction and were instead harmonically restrained by using a force constant of 1000 kcal mol-' A-2 to the equilibrium value of 2.5 A in order to keep the 01-05 interatom distance fixed. Thus, the in vacuo proton potential energy barrier, which parametrically depends on the 0 1 - 0 5 distance, is kept nearly constant. However, in the
Reviews of the thermodynamic integration (TI) methodology have been given el~ewhere,'~.'~ so we will only discuss the characteristics of the simulations. One molecule of HM was placed at the center of a truncated octahedron with side toopposite side distance of 26.7 A. The truncated octahedron was filled with 255 water molecules. A twin cutoff method was used in the simulations for treatment of nonbonding interactions with the cutoff radii of 8 and 1 1 A, respectively. The pair list of the inner shell and the interaction with the outer-shell atoms were updated every 10th time step. MD was carried out with coupling to a bath with constant temperature (298 K) and pressure (1 bar) with coupling constants of 0.05 ps. This implies that our free energies are in fact Gibbs free energies. The leapfrog algorithm with the SHAKE procedure to constrain the molecular geometryZo was applied for water molecules. In the equilibrium TI studies, the masses of all hydrogen atoms were increased to 10 amu in order to be able to use a time step of 2 fs for the integration of equations of motion. The TI runs were preceded by 100 ps of equilibration for the GROMOS and ab-initio charges, respectively. In the case of GROMOS charges, the starting point for the TI (correspondingto the state A) was the structure where the 01-H8 distance was 0.9 A. The final structure (corresponding to the state B) had a 01-H8 distance of 1.25 A. The TI was performed forward and after 50 ps of equilibration in the state B backward over 100 ps in each direction. By a similar procedure, the difference in free energy using the ab-initio charges was estimated,except that thestarting structurewasonecorresponding to the minimum of the free energy estimated in simulations using the GROMOS point charges and that in the state B only 20 ps of additional equilibration was performed prior to returning to state A. The free energy along the postulated reaction coordinate can be referred to as the potential of mean force (pmf). The pmf for the proton transfer is shown in Figure 3a using GROMOS symmetric atomic charges and in Figure 3b using asymmetric charges, The free energy differences are collected in Table 111. Half the absolute difference in free energy for TI in the forward and backward directions was taken as an estimate of the error. The free energy difference between barrier and well is 4.15 i 0.10 and 2.75 f 0.14 kcal mol-l for the ab-initio and GROMOS symmetric charges, respectively. The differences in free energy of hydration, i.e., with the in vacuo contribution subtracted, are 1.63 f 0.10 and 0.23 f 0.14 kcal mol-', respectively. We treated the free energy differences also with analytical formulas for the spherical dipoldielectric continuum interac-
13472 The Journal of Physical Chemistry, Vol. 97, No. 51, 1993
TABLE IIk Free Energies Calculated Using T b e r m o d y ~ d c Integration by the Method of Slow Crowthe set d(Ol-H8)i.iti.lb d(Ol-H8)fi,1~ t / p A@ If 0.90 1.25 100 2.615 lb 0.90 1.25 100 2.892 2.75 &0.14c 2f 0.977 1.25 100 4.25 2b 0.977 1.25 100 4.05 4.15*0.1od Free energy differencesof hydrogen malonate in CZ,symmetry and the minimum corresponding to the CIsymmetry. All energy values are given in kcal mol-l and distance in A. Couplingparameter X w a d " such that the value of 0 corresponds to the initial state A and the value of 1 to the final state B. In the simulations of set 1, the charges corresponding to the GROMOS symmetric values (@ in Table I) were not changed with the coupling parameter X. In the simulations of set 2, the charges preservingthe dipole moment were used and were changed (from to @ in Table I) by changing of the coupling parameter A. The error was estimated as half the absolute differences between the values obtained in forward and backward thermodynamic integration.
I
4.5 f
*"I 3.5
4.5 2
I
1.00
1.05
1.10
1.15
1.20
Mavri et al.
1.25
d(O1-H8)A
The calculated barriers can be compared to ab-initio HF/43 1G calculations9with inclusion of the solvent reaction field.I0 The in-vacuo barrier of 2.74 kcal mol-' increases under the influence of the solvent reaction field by 2.24 kcal mol-', which results in an overall barrier of 4.98 kcal mol-'. The barrier increase of 2.24 kcal mol-' is in reasonable agreement with the TI result with ab-initio charges of 1.63 k 0.10 kcal mol-I. We note that a 4-3 1G basis set exaggeratesdipole moments, and hence predicted differences in free energies of hydration are expected to be exaggerated, as well. One can rationalize that the predicted barrier increase using GROMOS symmetric charges is too low. Therefore, further applications will refer to the free energy predicted by the ab-initio charges. From the barrier in terms of free energy, AG,one can calculate the rate constant for proton transfer using classical TST:
Figure 3. Free energy as a function of the reaction coordinate when the
symmetric form of hydrogen malonate is transformed to an asymmetric form in the course of the thermodynamic integration by the slow growth method. (k,top) Charges correspondingto the original GROMOS force field were applied and did not change with the coupling parameter. The solid smooth line is the energy profile in vacuo, the solid noisy line corresponds to 100 ps of thermodynamic integration in the forward direction, and the dashed line corresponds to 100 ps in the backward direction. The starting point for the thermodynamic integration was the valueofthereactioncoordinated(O1-H8) = 0.9A. (b, bottom) Charges preserving the HF/DZPdipole moment were applied and weredependcnt on the coupling parameter. The solid smooth line is the energy profile in vacuo, the solid noisy line corresponds to 100 ps of thermodynamic integration in the forward direction, and the dashed line corresponds to 100 ps in the backward direction. The starting point for the thermodynamic integration was the value of the reaction coordinate d(01-H8) = 0.977 A, correspondingto the free energy minimum in the calculation using symmetric GROMOS charges. tion.21 Since the free energy differences are known from the TI and the change of the dipole moment associated by the transformation of asymmetric H M to its symmetric form can be calculated from the atomic charges, the radius of the spherical cavity can be estimated. From the free energy barrier (excluding the in vacuo contribution using ab-initio atomic charges) and the dielectricconstant oftheexternalcontinuum (e = 78.5), thecavity radius was estimated to be 2.17 A. We applied this cavity radius together with the change of the dipole moment calculated using GROMOS symmetric charges to estimate the free energy difference on the basis of a continuum model. The predicted difference in free energy is 0.16 kcal mol-], which is in good agreement with the TI difference of 0.23 f 0.14 kcal mol-'. Obviously, despite some unfavorable system characteristics for the free energy study (nonzero net charge, considerable dipole moment,strong hydrogen bonds) that usually result in large errors manifesting themselves in significant hysteresis, convergence of the free energy calculations seems to be satisfactory.
For the classical free energy barrier of 4.15 f 0.10 kcal mol-', the rate constant is (5.6 f 1.0) X lo9 s-l. TST requires the motion of the proton to be slow enough to be in equilibrium with the solvent for the whole range of the reaction coordinate and assumes a transmission coefficient equal to one. This is in contradiction with the common view of the fast proton movement in a hydrogen bond. In the next section, we consider the quantum mechanical tunneling motion of the proton. In section 5 , we shall return to the classical motion over the barrier and consider the nonadiabatic case.
4. Deosity Matrix Evolution
As pointed out already by Borgis and H y n e ~ , 2 the ~ - ~motion ~ of the proton in hydrogen bonds is considerably faster than the motion of the surrounding solvent. The appropriate physical picture of such a process would be that the proton moves in the in-vacuo potential distorted by the solvent. Since the motion of the light proton is considerablyfaster, the proton "feels" the solvent as being frozen, at least for the slow components of the polarizability. The configuration of the solvent changes with time, and the proton potential energy function changes with it. Tunneling of the proton is very likely to occur when the reactant (R) and product (P)levels are of about the same energy. The occurrence of a state in which R and P levels are of the same energy is a rare event; in other words, it is necessary to wait a considerable amount of time before a suitable solvent fluctuation occurs. From thesimulations,the probability densityof theenergy splitting can be evaluated, as well as the associated free energy. Thedifferencebetweenthesocalculated freeenergycorresponding to zerosplittingand freeenergycorrcapondingtothemost probable splitting equals the reversible work required to reorganize the solvent is such a way that tunneling is likely to occur. It is expected
Molecular Dynamics Simulation of Tunneling
The Journal of Physical Chemistry, Vol. 97, No. 51, 1993 13473
that the difference in free energy of hydration between the asymmetric and symmetric state obtained by thermodynamic integration in the previous section is a good estimate for the cost in free energy of solvent reorganization to a state with zero splitting. HM with the proton positioned between 0 1 and 05 favors solvent configurationswith zerosplitting, or, more precisely, such configurations are the most probable ones because of symmetry. In previous ~ ~ r kthe,rate ~ of~ tunneling , ~ ~ was considered only starting from a fully populated R state, with the solvent relaxed with respect to the R state, although in later simulation^^^ the proton was allowed to move. With strong solvent stabilization, the procedure may not yield efficient statistics since the sampling of states with high tunneling efficiency may be very sparse. Therefore, we extended our study to the case where R and P states are equally populated, which restricts the solvent phase space to configurations with zero or small splitting. The higher tunneling rate corrected by the free energy difference required for such solvent reorganization should yield the same result for the rate constant of tunneling. 1. Calculation of the Proton Potential Energy FMC~~OM Including Solvent Contributions. Borgis and HynesZ2introduced two quantities: splitting M ( t ) = ( + p l d + p ) - (+R(II+R),and coupling, c ( t ) = (+#&), where the +p and +R refer to the ground state wave functions of the proton in the product and in the reactant well, respectively. The splitting M(t)equals the classical potential energy difference M,lass(t)between the product and reactant states, assuming that only ground vibrational states are populated but corrected for zero point energy. Two MD runs were performed for 100 ps each, for fiied 05H 8 distances of 1.OO and 1.25 A, respectively. The equations of motions were integrated by a time step of 1 fs. Coordinates of theatomswerewrittentodiskevery2fs,suchthat50 000samples were available for the analysis. The analysis revealed that the stretching frequency associated with the 01-05 vibration had a value of about 1300 cm-l, which is a relatively high value even for such a strong intramolecular hydrogen bond. For each set of the coordinates, the proton potential was calculated by moving the proton from the R well to the P well for 21 equidistant points between d(OS-H8) = 0.7 and 1.8 A. In the proton potential energy calculation, we changed all the HM atomic charges as a function of the progressing reaction coordinate. The charges were interpolated for thevalues of the reaction coordinate between the R and P free energy minima, while values of the asymmetric ab-initio charges were used beyond these values. A parabolic interpolation was used-being the simplest even polynomial through three points. The in vacuo proton potential energy function, parametrically dependent on the 0 1 - 0 5 distance, was then added to the solvent contribution. The average classical splitting obtained was (MclBss) = 7.65 kcal mol-'. The change of all atom charges when moving the proton is essential, since otherwise the average value of the splitting appeared to be ( M )= 2.61 kcal mol-' compared to the value of 6.78 kcal mol-' obtained by changing the charges (see next subsection). In previous the problem was circumvented by increasing the charge of the proton. 2. Solving tbe Time-Imlepedent SchrWnger Equation. The time-independent Schrgdinger equation was solved in one dimension using a vibrational method with two Gaussians as basis functions, centered around both in vacuo minima and with a standard deviationof 0.0935 A. Kineticenergy, potential energy, and overlap integrals were calculated, and orthogonalization of the basis set was performed by the Schmidt procedure. The potential energy integral was calculated numerically by the Romberg integration procedure, using cubical splineinterpolation between the 21 points for which the proton potential energy was obtained as described in the previous subsection. After diagonalization of the Fock matrix, thevalues of coupling and splitting
o'
2
0
-2
6
4
8
10
12
14
16
& E kodlmol
Figure 4. Probability density of splitting w(AE) calculated from the
100-psimulation. The distributionwas calculated from 50 000samples. The dashed curve is a Gaussian approximation for the reactant state, while the dotted curve is the corresponding Gaussian for the product state.
-2
0
2
4
6
8
IO
I2
14
A € kcal/mol
Figure 5. Free energy profile corresponding to the probability density ) shown in Figure 4. AG(AE) = - k T l n ( w ( U ) ) kTln(w(( M ) ) (solid line). The h h e d curve correspondsto a harmonic approximation, while
+
the dotted curve corresponds to the product state. The latter curve was obtained by mirroring the reactant harmonic curve over the y axis.
were obtained simultaneously. Unlike in ref 22, it was not necessary to include a parametrical dependence of the coupling on the 0 1 - 0 5 distance. The average value of the splitting ( M ) obtained by solving the Schrgdinger equation was 6.78 kcal mol-', which can be compared to the classical value of the average splitting when zero point energies are not taken into account of (Mcl,)= 7.65 kcal mol-'. Borgis and Hynes22only took the classical splitting into account. The Schrgdinger equation was solved for each configuration as described above, and the values for coupling and splitting were tabulated for every 2 fs of the trajectory. For M(r),the normalized distribution w ( M ) was calculated (Figure 4). The associated free energy change AG = kBT ln(w) keT l n ( w ( ( M ) ) ) is shown in Figure 5. The corresponding P solvent curve was obtained by mirroring over the coordinate origin hE = 0. The R and P curves intersect at the height of AG = 1.6 f 0.1 kcal mol-'. For a harmonic solvent, the relation AG = (AE)/4 should hold,22 and the value for AG from this relation is 1.69 kcal mol-'. The difference in free energy of hydration between symmetric and asymmetric HM obtained from TI resulted in a reorganization free energy of AG = 1.63 fO.10kcalmol-I. Agreement betweenthefreeenergissofsolvent reorganization calculated by all three methods is quite good. This means that the restriction to solvent configurationsin equilibrium with symmetric HM indeed selects the relevant configurations for efficient tunneling with AE = 0.
+
Mavri et al.
13474 The Journal of Physical Chemistry, Vol. 97, No. 51. 1993
1; V
I
8
o
.
O
O
l
~
1
~
\-/------
-0.m
h
I I
v"
-0.001 -0.002
0.999995 0
10
20
30
40
50
1
-0.003'. 0
60
W . '
'
.
.
.
.
.
I
.
t fs
30
20
t fs
Figure 6. (z(r)) as a function of time, calculated from the simulation with the proton fixed at 0 1 (asymmetric position).
3. Solving the Time-Dependent Scbrijdinger Equation by Density MatrixEvolution. Thedensity matrixevolution approach is described in detail in the accompanyingarticle.I2 In this article, we will adopt the following form of the differential equations, suitable for numerical solution: (3) 1 j = *Im(pRp)= j-(2c(t)z + x U ( t ) )
........I........"
10
Figure 7. (x(r)) (solid line) and ( y ( r ) ) (dashed line) as functions of time, calculated from the simulation with the proton fixed at 01 (asymmetric position).
1 \
I
0.99999
h
=v
0.99998
(4) 0.99997
0.99996 0
Note that the quantities w, and wx in the preceding paper are related to the splitting M(t)and the coupling c ( t ) by
A E / h = -0,
(6)
2 c ( t ) / h = 0,
(7)
and 0, = 0. We note that diagonalization of the Fock matrix is not needed to obtain hE, but it has no numerical effect since the off-diagonal element is 3 orders of magnetude smaller than AI?. The vector (x,y,z) represents the essential components of the density matrix: z measures the relative population of the R and P states with z = 1 corresponding to the pure R state and z = -1 to the pure P state. Starting with the pure R state as initial condition (x,y,z)(t = 0) = (O,O,l), the phenomenological rate constant for the proton transfer by tunneling from the R well to the P well can be related to the averaged value of the initial decay rate of z ( t ) . As was shown in the preceding article,I2 the rate constant is related to the average decay rate of z(t) on a course-grained time scale: k = 1/2Az(t))/At. The factor '/2resultsfrom thedefinition of z. The system of three differential equations (eqs 3-5) was solved with a time step of 0.04 fs for times between 0 and 60 fs. Thevaluesof M ( t ) andc(t) were interpolatedbetween thevalues 2 fs apart, obtained in the previous subsection, using a cubical spline procedure. Averaging was performed over the entire set of M ( t ) and c ( t ) values obtained from the MD simulation by choosing every 4 fs of the trajectory as a new starting point for a density matrix evolution over 60 fs. The system of differential equations was integrated about 25 000 times using a fourth-order Runge-Kutta method. Thedifferential equationswere also solved with the same procedure using the data from the MD simulation with a symmetrically placed proton. The time evolution of the variables ( z ( f ) ) , ( x ( t ) ) , and ( y ( t ) ) for the asymmetric case is shown in Figures 6 and 7 and of ( z ( t ) ) for the symmetric case in Figure 8. The function ( z ( t ) ) in the
10
20
30
40
50
€0
t fs
Figure 8. ( z ( r ) ) as a function of time, calculated from the simulation with the proton fixed at the center of the line connecting the atoms 01 and 05 (symmetric position).
asymmetric case is characterized by a rapid decay in the first 6 fs, followed by a small increase and subsequently by an almost linear decay from 15 to 60 fs. A similar initial behavior of the probability was reported by CartlingZSfor a classical particle in a bistable potential coupled to an external bath by solving the Fokker-Planck equation. From the final slope of the decay, the rate constant was calculated. ( x ( t ) ) and ( y ( t ) )show an initial significant oscillation followed by constant values. In the case of symmetric HM, the decay of z ( t ) is considerably faster than that in the asymmetric case (Figure 8). In contrast to the asymmetric case, the proton in the symmetric case exhibits monotonic behavior of z ( t )at small values of time. The increased rate relative to the asymmetric case is easily understood, since the solvent in this case adopts configurations that are more favorable for tunneling. The thermodynamicweight of such states can be calculated via the free energy difference that is available from our equilibrium solvation study. Thus, we have two alternative ways of computing the quantum mechanical nonadiabatic proton transfer from R to P: by computing improbable transfer from the full solvation ensemble of asymmetric HM or by computing the probable transfer in a selected solvation ensemble of symmetric HM. In the next subsection, the results will be compared. 4. Rate Constant for the Proton Transfer by Tunneling. Rate constants of tunneling were calculated as one-half the coarsegrained derivative (z)with respect to time on the time interval between 20 and 60fs for the asymmetrical H M and on the interval between 10 and 60 fs for the symmetrical HM. The same results can be obtained graphically from Figure 6 and Figure 8. The calculated rate constants of tunneling are (1 -36f 0.05) X lo7 and (2.67 i 0.05) x 108 s-] for the asymmetric and symmetric
Molecular Dynamics Simulation of Tunneling
The Journal of Physical Chemistry, Vol. 97, No. 51, 1993 13475
16
/I
12 141
-2 0.8
1.0
1.1
1.2
1.3
d(01-HE)
1.4
1.5
1.6
B
Figure 9. Averaged nonadiabatic proton potential energy function
calculation which includes the zero point energy difference betweenthePandRstates. WeappliedTSTandused theaverage splitting as the energy difference, assuming that the P state is a transition state, which is an extremely crude approximation. In addition, we assumed that in all cases the proton remains in the P state if it reaches this state (which is not true) and obtained a rate constant of 6.6 X lo7 s-I. This value is somewhat larger than the rate of tunneling, but one can expect a low transmission coefficient and therefore a significant lowering of the TST rate constant. We emphasize at this point that in the above calculations of proton potential energy functions the electronic polarizability was not taken into account. Since it is expected that electronic polarizability follows the proton motion, some additional stabilization of the transition state and product state is lost by this neglect.
calculated from the simulation with the asymmetricallyplaced proton.
6. Discussion
case, respectively. The Boltzmann factor for the latter case is exp(-AG/keT) = (6.4 f 1.0) X where AG = 1.63 A 0.1 kcal mol-l is the free energy difference for the solvent reorganization. Multiplied by this Boltzmann factor, the rate constant of tunneling derived from the symmetric case is (1.7 f 0.3) X lo7 s-I. This is in good agreement with the value directly determined from the asymmetric case. In addition, we calculated the rate of tunneling from the analytical formula (eq 5.17 in ref 22):
where (c) refers to the average coupling in the asymmetric state and (AI?)refers to the average splitting. Using our values (c) = 4.69 X kcal mol-I and (AE)= 6.78 kcal mol-’, the calculated rate constant was found to be (7.3 f 2.4) X lo7 s-l, which can be compared to the DME value of (1.36 f 0.05) X lo7s-I. One should bear in mind that in deriving the analytical expression several approximations were made.22 5. Classical Nonadiabatic Proton Transfer
The DME approach gives only the contribution of tunneling to theoverall rateconstant. If the quantumparticleis transformed to a classical particle, e.g. by increasing its mass, this contribution vanishes. The proof is straightforward, since in eq 5 dz/dt is proportional to the coupling which rapidly approaches zero for species heavier than the proton. Therefore, additional calculations using corrected TST are necessary. Classical TST assumes that the reacting system is in equilibrium with the environment for all values of the reaction coordinate. Due to the fast motion of the proton, it is expected that the solvent will not be in equilibrium with all its components of polarization. This, even in the classical case, leads to an additional increase of the barrier height and to an asymmetric shape of the proton potential. We shall estimate the classical nonadiabatic rate of proton transfer and compare this rate with the DME value. Proton potential energy functions have been calculated in order to solve the one-dimensionalSchrBdinger equation, and the 50 000 functions have been stored on disk for the asymmetric case. Averaging was performed, and the obtained average proton potential energy function (V)is shown in Figure 9. Instead of by a transition state, the average proton potential energy function is characterized by an inflection point. In other words, the point at which the transition state is expected is more favorable than the product state. The average classical energy of the expected transition state is 6.13 kcal mol-I, while at the product site it is 7.65 kcal mol-’. Both values are relative to the R state. An average splitting of 6.78 kcal mol-’ is available from the DME
The present calculations show that it is possible to calculate the rate constant of proton tunneling by the method of DME. The calculated rate constant is 5 times lower than the analytical results of Borgis and Hynes (eq 5.17 in ref 22), which can be considered as a reasonable agreement in view of the approximations in the analytical theory. Nevertheless, there are three weak points of the present calculation. The first is that we did not take into account the vibrationally excited states of the proton, despite the large value of the splitting. The DME method would allow the inclusion of excited states at the expense of a higher dimensionality of the density matrix. The second weak point of the DME calculation of the rateof tunneling relative to thecumulant expansion method proposed by Borgis et al.22323 is that 0-0vibration is in the present case treated as classical oscillator. This is not likely to present a serious error. It could be redressed by including the 0-0 vibration in the quantum mechanical part of the calculation. The proton degrees of freedom orthogonal to the arc along the 0105 line, which is treated quantum mechanically, should in principle also be treated as quantum degrees of freedom. This implies solving the three-dimensional SchrGdinger equation to obtain the proton potential energy functions, instead of the onedimensional one. The third weak point is that the DMEdynamics is not appropriately coupled to the classical dynamics of the environment,and hence the dynamical response of the environment to a proton transfer is not included. In the present calculations, it was essential to change simultaneously atomic charges of all atoms when moving the proton in the calculation of the potential of mean force as well as in the calculation of nonadiabatic proton potential energy functions. Especially in the latter case, the derived quantities were qualitatively dependent on the values of the charges. For thecalculation ofthe rateconstant from MD, weevaluated for every snapshot the value of splitting and coupling by solving the one-dimensionalSchriidinger equation without any additional simplifications and assumptions. Solving the SchrBdinger equation for every snapshot is necessary, since it was shown that the zero point energies of the proton are not the same in the P and R well. Approximation of the splitting by a classical energy difference22is therefore too crude for the present case. In addition, using our approach no parametric approximation of the coupling was necessary, as in ref 22. We found that the difference in free energy of solvation (i.e., the free energy difference with the in vacuo contribution subtracted) between the symmetric and asymmetric positions of the proton has the meaning of a reorganization free energy of the solvent required to obtain zero splitting. Within thecomputational errors, we found this free energy in agreement with the estimated
13476 The Journal of Physical Chemistry, Vol. 97, No. 51, 1993
value of (hE)/4. Evaluation of the tunneling rates using a cumulant expansion requires a harmonic free energy associated with the probability density of splitting. In the present case, perfect harmonicitywas not observed (see Figure 5). One expects even less harmonicity when replacing the aqueous environment by the fluctuating environment of the macromolecule, e.&, enzyme. Rate determination by DME has, on the other hand, no such requirements with respect to the shape of the free energy curve. Classical proton transfer over the barrier is an alternative reaction pathway. Averagenonequilibriumproton-transfer c u m a calculated by assuming solvent to be frozen give evidence that proton transfer in HM is not a transition-state process. The roughly estimated TST rate constant was shown to be somewhat larger than the tunneling rate, although it is believed that the transmission coefficient would be small and reduce the TST rate considerably. HM in aqueous solution is an example of a highly polar and charged hydrogen-bonded system (as most of the biological systems are). We expect treatment of the related enzymatic reactions where proton transfer is the rate-limiting step, as well as further progress in the methodology of quantum simulations in the near future.
Acknowledgment. We are grateful to Dr. P.Th.van Duijnen for performing the ab-initio calculations. One of us (J.M.) is grateful for financial aid from NUFFIC and to the University of Groningen for the hospitality during his stay in the Netherlands in 1990 and 1992.
Refereaees and Notes (1) Chandrasekhar, J.; Smith, S.F.; Jorgensen, W. L. J . Am. Chem.Soc. 1985, 107, 154.
Mavri et al. (2) Chandler, D. J. Srat. Phys. 1986, 42, 49. (3) Bergsma, J. P.; Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J . Chem. Phys. 1987.86, 1356. (4) Kurz, J. L. J. Am. Chem. Soc. 1989,111,8631. (5) Daniel, D. H. Mc.; Brown, H. C. Science 1953, 118, 370. (6) DjinoviE, K.; GoliE. L.; Had%, D.; Orel, B. Croar.€him. Acra 1988, 61,405. (7) Maury, L.; Bardet, L. J. Raman Specrrosc. 1978, 7, 197. (8) KidriE, J.; Mavri, J.; Podobnik, M.;Had%, D. J. Mol. Srrucr. 1990, 237, 265. (9) Mavri, J.; Hoddkk, M.; Had%, D. J. Mol. Srruct. (THEOCHEM) 1990, 209,421. (10) Miertd, S.;Scrocco, E.; Tomasi, J. Chem. Phys. 1981,55, 117. (11) Vikingstad, E.; Saetersdal, H. J . Colloid Inrerj. Sci. 1980,77,407. (12) Berendsen, H. J. C.; Mavri, J. J . Phys. Chem., preceding paper in
this issue. (13) Bcrcndsen,H. J.C.;Postma, J.P.M.;Gunsteren, W. F.van;Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; p 331. (14) Gunsteren, W. F. van; Berendsen, H. J. C.; GROMOS,Groningen Molecular Simulation Program Package, Biomos b.v., Nijenborgh, 4, 9747 AG Groningen, The Netherlands, 1987. (15) Thole, B. T.; Duijnen, P.Th. van. Theor. Chim. A d a 1980,55,307. (16) Dupuis, M.; Watts, J. D.; Villar, H. 0.;Hurst, G.J. B. Comp. Phys. Commun. 1989,52.415. (17) Stouten, P. F. W.; Kroon-Batenburg, L. M. J.; boon, J. J . Mol. Srrucr. (THEWHEM) 1989, 200, 169. (18) Gunsteren, W. F. van. In Compurarion ofFree Energyfor Biomo-
lecular Sysrems; Gunsteren, W. F. van, Weiner, P. K., Eds.; Escom Science Publishers: Leiden, The Netherlands, 1989; p 27. (19) Mezei, M.;Beveridge, D. L. Free Energy Simulations. Ann. Acad. Sci. N.Y. 1986, 482, 1. (20) Ryckaert, J. P.; Ciccotti, G.;Berendsen, H. J. C. J. Comp. Phys. 1977, 23, 327. (21) Kirkwood, J. G.J . Chem. Phys. 1934, 2, 351. (22) Borgis, D. C.; Hynes, J. T. J . Chem. Phys. 1991, 94, 3619. (23) Borgis, D. C.; Lee, S.;Hynes, J. T. Chem. Phys. Lcrt. 1989,162.19, (24) Borgis, D. C.; Tarjus, G.;Azzouz, H. J. Chem. Phys. 1992,96,3188. (25) Cartling, B. J. Chem. Phys. 1989, 90,1819.
