Advanced Molecular Dynamics Lyopanov instability and the shadow - PowerPoint PPT Presentation

Advanced Molecular Dynamics Lyopanov instability and the shadow theorem (4.3.4) Lagrangian and Hamiltonians (Appendix A) MD at constant temperature (6.1) Constraints (15.1) Computer experiments (4.4)

Contents • Recap of the Verlet algorithms • Lyopunov instability and the shadow theorem • Lagrangian and Hamiltonian approach • MD at constant temperature: Thermostats • Constraints • Ensemble averaging and computing static and dynamic properties • Example of biomolecular simulation

Recap: The Verlet algorithms Original Verlet algorithm ) + Δ t 2 ( ) ≈ 2 r t ( ) − r t − Δ t ( ( ) r t + Δ t m f t Downside regular verlet algorithm: velocity is not known, worse accuracy. Velocity verlet (Andersen 1983): 2 t Δ ( ) ( ) ( ) ( ) r t t r t v t t f t + Δ ≈ + Δ + 2 m t Δ [ ] ( ) ( ) ( ) ( ) v t t v t f t t f t + Δ ≈ + + Δ + 2 m (Is based on Trotter decomposition of Liouville operator formulation, also basis of Multiple time steps). Velocity Verlet has the advantage of allowing multiple time stepping

Is Verlet a good algorithm? Verlet algorithm – is time reversible – does conserve volume in phase space – (is “ symplectic ” ) – does not suffer from energy drift ...but is it a good algorithm? i.e. does it predict the time evolution of the system correctly???

Molecular chaos Dynamics of “ well-behaved ” classical many-body system is chaotic. Consequence: Trajectories that differ very slightly in their initial conditions diverge exponentially ( “ Lyapunov instability ” )

Lyapunov instability The Lyapunov disaster in action... ( N N ) r 0 , p 0 ( ) ( ) ( N ) r 0 , p 0 ,  , p 0 , p 0  , p 0 ( ) ( ) ( ) ( ) ( ) + ε − ε 1 j i N 10 10 − ε =

Any small error in the numerical integration of the equations of motion, will blow up exponentially.... always... ...and for any algorithm!! SO: Why should anyone believe Molecular Dynamics simulations ???

Answers: 1. In fact, one should not … 2. Good MD algorithms (e.g. Verlet) can also be considered as good Monte Carlo algorithms –they therefore yield reliable STATIC properties ( “ Hybrid Monte Carlo ” ) 3. What is the point of simulating dynamics, if we cannot trust the resulting time-evolution??? 4. All is well (probably), because of... The Shadow Theorem.

Shadow theorem • For any realistic many-body system, the shadow theorem is merely a hypothesis. • It basically states that Good algorithms generate numerical trajectories that are “ close to ” a REAL trajectory of the many-body system. • Question: Does the Verlet algorithm indeed generate “ shadow ” trajectories? • Take a different look at the problem. – Do not discretize NEWTON ’ s equation of motion … – ...but discretize the ACTION

Lagrangian Classical mechanics d 2 x i ( t ) F i = m i • Newton: dt 2 • Lagrange: – Consider a system that is at a point r 0 at time t=0 x and at point r t at time t=t, then the system follows a trajectory r(t) such that: � t dt � L ( r ( t � )) S ≡ t=0 t 0 S<S is an extremum. The Lagrangian L is defined as: S<S L ( r ( t )) = K − U ( r ) kinetic energy

Langrangian For example, if we use cartesian coordinates: N 1 � r 2 L ( r ( t )) = 2 m i ˙ i − U ( r 1 , r 2 , . . . r N ) i =1 What does this lead to? Consider the “ true ” path R(t), with R(0)=r 0 and R(t)=r t . Now, consider a path close to the true path: r ( t � ) = R ( t � ) + δ r ( t � ) Then the action S is an extremum if ∂ S ∂ r ( t � ) = 0 for all t What does this lead to?

Discretized action � t 1 S cont = dt L ( t ) t 0 i max � S disc = ∆ t L ( t i ) L ( t i ) = K ( t i ) − U ( t i ) i =0 For a one dimensional system this becomes 2 m ∆ t ( x i +1 − x i ) 2 L ( t i ) ∆ t = 1 − U ( x i ) ∆ t ∆ t 2 i max � m ( x i +1 − x i ) 2 ⇥ ⇤ S disc = − U ( x i ) ∆ t 2 ∆ t i =1

Minimize the action Now do the standard thing: Find the extremum for small variations in the path, i.e. for small variations in all x i . ∂ S disc = 0 for all i ∂ x i This will generate a discretized trajectory that starts at time t 0 at X, and ends at time t at X t .

Minimizing the action i max � m ( x i +1 − x i ) 2 ⇥ ∂ S disc ∂ ⇤ = − U ( x i ) ∆ t ∂ x i ∂ x i 2 ∆ t i =1 = − m ( x i +1 − x i ) + m ( x i − x i − 1 ) − ∆ t ∂ U ( x i ) ∂ S disc ∂ x i ∆ t ∂ x i 2 x i − x i +1 − x i − 1 − ∆ t 2 � ∂ U ( x i ) ⇥ 0 = m ∆ t ∂ x i m

0 = 2 x i − x i +1 − x i − 1 − ∆ t 2 ∂ U ( x i ) ∂ x i m x i +1 = 2 x i − x i − 1 + ∆ t 2 m F ( x i ) • which is the Verlet algorithm! • The Verlet algorithm generates a trajectory that satisfies the boundary conditions of a REAL trajectory –both at the beginning and at the endpoint. • Hence, if we are interested in statistical information about the dynamics (e.g. time-correlation functions, transport coefficients, power spectra...) ...then a “ good ” MD algorithm (e.g. Verlet) is fine.

Contents • Recap of the Verlet algorithms • Lyopunov instability and the shadow theorem • Lagrangian and Hamiltonian approach • MD at constant temperature: Thermostats • Constraints • Ensemble averaging and computing static and dynamic properties • Example of biomolecular simulation

Lagrangian approach Lagrangian is sum of two terms r 2 r ) − U ( r ) = m ˙ L ( ˙ r, r ) = K ( ˙ + U ( r ) 2 ∂ L r = ∂ K ∂ L ∂ r = − ∂ U r = p ∂ r = F ∂ ˙ ∂ ˙ p = ∂ L ( ˙ r, r ) p = ∂ L ( ˙ r, r ) ˙ ∂ ˙ ∂ r r Newton : F=ma

Hamiltonian approach ( ) ( ) ∂ L r , ! ∂ L r , ! r r L r , ! ( ) r ! p = p = ∂ ! r ∂ r ( ) = ! rp − L r , ! ( ) H r , p r ( ) = d ! ( ) − dL r , ! ( ) dH r , p rp r ⎡ ⎤ ∂ L r , ! ( ) ∂ L r , ! ( ) r r = ! ( ) + p d ! ( ) − d ! r d p r dr + r ⎢ ⎥ ∂ ! ∂ r r ⎢ ⎥ ⎣ ⎦ r = ∂ H = − ! p dr + ! ( ) r d p ! ∂ p ⎛ ⎞ ⎛ ⎞ ) = ∂ H ⎟ dr + ∂ H { ( dH r , p ⎟ d p ⎜ ⎜ p = −∂ H ∂ r ∂ p ⎝ ⎠ ⎝ ⎠ ! Hamilton’s equations of motion ∂ p

Hamiltonian approach The Hamiltonian is defined as H ( p, r ) = p ˙ r − L ( ˙ r, r ) p 2 � H ( p N , r N ) = U ( r N ) + i 2 m i i Hamilton ’ s equations are then = p m = − ∂ U ( r N ) ∂ r Integrating equations of motion (by Verlet) conserves the Hamiltonian

Conservation of Hamiltonian dH ( p , r ) = ∂ H ∂ p dp + ∂ H ∂ r dr ∂ H ∂ H ∂ p = ˙ ∂ r = − ˙ r p dH ( p , r ) = ∂ H p + ∂ H ˙ r = ˙ ˙ p − ˙ ˙ ˙ r p r = 0 dt ∂ p ∂ r So a solution to the Hamiltonians equation conserves the TOTAL energy E = K + U MD generates the NVE ensemble How do we sample the canonical ensemble with MD?

Contents • Recap of the Verlet algorithms • Lyopunov instability and the shadow theorem • Lagrangian and Hamiltonian approach • MD at constant temperature: Thermostats • Constraints • Ensemble averaging and computing static and dynamic properties • Example of biomolecular simulation

Constant Temperature: a naïve approach Velocity scaling 3 1 N 1 2 k T mv ∑ = B i 2 N 2 i 1 = T req v v → i i T This is called the isokinetic thermostat Do we sample the canonical ensemble?

Partition function 1 N 2 N ( N ) Q dp exp p 2 m dr exp U r ⎡ ⎤ ∑ ⎡ ⎤ = − β − β ∫ ∫ ⎣ ⎦ NVT i ⎣ ⎦ 3 N h N ! Maxwell-Boltzmann velocity distribution 3 2 β ⎛ ⎞ 2 P p exp p 2 m ( ) ⎡ ⎤ = − β ⎜ ⎟ ⎣ ⎦ 2 m π ⎝ ⎠ 2 2 p dp P p p ( ) = ∫ 3 2 β ⎛ ⎞ 4 2 d 4 p p exp p 2 m ⎡ ⎤ = ⎠ ∫ π − β ⎜ ⎟ ⎣ ⎦ 2 m π ⎝ 3 m = β

3 2 β ⎛ ⎞ 1 N 2 1 P p exp p 2 m ( ) ⎡ ⎤ = − β 2 2 ⎜ ⎟ k T p m N p ∑ ⎣ ⎦ = = 2 m π B i ⎝ ⎠ 3 N 3 Nm i 1 = 3 2 3 m β ⎛ ⎞ 2 2 2 4 2 p dp P p p d 4 p p exp p 2 m ( ) ⎡ ⎤ = ∫ = 1 ∫ N π − β = ⎛ ⎞ ⎜ ⎟ 2 ⎣ ⎦ 2 k T p ( ) 2 m ∑ = π β ⎝ ⎠ ⎜ ⎟ B i 2 3 mN ( ) ⎝ ⎠ i 1 = 2 15 m ⎛ ⎞ 4 4 p dp P p p ( ) = ∫ = ⎜ ⎟ ⎛ ⎞ β N N N 1 ⎝ ⎠ ⎜ ⎟ 4 2 2 p p p ∑ ∑∑ = + i i j Fluctuations in the momentum squared: 2 ⎜ ⎟ 3 mN ( ) ⎜ ⎟ i 1 i 1 j 1 = = = ⎝ ⎠ j i ≠ 2 2 2 2 4 2 p p 1 15 m 3 m ( ) ( ) σ − 2 β − β ( 2 ) 4 2 2 N p N N ( 1 ) p p = − − = = = 2 ( 3 mN ) 2 2 2 3 3 m ( ) 2 2 p p β 2 2 ( ) 4 2 2 N p N N 1 p N p ( ) 2 − − − σ Fluctuations in the temperature k T B = 2 2 k T ( ) 2 N p 2 2 k T k T ( ) B 2 2 − σ B B k T = 2 B = 4 2 p p − 1 3 N 2 2 2 k T k T = = B B 2 N 3 N 2 p So, no the canonical ensemble is not isokinetic