经典分子动力学方法

D.Y.Sun and X.G.Gong, J.Phys.Condens.Matter 14 (2002)L487-L493
Parrinello-Rahman等压分子动力学
&2 H = H 0 + Q ∑ Tαβ PexV
1 2
r r 1 && = H F mG 1Gs && ms & QT& = ( P IP )V ( H 1 )T
首先由方程(A)初步预测新的位置rip(t+h),并且计算 力Fi(t+h)。然后根据方程(B)计算新的速度vi(t+h). 再由方程(C)计算最后新的位置ri(t+h)
原子间相互作用势V({ri})
两体势(Lennard-Jones,Morse) 惰性气体、简单金属。 原子内嵌势(Embed Atom Methods) 简单金属，过渡金属
M.P.Allen and D.J.Tildesley, Computer Simulation of Liquids, Oxford University Press(1987)
j i θijk k
f C (rij ) g (θ ijk ) exp[λ3 (rij rik ) 3 ] 3
C
k ≠i , j
∑f
(rik ) exp[λ (rij rik ) ]
c2 d2
g (θ ijk ) = 1 +
情参阅,Brenner
c2 d 2 [ h cos(θ ijk ) 2 ]
在实际的应用中，我们把上面的哈密顿方程化为下面的 牛顿方程，并且用位置ri和速度vi做为描述体系的参量。
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)]
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
*由前两个时刻的位置，根据方程(A)推得下一个时刻得位置 速度由方程(B)计算速度。 *需要连续记录两个时刻得位置。
积分牛顿方程的方法(II)
2.Verlet 速度法则
ri (t + h) ri (t ) F (t ) h.........( A) = vi (t ) + h 2mi
F (t + h) + F (t ) vi (t + h) = vi (t ) + h........( B) 2mi
有限体系的等压分子动力学
L
) P ≡ Pext + P = 0
L = L0 + PextV ({rij })
Pext
L0 = ∑
i
N
P
pi2 2 mi
φ ({ri })
L0
) P=
1 3V
1 3V
(mi vi2 ri .φ ri .PeV ) = 0 ∑
i Pe 3V
N
P=
N i
(mi vi2 ri .φ ) ∑
0
Baidu Nhomakorabea
r
2e
ρ 0 ( rij 1)
0
r
]
冲击势
吸引势
冲击(排斥)部分(阻止原子靠的过近）
原子间平衡距离
ac
rc
吸引势部分(约束原子不散开)
其它原子势尽管数学形势不同，基本上由这两部分组成。
原子内嵌势(EAM)
EAM,Johnson,Finnis-Sutton,Glue,Mei
E=
1 2
∑ φ (r ) ∑ F ( ρ )
经典分子动力学方法
中国科学院固体物理研究所 计算材料科学研究室 范巍
分子动力学基本原理
一个体系有N个原子 体系的状态由这N个原子的位置{ri}和动量{pi} 或速度{vi}来标志。 体系的能量为H({ri,pi}) 体系的运动方程为
pi = H ri t
ri = H pi t
分子动力学的主要目的是解上面的方程求得体系状态 相空间演化的轨迹{ripi}t0,{ripi}t1,{ripi}t2,{ripi}t3,……。 进而可计算我们感兴趣的物理量的值Q({ri,pi})。
i
N
=
n
∑ r .V = P
i i
N i
N
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 ≡ 常数
不需要引入质量参数
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)
vi (t + 1 h) vi (t 1 h) 2 2 mi = Fi (t ).........( A) h ri (t + h) ri (t ) = vi (t + 1 h).........( B ) 2 h
需要知道上一个时刻得位置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法则
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
H
ne
r r r = H + Α({ri }, { pi }) Γ(t ) (1)
非平衡态分子动力学(NEMD)
r & r=
r r r & p = F Ap Γ(t ) (3)
r p m
r + Ar Γ(t ) (2)
Ar = r A
AP = P A
在外场的作用下，系统内由粒子输运和热输运。系统 一般处于非平衡态下。可用含时哈密顿量来表示(1)。 运动方程由式(2)(3)表示。A是3Nx3N阶矩阵，Γ是3N维 的矢量.
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
Brenner C-H势与Tersoff势基本相同，
PRB42,9458(1990)
分子动力学中的系综
1.微正则系综 2.等温分子动力学 3.等压分子动力学 4.非平衡态分子动力学 5.最速下降(Steep Desend)
Nose-Hover等温分子动力学
& H = H 0 + Qs ( f + 1)k BT ln(s )
Anderson等压分子动力学
& 2 P V (1) H = H 0 + QV ex r r r && = F 1 2V& s (4) V& = ( P Pex ) (5) & & s 3V Q
1 2
mV 3
P = ∑( + ∑ rijα Fijα ) j >i r r α ,i r r & r = hs (2) v = hs (3) V = h 3
需要知道上一个时刻得位置，速度和力，首先由方程(A)计算 新得位置，然后计算新得力F(t+h)，再由由方程(B)计算新时 刻得速度，需要储存前一个时刻的位置，速度和力。
积分牛顿方程的方法(III)
3.Leap Frog法则
v(t + 1 h) + v(t 1 h) 2 2 v(t ) = ........(C ) 2
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势形式上完全一样，第一项代表离子 间的排斥项，第二项为电子气对原子的相互势，是吸 引项。再进行合金的模拟中，需要确定不同原子间的 紧束缚参数。
类原子内嵌势 Johnson, Mei potential, Glue Potential, Finnis-Sutton Potential 紧束缚势(Tight Binding) Tersoff,Brennerd,Stillinger-Weber Potential 半导体
两体势
Lennard-Jones
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)
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)
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)
1 2 2
r r & v = sr
r && = r
r F ms 2
r && 2 sr s
r r ( f +1) k BT & & Q&& = ∑ mr r s s s
i
引入了新的自由度s,定义与其相关的动能和势能。 与老系统的耦合体现在对速度的重新标定。 f为系统的自由度数，Q为Nose质量。 新的自由度的引入相当与引入了一个恒温热库，系统与 其耦合趋于和热库热平衡。系统的温度恒温热库相同而 保持恒温 Nose, Mol. Phys. 52, 255-68(1984)
1 Q= M
t (m)
∑ Q({r , v }
i
i t (m)
)....(m = 1,..., M )
积分牛顿方程的方法(I)
1.Verlet法则
mi ri (t + h) 2ri (t ) + ri (t h) = V ({ri }) + O(h 4 )......( A) h2 ri
ri (t + h) ri (t h) vi (t ) = + O(h 3 )......( B) 2h