Transcript Verlet-I/r-RESPA/Impulse is Limited by Nonlinear Instabilty

Nonlinear Instability in Multiple Time
Stepping Molecular Dynamics
Jesús Izaguirre, Qun Ma,
Department of Computer Science and Engineering
University of Notre Dame
and
Robert Skeel
Department of Computer Science and Beckman Institute
University of Illinois, Urbana-Champaign
SAC’03
March 10, 2003
Supported by NSF CAREER and BIOCOMPLEXITY grants
1
Overview

Background




Nonlinear instabilities





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
2
Overview

Background




Nonlinear instabilities





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
3
Classical molecular dynamics

Newton’s equations of
motion:
Mr ''  U (r)  F(r). - - - (1)



Atoms
Molecules
CHARMM potential
(Chemistry at Harvard Molecular
Mechanics)
Bonds, angles and torsions
4
The CHARMM potential terms
Bond
Angle
Dihedral
Improper
5
Energy functions
6
Overview

Background




Nonlinear instabilities





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
7
Multiple time stepping

Fast/slow force splitting



Bonded: “fast” (small periods)
Long range nonbonded: “slow” (large char. time)
Evaluate slow forces less frequently


Fast forces cheap
Slow force evaluation expensive
Fast forces, t
Slow forces, t
8
Verlet-I/r-RESPA/Impulse
Grubmüller,Heller, Windemuth and Schulten, 1991
Tuckerman, Berne and Martyna, 1992



The state-of-the-art MTS integrator
Fast/slow splitting of nonbonded terms via switching
functions
2nd order accurate, time reversible
half kick: vn 1/ 2  vn  t f slow (rn ) / 2
oscillate: update positions and momentum
using Verlet/leapfrog (t/2, much smaller time steps)
Algorithm 1. Half step discretization of Impulse integrator
9
Overview

Background




Nonlinear instabilities





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
10
Linear instability of Impulse
Linear instability: energy growth occurs unless
longest t < 1/2 shortest period.
Total energy(Kcal/mol) vs. time (fs)
Impulse
MOLLY - ShortAvg
MOLLY - LongAvg
11
Overview

Background




Nonlinear “instabilities” (overheating)





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
12
Nonlinear instability of Impulse
Ma, Izaguirre and Skeel (SISC, 2003)


Approach

Analytical: Stability conditions for a nonlinear model
problem

Numerical: Long simulations differing only in outer time
steps; correlation between step size and overheating
Results: energy growth occurs unless
longest t < 1/3 shortest period.

Unconditionally unstable 3rd order nonlinear resonance

Flexible waters: outer time step less than 3~3.3fs

Constrained-bond proteins w/ SHAKE: time step less
than 4~5fs
13
Overview

Background




Nonlinear instabilities





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
14
Nonlinear instability: analytical

Approach:


1-D nonlinear model problem, in the neighborhood of
stable equilibrium
MTS splitting of potential:
U (q)  (2q2 / 2)oscillate  ( Aq2 / 2  Bq3 / 3  Cq4 / 4)kick  O(q5 ).



Analyze the reversible symplectic map
Express stability condition in terms of problem
parameters
Result:



3rd order resonance stable only if “equality” met
4th order resonance stable only if “inequality” met
Impulse unstable at 3rd order resonance in practice
15
Nonlinear: analytical (cont.)

Main result. Let
   2 , where     i , where
  1  hs ' A /(2c ') if c '  0,   0 if c '  0, and
  1  hc ' A /(2s ') if s '/   0,   0 if s '/   0, and
s '  sin(h / 2) and c '  cos(h / 2).
1. (3rd order) Map stable at equilibrium if
B  0, C  0, and unstable if B  0.

Impulse is unstable in practice.

May be stable at the 4th order resonance.
2. (4th order) Map stable if C  0 or C  2hB2 s ' c '/ ,
and unstable if 0  C  2hB2 s ' c '/ .
16
Overview

Background




