粒子位置的Taylor展开式：
ri
(t
t)
ri
(t)
vi
(t)
t
1 2
ai
(t)
t 2
1 6
bi
(t)
t 3
+
ri
(t
t)
ri
(t)
vi
(t)
t
1 2
ai
(t)
t 2
1 6
bi
(t)
t 3
粒子位置 ： 粒子速度 ：
ri (t t) 2ri (t) ri (t - t) ai (t) t2
vi
利用此预测误差，对预测出的位置、速度、加速度等量进行校正：
rc (t t) r p (t) c0a(t t) vc (t t) v p (t) c1a(t t) ac (t t) a p (t) c2a(t t) bc (t t) b p (t) c3a(t t)
预测阶段运动方程的变换：
pi
mi
dri dt
mi vi
dpi
dt
N 1 i 1
N
F(rij
j i 1
)
N 1
i 1
N U (rij j i1 rij
)
r r (0) 初始条件： i t0 i
dri dt
t0 vi (0)
分子动力学方法特征：
分子动力学是在原子、分子水平上求解多体问题的重要的计 算机模拟方法，可以预测纳米尺度上的材料动力学特性。
④ 计算第n步的力 ⑤ 计算第n+1步的位置： 计算第n步的速度： 重复④至⑥
ri (t t) 2ri (t) ri (t - t) ai (t) t2
vi
(t)
ri
(t
t) ri 2t
(t
-
t)
Verlet算法程序：
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)
a p (t t) a(t) b(t)t
b p (t t) b(t)
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
) ) )
c1 c2 c3
r
r C0， C1， C2， C3的值以及
的形式：
取决于运动方程的阶数。
① 一阶运动方程：
9
rc
3
采用对比量：r* r / * 3 E* E /
U
* c
(r
*
)
4
1 r*
12
1 r*
6
U* lrc
8 3
πN
*
1 3
1 rc*
9
1 rc*
3
③ 内能：
内能由势能和动能组成：
E U EKLeabharlann 采用对比量：E*U*
E
* K
④ 压力： P
k BT
Wc V
定义一组矢量:
r0 r(t)
dr(t) r1 dt t
r2
d 2r(t) dt 2
1 t 2 2
r3
d 3r(t ) dt 3
1 t 3 6
r p (t t) r(t) v(t)t 1 a(t)t 2 1 b(t)t3
2
6
v p (t t) v(t) a(t)t 1 b(t)t 2 2
vi
(t
t)
vi
(t
1 2
t)
1 2
ai
(t
t)t
Velocity Verlet算法的表述：
算法启动
① 规定初始位置
② 规定初始速度 ③ 计算第n+1步的位置： ④ 计算第n+1步的力 ⑤ 计算第n+1步的速度： ⑥ 重复③至⑤
ri
(t
t)
ri
(t)
vi
(t)
t
1 2
ai
(t)
t
2
vi
(t
t)
vi
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)
1 2
t) 2
vi
(t
-
1 2
t)
vr
v
t-Δt/2 t t+Δt/2 t+Δt t+3Δt/2 t+2Δt
Leap-frog算法的表述：
算法启动
① 规定初始位置
② 规定初始速度
③ 扰动初始速度：
v(t/2) v(0) ai (0) t/2
④ 计算第n步的力 ⑤ 计算第n+1/2步的速度： ⑥ 计算第n+1步的位置： ⑦ 计算第n步的速度： ⑧ 重复④至⑦
v
i
(t
1 2
t)
v
i
(t
-
1 2
t)
ai
(t)
t
1 ri (t t) ri (t) vi (t 2 t) t
vi (t)
vi (t
1 2
t) 2
vi (t
-
1 2
t)
Leap-frog算法的优缺点：
优点： 1、提高精确度 2、轨迹与速度有关，可与热浴耦联
缺点： 1、速度近似 2、比Verlet算子多花时间
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 应该减小
(t)
1 2
[ai
(t)
ai
(t
t)]t
Verlet三种形式算法的比较：
Verlet
Leapfrog Velocity Verlet
四、预测－校正(Predictor-Corrector)格式算法：
1. 预测(Predictor)阶段：其基本思想是Taylor展开，
r p (t t) r(t) v(t)t 1 a(t)t 2 1 b(t)t 3
(t)
ri
(t
t) ri 2t
(t
-
t)
粒子加速度：
ai
(t)
Fi (t) mi
开始运动时需要r(t-Δt)：
r(t) r(0) vi (0) t
缺点：Verlet算法处理速度非常笨拙
Verlet算法的表述：
算法启动
① 规定初始位置
② 规定初始速度
③ 扰动初始位置：
r(t) r(0) vi (0) t
给定每个分子的初始位置ri(0)和速度vi(0)
计算每个分子的受力Fi和加速度ai
解运动方程并求出每个分子运动一个时间步 长后到达的位置所具有的速度
移动所有分子到新的位置并具有当前时刻的 速度
统计系统的热力学性质及其它物理量
No
Yes
统计性质不变？
打印结果，结束
微正则系综MD模拟程序F3讲解（LJ, NVE）：
RX(I) = RXNEWI RY(I) = RYNEWI RZ(I) = RZNEWI 100 CONTINUE
Verlet算法的优缺点：
优点： 1、精确，误差O(Δ4) 2、每次积分只计算一次力 3、时间可逆
缺点： 1、速度有较大误差O(Δ2) 2、轨迹与速度无关，无法与热浴耦联
二、蛙跳(Leap-frog)算法：半步算法
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 r2c (t t) r2p (t t)
Values
c0
c1
c2
3
0
1
1
4
1/6
5/6
1
通过求解所有粒子的运动方程，分子动力学方法可以用于模 拟与原子运动路径相关的基本过程。
在分子动力学中，粒子的运动行为是通过经典的Newton运动 方程所描述。
分子动力学方法是确定性方法，一旦初始构型和速度确定了， 分子随时间所产生的运动轨迹也就确定了。
分子动力学的算法：有限差分方法