Time-Step Methods

In an effort to lengthen the feasible simulation time scale of molecular simulations, Streett and co-workers introduced the multiple time-step method in 1978. These authors recognized that the components of the force that vary most rapidly, and hence require small time steps for numerical resolution, are typically associated with atom pair interactions at small separations. This spatial localization is important because each of the N particles in the simulation has such an interaction with only a few, say C N, neighboring particles. 

Success of the multiple time-step methodology depends on computing these kN C interactions at each step, i.e., at intervals of h Eqs. [3]-[6] and [8], while computing the remaining N N — k) pair interactions (which correspond to forces that vary more slowly) at longer time intervals ih. 

Street et al. originally presented the MTS method in the context of a distance truncated Lennard-Jones potential, so that the total number of computed interactions was somewhat less that N. For biological MD applications there is evidence that cutoffs can cause undesirable artifacts.  

Before discussing the implementation details, we need to state the general issues of multiple time-step numerical methods. The central objectives are (1) to devise a splitting of the systematic forces into a hierarchy of two or more force classes based on the time interval over which they vary significantly, and (2) to incorporate these force classes into a numerical method in a way that realizes enhanced computational efficiency and maintains stability and accuracy of the computed solution. 

---

A straightforward derivation (not reproduced here) shows that the effect of the diree successive steps embodied in equation (b3.3.7), with the above choice of operators, is precisely the velocity Verlet algorithm. This approach is particularly usefiil for generating multiple time-step methods. 

MD, one needs to use multiple time step methods to ensure proper handling of the sprmg vibrations, and there is a possible physical bottleneck in the transfer of energy between the spring system and the other degrees of freedom which must be handled properly [199]. In MC, one needs to use special methods to sample configuration space efficiently [200, 201]. 

In general, multiple-time-step methods increase computational efficiency in a way complementary to multipole methods The latter make use of regularities in space, whereas multiple-time-stepping exploits regularities in time. Figure 2 illustrates the general idea ... 

W. B. Streett, D. J. Tildesley, and G. Saville. Multiple time step methods in molecular dynamics. Mol. Phys., 35 639-648, 1978. 

M. Watanabe and M. Karplus. Dynamics of molecules with internal degrees of freedom by multiple time-step methods. J. Chem. Phys., 99 8063-8074, 1993. 

J. J. Biesiadecki and R. D. Skeel. Dangers of multiple-time-step methods. J. Comp. Phys., 109 318-328, 1993. 

B. Garcia-Archilla, J.M. Sanz-Serna, and R.D. Skeel. Long-time-step methods for oscillatory differential equations. SIAM J. Sci. Comp., 1996. To appear, [Also Tech. Kept. 1996/7, Dep. Math. Applic. Comput., Univ. Valladolid, Valladolid, Spain). 

Another way to overcome the step-size restriction fc < is to use multiple-time-stepping methods [4] or implicit methods [17, 18, 12, 3). In this paper, we examine the latter possibility. But for large molecular systems, fully implieit methods are very expensive. For that reason, we foeus on the general class of scmi-implicit methods depicted in Fig. 1 [12]. In this scheme. Step 3 of the nth time step ean be combined with Step 1 of the (n - - l)st time step. This then is a staggered two-step splitting method. We refer to [12] for further justification. 

Fig. 1. Schematic for the impulse multiple time stepping method.

Watanabe, M., Karplus, M. Dynamics of Molecules with Internal Degrees of Freedom by Multiple Time-Step Methods. J. Chem. Phys. 99 (1995) 8063-8074 Figueirido, F., Levy, R. M., Zhou, R., Berne, B. J. Large Scale Simulation of Macromolecules in Solution Combining the Periodic Fast Multiple Method with Multiple Time Step Integrators. J. Chem. Phys. 106 (1997) 9835-9849 Derreumaux, P., Zhang, G., Schlick, T, Brooks, B.R. A Truncated Newton Minimizer Adapted for CHARMM and Biomolecular Applications. J. Comp. Chem. 15 (1994) 532-555... 

Garcia-Archilla, B., Sanz-Serna, J.M., Skeel, R.D. Long-Time-Steps Methods for Oscillatory Differential Equations. SIAM J. Sci. Comput. (to appear)... 

Abstract. The overall Hamiltonian structure of the Quantum-Classical Molecular Dynamics model makes - analogously to classical molecular dynamics - symplectic integration schemes the methods of choice for long-term simulations. This has already been demonstrated by the symplectic PICKABACK method [19]. However, this method requires a relatively small step-size due to the high-frequency quantum modes. Therefore, following related ideas from classical molecular dynamics, we investigate symplectic multiple-time-stepping methods and indicate various possibilities to overcome the step-size limitation of PICKABACK. 

