NONLINEAR RESONANCE ARTIFACTS IN MOLECULAR-DYNAMICS SIMULATIONS

The intriguing phenomenon of resonance, a pronounced integrator-induce d corruption of a system's dynamics, is examined for simple molecular systems subject to the classical equations of motion. This source of t imestep limitation is not well appreciated in general, and certainly a nalyses of resonance patterns have been few in connection to biomolecu lar dynamics, Yet resonances are present in the commonly used Verlet i ntegrator, in symplectic implicit schemes, and also limit the scope of current multiple-timestep methods that are formulated as symplectic a nd reversible, The only general remedy to date has been to reduce the timestep. For this purpose, we derive method-dependent timestep thresh olds (e.g., Tables 1 and 2) that serve as useful guidelines in practic e for biomolecular simulations, We also devise closely related symplec tic implicit schemes for which the limitation on the discretization st epsize is much less severe. Specifically, we design methods to remove third-order, or both the third-and fourth-order, resonances. These sev ere low-order resonances can lead to instability or very large energie s. Our tests on two simple molecular problems (Morse and Lennard-Jones potentials), as well as a 22-atom molecule, N-acetylalanyl-N'-methyla mide, confirm this prediction; our methods can delay resonances so tha t they occur only at larger timesteps (EW method) or are essentially r emoved (LIM2 method), Although stable for large timesteps by this appr oach, trajectories show large energy fluctuations, perhaps due to the coupling with other factors that induce instability in complex nonline ar systems. Thus, the methods developed here may be more useful for co nformational sampling of biomolecular structures. The analysis present ed here for the blocked alanine model emphasizes that one-dimensional analysis of resonances can be applied to a more complex, multimode sys tem to analyze resonance behavior, but that resonance due to frequency coupling is more complex to pinpoint. More generally, instability, ap parently due to numerically induced has been observed in the applicati on of the implicit midpoint scheme to vibrating structures and could b e expected also in the simulation of nonlinear wave phenomena; in such applications it is adequate not to resolve the highest frequency mode s, so the proposed methods could be very useful. (C) 1998 Academic Pre ss.