经典分子动力学方法
中国科学院固体物理研究所
计算材料科学研究室
范巍
分子动力学基本原理
? 一个体系有N个原子
? 体系的状态 由这N个原子的位置{r
i
}和动量{p
i
}
或速度{v
i
}来标志。
? 体系的能量为H({r
i
,p
i
})
? 体系的运动方程为
H
r
p
t
i
i
?
?
?=
?
?
H
p
r
t
i
i
?
?
=
?
?
分子动力学的主要目的是解上面的方程求得体系状态
相空间演化的轨迹 {r
i
p
i
}
t0
,{r
i
p
i
}
t1
,{r
i
p
i
}
t2
,{r
i
p
i
}
t3
,……。
进而可计算我们感兴趣的物理量的值 Q({r
i
,p
i
})。
在实际的应用中,我们把上面的哈密顿方程化为下面的
牛顿方程,并且用位置 r
i
和速度 v
i
做为描述体系的参量。
})({
2
1
1
2
i
N
i
ii
rVvmH +=
∑
=
})({
2
2
i
i
i
dt
d
m rV
r
r
i
?
?
?=
V({r
i
})是原子间相互作用势,通过解上面的方
程我们可以得到体系在相空间得由轨迹,进而
求得物理量得平均值 [t(1),t(2),t(3),…t(M)]
),...,1....()},({
1
)(
)(
MmvrQ
M
Q
mt
mtii
==
∑
积分牛顿方程的方法(I)
1.Verlet法则
))......((})({
)()(2)(
4
2
AhOrV
rh
htrtrhtr
m
i
i
iii
i
+
?
?
?=
?+?+
))......((
2
)()(
)(
3
BhO
h
htrhtr
tv
ii
i
+
??+
=
*由前两个时刻的位置,根据方程 (A)推得下一个时刻得位置
速度由方程 (B)计算速度。
*需要连续记录两个时刻得位置。
积分牛顿方程的方法(II)
2.Verlet 速度法则
).........(
2
)(
)(
)()(
Ah
m
tF
tv
h
trhtr
i
i
ii
+=
?+
)........(
2
)()(
)()( Bh
m
tFhtF
tvhtv
i
ii
++
+=+
需要知道上一个时刻得位置,速度和力,首先由方程(A)计算
新得位置,然后计算新得力F(t+h),再由由方程(B)计算新时
刻得速度,需要储存前一个时刻的位置,速度和力。
积分牛顿方程的方法(III)
3.Leap Frog法则
)().........(
)()(
2
1
Bhtv
h
trhtr
i
ii
+=
?+
)().........(
)()(
2
1
2
1
AtF
h
htvhtv
m
i
ii
i
=
??+
)........(
2
)()(
)(
2
1
2
1
C
htvhtv
tv
?++
=
需要知道上一个时刻得位置 r
i
(t),中间时刻的速度
v
i
(t-(1/2)h)和力F
i
(t),首先由方程(A)计算
v
i
(t+(1/2)h),然后由方程(B)计算新时刻得位置
r
i
(t+h),再由方程(C)计算 v
i
(t)。
积分牛顿方程的方法(IV)
4.Gear Predictor-Corrector法则
).........(
2
)()()()(
B
tFhtF
h
tvhtv
m
iiii
i
++
=
?+
)().........(
2
)()(
Atv
h
htrhtr
i
i
p
i
=
??+
).........(
2
)()()()(
C
tvhtv
h
trhtr
iiii
++
=
?+
首先由方程(A)初步预测新的位置r
i
p
(t+h),并且计算
力F
i
(t+h)。然后根据方程(B)计算新的速度v
i
(t+h).
再由方程(C)计算最后新的位置r
i
(t+h)
原子间相互作用势 V({r
i
})
? 两体势 (Lennard-Jones,Morse)
惰性气体、简单金属 。
? 原子内嵌势 (Embed Atom Methods)
简单金属,过渡金属
? 类原子内嵌势
Johnson, Mei potential,
Glue Potential, Finnis-Sutton Potential
? 紧束缚势 (Tight Binding)
? Tersoff,Brennerd,Stillinger-Weber Potential
半导体
两体势
Lennard-Jones Morse
M.J.Weins,Surf. Sci. 31, 138(1972)
L.A.Girifalco and G.V.Weizer, Phys.Rev, 114,687(1959)
])()([4
612
ijij
rr
ji
V
σσ
ε ?=
∑
<
]2[
)1()1(2
0
0
0
0
??
<
??
∑
?=
r
ij
r
r
ij
r
eeV
ji
ρρ
ε
吸引势
冲击势
c
r
c
a
阻止原子靠的过近)部分排斥冲击 ()(
原子间平衡距离
)(约束原子不散开吸引势部分
部分组成。势不同,基本上由这两其它原子势尽管数学形
原子内嵌势 (EAM)
EAM,Johnson,Finnis-Sutton,Glue,Mei
? 嵌到电子气中原子所受
到的力,是原子所处位
置电子密度的函数 F(ρ
i
)
是 吸引势。
∑ ∑
≠
?=
iji
iij
FrE )()(
2
1
ρφ
∑
≠
=
ij
iji
rf )(ρ
? 两体势
? 一般情况下,带正
电的原子实间相互
作用,是 冲击势
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)
e
∑ ∑
≠
?=
iji
iij
FrE )()(
2
1
ρφ
⊕
⊕
紧束缚势
∑∑∑
????
?=
ij
q
ji
p
r
ij
r
r
ij
r
eeAE
2
1
00
}{
)1(2
2
,
)1(
αβ
αβ
αβ
αβ
αβαβ
ξ
紧束缚势与 EAM势形式上完全一样,第一项代表离子
间的排斥项,第二项为电子气对原子的相互势,是吸
引项。再进行合金的模拟中,需要确定不同原子间的
紧束缚参数。
PRB48,22(1993)
Tersoff&Brenner势 (I)
∑
<
+=
ji
ijAijijRijijC
rfbrfarfE )]()()[(
)exp(
2
rBf
A
λ??=
吸引势
)exp(
1
rAf
R
λ?=
冲击势
?
?
?
?
?
+>
+<<??
?<
=
?
DRr
DRrDR
DRr
rf
D
Rr
C
,0
],sin[
,1
)(
)(
22
1
2
1 π
Tersoff ?PRB37,6991(1988),PRB39,5566(1989)
Tersoff&Brenner势 (II)
三体效应包含在 a
ij
和 b
ij
中
n
n
ij
n
ij
b
2
1
)1(
?
+= ζβ
n
n
ij
n
ij
a
2
1
)1(
?
+= ηα
∑
≠
?=
jik
ikijijkijCij
rrgrf
,
33
3
])(exp[)()( λθ?
i
j
k
θ
ijk
∑
≠
?=
jik
ikijikCij
rrrf
,
33
3
])(exp[)( λη
])cos([
22
2
2
2
1)(
ijk
hd
c
d
c
ijk
g
θ
θ
??
?+=
Brenner C-H势与 Tersoff势基本相同,
情参阅 ,Brenner ?PRB42,9458(1990)
分子动力学中的系综
1.微正则系综
2.等温分子动力学
3.等压分子动力学
4.非平衡态分子动力学
5.最速下降 (Steep Desend)
Nose-Hover等温分子动力学
rsv
&
rr
=
)ln()1(
2
2
1
0
sTkfsQHH
B
+?+=
&
∑
+
??=
i
s
Tkf
B
srrmsQ
)1(
&
r
&
r
&&
s
rs
ms
F
r
&
r
&
r
&&
r
2
2
?=
引入了新的自由度s,定义与其相关的动能和势能。
与老系统的耦合体现在对速度的重新标定。
f为系统的自由度数,Q为Nose质量。
新的自由度的引入相当与引入了一个恒温热库,系统与
其耦合趋于和热库热平衡。系统的温度恒温热库相同而
保持恒温
Nose, Mol. Phys. 52, 255-68(1984)
Anderson等压分子动力学
)1(
2
2
1
0
??+= VPVQHH
ex
&
)4(
3
2
3
1
??= ss
V
V
mV
F &
r
&&
r
&
r
)5(
)(
?=
?
Q
PP
ex
V
&&
h
)2(?= shr
rr
)3(?= shv
&
rr
∑ ∑
>
+=
iij
ijij
m
pp
V
FrP
i
ii
,
3
1
)(
α
αα
αα
3
hV =
一个正立方体中的原子系统,立方体边长为h,s为原子相立
方体的坐标,现在把(V,s)均做为变量。和原先系统(r)
比自由度增加1.(1)中右边2,3项为与V相关的能量。可以
看出等压过程是通过调节立方体的体积与速度来实现的。
Q为Anderson质量,调节Q可以控制压力涨落的大小.
Andeson, H.C. J.Chem.Phys. 72, 2384-93(1980)
有限体系的等压分子动力学
0=+≡ PPP
ext
)
ext
P
L
P
})({
0 ijext
rVPLL +=
0
L
∑
?=
N
i
i
m
p
rL
i
i
})({
2
0
2
φ
0)..(
2
3
1
=????=
∑
N
i
eiiii
V
VPrrvmP φ
)
∑
??=
N
i
iii
V
rvmP ).(
2
3
1
φ
ext
N
i
i
V
P
PVr
e
?=?
∑
.
3
为常压的系统也为常数则为常数显然
0
,,. LPVr
N
i
i∑
?
是常数从而是常数则的齐次函数是若选取 PnVVrrVrVrV
N
i
i
n
,.),()(
∑
=?=λλ
常数则 ≡==
∑ e
N
i
r
i
PPV
ij
,)(
3
23
4
γ
π
=
不需要引入质量参数
D.Y.Sun and X.G.Gong, J.Phys.Condens.Matter 14 (2002)L487-L493
Parrinello-Rahman等压分子动力学
sTr
rr
?=VPTQHH
ex
?+=
∑
αβ
αβ
2
2
1
0
&
),,(
321
tttT
rrr
=
sGmGFHsm
&
&
r
&&
r
11 ??
?=
321
tttV
rrr
×?=
T
ex
HVIPPTQ )()(
1?
?=
&&
TTG
T
=
与 Anderson方法不同之处在于 ,不但能改变
元胞的大小又可以改变元胞的形状,附加
的自由度是 6个
Parrinello and Rahman Phys.Rev. Lett. 45, 1196-99(1980);
J. appl. Phys. 52, 7182-90(1981)
非平衡态分子动力学 (NEMD)
)1()(}){},({ ?Γ?Α+= tprHH
ii
ne
r
rr
)2()( ?Γ?+= tAr
r
m
p
r
&
r
r
)3()( ?Γ??= tAFp
p
rr
&
r
AA
rr
?=
AA
PP
?=
在外场的作用下,系统内由粒子输运和热输运。系统
一般处于非平衡态下。可用含时哈密顿量来表示 (1)。
运动方程由式 (2)(3)表示。 A是 3Nx3N阶矩阵, Γ是 3N维
的矢量 .
M.P.Allen and D.J.Tildesley, Computer Simulation of Liquids, Oxford University Press(1987)
Langevin分子动力学
)()( tGrmrVrm
iiiii
+???=
&
r
&&
r
γ
Gauess白噪声
)()()(
''
ttCtGtG ?>=< δ
通过调节噪声强度 C来控制系统的温度
K.W.Jacabsen, et.al Phys.Rev. B35, 7423(1987)
最速下降 (Steep Descend)
1
v
r
0
1
=v
r
2
2
1
112
ttvrr
m
F
δδ
r
rrr
++=
MD
0
1
=v
r
0
2
=v
r
0
4
=v
r
稳定构形
F
r
0
3
=v
r
2
2
1
12
trr
m
F
δ
r
rr
+=
Steep Descend
VF ??=
r
由于力是势能的梯度,方向垂直与势能面,并且沿着能量减小的方向,
因此沿着力的方向运动,势能总是减小,最后趋向一个极小点
物理量的计算
A.扩散系数的计算
∑∑∑
===
?+=?+=
M
m
N
i
mimi
NM
s
N
i
ii
N
sRstRsRstRtR
11
2
11
1
2
1
2
)]()([)]()([)(
rrrr
t
tR
t
D
6
)(
2
lim
??
∞→
=
.代表对初试点的平均
s
为了使均方位移与初始点 (s)的选取无关,通常选取
多个初始点 (M个 ),然后再取平均。
,6)(,
2
DttR =子系统对于处于平衡态的多原
α
DttR =)(
2
处于非平衡态的系统
)(1
),(sup,1
onsubdiffusi
nerdiffusio
是亚扩散
是超扩散
<
>
α
α
也可以称为扩散系数慢的标志
做为扩散快但我们仍然可以用
,
,1 D≠α
分布距离是飞行的特征,其跳跃的
其扩散轨迹具有对于团簇的扩散
Levy
Levy,
1<α
1=α
1>α
3',1
)(,)(
'
<<
==
??
μμ
μμ
CttPCllP
taillongLevy ?分布
玻尔兹曼分布
B.速度关联函数
∑∑
=
∑
∑
+
=
=
=
==
N
i
svsv
svstv
N
N
i
vv
vtv
N
M
m
mimi
M
m
mimi
ii
ii
tC
1
)()(
)()(
1
1
)0()0(
)0()(
1
1
1
)(
∑∑
=
∑
∑
+
=
=
=
==
N
i
svsv
svstv
N
N
i
vv
vtv
N
M
m
mimi
M
m
mimi
ii
ii
tC
1
)()(
)()(
1
1
)0()0(
)0()(
1
1
,,
1
,,
,,
,,
)(
βα
βα
βα
βα
αβ
速度关联函数是相同原子不同时间的关联 ,N是你所感兴趣的原子数,
可以是一个原子,或多个原子。我们也可以研究,速度
不同分量间的关联。 α=x, β=y
∫
= dtttCC
P
)cos()()(
2
ωω
π
声子普 (声子态密度)
∑
≠
??=
ji
ji
rrrrg |)|()(
0
δρ
1)(lim =
∞→
rg
r
0)(lim
0
=
→
rg
r
0)(lim =
∞→
rg
r
0)(lim
0
=
→
rg
r
)(4
2
rgrRCF π=
)(rg
r
1)( =rg
)(rg
r
C.对相关函数
无限系统
有限系统
径向分布函数
RCF是实验可以直接测量的量 (中子衍射, X射线衍射)
D.结构因子 s(k)
∑
???
=
ji
rrki
N
ji
ekS
,
)(
1
2
)(
rr
r
r
)(kS
r
1
有序程度
需要选择一系列倒格失 k进行计算或选择某一特殊的 k
0
3
2
1
.............
3
3
2
2
1
1
ρ=→
→
→
→
∞
∞
∞
V
N
V
N
V
N
V
N
r
r
r
r
?
?
?
?
?
?
?
?
11
Vr ?→
22
Vr ?→
33
Vr ?→
1
N?
2
N?
3
N?
∑
??=
ji
ji
N
rrtrtrG
,
1
)|)0()((|),(
2
δ
E.Van Hove 函数
相同原子间的关联,不同原子间的关联,时间演化效应
∑
??=
i
ii
N
rrtrtrG )|)0()((|),(
2
1
δ
Van Hove自相关函数
∑
≠
??=
ji
ji
N
rrtrtrG )|)0()((|),(
2
1
δ
随时间演化的对关联函数
)(rg
对时间求和得到对关联函数
∫
=
c
r
kr
kr
drtrGtkF
0
)sin(
),(),(
中子非弹性散射
The Intermediate Scattering Function
角分布函数
∑∑
=<
?
??Θ?Θ=
N
kji
rr
rr
kjki
N
kjki
kjki
c
rrrrrA
1
1
)(cos)()(),(cos
rr
θδθ
∑∑
=<
?Θ?Θ=
N
kji
kjkiC
rrrrN
1
)()(
2
1
),(cos rA θ
),(cos rA θ
0
2
1
θcos θcos
2
1
?
1?0
分布函数 P(Q)
∫
∑
==
=
dQQQPQQ
step
N
i
i
T
)(
1
1
∫
=1)( dQQP
如果分子动力学运行 N
step
步,那么它在相空间采样 N
step
个点 ,其中 P(Q)是这
Nstep个采样点在区间 Q+dQ中出现的概率。其中 Q可以是 能量 或其它你感兴趣
的物理量。可以从比较少的采样中获得系统的有用信息
∑
=
=?
N
i
ii
vmK
3
1
2
2
1
动能
)(EP
B
Nk
K
T
3
2
=?温度
∑∑
>=
?+=?
ji
ijij
N
i
ii
V
FrvmP )(
1
2
3
1
r
rr
压力是
E
1
E
2
E
一个假想系统的能量分布函
数,能量在 E
1
,E
2
附近有两个
峰值,预示系统有两个理想
的状态
2
2
2
)(
2
3
2
3
)1(
TkN
KK
C
Nk
B
B
?
=??
υ
比热
分
子
动
力
学
的
程
序
结
构
Verlet列表技术
f
r
skin
r
1
d
2
d
0
t
0
t
?+
0
t
?+
0
t
).(,
2121
dddd >中两个最大的是对初始时刻位置偏离
有变化。
本身没列表重新计算,但列表的变化大于
离原子相对于中心原子距没有变化。例如图中的
列表有可能尽管列表就全都重新计算,
整个对原子该条件被破坏,序设计中,只要任何一
在实际的程列表不需要重新计算。都小于
间的距离变化也就是说任何两个原子
。如果变化不会大于任何两个原子间的距离
Verletr
c
VerletVerlet
Verltr
rdd
dd
skin
skin
skin
,
,
,
21
21
<+
+
c
c
c
r
)(rV
r
算量大大增加。来的复杂性有可能使计会更加复杂性,由次带
计算列表的条件算列表的次数,但判断我们可以尽一步减小计
周期性边界条件
)int(
)int(
)int(
boz
z
ijij
boy
y
ijij
box
x
ijij
ij
ij
ij
anbozzz
anboyyy
anboxxx
×?=
×?=
×?=
)|(|)(
)|(|)(
)|(|)(
2
1
2
1
2
1
??=
??=
??=
ijijijij
ijijijij
ijijijij
zzsigzz
yysigyy
xxsigxx
θ
θ
θ
是质心
ccc
ijijij
boz
zz
ij
boy
yy
ij
box
xx
ij
zyx
zyx
z
y
x
cij
cij
cij
,,
),,(
2
1
2
1
?≥≥
?
?
?
?
?
?
2
3
2
1
2
1
||
1
1
≤
+=?<
?=>
x
xxx
xxx
要求
则如果
则如果
期性边界条件。
耗时比较少的周
0
2
1
2
1
?
box
boy
θ
t
tt δ+
Cii
N
i
i
N
C
Rrr
rR
r
rr
r
r
?=
=
∑
=
1
1
1
Cii
N
i
i
N
C
vvv
vv
rrr
rr
?=
=
∑
=
1
1
1
JI
vrmJ
N
i
iii
r
r
rr
r
?=
×=
?
=
∑
1
1
)(
ω
trrrr
rvvv
iiii
iiii
δω
ω
?×+=?
×+=?
)(
112
112
rrrrr
rrrr
质心移至原点 相对于质心速度
整体转动角速度
以后的位置和速度
扣除平动和转动
与以削除。
由度的影响可暂时体平动和转动等外部自
自由度感兴趣,整质时,我们只对其内部
究团簇等体系的性统的平动和转动。在研
差也有可能引起系分运动方程所带来的误
动和平动。另外积统始终会有一个整体转
,那么系速度不为果初始状态的角速度和
角动量均守恒,如在孤立系统中,动量、
0
∑
?+++=
N
i
iiiiii
xxxxxmI ])1()[(
2
3
2
2
2
1 βα
δ
αβαβ
αβ
δ
?
?
?
?
?
?
?
?
?
?
??
??
??
333231
232221
131211
III
III
III
扣除整体平动与转动