12 Pages

#### Tuckerman-RESPA-JCP92

Course: CS 653, Fall 2008

School: USC

Word Count: 8767

Rating:

###### Document Preview

scale Reversible multiple time molecular dynamics M. Tuckermar?) G. J. Martyna and B. J. Berne Department of Chemistry, Columbia University, New York, New York 10027 Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323 (Received 13 January 1992; accepted 17 March 1992) The Trotter factorization of the Liouville propagator is used to generate new reversible molecular...

##### Unformatted Document Excerpt
Coursehero >> California >> USC >> CS 653

Course Hero has millions of student submitted documents similar to the one
below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

Course Hero has millions of student submitted documents similar to the one below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

scale Reversible multiple time molecular dynamics M. Tuckermar?) G. J. Martyna and B. J. Berne Department of Chemistry, Columbia University, New York, New York 10027 Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323 (Received 13 January 1992; accepted 17 March 1992) The Trotter factorization of the Liouville propagator is used to generate new reversible molecular dynamics integrators. This strategy is applied to derive reversible reference system propagator algorithms (RESPA) that greatly accelerate simulations of systems with a separation of time scales or with long range forces. The new algorithms have all of the advantages of previous RESPA integrators but are reversible, and more stable than those methods. These methods are applied to a set of paradigmatic systems and are shown to be superior to earlier methods. It is shown how the new RESPA methods are related to predictorcorrector integrators. Finally, we show how these methods can be used to accelerate the integration of the equations of motion of systems with Nose thermostats. 1. INTRODUCTION Recently we have devised reference system algorithms for treating problems with separation in time scales (stiff oscillators in soft fluids, disparate masses) and/or forces that can be separated into short and long range components3 ). These algorithms run at much faster speeds than do the standard algorithms; sometimes with speedup factors on the order of 20 to 30. Some of these reference system methods are not reversible in time, and thus may exhibit numerical instabilities after very long times. Moreover, they can not be used in conjunction with hybrid Monte Carlo methods because these require strictly reversible molecular dynamics in order to satisfy detailed balance.4 In this paper we show how to make these reference system algorithms reversible in time. The starting point is the Trotter expansion of the classical Liouville propagate? and the reversable Trotter expansion. From this we derive several new integrators for solving Newton equations of s motion. These new integration schemes involve the same number of force evaluations as the methods previously introduced, and therefore require no more CPU time in their execution. Furthermore, we find that making these methods reversible in time increases their stability over the RESPA and NAPA methods previously introduced. In Sec. II we present a necessary review of the Liouville operator formalism. Trotter factorizations of the classical propagator can be used to derive simple integrators. In Sec. II A the velocity Verlet integrator as well as a new integrator, the position Verlet integrator is derived. In Sec. II C generalized Trotter factorizations are introduced. In Sec. III we treat the long range force problem and apply it to a particular example. In Sec. IV we treat the multiple time scale problem and apply it to the disparate mass mixtures in Sec. IV A and to systems with stiff oscillators in soft baths in Sec. IV B. In Sec. V we show how to combine these methods to treat problems with separation of time scales and long and short range forces. In Sec. VI we relate the Trotter expansion In partial fulfillment of the Ph.D. in the Department of Physics, Colum) bia University. method to the more conventional predictor-corrector methods. These methods are then extended in Sec. VII to treat constant temperature molecular dynamics simulations based on a Nose heat bath. II. THE TROTTER EXPANSION OF THE LIOUVILLE PROPAGATOR The Liouville operator L for a system off degrees of freedom is defined in Cartesian coordinates as / iL={--. , ~L,F~L ,Hl= C (2.1) j=l axj Pj where IY = {xj,pj} are the position and conjugate momenta of the system, Fj is the force on thejth degree of freedom, and c *-* -**) is the Poisson bracket of the system. L is a linear Hermitian operator on the space of square integrable functions of F. The classical propagator is then 1 U(t) = e and the state of the system at time t is given by r(t) = wt)r(o), (2.2) (2.3) and U(t) is a unitary operator; i.e., U( - t) = U - f). In ( the following we will decompose the Liouville operator L into two parts such that iL=iL, +iL,. + L,)f = [ei(L,+L2)f/P]P rL,(Ar/2)eiL,AreiL,(Ar/2) (2.4) For this decomposition the Trotter theorem yields ei(L, = [e 1P (2.5) where At = t/P. From this we define the discrete time propagator as + at 3/P2), =e LL,(Ar/Z)ei~Ar~iL,(Af/2) (2.6) 0021-9606/92/l 51990-l 2$06.00 @ 1992 American Institute of Physics 1990 J. Chem. Phys. 97 (3), 1 August 1992 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp Tuckerman, Berne, and Martyna: Reversible multiple time scale MD 1991 Because the three factors in G( At) are separately unitary it is easy to show that G(t) is unitary, and therefore that G - (t) = G +(t) = G( - t) . This means that any integrator based on this Trotter factorization will be reversible. The formal solution for any decomposition of the Liouvillian can be generated as follows: Let us define G(At) = e (Aht/2)F(x)d/dpeAr*d/d~e(Ar/2)F(x)d/d~ (2.15) Operating with this factorization on {x( 0) ,p( 0) ] and using the explicit property of any operator of the form eca that aq ecd dY~q) =f(q + ~1, (2.16) r,[Acr(o)] and = U,(At)r(O) = u,(At)r(o) (2.7) where c is independent of q, it is a routine matter to derive the integrator (expressed in terms of positions and velocities), x(At)=x(O) r,[Acr(o)] + Ark(O) +~F'x(o)l (2.8) to be, respectively, the state at time At when the system is propagated by U, (At) or U, (At) starting from the state JY (0). Then by applying the operators in Eq. (2.6) serially to r(0) gives m , (2.17) i(At) =i(O) +e{F[x(O,] +F[x(At)]}. This is identical to the velocity Verlet algorithms.1o Because G( - t) = G - l(t), it follows that this integrator is exactly time reversible. This guarantees stability of this integrator and also allows it to be used in the hybrid Monte Carlo algo&hrn.4. 1 . Another way to write this integrator is r(At) = r(At) =rl u, (9) rz [Atg+$(o)]], [+(o)]]). =2(O) +~F[x(OH, (2.10) x(At)=x(O)+ (+-,(At,r, (2.11) Ad (2.18) To generate this solution on a computer: (a) Start with the initial state r(O) and generate the motion under the propagator eiLIAt. This gives the state i(At) =f ( -$- +&F[x(At)]. > r, [At/2;r(o)]. (b) Use the state just generated as the initial state and generate the motion using the propagator eiL2Ar.This gives thestate r,{At;r, [Ar/2;r(O)]). (c) Start with state generated in (b) as the initial condition and generate the motion using eiL1. This gives the final state T(At) specified in Eq. (2.11). The solution for a more complicated breakup proceeds analogously. We note in passing that the Trotter expansion carried out to higher orders will yield higher order integrators. For example, to obtain an integrator good to 0( At 4, one would start with the Trotter expansion ei(L, -t- L,)At z-e iL,(At/2)eiL,(At/2)e Xei~'A'/2'eiL,'A'"' , _ iC(At /24) Starting from the initial condition {x(O),p( 0) } one computes the velocity at the half-step, then the position at the full step, and finishes by calculating the velocity at the full step. This looks like leapfrog, but in practice it is different. In the leapfrog method one initiates the process differently, does not calculate the velocities at the full steps, and only in some cases calculates the last velocity at the full step. It seems that the most common integrator used in simulations, the Verlet algorithm, cannot be derived using this formalism. However, it can be shown that standard Verlet gives exactly the same trajectory as velocity Verlet. The standard Verlet equations of motion are x(At)=2x(O)-x( f(O) = x(At) -At)+ +bWl, (2.19) (2.12) -XC - At) 2At where -iC= [(iL, +2iL,),[iL,,iL,]]. (2.13) If the initial condition is assumed to be x(0)-x( -At)=i(O) - 2mFb(0)l, At2 (2.20) Here [*.*;** ] denotes the Lie bracket or commutator. The disadvantage of using such a higher order integrator is that the commutator will bring in derivatives of the force with respect to position which may be dithcult to calculate. A. Examples of simple reversible integrators Consider the propagator generated by the subdivision, iL, =G?-; a dX iL, = F(x) a. JP This leads to the propagator then it can be shown by induction that the solution obtained by iterating these equations to time t, {x(t) ,i (t - At)} will be equal to that obtained from the velocity Verlet algorithm with initial conditions {x (0) ,i (0)). Equivalently, the standard Verlet equations of motion can be obtained from those of velocity Verlet using time reversal symmetry and eliminating the velocities. However, Eq. (2.20) is much more useful as it can be used to convert all the algorithms derived in this paper to the standard Verlet form, as is demonstrated in the Appendix. A similar analysis is possible for leapfrog Verlet, but this is much less widely used and we omit it. J. Chem. Phys., Vol. 97, No. 3,1 August 1992 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp 1992 Tuckerman, Berne, and Martyna: Reversible multiple time -3.0 scale MD I 1 I I I Interestingly, we can generate another reversible algorithm by taking (2.21) iL, = F(x) a. ap Applying Eq. (2.11) with this subdivision to lY( 0) gives iL, =.k-; C3X a -3.5 -4.0 f(At) x(At) =2(O) =x(O) + ArF x(0) +$5(O) [ + -$[f(O) +k(At)]. I , (2.22) 2 M 2 -4.5 This is an entirely new reversible integrator which we name the position Verlet integrator. We will compare it to the velocity Verlet integrator later. It should be clear from the foregoing that both of these algorithms are of order At 3 and that the Trotter factorization is an excellent method for inventing new algorithms. Higher order factorizations should generate higher order reversible algorithms.6 B. Numerical comparison Verlet integrators of the velocity and position -5.0 -5.5 0 1 2 3 At(x103 4 ) 5 6 FIG. 1. Comparison of the energy conservation dependence on time step as given by Eq. (2.1) for the velocity Verlet (dashed line) and position Verlet (bold line) integrators for the LJ( 12-6) potential. To compare the velocity and position Verlet integrators of Eqs. (2.17) and (2.22), respectively, we choose a pure Lennard-Jones ( 12-6) fluidFith 864 particles. The reduced temperature of the system is T = 2.0 and the reduced density is po 3 = 1.0. The potential is cutoff at r, = 3~. The comparison is made by studying the energy conservation as a function of time step At for a run whose real time is NAt. This energy conservation is defined to be A& At) = -$ $, / E(kA;(; =( ) /. (2.23) e if3 _ - i [ii uk(At/z)] u,,(At) k=l n-l : ,& U,,-,(At/2) +~(t3/P2), where U, (h) = eiLkhand At = t/P. Let us define n-l l-J U,(At/2) U,(At) k=l n-l k& U,,-,(ht/2) . II 1 (2.25) A value of NAt = 1.0 is chosen and kept the same for all runs. Thus for the two integrators, we solve the classical equations of motion and measure the energy conservation as a gnction of time step. In Fig. 1, we plot the dependence of hE on At for the two integrators. The result is very interesting. The curve shows that the velocity Verlet integrator is better at small time steps, but that it becomes unstable more quickly than the position Verlet integrator. For energy conservation tolerances acceptable for typical molecular dynamics calculations, it seems that the velocity Verlet integrator is superior. However, for stability at large time steps, the position Verlet seems better and would thus be a suitable integrator for hybrid Monte Carlo calculations. 1 (2.26) It is a simple matter to show that G( At) gives rise to a reversible algorithm given that for each U, (At) Ut, (At) = U, (At) = U, ( - At) (2.27) - At) = 1, and (i.e., each L, isHermitian). Thus G(At)G( G( At) generates reversible dynamics. Ill. REVERSIBLE REFERENCE SYSTEM PROPAGATOR ALGORITHM (RESPA) FOR SYSTEMS WITH SHORT AND LONG RANGE FORCES We have previously introduced RESPA3, for accelerating the integration of the equations of motion for systems with forces, F(x), that can be decomposed into short, F, (x) , and long range forces Fl (x) according to F(x) =F;(x) +F,(x). (3.1) C. Reversible algorithms factorizations based on generalized Trotter It is well known that the Trotter factorization can be generalized to a sum of operators. If the Liouvillian takes the form iL = i k=l iL, (2.24) then the propagator can be expressed as The short range force then determines the time step St. In standard methods the full force must be recomputed at the end of each time step. In the RESPA method, however, the short range force is computed after each time step and the long range force is computed every n time steps. This reduces the number of forces calculated and thereby reduces the CPU time. The decomposition of the forces is accomplished by introducing a switching function S(x) that switches from J. Chem. Phys., Vol. 97, No. 3,1 August 1992 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp Tuckerman, Berne, and Martyna: Reversible multiple time scale MD 1993 1 to 0 at some interatomic distance chosen to minimize the CPU time. Then F(x) = S(xV Cc) + [ 1 -S(x) l.F(x) (3.2) x(At) k(At) =x, Ar;x(O),i(O) =.ks At;x(O),i(O) ++F,[x(At)] . (3.8) and we make the identification F, (xl = WV (x), F,(x) = [ 1 -S(x)]F(x). (3.3) In these equations x,{h~x(O),k(O) + (At/2m)F, [x(O)]} and .ic,CAt;x(O),i(O) + (At /2m)I;; [x(O)]} denote the solutions of the short range force subject to the initial conditions{x(O),f(O) + (At/2m)F~[x(O)]}.Thefullfactorization above dictates that x, and jc, are to be determined numerically using the velocity Verlet integrator [see Eq. (2.17) ] iteratively for n small time steps of size St subject to the initial conditions {x(O),k(O) + (At/2m)Fl[x(0)]}. The factorization in Eq. (3.7) is implemented as follows: (a) Starting with the initial state {x(O),p(O)} generate the motion using the propagator ecAt Fdp. This gives the 2 In our original formulation we took the reference system to be the dynamical system with only short range forces F,(x). This was solved numerically for n time steps St~ht /n. The equations of motion for the deviations of the reference system from the real system was then solved only once for the time step At. This algorithm is not reversible in time. The Trotter factorization allows us to make this strategy reversible. The basic idea is to define the short range force F, (x) as a reference system force and to rewrite the Liouville operator as iL =i-$+Fs(x) a+Fl(x) JP ?JP state{x(OLpW + (At/2)F,[x(O) I}. (b) Using the final state of step (a) as the initial state, generate the motion using the middle propagator in Eq. (3.7). This is equivalent to using the velocity Verlet integrator iteratively n times on the system with only short range forces Fs. (c) Starting with state generated in (b) as the initial state, generate the motion using e(Ar 2)Fdp. This gives the final state Eq. (2.11). In Appendix A, we show to convert this scheme to one which is based on standard Verlet. In Appendix B, we show a few lines of FORTRAN code to illustrate how the algorithm of Eq. ( 3.8 ) is implemented in a typical program. The same treatment can be applied to the factorization given by Eq. (3.8) (3.9) but we shall not do so here because it leads to more computational effort in calculating the potential energy at each step than does the previous factorization. In addition, the position Verlet integrator Eq. (2.22) could have been used in place of the velocity Verlet integrator. Thus there are eight possible factorizations, of which we only discuss one. A system with long and short range forces. To compare reversible RESPA with nonreversible RESPA, we look at a system studied in our previous work.3 The system is a LJ( 12-6) fluid with a cutoff at r, = 3a in,a box containing 864 particles at a reduced temperature T= 1.0 and a reduced density PCY = 0.8. For the nonreversible RESPA as 3 outlined in Ref. 3, we study only the reference systems given in Eq. (2.2) of Ref. 3. That is, we let the reference force be simply F, (x) . In all comparisons, we make use of switching function discussed in our previous publications,3s with an 2 inner cutoff of 1.90. For reversible RESPA, we use the integration scheme of Eq. (3.8). These two methods are compared with each other and with ordinary velocity Verlet. We study the dependence of energy conservation on At according to Eq. (2.1) using NAt = 1.0 for all runs. St is fixed at 1.0~ 10 - 3. In Fig. 2 we show the plot of A2 vs At for the RESPA three integration methods. We note that reversible Gs,s (At) = eiL~(At/2)eArF~/d~eiL*(Ar/2), =iL, +F,(x)a. aP The propagator eiLArcan now be factorized as, Gls, (A0 = e (At/Z)F,(x)d/dpeiL~te(A~/Z)F,(x)d/dp (3.5) The middle propagator generates motion using the short range forces. Using the Trotter factorization (6r/2)F,(x)d/dped,d/d~e(dt/2)F,oa/d~ eiL& = [e 1n , (3.6) where St = At/n and n is chosen to guarantee stable dynamics. This gives a good approximation to the middle propagator. The propagators on each end of Eq. (3.5) involve the slowly varying long range force and thus will be well approximated as shown with large time step At = nSt. These factorizations give the overall propagator G,,, (At) = e (Ar/Z)F,(x)d/dp [e (df/Z)F,(x)d/dp Xedr*d/dxe(dl/2)F.(x)d/dp 1 ne(Af/2)F,Wd/dp (3.7) Examination of Eq. (3.7) reveals that because the long range force propogators correct only the velocities, the long range force as well as the short range force from the previous step carry over as input to the next step. Thus one long range force evaluation and n short range force evaluations will be required in a large time step At. There is therefore no more work required for this reversible propogator than for the corresponding nonreversible RESPA method previously introduced.3 When applied to the initial state {x(O),p(O)}, this propagator generates the solution (expressed in terms of positions and velocities) Downloaded 05 Jan 2006 to 128.125.52.12. J. Chem. Phys., Vol. 97, No. AIP license1992 Redistribution subject to 3,i August or copyright, see http://jcp.aip.org/jcp/copyright.jsp 1994 -4.0 ,-_-------- -- -~ 1 I I Tuckerman, Beme, and Martyna: Reversible multiple time scale MD I I *clock), d.al -. -A rnenibla RESPA _a-*_-I _--- -4.5 .--- /- _/--- e- time. As we shall see below these methods are not analytically reversible. The Trotter factorization of the propagator allows us to make this strategy reversible and higher order in time. A. Disparate mass systems -5.0 2 xl 0 _/-.-- -5.5 & _H-- _/-- -6.0 r Y Consider the case where the degrees of freedom of the system can be subdivided into fast and slow degrees of freedom labeled x and y, respectively. An example is2 that of a binary mixture in which one of the components consists of heavy particles (they subsystem) and the other component consists of light particles (the x subsystem). Then we can decompose the Liouvillian as iL = iL, + iL,, (4.1) where iL, =i---+F&,y) -6.5 1 I 1 I I 1 2 3 4 At(x103 5 ) 6 7 FIG. 2. Comparison of the energy conservation dependence on time step as given by Eq. (2.1) for straightforward velocity Verlet (dashed line), nonreversible RESPA (double dashed line), and reversible RESPA of Eq. (3.7) (bold line) for the long range force breakup. a ax -, a ah (4.2) (4.3) iL, =ja+Fy(xy) d. aY ah The propagator can thus be factorized as Gxyx (At) = ei~~(A /2)ei~P eiL~(A?/2~. and nonreversible RESPA agree well at small time steps, but that reversible RESPA is more stable at large time steps. This is to be contrasted with velocity Verlet for which the energy conservation degrades rapidly at large time steps. To quantify this, iz we wanted to conserve energy to within a tolerance of AE = 5 X 10 - 6, we see from the figure that a time step of 2 X 10 - 3 would be required for velocity Verlet, while nonreversible RESPA would use a time step of 6 X 10 - with n = 6, and for reversible RESPA, a time step larger than 8 X 10 - 3 with n around 8 could be used. From our previous work on long range forces,3 we predict that the CPU savings would be a factor of roughly 3 for nonreversible RESPA and around 4 for reversible RESPA. IV. REVERSIBLE REFERENCE SYSTEM PROPAGATOR ALGORITHM FOR SYSTEMS WITH MULTIPLE TIME SCALES We have previously introduced methods for accelerating the integration of the equations of motion for systems with two or more time scales. Examples are systems con, sisting of high frequency oscillators coupled to low frequency oscillators or slow baths, and systems consisting of massive (slow moving) particles interacting with light (fast moving) particles, the so-called disparate mass problem.2 These methods are based on defining a reference system for the fast subsystem, solving this either analytically or numerically for a sequence of small time steps, deriving equations of motions for the deviations from the exact trajectory, and then solving the equations of motion for this deviation and for the slow degrees of freedom numerically using a much larger time step. When the reference system is solved numerically the method is called RESPA (reference system propagator algorithm), whereas when it is solved analytically the method is called NAPA (numerical analytical propagator algorithm). These methods reduce the number of forces that must be calculated and thus reduce the CPU (4.4) For this GX,,Xfactorization the end propagators involve the fast motion whereas the middle propagator involves the slow motion. Although is it not as obvious as in the long range force case [ Eq. (3.7) 1, these factorization involves the same number of force evaluations as does the nonreversible disparate mass RESPA method previously introduced.2 This follows from the fact that FXY = - Fyx, so that this part of the force need only be calculated while advancing the fast coordinate. The fast propagator can be further factorized e iL,(At/Z) = (s1/2)F~/ap,~s~~aaa~~(6f/2)F~/ap, nn [e 19 (4.5) where St = At/n. The middle propagator in GX,,X also be can factorized, but because it involves only the slow motion we can use ei=Ar = e(ht/2,F~/aP~eArLa/ayecar/2)F~/ap, (4.6) If these factorizations are used in Eq. (4.4) the fast degrees of freedom and their conjugate momenta {x,p,) are to be determined numerically using the velocity Verlet integrator [see Eq. (2.17) ] iteratively for n small time steps St subject to whatever the initial conditions might be, whereas the slow degrees of freedom are to be determined using the velocity Verlet integrator for only one large step At. We can summarize the numerical procedure that follows from application of GXYX follows: as (a) Use the velocity Verlet integrator for n/2 time steps St = At /n to generate the state at time At /2 under the action of the the propagator eiLxAr starting from the initial state {x(0)s(O),pX (O),pu(O)}. Note that during this step &p,} is unchanged. (b) Use the velocity Verlet integrator for one time step At to advance the system from the final state of step (a) using the slow propagator Eq. (4.6). (c) Use the velocity Verlet integrator for n/2 time steps J. Chem. Phys., Vol. 97, No. 3,l August 1992 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp Tuckerman, Berne, and Martyna: Reversible multiple time scale MD 1995 I I I 6t = At/n to generate the state at time At /2 under the action of the propagator eiLPrn , starting from the final state in (b). Note that during this step fj,pY} is again held fixed. The resulting state is F( At). In the foregoing the forces on the slow coordinates are recalculated only once every time step, At = n&. If the dimensionality of the fast subsystem is small compared to the dimensionality of the whole system, the forces on the slow coordinates will thus be calculated much less frequently than would be required in the standard methods. There are many variations on the above theme. For example, we could have adopted the Trotter factorization (4.7) and then applied the velocity Verlet breakups. Alternatively, we could use the position Verlet integrator instead of the velocity Verlet factorizations of the different factors in GXYX or Gy.ry.There are eight possibilities in all. For the purposes of the subsequent discussion, we look only at GX,,Xwith the velocity Verlet [see Eq. (2.17) ] breakup applied to each of the factors in GX,,X. Numerical comparison for the disparate mass system. The comparison of reversible RESPA to nonreversible RESPA proceeds now on a system studied in our previous work.2 Again, we look at a sysf\em of LJ( 12-6) particles this time at reduced temperature T = 0.67 and reduced density pa 3 = 0.86, which correspond to the triple point. Again, the potential is cutoff at r, = 3~. The system is comprised of a mixture of 40 light (m = 1) particles and 824 heavy (M = 100) particles. The nonreversible RESPA simulations are carried out as described in Ref. 2 in which the reference force is taken to be Fr (x) = FX [x,y( 0) 1, whereas the reversible RESPA simulations are done using the propagator of Eq. (4.4). When this propagator is used the explicit algorithm for advancing the positions and velocities is given by At ;xo A0 PYO 3 > At (2 i( At) = i:, 30 20 aY0 9 > (2 x(At) =x, y(At) i(At) where x0 =x, 5zo= ir -$;xmf(oLY(o)] , ) (4.9) =v,[k;yUWUko], =Y,[WUW(OL~~], (4.8) Gyxy =ei~,(Ar/2)eiLpfeiL,(Az/Z) -3.5 ,-------------_--_ I cloclt, RmpA ar,e, ~ m.eralble FSSPA -4.0 / ,/ ,.I ,/ ,/ ,I -4.5 5 M 2 ,**- / ,/ -5.0 // ,- /- / / /- 0. /. / -5.5 --fi-. 0 0.5 A 1.0 2.0 1.5 At(xl0 ) 2.5 3.0 FIG. 3. Comparison of the energy conservation dependence on time step as given by Eq. (2.1) for straightforward velocity Verlet (double dashed line), nonreversible RESPA (dashed line), and the disparate mass breakup of Eq. (4.4) (bold line). versible RESPA is consistently better than velocity Verlet, reversible RESPA is dramatically better and remains remarkably stable even at large time steps. While velocity Verlet becomes unstable at At = 0.04, the energy conservation for reversible RESPA is better than 10 - 5, and for nonreversible RESPA it is better than 5 x 10 - 5. In Fig. 4, we plot the ,^ 2 0 x 3 a P 0 a I I 0 2 I 4 I 6 I 6 10 I 12 -2 11 0 , 2 I 4 I 6 I 6 I 10 I 12 t [ [ -$x(0),*(0)9(0)] x, and y, refer to the evolution ofx and y under the actions of ei=*r and ei=Y(,respectively. The velocity Verlet breakup of Eq. (2.17) is applied to the propagators eiLXrand eiL9. In all simulations, the light (x) particles are integrated using 6t = At / 10. The eneri-y conservation is measured as before using Eq. (2.1) an> T is chosen equal to 1.0. In Fig. 3, we show the plot of AE vs At. Figure 3 shows that while nonre- FIG. 4. (a) The instantaneous fluctuations of the energy about E(0) for velocity Verlet in the disparate mass system at a time step At = 2~ lo- . (b) The same except for nonreversible RESPA. (c) The same except for reversible RESPA. J. Chem. Phys., Vol. 97, No. 3, 1 August 1992 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp IQ96 Tuckerman, Berne, and Martyna: Reversible multiple time scale MD instantaneous fluctuations of the energy about the initial value. That is, we plot AE(Q = E(t) -E(O) (4.10) E(O) for each of the three integrators for a time step of At = 2 X 10 - *. In Fig. 3, this corresponds to an average energy conservation of about 3 x 10 - for velocity Verlet, 1 X 10 - 5 for nonreversible RESPA, and 3 x 10 - 6 for reversible RESPA. From these plots, we can determine the maximum of the fluctuations, and discern whether or not the system has a systematic drift. Although the behavior of the instantaneous fluctuations is different in each of the three cases,we see that they are all stable at this time step. B. A stiff oscillator dissolved in a soft liquid G( At) = e (Ar/Z) F+sf Xe iLpte (u/2) aey e (A~/~)Ax.Y) pajay a/ap, e (Ac/z))a,ay flwa/ap, (4.17) e(Ar/2) Xe (at/z) qda/ap, Using this factorization of the propagator, the integration procedure generated in terms of positions and velocities is At;x(O),k(O) + $0, , I At;x(O),f(o) +$(o) v(At) I +$(At) , Another kind of system of interest is that of a stiff oscillator with vibrational coordinate labeled x buried in a soft bath with coordinates labeled y. For convenience we will also consider the center of mass coordinate of the oscillator as well as its orientational coordinates as part of the bath. It is possible to break the force F, (xg) on the vibrational degree of freedom into a strong, possibly analytically solvable force, and a weak correction. The basic idea is to define a reference system force F,(x) and to rewrite the Liouville operator as iL =k-$-+F,,(x) d aPX 1 a+Ld+F,(x,y) JP, av -% =y(O) At2 + A@(O)+2m F y (0), L(At) =j(O) +& [Fy(0) + F,(At)]. (4.18) +m,Y JPY (4.11) where&y) = F, (X,-Y) - F.,(x). f contains forces on x due to the solvent atoms. This allows us to subdivide the Liouvillian into parts iL, =is+F,(x) and -!JPX (4.12) This particular factorization has the obvious simplicity that the solvent is integrated using velocity Verlet. This breakup is a generalization of the Trotter form as discussed in Sec. II C. The same treatment can be applied to the factorization given by Eq. (4.15) but we shall not do so here. Numerical comparisons for a sti#oscillator in a soft& id. The comparison between reversible and nonreversible RESPA proceeds, as before, on a system considered in our previous work. We look at a system of 62 LJ( 12-6) particles in which is imbedded a high frequency harmonic oscillat,r. The reduced temperature and density of the liquid are T = 2.5 and pa 3 = 1.05. The frequency of the oscillator is taken to be w = 300. Since the peak of the spectral density of the neat liquid at this density and temperature is around w = 20, the oscillator frequency w = 300 is an extremely high frequency. Since the reference force F, (x) = - mw22 is harmonic, we can evaluate the action of the reference system propagator analytically in the usual way as e iL,At iL, =B$j+FyW) ~+~(x,JJ) JPY Gy, tAtj = eiL~(At/2)eiLPt@iL~(d R x = x cos oAt + * sin @At, x = 2 cos wht - wx sin wht. . -% ap, (4.13) e iL& (4.19) The propagator eiLA*can now be factorized as (4.14) or G,(At) = e~,,cAt/2)e e LP LAt/Z). (4.15) Factorization of the middle propagator in Eq. (4.14) as iL$r = [e (b~/z)F~/ap~~g~jd,a~~(~f~)F?/ap, e (4.16) I9 where 6t = At/n. If the reference system can be solved analytically (e.g., a harmonic, cubic, or Morse oscillator) then the above factorization is not necessary. A simple way to include the solvent effects into the full propagator yields the following result for GY,.,, (At) : We follow the evolution of the system for 500 periods of the oscillator using ordinary velocity Verlet, NAPA, and reversible NAPA. The result of the comparison of the three integrators is shown in Fig. 5 in which is plotted the energy conservation vs At as in Eq. (2.1) . We see that the reversible NAPA as given by Eq. (4.18) is highly accurate at all time steps, while the nonreversible NAPA tends to become significantly less accurate than reversible NAPA for time steps greater than about 10 - 3. This is due to the fact that nonreversible NAPA, being less stable than reversible NAPA. However, since both NAPA methods are based on analytically solvable stiff reference systems we expect them to agree well at smaller time steps, which they do. The extreme stiffness of the reference system illustrates the advantage of using a reference system algorithm over ordinary velocity Verlet J. Chem. Phys., Vol. 97, No. 3.1 August IQ92 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp Tuckerman, Berne, and Martyna: Reversible multiple time scale MD 1997 _.. / C The middle propagator in Eq. (5.2) involves the fast degrees of freedom and the short range forces acting on them and can thus be further factorized as in Eq. (4.5), whereas the middle propagator in Eq. (5.3) involves the slow degrees of freedom propagating under the short range forces and can thus be factorized as in Eq. (4.6). The overall propagator then generates a reversible double RESPA algorithm which can lead to enormous savings in CPU times. The extension of this idea to the stiff oscillator problem is straightforward. VI. REVERSIBLE ALGORITHMS FROM PREDICTORCORRECTORINTEGRATORS The use of the Liouville operator formalism to derive integrators is a somewhat novel approach. It is, however, satisfying to know that even the standard integrators of molecular dynamics can be obtained from a first principles approach. Nevertheless, it is useful to draw the connection between the Liouville operator formalism and the more conventional predictor-corrector scheme.8S In this section we 3 address this issue. The usual approach in a predictor-corrector scheme is to predict the state of the system based on the initial conditions for the current time step, and then to correct the prediction based on the forces evaluated at the predicted positions. We show that the Trotter breakup ofthe propagator is equivalent to a predictor-corrector scheme. Consider the Trotter breakup of Eq. (2.15) which yields the velocity Verlet (6.1) Mathematically, this breakup is equivalent to the following set of equations $=iF(x), dx dt = v, (6.2) eiLAf = e(At/2)F(x)a/apeArjia/ax~(A2/2)F(x)a/ap -5.0 0 I I I I I I 2 4 6 At(xlO+ 8 ) 10 12 14 FIG. 5. Comparison of the energy conservation dependence on time step as given by Eq. (2.1) for straightforward velocity Verlet (double dashed line), nonreversible NAPA (dashed line), and the reversible NAPA scheme of Eq. (4.17) (bold line) for the stiff oscillator dissolved in a soft solvent. which is significantly less accurate than NAPA or reversible NAPA for reasonable energy conservation tolerances (10-3-10-40rbetter). V. SYSTEMS WITH SEPARATION OF TIME SCALES AND SHORT AND LONG RANGE FORCES In the foregoing we have treated systems in which forces can be subdivided into short and long range components. Reference system methods lead to a considerable saving in CPU time. We have also discussed systems in which there is a large separation in time scales such as stiff oscillators embedded in soft fluids or disparate mass systems. In systems with separation of time scales and with force fields that can be decomposed into short and long range forces it is possible to devise reference system methods that can accelerate the solution of the equations of motion. In a previous publication we were able to apply a double reference system approach (double RESPA) to greatly accelerate such systems. This method was not reversible. Here we show how one can generate reversible multiple RESPA procedures to handle systems like these. The propagator for the disparate mass system given in Eq. (4.4) can be expressed as (6.3)$= +F(x), (6.4) G,,, (At) = GI; where Gi$(At) = e (Af/2)F,(x.~)a/aP~~iL~f Xe and G$ (At) = e (At/2)F~(~y)a/aP~~iL,~t Xe (At/2vy,w)a/apy (Ar/2)F,,kYv)a/ap, (5.1) (5.2) with the prescription that Eq. (6.2) is solved for a half time step At /2, with initial conditions {x(O), v(O)} and holding This fixed. the yields predicted value x(O) v(At/2) = v(0) + (At/2m)F(x).ThenEq. (6.3) issolved for a full time step using the predicted value of v and x( 0) as the initial conditions. This yields the true position x(At). Finally, we correct the velocities using Eq. (6.4) with initial conditions x( At) and v( At /2). Following this procedure will yield the scheme of Eq. (2.18). Reversibility comes from having a symmetric set of equations. The middle step in which x is integrated can be thought of as a reference system integration done analytically, since eA* ax would be the a propagator for a reference system with F, (x) = 0. With this in mind, we can easily extend this idea to reversible RESPA. Now consider the example of a system with long and short range force components. The Trotter expansion of the propagator takes the form given in Eq. (3.7) @At = e(ACn)Fl(X)a/aP=iLpre(ht/2)Fl(x,a/aP , (5.3) J. Chem. Phys., Vol. 97, No. 3,i August 1992 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp 1998 Tuckerman, Berne, and Martyna: Reversible multiple time scale MD where iL, = ka /ax + Fs (x>a/ap. this into the equivalent form $=+F,(x), dX dt=vs As above we can cast (6.5) Here we have chosen to use the variable ~7= 1 dt instead s of c as introduced by Hoover in order to make the equation for 7 second order. In order to derive a reversible integrator for the NoseHoover equations, we recognize that these equations can be derived using the Liouville operator iL=zi$+F,.(x)a+f(x)d ap -+&+&+F,(p) ap a7 ap d, a& dv -=LF,(x), (6.6) dt m dv -=-l-F,(x). (6.7) dt m The scheme is now transparent. Eq. (6.5) is integrated for a half step At /2 with initial conditions x (0), v( 0). This yields the predicted value v(At/2) = v(0) + (At/2m)F,(O). Then Eqs. (6.6) are integrated for a full step At using the predicted value v(At /2) and x(O) as initial conditions. This yields a reference system trajectory in which the position is the true position x, (At) = x(At). Then the reference velocity is corrected by Eq. (6.7) using as initial conditions x( At) and v( At /2). As before, reversibility comes from having a symmetric set of equations. It should now be clear how to generate predictor-corrector schemes for other systems with reference forces. VII. A REVERSIBLE ALGORITHM FOR NOSk DYNAMICS We now discuss integration of systems with multiple time scales obeying Nose dynamics. 5 Recall that Nose 4P s equations of motion for a system of N particles in virtual variables can be generated from the Hamiltonian where F,,(p) = Z(p2/m) - 3NkT or F,,(i) = Z(m5z2) - 3NkT. We have introduced the reference system force F,. (x) and the deviation of the true force from the reference force f(x) = F(x) - F,(x). It should be noted that although the conserved energy is H =x I!-+ v(x) +&+ (3N+ l)kTvt (7.5) as can be verified by computing iLH using Eqs. (7.4) and (7.5), this energy is not a Hamiltonian, and the Liouville operator of Eq. (7.4) cannot be derived from it. However, it can be used in conjunction with Eq. (2.1) to measure how well the integration algorithm conserves the energy. Definpropagator as ing reference system the iL, = S/ax + F,(x)a/ap, the propagator eiU* is factorized as G(At) = e (Af/2) F,wm, e (At/4)flx) ia/* (At/4)f(x)a/ap,-~t/z) a/ap e(Af/2) F,&P)a/h, a/ap e (~t/2) *pa/ap H=E$+--H- pf + V(x) +%+ (3Nf l)krlns. (7.1) X,cAt/4~fcx~a/apecAr/2~ Xec~t/2~f7a/*e efLPr *palap The equations of motion derived from this Hamiltonian are 1~ ap -X7 p= -==aH Ps e l)kT] . (7.2) av .--$-... Xe (~t/4)fl~) (7.6) state (7.7) W-e act with [SJJj,p(Oj<bj,p; e---. cAt/2ma/ap~(p) this propagator on the initial (O)], and use the-fact that = 4cpecAt/z)i). . aH s=-.--=-, ap, _. The above identity as well as Eq. (2.16) are special cases of the general identity for the function eca a*c@f(q), ecammf(q) = ,caa3mfcg- 1lgcq) 1) By introducing the time scaling dt- dt /s and the change of variables p-p/s, p, -pJs, and introducing as a dynamical variable 77= In s, one obtains Nose-Hoover equations16 jc=--, P m =fCg- lk(q) -I- cl II, (7.8) where c is independent of q. The result is a reversible integrator for Nose dynamics, x(At) ItAt) =x,[Aht;x(O),f,], = x,[At;x(O),x:,] I Xe - (Ar/2)+70 ;; + $f(At) f(At), I p= -av-pq ax ePy j=p, Q 3NkT. (7.3) q(At) = ~(0) + At*(O) +$$Fq jj,, =C$- [f(O) I, J. Chem. Phys., Vol. 97, No. 3.1 August 1992 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp Tuckerman, Berne, and Martyna: Reversible multiple time scale MD 1999 MAO = j1(0) + $-$ where [ii-(O) I + Fv [i(At (7.9) jr, = f(0) +LLf(o) [ rio =tw +$F,MoH. 1 e- cat + 2)*Q$p done by iterating the velocity integration step starting the &prox ( At) = i ( 0 ) with approximation + At/m [F(O) - e(O)k(O) 1. Ten iterations ofthis approximation are used to obtain the final velocity 5 (At). Figure 6 shows the comparison of energy conservation defined by Eq. (2.1). We see that the reversible RESPA method for Nose dynamics is a dramatic improvement over straight velocity Verlet and is remarkably stable at large time steps. VIII. CONCLUSION Dynamical systems with multiple time scales pose a major problem in simulations because the small time steps required for stable integration of the fast motions lead to large numbers of time steps required for the observation of the slow degrees of freedom and thus to the need to compute a large number of forces. As we have shown in several previous publications, reference system methods (RESPA) greatly reduce the CPU time for the simulation of such sysIn tems. 1-3*12,18 this paper we have used reversible forms of RESPA that are superior to the previous RESPA methods both with respect to the order and stability of the integrators. Moreover, because these are reversible they can also be used in the hybrid Monte Carlo methods.4,5 We have limited this discussion to a few obvious factorizations of the propagator based on the Trotter factorization, and we have shown how one can use a more general breakup. This procedure is put in the context of predictor-corrector methods and is applied to several different dynamical systems. It is shown how simple factorizations give rise to the velocity Verlet integrator and to a new integrator which we call the position Verlet integrator. It is shown how the problem with multiple time steps like the disparate mass problem and the problem of stiff oscillators in soft baths can be treated. A very important class of problems in which the forces can be subdivided into short and long range forces is treated by these methods. It is shown how to combine these latter two classes of problems with very large speedups. Finally, we incorporate these RESPA methods in the treatment of systems interacting with a Nose heat bath for constant temperature molecular dynamics.The approach outlined here is easily generalizable to many more integration strategies. ACKNOWLEDGMENTS This work was supported by a grant from the National Science Foundation (NSF CHE-87-00522 and by a grant from the National Institutes of Health (GM43340-OlAl). One of us (B.J.B.) would like to thank Dr. J. Sexton and Dr. D. Weingarten for discussions about the time reversibility of the Trotter factorization. APPENDIX A In this Appendix, we show how to obtain reference system integration schemes based on the standard Verlet algorithm. The strategy is the same as that employed in Eq. (2.20). Thus the starting point for these derivations is the velocity Verlet based integrators. Consider, first, the case of the long range force integrator discussed in Sec. III. The velocity Verlet based algorithm gave the positions and velocities at each time step to be The simplicity of this integrator is that the bath is integrated with ordinary velocity Verlet, and the dissipative effect of the bath on the system is handled analytically in the e - (Ar terms. We see that this factorization is of a gener2)il, alized Trotter form as discussed in Sec. II C. In the Appendix, we show how to convert this algorithm to one based on the standard Verlet. Note that it is possible to incorporate the both variables 7 and P, into the reference system and define the Liouville operator as + iLr,n, JP where the definition of iLr,n can be inferred from Eq. (7.4). The breakup of the propagator is then G(At) = e LAf/2)Rx)d/dpeiLr.nAI~(Af/2)flx)d/dp (7.12) Although this tends to be a slower algorithm, it proves to be more stable when multiple Nose baths are involved and can handle large fluctuations in 17more stably. Numerical comparisons for No& dynamics. The integrator in Eq. (7.9) is compared to straight velocity Verlet on the system considered in Sec. III. That is, we consider a system of 864 LJ( 12-6) particles at a temperature 1.0 and density 0.8 using the long range force breakup to define the reference system. The mass of the Nose particle is Q = 2.0 chosen according to the formula Q = 3NkTr 2, where 7 is some characteristic time in the system. Since the forces are velocity dependent, the straight velocity Verlet simulations must be iL =f(x) - a -3.0 I I 1 I I ---_-- vdodt)r erlet rarenib,* aE.sP* I .-----3.5 /.*.- 1 I .-- 3 -5 - -4.0 -4.5 -_-_e--/ *__--- /, ,,= _.-I - _,,/- -e-_o-_e-- -5.0 -5.5 ,/------q 1 -6.0 1 1 2 1 3 I 4 At(xl0 I 5 ) I 6 I 7 FIG. 6. Comparison of the energy conservation dependence on time step as given by Eq. (2.1) for straightforward velocity Verlet (dashed line), and the reversible RESPA scheme for Nos.6 dynamics of Eq. (7.6). Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp J. Chem. Phys., Vol. 97, No. 3, 1 August 1992 2000 Tuckerman, Berne, and Martyna: Reversible multiple time scale MD x(At) =x, At;x(O),W) + $Jx(O)l] x(At) =xX At;x(O),f(O) 1 +&F,[x(O)]] , 2At FinaIly, to evolve the system coordinate x, we define x( -At) =e (AU~)+~C,VJ)~,~ -At> i(O) = v(At) - vI( - At) P-12) ++x(At)l. (Al) +[I--e (At/2)ioVJ) 1 In standard Verlet, we start with the initial conditions {x(0),x( -St)), where St=At/n. Using the condition in Eq. (2.20), the standard V er 1et SC eme is obtained as folh lows: First, define x -St) ( and x (St) =2x(O) -x ( -St) +$FJx(o)l. ( 43) =x( - St) - (Ar2)Ast) F, [x(O)] (AZ) x x(O)+zAt2 x(At) F(0) 1 (Al3) from which, we obtain the position x( At) = 2x(O) - x ( - At) + -At2 F(0). m In the above, v. is defined in Eq. (7.10). APPENDIX B The following is a piece of FORTRAN code which shows how to implement the integration algorithm expressed in Eq. (3.8). These lines of code assume that the long and short range force components are known from a previous integration step. dt-long is the big time step At in Eq. (3.8), while dt-short corresponds to the small time step St with which the reference system is integrated. Ninner is the number of small time steps per big time step. DO i = l,Natoms VX(i) = VX(i) + OS*dt.Zong*FX-long(i) VY(i) = W(i) -ENDDO DO inner = 1, Ninner DO i = 1,Natoms X(i) = X(i) + dt-short * VX( i) + 0.5 *dt 2-short *FX-short(i) Y(i) = Y(i) + dt-short *VY(i) Z(i) = Z(i) + 0.5 i) (i) *dt 2-short *FY-short( *dt 2~short *FZ-short + OS*dt-Zong*FY-long(i) VZ( i) = VZ( i) + OS*dt. Zong*FZ- Zong( i) (-414) Then, we compute the velocity x(0) from f(O) = x (St) Then, by defining x( -St) = x -St) ( - (A;;st) F,[x(O)], (AS) - x - 60 ( 2St I .I (-44) we obtain the positions at x( At - St) and x( At) x(At-St) x(At) =x,[At-St,x(O),x( -St)], -St)], (A61 (-47) where x, (7) denotes the solution of the reference system using the standard Verlet algorithm for p steps where p = r/St subject to the initial conditions {x(0),x ( - St)}. The need for both x(At - St) and x(At) is cIear from the initial conditions required to initiate each step. As a final example, we show how to obtain the standard Verlet scheme for Nose dynamics as discussed in Sec. VII. For simplicity, we consider the case of F,(x) = F(x) although generalization to a different reference force is straightforward. The initial conditions are given as {x(0),x( - At),v(O),v( - At)}. Again, using the condition in Eq. (2.20), we define x ( =x,[At,x(O),x( + dt-short * VZ( i) + 0.5 *FX-short(i) *FY-short(i) *FZ-short(i) -At) =e - (Ar/2)tio(+ [l-e- Ux( (Ar/l)i~( _ At) -A0 1 I VX( i) = VX( i) + 0.5*dt-short x x(0) + ----HO) At2 2m and x (At) = 2x(O) C-48) W(i) = W(i) + 0.5*dt-short KZ( i) = VZ( i) + 0.5*dt-short ENDDO - x - At) + -At2 F(0). ( m (-49) CALL FORCEs_short DO i = 1,Natoms VX(i) = VX(i) + 0.5*dt-short *FX-short(i) *FY-short(i) *FZ-short(i) The velocity x( 0) is then obtained from f(O) = x (At) - x - At) ( ~-- 2At = The evolution of the bath coordinate is given by the usual Verlet scheme (AlO) VY(i) = VY( i) + 0.5*dt-short m(i) = VZ(i) + 0.5*dt-short ENDDO ENDDO CALL FORCEs_Zong rl(At) = 2q~(O) - v7(-At) + +$Fv (01, (All) J. Chem. Phys., Vol. 97, No. 3, 1 August 1992 Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp Tuckerman, Berne, and Martyna: Reversible multiple time scale MD 2001 DO i = 1,Natoms V.(i) n(i) = VX(i) + 0.5*dt-Zong*FX-long(i) + 0.5*dt-Zong*FY-long(i) + 0.5*dt-Zong*FaZong(i) (Bl) = vZ(i) F Y(i) = W(i) ENDDO In the above, the first loop uses the long range forces from the previous step to obtain the initial conditions on the velocities for the current step. The middle loop runs over the number of inner time steps and computes the reference system positions and velocities from the short range forces. The final loop corrects the velocities from the current values of the and long range force. The routines FORCES-short FORCES-Zong are assumed to update the short and long range components of the force, respectively, dt-short dt 2-short and = dt- Zong/FLO-4 T( Ninner), = (dtwshort)2. M. E. Tuckerman, G. J. Martyna, and B. J. Beme, J. Chem. Phys. 93, 1287 (1990). M. E. Tuckerman, B. J. Berne, and A. Rossi, J. Chea Phys. 94, 1465 (1990). M. E. Tucker-man and B. J. Beme, J. Chem. Phys. 94,681l ( 1991). Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth, Phys. Lett. B S. 195,216 (1987). M. Creutz and A. Goksch, Phys. Rev. Lett. 63,9 ( 1989). De Raedt and B. De Raedt, Phys. Rev. A 28,3575 (1983). H. Takahashi and M. Imada, J. Phys. Sot. Jpn. 53, 3765 (1984). M. aC. W. Gear, NumericalInitial ValueProblems Ordinay Dtrerential in Equations (Prentice-Hall, Englewood Cliffs, 197 1) . 9 H. F. Trotter, Proc. Am. Math Sot. 10, 545 (1959). W C Swope H. C. Andersen, P. H. Berens, and K. R. Wilson, J. Chem. Phys: 76,63+ ( 1982). J. C. Sexton and D. H. Weingarten, Nucl. Phys. B (in press); Nucl. Phys. B Proc. Suppl. 26,613 (1992). M. E. Tuckerman and B. J. Beme, J. Chem. Phys. 95,8362 (1991). I3 M. P. Allen and D. J. Tildesley, ComputerSimulation ofLiquids (Oxford University, Oxford 1989). 14S.Nest, Mol. Phys. 52, 255 (1984). Is S. No&, J. Chem. Phys. 81,5 11 ( 1984). 16W. G. Hoover, Phys. Rev. A 31, 1695 (1985). D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids (Academic, San Diego, 1990). I M. E. Tuckerman and B. J. Beme, J. Chem. Phys. 9.5,4389 (1991). Downloaded 05 Jan 2006 to 128.125.52.12. Redistribution subject to No. 3,1 August 1992 J. Chem. Phys., Vol. 97, AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp Find millions of documents on Course Hero - Study Guides, Lecture Notes, Reference Materials, Practice Exams and more. Course Hero has millions of course specific materials providing students with the best way to expand their education. Below is a small sample set of documents: USC - CS - 653 Computer Physics CommunicationsELSEVIERComputer Physics Communications 83 (1994) 197214Multiresolution molecular dynamics algorithm for realistic materials modeling on parallel computersAiichiro Nakano, Rajiv K. Kalia, Priya VashishtaConcurrent USC - CS - 653 USC - CS - 653 Message Passing Interface (MPI) ProgrammingMPI (Message Passing Interface) is a standard message passing system that enables us to write and run applications on parallel computers. In 1992, MPI Forum was formed to develop a portable message passing USC - CS - 653 OpenMP ProgrammingAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern California Email: ana USC - CS - 653 Hybrid MPI+OpenMP Parallel MDAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern California USC - CS - 653 Parallel Molecular DynamicsThis chapter explains the example parallel MD program, pmd.c, in detail.Spatial Decomposition Spatial decomposition: The physical system to be simulated is partitioned into subsystems of equal volume. Processors in a pa USC - CS - 653 Parallel Molecular DynamicsAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern California E USC - CS - 653 Scalability Metrics for Parallel Molecular DynamicsParallel EfficiencyWe define the efficiency of a parallel program running on P processors to solve a problem of size W. Let T(W, P) be the execution time of this parallel program. Speed of the prog USC - CS - 653 Anton, a Special-Purpose Machine for Molecular Dynamics SimulationDavid E. Shaw, Martin M. Deneroff, Ron O. Dror, Jeffrey S. Kuskin, Richard H. Larson, John K. Salmon, Cliff Young, Brannon Batson, Kevin J. Bowers, Jack C. Chao, Michael P. Eastwood, USC - CS - 653 A Fast, Scalable Method for the Parallel Evaluation of Distance-Limited Pairwise Particle InteractionsDAVID E. SHAWD. E. Shaw Research and Development, LLC and Center for Computational Biology and Bioinformatics, Columbia University, 120 W. 45th S USC - CS - 653 Divide-and-Conquer Parallelization ParadigmDivide-and-Conquer Simulation Algorithms Divide-and-conquer (DC) algorithms: Recursively partition a problem into subprogram of roughly equal size. If subprogram can be solved independently, there is a pos USC - CS - 653 Divide-&amp;-Conquer ParallelismAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern California USC - CS - 653 Computational Materials Science 38 (2007) 642652 www.elsevier.com/locate/commatsciA divide-and-conquer/cellular-decomposition framework for million-to-billion atom simulations of chemical reactionsAiichiro Nakano a,*, Rajiv K. Kalia a, Ken-ichi No USC - CS - 653 DE NOVO ULTRASCALE ATOMISTIC SIMULATIONS ON HIGH-END PARALLEL SUPERCOMPUTERSAiichiro Nakano Rajiv K. Kalia1 Ken-ichi Nomura1 Ashish Sharma1, 2 Priya Vashishta1 Fuyuki Shimojo1,3 4 Adri C. T. van Duin 4 William A. Goddard, III 5 Rupak Biswas Deepak S USC - CS - 653 Load BalancingAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern California Email: anakano USC - CS - 653 Lanczos Method for EigensystemsAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern Californ USC - CS - 653 Supplementary Derivations for the Lanczos-Algorithm LectureSpectral representation The eigenvalues and eigenvectors satisfyn n$ Aij q (&quot; ) = #&quot; qi(&quot; ) = $qi(&quot; ) #%&amp;%&quot; , jj=1%=1()(1)where &quot;#$ = 1 ($= #); 0 ($ % #). Define an orthogona
USC - CS - 653
C3P 913June 1990Performance of Dynamic Load Balancing Algorithms for Unstructured Mesh CalculationsRoy D. WilliamsConcurrent Supercomputing Facility California Institute of Technology Pasadena, CaliforniaAbstract If a nite element mesh has a
USC - CS - 653
SIAM J. SCI. COMPUT. Vol. 20, No. 1, pp. 359392c 1998 Society for Industrial and Applied MathematicsA FAST AND HIGH QUALITY MULTILEVEL SCHEME FOR PARTITIONING IRREGULAR GRAPHSGEORGE KARYPIS AND VIPIN KUMAR Abstract. Recently, a number of researc
USC - CS - 653
USC - CS - 653
Applied Numerical Mathematics 52 (2005) 133152www.elsevier.com/locate/apnumNew challenges in dynamic load balancingKaren D. Devine a, , Erik G. Boman a , Robert T. Heaphy a , Bruce A. Hendrickson a , James D. Teresco b,1 , Jamal Faik c , Joseph E
USC - CS - 653
Hypergraph-based Dynamic Load Balancing for Adaptive Scientic ComputationsUmit V. Catalyurek, Erik G. Boman, Karen D. Devine, Doruk Bozda , Robert Heaphy, g and Lee Ann Riesen Ohio State University Sandia National Laboratories Dept. of Biomedical
USC - CS - 653
Optimizing Molecular DynamicsAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern California
USC - CS - 653
Improving Memory Hierarchy Performance for Irregular Applications*John Mellor-Crummeyt, T Department of Computer Science, MS 132 Rice University 6100 Main Houston, TX 77005 David Whalleyz, Ken Kennedy?Cjohnmc,ken}@cs.rice.edu AbstractThe gap betw
USC - CS - 653
Metrics and Models for Reordering TransformationsMichelle Mills StroutMathematics and Computer Science Division Argonne National Laboratory Argonne, IL 60439 USAPaul D. HovlandMathematics and Computer Science Division Argonne National Laboratory
USC - CS - 653
SIAM REVIEW Vol. 46, No. 1, pp. 345c 2004 Society for Industrial and Applied MathematicsRecursive Blocked Algorithms and Hybrid Data Structures for Dense Matrix Library SoftwareErik Elmroth Fred Gustavson Isak Jonsson Bo K gstrom aAbstract. Ma
USC - CS - 653
Lattice-Boltzmann (LB) Fluid Simulation on a Playstation3 (PS3) ClusterKen-ichi Nomura, Liu Peng &amp; Richard SeymourCollaboratory for Advanced Computing &amp; Simulations Dept. of Computer Science, Dept. of Physics &amp; Astronomy, Dept. of Chemical Engineer
USC - CS - 653
SCOP3A Rough Guide to Scientic Computing On the PlayStation 3Technical Report UT-CS-07-595 Version 0.1by Alfredo Buttari Piotr Luszczek Jakub Kurzak Jack Dongarra George Bosilca Innovative Computing Laboratory University of Tennessee Knoxville Ap
USC - CS - 653
Parallel Lattice Boltzmann Flow Simulation on a Low-cost PlayStation3 ClusterInternational Journal of Computational Science 1992-6669 (Print) 1992-6677 (Online) www.gip.hk/ijcs 2008 Global Information Publisher (H.K) Co., Ltd. All rights reserved.
USC - CS - 653
CTWatchISSN 1555-9874 VOLUME 3 NUMBER 1 FEBRUARY 2007Coming Multicore Revolution and its ImpactJACK DONGARRAGUEST EDITORthe Promise and Perils of theCTWatch Quarterly February 2007IntroductionOver the past few years, the familiar idea t
USC - CS - 653
Teraflops Research ChipMoores Law Motivates Multi-Core10945nm 45nm2007107More, better transistors Continued benefits from Moores Law+ More cores1051032Teraflops of performance operating on Terabytes of dataEntertainmentWhat is
USC - CS - 653
The Landscape of Parallel Computing Research: A View from BerkeleyKrste Asanovic Ras Bodik Bryan Christopher Catanzaro Joseph James Gebis Parry Husbands Kurt Keutzer David A. Patterson William Lester Plishker John Shalf Samuel Webb Williams Katheri
USC - CS - 653
Accelerating Molecular Modeling Applications with Graphics ProcessorsJOHN E. STONE,1* JAMES C. PHILLIPS,1* PETER L. FREDDOLINO,1,2* DAVID J. HARDY,1* LEONARDO G. TRABUCO,1,2 KLAUS SCHULTEN1,2,32Beckman Institute, University of Illinois at Urbana-
USC - CS - 653
Harvesting graphics power for MD simulationsarXiv:0709.3225v1 [cond-mat.other] 20 Sep 2007J.A. van Meel , A. Arnold , D. Frenkel , S.F. Portegies Zwart , R.G. Belleman February 2, 2008Abstract We discuss an implementation of molecular dynamics (M
USC - CS - 653
1 Quantum Dynamics BasicsIn this chapter, we will simulate the dynamics of a particle, such as an electron, which follows the law of quantum mechanics [1]. Basics of the quantum-dynamics (QD) method [2-5] are described, along with corresponding data
USC - CS - 653
Quantum Dynamics BasicsSpectral MethodIn this chapter, we will solve the time-dependent Schrdinger equation using another numerical technique, i.e., the spectral method, which is based on Fourier transformation.1.Discrete Fourier TransformCons
USC - CS - 653
Parallelizing Quantum DynamicsAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern Californi
USC - CS - 653
Multiresolution MethodsAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern California Email
USC - CS - 653
13.10 Wavelet Transforms591The splitting point b must be chosen large enough that the remaining integral over (b, ) is small. Successive terms in its asymptotic expansion are found by integrating by parts. The integral over (a, b) can be done usi
USC - CS - 653
Multiresolution Analysis Using WaveletsHAAR BASIS Consider a one dimensional image on 2 pixels: I[2] = (6, 2). We decompose this information into a smooth and a detailed components. The smooth component is an average of the two intensities:s= (6
USC - CS - 653
19.6 Multigrid Methods for Boundary Value Problems871standard tridiagonal algorithm. Given u n , one solves (19.5.36) for un+1/2 , substitutes on the right-hand side of (19.5.37), and then solves for u n+1 . The key question is how to choose the
USC - CS - 653
Computer Physics CommunicationsELSEVIER Computer Physics Communications 83 (1994) 181196Massively parallel algorithms for computational nanoelectronics based on quantum molecular dynamicsAiichiro Nakano, Priya Vashishta, Rajiv K. KaliaConcurrent
USC - CS - 653
Computer Physics Communications 167 (2005) 151164 www.elsevier.com/locate/cpcEmbedded divide-and-conquer algorithm on hierarchical real-space grids: parallel molecular dynamics simulation based on linear-scaling density functional theoryFuyuki Shi
USC - CS - 653
! ! !43 2 1
USC - CS - 653
Hybrid Particle-Continuum SimulationAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Department of Computer Science Department of Physics &amp; Astronomy Department of Chemical Engineering &amp; Materials Science University of Southern Cal
USC - CS - 653
RAPID COMMUNICATIONSPHYSICAL REVIEW B, VOLUME 64, 161102 RLinear scaling relaxation of the atomic positions in nanostructures Stefan Goedecker, Frederic Lancon, and Thierry Deutsch ` Departement de Recherche Fondamentale sur la Matiere Conden
USC - CS - 653
Computer Physics Communications 138 (2001) 143154 www.elsevier.com/locate/cpcHybrid nite-element/molecular-dynamics/ electronic-density-functional approach to materials simulations on parallel computersShuji Ogata a , Elefterios Lidorikis b , Fuyu
USC - CS - 653
PerspectiveEquation-Free: The Computer-Aided Analysis of Complex Multiscale SystemsIoannis G. KevrekidisDept. of Chemical Engineering, and PACM and Mathematics, Princeton University, Princeton, NJ 08544C. William GearDept. of Chemical Engineer
USC - CS - 653
ARTICLE IN PRESSJournal of the Mechanics and Physics of solids 53 (2005) 16501685 www.elsevier.com/locate/jmpsMultiscale modeling of the dynamics of solids at nite temperatureXiantao Lia, Weinan EbaProgram in Applied and Computational Mathema
USC - CS - 653
Comput. Methods Appl. Mech. Engrg. 196 (2007) 908922 www.elsevier.com/locate/cmaGeneralized mathematical homogenization of atomistic media at nite temperatures in three dimensionsJacob Fish *, Wen Chen, Renge LiRensselaer Polytechnic Institute, T
USC - CS - 653
VOLUME 93, NUMBER 17PHYSICA L R EVIEW LET T ERSweek ending 22 OCTOBER 2004Learn on the Fly: A Hybrid Classical and Quantum-Mechanical Molecular Dynamics Simulation Gabor Csa nyi,1 T. Albaret,3 M. C. Payne,1 and A. De Vita 2,3Cavendish Laborat
USC - CS - 653
Center for Turbulence Research Annual Research Briefs 200597A python approach to multi-code simulations: CHIMPSBy J. U. Schlter, X. Wu, E. v. d. Weide, S. Hahn, u M. Herrmann, J. J. Alonso A N D H. Pitsch1. IntroductionSimulations of real wor
USC - CS - 653
Visualizing Molecular Dynamics IWindowingGraphic LibraryOpenGL OpenGL: A standard, hardware-independent interface to graphics hardware (e.g., SGI InfiniteReality2 Graphics Engine). A set of graphics routines defined on the OpenGL virtual graphics
USC - CS - 653
Scientic Visualization BasicsAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Dept. of Computer Science, Dept. of Physics &amp; Astronomy, Dept. of Chemical Engineering &amp; Materials Science University of Southern California Email: anak
USC - CS - 653
Massive Dataset VisualizationAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Dept. of Computer Science, Dept. of Physics &amp; Astronomy, Dept. of Chemical Engineering &amp; Materials Science University of Southern California Email: anak
USC - CS - 653
Ashish Sharma sharmaa@usc.edu Collaboratory for Advanced Computing and Simulations Dept. of Computer Science Aiichiro Nakano anakano@usc.edu Dept. of Computer Science Rajiv K. Kalia rkalia@usc.edu Dept. of Physics &amp; Astronomy Priya Vashishta priyav@u
USC - CS - 653
From Mesh Generation to Scientic Visualization: An End-to-End Approach to Parallel SupercomputingTiankai Tu Hongfeng Yu Leonardo Ramirez-Guzman Jacobo Bielak Omar Ghattas Kwan-Liu Ma David R. OHallaronAbstractParallel supercomputing has tradition
USC - CS - 653
A Spatially Decomposed Parallel Visualization AlgorithmInternational Journal of Computational Science 1992-6669 (Print) 1992-6677 (Online) Global Information Publisher 2007, Vol. 1, No. 4, 407-421ParaViz: A Spatially Decomposed Parallel Visualiza
USC - CS - 653
Singular Value Decomposition &amp; Data MiningAiichiro NakanoCollaboratory for Advanced Computing &amp; Simulations Dept. of Computer Science, Dept. of Physics &amp; Astronomy, Dept. of Chemical Engineering &amp; Materials Science University of Southern California
USC - CS - 653
Mining Scientic DataNaren Ramakrishnan Department of Computer Science Virginia Tech, VA 24061 Tel: (540) 231-8451 Email: naren@cs.vt.edu Ananth Y. Grama Department of Computer Sciences Purdue University, IN 47907 Tel: (765) 494-6964 Email: ayg@cs.pu
USC - CS - 653
DATAMININGGraph-Based Data MiningDiane J. Cook and Lawrence B. Holder, University of Texas at ArlingtonLARGE AMOUNT OF DATA collected today is quickly overwhelming researchers abilities to interpret the data and discover interesting patterns i