分子动力学方法

1/ 2
vin 1

（6） 返回到步骤（3），开始第n+2步的模拟计算。
3.4 MD在材料科学中的应用


3.4.1 高分子链动力学模拟
计明娟等用高温淬火分子动力学模拟方法研究了甲硫氨酸 一脑啡肽(Met-enkephalin)在真空中的构象性质，经聚类分 析和能量优化得到了十三个低能构象，并与吗啡进行了空 间拟合。
( n) ( n1) ( n1) vi ri ri / 2h


（5）返回到步骤（2），开始下一步的模拟计算。
3.3 平衡态分子动力学模拟
Verlet算法的速度形式： (1) r （1）给定初始空间位臵 i 。 (1) v （2）给定初始速度 i 。 （3）利用公式：
3.2 分子动力学模拟的基本步骤

1．模拟模型的设定 硬球势
, U r 0,
如果 r , 如果 r .
4 3 2 1
位势V(r)

Lennard-Jones型势
U r 4 r r


元胞 周期性边界条件,
最小像力约定 rij min ri r j n L


分子动力学模拟的最小像力约定示意图
3.2 分子动力学模拟的基本步骤
2.给定初始条件 给定粒子的初始位臵和速度的数值： （1）令初始位臵在差分网格格子上，初始速度从玻尔兹曼 分布随机抽样得到。 （2）令初始位臵随机地偏离差分网格格子，初始速度为零。 （3）令初始位臵随机地偏离差分网格格子，初始速度从玻 尔兹曼分布随机抽样得到。 3．趋于平衡 使系统达到平衡，模拟中需要一个趋衡过程。在这个过 程中，我们增加或从系统中移出能量，直到系统具有所 要求的能量。
(1) (2)
(3) (4) (5) (6)
3.1.1粒子运动方程的数值求解
其他求解算法： Leap-frog algorithm r(t+∆t) = r(t) + v(t+(1/2)∆t) ∆t v(t+(1/2)∆t) = v(t-(1/2)∆t) + a(t) ∆t Beeman’s algorithm r(t+∆t) = r(t) + v(t)∆t + (2/3)a(t)∆t2 – (1/6)a(t-∆t)∆t2 v(t+∆t) = v(t) + v(t)∆t + (1/3)a(t)∆t + (5/6)a(t)∆t– (1/6)a(t∆t)∆t
u P T rij PC 6 N i j rij
k B


长程改正项为
u 3 PC g ( r ) 4 r dr r 6 c r

2
3.2 分子动力学模拟的基本步骤

分子动力学模拟的实际步骤可以划分为四步：

设定模拟所采用的模型； 给定初始条件； 趋于平衡的计算过程； 宏观物理量的计算。
1 1 2 m i k B T 2 2
3.1.2 能量均分定理

由于系统的每个粒子有三个自由度，因此
3 E k Nk B T 2 假定位势在处被截断。系统内部的位形能量的轨道平均值 为 n 1 v U U n n0 vn0


其中
U v u(rijv )
i j
3.1.2 能量均分定理
3.3 平衡态分子动力学模拟
1. 微正则系综的分子动力学模拟
粒子数恒定、体积恒定、能量恒定、整个系统的总动量恒 等于零。 分子动力学模拟步骤如下（Verlet算法）： （1）给定初始空间位臵。 （2）计算在第n步时粒子所受的力。 （3）计算在第n +1步时所有粒子所处的空间位臵。 ( n1) (n) (n1) ( n) 2 ri 2ri ri Fi h / m （4）计算第n步的速度。


（5） 返回到步骤（3），开始第n+2步的模拟计算。
3.3 平衡态分子动力学模拟


Байду номын сангаас

一般对于能量确定的系统不可能给出精确的初始条件。可 以将初始位臵固定在格子的格点上，而初始速度则由波尔 兹曼分布得出。 以上述方法建立的系统一般不会具有所需要的能量，需要 在模拟过程中逐渐调节系统能量达到给定值。步骤为： (1)对运动方程求解若干步； (2)计算动能和势能； (3)若能量不等于所需的值，对速度进行标度； (4)从第一步开始重复，直至系统达到平衡为止。
12 6
力F(r)
0 -1 -2 -3 0.6
排斥力
0.8 1.0 1.2 1.4 1.6
吸引力
1.8 2.0 2.2 2.4 2.6
3.2 分子动力学模拟的基本步骤
根据经典物理学的规律，可以知道在系综模拟中的守恒 量。 微正则系综的模拟中能量、动量和角动量均为守恒量。 在此系综中他们分别表示为： 2 1 V r P pi M ri pi E m r i i i i i 2
3.3 平衡态分子动力学模拟
正则系综分子动力学的模拟具体步骤: (Verlet算法的速度形式)
（1）给定初始空间位臵， （2）给定初始速度， （3）利用公式：
( n1) ( n) ( n) ( n) 2 ri ri hvi Fi h / 2m
( n1 ) ( n ) ( n1 ) ( n ) vi vi h Fi Fi / 2m
3.1.2 热力学量的计算

