Computer Physics Communications 176 (2007) 271–291 www.elsevier.com/locate/cpc Error and timing analysis of multiple time-step integration methods for molecular dynamics Guowen Han a , Yuefan Deng a,b,∗ , James Glimm a , Glenn Martyna c a Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, NY 11794-3600, USA b Institute of Scientific Computing, Nankai University, Tanjin 300071, PR China c IBM Research, Physical Sciences Division, IBM T.J. Watson Laboratory, P.O. Box 218, Yorktown Heights, NY 10598, USA Received 15 September 2006; accepted 4 October 2006 Available online 16 November 2006 Abstract Molecular dynamics simulations of biomolecules performed using multiple time-step integration methods are hampered by resonance instabili- ties. We analyze the properties of a simple 1D linear system integrated with the symplectic reference system propagator MTS (r-RESPA) technique following earlier work by others. A closed form expression for the time step dependent Hamiltonian which corresponds to r-RESPA integration of the model is derived. This permits us to present an analytic formula for the dependence of the integration accuracy on short-range force cutoff range. A detailed analysis of the force decomposition for the standard Ewald summation method is then given as the Ewald method is a good candidate to achieve high scaling on modern massively parallel machines. We test the new analysis on a realistic system, a protein in water. Under Langevin dynamics with a weak friction coefficient (ζ = 1 ps −1 ) to maintain temperature control and using the SHAKE algorithm to freeze out high frequency vibrations, we show that the 5 fs resonance barrier present when all degrees of freedom are unconstrained is postponed to ≈ 12 fs. An iso-error boundary with respect to the short-range cutoff range and multiple time step size agrees well with the analytical results which are valid due to dominance of the high frequency modes in determining integrator accuracy. Using r-RESPA to treat the long range interactions results in a 6× increase in efficiency for the decomposition described in the text.  2006 Elsevier B.V. All rights reserved. Keywords: Multiple time-step; Ewald methods; Coulomb forces; Molecular dynamics; Error analysis 1. Introduction Classical molecular dynamics (MD) simulation studies have provided invaluable insight into the properties of biomolecular sys- tems. Unfortunately, biomolecular simulations are computationally intensive due to the large number of particles and the complex interactions present. In a MD calculation, the trajectory of a group of interacting atoms is determined through the approximate solution of a Hamiltonian set of equations of motion based on an empirical potential function or force field. A typical force field contains bonded interactions through bond, angle, and torsion potentials, and non-bonded interactions through van der Waals and the electrostatic potential terms. The first four terms are short-range in nature, while the last term is long-range in nature. The efficient computation of the long range interactions is recognized to be a challenging task in atomistic computer simulation studies. Intuitively, the long-range interactions can be calculated for each pair of particles in the system. However, the computing time complexity would then be on the order of O(N 2 ). This naïve method becomes impossibly computationally intensive when periodic boundary conditions are applied and/or systems are large. Historically, long-range forces were spherically truncated in macromolecular simulations because of insufficient computer power. However, such a scheme has been shown to produce severe artifacts in the results of the simulation [39,47,48]. Thus the need for a correct and efficient treatment of long-range electrostatic * Corresponding author. E-mail address: Yuefan.Deng@StonyBrook.edu (Y. Deng). 0010-4655/$ – see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2006.10.005