J. Phys. Chem. 1993,97, 6907-6913
6907
Stochastic Dynamics Simulations of the Alanine Dipeptide Using a Solvent-Modified Potential Energy Surface Paul E. Smith and B. Montgomery Pettitt' Department of Chemistry, University of Houston, 4800 Calhoun Road, Houston, Texas 77204-5641
Martin Karplus' Department of Chemistry, Harvard University, I2 Oxford Street, Cambridge, Massachusetts 02138 Received: November 18, 1992; In Final Form: March 31, I993
Langevin dynamics on a potential of mean force energy surface is used to model the effects of aqueous solvent on the structure and dynamics of the alanine dipeptide. Conformational transition rates obtained by correlation function analysis and hazard plots from several simulations are compared. In particular, averages obtained over three 25-11s runs are compared with a single run of 100 ns. Rate constants for selected conformational transitions are also examined. For the conformational processes considered here (ato 0 and the reverse), we observe two rates, one of which is only significant for simulations of more than 25 ns. On further decomposition of the rate process, it is shown that the two observed rates correspond to the individual rates for rotation around the 9 and rl. dihedrals. It is also shown that, for the conformational transitions investigated here, the transition process corresponds to an essentially uncorrelated motion of the two dihedral angles.
I. Introduction A number of advances in the understanding of the chemical function of biological molecules have come directly from the analysis of dynamical simulation^.^-^ However, many important events (e.g.,conformationaltransitions) are characterized by long time scales which require trajectories beyond those currently available for large systems (e.g.,a protein surrounded by water). In an effort to access longer time processes it is necessary to reduce the complexity of the calculation, and therefore early work typically dealt with the dynamics of macromolecules in vacuo. Although recently it has become common to simulate the system of interest surrounded by an explicit solvent en~ironment,~.~ most of the simulationshave been limited to subnanosecondtime scales. Polyatomic solutes with internal degrees of freedom display a rich variety of solvation behavior, especially in a strongly associating liquid such as water. There are two types of solvent effects. One is an equilibrium effect that corresponds to altering the thermodynamics of the solute by modifying the effective potential energy or potential of mean force. Such a change in the potential also alters the dynamics. For peptides there is a pronounced shift in the relative conformational populations on changing solvent conditions due to the competition between hydrogen bonding and packing requirements.5 Hence, it is necessary to include these solvation effects in a reliable way. The second solvent effect is a purely dynamical effect on the solute due to collisions with surrounding solvent molecules. If the solvent is to be treated explicitly,enough solvent molecules must be included to ensure that the solute molecule is out of range of its own images when using periodic boundary conditions. Ideally, it is desirable to study the solute internal structure and dynamics in a solution sufficiently dilute that one may safely ignore finite concentration phenomena. Computationally, this has meant surrounding the solute molecule with a minimum of 4-5 layers of solvent,2 although such a system is still of the order of 0.05 M, which is large compared to experimental conditions. Since the complexity of a dynamic simulation increases rapidly with the number of sites, a simulation with a multisite model for water as the solvent (TIPS: SPC? etc.) can be many times (101103) more expensive than a simulation without explicit solvent. In this paper, we present a methodology for reducing a dynamical many-body solvation problem to an effective few-
body simulation and apply this technique to study conformational transitions of the blocked alanine dipeptide. The technique is similar in spirit (although different in detail) to an earlier approach used to study the trans to gauche isomerization rate constant for butane, and other simple alkanes, at infinite dilution in a semiempirical waterlike monatomic solvent? In that study, a stochastic dynamics formulation in the diffusive limit was used to model the many-body kinetic effects of the solvent, and a MacMillan-Mayer level potential of mean force (PMF) from thePratt-Chandler theory9was used to model the effectivesolventmediated internal potential. Here, we use a Langevin dynamics scheme that is not restricted to the diffusive regime; that is, we do not ignore inertial terms in the equations of motion. The intramolecular potential of mean force is derived by use of the superposition approximation from solutionsof Ornstein-Zernikelike (extended RISM) integral equations for the equilibrium structure of charged solute sites in an associated polar molecular fluid."J As the model solute, we have chosen the blocked alanine dipeptide, N-acetylalanyl-N'-methylamide,which has sufficient complexity to include solvation effects common to larger systems. In particular, it has the polar groups (C=O and N-H) of the peptide bond and methyl groups representative of nonpolar side chains. This molecule has previously been studied by molecular dynamicssimulationsin vacuum11J2and aqueous solution.l3Also, studies of the effective potential for the variation of the central 4 and 9 angles (Ramachandran plot) have been made in vacuum14J5 and in aqueous ~ o l u t i o n . ~Thus, ~ J ~ we shall be in a position to compare the results of our technique with previous theoretical results. It is not the purpose of this paper to compare the current potential model with the experimental data for this system. Also, this is not possible because equilibriumand dynamic information for this model system are very limited.'* We consider rate processes into and out of regions in the space that are related to the idealized helical and pleated sheet regions for polypeptides which we call the CY and 6 states, respectively. For our system, and indeed for any polypeptide system of less than 10 or 12 residues, such a classification is artificial as helices and sheets rarely are the stable state forming. However arbitrary, the state designations for short peptides are often used and provide a conceptual framework.
0022-3654/93/2097-6907$04.00/0 0 1993 American Chemical Society
Smith et al.
6908 The Journal of Physical Chemistry, Vol. 97, No. 26, 1993
In section I1 we describe the method used to reduce the manybody problem to an effective few-body problem. Section I11 describes the molecular model. In sections IV and V respectively we present the results and conclusions of our calculations.
11. Method Consider an infinite number of particles interacting through a pair potential. (The generalization to three-body and higherorder potentials is straightforward but not necessary at this point.) The dynamical equations of motion for the system may be written as an infinite set of coupled differential equations,
j= 1
where the prime on the sum indicates that i = j is excluded. Here V,,(r,,) is the pair potential, and 7, and mrare the acceleration and mass of particle i, respectively. If we let i = 1, N refer to the N particles of the solute and i = N 1, refer to solvent, we have the exact equations of motion for an infinitely dilute solute composed of N particles. We can formally reduce the number of equations to a finite set by introducing a Langevin-type approach. With W representing the solvent-mediated effective intramolecular potential, or potential of mean force (see below), we may write the solute dynamics in terms of a generalized Langevin equation (see ref 19 and references therein),
+
N
mir, = - x t V W , j ( r , j ) - m,Kdt’K,((r],r? h,(r-t?
+
j=1
AiW,r)
i = 1, N (2)
where A,((r),t)is a position and time-dependent random force that has the characteristics of the collisional input of energy into the solute from the solvent, it is the velocity of particle i, and the integral characterizes the collisional sink for energy provided by the solvent. The factor K,((r),?’)in the integral is the generalized friction kernel, where (3represents the set of position vectors of all particles. The generalized friction kernel and the random force are related by the fluctuation-dissipation theorem.20 Since it depends on the positions of all the solvent particles through (r), and their past history, the full generalized Langevin equation requires solving the complete set of equations of motion given by eq 1 and so does not directly simplify the problem. To avoid this difficulty, various approximations for the friction kernal have been introduced (see ref 21 and references therein). In this study we have adopted one of the more simpler approaches by choosing to model the solvent effect on the kinetics of the solute by means of the Langevin equation, N
i = l , N (3) where y, is the friction constant (Le., the generalized friction kernal is taken to be delta correlated in time and independent of atomic positions) and At(?)is the random forceon a given particle, also assumed independent of atomic positions. In our case, Ai(?) is assumed to have a Gaussian distribution with moments, ( A , ( t ) )= 0
(4)
and where i andjlabel different solute sites and kBTis the Boltzmann constant multiplied by the absolute temperature. We are interested in the intramolecular structure and dynamics of a solute molecule in a given solvent. Since we are not primarily concerned with overall translational diffusion, we have neglected hydrodynamic interactions in the model. For the choice of friction
TABLE I: Parameters Used in the Stochastic D ~ M I U ~ C S Simulations atom name q, le1 c, kl/mol u, A y, ps-1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
CL HLl HL2 HL3 CLP OL NL HL CA HA CB HB1 NB2 NB3 CRP OR NR HR CR HRl HR2 HR3
-0.0480 0.0246 0.0246 0.0246 0.4800 -0.5410 -0.3210 0.2960 0.0700 0.0246 -0.0480 0.0246 0.0246 0.0246 0.4800 -0,5410 -0,3210 0.2960 -0,0480 0.0246 0.0246 0.0246
0.3768 0.0188 0.0188 0.0188 0.3768 0.9630 0.6699 0.0188 0.3768 0.0188 0.3768 0.0188 0.0188 0.0188 0.3768 0.9630 0.6699 0.0188 0.3768 0.01 18 0.01 18 0.0118
3.208 2.616 2.616 2.616 3.208 2.640 2.770 1.604 3.208 2.616 3.208 2.616 2.616 2.616 3.208 2.640 2.770 1.604 3.208 2.616 2.616 2.616
84.38 0.00 0.00 0.00 40.25 45.89 45.73 0.00 50.19 0.00 84.38 0.00 0.00 0.00 40.25 45.89 45.73 0.00 84.38 0.00 0.00 0.00
coefficients we have used a slight modification of the Pastor and Karplus accessible surface area (ASA) frictional drag In this picture,21as in a number of other^,^^^^^ the translational drag coefficientsare calculated using Stokes’law for slip boundary conditions, Yi = 4 w , / m ,
(6)
where 7 is the viscosity of the pure solvent, mi is the mass of the solute particle, and a, is the hydrodynamic radius. In the ASA model, the hydrodynamic radius ai is chosen to reproduce the exposed surface area (A,) of the ith site in the solute, Le., a, = (A,/47r)”2
(7)
This method has been shown to produce reliable translational diffusion coefficients for a number of different molecules.21For the traditional extended atom model of hydrocarbon chains, it has been shown to be equivalent to the bead model.21 However, this picture must be modified when hydrogen atoms are treated explicitly rather than being incorporated into the heavy atoms to which they are attached. Our approach is to model the hydrogen atoms in the rotational slip limit, Le., YH = 0. This has been shown to be reasonablefor hydrogens involved in rotational motion around high-symmetry axes, such as the rotation of methyl hydrogens around the C3 axis or the rotation of benzene ring hydrogens around the Ca symmetry axis.21 We use the ASA calculated from an extended atom model25 and assign all of the drag to the heavy atom, leaving the explicit hydrogens with zero drag. The resulting drag coefficients for each of the atoms in the dipeptide are shown in Table I. Next we must consider the problem of determining the effective solvent modified intramolecular potential, W. To simplify the analysis, we assume spherically symmetric pairwise additive atomic site-site potentials; i.e., we approximate the solute n,point internal distribution function by a superposition of atomic pair terms. This is similar in spirit to Kirkwoods’ superposition approximation to the triplet distribution function.2* That a p proximation has been shown to give good qualitative results when compared to ~imulations.2~ However, it has been demonstrated that the superposition approximation breaks down for highly overlapping ~ p h e r e s .This ~ can be rationalized in terms of an overcounting of the mutually excluded volume. In our approximation we only use cavity correction terms for the interactions between nonoverlapping spheres (1,4 neighbors or higher) in an
The Journal of Physical Chemistry, Vol. 97, NO. 26, 1993 6909
Stochastic Dynamics of the Alanine Dipeptide effort to alleviate this overcounting problem. This is consistent with the use of a pair potential that excludes nonbonded interactions for 1,2 and 1,3 neighbors that are connected by stiff bonding terms. We use the notational convention that 1~ represents the set of vectors (r& ...,r$"$) where n M is the total number of sites in a molecule of type M and r):) denotes the position of the ath site in the lth molecule of type M. As previously we approximate the one molecule (M) nu-point cavity distribution function (exponential of the solvent contribution to the PMF) by a superposition of two-site, atomic particle, cavity distribution functions, ~*,(l,)
=
yaMyM(r)
180.
120.
60.
*
-60.
(8)
- 120.
amM
where y,,,(r) is the cavity distribution function matrix element for sites a and y in molecule M . The product is taken for intramolecular pairs separated by three bonds or more. To complete this picture, we must have a theory to compute the relevant distribution functions, yaH,(r), for the general case of charged particles in a polar molecular solvent. This is equivalent to finding the PMF for each pair of sites in the that are at least three bonds apart. This process would be eXtrf"ly tedious by molecular dynamicssimulations. The availabilityof an analytic technique that is readily solved by elementary numerical methods is therefore important. Specifically,we obtain the necessarycavity distribution functions from the Ornstein-Zernike-like RISM equation (XRISM) in renormalized form with a hypernetted chain closure. The details of this procedure have been published elsewhere, and we refer the interested reader to this i i t e r a t ~ r e . ~ The ~ . ~results ~ . ~ ~ for the effective solvent-mediated potential surface have been shown to be in reasonable accord with explicit simulations of the free energy solvation surface for the alanine dipeptide.30 This approximationto the intramolecular pair potential of mean force, when coupled with an effective dynamical theory such as Langevin dynamics, represents an approximate reduction of a many-body Hamiltonian problem to an effectivefew-body problem and thus a substantial computational saving.
0.
-180. -180.
-120.
-60.
0.
60.
120.
180.
9 surfam for the alanine dipeptide, Contoun are separated by 1 kcal/mol, The lowest antour is darkend, the highest antoum are dashed. Figure
Solvent-modificd
TABLE II: Averages from the Stochastic Dyllpdcs Runs' run 1 length of run total energy potential kinetic bond angle dihedral Urey-Bradley
AW translation rotation
run 2
25 65.4 f 21.3 65.2 & 21.3 -24.5 13.9 -24.7 13.8 82.5 f 14.4 82.4 & 14.3 45.0 f 10.8 44.9 f 10.8 32.8 8.4 32.9 8.4 12.5 f 5.0 12.4 4.9 -95.2 10.0 -95.2 10.0 -19.6 3.0 -19.7 3.0 3.7 3.1 3.7 3.1 3.1 f 3.1 3.7 f 3.1 300.5 f 52.3 300.3 f 52.2
* * * *
* * * * *
run 3
run 4
25 25 65.8 f 21.2 65.6 & 21.2 -24.2 f 13.8 -24.3 f 13.8 82.5 f 14.4 82.5 14.4 45.1 f 10.8 45.0 f 10.8 32.7 8.3 32.8 8.4 12.6 5.0 12.5 5.0 -95.0 10.0 -95.1 f 10.0 -19.5 2.8 -19.5 3.0 3.7 f 3.1 3.7 3.1 3.1 f 3.1 3.7 f 3.1 300.7 f 52.3 300.5 f 52.4
* * * *
* * * * *
temperature a Energies and rms fluctuations are in kJ/mol, temperatures in K, and length of runs in ns.
Values of the friction constants were obtained from ASA calculations performed on the C p conformation of the dipeptide. Following the approach of Pastor and Karplus,22the radius of a The system chosen for this study was the blocked alanine sphere with an area equal to the accessible surface area38 was dipeptide at infinite dilution in water at 300 K. The water model calculated and then used in Stokes' law with slip boundary used was the simple three-site TIPS model of Jorgensen? as adapted for analytic theories by Pettitt and R0ssky.3~ The conditions (see eq 6). The resulting atomic friction constants are shown in Table I, together with the relevant nonbonded paramXRISM pair correlations and properties for this potential model eters. These friction constants compare well with previously were the focus of a previous s t ~ d y . 3It~ has been demonstrated calculated values.21 that the XRISM pair correlations yield adequate agreement with For each simulation the initial conformation of the dipeptide simulationresults both for peak placement of the radial probability corresponded to the C709 conformation (4 = -67', $ = 68"), distribution function and for coordination numbers.10.32,33 which is the global energy minimum on the potential energy For the internal intramolecular potential, we use an empirical surface, U(1~).28 Each of the simulations differed only in the molecular mechanics potential with two-, three-, and four-body terms to mimic the gas-phase (vacuum)potential surface, U(l ~ ) . ~initial ~ random number seeds. Four parallel simulations were performed of 25 ns each. In addition, one of the runs was extended For convenience the nonbonded parameters are listed in Table for a total of 100 ns. The averages from each of the trajectories I. To this potential the liquid-state correction terms (AW(1,)) are presented in Table 11. Each nanosecond required 24 node are added to yield the total PMF, given by hours on an Intel iPSC/860. W1,) = U(1,) + A w l , ) (9) IV. Results and Discussion An all-atom model of the alanine dipeptide (22 atoms) was used througho~t.l3-3~*~6 The free energy surface W(),1 for the alanine Figure 2 shows the time histories of the 4 and $ dihedrals over dipeptide is displayed in Figure 1.28 the 100 ns of run 1. As one can see from the figure, rotation around $ is more facile than for 4. A fraction of the cogigurations The Langevin equation was integrated using a leapfrog generated from run 1 are also displayed as a I$,$map in Figure algorithm developed for stochastic dynamic~,3~ with a time step 3. It is apparent that the dipeptide sampled all of the minima of 0.5 fs. No bond or angleconstraints wereused. All nonbonded on the potential surface. Also shown in Figure 3 are the two interactions were calculated from the relevant atom-atom PMFs regions, a and 8, for which we were interested in calculating utilizing a grid look-up procedure. All nonbonded interactions transition rate constants. The a and 8 regions are so named were included. (No cutoffs were employed.) 111. Model
Smith et al.
6910 The Journal of Physical Chemistry, Vol. 97, No. 26, 1993 540
I
1
i
360
180
s
0.8
o -180 -360 IS1
540
I
1
0.2
' 0.
I
-540
20.
I
80.
60.
40.
I
I
V."
r ~
0.
100.
20.
40.
60.
80.
100.
Time (ns)
Time (ns)
Figure 2. Time histories of the Q and +dihedralsover the 100 ns of run 1.
Figwe 4. Average number of a states during each simulationasa function of time: run 1 (solid),run 2 (dashed),run 3 (dotted), run 4 (dot-dashed). The straight line is the Boltzmann-weighted average (0.56) from the solventmodified surface.
the transition a
-
j3 is a first-order process, then we have
Hence, the decay of Nuwill be exponential with a time constant of (kd kpJ-1. To determine the time constants, one can model theshort-time behavior of the transition counting autocorrelation function by
+
CN(~)= e x p ( - t ) where -120
-180 -180
-120
-60
0
60
120
180
CP
Figure 3. A $,+plot of the configurationsgenerated every 100 p during run 1. This is 1% of the configurationsused in the analysis. The a and B regions used in the rate analysis are also shown (seetext for the exact
definitions).
because they include the predominantly a-helix and @-sheet portions of the standard peptide 4,$ but they are extended beyond the usually defined regions. This was done to obtain a surface divided into two parts to simplify the transition analysis. The exact definition of the two regions used here are j3 if -240 < 4 < 0 and 0 < J/ < 240; otherwise a. The dividing lines between the two regions correspond roughly to maxima on the solventmodified potential surface.36 In this study we have concentrated on transitions of biochemical significance rather than kinetic interest, hence the seemingly arbitrary division between a and j3 states.
To calculate the rates for transitions between the a and j3 regions, we adopt the following definitions.12J9 Let N,=l if the dipeptide is in the a region of the 4,J/map; otherwiseNp0. The transition counting autocorrelation function is then defined as
where6Nu(r) =N,(t)- (Na),and ( N a )is thefractionofastates over the whole trajectory. With this technique recoil transitions are not problematic, unlike for simple counting approaches.12If
TN = l/(kap + k@a) (13) The conditions under which the time constant may be obtained from stochastic dynamics simulations have been investigated by ZhouM and are dependent on the barrier height and friction constant. For barriers of more than 3 or 4 ksT, which are representative of the present surface,36 the decay of Cdt) is expected to be exponential for all finite values of the friction constant.40 At equilibrium we also have the following relationship
from which we obtain
We have applied these techniques to the present simulation of the alanine dipeptide to calculate rates for transitions between the a and /3 regions starting from the global energy minimum which is a @ state as defined by Figure 3. In addition, we have also modeled the number correlation function assuming a biexponential decay of the form CJt) = A, exp(-k,t)
+ A, exp(-k,t)
(16) where kl and kz are rate constants. This allows for the possibility of a more complicated rate process than that of a purely first order decay from state a. In Figure 4 the average number of a states as a function of time is displayed. Also shown in the figure is the Boltzmannweighted population of a states calculated from the potential surface.36 This surface predicts the fraction of a states to be 0.56, and although this value refers to an adiabatic free energy surface, it should be a good estimate of the infinite time average
The Journal of Physical Chemistry, Vol. 97, No. 26, 1993 6911
Stochastic Dynamics of the Alanine Dipeptide
n ~
0.6
dynamics. The backward (kb,) rates are, with the exception of run 3, consistently slightly larger in magnitude (cf. k,in Table I1 of ref 12). The rates obtained assuming a biexponential decay span the range of values obtained here. Again, the rates obtained for run 3 are probably distorted by an inadequate sampling of the a state. The relaxation times observed for the present simulation are an order of magnitude larger than those obtained for the vacuum surface at 800 K. Also, in the latter simulation, the B state is more heavily favored with (Nb) 0.95. The present results for a solvent-modified potential surface at room temperature are likely to be more appropriate for use in protein folding studies. We have also analyzed the transition process by performing a hazard analysis of the dynamics traje~tory.~' In this procedure one calculates a set of n first passage times, which are defined as the time between arrival in the reactant well and a subsequent transition into the product well. The calculated first passage times are then arranged in ascending order and plotted against the cumulative hazard, which for point k is given by
-
\
0.0
0.
.-,500.
I
I
1000.
1500.
2000.
t (in ps) FigureS. Number correlationfunctionforescapefromtheastateobtained from the 100-ns simulation (solid). Also shown are the fitted singleexponential (dashed) and biexponential (dot-dashed) decays as given in
Table 111.
TABLE IIk Averages, Time Constants in ps, and Rates in ns-1 from the Stochastic Dynamics Runs run1 run2 run3 run4 0.62 0.72 0.37 0.54 Cfvu ) single exponential 123.6
114.4
77.9
kag ke,,
3.06 5.03
2.43 6.31
8.05 4.79
A1 Az
0.437 0.563
0.407 0.593
0.587 0.413
20.0
25.2 1.0
22.0 2.1
7N
biexponential ki
kz hazard kiH kzH
0.7
113.8 4.05
4.74 0.470 0.530 17.1 1.1
16.0 0.4
for our flexible model (potential of mean force surface). As one can see in Figure 4, the average number of a! states has not fully convergedin the simulations;even for the 100-11s run, the deviation from the calculated value is -10%. What is most striking is the slow rate of convergence to the expected value. This suggests that even for small systems, such as the alanine dipeptide, simulations on the order of microseconds are required in order to sample all the states with their correct Boltzmann weights. Larger systems would clearly exacerbate this problem. It is possible that an extended series of short runs (25 ns) would be the best way to obtain a converged value. In Figure 5 we display the short-time behavior of the number correlation functioncalculated from the 1OO-nsrun. It does indeed display a fast initial exponential decay, followed by a much slower exponential decay. The time constants and the corresponding rates calculated assuming both a single-exponential and a biexponential decay are presented in Table 111. The number correlation function obtained from run 3 is markedly different to that of the other three runs. This difference is related to the small value of (Nu)=0.37, indicating that during this run the dipeptide preferentially sampled the B region. This was also apparent from the time history of the 4 dihedral which displayed relatively few transitions compared to the other three runs (data not shown). The small number fluctuations of a! states enters into the numerator of the correlation function, causing it decay faster than the functionscharacteristic of runs with a more normal distribution (and consequently fluctuations between) of a! and @ states. The prime mechanistic reason for this is in the low number of 4 transitions for that run. The calculated rates are to be compared with those obtained in a recent high-temperature (800 K)stochastic dynamics study of the alanine dipeptide on a vacuum potential energy surface;12 the latter study had as its primary objective a search of configurational space rather than an examination of dipeptide
For a Poisson process the slope of this plot is equal to the correspondingrate constant. For cooperativeprocesses involving rotation around two or more dihedrals the rate constant obtained is a composite of the individual rates,42 while for noncooperative rotations a more complicated hazard plot is obtained. Due to the large friction constants of atoms in the present simulation (-60 ps-I), errors due to recoiling should be negligible.12 A hazard plot generated from the 100-ns run is shown in Figure 6a. It is interestingin that there appears to be at least two distinct processes occurring on this time scale, correspondingto two different rates, which we label klH and k2H. The two rates obtained from this plot are also presented in Table 111. Both of these rate constants agree well with the rates obtained from the correlation function approach assuming a biexponential decay of C d t ) . To further investigate the nature of the two rates, we have decomposed the hazard plot to determine the contribution from the 4 and $ dihedrals. The results are shown in Figure 6b. It is evident that the faster rate corresponded to transitions around the $ dihedral and the slower rate to transitions around the 4 dihedral. This is in accord with the time histories of 4 and $ plotted in Figure 2. It is also noticeable that these two rates did not appear to be significantly coupled; Le., rotation around 4 is essentially independent of rotation around $.
V. Conclusions We have implemented a method based on an implicit model for studying certain aspects of dynamic and equilibrium behavior in a condensed-phaseenvironment. Theeffect of solventcollisions are modeled by Langevin dynamics, and the vacuum nonbonded interactions (Lennard-Jones and Coulombic) are replaced by an effective potential of mean force for each unique atom pair, which includesthe effect of solvent correlations on the potential surface. This large reduction in the number of degrees of freedom allows the possibility of performing many nanosecond trajectories to study transitions and rate constants in sufficient detail to propose mechanisms based on the details of an atomic model. The approach is based on a desire to retain molecular details for selected atoms or degrees of freedom of the solute, while the solvent degrees of freedom, assumed to be of less direct interest, are projected (integrated) out of this particular picture of a manybody system. A number of ad hoc procedures (such as a distancedependent dielectric constants) for the approximate inclusion of solvent have been used relatively successfully for some solvation properties, but they do not contain the effects of the atomic nature of the solvent such as packing and explicit hydrogen bonding. To
6912 The Journal of Physical Chemistry, Vol. 97, No. 26, 1993 8.
a
8
8.
4 ;
