2.4 Equations of motion
分子动力学模拟
为了在计算机上解运动方程，必须为微分方程建立一个 有限差分格式，从差分方程中再导出位置和速度的递推关系 式。这些算法是一步一步执行的，先算t 时刻的位置和速度， 然后在此基础上计算t+1时刻的位置和速度。
微分方程最为直接的离散化格式来自泰勒展开: r(th)r(t)n i 1 1hi!ir(i)(t)Rn
1.5
1
间间
0.5
rij 6 2
0
-0.5
-1
0.8
1
1.2 1.4 1.6 1.8 间间
2
2.2 2.4 2.6
对势能的最大贡献来自于粒子的近邻区域，位势截断
常用的方法是球形截断，截断半径一般取2.5σ或3.6 σ，对
截断距离之外分子间相互作用能按平均密度近似的方法进
行校正。
分子动力学模拟
The disk processed after the simulation is finished. It contains at least all the positions and velocities of all particles. This information is sufficient to calculate all the properties of the system. However, it is more economical to calculate properties during the simulation and store them in the than reading the calculating them afterwards.
➢二、分子动力学方法
分子动力学模拟
2.1 Newtonian mechanics
In the MD method, every new distribution is derived from the previous one by using the interactions between the particles. The interactions depends on the position of the particles.
分子动力学模拟
2.3 Calculations of force, velocity, position
The initial distribution of the Molecular dynamics simulation is generated in a random distribution.
粒子数密度： 温度： 能量：
* 3
T*kBT/
E*E/
1 m3
TK
EJ
压强：
p* p3 /
p N/m2
时间： 力：
t*t(48/m2)12 t s
f*f/
f N
➢四、参量的计算
分子动力学模拟
A simulation run produces the raw data in form of a very large disk file, The the complete state of the system at each step.
1、有限差分方法-预测校正法
rp(tt)r(t)tv(t)t2a(t)/2t3b(t)/6 vp(tt)v(t)ta(t)t2b(t)/2 ap(tt)a(t)tb(t) bp(tt)b(t)
2、有限差分方法-Verlet算法
分子动力学模拟
①、Verlet算法的一般形式
fma 为了用数值方法求解微分方程， i
有效两体势 VV1(ri)V2eff(rij) V (rij )
i
i ji
它包含多体效应，可很好地反映系统粒子间的相互作用。
分子动力学模拟 下面仅对简单系统的相互作用模型给予简介
1、Lennard-Joans势
Lennard-Jones 间间间间 2
1.5
1
间间
0.5
rij 6 2
0
-0.5
-1
①、简单立方晶格
②、体心立方晶格 ③、面心立方晶格
c b
a
c b
a
c c
分子动力学模拟
c c
c ca a
bc bca a
a a
b ba a
b b
b b
初始速度
分子动力学模拟
模拟时，各粒子的初始速度按麦克斯韦速度分布取样。
Maxwell’s distribution law of velocity
分子动力学模拟
本章主要内容
分子动力学模拟
➢一、系综理论
➢二、分子动力学方法
➢三、模拟细节
➢四、参量的计算 ➢五、液态水的MD模拟
➢六、误差分析 ➢七、分子动力学模拟方法的应用
➢一、系综理论
分子动力学模拟
分子动力学模拟（molecular dynamics simulation，简称MD） 方法首先是由Alder和Wainwright提出的，现已逐渐成为预测系统特性、 验证理论和改进模型的计算工具。
在计算势能时必须考虑相互作用的力程，对粒子数为 N的模拟系统，原则上任何两个粒子间都有相互作用，在 计算体系势能时须进行N（N-1）/2次运算。
为了提高计算效率，实际模拟 过程中通常在一个方便的力程上 将位势截断，用以减少计算势能 所消耗的时间。
Lennard-Jones 间间间间 2
分子动力学模拟
Each particle is also assigned an initial velocity vi
In simulation:
fx48 2(xi xj)[ ri(j)141 2( rij)8] fy48 2(yi yj)[ ri(j)141 2( rij)8] fz4 82(zi zj)[ ri(j)141 2( rij)8]
对于正则系综，体系趋衡是通过增减体系能量来实现的。如果 在进行足够多的运算步数以后，系统能持续给出确定的平均动能和 平均势能的数值，那么就认为平衡已经建立。
间间间间 间间间间
间1 间间间间间间间间间间
-2900
-3000
-3100
-3200
-3300
-3400
T=293K
-3500
0
1
2
3
4
5
6
7
8
ri(th)2ri(t)ri(th)m 1h2fi(t)
rin12rinrin1m 1h2fin
速度对轨道计算没有关系，但对 估算动能及速度相关函数（用来研 究粒子的输运特性）非常有用。速 度形式为
vin
rin1 rin1 2h
②、Verlet 蛙跳格式（leap-frog）
分子动力学模拟
The velocity varies so one has to choose a reasonable average value to be used. The velocity at the middle of the step ought to be a good compromise,
The force causes an acceleration
fi mai
分子动力学模拟
Which in turn modifies the initial velocity vi as
v i' v i a i t v i a ih
And modifies the initial position ri as
对一个由N个原子构成的简单系统，其势能项由下式给出
V V 1 ( r i) V 2 ( r i,r j) V 3 ( r i,r j,r k )
i
ij i
ij ik j i
式中右端第一项是外场（如电场、 磁场、声场等）对系统的作用；第二项 是两体势即系统中每两个粒子间的相互 作用；第三项是三体势，表示系统中每 三个粒子间的相互作用……
0.8
1
1.2 1.4 1.6 1.8 间间
2
2.2
2.4 2.6
2. 硬球势（Hard-Sphere）
VHS(r) 0
r r
分子动力学模拟
3. 软球势（Soft-Sphere）
VSS(r)() r
r
通常，v 是为整数的参数。
4. 方阱势（Square-Well）
VSW(r)
0
r 1 1 r2 r 2
61
61
61
5D 2
5C 2
5B 2
3 4
3 4
3 4
61
61
61
5E
2
2 5基本元胞
5A
2
3 4
3 4
3 4
61
61
61
5F
2
5G 2
5H 2
3 4
3 4
3 4
L
3.3 Calculation of interactions 分子动力学模拟
i , j 两个离子的相互作用势
体系总势能可表示为：
U u(rij) ij
指平衡状态下理想气体分子 速度分布的统计规律。在平衡状 态下，分布在任一速率区间内的 分子数与总分子数的比率为
f(v)dv 42m kT 3/2exp2 m k2T vv2dv
分子动力学模拟
3.2 Periodic boundary conditions
通常选取小数体系（几十个到数千个分子）作为研究 对象，但由于位于表面和内部的分子受力差别较大，将会 产生表面效应。
指定的体系温度，T*是每一时间步的瞬时动能温度。
Ekin12mi
vi2
3Nk*T 2