4. Molecular dynamics Understanding Molecular Simulation Molecular - PowerPoint PPT Presentation

4. Molecular dynamics Understanding Molecular Simulation

Molecular Simulations ➡ Molecular dynamics : solve equations of motion r 1 r 2 r n ➡ Monte Carlo : importance sampling r 1 r 2 r n Understanding Molecular Simulation

Molecular Dynamics 4. Molecular Dynamics 4.1.Introduction 4.2.Basics 4.3.Liouville formulation 4.4.Multiple time steps Understanding Molecular Simulation

4. Molecular dynamics 4.2 Basics Understanding Molecular Simulation

“Fundamentals” Theory: F = md 2 r dt 2 • Compute the forces on the particles • Solve the equations of motion • Sample after some # of timesteps Understanding Molecular Simulation

4. Molecular dynamics 4.3 Some practical details Understanding Molecular Simulation

Molecular Dynamics Initialization Total momentum should be zero (no external forces) • Temperature rescaling to desired temperature • Particles start on a lattice • Force calculations Periodic boundary conditions • Order NxN algorithm, • Order N: neighbor lists, linked cell • Truncation and shift of the potential • Integrating the equations of motion Velocity Verlet • Kinetic energy • Understanding Molecular Simulation

Molecular Dynamics Algorithm 3 (A Simple Molecular Dynamics Program) simple MD program program md initialization call init t=0 MD loop do while (t.lt.tmax) determine the forces call force(f,en) integrate equations of motion call integrate(f,en) t=t+delt sample averages call sample enddo stop end Understanding Molecular Simulation

3. Molecular dynamics: practical details 3.3.1 Initialization Understanding Molecular Simulation

Algorithm 4 (Initialization of a Molecular Dynamics Program) initialization of MD program subroutine init sumv=0 sumv2=0 do i=1,npart place the particles on a lattice x(i)=lattice pos(i) give random velocities v(i)=(ranf()-0.5) velocity center of mass sumv=sumv+v(i) kinetic energy sumv2=sumv2+v(i)**2 enddo velocity center of mass sumv=sumv/npart mean-squared velocity sumv2=sumv2/npart scale factor of the velocities fs=sqrt(3*temp/sumv2) set desired kinetic energy and set do i=1,npart velocity center of mass to zero v(i)=(v(i)-sumv)*fs position previous time step xm(i)=x(i)-v(i)*dt enddo return end Understanding Molecular Simulation

3. Molecular dynamics: practical details 3.3.2 Force calculation Understanding Molecular Simulation

Understanding Molecular Simulation

Periodic boundary conditions Understanding Molecular Simulation

The Lennard-Jones potential s • The Lennard-Jones potential ⎡ ⎤ 12 6 ⎛ ⎞ ⎛ ⎞ σ − σ ( ) = 4 ε U LJ r ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ r r ⎢ ⎥ ⎣ ⎦ • The truncated Lennard-Jones potential ( ) ⎧ U LJ r r ≤ r c ⎪ ( ) = LJ r U TR ⎨ r > r c 0 ⎪ ⎩ • The truncated and shifted Lennard-Jones potential ( ) − U LJ r c ( ) ⎧ U LJ r r ≤ r c ⎪ ( ) = LJ U TR − SH r ⎨ r > r c 0 ⎪ ⎩ Understanding Molecular Simulation

The Lennard-Jones potentials Understanding Molecular Simulation

Saving CPU-time Cell list Verlet-list Understanding Molecular Simulation

3. Molecular dynamics: practical details 3.3.3 Equations of motion Understanding Molecular Simulation

Understanding Molecular Simulation

Equations of motion We can make a Taylor expansion for the positions: ( ) ( ) d 2 r t dr t Δ t 2 ( ) ( ) = r t ( ) + r t + Δ t Δ t + 2! + O Δ t 3 dt 2 dt The simplest form (Euler): ( ) ( ) = r t ( ) + v t ( ) Δ t + O Δ t 2 r t + Δ t ( ) df t ( ) = v t ( ) + m v t + Δ t Δ t dt We can do better! Understanding Molecular Simulation

