分子动力学模拟方法简介

∆r
=
rc 1
(t
+
∆t)
−
rp 1
(t
+
∆t)
Values c0
c1
3
5/12
1
4
3/8
1
5 251/720 1
6 95/288
1
c2
c3
c4
c5
1/2
3/4
1/6
11/12 1/3
1/24
25/24 35/72 5/48 1/120
② 二阶运动方程之一：
&r& = f (r)
∆r
=
rc 2
(t
计算当前时刻的速度：
vi (t)
=
vi (t
+
1 2
∆t) + 2
vi (t
-
1 2
∆t)
vr
t-∆t/2 t
v
t+∆t/2 t+∆t t+3∆t/2 t+2∆t
Leap-frog算法的表述：
算法启动
① 规定初始位置
② 规定初始速度
③ 扰动初始速度：
v(−∆t/2) = v(0) − ai (0) ∆t/2
② 规定初始速度 ③ 计算第n+1步的位置： ④ 计算第n+1步的力 ⑤ 计算第n+1步的速度： ⑥ 重复③至⑤
ri
(t
+
∆t)
=
ri
(t)
+
vi
(t)
∆t
+
1 2
ai
(t)
∆t
2
vi
(t
+
∆t)
=
vi
(t)
+
1 2
[ai
(t)
+
ai
(t
+
∆t)]∆t
Verlet三种形式算法的比较：
Verlet
Leapfrog Velocity Verlet
④ 计算第n步的力 ⑤ 计算第n+1/2步的速度： ⑥ 计算第n+1步的位置： ⑦ 计算第n步的速度： ⑧ 重复④至⑦
vi
(t
+
1 2
∆t)
=
v
i
(t
-
1 2
∆t)
+
ai
(t)
∆t
ri
(t
+
∆t)
=
ri
(t)
+
vi
(t
+
1 2
∆t)
∆t
vi (t)
=
vi (t
+
1 2
∆t) + 2
vi (t
-
1 2
无因次量：
r* = r /σ
ρ* = ρσ 3
E* = E / ε
f * = fσ / ε
t* =
t
(mσ 2 / ε )1/ 2
v* =
v
(ε / m)1/ 2
P* = Pσ 3 / ε
T* = kT / ε
MD模拟中几个热力学量的计算：
对于由N个单原子组成的系统：
① 动能和温度：
∑ E K
=
1 2
四、预测－校正(Predictor-Corrector)格式算法：
1. 预测(Predictor)阶段：其基本思想是Taylor展开，
r
p
(t
+
∆t)
=
r(t)
+
v(t)∆t
+
1 2
a(t)∆t
2
+
1 6
b(t)∆t 3
+L
v
p
(t
+
∆t)
=
v(t)
+
a(t)∆t
+
1 2
b(t)∆t
2
+L
a p (t + ∆t) = a(t) + b(t)∆t + L
r0 = r(t)
r1
=
dr(t) dt
∆t
r2
=
d 2r(t) dt 2
1 2
∆t 2
r3
=
d
3
r
(t
)
dt 3
1 6
∆t 3
r
p
(t
+
∆t)
=
r(t)
+
v(t)∆t
+
1 2
a(t)∆t 2
+
1 6
b(t )∆t 3
+L
v
p
(t
+
∆t)
=
v(t)
+
a(t)∆t
+
1 2
b(t)∆t 2
+
L
a p (t + ∆t) = a(t) + b(t)∆t +L
RX(I) = RXNEWI RY(I) = RYNEWI RZ(I) = RZNEWI 100 CONTINUE
Verlet算法的优缺点：
优点： 1、精确，误差O(Δ4) 2、每次积分只计算一次力 3、时间可逆
缺点： 1、速度有较大误差O(Δ2) 2、轨迹与速度无关，无法与热浴耦联
二、蛙跳(Leap-frog)算法：半步算法
+
∆t)
−
rp 2
(t
+
∆t)
Values
c0
c1
c2
3
0
1
1
4
1/6
5/6
1
5
19/120
3/4
1
6
3/20 251/360
1
c3
c4
c5
1/3 1/2 11/18
1/12 1/6
1/60
② 二阶运动方程之二：
&r& = f (r,r&)
∆r
=
rc 2
(t
+
∆t)
−
rp 2
(t
+
∆t)
Values
： 粒子速度
：
vi
(t)
=
ri
(t
+
∆t) − ri 2∆t
(t
-
∆t)
粒子加速度：
a i
(t)
=
F i
(t)
mi
开始运动时需要r(t-∆t)：
r(−∆t) = r(0) − vi (0) ∆t
缺点：Verlet算法处理速度非常笨拙
Verlet算法的表述：
算法启动
① 规定初始位置
② 规定初始速度
③ 扰动初始位置：
分子动力学introduction
分子动力学简史
•1957年：基于刚球势的分子動力学法（Alder and Wainwright） •1964年：利用Lennard-Jone势函数法对液态氩性质的模拟（Rahman） •1971年：模拟具有分子团簇行为的水的性质（Rahman and Stillinger） •1977年：约束动力学方法（Rychaert, Ciccotti & Berendsen; van Gunsteren） •1980年：恒压条件下的动力学方法（Andersen法、Parrinello-Rahman法） •1983年：非平衡态动力学方法（Gillan and Dixon） •1984年: 恒温条件下的动力学方法(Berendsen et al.) •1984年：恒温条件下的动力学方法（Nosé-Hoover法） •1985年：第一原理分子動力学法（→Car-Parrinello法） •1991年：巨正则系综的分子动力学方法（Cagin and Pettit）
c0
c1
c2
3
0
1
1
4
1/6
5/6
1
5
19/90
3/4
1
6
3/16 251/360
1
c3
c4
c5
1/3 1/2 11/18
1/12 1/6
1/60
五、积分时间步长∆t的选择：
太长的时间步长会造成分子间的激烈碰撞，体系数据溢 出;
太短的时间步长会降低模拟过程搜索相空间的能力
室温下， ∆ t ≈ 1 fs (femtosecond 10-15s)， 温度越高，∆ t 应该减小
r(−∆t) = r(0) − vi (0) ∆t
④ 计算第n步的力 ⑤ 计算第n+1步的位置：
计算第n步的速度： 重复④至⑥
ri (t + ∆t) = 2ri (t) − ri (t - ∆t) + ai (t) ∆t2
vi
(t)
=
ri
(t
+
∆t) − ri 2∆t
(t
-
∆t)
Verlet算法程序：
VXI = ( RXNEWI – RXOLD(I) ) / DT2 VYI = ( RYNEWI – RYOLD(I) ) / DT2 VZI = ( RZNEWI – RZOLD(I) ) / DT2
RXOLD(I) = RX(I) RYOLD(I) = RY(I) RZOLD(I) = RZ(I)
粒子位置的Taylor展开式：
r i
(t
+
∆t)
=
r i
(t)
+
vi (t)
∆t
+
1 2
ai (t)
∆t 2
+
1 6
bi (t)
∆t 3
+
r i
(t
−
∆t)
=
r i
(t)
−
v i
(t)
∆t
+
1 2
a i
(t)
∆t 2
−
1 6
b i
(t)
∆t 3
粒子位置 ri (t + ∆t) = 2ri (t) − ri (t - ∆t) + ai (t) ∆t2
(t)
+
ai
(t
+
∆t)]∆t
等价于
vi
(t
+
1 2
∆t)
=
vi
(t)
+
1 2
ai
(t) ∆t
优点：速度计算更加准确
ri
(t
+
∆t)
=
ri
(t)
+
vi
(t
+
1 2
∆t)
∆t
vi
(t
+
∆t)
=
vi
(t
+
1 2
∆t)
+
1 2
ai
(t
+
∆t)∆t
Velocity Verlet算法的表述：
算法启动
① 规定初始位置
1. 首先利用当前时刻的加速度，计算半个时间步长后的速度：
vi
(t
+
1 2
∆t)
=
vi
(t
-
1 2
∆t)
+
ai
(t)
∆t
开始运动时需要v(-∆t/2)：
2. 计算下一步长时刻的位置：
ri
(t
+
∆t)
=
ri
(t)
+
vi
(t
+
1 2
∆t)
∆t
v(−∆t/2) = v(0) − ai (0) ∆t/2
3.
b p (t + ∆t) = b(t) +L
r0p (t + ∆t) 1 1 1 1r0 (t)
r1p r2p r3p
(t (t (t
+ + +
∆t ∆t ∆t
) ) )
=
0 0 0
1 0 0
2 1 0
3 r1(t)
13
r2 r3
(t (t
) )
校正阶段运动方程的变换：
r0c (t + ∆t) r0p (t + ∆t) c0
r1c r2c r3c
(t (t (t
+ + +
∆t) ∆t) ∆t)
=
r1p r2p r3p
(t (t (t
+ ∆t
C0，
+ ∆t
+ ∆t
) ) )
+
c 1
c 2
c3
∆r
∆r C0， C1， C2， C3的值以及
的形式：
取决于运动方程的阶数。
① 一阶运动方程：
r& = f (r)
∆t)
Leap-frog算法的优缺点：
优点： 1、提高精确度 2、轨迹与速度有关，可与热浴耦联
缺点： 1、速度近似 2、比Verlet算子多花时间
三、Velocity Verlet算法：
ri
(t
+
∆t)
=
ri
(t)
+
vi
(t)
∆t
+
1 2
ai
(t)
∆t 2
vi
(t
+
∆t)
=
vi
(t)Hale Waihona Puke Baidu
+
1 2
[ai
微正则系综分子动力学(NVE MD)
它是分子动力学方法的最基本系综 具有确定的粒子数N，能量E和体积V 算法： ① 规定初始位置和初始速度
② 对运动方程积分若干步 ③ 计算势能和动能 ④ 若能量不等于所需要的值，对速度进行标度 ⑤ 重复②至④，直到系统平衡
微正则系综(NVE)MD模拟算法的流程图：
Do 100 I = 1, N RXNEWI = 2.0 * RX(I) − RXOLD(I) + DTSQ * AX(I) RYNEWI = 2.0 * RY(I) − RYOLD(I) + DTSQ * AY(I) RZNEWI = 2.0 * RZ(I) − RZOLD(I) + DTSQ * AZ(I)
给定每个分子的初始位置ri(0)和速度vi(0) 计算每个分子的受力Fi和加速度ai
解运动方程并求出每个分子运动一个时间步 长后到达的位置所具有的速度
移动所有分子到新的位置并具有当前时刻的 速度
统计系统的热力学性质及其它物理量
No
Yes
统计性质不变？
打印结果，结束
微正则系综MD模拟程序F3讲解（LJ, NVE）：
b p (t + ∆t) = b(t) + L
2. 校正(Corrector)阶段：
根据新的原子位置rp，可以计算获得校正后的ac(t+∆t),定义预测误差:
∆a(t + ∆t) = ac (t + ∆t) − a p (t + ∆t)
利用此预测误差，对预测出的位置、速度、加速度等量进行校正：
rc
(t
+
∆t)
=
r
p
(t)
+
c 0
∆a(t
+
∆t
)
vc (t + ∆t) = v p (t) + c1∆a(t + ∆t)
ac (t + ∆t) = a p (t) + c2∆a(t + ∆t)
bc (t + ∆t) = b p (t) + c3∆a(t + ∆t)
预测阶段运动方程的变换：
定义一组矢量:
N i=1
mv2 i
=
3 2
Nk T B
采用对比量：
v* =
v
(ε / m)1/ 2
T* = kT / ε
E* = E / ε
∑ E* K
=
1 2
N i =1
v*2 i
=
3 2
NT *
② 势能：
∑ ∑ ∫ U
=Uc
+ Ulrc
=
N −1 i =1
N
Uc (rij ) +
j =i +1
Nρ
2
∞
dr
u(r)4πr
2
rc
对于LJ流体：
U
(r) = 4ε σ
12 − σ
6
c
r r
U lrc
=
8 3
πNρσ
3ε
1 3
σ
r c
9
−
σ
课程讲解内容：经典分子动力学 (Classical Molecular Dynamics)
粒子的运动取决于经典力学 (牛顿定律（F=ma）
分子动力学方法基础：
原理： 计算一组分子的相空间轨道，其中每个分子各自服从 牛顿运动定律：
∑ ∑ ∑ H
=
1 2
N i =1
pi2 mi
N −1 N
+
U (rij )
通过求解所有粒子的运动方程，分子动力学方法可以用于模 拟与原子运动路径相关的基本过程。
在分子动力学中，粒子的运动行为是通过经典的Newton运动 方程所描述。
分子动力学方法是确定性方法，一旦初始构型和速度确定 了，分子随时间所产生的运动轨迹也就确定了。
分子动力学的算法：有限差分方法
一、Verlet算法
i=1 j =i+1
pi
=
m i
dri dt
=
m i
vi
∑ ∑ ∑ ∑ dpi
dt
=
N −1 i =1
N
F(rij
j =i +1
)
=
N −1
−
i =1
N j =i +1
∂U (rij ∂rij
)
r = r (0) 初始条件： i t=0 i
dr
i
dt
t=0 = vi (0)
分子动力学方法特征：
分子动力学是在原子、分子水平上求解多体问题的重要的计 算机模拟方法，可以预测纳米尺度上的材料动力学特性。