在物理系统的计算机模拟中，系综平均必须用时间平均代 替，在通常的模拟中，粒子数N和体积V是固定的。给定 初始位臵rN(0)和 初始动量pN(0)后，一个MD算法将从运 动方程生成轨道(rN(t), pN(t))，轨道平均的定义为
A lim(t 't 0 )
t '
E lim(t 't ) 1 t ' E ( (t ))dt, 0 t0 k k t ' t' 1 U lim(t 't 0 ) U (r (t ))dt, t0 t '
1 Ek n n0
v n0
v E k
n
3.1.2 热力学量的计算
g (r ) V n( r ) N 4r 2 r
3.1.2 能量均分定理


其他的量也需要尾部改正，以压强的计算为例，这时维里 (virial)状态方程成立。 2 u k 3 P BT g ( r ) 4 r dr 0 6 r 至于势能的计算，可以把积分分成两项，一项是由相互作 用力程之内的贡献引起的，一项是对位势截断的改正项：

生成的动能路径是不连续的，必须在各个时间间断点上计 算动能的值以求平均:
1 Ek n n0

其中
v k
v n0
v E k
n
1 E m( i2 ) v i 2

从平均动能可以计算系统的温度，温度是一个重要的物理 量，需要加以监测，特别是在模拟的起始阶段。
3.1.2 能量均分定理
3.1 基本原理
1.
分子动力学是在原子、分子水平上求解多体问题
的重要的计算机模拟方法。
2.
通过求解所有粒子的运动方程，分子动力学方法
可以用于模拟与粒子运动路径相关的基本过程。
3.
在分子动力学中，粒子的运动行为是通过经典的
Newton运动方程所描述。
3.1.1粒子运动方程的数值求解
粒子体系的运动方程Lagrangian方程



在热力学极限下，一切系综都是等同的，并且可以应用能 量均分定理。 热力学极限是指粒子数（或体积）趋向无穷大时的极限。 一般宏观物体包含了1023个粒子，可以认为是满足热力学 极限的。相变只有在热力学极限下才会发生。 系统的哈密顿量为

则：
Pi 2 1 u rij 2 i mi i j
3.3 平衡态分子动力学模拟
速度标度因子:
* kT ( N 1 ) 16 m v 2 i i i
1/ 2
2.正则系综的分子动力学模拟 速度标度因子:
( 3 N 4 ) kT m vi2 i
1/ 2
1

t'
t0
dtA(r N (t ) p N (t );V (t ))

假定能量守恒，并且轨道在一切具有同一能量的相同休积 内经历相同的时间，则轨道平均等于微正则系综平均
A A
NVE
3.1.2 热力学量的计算


孤立系统的总能量是一个守恒量，沿着分子动力学模拟生 成的任何一条轨道，能量应保持不变，即 E E 。 孤立系统的动能和势能不是守恒量，它们的大小沿着生成 的轨道逐点变化
k
qk
k p
H qk
3.1.1粒子运动方程的数值求解

单原子的牛顿运动方程
pi ri mi
.
m i r i f i ri V
..
3.1.1粒子运动方程的数值求解

Verlet 算法 r(t+∆t) = r(t) + v(t)∆t + (1/2)a(t)∆t 2 r(t-∆t) = r(t) – v(t)∆t + (1/2)a(t)∆t2 将上面两式相加，得到： r(t+∆t) = 2r(t) – r(t- ∆t) + a(t)∆t2 v(t+∆t) = v(t) + a(t)∆t + (1/2)b(t)∆t2 a(t+∆t) = a(t) + b(t)∆t 将 (5) 式的b(t) 代入 (4) ，得到： v(t+∆t) = v(t) + (1/2)[a(t) + a(t+∆t)] ∆t
1.
Lagrangian函数的定义为
则运动的Lagrangian方程为
d L L 0 k qk dt q
3.1.1粒子运动方程的数值求解
粒子的运动方程-Lagrangian方程
L pk k q H pk
) k pk L(q, q H (p, q) q
分子动力学方法
3.1 基本原理
计算机模拟分类：
(1)随机模拟方法。 优点: 随机模拟方法计算的程序简单，占内存少，但 是该方法难于处理非平衡态的问题。
(2)分子动力学方法（Molecular Dynamics或简称 MD） 。 可以处理非平衡态问题。但是使用该方法的程序较 复杂，计算量大，占内存也多。 分子动力学方法利用牛顿古典力学来计算许多分子 在相空间中的轨迹。
3.1 基本原理


分子动力学(MD)方法的出发点是物理系统的确定的微 观描述（哈密顿描述方程、拉格朗日方程或牛顿运动 方程）。因此，分子动力学方法是用运动方程来计算 多体或少体系统的性质，结果得到的既有系统的静态 特性，也有动态特性。蒙特卡罗方法只能得到系统的 位形特性。 MD方法的具体做法是在计算机上求运动力程的数值解。 为此，通过适当的格式对方程进行近似，使之适于在 计算机上求数值解。其实质是计算一组分子的相空间 轨道，其中每个分子各自都服从经典运动定律。这里 的系统不仅是点粒子系统，也包括具有内部结构的粒 子组成的系统。
Ek 1 m vi( n 1) 2 i
计算在第n+1步时所有粒子所处的空间位臵， （4） 计算在第n+1步时所有粒子的速度：

2

动能和速度标度因子：


（5） 计算将速度乘以标度因子的值，并让该值作为下一次计算时， 第n+1步粒子的速度：
v
n 1 i
(3N 4)kT m(v ( n 1) ) 2 i i

由于位势被截断，总能量和势能含有误差，为了估计必须 作出改正。所有的内部位形能都加到截止距离为止，尾部 改正可以取
U c 2 u(r ) g (r )r 2 dr
rc


其中是g(r)对关联函数，它是在原点r=0处有一个粒子时,在 r周围的体积元dr内找到一个粒子的概率，令n(r)为离一个 给定粒子的距离在r和r+dr之间的平均粒子数，于是
( n1) ( n) ( n) ( n) 2 ri ri hvi Fi h / 2m
r 计算在第n+1步时所有粒子所处的空间位臵
( n 1) i
。
（4） 计算在第n+1步时所有粒子的速度： (n1) (n) (n1) ( n) vi vi h Fi Fi / 2m