分子动力学简介

ri
( n)
hv
( n) i
1 2 (n) h Fi 2m
微正则系综
③计算n+1步所有粒子的速度 vi( n 1) vi( n 1) vi( n ) h( Fi n1 Fi n ) / 2m ④返回到步骤3，进行下一步的模拟计算。 这样的优点是成功的得到了同一时间步长上的空间位置和 速度，另外，数值计算的稳定性也加强了。
2 r2 d 2 P Rn 2 Rn 2 Dt D 2 d P r t r2 是一个与晶格属性及温度有关的量，与扩散系数D d 2P 无关，我们记作K，则：
2 1 Rn D K t
固体扩散中，扩散原子并不是沿着直线扩散，二是呈折线 无规则地在晶体的间隙或者空位之间跳动。
均方位移与扩散系数关系式推导
由于原子跳动方向有正负，会出现抵消，因此只考虑纯位 移就失去了意义因此提出了均方位移。
R 2 nr 2 n
均方位移是总位矢的平方，可代表原子经过n次跳动之后 的净位移。 上式中r代表原子每次跳动的距离。 •
微正则系综
• 从牛顿运动方程出发：
d 2 ri (t ) 1 Fi (rij ) 2 dt m i j
i=1,2,…,N
采用Verlet方法（中心差分格式）求上式： ri (t h) 2ri (t ) ri (t h) Fi (t )h 2 / m
i j
i=1,2,…,N

