《第6章 分子动力学方法.docx》由会员分享，可在线阅读，更多相关《第6章 分子动力学方法.docx（24页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第6章分子动力学方法第6章分子动力学方法经典分子动力学方法无疑是材料，尤其是大分子体系和大体系模拟有效的方法之一。分子动力学能够用于NPT，NVE，NVT等不同系综的计算，是一种基于牛顿力学确定论的热力学计算方法。与蒙特卡罗法相比在宏观性质计算上具有更高的准确度和有效性，能够广泛应用于物理，化学，生物，材料，医学等各个领域。本章在介绍分子动力学的基本概念的基础上，简单介绍了分子动力学的基本思想，势函数分类和基本方程。然后介绍了分子动力学的常用系综和典型的NPT，NVE，NVT系综基本方程。结合材料建模中的基本简化方法和技巧，阐述了边界条件和时间积分的数值处理技巧。最后，利用统计力学的基本概念给
2、出分子动力学的计算信息的解析方式。并且结合MaterialsExplore软件计算分析了CNT的几何构造稳定性。6.1引言分子动力学方法(MolecularDynamics,MD)方法是一种按该体系内部的内禀动力学规律来计算并确定位形的变化确实定性模拟方法。首先需要在给定的外界条件下建立一组粒子的运动方程，然后通过直接对系统中的一个个粒子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计力学方法得到多体系统的静态和动态特性，进而获得系统的宏观性质。能够看出，分子动力学方法中不存在任何随机因素，这个也是分子动力学方法和后文要提到的蒙特卡洛方法的区别之一。在分
3、子动力学方法的处理经过中，方程组的建立是通过对物理体系的微观数学描绘给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学定律(或者是拉格朗日方程)。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日函数来描绘，可以以直接用牛顿运动方程来描绘。确定性方法是实现玻尔兹曼的统计力学途径。这种方法能够处理与时间有关的经过，因此能够处理非平衡态问题。但是分子动力学方法的计算机程序相对蒙特卡罗较复杂，其计算成本较高。分子动力学方法发展历史改革经历了近60年，分子动力学方法是20世纪50年代后期由AlderBJ和WainwrightTE创造发展的。Alder和Wainwright在19
4、57年利用分子动力学模拟，发现了“刚性球组成的集合系统会发生由其液相到结晶相的相转变，后来人们称这种相变为Alder相变。其结果表明，不具有引力的粒子系统也具有凝聚态。到20世纪70年代，产生了刚性体系的动力学方法，成功地被应用于水和氮等分子性溶液体系的处理；1972年，LessAW和EdwardsSF等发展了该方法并扩展到了存在速度梯度，即处于非平衡状态的系统。之后，此方法被GinanMJ等推广到了具有温度梯度的非平衡系统，进而构造并构成了所谓的非平衡分子动力学方法体系。进人20世纪80年代之后，出现了在分子内部对一部分自由度施加约束条件的新的分子动力学方法，进而使分子动力学方法可适用于类似
5、蛋白质等生物大分子的解析与设计。分子动力学方法真正作为材料科学领域的一个重要研究方法，开场于由Andersen,Parrinello和Rahman等创立恒压分子动力学方法和Nse等完成恒温分子动力学方法的建立及在应用方面的成功。后来，针对势函数模型化比拟困难的半导体和金属等，1985年人们又提出了将基于密度泛函理论的电子论和分子动力学方法有机统一起来的所谓Car-Parrinello方法，亦即第一性原理分子动力学方法。这样，分子动力学的方法进一步得到发展和完善，它不仅能够处理半导体和金属的间题，同时还可应用于处理有机物和化学反响。关于Car-Parrinello分子动力学方法将在第7章重点学习
6、讨论。本章重点讨论经典分子动力学方法的基本原理和计算方法。6.2分子动力学的计算框架6.2.1基本思想分子动力学方法是沿用牛顿运动方程或拉格朗日方程、哈密顿方程，通过考察粒子的运动来研究多粒子系统的物理性质。在处理孤立粒子或原子团簇(不考虑与其他粒子的互相作用)时，单纯地对牛顿运动方程进行积分求解就行了；而在处理凝聚系统时，能够以为将所考虑的N个粒子放入具有一定体积的容器中，在这样组成的封闭系中，建立直角坐标(x,y,z)，每个粒子的位置就由三个坐标分量决定，通过求解3N个联立方程组就能够得到保守系的总能量。分子动力学的基本思想在于通过设定原子之间的互相作用(势函数)和相关的系综(亦即作用对象
7、和条件)，确定其基本的模拟范畴。在给定的牛顿方程、拉格朗日方程或者是哈密顿方程进行时间的迭代，在到达指定的力学收敛条件之后得到一个最终的位置坐标，然后通过相关的位形信息和运动的速度和加速度等信息通过热力统计学的方法，给出模拟体系的统计信息，主要包括复杂的光学，电学和一些其他的物理信息。其基本思想，总结如图6-1所示。一次信息二次信息动力学方程初始化条件图6-1分子动力学的一般思想6.2.2计算流程总体来讲，分子动力学的基本计算有如下四个步骤：1.模型的设定，也就是势函数的选取势函数的研究和物理系统中对物质的描绘研究息息相关。分子动力学的模拟最早使用是硬球势，即小于临界值时无穷大，大于等于临界值
8、时为零。常用的是Lennard-Jones势函数，还有嵌入原子EAM势函数，不同的物质状态描绘用不同的势函数。模型势函数一旦确定，就能够根据物理学规律求得模拟中的守恒量。关于势函数选取将在6.4节进一步讨论。2.给定初始条件运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如Verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。一般意义上讲系统的初始条件不可能知道，实际上也不需要准确选择代求系统的初始条件，由于当模拟时间足够长时，系统就会忘掉初始条件(对于无记忆的体系而言)。当然，合理的初始条件能够加快系统趋于平衡的
9、时间和步程，获得好的精度。常用的初始条件能够选择为，令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令3.趋于平衡计算在边界条件和初始条件给定后就能够解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。为使得系统平衡，模拟中设计一个趋衡经过，即在这个经过中，增加或者从系统中移出能量，直到持续给出确定的能量值，称此时系统已经到达平衡，到达平衡的时间称为弛豫时间。分子动力学中，积分步程的大小选择特别重要，这决定了模拟所需要的时间。为了减小误差，步长要小，但过
10、小系统模拟的弛豫时间就长。因而根据经历选择适当的步长。如，对一个具有几百个氩气分子的体系，采用Lennard-Jones势函数，发现取h为0.01量级，能够得到很好的结果。这里选择的h是没有量纲的，实际上这样选择的h对应的时间在10-14s的量级。假如模拟1000步，系统到达平衡，弛豫时间只要10-11s。4.宏观物理量的计算实际计算宏观的物理量往往是在模拟的最后阶段进行的，它是沿相空间轨迹求平均来计算得到的，根据各态历经假讲，时间平均代替系综平均。关于各态历经假讲将在第8章的蒙特卡洛方法中深化讨论。6.2.3经典分子动力学中近似处理实际的计算体系由于其外在条件的复杂性，往往需要对计算的系统进
11、行一定的近似处理，在经典分子动力学计算中，引进下面近似处理：?经典粒子互相作用，不考虑电子互相作用量子效应。?力的作用形式，由参数可调的互相作用势函数决定，并经过实验测定来验证。?模拟体系与实际体系相差较大，一般需要采用周期边界来扩展计算体系。时间平均是在有限时间内完成。【练习与考虑】6-1.查找文献，根据上述的分子动力学的处理流程图，编写实现分子动力学的简单程序，可参考DaanF和BerendS编著的（分子模拟）一书。6.3分子动力学的系综系综(ensemble)是指在一定的宏观条件下(约束条件)，大量性质和构造完全一样的、处于各种运动状态的、各自独立的系统集合，或称为统计系综。系综是用统计
12、方法描绘热力学系统的统计规律性时引入的一个基本概念；系综是统计理论的一种表述方式，系综理论使统计物理成为普遍的微观统计理论；系综并不是实际的物体，构成系综的系统才是实际物体。约束条件是由一组外加宏观参量来表示。在平衡统计力学范畴下，能够用来处理稳定系综。6.3.1常用系综分类1.正则系综(canonicalensemble)全称应为“宏观正则系综，简写为NVT，即表示具有确定的粒子数(N)、体积(V)、温度(T)。正则系综是分子动力学方法模拟处理的典型代表。假定N个粒子处在体积为V的盒子内，将其埋入温度恒为T的热浴中。此时，总能量(E)和系统压强(P)可能在某一平均值附近起伏变化。平衡体系代表
13、封闭系统，与大热源热接触平衡的恒温系统。正则系综的特征函数是亥姆霍兹自由能(,)FNVT。2.微正则系综(micro-canonicalensemble)简写为NVE，即表示具有确定的粒子数(N)、体积(V)、总能量(E)。微正则系综广泛被应用在分子动力学模拟中。假定N个粒子处在体积为V的盒子内，并固定总能量(E)。此时，系综的温度(T)和系统压强(P)可能在某一平均值附近起伏变化。平衡体系为孤立系统，与外界即无能量交换，也无粒子交换。微正则系综的特征函数是熵(,)SNVE。3.等温等压(constant-pressure,constant-temperature)简写为NPT，即表示具有确定
14、的粒子数(N)、压强(P)、温度(T)。其总能量(E)和系统体积(V)可能存在起伏。体系是可移动系统壁情况下的恒温热浴。特征函数是吉布斯自由能GNPT。(,)4.等压等焓(constant-pressure,constant-enthalpy)简写为NPH，即表示具有确定的粒子数(N)、压强(P)、焓(H)。由于HEPV=+，故在该系综下进行模拟时要保持压力与焓值为固定，其调节技术的实现也有一定的难度，这种系综在实际的分子动力学模拟中已经很少碰到了。5.巨正则系综(grandcanonicalensemble)简写为VT，即表示具有确定的粒体积(V)、温度(T)和化学势()。巨正则系综通常是蒙
15、特卡罗模拟的对象和手段。此时、系统能量(E)、压强(P)和粒子数(N)会在某一平均值附近有一个起伏。体系是一个开放系统，与大热源大粒子源热接触平衡而具有恒定的温度(T)。特征函数是马休(Massieu)函数J(，V，T)。对于不同的系综，由于其选择的运动方程的求解方法亦有不同。针对不同的计算体系详细的细节比拟复杂，下文将以NEV和NTP这两类简单的系综运用两类不同的条件方法下求解的不同方程做扼要的介绍。6.3.2NEV系综基本方程NEV系综是一类较为常见的系综。对于该系综，经典分子动力学方法一般采用差分近似法求解牛顿运动方程，并追踪系统的时间变化。考虑由N个粒子构成的系统，假设第i个粒子在时刻
16、()ixt位置，速度为()ivt遭到的作用力为()iFt，为求出该粒子(i)在时刻tt+?的位置，经常使用的方法是Verlet算法(对于该算法的详细细节在下文将有具体的推导)。根据系统的动能及统计物理学中的能量均分定理，则系统的平衡温度可表达为2212iiBfipTkNm?=?(6.1)根据此方法，在晶体中粒子的配置与排布、晶形变化等构造相变可以以用计算机模拟进行研究。根据物理化学的基本气体的基本方程能够知道，体系的NEV系综下的计算体系主要变化的量有两个压力和温度。下文就从恒温和恒压两个条件下对NEV系综的基本方法做简单的推导。1.恒温方法NseS和HooverWG引入与热源相关的参数?，尝
17、试构造恒温状态，亦即具有确定N，V，T值的系统，对此可设想为与大热源接触而到达平衡的系统。由于热源很大，交换能量不会改变热源的温度。在两者建立平衡后，系统将与热源具有一样的温度，系统与热源合起来构成一个复合系统，如图6-3所示。这样，系统中微观粒子的动力学方程：iiidqpdtm=(6.2)iiidppdtdq?=-?(6.3)23222iBexiipdNkTdtmQ?=-?(6.4)式中为系统的势函数，Q为有效质量，?表示热力学摩尔系数。式(6.2)至式(6.4)同往常的分子动力学方法的区别体如今式(6.3)中，即增加了与热源的互相作用相关的并与力的量纲一样的一项(ip?)。与热源相关的变化
18、参数?的运动方程表明，当系统的总能量大于32BNkT时，?是增加的，进而显示出使粒子速度减小那样的作用；反之则显示出使粒子速度增大；Q是表示与温度控制有关的常数。图6-3采取扩展系方法恒温分子动力学方法原理示意图该方法的特征是系统处于微观状态的分布服从正则分布，若在推导式(6.2)、式(6.3)、式(6.4)中，引人由公式lndsdt?=定义的变量，则系统的哈密顿量(即能量)能够表示为203()ln22BQHHNkTs?=+?+(6.5)式中第二项和第三项对应于与热源相联络的部分。这讲明，由于系统与热源存在热接触，二者能够交换能量，因而粒子的能量发生变化(H0H)，同时也讲明系统可能的微观状态
19、可具有不同的能量值。2.恒压方法为了进行压力的控制，考虑如图6-4所示的方法，即利用活塞调控系统的体积变化。对于质量一定的理想气体系统，体积变大，其压力就变小，反之亦然。利用这样的原理，能够通过改变体积使压力到达某一个预定值。就实际的模拟而言，可控制体积大小均匀(三维)地变化，将粒子的坐标、动量表示成如下形式1313iiiiqVqpVp-?=?=?(6.6)进而将粒子的坐标、动量与对分子动力学单元(边长1/3LV=)归一化的参变量iq和ip联络起来。应用归一化的参变量并根据正则方程可得到粒子运动方程。图6-4恒压分子动力学方法原理示意图(实际上，活塞为三维运动)粒子系和活塞组成复合系统，其哈密
20、顿量可由下式给出202vexpHHPVW=+(6.7)212330()2qiiipHVVm-?=+?(6.8)式中，exP是外压；vp是体积V对应的共轭动量；W是决定于体积变化速度的因子。若将由上述哈密顿量导出的运动方程用粒子的实际坐标直接写出来，则有3iiiidqpqdVdtmdtV=+?(6.9)()3iiidppdVdtdtVq?=-?(6.10)222()()3iiiiiiexpqmqdVWPVdt?-?=-(6.11)式(6.9)中右边第一项是根据维里(Virial)定理求出的系统的霎时压P。在此处理中，十分重要的是式(6.9)右边第一项与统计力学压力表达式一样。6.3.3NTP系综
21、质点系的基本方程前面已讨论了关于NEV系综的质点系方程。鉴于NTV和NHP系综的基本方程包含于NTP系综的基本方程之中，在此只对NTP系综的基本方程作一些讲明。NTP系综是在宏观上具有确定的粒子数N、温度T和压力P的系综。在此系综中，T和P可作为外部参数给出。那么，如何才能知足N，T，P恒定的条件呢？显然，粒子数N保持一定是最简单的，由于只要使所考察的基本单元内的粒子数恒定就行了。1.压力恒定体系为简单起见，以由圆柱容器内充满气体和调节活塞构成的系统为例，如图6-5所示。设活塞的重量为W，则在平衡状态下，由活塞的重量产生的外部压力(exP)与由气体产生的内压平衡。假如外压exP变大，为了保持平
22、衡，则气体部分的体积将收缩，进而内压变大。图6-5活塞与气体容器构成的示意图因而，把体积作为系统的动力学参数建立其运动方程，就好似维持压力恒定。其方程可写为22()exdVWPVdt?=-?(6.12)式中，V为体积，是由维里(Virial)定理求出的内压，W是对应于V的假想质量。图6-5为活塞与气体容器构成的示意图。将上述想法推广到一般情况下的晶格系统，如图6-6所示，对一般的晶体，其基本单元可由三个平移矢量(,)xyzaaaa=,(,)xyzbbbb=,(,)xyzcccc=决定其形状和大小。将,abc的三个分量写成矩阵h，该矩阵规定了晶体基本单元的形状和大小，这样h也成为原子的实际坐标i
23、r和晶格矢量is之间的变换矩阵，亦即xyzxyzxyzaaahbbbccc?=?(0.1)iirhs= (0.2)ca图6-6基于单元形状示意图在式(0.2)中，体积V作为动力学参数，而在一般晶体系统中，h应该作为动力学参数。若给出此动力学系统的哈密顿量，则有1121()1(,)()222tTNiiNexiiGtrHrrrPtrGmW-=+LL(6.13)式中，i是与格矢is相应的原子共扼动量，TGhh=为数值行列式，中是多粒子体系的势函数，是h中与体积V相应的共扼动量，是基本单元的体积，且由()()habc=?=?给出，W为相应的质量，exP为外部压强，代表各向异性应力张量。2.温度恒定体系
24、接下来讨论知足温度恒定的基本方法。恒温体系比定压体系要抽象，简单地给出视觉化图像是困难的。势能在原子坐标和动量之外引入了附加自由度f，并由下式用f把现实系统与具有能量确定的(微正则系统)假想系统联络起来。dtdtf=(6.14)式中，t为现实系统的时间，t为假想系统的时间。该假想(力学)系统的哈密顿量H能够写成221221()(,)(ln)22NfiNBexiiPpHrrrgkTfMmf=+?L(6.15)式中，iP是在假想系统中原子的动量，fP是f对应的共扼动量，g是由公式1gN=+给出的假想系统的自由度数目，exT二是热浴的温度，Bk为玻尔兹曼常数。显然，式错误!未找到引用源。中右边第一项
25、代表粒子的动能，第二项为粒子系统的势能，第三项是热浴的假想动能，第四项是热浴的势能；由于在假想系统中，力学系综是微正则分布的，所以其配分函数NEWZ变成具有总能量E的状态总数12121212()NNfNEWNNfHEdrdrdrdpdpdpdfdPdrdrdrdpdpdpdfdP-Z=?LLLLLL(6.16)12121212()NNfNEWNNfHEdrdrdrdpdpdpdfdPdrdrdrdpdpdpdfdP-Z=?LLLLLL(6.17)若将此配分函数积分变换到现实系统中，则可得到在恒定温度T下的正则分布(也称为玻尔兹曼分布)的配分函数NTVZ12121212exp()NNBNTVNN
26、HdrdrdrdpdpdpkTdrdrdrdpdpdp-Z=?LLLLLL(6.18)式中，H是由下式给出的现实系统的哈密顿量2121(,)2NiNiipHrrrm=+L(6.19)H是假想系统的守恒量，而H对真实系统来讲，由于存在热交换将不是守恒量。为了判定计算结果能否合理，只要简单地检验H能否保持恒定(守恒)就行了。3.NTP系综到此，已经导出了压力恒定系综和温度恒定系综的哈密顿量，将两个系综的哈密顿量组合起来就能够得到NTP系综的哈密顿量，其模型如图6-7所示。NTP系综的哈密顿量为1122212()(,)221()(ln)22tTNiiNiifexBexGtrHrrrmfWfPPtrG
27、gkTfM-=+?+?+?L(6.20)基本单元图6-7NTP系统示意图若根据此哈密顿量导出正则方程式，并用速度置换动量，则可得到各变数的运动方程2121()NiiijijiijdsfmsmGGsfdt-=-+&(6.21)22()exdhWfhWPhfdt=-&(6.22)2221()()()()NTTiiiBexidffMmhshsWtrhhgkTfMfdt=+-?+&(6.23)式中，变量上方的圆点表示对时间的微分，而,x分别由下列各式给出1ijijijijdrdr=-(6.24)111()()()NNNiiiijijijiijimhshsxrr=?+?&(6.25)1h-=(6.26)
28、4.NPT系综分子动力学求解若在三维周期边界条件下求解上述方程，则可模拟固体的构造相变、溶解现象和晶化经过。各变量的运动方程式(6.21)，式(6.22)和式(6.23)是非线性常微分方程，不能用贝鲁勒(Berrele)法进行求解。因而，对这些方程通常用牛顿迭代方法求解。而问题的关键在于确定式(6.22)和式(6.23)中的假想质量。所考察基本单元内的粒子体系存在着由粒子作用势以及粒子数与粒子质量共同决定的固有压力和温度等的起伏变化。若想再现这些波动周期那样给出M和W，则由体积变化或温度变化引起的粒子系统对平衡态的偏离将迅速回复而变成稳定状态。因而，若设此周期分别为Pt和Tt，则M和W能够表示
29、为26()2TBtMNkT=(6.27)23()()2PtLWk=(6.28)式中，T为温度，Bk为玻尔兹曼常数，N为原子数，L为基本单元的线尺度，k为压缩率，Pt和Tt分别根据对NEV系综的计算结果而定。【练习与考虑】6-2.推演NEV系综的牛顿方程，哈密顿方程，讲明各个物理量的意义，在那些实际的计算体系中得到应用。6-3.根据NPT系综的推导方法，对NPH系综的哈密顿方程做简单的推导。6-4.具体的分析不同的系综的使用范畴，给出各个系综的使用条件，列表进行简单的比拟。6.4原子势函数和分子力场构造分子动力学中存在着两个基本的概念，一是在前面提到的系综的问题，它所确定的是模拟系统的大小，规格
30、和相关的模拟条件等信息；二是本节需要讲述的势函数。简单的讲，势函数就是用来描绘原子和原子之间互相作用的。随着势函数的不断发展，现有的势函数主要分成了两个大的部分。即是经典的势函数和电子理论的势函数。如图6-8所示，经典分子动力学的基本的势函数的一个简单分类。图6-8经典的分子动力学互相作用势分类作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、
31、Morse势等二体势模型，对于金属计算，主要采用Morse势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因而在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是相对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大的困难，因而在一般应用中，通过第一性原理计算结果拟合势函数的Lennard-Jones，Morse等势模型的应用仍然非常广泛。6.4.1经典理论的原子势函数在分子动力学的框架之下，早期的势函数是主要以描绘原子和原子之间的互相作用而产生和发展起来的。早在1903年的时候，MieG就研究了两个粒子之间的互相作用势，他给出的势函数是由两项构成的，一项代表原子之间的互相吸引，另一项则代表着原子之间的相