Nonlinear instabilities





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
17
Nonlinear resonance: numerical
Fig. 1: Left: Flexible water system. Right: Energy drift from 500ps MD simulation of flexible
water at room temperature revealing 3:1 and 4:1 nonlinear resonance (3.3363 and 2.4 fs)
18
Nonlinear resonance: numerical
Fig. 2. Energy drift from 500ps MD simulation of flexible
water at room temperature revealing 3:1 (3.3363)
19
Nonlinear: numerical (cont.)
Fig. 3. Left: Flexible melittin protein (PDB entry 2mlt). Right: energy drift from 10ns
MD simulation at 300K revealing 3:1 nonlinear resonance (at 3, 3.27 and 3.78 fs).
20
Overview

Background




Nonlinear instabilities





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
21
Concluding remarks

MTS restricted by a 3:1 nonlinear
resonance that causes overheating


Important for long MD simulations due to:



Longest time step < 1/3 fastest normal mode
Faster computers enabling longer simulations
Long time kinetics and sampling, e.g., protein
folding
Use stochasticity for long time kinetics

For large system size, NVE  NVT
22
Overview

Background




Nonlinear instabilities





Classical molecular dynamics (MD)
Multiple time stepping integrator
Linear instability
Analytical approach
Numerical approach
Concluding remarks
Acknowledgements
Key references
23
Acknowledgements

People





Resources



Dr. Thierry Matthey
Dr. Ruhong Zhou, Dr. Pierro Procacci
Dr. Andrew McCammon hosted JI in May 2001 at UCSD
Dept. of Mathematics, UCSD, hosted RS Aug. 2000 –
Aug. 2001
Hydra and BOB clusters at ND
Norwegian Supercomputing Center, Bergen, Norway
Funding


NSF CAREER Award ACI-0135195
NSF BIOCOMPLEXITY-IBN-0083653
24
Key references
[1] Overcoming instabilities in Verlet-I/r-RESPA with the mollified impulse
method. J. A. Izaguirre, Q. Ma, T. Matthey, et al.. In T. Schlick and H. H. Gan,
editors, Proceedings of the 3rd International Workshop on Algorithms for
Macromolecular Modeling, Vol. 24 of Lecture Notes in Computational Science
and Engineering, pages 146-174, Springer-Verlag, Berlin, New York, 2002
[2] Verlet-I/r-RESPA/Impulse is limited by nonlinear instability. Q. Ma, J. A.
Izaguirre, and R. D. Skeel. Accepted by the SIAM Journal on Scientific
Computing, 2002. Available at http://www.nd.edu/~qma1/publication_h.html.
[3] Targeted mollified impulse – a multiscale stochastic integrator for molecular
dynamics. Q. Ma and J. A. Izaguirre. Submitted to the SIAM Journal on
Multiscale Modeling and Simulation, 2003.
[4] Nonlinear instability in multiple time stepping molecular dynamics. Q. Ma, J.
A. Izaguirre, and R. D. Skeel. In Proceedings of the 2003 ACM Symposium on
Applied Computing (SAC’03), pages 167-171, Melborne, Florida. March 9-12,
2003
25
Key references
[5] Long time step molecular dynamics using targeted Langevin
Stabilization. Q. Ma and J. A. Izaguirre. In Proceedings of the 2003
ACM Symposium on Applied Computing (SAC’03), pages 178-182,
Melborne, Florida. March 9-12, 2003
[6] Dangers of multiple-time-step methods. J. J. Biesiadecki and R. D.
Skeel. J. Comp. Phys., 109(2):318–328, Dec. 1993.
[7] Difficulties with multiple time stepping and the fast multipole
algorithm in molecular dynamics. T. Bishop, R. D. Skeel, and K.
Schulten. J. Comp. Chem., 18(14):1785–1791, Nov. 15, 1997.
[8] Masking resonance artifacts in force-splitting methods for
biomolecular simulations by extrapolative Langevin dynamics. A.
Sandu and T. Schlick. J. Comut. Phys, 151(1):74-113, May 1, 1999
26
THE END.
THANKS!
27
Nonlinear: numerical (cont.)
Fig. 4. Left: Melittin protein and water. Right: Energy drift from 500ps SHAKEconstrained MD simulation at 300K revealing combined 4:1 and 3:1 nonlinear resonance.
28