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
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)
properties
Original Verlet algorithm Downside regular verlet algorithm: velocity is not known, worse accuracy. Velocity verlet (Andersen 1983): (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
t t t m t t t t t m t t t t t t f f v v f v r r + Δ + Δ + ≈ Δ + Δ + Δ + ≈ Δ + 2 2
2
r t + Δt
m f t
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???
Dynamics of “well-behaved” classical many-body system is chaotic. Consequence: Trajectories that differ very slightly in their initial conditions diverge exponentially (“Lyapunov instability”)
The Lyapunov disaster in action...
1 10
0 , 0 , 0 , , , , 10
N N N j i N
ε ε ε
−
+ − = r p r p p p p
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:
Monte Carlo algorithms –they therefore yield reliable STATIC properties (“Hybrid Monte Carlo”)
time-evolution???
The Shadow Theorem.
hypothesis.
that are “close to” a REAL trajectory of the many-body system.
trajectories?
– Do not discretize NEWTON’s equation of motion… – ...but discretize the ACTION
– Consider a system that is at a point r0 at time t=0 and at point rt at time t=t, then the system follows a trajectory r(t) such that: is an extremum. The Lagrangian L is defined as:
L(r(t)) = K − U(r)
kinetic energy
For example, if we use cartesian coordinates: What does this lead to? Consider the “true” path R(t), with R(0)=r0 and R(t)=rt. Now, consider a path close to the true path: Then the action S is an extremum if What does this lead to?
r(t) = R(t) + δr(t)
∂S ∂r(t) = 0 for all t
L(r(t)) =
N
1 2mi ˙ r2
i − U(r1, r2, . . . rN)
For a one dimensional system this becomes
t0
L(ti) = K(ti) − U(ti)
imax
L(ti)∆t = 1 2m∆t(xi+1 − xi)2 ∆t2 − U(xi)∆t
Sdisc =
imax
⇤
i=1
m(xi+1 − xi)2 2∆t − U(xi)∆t ⇥
Now do the standard thing: Find the extremum for small variations in the path, i.e. for small variations in all xi.
This will generate a discretized trajectory that starts at time t0 at X, and ends at time t at Xt.
∂Sdisc ∂xi = ∂ ∂xi
imax
⇤
i=1
m(xi+1 − xi)2 2∆t − U(xi)∆t ⇥
∂Sdisc ∂xi = −m(xi+1 − xi) + m(xi − xi−1) ∆t − ∆t∂U(xi) ∂xi
0 = m ∆t
m ∂U(xi) ∂xi ⇥
conditions of a REAL trajectory –both at the beginning and at the endpoint.
(e.g. time-correlation functions, transport coefficients, power spectra...) ...then a “good” MD algorithm (e.g. Verlet) is fine.
0 = 2xi − xi+1 − xi−1 − ∆t2 m ∂U(xi) ∂xi
properties
Lagrangian is sum of two terms
∂L ∂ ˙ r = ∂K ∂ ˙ r = p ∂L ∂r = −∂U ∂r = F ˙ p = ∂L( ˙ r, r) ∂r p = ∂L( ˙ r, r) ∂ ˙ r
Newton : F=ma
L r, ! r
p = ∂L r, ! r
∂! r H r, p
rp − L r, ! r
dH r, p
rp
r
= ! rd p
r
∂L r, ! r
∂r dr + ∂L r, ! r
∂! r d! r ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ! p = ∂L r, ! r
∂r = − ! pdr + ! r dp
dH r, p
∂r ⎛ ⎝ ⎜ ⎞ ⎠ ⎟dr + ∂H ∂p ⎛ ⎝ ⎜ ⎞ ⎠ ⎟dp ! r = ∂H ∂p ! p = −∂H ∂p
Hamilton’s equations of motion {
The Hamiltonian is defined as Hamilton’s equations are then Integrating equations of motion (by Verlet) conserves the Hamiltonian
i
= p m
= −∂U(rN) ∂r
So a solution to the Hamiltonians equation conserves the TOTAL energy
dH(p,r) = ∂H ∂p dp + ∂H ∂r dr ∂H ∂p = ˙ r ∂H ∂r = − ˙ p dH(p,r) dt = ∂H ∂p ˙ p + ∂H ∂r ˙ r = ˙ r ˙ p − ˙ p ˙ r = 0
MD generates the NVE ensemble
How do we sample the canonical ensemble with MD?
properties
Velocity scaling Do we sample the canonical ensemble?
2 1
N B i i
=
req
i i
This is called the isokinetic thermostat
Maxwell-Boltzmann velocity distribution
3 2 2
exp 2 2 P p p m m β β π ⎛ ⎞ ⎡ ⎤ = − ⎜ ⎟ ⎣ ⎦ ⎝ ⎠
2 3
N N N NVT i N
Partition function
2 2
dp p P p p =∫
3m β =
3 2 4 2
d 4 exp 2 2 p p p m m β π β π ⎛ ⎞ ⎡ ⎤ = − ⎜ ⎟ ⎣ ⎦ ⎝ ⎠ ∫
Fluctuations in the momentum squared:
3 2 2
exp 2 2 P p p m m β β π ⎛ ⎞ ⎡ ⎤ = − ⎜ ⎟ ⎣ ⎦ ⎝ ⎠
3 2 2 2 4 2
dp d 4 exp 2 2 p P p p p p p m m β π β π ⎛ ⎞ ⎡ ⎤ = = − ⎜ ⎟ ⎣ ⎦ ⎝ ⎠
3m β =
2 4 4
dp 15 m p P p p β ⎛ ⎞ = = ⎜ ⎟ ⎝ ⎠
2
2 2 4 2 2 2 2 2 p
p p p p σ − =
2 2 2
15 3 2 3 3 m m m β β β − = =
Fluctuations in the temperature
2 2 2 2 2
B
B B k T B B
k T k T k T k T σ − = 2 3N =
2 2 1
1 1 3 3
N B i i
k T p m N p N Nm
=
= =
( ) ( )
2 2 2 2 1
1 3
N B i i
k T p mN
=
⎛ ⎞ = ⎜ ⎟ ⎝ ⎠
( )
4 2 2 2 1 1 1
1 3
N N N i i j i i j j i
p p p mN
= = = ≠
⎛ ⎞ ⎜ ⎟ = + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠
( ) ( )
2 4 2 2
1 1 3 N p N N p mN = − −
( )
2 2 4 2 2 2 2 2 2
1
B
k T B
N p N N p N p k T N p σ − − − =
2 4 2 2 2
1 2 3 p p N N p − = =
So, no the canonical ensemble is not isokinetic
– Nose-Hoover
– Andersen – Langevin – Bussi – Nose Hoover- Langevin All of these alter the velocities such that the trajectory samples the canonical NVT ensemble, and the partition function becomes These thermostats differ in how they achieve this
to collide with the Andersen demon
new velocity
uncorrelated (Poisson distribution)
3 2 2
exp 2 2 P v mv m β β π ⎛ ⎞ ⎡ ⎤ = − ⎜ ⎟ ⎣ ⎦ ⎝ ⎠
; exp P t v vt ν = −
) (
i 2 i 2 i
p p r
L t L t L t
Δ Δ Δ
Velocity Verlet:
e
iLpΔt 2
( ) : v t
2m ! f 0
e
iLrΔt
( ) :r t + Δt
v t + Δt 2
e
iLpΔt 2
( ) : v t + Δt
2m ! f t
Couple to the heat bath
to the dynamics.
– Evolve a single time step with a symplectic integrator e.g. velocity Verlet – Calculate the kinetic energy Kt . – Evolve the kinetic energy Kt =Kt + dK for a time corresponding to a single time step using an auxiliary continuous stochastic dynamics: – Rescale the velocities so as to enforce this new value of the kinetic energy Kt.
relaxation time parameter Wiener noise term
Solving the differential equation gives Where the Ri are independent Gaussian random numbers (Bussi et al 2007)
goal: compute MD trajectory sampling NVT ensemble. Take kinetic energy out of the system and put it back in via a ‘piston’. piston can be seen as additional variable s storing kinetic energy Approach: extended Lagrangian
extended variable
+
effective mass constant to be set
now define then it is possible to show that the partition function Znose is for g=3N+1 the system samples the canonical distribution if p’ is interpreted as the real momentum
2 Nose 2 1
ln 2 2
N N i s i
p p g U r s ms Q β
=
= + + +
H
Znose ∝ 1 N! ⇤ dpNdrN exp
g H(p, r) ⇥
H(p, r) =
N
p
i 2
2mi + U(rN)
p = p/s
= 1 N! d ps
d p'N dr N
ds
s3Nδ H p',r
2
2Q + g β ln s− E # $ % % & ' ( (
2
QNose = 1 N! d ps
d pN dr N
ds
δ HNose − E
= 1 N! d ps
d p'N dr N
ds
βs3N+1 g δ s−e
β L E−H p',r
( )− p2s
2Q # $ % & ' (
) * + + ,
. = 1 N! d ps
d p'N dr N
β g e
β 3N+1
( )
g E−H p',r
( )− p2s
2Q # $ % & ' (
= C N! d p'N dr N
e
− β 3N+1
( )
g H p',r
( )
= C N! d p'N dr N
e
−βH p',r
( )
g = 3N +1
Gaussian integral Constant plays no role in thermodynamics Recall MD MC
QNVE = 1 N! d pN dr N
δ E − H r N , pN
QNVT = 1 N! d pN dr N
exp −βH r N , pN
# $ % & h s
2
2Q + g β ln s− E h' s
β 1 s s0 = exp β g E − H p',r
2
2Q " # $ $ % & ' ' ( ) * + * ,
. * δ H p',r
2
2Q + g β ln s− E " # $ $ % & ' ' βs g δ s−e
β g E−H p',r
( )− ps
2
2Q " # $ $ % & ' '
( ) * * * + ,
equations of motion follow from Hamilton's equations.
2 Nose 2 1
ln 2 2
N N i s i
p p g U r s ms Q β
=
= + + +
H
Nose 2
d d
i i i
r p t p ms ∂ = = ∂ H
Nose
d d
N i i i
U r p t r r ∂ ∂ = − = − ∂ ∂ H
Nose
d d
s s
p s t p Q ∂ = = ∂ H
2 Nose 2
d 1 d
s i
p p g t s s ms β ⎛ ⎞ ∂ = − = − ⎜ ⎟ ∂ ⎝ ⎠
H
NH equation of motion can be rewritten as ( Hoover 1984) Where now denotes a kind of ‘friction’ term Mass Q determines the damping
ξ = ps/Q
i
i /mi − 3N/β
Lennard-Jones fluid mean square displacement temperature relaxation
˙ = ˙ = ( )−ξ ˙ ξ = ( − )/µ−γξ + σ ˙
−∆ ˆ
≈
−(∆ / )ˆ
+
−(∆ / )ˆ + −(∆ / )ˆ + −∆ ˆ
+ −(∆ / )ˆ +
−(∆ / )ˆ + −(∆ / )ˆ
ˆ = ˆ + ˆ + ˆ + ˆ ( , ; + ∆ ) =
−∆ ˆ ( , ; )
µ is mass Friction and noise term
s
= (−γ ∆ ) = s µ ( − ) ξ := ξ + :=
−ξ∆ /
:= + ∆ ( ) := + ∆ ξ := ξ + ∆ ( − )/µ := + ∆ := + ∆ ( ) :=
−ξ∆ /
ξ := ξ +
Note the similarity to Velocity rescaling (Bussi)
(Leimkuhler 2012)
Lagrangian formalism as the Nose Hoover thermostat.
– Andersen barostat, – Martyna-Tuckerman-Klein barostat, – Nosé-Hoover Langevin barostat
– Parrinello Rahman barostat
properties
CH bonds by fixing the the bond length.
m1 r
1 = f1 + g1
m2 r
2 = f2 + g2
m3 r
3 = f3 + g3
r1 r2 r3 d12 d23
ga = 1 2 λ12∇ra χ12 + 1 2 λ23∇ra χ23
χ12 = r
12 2 (t)− d12 = 0
χ23 = r
23 2 (t)− d23 = 0
Lagrange multiplier
Newton’s 3rd law (action = - reaction)
ma ra = fa + ga ≈ fa + ga
(r)
r
a(t + Δt) = r a
"(t + Δt)+(Δt2 / ma)aga
(r)(t)
g1
(r) = λ12r 12
g2
(r) = λ23r 23 − λ12r 12
g3
(r) = −λ23r 23
time averages over a NVT MD trajectory ensemble average Ergodicity theorem states that for an ‘ergodic system’
⇥A⇤ =
Ensemble averages follow from time averages Temperature follows from equipartition: ½ kBT per d.o.f. Where f is number of degrees of freedom
T = 2K kB f K = 1 2 mvi
2 i=1 N
¯ A = 1 T T A(t)dt hAi =
P = NkBT V + 1 3V fijr
ij i<j N
Pressure follows from virial expression
Statistical mechanics gives us Project free energy on collective variable q The result is a Landau free energy This can be generalized to multiple dimensions
F(q) = −kBT ln Z drNδ(q − q(rN))e−βU(rN)
Diffusion equation (Fick’s second law)
Solution for an initial c(x,0)=δ(0): all molecules at origin
hx2(t)i = R dxx2c(x, t) R dxc(x, t) = Z dxx2e− x2
4Dt
p 4πDt = 2Dt
Mean square displacement of the molecules Time derivative gives
General c(x,t)
D = 1 6 lim
t→∞
dhr2(t)i dt
Diffusion in 3 dimensions
Relation to velocity
∂hx2(t)i ∂t = 2 ∂ ∂t Z t dt0 Z t0 dt00hvx(t0)vx(t00)i
Define τ = t –t’’ Green –Kubo relation Also exists for other transport coefficients, such as viscosity and conductivity
V(r) = kr(r − r
eq)2 bonds
+ kθ (θ −θeq)2
angles
+ 1 2 vn(1+ cos(nφ − φ0)
dihedrals
+ aij r
ij 12 − bij
r
ij 6 + qiq j
εr
ij
' ( ) ) * + , ,
i< j
1 2
3
4 2,3 4 1
vdW interactions only between non-bonded |i-j|>4
UCoul = 1 2 qiφ(r
i) i=1 N
φ(r
i) =
qj r
ij + nL j,n
screened potential of the opposite sign, such each charge q is canceled
converges much quicker.
Ucoul = 1 2V 4π k2
k≠0
ρ(k)
2 exp −k2 / 4α
π $ % & ' ( )
1/2
qi
2 + 1
2 qiqjerfc αr
ij
r
ij i≠j N
1=0 N
charge (converges quickly)
by a Fourier series)
(MacKerrel et 96)
(Cornell et al. 95)
(Berendsen et al 87)
(Jorgensen et al 95)
(Levitt et al 83)
atom rep.
– TIP3P, TIP4P – SPC/E
Head N-Terminal C-Terminal Arm1 Arm2 Linker
– Ribosome exit tunnel – Downstream in cytosol
– N-terminal: Tail – C-terminal: Body + Arms – PPIase domain: Head
Biochimica et Biophysica Acta 1803 (2010), 650
MD 250 ns trajectories in 50 mM salt solution Amber FF System with explicit water: 200000 atoms
Probability Histogram gives free energy landscape Caveat: MD not long enough to be converged
Head --- Arm1 N-Terminal --- Arm2
Vreede, A. Mashaghi, S. J. Tans, P. G Bolhuis, PLoS One, 8, e59683 (2013)
βF(x, y) = − ln P(x, y)
– crystal structure collapses – first hydrophobic collapse of Head-Arm1 – followed by hydrophilic interaction of Nterm-Arm2
caused by high barriers Next Thursday: rare events!