1、第三章第三章 经典分子动力学方法经典分子动力学方法第1页，共39页。3.1 引言引言分子动力学(Molecular Dynamics,简写为MD)方法是确定性模拟方法，这方法是按该体系内部的内禀动力学规律来计算确定位形的转变。首先需要建立一组分子的运动方程，然后通过直接对系统中的每一个原子/分子运动方程进行数值求解，得到每个时刻每个原子/分子的坐标与动量(速度)，即在相空间的运动轨迹，再利用统计方法得到多体系统的静态和动态特性，从而得到系统的宏观性质。在MD方法的处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的，在这个微观的物理体系中，每个原子/分子都各自服从经典的牛顿力学定律。M
2、D方法是实现玻尔兹曼的统计力学途径，可以处理与时间有关的过程，因而可以处理非平衡态问题，但是该方法的计算机程序较复杂，计算量大，占内存也多。第2页，共39页。MD方法的发展史方法的发展史 MD方法是20世纪50年代后期由B.J Alder和T.E.Wainwright创造发展的。他们在1957年利用MD方法，发现了早在1939年根据统计力学预言的“刚性球组成的集合系统会发生由其液相到结晶相的相转变”。20世纪70年代，产生了刚性体系的动力学方法被应用于水和氮等分子性溶液体系的处理，取得了成功。1972年，A.W.Less和S.F.Edwards等人发展了该方法，并扩展到了存在速度梯度(即处于非
3、平衡状态)的系统。之后，此方法被M.J.Gillan等人推广到了具有温度梯度的非平衡系统，从而构造并形成了非平衡MD方法体系。第3页，共39页。MD方法的发展史方法的发展史 MD方法真正作为材料科学领域的一个重要研究方法，开始于恒压MD方法(1980)和恒温MD方法(1984)的建立及在应用方面的成功。1985年人们又提出了将电子论和分子动力学方法有机统一起来的所谓CarParrinello方法，即第一性原理MD方法。它不仅可以处理半导体和金属的问题，同时还可应用于处理有机物和化学反应。1991年有人进一步提出了巨正则系综MD方法，从而又可适用于吸附问题的处理等，该方法还在进一步发展之中。分子
4、动力学方法的主要发展可见表3.1。第4页，共39页。年代创立者创造内容工作（MD分类名称）1957B.J Alder&T.E.Wainwright刚性球MD方法1963A.Rahman质点系MD方法1971Rahman&F.H.Stillinger刚性系统MD方法1972W.Lees&S.F.Edwards平衡系统MD方法（存在速度梯度）1977J.P.Rychaert et al.约束系统MD方法1980Andersen,Parrinello&Rahman恒压恒压MD方法方法1983N.J.Gillan&M.Dixon非平衡MD方法（存在温度梯度）1984 S.Nos恒温恒温MD方法方法19
5、85R.Car&M.Parrinello第一性原理第一性原理MD方法方法(Car-Parrinello方法方法)1991Cagin&Pettitt巨正则系统MD方法表3.1 MD方法的里程碑工作方法的里程碑工作 第5页，共39页。3.2 MD方法计算初步方法计算初步 在计算机出现以前，作为根据原子间相互作用力等微观信息了解多原子或分子团的结构、性质的方法，所采用的是基于统计理论的数学解析法。然而，原子间相互作用力稍微复杂一些，不用说求解统计理论严格方程解，就是进行数值求解也是一件很困难的事。MD方法就是数值求解多体系统的确定性运动方程，并根据对所求结果进行统计处理，决定粒子的轨迹，从而给出物性
6、预测和微观结构信息的一种模拟方法。第6页，共39页。内能比热容运动方程温度、压力相互作用原子位置坐标3维结构原子坐标、速度原子运动热力学性质动力学性质光学性质（输出信息）（二次信息）扩散系数粘滞系数电导率红外吸收图图3-1 MD方法信息输入输出信息方框图方法信息输入输出信息方框图第7页，共39页。MD这种方法并不严格。因此，必须根据情况，检验改变所模拟的基本单元尺寸所得结果是否会改变，直到所得结果不随基本单元尺寸变化而变化。通常这样的处理在很多情况下是有效的。第8页，共39页。对于基本单元内的原子、分子运动方程，使用什么样的形式合适，要具体问题具体分析。若是考虑具有确定的粒子数N，体积V和能量
7、E的NEV系综(称为微正则系综，Micro-Canonical Ensemble)，则其运动方程可以表达成式(3-2-1)所示的普通牛顿方程的形式22iiid rmFdt (3-2-1)式中mi为所考察的原子质量，ri为原子的位置坐标，Fi为作用在原子上的原子相互作用的合力，它由下式给出1NiiijjF (3-2-2)其中，ij是原子和原子j之间的势函数(有时亦称为力场)第9页，共39页。例如，由氩原子等组成的稀薄气体，其势函数可采用LennardJones势，1264ijijijijrr (3-2-3)式中，r 是原子间距，是结合强度参数，是表示原子半径的参数。ijij第10页，共39页。在
8、t时间内，对系统内的所有粒子解运动方程111()()()()ininininnnr tr tvtvtttt 1nmaxttYesNo启动计算设定坐标、速度初始值计算作用在原子上的力Fi计算要求的物理量，将数据写入轨迹文件 ttmax输出计算结果，并结束计算对(3-2-1)可用数值积分法求解，其数据处理流程图见图3-2 图图3-2 MD数据处理流程图数据处理流程图第11页，共39页。MD方法NEV能量恒定NTV恒温NHP恒压NTP恒温恒压恒温恒压VT巨正则系VL恒化学势弹性力学（原子分子）质点力学质点力学（原子分子）刚体力学（分子）约束力（分子和晶体）系）动力学模型目标系统团簇块体材料块体材料表
9、面界面MBE/CVDBerrele法Green法多重时间刻度法数值积分法边界条件到目前为止已经确立的MD方法的方法的主要技术体系 第12页，共39页。统计系综统计系综系综是一个巨大的系统，由组成、性质、尺寸和形状完全一样的全同体系构成数目极多的系统的集合。不同的系综，MD方法的基本方程有所不同。目前除微正则系综（NEV系综）外，已完成了正则系综（NTV系综），等温等压系综（NTP系综），等压等焓系综（NHP系综），巨正则系综（VL系综），恒定化学势系综（VT系综）等五个系综的MD方法的基本方程的确立。现在已经能够处理许多体系，例如：孤立宏观团簇的模拟(用NEV或NTV系综)固体的结构相变，玻璃
10、转变，晶化过程的模拟(用NTP系综)固体(晶体)表面的原子、分子吸附现象的模拟(用VT系综)第13页，共39页。力学条件力学条件 已建立了弹性力学、质点力学、刚体力学、约束力学等不同力学条件下的四种体系的MD方法。弹性力学方法是将所考察的原子分子看作刚性球来处理，建立完全弹性碰撞方程，借以求解出原子、分子的运动规律。这种处理可以在液晶的模拟中使用。质点力学模型质点力学模型是将原子、分子作为质点处理，粒子间的相互作用力采用坐标的连续函数。这种力学体系的应用对象非常多，可以用于处理陶瓷、金属、半导体等无机化合物材料以及有机高分子、生物大分子等几乎所有的材料。第14页，共39页。力学条件力学条件 刚
11、体力学方法是把分子作为刚体处理，建立对于刚体的欧拉方程和对于质点系的牛顿方程，联立求解所研究对象的运动问题，以前主要是在处理像水和四氯化碳那样的低分子量体系，现已用于研究晶体的相变。约束力学是冻结粒子体系的一部分自由度，进而求解因此而生成约束条件下的质点力学运动方程。对于有机分子来讲，因为键长、键角的振动变化对体系的结构影响不大，将这些自由度冻结是合适的。另外也有采用固定晶体结构的考虑。第15页，共39页。边界条件问题边界条件问题 在处理原子、分子的聚合体问题时，就MD方法而言，能处理的原子(分子)数目要受到计算机运行速度和能力的限制。目前报导的最好水平是能处理109量级的原子数目。这与现实物
12、质含有1023个原子或分子的差距还很大，导致模拟系统原子数少于真实系统的所谓“尺寸效应”的问题。为了减小“尺寸效应”而又不至于使计算工作量过大，对于平衡态MD模拟采用“周期性边界条件”。第16页，共39页。周期性边界条件 图3-4假设现实物质中的一部分原子(通常为102105个)，被取出配置在所谓基本单元的箱中，由于基本单元周围的原子、分子变成表面，从而不同于本来要处理的体状态。为了防止这种情况，在基本单元周围配置其复制品(图3-4)。第17页，共39页。周期性边界条件 对于在基本单元周围边界的原子、分子所受到的作用力，不仅有来自基本单元之内部的原子或分子的作用力的贡献，还要考虑来自其近邻复制
13、单元的原子或分子的作用力的贡献。这样，把来自距粒子某一距离（截断距离）内的j粒子的贡献（这里与粒子j是否处在基本单元内或复制单元内无关）求和。给出力的方法称为最小镜像距离法（Minimum Image Distance）。第18页，共39页。周期性边界条件如果在基本单元内的原子的位置为,周期边界条件会产生该原子的镜像。它的位置在 这儿a、b、c是基本单元的三个边长，l，m，n是整数，取值范围为-到+。在基本单元内的原子不仅与在本单元内的其他原子有相互作用，还与在相邻单元内的镜像原子有作用。imageiirrlambnc (3-2-4)第19页，共39页。基本单元大小的选择基本单元大小的选择 基
14、本单元的大小必须大于2Rcut（Rcut是相互作用势的截断距离）或Rcut1012 m2/s时适用。2211NiiiMSDr tr tr oN (3-3-10)(3-3-9)第31页，共39页。GreenKubo formula 扩散系数D也可用速度自相关函数通过下面的GreenKubo 公式表示 （3311）GreenKubo公式和Einstein 关系是等价的。101133NiiiDv tv odtv tvo dtN第32页，共39页。Partial diffusion 在多种原子组成的体系中（例如AB2），可以对某一种原子求和，得到Partial diffusion，从而得到原子A或原子
15、B的扩散系数DA、DB。第33页，共39页。非均匀体系中的扩散非均匀体系中的扩散 对一非均匀体系，例如一块材料，有体和晶界部分。可分别对某一类原子做平均，如对晶界原子，或对体中的原子，得到不同部分原子不同的扩散。第34页，共39页。与温度有关的扩散与温度有关的扩散 在一个系统里，若原子的扩散是通过跳跃进行的，可将扩散用Arrhenius关系 来表达，式中是prefactor，是原子迁移率或激发能(activation energy)。对（3312）两边求对数有0expaBEDDk T(3-3-12)TkEDDBa1lnln0(3-3-13)第35页，共39页。与温度有关的扩散与温度有关的扩散
16、作lnD1/T图(图3.8)，从图中直线的斜率可以得到Ea，从截距可以得到D0。同理也可以得到表面原子、晶界原子的扩散系数(图3.9)。图3.8图3.9TkEDDBa1lnln0第36页，共39页。MD模拟材料的溶解模拟材料的溶解 固态到液态的相变（溶解）的过程发生在几ps时间，在MD模拟中可通过以下手段知道是否发生了固-液相变：观察不同时间段系统的原子位置情况（snapshots）能量从动能向势能转移 总能量和势能随温度的变化有跳跃 扩散系数有变化 径向分布函数有变化 在做有固定体积的模拟时，压力有突变第37页，共39页。MD模拟材料的溶解模拟材料的溶解 溶点对应于固-液共存状态，这时在computational cell 中应含有固态、液态及固-液界面。由图310可以得到溶点温度。图310第38页，共39页。体膨胀系数可以由VT曲线（图311）得到 PTVV)(1图311对不同的压强P做模拟，可以得到（T，P）。第39页，共39页。