Here we suggest a different approach that propagates the system using multiple step-sizes, i.e., few steps with step-size At are taken in the slow classical part whereas many smaller steps with step-size 5t are taken in the highly oscillatory quantum subsystem (see, for example, [19, 4] for symplectic multiple-time-stepping methods in the context of classical molecular dynamics). Therefore, we consider a splitting of the Hamiltonian H = Hi +H2 in the following way ... 

We have derived time-reversible, symplectic, and second-order multiple-time-stepping methods for the finite-dimensional QCMD model. Theoretical results for general symplectic methods imply that the methods conserve energy over exponentially long periods of time up to small fluctuations. Furthermore, in the limit m —> 0, the adiabatic invariants corresponding to the underlying Born-Oppenheimer approximation will be preserved as well. Finally, the phase shift observed for symmetric methods with a single update of the classical momenta p per macro-time-step At should be avoided by... 

A different long-time-step method was previously proposed by Garci a-Archilla, Sanz-Serna, and Skeel [8]. Their mollified impulse method, which is based on the concept of operator splitting and also reduces to the Verlet scheme for A = 0 and admits second-order error estimates independently of the frequencies of A, reads as follows when applied to (1) ... 

Table 7.1 presents us with something of a dilemma. We would obviously desire to explore i much of the phase space as possible but this may be compromised by the need for a sma time step. One possible approach is to use a multiple time step method. The underlyir rationale is that certain interactions evolve more rapidly with rime than other interaction The twin-range method (Section 6.7.1) is a crude type of multiple time step approach, i that interactions involving atoms between the lower and upper cutoff distance remai constant and change only when the neighbour list is updated. However, this approac can lead to an accumulation of numerical errors in calculated properties. A more soph sticated approach is to approximate the forces due to these atoms using a Taylor seri< expansion [Streett et al. 1978] ... 

The described method can generate a first-order backward or a first-order forward difference scheme depending whether 0 = 0 or 0 = 1 is used. For 9 = 0.5, the method yields a second order accurate central difference scheme, however, other considerations such as the stability of numerical calculations should be taken into account. Stability analysis for this class of time stepping methods can only be carried out for simple cases where the coefficient matrix in Equation (2.106) is symmetric and positive-definite (i.e. self-adjoint problems Zienkiewicz and Taylor, 1994). Obviously, this will not be the case in most types of engineering flow problems. In practice, therefore, selection of appropriate values of 6 and time increment At is usually based on trial and error. Factors such as the nature of non-linearity of physical parameters and the type of elements used in the spatial discretization usually influence the selection of the values of 0 and At in a problem. 

After application of the 6 time-stepping method (see Chapter 2, Section 2.5) and following the procedure outlined in Chapter 2, Section 2.4, a functional representing the sum of the squares of the approximation error generated by the finite element discretization of Equation (4.118) is formulated as... 

The remaining terms in equation set (4.125) are identical to their counterparts derived for the steady-state case (given as Equations (4.55) to (4.60)). By application of the 9 time-stepping method, described in Chapter 2, Section 2.5, to the set of first-order ordinary differential equations (4.125) the working equations of the solution scheme are obtained. The general form of tliese equations will be identical to Equation (2.111) in Chapter 2,... 

Application of the previously described 6 time-stepping method to Equation (5.11) gives... 

Constraint and multiple time step methods—Section VIII, this chapter. 

Hamiltonian systems. Thus, one has to treat this non-volume-preserving piece of the integrator a bit more carefully. To ensure numerical stability, higher order reversible integration schemes in conjunction with multiple time step methods are preferred. The details of implementing this scheme are provided in Ref. 28. 

The simulation uses a variable-time step method in order to simulate the kinetics over different Pd and PdAu surfaces [85]. The temporal behavior of all intermediates is explicitly tracked throughout the simulation. All atop, bridge, and 3-and 4-fold hollow sites are specifically followed as a function of time. The simulation follows all lateral and through-space interactions between coadsorbed intermediates within a cut off of two nearest-nearest neighbors. 

The major limitation of the approaches to multiscale modeling discussed thus far is the timescale. In each of these examples, there are atomic vibrations (on the order of 10 seconds) that need to be followed. This pins down the total simulation time to 0(10 seconds for reasonable calculations. There are many clever multiple time step methods for improving efficiency (e.g., Nakano 1999) by using a quatemion/normal mode representation for atoms that are simply vibrating or rotating, but this buys only a factor of 0(10). 

---