经典分子动力学方法

• 类原子内嵌势 Johnson, Mei potential, Glue Potential, Finnis-Sutton Potential • 紧束缚势(Tight Binding) • Tersoff,Brennerd,Stillinger-Weber Potential 半导体
两体势
Lennard-Jones
t (m)
∑ Q({r , v }
i
i t (m)
)....(m = 1,..., M )
积分牛顿方程的方法(I)
1.Verlet法则
mi ri (t + h) − 2ri (t ) + ri (t − h) ∂ 4 = − V ({ r }) + O ( h )......( A) i 2 h ∂ri
−
r & &r 2s s
r r ( f +1) k BT & & Q& s& = ∑ mr ⋅ r s − s
i
引入了新的自由度s,定义与其相关的动能和势能。 与老系统的耦合体现在对速度的重新标定。 f为系统的自由度数，Q为Nose质量。 新的自由度的引入相当与引入了一个恒温热库，系统与 其耦合趋于和热库热平衡。系统的温度恒温热库相同而 保持恒温 Nose, Mol. Phys. 52, 255-68(1984)
首先由方程(A)初步预测新的位置rip(t+h),并且计算 力Fi(t+h)。然后根据方程(B)计算新的速度vi(t+h). 再由方程(C)计算最后新的位置ri(t+h)
原子间相互作用势V({ri})
• 两体势(Lennard-Jones,Morse) 惰性气体、简单金属。 • 原子内嵌势(Embed Atom Methods) 简单金属，过渡金属
j i θijk k
3 f C (rij ) g (θ ijk ) exp[λ3 ( r − r ) ] ik 3 ij
k ≠i , j
∑f
C
(rik ) exp[λ (rij − rik ) ]
c2 d2
g (θ ijk ) = 1 +
−
c2 d 2 −[ h − cos(θ ijk ) 2 ]
Brenner C-H势与Tersoff势基本相同， 情参阅,BrennerÆPRB42,9458(1990)
ri (t + h) − ri (t − h) vi (t ) = + O(h 3 )......( B) 2h
*由前两个时刻的位置，根据方程(A)推得下一个时刻得位置 速度由方程(B)计算速度。 *需要连续记录两个时刻得位置。
积分牛顿方程的方法(II)
2.Verlet 速度法则
ri (t + h) − ri (t ) F (t ) h.........( A) = vi (t ) + h 2mi
ri p (t + h) − ri (t − h) = vi (t ).........( A) 2h vi (t + h) − vi (t ) Fi (t + h) + Fi (t ) mi = .........( B) h 2 ri (t + h) − ri (t ) vi (t + h) + vi (t ) = .........(C ) h 2
1 3V p iα p iα mi
h
一个正立方体中的原子系统，立方体边长为h，s为原子相立 方体的坐标，现在把(V,s)均做为变量。和原先系统(r) 比自由度增加1.(1)中右边2,3项为与V相关的能量。可以 看出等压过程是通过调节立方体的体积与速度来实现的。 Q为Anderson质量，调节Q可以控制压力涨落的大小. Andeson, H.C. J.Chem.Phys. 72, 2384-93(1980)
F (t + h) + F (t ) vi (t + h) = vi (t ) + h........( B) 2mi
需要知道上一个时刻得位置，速度和力，首先由方程(A)计算 新得位置，然后计算新得力F(t+h)，再由由方程(B)计算新时 刻得速度，需要储存前一个时刻的位置，速度和力。
积分牛顿方程的方法(III)
PRB48,22(1993)
Tersoff&Brenner势(I)
E = ∑ f C (rij )[aij f R (rij ) + bij f A (rij )]
i< j
f R = A exp(−λ1r ) f A = − B exp(−λ2 r )
冲击势 吸引势
1, r < R − D 1 1 π (r −R) f C (r ) = 2 − 2 sin[ 2 D ], R − D < r < R + D 0, r > R + D
j ≠i ij i i
• 两体势 • 一般情况下，带正 电的原子实间相互 作用，是冲击势
• 嵌到电子气中原子所受 到的力，是原子所处位 置电子密度的函数F(ρi) 是吸引势。 ρi = ∑ f (rij )
j ≠i
EAM. S.W.Daw and M.I.Baskes, PRL50,1285(1983);PRB29,6443(1984) S.M.Foiles, M.I.Bakes and S.W.Daw, PRB33,7983(1986) Johnson. PRB37,3924(1988);PRB39,12554(1989). Finnis-Sutton. Philos. Mag.50,45(1984),Philos.Mag.Lett,61,139(1990) Glue. F.Ercolessi and J.B.Adams, EuroPhys.Lett.26,583(1994) Mei. J.Mei etal Phys.Rev.B43,4653(1991)
TersoffÆPRB37,6991(1988),PRB39,5566(1989)
Tersoff&Brenner势(II)
三体效应包含在aij和bij中
aij = (1 + α η )
n
n − 21n ij
ς ij =
ηij =
k ≠i , j
∑
bij = (1 + β ζ )
n
3 3 3
n − 21n ij
M.J.Weins,Surf. Sci. 31, 138(1972)
Morse
L.A.Girifalco and G.V.Weizer, Phys.Rev, 114,687(1959)
V = ∑ 4ε [( rij ) − ( rij ) ]
σ 12 σ
6 i< j
V = ∑ ε [e
i< j
− 2 ρ 0 ( rij −1)
H
ne
r r r = H + Α({ri }, { pi }) ⋅ Γ(t ) ⇒ (1)
非平衡态分子动力学(NEMD)
r &= r
r r r & = F − A ⋅ Γ(t ) ⇒ (3) p p
r p m
r + Ar ⋅ Γ(t ) ⇒ (2)
3.Leap Frog法则
1 vi (t + 1 h ) − v ( t − i 2 2 h) mi = Fi (t ).........( A) h ri (t + h) − ri (t ) = vi (t + 1 2 h).........( B ) h
1 v(t + 1 h + v t − ) ( 2 2 h) v(t ) = ........(C ) 2
0
r
− 2e
− ρ 0 ( rij −1)
0
r
]
冲击势
吸引势
冲击(排斥)部分(阻止原子靠的过近）
原子间平衡距离
ac
rc
吸引势部分(约束原子不散开)
其它原子势尽管数学形势不同，基本上由这两部分组成。
原子内嵌势(EAM)
EAM,Johnson,Finnis-Sutton,Glue,Mei
E=
1 2
∑ φ (r ) − ∑ F ( ρ )
需要知道上一个时刻得位置ri(t)，中间时刻的速度 vi(t-(1/2)h)和力Fi(t)，首先由方程(A)计算 vi(t+(1/2)h),然后由方程(B)计算新时刻得位置
ri(t+h)，再由方程(C)计算vi(t)。
积分牛顿方程的方法(IV)
4.Gear Predictor-Corrector法则
1 N 2 H = ∑ mi vi + V ({ri }) 2 i =1
d2 mi 2 r dt i
∂ = − V ({ri}) ∂ri
V({ri})是原子间相互作用势，通过解上面的方
程我们可以得到体系在相空间得由轨迹，进而 求得物理量得平均值[t(1),t(2),t(3),…t(M)]
1 Q= M
分子动力学的主要目的是解上面的方程求得体系状态 相空间演化的轨迹{ripi}t0,{ripi}t1,{ripi}t2,{ripi}t3,……。 进而可计算我们感兴趣的物理量的值Q({ri,pi})。
在实际的应用中，我们把上面的哈密顿方程化为下面的 牛顿方程，并且用位置ri和速度vi做为描述体系的参量。
有限体系的等压分子动力学
L
) P ≡ Pext + P = 0
L = L0 + PextV ({rij })
Pext
L0 = ∑
i
N
P
pi2 2 mi
− φ ({ri })
L0
) P=
1 3V 2 ( m v ∑ i i − ri .∇φ ) i N
1 3V
2 ( m v ∑ i i − ri .∇φ − ri .Pe∇V ) = 0 i N i
e
⊕
⊕
E=
1 2
∑ φ (r ) − ∑ F ( ρ )
j ≠i ij i i
紧束缚势
E = ∑ Aαβ e
i, j − pαβ (
αβ r
0 rij
−1)
− ∑ {∑ ξαβ e
2 i j
− 2 qαβ (
αβ r
0
rij
−1)
}
1 2
紧束缚势与EAM势形式上完全一样，第一项代表离子 间的排斥项，第二项为电子气对原子的相互势，是吸 引项。再进行合金的模拟中，需要确定不同原子间的 紧束缚参数。
经典分子动力学方法
中国科学院固体物理研究所 计算材料科学研究室 范巍
分子动力学基本原理
• 一个体系有N个原子 • 体系的状态由这N个原子的位置{ri}和动量{pi} 或速度{vi}来标志。 • 体系的能量为H({ri,pi}) • 体系的运动方程为
∂ ∂ pi = − H ∂ri ∂t
∂ ∂ ri = H ∂pi ∂t
不需要引入质量参数
Fra Baidu bibliotek
D.Y.Sun and X.G.Gong, J.Phys.Condens.Matter 14 (2002)L487-L493
Parrinello-Rahman等压分子动力学
&2 − P V H = H 0 + Q∑ T αβ ex
1 2
r r − 1 −1 & & & & ms = H F − mG Gs && = ( P − IP )V ( H −1 )T QT
ex
αβ
r r r T = (t1 , t2 , t3 ) r r r V = t1 ⋅ t2 × t3
G = T TT
r r r =T ⋅s
与Anderson方法不同之处在于,不但能改变 元胞的大小又可以改变元胞的形状，附加 的自由度是6个
Parrinello and Rahman Phys.Rev. Lett. 45, 1196-99(1980); J. appl. Phys. 52, 7182-90(1981)
Anderson等压分子动力学
2 & H = H 0 + QV − PexV ⇒ (1) r r & r 2V F & & & ⇒ ( 4) V && = ( P − Pex ) ⇒ (5) s = 1 − 3V s Q 1 2
mV 3
P = ∑( + ∑ rijα Fijα ) j >i r r α ,i r r & ⇒ (3) V = h 3 r = hs ⇒ ( 2) v = hs
分子动力学中的系综
1.微正则系综 2.等温分子动力学 3.等压分子动力学 4.非平衡态分子动力学 5.最速下降(Steep Desend)
Nose-Hover等温分子动力学
& − ( f + 1)k BT ln(s ) H = H 0 + Qs
1 2 2
r r & v = sr
r & r& =
r F ms 2
N
P=
N i
=
n
Pe 3V
∑ r .∇V = − P
i
N i
ext
显然∑ ri .∇V为常数, 则P也为常数, 系统L0为常压的
若选取V是r的齐次函数V (λr ) = λ V (r ), 则∑ ri .∇V = nV 是常数, 从而P是常数
V=
4π 3
∑γ (
i
N
rij 3 i 2
) , 则P = Pe ≡ 常数