6
具体步骤
MD中考虑粒子间相互作用时，通常采用最小像力 的约定： 在无穷重复的分子动力学基本元胞中，只考 虑每个粒子只同它所在的基本元胞内的另外N-1个 粒子或其最邻近的影像粒子相互作用； 既有， rij min(| ri rj nL |) 式中rij为位于ri处的粒子i同ri处的粒子j之间的距离。
分子之间相互作用势：
势函数的选取是模型设定步骤中的关键 MD中通常是基于经验数据和规律来选取势函数 1、硬球势（最先被应用的势函数） ∞， 如果r<σ V(r)= 0， 如果r≥ σ 其中σ为相互作用势为0时的两分子之间的距离。
具体步骤
2、Lennard-Jones位势（最简单和最常用）
V (r ) 4 [( )12 ( )] r r 其中，-ε是位势的最小值（即势阱最低处的V(r) 值）， ε可确定能量的单位；σ为V(r)等于零时的r 值， σ可以确定长度的单位。 周期性边界条件 A(x)=A(x+nL) A为任意的可观测量
对于两个相邻平行晶面
面密度是n1，n2；晶面距离为d，原子跳动的频率是f，原 子跳动概率的平方P
均方位移与扩散系数关系式推导
那么在△t内，面1跳向面二和面二跳向面1的原子数为：
N12 n1 Pft
两式相减并利用扩散通量J的定义有：
N 21 n2 Pft
N12 N 21 J Pf (n1 n2 ) t
• 模拟步骤
(1) ①给定粒子初始空间位置 ri ；
②给定粒子初始速度 vi(1) ；
正则系综
ri ( n 1) ③计算出n+1时间步时所有粒子所处的空间位置
ri
( n 1)
ri
( n)
hv
( n) i
1 2 (n) h Fi 2m
④计算出n+1时间步时所有粒子速度，
vi( n 1) vi( n ) h( Fi n1 Fi n ) / 2m ⑤计算相应的动能 1 ( n 1)
ri( n ) ri( n1) hzi( n 1)
zi( n ) zi( n 1) h (n) Fi m
Verlet算法的速度形式模拟步骤可写为： ri(1) ，及初始速度 vi(1) ①给定粒子的初始空间位置 ri ( n 1) ②计算第n+1步所有粒子所处的位置
ri
( n 1)
ri(1) ，
②在n步时计算粒子所受的力： Fi ( n ) Fi (tn ) ③计算粒子第n+1步的位置 ri ( n 1)
ri( n1) 2ri( n ) ri( n1) Fi ( n ) h 2 / m
④计算第n步的速度， ( n ) (r ( n 1) r ( n 1) ) / 2h Vi i i ⑤返回步骤2，开始下一次模拟计算。
• 改进：
把N个粒子的初始位置放置在网格的格点上，然后加以扰 动，给出的初始条件是粒子的空间位置和运动速度，可用 如下公式计算粒子位置： ri(1) 2ri( 0) hvi0 Fi ( 0) h 2 / 2m
微正则系综
• 上述Verlet算法可以作变形改进，使数值计算的稳定性得 到改善，即为Verlet算法的速度形式： 令， zi( n ) (ri( n1) ri( n ) ) / h 则前面提及的差分方程组可写为：
具体步骤
• 数值求解运动方程和趋衡计算
模拟计算首先从数值求解运动方程出发，方法有： 一步法： 采用有限差分法，将微分方程变为有限差分方程， 取差分计算的时间步长为h,运用到展开的一阶泰 勒公式：
具体步骤
整理后得到微分方程的欧拉算法：
这是一组递推公式。 Verlet方法（二步法的一种） 将泰勒公式二次展开
平衡态分子动力学模拟
平衡态MD 经典MD 非平衡态 MD 微正则系综 （NVE） 正则系综 （NVT）
等温等压系 综（NPT）
等焓等压系 综（NPH）
量子MD
• 微正则系综（NVE）
系统的能量、粒子数和体积均被处理成恒定量。
微正则系综
• 前提条件
对微正则系综进行MD模拟时首先要确定系统中粒子相互 作用的模型，故假定一个孤立的多粒子体系，粒子间的相 互作用位势是球对称的，则系统的哈密顿量为： （rij第i个粒子和第j个粒子 1 1 2 H pi u (rij ) 之间距离） 2 i m i j 时间不显示地出现在哈密顿量中，系统受一下约束： ①系统能量E守恒； ②系统的粒子数恒定； ③系统的体积不变； ④整个系统并未运动，系统的总动量P恒等于零。
• 总述
一般来说，一个给定的系统并不知道其精确的初始条件， 需要给出一个合理的初始条件，然后在模拟过程中对能量 进行增减调节。具体，先算出若干步的动能和势能，如果 不符合给定的恒定量，则乘以一个标度因子，再回到第一 步。
正则系综
• 概述
很多情况下，所研究系统的温度恒定，如同将系统置于热 浴中，系统能量可能有涨落，这就是正则系综（NVT）， 粒子数N，体积V，温度T和总动量P为守恒量。 1 A mvi2 常量 正则系综MD模拟的约束为： 2 i 平衡化过程中，用以上约束 进行标度调整。
可知，由t和t-h的位置及t时刻的受力能求出t+h时刻粒子的位置； 现令， t nh F ( n ) F (t ) r ( n ) r (t )
nቤተ መጻሕፍቲ ባይዱ
i i n
i i n
上式可以改写为：ri
( n 1)
2ri( n) ri( n1) Fi ( n) h 2 / m , i=1,2,…,N
分子动力学方法
概述
• 首先，分子动力学方法是一种确定性方法。 • 分子动力学方法是按照该体系内部的动力 学规律来确定位置和速度随时间的演化的 一种模拟方法。 • 跟踪系统中每个粒子的个体运动。 • 根据统计物理规律，给出微观量（座标、 速度）与宏观可观测量（温度、压力、比 热、弹性模量等）的关系，从而研究物质 和材料的性能。
i j
即为计算粒子坐标随时间演化的递推公式； 一组初始空间位置 随后粒子的所有位置； 粒子运动速度： ( n ) ( n 1) ( n 1)
Vi
(ri
ri
) / 2h
由n+1步位置算出n步的速度，可见动能的计算比势能落后一步。
微正则系综
• 具体模拟步骤：
①给定初始空间位置： ri
( 0)
把面密度n1，n2改为体密度C1，C2
J (n1 n2 ) Pf (C1d C2 d ) Pf C 2 d Pf x
均方位移与扩散系数关系式推导
将上式与Fick第一定律比较有：D d 2 Pf d和P是由晶体结构类型决定的，f除了与晶体结构类型有 关外，和温度也有极大的关系。 2 上式与 Rn ftr 2 相结合有：
均方位移与扩散系数关系式推导
• 扩散定律
Fick第一定律：单位时间内通过垂直扩散方向的单位截面 的物质量（扩散通量）与该物质在该面积处的浓度梯度成 C 正比，
J D x
Fick第二定律：考虑到时间对浓度的影响，同时近似将扩 散系数D看做常数， C 2C D 2 t x
• 均方位移法（MSD）
Ek 2
m(vi
i
)2
和标度因子β
[
(3N 4)kT 1/ 2 ] m(vi( n 1) ) 2
i
⑥对速度进行标度调整，作为下一次的计算值，
vi( n 1)
vi( n 1)
正则系综
⑦返回步骤3，开始n+2时间步的模拟计算。 按上述步骤，对时间逐步递进模拟计算，知道系统趋至平 衡，完成正则系综的分子动力学模拟。
具体步骤
对于短程力采用截断半径法
相互作用势的确定十分关键，然而势函数的确定 往往相当困难，并不能精确知道。
具体步骤
• 给定初始条件
一般来说MD中初始条件是不知道的，而且初始 条件往往不重要； 常用的初始条件选择：
1、令初始位置处在离散网格的格点上，初始速度则从波 尔兹曼分布随机抽样确定； 2、令初始位置随机地偏离网格的格点，初始速度取为零； 3、初始位置与速度都随机确定：初始位置随机偏离网格 格点，初始速度有波尔兹曼分布随机抽样确定。
ii
iii
iv
具体步骤
• 设定模拟物理模型 首要步骤
对经典MD而言，每个分子遵循牛顿运动规律 mai(t)=Fi(r), i=1,2,…,N Fi 是体系内N个经典粒子作用在第i个粒子上力的 总和。
Fi (r ) Fi (| ri r j |) r ij
j 1 i j
N
^
具体步骤
具体步骤
蛙跳Leap-Frog算法：
具体步骤
具体步骤
另外还有速度Verlet算法 可以同时给出位置、速度与加速度，并且不牺牲 精度。 预测—校正算法 时间步长h的选取 过大误差太大，过小平衡时间太长，同时增加累 加误差。
具体步骤
• 宏观物理量的计算 最终目的 通过相对轨迹求平均得到。
分子动力学方法的分类
• 按所遵循的运动规律或研究的对象分类：
平衡态MD
经典MD 非平衡态 MD 微正则系综 （NVE）
正则系综 （NVT）
等温等压系 综（NPT） 等焓等压系 综（NPH）
量子MD
分子动力学模拟的基本步骤
i
• 设定模拟的物理模型，建立运动微分方程 • 给定初始条件，数值求解运动微分方程 • 趋于平衡的计算过程 • 宏观物理量的计算