We can make a Taylor expansion for the positions: ( ) ( ) ( ) d 2 r t d 2 r t dr t Δ t 2 Δ t 3 ( ) ( ) = r t ( ) + r t + Δ t Δ t + 2! + 3! + O Δ t 4 dt 2 dt 2 dt ( ) ( ) ( ) d 2 r t d 2 r t dr t Δ t 2 Δ t 3 ( ) ( ) = r t ( ) − r t − Δ t Δ t + 2! − 3! + O Δ t 4 dt 2 dt 2 dt When we add the two: ( ) d 2 r t ( ) ( ) + r t − Δ t ( ) = 2 r t ( ) + Δ t 2 + O Δ t 4 r t + Δ t dt 2 ( ) Δ t 2 ( ) ( ) = 2 r t ( ) − r t − Δ t ( ) + f t Verlet algorithm r t + Δ t m + O Δ t 4 no need for velocities numerically not ideal Understanding Molecular Simulation

( ) Δ t 2 ( ) ( ) = 2 r t ( ) − r t − Δ t ( ) + f t Verlet algorithm: r t + Δ t m + O Δ t 4 Velocity Verlet algorithm ( ) Δ t 2 ( ) ( ) = r t ( ) + v t ( ) Δ t + f t r t + Δ t 2 m + O Δ t 4 ( ) + Δ t ( ) = v t ( ) + f t ( ) ⎡ ⎤ v t + Δ t 2 m f t + Δ t ⎣ ⎦ to see the equivalence: ) Δ t 2 ( ) = r t + Δ t ( ) + v t + Δ t ( ) Δ t + f t + Δ t ( r t + 2 Δ t 2 m ( ) Δ t 2 ( ) = r t + Δ t ( ) − v t ( ) Δ t − f t r t 2 m adding the two Δ t 2 ( ) = 2 r t + Δ t ( ) − r t ( ) + v t + Δ t ( ) − v t ( ) ( ) − f t ( ) ⎡ ⎤ ⎡ ⎤ r t + 2 Δ t ⎦Δ t + f t + Δ t ⎣ ⎣ ⎦ 2 m ( ) + Δ t ( ) = v t ( ) + f t ( ) with ⎡ ⎤ v t + Δ t 2 m f t + Δ t ⎣ ⎦ ) Δ t 2 ( ) = 2 r t + Δ t ( ) − r t ( ) + f t + Δ t ( r t + 2 Δ t m Understanding Molecular Simulation

Lyaponov instability MD: reference trajectory ( ) ( ) , ! , r N 0 ( ) , p 1 0 ( ) , ! , p N 0 ( ) r 1 0 with initial condition: ( ) ( ) , ! , r N 0 ( ) , p 1 0 ( ) , ! , p i 0 ( ) + ε , p j 0 ( ) − ε , ! , p N 0 ( ) MD: compare: r 1 0 ε = 10 − 10 Understanding Molecular Simulation

4. Molecular dynamics: 4.4 Liouville Formulation Understanding Molecular Simulation

the dot above, ḟ , Liouville formulation implies time derivative Let us consider a function that f which depends on the ( ) positions and momenta of the particles: f p N , r N We can “solve” how f depends on time: ⎛ ⎞ ⎛ ⎞ ∂ f ∂ f ! f = r + ⎟ ! ⎟ ! p ⎜ ⎜ ∂ r ∂ p ⎝ ⎠ ⎝ ⎠ Define the Liouville operator: ⎛ ⎞ ⎛ ⎞ ∂ ∂ iL ≡ ! ⎟ + ! r p ⎜ ⎜ ⎟ ∂ r ∂ p ⎝ ⎠ ⎝ ⎠ df the time dependence follows from: dt = iLf with solution: beware: the solution is ( ) f = e iLt f 0 equally useless as the differential equation Understanding Molecular Simulation

In an ideal world it would be less useless: ⎛ ⎞ ⎛ ⎞ ∂ ∂ iL ≡ ! ⎟ + ! r p ⎜ ⎜ ⎟ ∂ r ∂ p ⎝ ⎠ ⎝ ⎠ ⎛ ⎞ ∂ Let us look at half the equation iL r ≡ ⎟ ! r ⎜ ∂ r ⎝ ⎠ which has as solution: ( ) iL r t f 0 f = e e x = 1 + x + x 2 2! + x 3 Taylor expansion: 3! + ! 2 + 1 3 + … ⎡ ⎤ ( ) = 1 + iL r t + 1 ( ) ( ) ( ) iL r t f 0 e 2 iL r t 3! iL r t ⎥ f 0 ⎢ ⎣ ⎦ ⎡ ⎤ 2 ⎛ ⎞ ⎛ ⎞ ( ) ∂ ∂ ⎟ + 1 ( ) = 1 + ! ( ) t ( ) t ( ) the operator iL r 2 ⎢ ⎥ iL r t f 0 + … ! e r 0 r 0 f 0 ⎜ ⎜ ⎟ gives a shift of ∂ r ∂ r ⎝ ⎠ ⎝ ⎠ 2 ⎢ ⎥ ⎣ ⎦ the positions ( ) ( ) 2 ⎛ ⎞ ⎛ ⎞ ∂ f 0 2 ∂ f 0 ( ) = f 0 ( ) ⎟ + 1 ( ) t ( ) + ! ( ) t ( ) t f 0 + ! + ! ! r 0 r 0 r 0 ⎜ ⎜ ⎟ ∂ r ∂ r 2 ⎝ ⎠ ⎝ ⎠ ( ) ( ) = f 0 + ! ( ) t Hence: iL r t f 0 e r 0 Understanding Molecular Simulation

The operation iL r gives a shift of the positions ⎛ ⎞ ⎛ ⎞ ∂ ∂ iL ≡ ! ⎟ + ! r p ⎜ ⎜ ⎟ ∂ r ∂ p ⎝ ⎠ ⎝ ⎠ ⎛ ⎞ ∂ Similarly for the operator iL p iL p ≡ ⎟ ! p ⎜ ∂ p ⎝ ⎠ which has as solution: ( ) iL p t f 0 f = e Taylor expansion: ⎡ ( ) ( ) ⎤ ( ) = 1 + iL p t + 1 + 1 ( ) 2 3 iL p t f 0 + … e 2 iL p t 3! iL p t ⎥ f 0 ⎢ ⎣ ⎦ ⎡ ⎤ 2 ⎛ ⎞ ⎛ ⎞ ( ) ∂ ∂ ⎟ + 1 ( ) = 1 + ! ( ) t ( ) t ( ) the operator iL p 2 iL p t f 0 ⎢ ⎥ + … ! e p 0 p 0 f 0 ⎜ ⎜ ⎟ gives a shift of ∂ p ∂ p ⎝ ⎠ ⎝ ⎠ 2 ⎢ ⎥ ⎣ ⎦ the momenta ( ) ( ) 2 ⎛ ⎞ ⎛ ⎞ ∂ f 0 2 ∂ f 0 ( ) = f 0 ( ) ⎟ + 1 ( ) t ( ) + ! ( ) t ( ) t f 0 + ! + ! ! p 0 p 0 p 0 ⎜ ⎜ ⎟ ∂ p ∂ p 2 ⎝ ⎠ ⎝ ⎠ ( ) ( ) = f 0 + ! ( ) t Hence: iL p t f 0 e p 0 Understanding Molecular Simulation

The operation iL r gives a shift of the ( ) ( ) = f 0,0 + ! ( ) t positions: iL r t f 0,0 e r 0 … and the operator iL p a shift of the ( ) ( ) = f 0 + ! ( ) t ,0 iL p t f 0,0 momenta: e p 0 This would have been useful if the operators would commute ( ) t f 0,0 ( ) = e ( ) ≠ e ( ) iL r + iL p iL p t f 0,0 iL r t e e iLt f 0,0 Trotter expansion: we have the non-commuting operators A and B: e A + B ≠ e A e B then the following expansion holds: ( ) P e A + B = lim P →∞ e A B A 2 P e P e 2 P Understanding Molecular Simulation

( ) ( ) = f 0,0 + ! ( ) t iL r t f 0,0 e r 0 ( ) ( ) = f 0 + ! ( ) t ,0 iL p t f 0,0 e p 0 We can apply the Trotter expansion: ( ) P e A + B = lim P →∞ e A B A 2 P e P e 2 P iL p t Δ t = t iL r t Δ t P = iL r Δ t 2 P = iL p P 2 These give as operations: ( ) = f p t ( ) ( ) , r t ( ) ( ) , r t ( ) + ! ( ) Δ t iL r Δ t f p t e r t gives us a shift of the position: ( ) → r t ( ) + ! ( ) Δ t r t + Δ t r t ⎛ ⎞ ( ) = f p t ( ) Δ t ( ) , r t ( ) ( ) + ! ( ) iL p Δ t 2 f p t e p t 2 , r t ⎜ ⎟ ⎝ ⎠ gives us a shift of the momenta: ( ) Δ t ( ) → p t ( ) + ! p t + Δ t p t 2 Understanding Molecular Simulation