7330
J. Phys. Chem. 1995, 99, 7330-7338
Subspace Method for Long Time Scale Molecular Dynamics Attila Askar Mathematics Department, Koc University, Istinye-Istanbul 80860, Turkey
Brian Space and Herschel Rabitz* Department of Chemistry, Princeton University, Princeton, New Jersey 08544 Received: July 12, 1994; In Final Form: October 12, 1994@
This article presents the dual foundations of an approach to the long time integration of molecular dynamics equations. This work demonstrates that the long time scale molecular dynamics takes place in a relatively low dimensional subspace spanned by a set of collective eigenmodes and that the molecule remains in this subspace for long spans of time, at least on the order of picoseconds or greater. Both of these points are central to constructing an efficient numerical algorithm for long time scale dynamics. The first point allows for automatic filtering of the high frequency components of the dynamics and thus enables the stable use of very large time steps. The second observation allows for possibly using the same basis set for a long time span. Calculations are presented to substantiate these points. The model molecule consists of a 32 atom chain for which the equilibrium configuration is a helix. The chain interactions are of a two, three, and four body nature, allowing respectively for stretch, bending, and torsional motions. The molecule is intentionally chosen to undergo extreme dynamical changes for a severe test of both assertions 1 and 2 and the resulting algorithm. Moreover, one initial state is chosen very far from thermal equilibrium, with all of the energy of the molecule residing in a single local normal mode of the equilibrium configuration. Another initial condition, with a thermal distribution of kinetic energy, is explored. The methods presented are shown to be capable of handling both diverse initial conditions. The dynamical results permit a detailed analysis of the spectral aspects of the problem and provide further support for the subspace concept and the methodology based on it. Stable subspace dynamic integration for time steps of 6t = 100 fs were executed, and the results agree well with the full dynamics (for which the maximum allowable time step using standard molecular dynamics would be 6t M 1 fs); the dynamics included torsional barrier crossings.
1. Introduction There is presently great interest in developing long time molecular dynamics (MD) algorithms for use in polymer and biopolymer systems.’-ll The highest frequency motion of the fully coupled set of Newton’s equations determines the largest stable integration time step for the system. Finding an a priori method to systematically filter out the high frequency portion of the collective dynamics would permit the use of a very large time step for integration out to much longer times. Although the high frequency motions of a system may not be of much interest in the context of long time dynamics, it is not obvious how to directly access just the low frequency motions. Also, it is desirable to avoid eliminating degrees of freedom from a system by intuition, as is done in traditional constraint metho d ~ . ~ %For ’ * example, gas phase stretching or bending motions are often removed in constructing the dynamics model itself, with no allowance for possible couplings between removed and retained degrees of freedom which might be important in the condensed phase dynamics. The subspace molecular dynamics method (SMD) is formulated in a manner in which one can naturally identify and represent a continuously variably sized set of collective low frequency coordinates of a condensed phase system. For the purposes of this paper, the chosen coordinates are a basis of low frequency local normal modes of the coupled collective system. The dynamics is performed by updating the local normal coordinates each time after the numerical solutions to the projected equations of motion are advanced by a dt *Abstract published in Advance ACS Abstracts, May 1, 1995.
0022-3654/95/2099-7330$09.00/0
appropriate for the set of modes included in the subspace. These low frequency modes are shown to produce accurate long time dynamics. Further, the error in the dynamics caused by eliminating the high frequency motions of the system appear to be primarily local in time and do not seem to accumulate. Therefore, the approximate low frequency dynamics tracks the full molecular dynamics over times that are large compared to local dynamical deviations. These observations are consistent with recent works1OJ1which project full MD simulations onto low frequency collective normal coordinates and show that the “essential dynamics” is spanned by these low frequency Coordinates. This work emphasizes the feasibility of SMD for polymer dynamics. Computationally, it is not reasonable to calculate even a relatively small set of eigenvectors of the system force constant matrix at every time step.13 But, in SMD it is not necessary to use the eigenvalues and eigenvectors of the force constant matrix to analytically advance the equations of motion. Rather, a set of local eigenvectors can be used as a basis for nonlinear subspace dynamics. A computationally efficient version of the SMD method was previously successfully applied to highly anharmonic crystals and g1a~ses.l~ However, before tackling technical obstacles inherent in the SMD approach as applied to complex systems such as proteins, it is necessary to establish that the low frequency subspace dynamics is representative of the long time behavior in a polymer system. These systems exhibit complex and dynamically coupled torsional, bending, and stretching motions. Encouragingly, it will be shown below that subspace dynamics with only 24 of 96 degrees of freedom represented (which corresponds to a frequency cutoff 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 19, 1995 7331
Long Time Scale Molecular Dynamics of about 120 cm-’ in a model polymer with a dense spectrum and a maximum frequency of 1850 cm-’) provides an excellent and robust dynamical model of the system. This suggests that as long as one can span the low frequency dynamics of a complex condensed phase system, then one can integrate only the low frequency dynamics out to long times using the SMD formulation. It will also be suggested that if the computational obstacles associated with using a basis of low frequency modes in internal, as opposed to Cartesian, coordinates can be overcome,14J5then subspaces can be expected to be valid for times on the order of picoseconds. To perform SMD on large protein systems, even partial diagonalization of the force constant matrix, which still requires the computation of the entire force constant matrix (an operation which scales as the square of the number of degrees of freedom in the problem), will become undesirable. For systems with thousands of degrees of freedom it will become desirable to break the system up into subsystems on which we can perform SMD, while treating the subsystem-subsystem interactions with use of multiple time scale MD methods’ or by more approximate methods. The incorporation of substructuring ideas into SMD is currently under development. In section 2 below we summarize the subspace algorithm. The model polymer is described in section 3. The results obtained with this model are presented in section 4. Section 5 presents conclusions.
2. Subspace Molecular Dynamics Method The basic idea behind SMD is to judiciously employ the vibrational spectrum of the molecular system to identify the subspace of the collective vibrations where the low frequency dynamics “mostly” takes place. A subspace is selected at the beginning of the calculations from the eigenvectors of the locally linearized system: based on a high frequency cutoff of the spectrum. Using these eigenvectors, the full dynamical equations are projected into this subspace leading to a smaller number of effective dynamical equations. Finally, this reduced set of equations is integrated by a suitable scheme to evolve the molecular system through a large time step 6t. In this scheme, after progressing for some time ( i e . , possibly many time steps), the subspace may need updating to continue again for further time. An extreme is that of using the same subspace throughout the whole process. Indeed, such an approach has proven to be successful for studying the thermodynamic properties of anharmonic crystals and g1a~ses.I~In crystals and glasses, the solid nature of the system does not permit large deformations and the dynamics consists of oscillations about a mean configuration. In another extreme, that of the present results, the polymer chain undergoes extreme deformations and the time steps are taken to be very large with the subspace being updated frequently or even at every time step. The subspace is selected at each step by constructing again the relevant low frequency eigenvectors of the locally linearized system. In this latter case, the linearized equations are integrated numerically but are updated frequently. To elaborate further on the method, consider Newton’s equations for a molecule with n atoms.
(2.3) Here q is a 3n vector of mass weighted Cartesian positions with qi = mil!xi,where mi is the mass of the atom corresponding to
the ith Cartesian coordinate, xi. g is the vector of system forces in mass weighted coordinates, and gi = -aV(x)/aqi, where V(x) is the potential energy of the system. qo and 40 denote respectively the initial conditions for the position and velocity vectors. In order to obtain a reduced set of dynamical equations, consider the local force constant matrix at 90: (2.4) and the coordinate transformation:
q = qo
+ Pz
(2.5)
where P is a rectangular matrix with its columns made up of a subset of m of the 3n eigenvectors of F and z is the reduced vector of dimension m. It must be noted that the product P’P = I,, forms the m x m identity matrix of the subspace, while at the same time the product PPTis a projection matrix into the m dimensional subspace. Substituting the transformation of q in eq 2.5 into the original in eq 2.1 and premultiplication by PT yields:
z = PTg(q) = h(z) The key assumption introduced here is that Pz % PPT(q - 90) q - qo, implying that the essential dynamics of z is mostly contained in the retained low frequency subspace. Here the effective force h in the subspace is now a function of the projected coordinate z. It should be noted that while the equations in eqs 2.1 and 2.6 are in the same form, the former is in terms of the 3n dimensional vector q and the latter is in terms of the m dimensional vector z projected into the subspace, with m 1
(3.8)
This simple model potential does not include nonbonded interactions, but their absence will not qualitatively effect the spectrum of the polymer and is therefore not a key issue in the efficacy of SMD. In fact, SMD successfully described anharmonic glasses and crystals where the potential interactions were exclusively n~nbonded,’~ and there was no obvious frequency separation provided by the potential interactions. Further, nonbonded interactions are simple distance dependent interactions which change slowly with respect to a fixed reference frame since they do not depend on the orientation of atoms with respect to the fixed reference frame used in the present version of SMD. In the calculations, the mass m of all the polymer atoms is 2 x g, and the initial configuration, corresponding to mechanical equilibrium, is a helix defined as
rk = (a cos kd, a sin k h , hk)
k = 1,2, ..., n (3.9)
For a = h and A = d 3 , the equilibrium values for the bond length and the bending and torsion angles are found to be
Figure 1. Notation and definitions of the kinematic variables: bond
length, bending angle, and torsion angle. method. The kinematic variables are defined in terms of the 3n Cartesian coordinates {x} (each atom located at the position vector ri, i = 1, 2, ...,n), for the atoms forming the polymer as
Above, det denotes the determinant equivalent to the triple product of the vectors appearing as its arguments and Cjk is the unit vector along the bonds between the particles at rj and rk:
Figure 1 illustrates these internal variables. For a chain having n particles, the internal variables consist of n - 1 bond lengths, n - 2 bending angles, and n - 3 torsion angles. In the transformation from the 3n Cartesian coordinates of the n particles, the remaining six variables consist of the center of mass coordinate, X = &m,rfi,mj and the rigid rotation variable Y = &X x rj. It must be remembered that the interaction potential is cyclic in these latter variables as a consequence of its invariance under rigid body translation and rotation. The potential energy is the sum of the stretching, bending, and torsional contributions correspondingly as
bj” = h a
(3.10)
6; = 41.410a
(3.11)
4; = 44.415a
(3.12)
a = n/(180 rad) The numerical values of the parameters used in the calculations are a=h=lA
(3.13)
Ki = 25.00 (kcaYmol)/A’
(3.14)
KE,= 3.625 (kcal/mol)/rad*
(3.15)
KTi= 0.0550 kcaYmol
(3.16)
They are chosen to give a frequency spectrum that resembles a helix without explicit hydrogens, with a maximum frequency of about 1850 cm-’. The ratio of the potential parameters for stretches, bends, and torsions are also chosen to resemble a polymer system. One initial condition is chosen with the molecule at its mechanical equilibrium configuration subjected to an initial velocity in the shape of the 8th mode of the linearized system around equilibrium with an amplitude of 9127 K. Another initial condition has the molecule at mechanical equilibrium with a Maxwell-Boltzman velocity distribution corresponding to an initial temperature of 300 K. These choices of initial conditions are rather arbitrary, and the physical conclusions regarding the SMD method are not sensitive to the initial conditions. 4. Results and Discussion
with
Si = ( 1/2)Ki(bi- bp)’
(3.6)
Bi = (1/2)KEi(ei- 6:)’
(3.7)
The results of the calculations are shown in Figures 2-6. The figures have a dual purpose: (1) to provide a basic understanding of the large deformation .dynamics in terms of its modes and (2) to test the SMD method by comparing with the results from a full dynamics calculation. The SMD and the full dynamics are performed using the local eigenvector^,'^ although other methods could be used as well. The chain
J. Phys. Chem., Vol. 99, No. 19, 1995 7333
Long Time Scale Molecular Dynamics
z
30
30
25
25
20
20
(i) 15
z
(A)
f
t
0
20
15
10
10
5
5
0
0 40
Y 30 25
20
20 -
15
15.
10
10.
0I
I
0
20
40
Y
60
I
80
100
(4
i
I i
0I
3
L
(I
i
I
80
(4
i ii
-
25
5
t
60
A
/
0
20
40
Y
60
80
1
(4
Figure 2. Projections of the polymer configuration as calculated with (a, top, left) the full 96 modes; (b, top, right) 48 modes; (c, bottom, left) 24 modes, with dt = 1 fs; and (d, bottom, right) 24 modes, with 6t = 100 fs. Each snapshot is shifted by 0.5 ps corresponding to 17.5 8, along the
abscissa consists of n = 32 atoms with the full dynamics involving 96 modes (which is equivalent to standard MD), and the SMD is evaluated at subspace dimensions (modes) of 48 and 24. The present model chain is more flexible than a protein due principally to the absence of hydrogen bonds along the helix. A synopsis of the dynamical results is given below. For the first initial condition, the impulsive loading and the large initial energy of the system are chosen in such a way as to force the chain to extreme deformations far away from both mechanical and thermodyamical equilibrium. (i) Projected snapshots of the dynamics in time are shown in Figure 2. The time evolution of the various geometric and energetic measures of the molecule are shown in Figure 3a-e. The figures demonstrate the extreme deformations that the polymer configuration underwent. The overall dynamics is well described by the SMD, even out to the longest time shown. Figure 3e shows the time dependence of the system potential energies for the full space and subspace runs. We can see that both subspace potential energy surfaces follow the potential energy of the full dynamics but eliminate some of the high frequency “noise” associated with the constrained higher frequency modes. It is important that local deviations in the potential do not seem to accumulate, and the potential energies
track well throughout the run. Figure 3e also shows the total energy of the system for the various runs. The total energy is necessarily conserved for the full space calculation, while it is only conserved for subspace runs until an update is performed. There is a slight degradation in the energy for the subspace runs. This is due to imperfect overlap in subspaces during updating from time step to time step and is in some way a measure of the error made in the calculation of the low frequency dynamics.13 (ii) Parts a-c of Figure 4 display time averaged results respectively for the position of the particles, and bending and torsion angles along the chain. It is seen that the particles have undergone large excursions; the bending angles did not change much, although the dihedral angles did. In fact some of the dihedral angles have made transitions over their barriers. An equally important observation concerns the slow changes in going from one global polymer configuration to its neighbor. This is indicative of the dominance of global modes. Most importantly, the high frequency dynamics averages out and does not significantly influence the overall motion. This observation is key supporting evidence for the SMD method. (iii) The partitioning of the internal energy of the molecule (not shown here) between the various types of interactions shows
Askar et al.
7334 J. Phys. Chem., Vol. 99, No. 19, 1995 70
\
-
25
60
50
\
2ot
40
x (ded 30 20 10 1
0.5
0;
1
1.5
s
2
time (ps)
time (ps) 4500
1000,
G
1
900 800 -
700. 600 IJm
(Aa) 500 400 300 -
2oo 100
t/
0;
time (ps)
0.5
I 1
1.5
2
2.5
time (ps)
10000 9000 -
..........................................
8000 -
7000 -
6000 5000 -
4000 -
3000 -
~
0.5
*O
1
1.5
2
2.5
time (ps) Figure 3. Time evolution of the geometry of the chain as calculated with the full 96, 48, and 24 modes. In the remaining figures, the thin line represents the full 96 mode space, the dotted line the lowest frequency 48 mode space, and the thick line the lowest 24 mode space; (a, top, left) the distance between the f i s t and last particles; (b, top, right) the azimuthal angle, d e f i e d as the angle between the vector defined by the f i s t Z;)”’; (d, middle, right) .,Z In parts c and d, ,Z = &n,[(yL - ycm)’ (z, - am)’], Z, and last particles and the z-axis; (c, middle, left) ,Z = = C,m[(x,- xC,)* (z, - z,,)~] and Za = C,m,[(x,- x,)’ (y, - Y,,)~], are the moments of inertia of the body with respect to the initial center of the mass. Here, the subscript cm refers to center of mass, and x , y , z refer to a Cartesian component of the position vector. The moments of inertia reported are normalized by the mass as the masses are equal for the polymer atoms; (e, bottom, left) is the total and potential energy (kelvin) of the system as a function of time.
+
(ZL +
+
x,
+
Long Time Scale Molecular Dynamics
J. Phys. Chem., Vol. 99, No. 19, 1995 7335 2 0 0 0 1 ' I
I
"
I
I
I
,
,
1500 -
1000 E
(m-')
500 5-
0-50
5
10
15
20
25
30
-500
Particle #
Made Number
50 1
2o0
1
5
10
15
4 10 20 30 40 50 60 70 80 90
20
25
2.5
30
Bend # Mode Number
Figure 5. Time averaged values of the frequencies of the chain 5s calculated with the full 96, 48, and 24 modes: (a, top) mean value, E (b, bottom) standard deviation, ulE. For mode numbers where the average frequency is imaginary, it is shown as a negative frequency for convenience.
-100
5
10
15
20
25
I 30
Torsion # Figure 4. Time averaged values over the interval of 2.5 ps as calculated with the full 96, 48, and 24 modes: (a, top) @e particle Cartesian coordinates Z, J , Z; (b,-middle) bending angles, 8;(c, bottom) torsion (Le. dihedral) angles, 9.
that, as expected, most of the energy of the system is in the torsional interactions with a moderate contribution in bending and less in the stretches. Indeed, the latter two become significant only when the dihedral angles go through a barrier transition. This is especially interesting since, in traditional
constraint methodologies, the stretching, and possibly the bending, motions would be excluded from the dynamic^.^ Yet, these intemal coordinates are part of the low frequency motion during the important barrier crossing dynamics and are represented in the SMD method. (iv) The plots in Figure 5a,b show time averaged values and the standard deviations of the dynamical eigenvalues. Although the SMD calculations used a subset of these, the whole set was calculated even for the SMD calculations in order to monitor the unused ones and compare them with their counterparts in either the subspace or the full space calculations. The eigenvalues and eigenvectors are calculated from the force constant matrix evaluated at PPT(q- qo), and only for m = 96 is this exactly (q - qo). The standard deviations provide a measure of the changes that have taken place. The behavior in Figure 5a,b indicates that although the polymer has undergone very large deformations, the spectrum has barely changed. In particular, the higher eigenvalues have remained the same to about one part in a thousand. This invariance in the spectrum is a very important observation in relation to the viability of the SMD method. AI1 the essential aspects of the dynamics are captured by the subspace motions, including the regeneration of the high frequency modes at any time. The latter behavior
Askar et al.
7336 J. Phys. Chem., Vol. 99, No. 19, 1995
I
I
*
0
40
20
60
80
100
as it shows that one may return to a full picture of the dynamics at any future time after propagation in the subspace. (v) The range of the periods of the modes in the molecule varies from picoseconds to femtoseconds. The long time subspace is dominated by the low end of the spectrum, and it was possible to conduct the dynamical calculations with a time step of dt = 100 fs. The results shown in Figure 2a-c for this initial condition cover 2.5 ps of dynamics with time steps of dt = 1 fs, with the subspace updated at every time step. Other 48 and 24 mode subspace runs were performed using larger time steps. The largest stable time step was found to be roughly proportional to the ratio of the frequencies of the highest frequency mode of the system to the highest frequency in the included subspace. The results were found to be virtually insensitive to the time step. Figure 2d has the same initial condition as Figure 2a-c, but using an extreme time step of dt = 100 fs. This time step is about 6 times larger than the ratio of the maximum polymer frequency (about 1850 cm-') to the cutoff frequency (about 120 cm-'). This demonstrates that the subspace dynamics is stable and can give a reasonable description of the dynamics, even with extremely large time steps. (vi) The eigenvectors are difficult to display graphically. However, we can report some quantitative observations about their behavior. We monitored the product of QT(t).Q(t dt), where Q is the set of all 3n eigenvectors, for the full and subspace calculations. All 96 modes were calculated even for the smaller subspace calculations for the purpose of comparison. If the eigenvectors had not rotated, we should get precisely the identity matrix for the product above generated by the subspace, as each of the vectors would remain orthogonal. Consequently, the deviation of the off diagonal terms from zero in these overlap matrices is a direct measure of the rotation of the subspace from step to step as the system evolves in time. These products indicated a strong overlap for the low frequency modes along with almost no overlap for the high frequency modes. This observation is a further verification of the hypothesis of the invariance of the low frequency subspace even under extreme dynamical changes in the molecule. It should be noted that a rotation among the low frequency subspace vectors would still be fully consistent with all the underlying assumptions of SMD. (vii) Both of the observations on the spectrum and the associated eigenvectors relating to the high frequency modes indicate that these modes are almost linear oscillations. Furthermore, a simple theoretical consideration for harmonic oscillations indicates that the amplitudes for the rapid modes ought to be smaller than those of the slow modes, even under an equipartition of the energy. In fact, it is clear from the relative invariance of the higher eigenvalues that the subspace spanned by low frequency collective modes, in internal coordinates, would be nearly invariant. The rotation of the vectors seen in these calculations is caused by rotation of the vectors with respect to a fixed coordinate frame.13 Performing the subspace calculation in internal coordinates would eliminate this problem but would introduce computationally prohibitive matrix operations,14 although methods to possibly circumvent this problem are under deve10pment.l~ (viii) In order to understand the nature of the approximation with the subspace it is useful to compare it to restricted dynamical schemes which freeze some of the degrees of freedom and again obtain a lower dimensional system as compared with the original one. In this sense, the SMD method also may be qualified as a restricted dynamics scheme. However, the restriction here is done in a more subtle way. The physical degrees of freedom are not completely dropped in SMD. Indeed, the stretching and bending modes are active at critical
+
1
0
30 -
25 20 15 10 50 -
i! I
40
20
i
ii
60
80
i
i J
i
B
I
100
i
Long Time Scale Molecular Dynamics
J. Phys. Chem., Vol. 99, No. 19, I995 7337
instances of the evolution of the system, for example, while malung a transition in the dihedral angle from one well in the potential to the next. This feature of the subspace method is most attractive in that it keeps the contributions of all types of interactions at important phases of the dynamics, while still permitting large integration time steps. (ix) Further insight and improvement of the SMD method is suggested by the following analysis. Toward building a better understanding of the nature of the approximation, consider the exact formulation by carrying the complement, R, to the projection of q - qo defined in eq 2.5. Hence we define
q = qo
+ Pz,+ Rz,
(4.1)
Above, zs and z, represent respectively the slow and rapid parts of the dynamics and are correspondingly m and 3n - m dimensional vectors. Furthermore, the projection operators P and R have as their columns respectively the m and 3n - m eigenvectors corresponding to the slow and rapid modes. They have the properties
+
PTR= 0
(4.2)
RTP= 0
(4.3)
+
and (PT RT)(P R) is the unit matrix in the full 3n dimensional space. By definition z, and z, satisfy the relations
PTPz,= 2,
(4.4)
RTRzr= z,
(4.5)
The substitution of the transformation in eq 4.1 into the dynamical equation in eq 2.1, linearization of g about 9 0 , and the use of the orthogonality properties of the operations in eqs 4.2 and 4.3 yield
R ~ F P= J,,, R ~ F R= A,
(4.10)
The matrices A, and A, are diagonal, and the coupling matrices, J, are expected to be weak. The equations above have the form of two weakly coupled systems. In the slow mode, eq 4.6, the rapid modes act as an external agent, whereas in the rapid mode, eq 4.7, the slow components are the driving external force. For sufficiently large time steps the amplitudes of the rapid oscillations corresponding to z, are approximately given by sines and cosines, evaluated at sparsely separated values of the arguments, and hence will be almost random. This observation draws attention to the Langevin equation.16 Indeed, this latter equation and its generalizations may be derived by starting with the splitting of the dynamical system into slow and rapid modes as above. Hence, it is seen that the rapid modes act as a weak external random load on the slow dynamics. The SMD method neglects this weak effect. Regarding the rapid modes, their dynamics is controlled by two phenomena: rapid self-oscillations of small amplitude with an almost zero mean and a slowly driven contribution. An
important observation shown in Figure 5 is that the eigenvalues of the rapid mode, Le. Ar, are almost fixed along with the corresponding subspace. This suggests a three level hierarchy of approximations based on integration of the rapid equations in their linearized form along with an exact numerical integration in the subspace of the slow modes, as described in eq 2.6. The simplest scheme is that used in obtaining the results of this paper by dropping the rapid modes altogether. The next level approximation would drop the high frequency self-oscillations and keep only the particular solutions, i.e. the motion driven by the slow components in eq 4.7. The most sophisticated level would also include the self-oscillations of the rapid modes, either in exact terms or as a random term chosen from a distribution of white noise. In all the levels, the main approximation is the use of the invariance of both the spectrum and the subspace corresponding to the rapid modes. The amount of computer effort in the latter two improved levels is increased only by a one time calculation of the rapid eigenmodes and thus does not necessitate a reduction of the large time steps for the slow motions. (x) In all of the above results, it is seen that the subspace calculations produce reasonable approximations. This is all the more encouraging when considering the extreme deformations the polymer has undergone as well as the initial state of the system being extremely far away from dynamical equilibrium. For larger systems the subspace approximation is expected to perform even better. (xi) Parts a-c of Figure 6 show projected snapshots of 2.5 ps of dynamics which began with the atoms at their mechanical equilibrium and the velocities chosen from a thermal distribution, with an initial temperature of 300 K. Figure 6a shows a full space run, with dt = 1 fs and updating at every step. Figure 6b is the lowest 24 mode subspace with a dt = 1 fs and updating at every step. Figure 6c is the lowest 24 mode subspace with a dt = 16.7 fs. The full and subspace dynamics produced by the thermal initial condition agree well. Both subspace runs use a 24 mode basis, but different time steps. In general, SMD will be run from thermalized initial conditions, since initial conditions that localize the energy in a mode or set of modes would, at long times, lead to systems at different temperatures. This happens as the energy distributes itself among the available degrees of freedom, which is the number of modes included in the subspace. The dynamics is again well reproduced by a 24 mode subspace. Further, one can get a good representation of the system by taking full advantage of the larger time step permitted by the removal of the high frequency modes. 5. Conclusions
This paper presented the concepts and evidence that a well defined slow subspace exists for long time scale molecular dynamics. The existence of such a subspace is physically reasonable, and the essential point is that this subspace may be exploited to develop a very large time step integration method for molecular dynamics. The present paper does not present a full algorithm for such a SMD method but rather lays the foundation for the algorithm. Furthermore, a multilevel hierarchy of SMD algorithms was also suggested. The development of algorithms based on these concepts could significantly enhance our ability to access the physically important regime of long time scale dynamical processes in all aspects of material science.
Acknowledgment. We acknowledge support from the Office of Naval Research and Bristol Myers Squibb.
7338 J. Phys. Chem., Vol. 99, No. 19, 1995 References and Notes (1) Tuckerman, M.; Beme, B. J.; Martyna, G. J. J . Chem. Phys. 1992, 97, 1990. (2) Teeter, M. M.; Case, D. A. J . Phys. Chem. 1990, 94, 8091. (3) Elber, R.; Karplus, M. Science 1987, 235, 318.
(4) Brooks, C. L., III;Karplus, M.; Pettitt, B. Adv. Chem. Phys. 1988, 71. (5) Gibson, K.; Scheraga, H. J . Comput. Chem. 1990, 11, 468. (6) Beme, B.; Tuckeman, M.; Straub, J.; Bug, A. J . Chem. Phys. 1990, 93, 5084.
(7) Schlick, G.; Schlick, T. J . Comput. Chem. 1993, 14, 1212. (8) Tumer, J.; Chun, H.; Lupi, V.; Weiner, P.; Gallion, S.; Singh, C. Chem. Des. Automat. News 1992, 7, 16. (9) Watanabe, M.; Karplus, M. J. Chem. Phys. 1993, 99, 8063.
Askar et al. (10) Amadei, A.; Linssen, A.; Berendsen, H. Proteins: Struct., Funct., Genet. 1993, 17, 412. (11) Horiuchi, T.; Go, N. Proteins: Struct., Funct., Genet. 1991, I O , 106. (12) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1989. (13) Space, B.; Rabitz, H.; Askar, A. J . Chem. Phys. 1993, 99, 9071. (14) Wilson, E. B.; Decius, J.; Cross, P. Molecular Vibrations; Dover Publications, Inc.: New York, 1955. (15) Jain, A.; Vaidehi, N.; Rodriguez, G. J . Comput. Chem. 1993, 106, 258. (16) Askar, A.; Owens, R. G.; Rabitz, H. J. Chem. Phys. 1993,99,5316. E2417438
