Lecture 3: Integration Algorithms Professor A. Martini Purdue - PowerPoint PPT Presentation

Short Course on Molecular Dynamics Simulation Lecture 3: Integration Algorithms Professor A. Martini Purdue University

High Level Course Outline MD Basics 1. Potential Energy Functions 2. Integration Algorithms 3. Temperature Control 4. Boundary Conditions 5. Neighbor Lists 6. Initialization and Equilibrium 7. Extracting Static Properties 8. Extracting Dynamic Properties 9. 10. Non-Equilibrium MD Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms � Potential energy (force) is a function of 3N atomic positions � There is no analytical solution to the equations of motion which must be solved numerically v v v v = d r v = v d v = F m a v i a i i i i i i dt dt Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms � General Rules: – Conservation of energy – Reversible – Computational efficient – Enable a “long” integration time step – Only one force evaluation per time step � Commonly used integrators: – Verlet – Velocity Verlet – Predictor-Corrector – Gear Predictor-Corrector Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms � Error: – Round off error vs. truncation – Local vs. global Round-Off Truncation Error Error Dominant Dominant Total Global Error Integration Time Step Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms � Verlet Algorithm: – Derived from two Taylor expansions dr t d r t d r t 2 3 ( ) 1 ( ) 1 ( ) + δ = + δ + δ + δ + δ r t t r t t t t O t 2 3 4 ( ) ( ) ( ) dt dt dt 2 3 2 3 ! dr t d r t d r t 2 3 ( ) 1 ( ) 1 ( ) − δ = − δ + δ − δ + δ r t t r t t t t O t 2 3 4 ( ) ( ) ( ) dt dt dt 2 3 2 3 ! – Add together and simplify d r t 2 ( ) + δ = − − δ + δ + δ r t t r t r t t t O t 2 4 ( ) 2 ( ) ( ) ( ) dt 2 Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms � Notes on Verlet: – Velocities not explicitly solved, calculated typically from first order central difference + δ − − δ r t t r t t ( ) ( ) = v t ( ) δ t 2 – Position vector at t+ δ t requires positions previous two time steps; a two-step method; not self starting – Advantages: simplicity and good stability – Global error O ( δ t 2 ) Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms � Velocity Verlet Algorithm: – Improved accuracy compared to standard Verlet – Start with position and velocity expansions 1 + δ = + δ + δ + r t t r t v t t a t t 2 ( ) ( ) ( ) ( ) ... 2 1 + δ = + δ + + δ + v t t v t t a t a t t ( ) ( ) [ ( ) ( ) ... 2 Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms Velocity Verlet Algorithm: � – Each integration cycle δ t 1 + = + δ v t v t a t t 1. Calculate velocities at mid- ( ) ( ) ( ) 2 2 step δ ⎛ + ⎞ t + δ = + δ ⎜ ⎟ r t t r t v t t ( ) ( ) 2. Calculate positions at the ⎝ ⎠ 2 next step 3. Calculate accelerations at next step from the potential δ ⎛ + ⎞ t 1 + δ = + + δ δ ⎜ ⎟ v t t v t a t t t 4. Update the velocities ( ) ( ) ⎝ ⎠ 2 2 Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms Predictor-Corrector Algorithms: � 1. Predict positions and velocities at the end of the next timestep 2. Evaluate forces at the next time step using the predicted positions 3. Correct the predicted positions and velocities Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms Predictor-Corrector Algorithms: � 1. Predict the system configuration at the end of the next timestep using Taylor expansion + δ = + r t t r t ( ) ( ) ... + δ = + v t t v t ( ) ( ) ... + δ = + a t t a t ( ) ( ) ... δ + = + b t t b t ( ) ( ) ... Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms Predictor-Corrector Algorithms: � 2. Evaluate forces at the next time step using the predicted system state; difference between the predicted and newly calculated acceleration is the error Δ + δ = + δ − + δ c p a t t a t t a t t ( ) ( ) ( ) Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms Predictor-Corrector Algorithms: � 3. Use the error calculated in the previous step to correct all next step values + δ = + δ + Δ + δ c p r t t r t t c a t t ( ) ( ) ( ) 0 + δ = + δ + Δ + δ c p v t t v t t c a t t ( ) ( ) ( ) 1 + δ = + δ + Δ + δ c p a t t a t t c a t t ( ) ( ) ( ) 2 + δ = + δ + Δ + δ c p b t t b t t c a t t ( ) ( ) ( ) 3 – Coefficients maximize stability and are dependant on the specific algorithm chosen Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms Gear Predictor-Corrector Algorithms: � Predict using 5 th order Taylor series – – So need five derivatives of position at each timestep – Coefficients are tabulated for q-order predictors: • For example, with q=3 C 0 =1/6 – C 1 =5/6 – C 2 =1 – C 3 =1/3 – – Error is O ( δ t q+1 ) Molecular Dynamics Simulation A. Martini Last Updated 8/2009

Integration Algorithms Choosing a time step � Too small � trajectory covers only a limited part of the – phase space Too large � numerical instability – – Timestep should be ~ 1 order of magnitude smaller than the shortest motion time scale – In general: System Types of Motion Time Step (s) 10 -14 Atoms Translation 5x10 -15 Rigid molecules Translation and Rotation Flexible molecules, Translation, Rotation and 2x10 -15 rigid bonds Torsion 10 -15 or 5x10 -16 Flexible molecules, Translation, Rotation, flexible bonds Torsion and Vibration Molecular Dynamics Simulation A. Martini Last Updated 8/2009