分子动力学简介

1、分子动力学方法概述 首先，分子动力学方法是一种确定性方法。 分子动力学方法是按照该体系内部的动力学规律来确定位置和速度随时间的演化的一种模拟方法。 跟踪系统中每个粒子的个体运动。 根据统计物理规律，给出微观量（座标、速度）与宏观可观测量（温度、压力、比热、弹性模量等）的关系，从而研究物质和材料的性能。分子动力学方法的分类 按所遵循的运动规律或研究的对象分类： 经典MD平衡态MD微正则系综（NVE）正则系综（NVT）等温等压系综（NPT）等焓等压系综（NPH）非平衡态MD量子MD分子动力学模拟的基本步骤i 设定模拟的物理模型，建立运动微分方程ii 给定初始条件，数值求解运动微分方程iii 趋于平

2、衡的计算过程iv 宏观物理量的计算具体步骤 设定模拟物理模型 首要步骤对经典MD而言，每个分子遵循牛顿运动规律 mai(t)=Fi(r), i=1,2,N Fi 是体系内N个经典粒子作用在第i个粒子上力的总和。 ijjiiNjijirrrFrF1|)(|)(具体步骤分子之间相互作用势： 势函数的选取是模型设定步骤中的关键 MD中通常是基于经验数据和规律来选取势函数1、硬球势（最先被应用的势函数） ， 如果r V(r)= 0， 如果r 其中为相互作用势为0时的两分子之间的距离。具体步骤2、Lennard-Jones位势（最简单和最常用） 其中，-是位势的最小值（即势阱最低处的V(r)值）， 可确

3、定能量的单位；为V(r)等于零时的r值， 可以确定长度的单位。周期性边界条件 A(x)=A(x+nL) A为任意的可观测量)()(4)(612rrrV具体步骤MD中考虑粒子间相互作用时，通常采用最小像力的约定： 在无穷重复的分子动力学基本元胞中，只考虑每个粒子只同它所在的基本元胞内的另外N-1个粒子或其最邻近的影像粒子相互作用； 既有， 式中rij为位于ri处的粒子i同ri处的粒子j之间的距离。 |)min(|nLrrrjiij具体步骤对于短程力采用截断半径法 相互作用势的确定十分关键，然而势函数的确定往往相当困难，并不能精确知道。具体步骤 给定初始条件 一般来说MD中初始条件是不知道的，而且

4、初始条件往往不重要； 常用的初始条件选择： 1、令初始位置处在离散网格的格点上，初始速度则从波尔兹曼分布随机抽样确定； 2、令初始位置随机地偏离网格的格点，初始速度取为零； 3、初始位置与速度都随机确定：初始位置随机偏离网格格点，初始速度有波尔兹曼分布随机抽样确定。具体步骤 数值求解运动方程和趋衡计算 模拟计算首先从数值求解运动方程出发，方法有：一步法： 采用有限差分法，将微分方程变为有限差分方程， 取差分计算的时间步长为h,运用到展开的一阶泰勒公式：具体步骤整理后得到微分方程的欧拉算法：这是一组递推公式。Verlet方法（二步法的一种） 将泰勒公式二次展开 具体步骤蛙跳Leap-Frog算法

5、： 具体步骤具体步骤另外还有速度Verlet算法 可以同时给出位置、速度与加速度，并且不牺牲精度。 预测校正算法时间步长h的选取 过大误差太大，过小平衡时间太长，同时增加累加误差。具体步骤 宏观物理量的计算 最终目的 通过相对轨迹求平均得到。平衡态分子动力学模拟 微正则系综（NVE） 系统的能量、粒子数和体积均被处理成恒定量。经典MD平衡态MD微正则系综（NVE）正则系综（NVT）等温等压系综（NPT）等焓等压系综（NPH）非平衡态MD量子MD微正则系综 前提条件 对微正则系综进行MD模拟时首先要确定系统中粒子相互作用的模型，故假定一个孤立的多粒子体系，粒子间的相互作用位势是球对称的，则系统的

6、哈密顿量为： （rij第i个粒子和第j个粒子 之间距离） 时间不显示地出现在哈密顿量中，系统受一下约束： 系统能量E守恒； 系统的粒子数恒定； 系统的体积不变； 整个系统并未运动，系统的总动量P恒等于零。 jiijiirupmH)(1212微正则系综从牛顿运动方程出发： i=1,2,N 采用Verlet方法（中心差分格式）求上式： i=1,2,N 可知，由t和t-h的位置及t时刻的受力能求出t+h时刻粒子的位置； 现令， 上式可以改写为： , i=1,2,N 即为计算粒子坐标随时间演化的递推公式； 一组初始空间位置 随后粒子的所有位置； 粒子运动速度： 由n+1步位置算出n步的速度，可见动能的

7、计算比势能落后一步。jiijiirFmdttrd)(1)(22jiiiiimhtFhtrtrhtr/)()()(2)(2nhtn)()(ninitrr)()(ninitFFjininininimhFrrr/22)()1()()1(hrrVninini2/ )()1()1()(微正则系综 具体模拟步骤： 给定初始空间位置： ， 在n步时计算粒子所受的力： 计算粒子第n+1步的位置 计算第n步的速度， 返回步骤2，开始下一次模拟计算。 改进： 把N个粒子的初始位置放置在网格的格点上，然后加以扰动，给出的初始条件是粒子的空间位置和运动速度，可用 如下公式计算粒子位置：)0(ir)1( nir)()(

8、ninitFF)1(irmhFrrrnininini/22)()1()()1(hrrVninini2/ )()1()1()(mhFhvrriiii2/22)0(0)0()1 (微正则系综 上述Verlet算法可以作变形改进，使数值计算的稳定性得到改善，即为Verlet算法的速度形式： 令， 则前面提及的差分方程组可写为： Verlet算法的速度形式模拟步骤可写为： 给定粒子的初始空间位置 ，及初始速度 计算第n+1步所有粒子所处的位置hrrzninini/ )()()1()()1()1()(nininihzrr)()1()(nininiFmhzz)1(ir)1(iv)1( nir)(2)()(

9、)1(21ninininiFhmhvrr微正则系综 计算n+1步所有粒子的速度 返回到步骤3，进行下一步的模拟计算。 这样的优点是成功的得到了同一时间步长上的空间位置和速度，另外，数值计算的稳定性也加强了。 总述 一般来说，一个给定的系统并不知道其精确的初始条件，需要给出一个合理的初始条件，然后在模拟过程中对能量进行增减调节。具体，先算出若干步的动能和势能，如果不符合给定的恒定量，则乘以一个标度因子，再回到第一步。)1( nivmFFhvvnininini2/ )(1)()1(正则系综 概述 很多情况下，所研究系统的温度恒定，如同将系统置于热浴中，系统能量可能有涨落，这就是正则系综（NVT），

10、粒子数N，体积V，温度T和总动量P为守恒量。 正则系综MD模拟的约束为： 平衡化过程中，用以上约束 进行标度调整。 模拟步骤 给定粒子初始空间位置 ； 给定粒子初始速度 ；iimvA常量221)1(ir)1(iv正则系综 计算出n+1时间步时所有粒子所处的空间位置 计算出n+1时间步时所有粒子速度， 计算相应的动能 和标度因子 对速度进行标度调整，作为下一次的计算值， )1( nir)(2)()()1(21ninininiFhmhvrrmFFhvvnininini2/ )(1)()1(inikvm2)1()(21E2/12)1()()43(inivmkTN)1( niv)1( niv正则系综

11、返回步骤3，开始n+2时间步的模拟计算。 按上述步骤，对时间逐步递进模拟计算，知道系统趋至平衡，完成正则系综的分子动力学模拟。均方位移与扩散系数关系式推导 扩散定律 Fick第一定律：单位时间内通过垂直扩散方向的单位截面的物质量（扩散通量）与该物质在该面积处的浓度梯度成正比， Fick第二定律：考虑到时间对浓度的影响，同时近似将扩 散系数D看做常数， 均方位移法（MSD） 固体扩散中，扩散原子并不是沿着直线扩散，二是呈折线 无规则地在晶体的间隙或者空位之间跳动。xCDJ22CxCDt均方位移与扩散系数关系式推导 由于原子跳动方向有正负，会出现抵消，因此只考虑纯位移就失去了意义因此提出了均方位移。 均方位移是总位矢的平方，可代表原子经过n次跳动之后 的净位移。 上式中r代表原子每次跳动的距离。 对于两个相邻平行晶面 面密度是n1，n2；晶面距离为d，原子跳动的频率是f，原子跳动概率的平方P 22Rnrn均方位移与扩散系数关系式推导 那么在t内，面1跳向面二和面二跳向面1的原子数为：两式相减并利用扩散通量J的定义有：把面密度n1，n2改为体密度C1，C2 tPfnN121tPfnN212)(211221nnPfJtNNPf