分子动力学简介

分子动力学方法概述•首先，分子动力学方法是一种确定性方法。•分子动力学方法是按照该体系内部的动力学规律来确定位置和速度随时间的演化的一种模拟方法。•跟踪系统中每个粒子的个体运动。•根据统计物理规律，给出微观量（座标、速度）与宏观可观测量（温度、压力、比热、弹性模量等）的关系，从而研究物质和材料的性能。分子动力学方法的分类•按所遵循的运动规律或研究的对象分类：经典MD平衡态MD微正则系综（NVE）正则系综（NVT）等温等压系综（NPT）等焓等压系综（NPH）非平衡态MD量子MD分子动力学模拟的基本步骤i•设定模拟的物理模型，建立运动微分方程ii•给定初始条件，数值求解运动微分方程iii•趋于平衡的计算过程iv•宏观物理量的计算具体步骤•设定模拟物理模型首要步骤对经典MD而言，每个分子遵循牛顿运动规律mai(t)=Fi(r),i=1,2,…,NFi是体系内N个经典粒子作用在第i个粒子上力的总和。ijjiiNjijirrrFrF^1|)(|)(具体步骤分子之间相互作用势：势函数的选取是模型设定步骤中的关键MD中通常是基于经验数据和规律来选取势函数1、硬球势（最先被应用的势函数）∞，如果rσV(r)=0，如果r≥σ其中σ为相互作用势为0时的两分子之间的距离。具体步骤2、Lennard-Jones位势（最简单和最常用）其中，-ε是位势的最小值（即势阱最低处的V(r)值），ε可确定能量的单位；σ为V(r)等于零时的r值，σ可以确定长度的单位。周期性边界条件A(x)=A(x+nL)A为任意的可观测量)]()[(4)(612rrrV具体步骤MD中考虑粒子间相互作用时，通常采用最小像力的约定：在无穷重复的分子动力学基本元胞中，只考虑每个粒子只同它所在的基本元胞内的另外N-1个粒子或其最邻近的影像粒子相互作用；既有，式中rij为位于ri处的粒子i同ri处的粒子j之间的距离。|)min(|nLrrrjiij具体步骤对于短程力采用截断半径法相互作用势的确定十分关键，然而势函数的确定往往相当困难，并不能精确知道。具体步骤•给定初始条件一般来说MD中初始条件是不知道的，而且初始条件往往不重要；常用的初始条件选择：1、令初始位置处在离散网格的格点上，初始速度则从波尔兹曼分布随机抽样确定；2、令初始位置随机地偏离网格的格点，初始速度取为零；3、初始位置与速度都随机确定：初始位置随机偏离网格格点，初始速度有波尔兹曼分布随机抽样确定。具体步骤•数值求解运动方程和趋衡计算模拟计算首先从数值求解运动方程出发，方法有：一步法：采用有限差分法，将微分方程变为有限差分方程，取差分计算的时间步长为h,运用到展开的一阶泰勒公式：具体步骤整理后得到微分方程的欧拉算法：这是一组递推公式。Verlet方法（二步法的一种）将泰勒公式二次展开具体步骤蛙跳Leap-Frog算法：具体步骤具体步骤另外还有速度Verlet算法可以同时给出位置、速度与加速度，并且不牺牲精度。预测—校正算法时间步长h的选取过大误差太大，过小平衡时间太长，同时增加累加误差。具体步骤•宏观物理量的计算最终目的通过相对轨迹求平均得到。平衡态分子动力学模拟•微正则系综（NVE）系统的能量、粒子数和体积均被处理成恒定量。经典MD平衡态MD微正则系综（NVE）正则系综（NVT）等温等压系综（NPT）等焓等压系综（NPH）非平衡态MD量子MD微正则系综•前提条件对微正则系综进行MD模拟时首先要确定系统中粒子相互作用的模型，故假定一个孤立的多粒子体系，粒子间的相互作用位势是球对称的，则系统的哈密顿量为：（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步的速度，可见动能的计算比势能落后一步。jiijiirFmdttrd)(1)(22jiiiiimhtFhtrtrhtr/)()()(2)(2nhtn)()(ninitrr)()(ninitFFjininininimhFrrr/22)()1()()1(hrrVninini2/)()1()1()(微正则系综•具体模拟步骤：①给定初始空间位置：，②在n步时计算粒子所受的力：③计算粒子第n+1步的位置④计算第n步的速度，⑤返回步骤2，开始下一次模拟计算。•改进：把N个粒子的初始位置放置在网格的格点上，然后加以扰动，给出的初始条件是粒子的空间位置和运动速度，可用如下公式计算粒子位置：)0(ir)1(nir)()(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)()()1(21ninininiFhmhvrr微正则系综③计算n+1步所有粒子的速度④返回到步骤3，进行下一步的模拟计算。这样的优点是成功的得到了同一时间步长上的空间位置和速度，另外，数值计算的稳定性也加强了。•总述一般来说，一个给定的系统并不知道其精确的初始条件，需要给出一个合理的初始条件，然后在模拟过程中对能量进行增减调节。具体，先算出若干步的动能和势能，如果不符合给定的恒定量，则乘以一个标度因子，再回到第一步。)1(nivmFFhvvnininini2/)(1)()1(正则系综•概述很多情况下，所研究系统的温度恒定，如同将系统置于热浴中，系统能量可能有涨落，这就是正则系综（NVT），粒子数N，体积V，温度T和总动量P为守恒量。正则系综MD模拟的约束为：平衡化过程中，用以上约束进行标度调整。•模拟步骤①给定粒子初始空间位置；②给定粒子初始速度；iimvA常量221)1(ir)1(iv正则系综③计算出n+1时间步时所有粒子所处的空间位置④计算出n+1时间步时所有粒子速度，⑤计算相应的动能和标度因子β⑥对速度进行标度调整，作为下一次的计算值，)1(nir)(2)()()1(21ninininiFhmhvrrmFFhvvnininini2/)(1)()1(inikvm2)1()(21E2/12)1(])()43([inivmkTN)1(niv)1(niv正则系综⑦返回步骤3，开始n+2时间步的模拟计算。按上述步骤，对时间逐步递进模拟计算，知道系统趋至平衡，完成正则系综的分子动力学模拟。均方位移与扩散系数关系式推导•扩散定律Fick第一定律：单位时间内通过垂直扩散方向的单位截面的物质量（扩散通量）与该物质在该面积处的浓度梯度成正比，Fick第二定律：考虑到时间对浓度的影响，同时近似将扩散系数D看做常数，•均方位移法（MSD）固体扩散中，扩散原子并不是沿着直线扩散，二是呈折线无规则地在晶体的间隙或者空位之间跳动。xCDJ22CxCDt均方位移与扩散系数关系式推导由于原子跳动方向有正负，会出现抵消，因此只考虑纯位移就失去了意义因此提出了均方位移。均方位移是总位矢的平方，可代表原子经过n次跳动之后的净位移。上式中r代表原子每次跳动的距离。•对于两个相邻平行晶面面密度是n1，n2；晶面距离为d，原子跳动的频率是f，原子跳动概率的平方P22Rnrn均方位移与扩散系数关系式推导那么在△t内，面1跳向面二和面二跳向面1的原子数为：两式相减并利用扩散通量J的定义有：把面密度n1，n2改为体密度C1，C2tPfnN121tPfnN212)(211221nnPfJtNNPfdxCPfdCdCPfnnJ22121)()(均方位移与扩散系数关系式推导将上式与Fick第一定律比较有：d和P是由晶体结构类型决定的，f除了与晶体结构类型有关外，和温度也有极大的关系。上式与相结合有：是一个与晶格属性及温度有关的量，与扩散系数D无关，我们记作K，则：PfdD222ftrRntRrPdDDtPdrRnn222222Pdr22tRKDn21