Get documents about "Advanced Molecular Dynamics
Shared by: dandanhuanghuang
-
Stats
- views:
- 2
- posted:
- 11/22/2011
- language:
- English
- pages:
- 43
Document Sample
Advanced Molecular Dynamics
Velocity scaling
Andersen Thermostat
Hamiltonian & Lagrangian Appendix A
Nose-Hoover thermostat
Multiple Timesteps
Car-Parrinello Molecular Dynamics
Constant Temperature
Naïve approach
Velocity scaling
N
3 1 1 2
2
k BT =
N
∑ 2 mvi
i =1
Treq
vi → vi
T
Do we sample the canonical ensemble?
Partition function
QNVT
1
= 3 N ∫ dp N exp − β ∑ pi2 2m ∫ dr N exp − βU r N
h N! ( )
Maxwell-Boltzmann velocity distribution
β
32
P ( p) = exp − β p 2 2m
2π m
p 2 = ∫ dpP ( p ) p 2
β
32
= ∫ dp 4π p 4 exp − β p 2 2m
2π m
3m
=
β
β
32
P ( p) = exp − β p 2 2m
2π m
β
32
p 2 = ∫ dpP ( p ) p 2 =
3m
dp 4π p 4 exp − β p 2 2m =
2π m ∫
β
2
m
p = ∫ dpP ( p ) p = 15
4 4
β
Fluctuations in the momentum:
15 ( m β ) − ( 3 m β )
2 2
σ 2
p − p
4 2 2
2
p2
= = =
(3 m β )
2
2 2 2 2 3
p p
Fluctuations in the temperature
( k BT ) − k BT
2 2
σ 2
2
k BT
= =
k BT
2
k BT
2 3N
N
1 1
k BT =
3N
∑
i =1
p m =2
i
3Nm
N p2
2
1 N
( k BT ) = 2 ∑
2
pi2
( 3mN ) i =1
N N N
1 p4 + pi2 p 2
2 ∑ i ∑∑ j
=
( 3mN ) i =1
i =1 j =1
j ≠i
=
(
1
3mN )
2 (
N p 4 − N ( N − 1) p 2
2
)
N p − N ( N − 1) p ( )
2
2 2
σ 2 4
− N p 2
=
k BT
(N )
2 2
k BT p 2
2 2
1 p − p
4
2
= =
N 2 2 3N
p
Andersen thermostat
Every particle has a fixed probability to collide
with the Andersen demon
After collision the particle is give a new velocity
β
32
P (v) = exp − β mv 2 2
2π m
The probabilities to collide are uncorrelated (Poisson distribution)
P ( t ;ν ) = v exp [ −vt ]
( iL p ∆t 2 ) ∆t &
e : v (t ) → v (t ) + f ( 0)
2m
e (iL r ∆ t ) : r ( t + ∆ t ) → r ( t ) + ∆ t v ( t + ∆ t 2 )
&
(i L p ∆ t ) ∆t &
: v (t + ∆ t ) → v (t + ∆ t 2 ) + f (t )
2
e
2m
Velocity Verlet:
(iLp ∆t 2) ( iLr ∆t ) ( iLp ∆t 2)
e e e
Andersen thermostat:
static properties
Andersen thermostat:
dynamic properties
Hamiltonian & Lagrangian
The equations of motion give the path that starts at t1 at position
x(t1) and end at t2 at position x(t2) for which the action (S) is
the minimum
t2
S = ∫ dt U k − U p
t1
S<S
x
S<S
t1 t2 t
Example: free particle t2 t2
Consider a particle in vacuum: vav ≡ ∫ dtv ( t ) ∫ dt
t1 t1
Up = 0 v(t)=vav t
1 2
1 2
U k = mv
= ∫ dt vav + η ( t )
t2 − t1 t1
2
t2 t2 t
1 2 1 2
S = ∫ dt U k − U p = ∫ dt mv
2
vav = vav + ∫ dtη ( t )
t2 − t1 t1
v(t ) = vav + η ( t )
t1 t1
t2
t2
1
S = m ∫ dt vav + η ( t )
2 t1
2
∫ dtη ( t ) = 0
t1
t2 t2
1
= Sav + m ∫ dtvavη ( t ) + m ∫ dt η ( t )
2
t1
2 t1 η(t)=0 for all t
t
1 2
= Sav + m ∫ dt η ( t )
2
2 t1 Always > 0!!
Lagrangian
Cartesian coordinates ( x, x ) (Newton) →
&
Generalized coordinates ( q, q )
& (?)
Lagrangian
L ( x, x ) = U k ( x ) − U p ( x )
& &
L ( q, q ) = U k ( q ) − U p ( q )
& &
The true path plus deviation
Action t2 t2
S = ∫ dtL ( x, x ) = ∫ dtL ( q, q )
& &
t1 t1
q (t ) = q (t ) +η (t )
d ∂L ( q, q ) ∂L ( q, q )
t2
& &
S[q+η]= S[q] + ∫ dt −
=S + η ( t )
t1 dt ∂q& ∂q
Should be 0 for all paths
t2
d ∂L ( q, q ) ∂L ( q, q )
& &
S[q+η]= S[q] + ∫ dt −
=S + η ( t )
t1 dt ∂q& ∂q
d ∂L ( q, q ) ∂L ( q, q )
& & Equations of motion
=
dt ∂q& ∂q
Lagrangian equations of motion
Conjugate momentum
∂L ( q, q )
& ∂L ( q, q )
&
pq = pq =
&
∂q& ∂q
Newton?
d ∂L ( q, q ) ∂L ( q, q )
& &
L ( q, q ) = U k ( q ) − U p ( q )
& & =
dt ∂q& ∂q
Valid in any coordinate system: Cartesian
1 2
L ( x, x ) = mx − U p ( x )
& &
2
Conjugate momentum
∂L ( x, x )
&
px = = mx
&
∂x&
∂L ( x, x )
& ∂U p ( x )
px =
& =− =F
∂x ∂x
Lagrangian dynamics
With these variables we can do
We have: statistical thermodynamics
2nd order differential equation
( q, q ) → q = L
& &&
Two 1st order differential equations
( q, p ) → q = L ∧ p = L
& &
Change dependence:
( q, q ) → ( q , p )
&
L ( q, q )
& ( q, q ) → ( q , p )
& Hamiltonian
∂L ( q, q )
& ∂L ( q, q )&
p= p=
&
∂q& ∂q
H ( q, p ) = qp − L ( q, q )
& &
dH ( q, p ) = d ( qp ) − dL ( q, q )
&
∂L ( q, q )
& ∂L ( q, q )
&
= qd ( p ) + pd ( q ) −
& & dq + &
dq
∂q ∂q&
= − pdq + qd ( p )
& & ∂H
q=
&
∂H ∂H ∂p
dH ( q, p ) =
∂q
dq +
∂p
dp { p=−
&
∂H
Hamilton’s equations of motion ∂q
Newton?
1 2
L ( x, x ) =
& mx − U p ( x )
&
2
Conjugate momentum
∂L ( x, x )&
p= = mx
&
∂x&
Hamiltonian
H ( x, p ) = xp − L ( x, x ) ∂H p
& & x=
& =
p2 1 2 ∂p m
= − mx − U p ( x )
&
m 2 ∂H ∂U p ( x )
p=−
& =−
=
1 p2
+U p ( x)
∂x ∂x
2 m
Lagrangian Nosé thermostat
( )
N
1 2 2 1 2 g
L Nose = ∑ ms ri − U r − Qs − ln s
& N
&
i =1 2 2 β
Hamiltonian Extended system 3N+1 variables
N
HNose = ∑ ri pi + sps − L ( x, x )
& & &
i =1 Associated mass
Conjugate momentum
∂L ∂L
pi = = ms ri2
& ps = = Qs&
∂ri& ∂s&
pi2
( )
N
ps g
HNose = ∑ + + U r + ln s
N
i =1 2ms
2
2Q β
p' = p s
ps g
= HNose ( p ', r ) + + ln s
2Q β
Nosé and thermodynamics
ps 2 L
HNose = HNose ( p ', r ) + + ln s
2Q β
1
dps ∫ dp N ∫ dr N ∫ dsδ ( HNose − E )
N!∫
QNose =
1 ps 2 L
= ∫ dps ∫ dp ' ∫ dr ∫ dss δ H ( p ', r ) + 2Q + β ln s − E
N N 3N
N!
1 βs 3 N +1 β
E − H ( p ', r ) − s
p
= ∫ dps ∫ dp ' ∫ dr ∫ ds
N N
δ s−e L 2Q
N! L
Gaussian integral
β ( 3 N +1) ps
N β
1 E − H ( p ', r ) −
= ∫ dps ∫ dp ' ∫ dr L e
N L 2Q
N!
β ( 3 N +1)
C − H ( p ', r ) L = 3N + 1
= ∫ dp ' ∫ dr e
N N L
N!
C
dp 'N ∫ dr N e − β H ( p ',r )
N!∫
=
Constant plays no role in thermodynamics
Recall
MD QNVE =
1
N!∫
( (
dp N ∫ dr N δ E − H r N , p N ))
MC QNVT =
1
N!∫ (
dp N ∫ dr N exp − β H r N , p N
)
Delta functions
∫ dsh ' ( s )δ ( h ( s ) ) = ∫ d h ( s )δ ( h ( s ) )
= ∫ dsδ ( s − s0 )
h ( s0 ) = 0
∫ dsh ' ( s )δ ( h ( s ) ) = ∫ dsδ ( s − s )0
δ (s − s )
δ ( h ( s )) = 0
h '( s)
δ ( s − s0 )
δ ( h ( s )) =
h '( s)
ps L
δ H ( p ', r ) + + ln s − E
2Q β
ps L
h ( s ) = H ( p ', r ) + + ln s − E
2Q β
g1 β ps
h '( s) = s0 = exp E − H ( p ', r ) −
β s L 2Q
βs β
E − H ( p ', r ) −
ps
δ s−e L 2Q
L
Lagrangian
Equations of Motion
( )
N
1 2 2 1 g
L Nose = ∑ ms r &i − U r N − Qs 2 − ln s
&
i =1 2 2 β
Hamiltonian
pi2 ps 2
( )
N
g
HNose =∑ + + U r + ln s
N
i =1 2ms
2
2Q β
Conjugate momenta
∂L ∂L
pi = = ms 2 ri
& ps = = Qs
&
∂ri
& ∂s&
Equations of motion:
dri ∂HNose
=
pi
= 2 dpi
=−
∂HNose
=−
( )
∂U r
N
dt ∂pi ms dt ∂ri ∂ri
ds ∂HNose ps dp s ∂HNose 1 2
pi g
= = =− = ∑ 2 −
dt ∂ps Q dt ∂s s ms β
Nosé Hoover
Multiple Timesteps
Time evolution
Liouville formulation
f (pN ,r N )
€
Beware: this solution is
equally useless as the
differential equation!
Solution
21
Time evolution
∂ ∂
˙
iL ≡ iLr + iLp = r ˙
+p
∂r ∂p
f (t) = exp(iLr t)f (0)
∂
˙
= exp r(0)t f (0)
∂r
∞
(˙ (0)t)n ∂ n
r
= f (0)
n=0
n! ∂rn
= f pN (0), (r(0) + r(0)t)N
˙
Shift of coordinates
˙
r(0) → r(0) + r(0)t
22
Time evolution
∂ ∂
˙
iL ≡ iLr + iLp = r ˙
+p
∂r ∂p
f (t) = exp(iLr t)f (0)
∂ f (t) = exp(iLp t)f (0)
˙
= exp r(0)t f (0)
∂r ∂
∞
(˙ (0)t)n ∂ n
r ˙
= exp p(0)t f (0)
= f (0) ∂p
n! ∂rn ∞
n=0 ˙
(p(0)t)n ∂ n
= f (0)
= f pN (0), (r(0) + r(0)t)N
˙ n=0
n! ∂pn
= f (p(0) + p(0)t)N , rN (0)
˙
Shift of coordinates Shift of momenta
r(0) → r(0) + r(0)t
˙ p(0) → p(0) + p(0)t
˙
22
Time evolution
Multiple time steps
• What to use for stiff potentials:
–Fixed bond-length: constraints (Shake)
–Very small time step
27
Time evolution
∂ F ∂
Multiple
iL ≡ iLr + iLp = v +
∂r m ∂v Time steps
iL ≡ iLshort + iLlong
Fshort ∂ Flong ∂
iLshort = iLlong =
m ∂v m ∂v
Trotter expansion:
ei(Llong +Lshort +Lr )∆t ≈ eiLlong ∆t/2 ei(Lshort +Lr )∆t eiLlong ∆t/2
Introduce: δt=Δt/n
n
iLlong ∆t/2 iLshort δt/2 iLr δt iLshort δt/2
≈e e e e eiLlong ∆t/2
28
Time evolution
Car-Parrinello Molecular Dynamics
Another Extended Lagrangian Method
Ab Initio Molecular Dynamics
Born-Oppenheimer
• Instantaneous relaxation to electronic ground state
• No coupling ionic and true electronic dynamics
¨ ∂
MI RI = − E(RN ) + Vnn(RN )
∂RI
• Ionic forces from electronic structure calculation
• Knowledge of electronic properties and its time evolution
• Electronic structure methods:
- DFT
- Hartree Fock, MCSCF
- Tight binding, semi-emperical
Beyond Born-Oppenheimer
• Surface hopping
• Time dependent DFT
Born Oppenheimer Molecular Dynamics (DFT)
• Kohn-Sham expression for electronic energy
E(RN ) = min E KS [n(r), RN ]
n
1 Determine E KS by direct minimization or self-consistent diagonalization
of HKS
2 Evaluate force in ground state using Hellman-Feynman theorem∗
N
n(r) = no(r) = i=1 ψi,0|ψi,0
∂ ∂
− E(RN ) = dr n0(r) Vext(r, RN )
∂RI Ω ∂RI
3 Propagate in time
4 Repeat from 1 on
• Verify that dynamics is performed properly by checking constants of mo-
tion (Energy).
Born-Oppenheimer MD issues
• Preserving constant of motion requires high accuracy of EKS
• Competition between accuracy and computational cost
Car–Parrinello MD
Define dynamical system with both nuclei and Kohn-Sham orbitals as degrees
of freedom.
Car-Parrinello Lagrangian LCP
˙ ˙
LCP (RN , RN , ψ N , ψ N ) =
1 ˙2 ˙ ˙ N
I 2 MI RI + i µ ψi |ψi − EKS [n, R ] − ij Λij ( ψi|ψi − δij )
N
n(r) = i=1 ψi|ψi
Equations of motion
d ∂LCP ∂LCP d δLCP δLCP
= =
˙
dt ∂ RI ∂RI ˙
dt δ ψi | δ ψi |
¨ ∂E KS ∂
MI RI = − + Λij ψi | ψj
∂RI ij
∂RI
¨ δE KS
µ | ψi = − + Λij | ψj
δ ψi | j
Car–Parrinello MD: Characteristics
˙ ˙
LCP (RN , RN , ψ N , ψ N ) =
1 ˙2 ˙ ˙ N
I 2 MI RI + i µ ψi |ψi −EKS [n, R ]− ij Λij ( ψi|ψi −δij ) =
Kn + Ke − EKS − constraints.
• µ fictitious “electronic mass”
• Λij Lagrange multipliers to ensure orthonormality of orbitals
• Time propogation of ionic positions and orbitals simulateously (in BOMD
subsequently).
• System not exactly in the ground state -> EKS not ground state energy
• HCP = Kn + Ke + EKS is constant of motion
• if Ke is ∼constant, Kn+EKS is approximately constant.
• if Ke is small, H = Kn+EKS is near to BO surface:
“low electronic temperature” or “cold electrons”.
CPMD: Conservation of Energy
Car-Parrinello Total energy rigorlously conserved by construction
• Forces on nuclei
∂L ∂EKS
=−
∂RI ∂RI
• Forces on electrons
δL
= −HKS ψi
δψi∗
• Imposing constraints
So Car-Parrinello total energy HCP conserved provided the time propagation
algorithm is accurate
Does Car-Parrinello MD work?
• Adiabatic decoupling
Imposing “cold electrons” and Ke is ∼constant decoupling of the dynamics
of the ions and the orbitals (electrons).
There is no net energy transfer from ions to “kinetic energy” of orbitals.
Adiabatic decoupling achieved by non-overlap of frequency spectrum of
ionic and orbital motion.
• Pictorial view
Example∗: Conservation of Energy
a
Si-crystal: Time evolution of various components of the energy
a∗
Pastore, Smargiassi, and Buda, Phys.Rev. A44 (1991)
Example∗: Deviation of BO surface
Si-crystal: Deviations of forces of BO surface are small and oscillating
Car-Parrinello Method: Summary
• Car-Parrinello method can yield very stable dynamical trajectories, pro-
vided the electrons and ions are adiabatically decoupled
• The method is best suited for e.g. liquids and large molecules with an
electronic gap
• The speed of the method is comparable or faster than using Born-
Oppenheimer dynamics — and still more accurate (i.e. stable)
