第三章_经典分子动力学方法

第三章经典分子动力学方法3.1引言•分子动力学(MolecularDynamics,简写为MD)方法是确定性模拟方法，这方法是按该体系内部的内禀动力学规律来计算确定位形的转变。•首先需要建立一组分子的运动方程，然后通过直接对系统中的每一个原子/分子运动方程进行数值求解，得到每个时刻每个原子/分子的坐标与动量(速度)，即在相空间的运动轨迹，再利用统计方法得到多体系统的静态和动态特性，从而得到系统的宏观性质。•在MD方法的处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的，在这个微观的物理体系中，每个原子/分子都各自服从经典的牛顿力学定律。•MD方法是实现玻尔兹曼的统计力学途径，可以处理与时间有关的过程，因而可以处理非平衡态问题，但是该方法的计算机程序较复杂，计算量大，占内存也多。MD方法的发展史•MD方法是20世纪50年代后期由B.JAlder和T.E.Wainwright创造发展的。他们在1957年利用MD方法，发现了早在1939年根据统计力学预言的“刚性球组成的集合系统会发生由其液相到结晶相的相转变”。•20世纪70年代，产生了刚性体系的动力学方法被应用于水和氮等分子性溶液体系的处理，取得了成功。1972年，A.W.Less和S.F.Edwards等人发展了该方法，并扩展到了存在速度梯度(即处于非平衡状态)的系统。•之后，此方法被M.J.Gillan等人推广到了具有温度梯度的非平衡系统，从而构造并形成了非平衡MD方法体系。MD方法的发展史•MD方法真正作为材料科学领域的一个重要研究方法，开始于恒压MD方法(1980)和恒温MD方法(1984)的建立及在应用方面的成功。•1985年人们又提出了将电子论和分子动力学方法有机统一起来的所谓Car－Parrinello方法，即第一性原理MD方法。它不仅可以处理半导体和金属的问题，同时还可应用于处理有机物和化学反应。•1991年有人进一步提出了巨正则系综MD方法，从而又可适用于吸附问题的处理等，该方法还在进一步发展之中。•分子动力学方法的主要发展可见表3.1。年代创立者创造内容工作（MD分类名称）1957B.JAlder&T.E.Wainwright刚性球MD方法1963A.Rahman质点系MD方法1971A.Rahman&F.H.Stillinger刚性系统MD方法1972A.W.Lees&S.F.Edwards平衡系统MD方法（存在速度梯度）1977J.P.Rychaertetal.约束系统MD方法1980Andersen,Parrinello&Rahman恒压MD方法1983N.J.Gillan&M.Dixon非平衡MD方法（存在温度梯度）1984S.Nosé恒温MD方法1985R.Car&M.Parrinello第一性原理MD方法(Car-Parrinello方法)1991Cagin&Pettitt巨正则系统MD方法表3.1MD方法的里程碑工作3.2MD方法计算初步•在计算机出现以前，作为根据原子间相互作用力等微观信息了解多原子或分子团的结构、性质的方法，所采用的是基于统计理论的数学解析法。然而，原子间相互作用力稍微复杂一些，不用说求解统计理论严格方程解，就是进行数值求解也是一件很困难的事。•MD方法就是数值求解多体系统的确定性运动方程，并根据对所求结果进行统计处理，决定粒子的轨迹，从而给出物性预测和微观结构信息的一种模拟方法。内能比热容运动方程温度、压力相互作用原子位置坐标3维结构原子坐标、速度原子运动热力学性质动力学性质光学性质（输出信息）（二次信息）扩散系数粘滞系数电导率红外吸收图3-1MD方法信息输入输出信息方框图•MD这种方法并不严格。因此，必须根据情况，检验改变所模拟的基本单元尺寸所得结果是否会改变，直到所得结果不随基本单元尺寸变化而变化。通常这样的处理在很多情况下是有效的。对于基本单元内的原子、分子运动方程，使用什么样的形式合适，要具体问题具体分析。若是考虑具有确定的粒子数N，体积V和能量E的NEV系综(称为微正则系综，Micro-CanonicalEnsemble)，则其运动方程可以表达成式(3-2-1)所示的普通牛顿方程的形式22iiidrmFdt(3-2-1)式中mi为所考察的原子质量，ri为原子的位置坐标，Fi为作用在原子上的原子相互作用的合力，它由下式给出1NiiijjF(3-2-2)其中，Φij是原子和原子j之间的势函数(有时亦称为力场)例如，由氩原子等组成的稀薄气体，其势函数可采用Lennard－Jones势，1264ijijijijrr(3-2-3)式中，r是原子间距，是结合强度参数，是表示原子半径的参数。ijij在Δt时间内，对系统内的所有粒子解运动方程111()()()()ininininnnrtrtvtvtttt1nmaxttYesNo启动计算设定坐标、速度初始值计算作用在原子上的力Fi计算要求的物理量，将数据写入轨迹文件ttmax输出计算结果，并结束计算对(3-2-1)可用数值积分法求解，其数据处理流程图见图3-2图3-2MD数据处理流程图MD方法NEV能量恒定NTV恒温NHP恒压NTP恒温恒压VT巨正则系VL恒化学势弹性力学（原子分子）质点力学（原子分子）刚体力学（分子）约束力（分子和晶体）系）动力学模型目标系统团簇块体材料表面界面MBE/CVDBerrele法Green法多重时间刻度法数值积分法边界条件到目前为止已经确立的MD方法的主要技术体系①统计系综•系综是一个巨大的系统，由组成、性质、尺寸和形状完全一样的全同体系构成数目极多的系统的集合。不同的系综，MD方法的基本方程有所不同。•目前除微正则系综（NEV系综）外，已完成了正则系综（NTV系综），等温等压系综（NTP系综），等压等焓系综（NHP系综），巨正则系综（μVL系综），恒定化学势系综（μVT系综）等五个系综的MD方法的基本方程的确立。•现在已经能够处理许多体系，例如：–孤立宏观团簇的模拟(用NEV或NTV系综)–固体的结构相变，玻璃转变，晶化过程的模拟(用NTP系综)–固体(晶体)表面的原子、分子吸附现象的模拟(用μVT系综)②力学条件已建立了弹性力学、质点力学、刚体力学、约束力学等不同力学条件下的四种体系的MD方法。•弹性力学方法是将所考察的原子分子看作刚性球来处理，建立完全弹性碰撞方程，借以求解出原子、分子的运动规律。这种处理可以在液晶的模拟中使用。•质点力学模型是将原子、分子作为质点处理，粒子间的相互作用力采用坐标的连续函数。这种力学体系的应用对象非常多，可以用于处理陶瓷、金属、半导体等无机化合物材料以及有机高分子、生物大分子等几乎所有的材料。②力学条件•刚体力学方法是把分子作为刚体处理，建立对于刚体的欧拉方程和对于质点系的牛顿方程，联立求解所研究对象的运动问题，以前主要是在处理像水和四氯化碳那样的低分子量体系，现已用于研究晶体的相变。•约束力学是冻结粒子体系的一部分自由度，进而求解因此而生成约束条件下的质点力学运动方程。对于有机分子来讲，因为键长、键角的振动变化对体系的结构影响不大，将这些自由度冻结是合适的。另外也有采用固定晶体结构的考虑。③边界条件问题•在处理原子、分子的聚合体问题时，就MD方法而言，能处理的原子(分子)数目要受到计算机运行速度和能力的限制。目前报导的最好水平是能处理109量级的原子数目。这与现实物质含有1023个原子或分子的差距还很大，导致模拟系统原子数少于真实系统的所谓“尺寸效应”的问题。•为了减小“尺寸效应”而又不至于使计算工作量过大，对于平衡态MD模拟采用“周期性边界条件”。周期性边界条件图3-4假设现实物质中的一部分原子(通常为102－105个)，被取出配置在所谓基本单元的箱中，由于基本单元周围的原子、分子变成表面，从而不同于本来要处理的体状态。为了防止这种情况，在基本单元周围配置其复制品(图3-4)。周期性边界条件•对于在基本单元周围边界的原子、分子所受到的作用力，不仅有来自基本单元之内部的原子或分子的作用力的贡献，还要考虑来自其近邻复制单元的原子或分子的作用力的贡献。这样，把来自距粒子某一距离（截断距离）内的j粒子的贡献（这里与粒子j是否处在基本单元内或复制单元内无关）求和。给出力的方法称为最小镜像距离法（MinimumImageDistance）。周期性边界条件•如果在基本单元内的原子的位置为,周期边界条件会产生该原子的镜像。它的位置在这儿a、b、c是基本单元的三个边长，l，m，n是整数，取值范围为-∞到+∞。•在基本单元内的原子不仅与在本单元内的其他原子有相互作用，还与在相邻单元内的镜像原子有作用。imageiirrlambnc(3-2-4)基本单元大小的选择•基本单元的大小必须大于2Rcut（Rcut是相互作用势的截断距离）或Rcut1/2基本单元的大小。这保证了任何原子只与原子的一个镜像有相互作用，不与自己的镜像作用。这个条件称为“minimumimagecriterion”•在我们所研究的体系内的任何结构特性的特征尺寸或任何重要的效应的特征长度必须小于基本单元的大小。•为了检验不同基本单元大小是否会引入“人为效应”，必须用不同的基本单元尺寸做计算，若结果能收敛，则尺寸选择是合适的。边界条件•边界条件的问题比较复杂，因所考虑的具体情况而定，概括地说，作为处理原子、分子团簇的边界条件可进行推广。例如，适用于块体状态的三维周期边界条件，可扩展处理表面重构和表面吸附、单层膜相变的二维周期边界条件，以及处理异质晶体界面的边界系统等，同时还开发了模拟非平衡的方法。诸如针对分子束外延（MBE）的边界条件，人们进行了很多研究。•以二维边界条件为例。只在x-y平面配置基本单元的复制品，使用周期性边界条件，同时在z方向不赋予周期性边界条件，固定两端的数层原子。由于采用了这样的人工边界条件，所以在z轴方向上的原子层数要有适当的数目，一般考虑的标准线度为4－5纳米。采用二维边界条件，就可以用MD在原子尺度上研究界面的结构。④数值积分法•MD方法的基本方程是线性或非线性的二阶常微分方程。对此人们研究了许多求解方法。例如，贝鲁勒（Berreele）法，阿达姆斯（Adams）法，龙格-库塔（Runge-Kuta）法，追赶法等。•最近，M.Tuckerman和B.J.Berne提出了所谓多重时间宽度计算方法。即在质量相差很大的多原子体系，力学常数有很大不同的体系以及含有长、短程力的体系等情况下，必须结合具有最小振动周期的自由度选择积分的时间标度（）.这方法缩短了计算用的时间，有望成为MD方法中有前途的数值积分法。•通常～10-15s=1fstt3.3MD模拟的数据分析MD模拟给出的直接数据是坐标和速度，对它们进行加工，可得到各种所要的信息。这些信息大致可分为3类：•·所模拟系统的平衡态性质•·系统处于亚稳态时的结构与性质•·系统在远离平衡态时的动力学过程平衡态•一个物理量A随时间变化而趋于平衡态的值A0的过程可用下式表示其中是驰豫时间。–对于小的，我们只需等待系统达到平衡态后，再开始收集系统在平衡态时的参量。–对很大的，要达平衡态，需要模拟很长的时间，这时MD方法是不合适的。–对中等大小的，如果我们不能直接得到，但可以估计A0。•在许多情况下，我们并不需要达到平衡。如果我们的目的是研究非平衡的过程。)/exp()(0tBAtA(3-3-1)在MD模拟中可观察量的测量•可把可观察量表示为一个位置和速度的函数对上式做时间平均，得•我们更关心的是整个系统的平均效应。(,)iiAtfrtvt(3-3-2)(3-3-3)stepsstartNNtstartstepstANNA)(1为什么要计算能量？•验证体系是否总能量守恒•若有能量从动能转移到势能，表示体系发生了相变•若在E～T曲线上有跳跃，表示有一个一阶相变。能量计算•每个原子的势能为（3－3－4）•每个原子的动能为（3－3－5