《第6章-分子动力学方法汇总.doc》由会员分享，可在线阅读，更多相关《第6章-分子动力学方法汇总.doc（28页珍藏版）》请在三一办公上搜索。
1、第6章 分子动力学方法经典分子动力学方法无疑是材料，尤其是大分子体系和大体系模拟有效的方法之一。分子动力学可以用于NPT，NVE，NVT等不同系综的计算，是一种基于牛顿力学确定论的热力学计算方法。与蒙特卡罗法相比在宏观性质计算上具有更高的准确度和有效性，可以广泛应用于物理，化学，生物，材料，医学等各个领域。本章在介绍分子动力学的基本概念的基础上，简单介绍了分子动力学的基本思想，势函数分类和基本方程。然后介绍了分子动力学的常用系综和典型的NPT，NVE，NVT系综基本方程。结合材料建模中的基本简化方法和技巧，阐述了边界条件和时间积分的数值处理技巧。最后，利用统计力学的基本概念给出分子动力学的计算
2、信息的解析方式。并且结合Materials Explore软件计算分析了CNT的几何结构稳定性。6.1 引言分子动力学方法(Molecular Dynamics, MD)方法是一种按该体系内部的内禀动力学规律来计算并确定位形的变化的确定性模拟方法。首先需要在给定的外界条件下建立一组粒子的运动方程，然后通过直接对系统中的一个个粒子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计力学方法得到多体系统的静态和动态特性，从而获得系统的宏观性质。可以看出，分子动力学方法中不存在任何随机因素，这个也是分子动力学方法和后文要提到的蒙特卡洛方法的区别之一。在分子动力学方
3、法的处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学定律(或者是拉格朗日方程)。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日函数来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现玻尔兹曼的统计力学途径。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是分子动力学方法的计算机程序相对蒙特卡罗较复杂，其计算成本较高。分子动力学方法发展历史改革经历了近60年，分子动力学方法是20世纪50年代后期由Alder B J 和Wainwright T E创造发展的。Alder和Wainwright在19
4、57年利用分子动力学模拟，发现了“刚性球组成的集合系统会发生由其液相到结晶相的相转变”，后来人们称这种相变为Alder相变。其结果表明，不具有引力的粒子系统也具有凝聚态。到20世纪70年代，产生了刚性体系的动力学方法，成功地被应用于水和氮等分子性溶液体系的处理；1972年，Less A W和 Edwards S F 等发展了该方法并扩展到了存在速度梯度，即处于非平衡状态的系统。之后，此方法被 Ginan M J等推广到了具有温度梯度的非平衡系统，从而构造并形成了所谓的非平衡分子动力学方法体系。进人20世纪80年代之后，出现了在分子内部对一部分自由度施加约束条件的新的分子动力学方法，从而使 分子
5、动力学方法可适用于类似蛋白质等生物大分子的解析与设计。分子动力学方法真正作为材料科学领域的一个重要研究方法，开始于由Andersen, Parrinello和Rahman等创立恒压分子动力学方法和Nse等完成恒温分子动力学方法的建立及在应用方面的成功。后来，针对势函数模型化比较困难的半导体和金属等，1985 年人们又提出了将基于密度泛函理论的电子论和分子动力学方法有机统一起来的所谓Car-Parrinello 方法，亦即第一性原理分子动力学方法。这样，分子动力学的方法进一步得到发展和完善，它不仅可以处理半导体和金属的间题，同时还可应用于处理有机物和化学反应。关于Car-Parrinello分子
6、动力学方法将在第7章重点学习讨论。本章重点讨论经典分子动力学方法的基本原理和计算方法。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、趋于平衡的时间和步程，获得好的精度。常用的初始条件可以选择为，令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。图6-2 分子动力学计算流程3. 趋于平衡计算在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，增加或者从系统中移出能量，直到持续给出确定的能量值，称此时系统已经达到平衡，达到
10、平衡的时间称为弛豫时间。 分子动力学中，积分步程的大小选择十分重要，这决定了模拟所需要的时间。为了减小误差，步长要小，但过小系统模拟的弛豫时间就长。因此根据经验选择适当的步长。如，对一个具有几百个氩气分子的体系，采用Lennard-Jones势函数，发现取h为0.01量级，可以得到很好的结果。这里选择的h是没有量纲的，实际上这样选择的h对应的时间在10-14s的量级。如果模拟1000步，系统达到平衡，弛豫时间只有10-11s。4. 宏观物理量的计算实际计算宏观的物理量往往是在模拟的最后阶段进行的，它是沿相空间轨迹求平均来计算得到的，根据各态历经假说，时间平均代替系综平均。关于各态历经假说将在第
11、8章的蒙特卡洛方法中深入讨论。6.2.3 经典分子动力学中近似处理实际的计算体系由于其外在条件的复杂性，往往需要对计算的系统进行一定的近似处理，在经典分子动力学计算中，引进以下近似处理：v 经典粒子相互作用，不考虑电子相互作用量子效应。v 力的作用形式，由参数可调的相互作用势函数决定，并经过实验测定来验证。v 模拟体系与实际体系相差较大，一般需要采用周期边界来扩展计算体系。时间平均是在有限时间内完成。【练习与思考】6-1. 查找文献，根据上述的分子动力学的处理流程图，编写实现分子动力学的简单程序，可参考Daan F和Berend S编著的分子模拟一书。6.3 分子动力学的系综系综(ensemb
12、le)是指在一定的宏观条件下(约束条件)，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统集合，或称为统计系综。系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念；系综是统计理论的一种表述方式，系综理论使统计物理成为普遍的微观统计理论；系综并不是实际的物体，构成系综的系统才是实际物体。 约束条件是由一组外加宏观参量来表示。在平衡统计力学范畴下，可以用来处理稳定系综。6.3.1 常用系综分类 1. 正则系综(canonical ensemble)全称应为“宏观正则系综”，简写为NVT，即表示具有确定的粒子数(N)、体积(V)、温度(T)。正则系综是分子动力学方法模拟处理的
13、典型代表。假定N个粒子处在体积为V的盒子内，将其埋入温度恒为T的热浴中。此时，总能量(E)和系统压强(P)可能在某一平均值附近起伏变化。平衡体系代表封闭系统，与大热源热接触平衡的恒温系统。正则系综的特征函数是亥姆霍兹自由能。 2. 微正则系综(micro-canonical ensemble)简写为NVE，即表示具有确定的粒子数(N)、体积(V)、总能量(E)。微正则系综广泛被应用在分子动力学模拟中。假定N个粒子处在体积为V的盒子内，并固定总能量(E)。此时，系综的温度(T)和系统压强(P)可能在某一平均值附近起伏变化。平衡体系为孤立系统，与外界即无能量交换，也无粒子交换。微正则系综的特征函数
14、是熵。3. 等温等压(constant-pressure, constant-temperature)简写为NPT，即表示具有确定的粒子数(N)、压强(P)、温度(T)。其总能量(E)和系统体积(V)可能存在起伏。体系是可移动系统壁情况下的恒温热浴。特征函数是吉布斯自由能。 4. 等压等焓(constant-pressure, constant-enthalpy)简写为NPH，即表示具有确定的粒子数(N)、压强(P)、焓(H)。由于，故在该系综下进行模拟时要保持压力与焓值为固定，其调节技术的实现也有一定的难度，这种系综在实际的分子动力学模拟中已经很少遇到了。 5. 巨正则系综(grand ca
15、nonical ensemble)简写为VT，即表示具有确定的粒体积(V)、温度(T)和化学势()。巨正则系综通常是蒙特卡罗模拟的对象和手段。此时、系统能量(E)、压强(P)和粒子数(N)会在某一平均值附近有一个起伏。体系是一个开放系统，与大热源大粒子源热接触平衡而具有恒定的温度(T)。特征函数是马休(Massieu)函数J(，V，T)。对于不同的系综，由于其选择的运动方程的求解方法亦有不同。针对不同的计算体系具体的细节比较复杂，下文将以NEV和NTP这两类简单的系综运用两类不同的条件方法下求解的不同方程做简要的介绍。6.3.2 NEV系综基本方程NEV系综是一类较为常见的系综。对于该系综，经
16、典分子动力学方法一般采用差分近似法求解牛顿运动方程，并追踪系统的时间变化。考虑由N个粒子构成的系统，假设第i个粒子在时刻位置，速度为受到的作用力为，为求出该粒子(i)在时刻的位置，经常使用的方法是Verlet算法(对于该算法的具体细节在下文将有详细的推导)。根据系统的动能及统计物理学中的能量均分定理，则系统的平衡温度可表达为(6.1)根据此方法，在晶体中粒子的配置与排布、晶形变化等结构相变也可以用计算机模拟进行研究。根据物理化学的基本气体的基本方程可以知道，体系的NEV系综下的计算体系主要变化的量有两个压力和温度。下文就从恒温和恒压两个条件下对NEV系综的基本方法做简单的推导。1. 恒温方法N
17、se S和Hoover W G引入与热源相关的参数，尝试构造恒温状态，亦即具有确定N，V，T值的系统，对此可设想为与大热源接触而达到平衡的系统。由于热源很大，交换能量不会改变热源的温度。在两者建立平衡后，系统将与热源具有相同的温度，系统与热源合起来构成一个复合系统，如图6-3所示。这样，系统中微观粒子的动力学方程：(6.2)(6.3)(6.4)式中为系统的势函数，为有效质量，表示热力学摩尔系数。式(6.2)至式(6.4)同往常的分子动力学方法的区别体现在式(6.3)中，即增加了与热源的相互作用相关的并与力的量纲相同的一项()。与热源相关的变化参数的运动方程表明，当系统的总能量大于时，是增加的，
18、从而显示出使粒子速度减小那样的作用；反之则显示出使粒子速度增大；是表示与温度控制有关的常数。图6-3 采取扩展系方法恒温分子动力学方法原理示意图该方法的特征是系统处于微观状态的分布服从正则分布，若在推导式(6.2)、式(6.3)、式(6.4)中，引人由公式定义的变量，则系统的哈密顿量(即能量)可以表示为(6.5)式中第二项和第三项对应于与热源相联系的部分。这说明，由于系统与热源存在热接触，二者可以交换能量，因此粒子的能量发生变化(H0H)，同时也说明系统可能的微观状态可具有不同的能量值。2. 恒压方法为了进行压力的控制，考虑如图6-4所示的方法，即利用活塞调控系统的体积变化。对于质量一定的理想
19、气体系统，体积变大，其压力就变小，反之亦然。利用这样的原理，可以通过改变体积使压力达到某一个预定值。就实际的模拟而言，可控制体积大小均匀(三维)地变化，将粒子的坐标、动量表示成如下形式(6.6)从而将粒子的坐标、动量与对分子动力学单元(边长)归一化的参变量和联系起来。应用归一化的参变量并根据正则方程可得到粒子运动方程。图6-4 恒压分子动力学方法原理示意图(实际上，活塞为三维运动)粒子系和活塞组成复合系统，其哈密顿量可由下式给出(6.7)(6.8)式中，是外压；是体积对应的共轭动量；是决定于体积变化速度的因子。若将由上述哈密顿量导出的运动方程用粒子的实际坐标直接写出来，则有(6.9)(6.10
20、)(6.11)式(6.9)中右边第一项是根据维里(Virial)定理求出的系统的瞬间压P。在此处理中，特别重要的是式(6.9)右边第一项与统计力学压力表达式相同。6.3.3 NTP 系综质点系的基本方程前面已讨论了关于NEV系综的质点系方程。鉴于NTV和NHP系综的基本方程包含于NTP系综的基本方程之中，在此只对NTP系综的基本方程作一些说明。NTP系综是在宏观上具有确定的粒子数N、温度T和压力P的系综。在此系综中，T和P可作为外部参数给出。那么，怎样才能满足N，T，P恒定的条件呢？显然，粒子数N保持一定是最简单的，因为只要使所考察的基本单元内的粒子数恒定就行了。1. 压力恒定体系为简单起见，
21、以由圆柱容器内充满气体和调节活塞构成的系统为例，如图6-5所示。设活塞的重量为W，则在平衡状态下，由活塞的重量产生的外部压力()与由气体产生的内压平衡。如果外压变大，为了保持平衡，则气体部分的体积将收缩，从而内压变大。图6-5 活塞与气体容器构成的示意图因此，把体积作为系统的动力学参数建立其运动方程，就好像维持压力恒定。其方程可写为(6.12)式中，为体积，是由维里(Virial)定理求出的内压，是对应于的假想质量。图6-5为活塞与气体容器构成的示意图。将上述想法推广到一般情况下的晶格系统，如图6-6所示，对一般的晶体，其基本单元可由三个平移矢量 , 决定其形状和大小。将的三个分量写成矩阵，该
22、矩阵规定了晶体基本单元的形状和大小，这样也成为原子的实际坐标和晶格矢量之间的变换矩阵，亦即(0.1)(0.2)图6-6 基于单元形状示意图在式(6.15)中，体积作为动力学参数，而在一般晶体系统中，应该作为动力学参数。若给出此动力学系统的哈密顿量，则有(6.13)式中，是与格矢相应的原子共扼动量，为数值行列式，中是多粒子体系的势函数，是中与体积V相应的共扼动量，是基本单元的体积，且由给出，为相应的质量，为外部压强，代表各向异性应力张量。2. 温度恒定体系接下来讨论满足温度恒定的基本方法。恒温体系比定压体系要抽象，简单地给出视觉化图像是困难的。势能在原子坐标和动量之外引入了附加自由度，并由下式用
23、把现实系统与具有能量确定的(微正则系统)假想系统联系起来。(6.14)式中，为现实系统的时间，为假想系统的时间。该假想(力学)系统的哈密顿量可以写成(6.15)式中，是在假想系统中原子的动量，是对应的共扼动量，是由公式给出的假想系统的自由度数目，二是热浴的温度，为玻尔兹曼常数。显然，式(6.18)中右边第一项代表粒子的动能，第二项为粒子系统的势能，第三项是热浴的假想动能，第四项是热浴的势能；因为在假想系统中，力学系综是微正则分布的，所以其配分函数变成具有总能量 E 的状态总数(6.16)(6.17)若将此配分函数积分变换到现实系统中，则可得到在恒定温度T下的正则分布(也称为玻尔兹曼分布)的配分
24、函数 (6.18)式中，是由下式给出的现实系统的哈密顿量(6.19)是假想系统的守恒量，而对真实系统来讲，由于存在热交换将不是守恒量。为了判断计算结果是否合理，只要简单地检验是否保持恒定(守恒)就行了。3. NTP系综到此，已经导出了压力恒定系综和温度恒定系综的哈密顿量，将两个系综的哈密顿量组合起来就可以得到NTP系综的哈密顿量，其模型如图6-7所示。NTP系综的哈密顿量为 (6.20)图6-7 NTP系统示意图若根据此哈密顿量导出正则方程式，并用速度置换动量，则可得到各变数的运动方程(6.21)(6.22)(6.23)式中，变量上方的圆点表示对时间的微分，而分别由下列各式给出(6.24)(6
25、.25)(6.26)4. NPT系综分子动力学求解若在三维周期边界条件下求解上述方程，则可模拟固体的结构相变、溶解现象和晶化过程。各变量的运动方程式(6.21)，式(6.22)和式(6.23)是非线性常微分方程，不能用贝鲁勒(Berrele)法进行求解。因此，对这些方程通常用牛顿迭代方法求解。而问题的关键在于确定式(6.22) 和式(6.23)中的假想质量。所考察基本单元内的粒子体系存在着由粒子作用势以及粒子数与粒子质量共同决定的固有压力和温度等的起伏变化。若想再现这些波动周期那样给出和，则由体积变化或温度变化引起的粒子系统对平衡态的偏离将迅速回复而变成稳定状态。因此，若设此周期分别为和，则和
26、可以表示为(6.27)(6.28)式中，为温度，为玻尔兹曼常数，为原子数，为基本单元的线尺度，为压缩率，和分别根据对NEV系综的计算结果而定。【练习与思考】6-2. 推演NEV系综的牛顿方程，哈密顿方程，说明各个物理量的意义，在那些实际的计算体系中得到应用。6-3. 根据NPT系综的推导方法，对NPH系综的哈密顿方程做简单的推导。6-4. 详细的分析不同的系综的使用范畴，给出各个系综的使用条件，列表进行简单的比较。6.4 原子势函数和分子力场构造分子动力学中存在着两个基本的概念，一是在前面提到的系综的问题，它所确定的是模拟系统的大小，规格和相关的模拟条件等信息；二是本节需要讲述的势函数。简单的
27、讲，势函数就是用来描述原子和原子之间相互作用的。随着势函数的不断发展，现有的势函数主要分成了两个大的部分。即是经典的势函数和电子理论的势函数。如图6-8所示，经典分子动力学的基本的势函数的一个简单分类。图6-8 经典的分子动力学相互作用势分类作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、Morse势等二体势模型，对于金属计算，主要采用Mor
28、se势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是相对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大的困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的Lennard-Jones，Morse等势模型的应用仍然非常广泛。6.4.1 经典理论的原子势函数在分子动力学的框架之下，早期的势函数是主要以描述原子和原子之间的相互作用而产生和发展起来的。早在1903年的时候，Mie G就研究了两个粒子之间的相互作用势，他给出的势函数是
29、由两项构成的，一项代表原子之间的相互吸引，另一项则代表着原子之间的相互排斥。这就是非常著名的Lennard-Jones势函数的原始思考。对于其他有关势函数的发展在这里就不做过多的赘述了。从图6-9可以清楚的知道原子之间的相互作用可以分为四类：对势、对泛函势、组合势、组合泛函势。下面就对相关的势函数的基本的发展做一个简单的总结。表6-1 原子中相互作用势函数势函数作用因素典型代表对势仅仅决定于原子的坐标L-J势、BMH势、GK势对泛函势多体相互作用EAM势有效介质理论组合势共价键的相互作用SW势组合泛函势Abell-Tersoff势CheliKowsky-Phillips势6.4.2 分子力场分
30、子力场根据量子力学的波恩奥本海默近似，一个分子的能量可以近似看作构成分子的各个原子的空间坐标的函数，简单地讲就是分子的能量随分子构型的变化而变化，而描述这种分子能量和分子结构之间关系的就是分子力场函数。分子力场函数为来自实验结果的经验公式，可以讲对分子能量的模拟比较粗糙，但是相比于精确的量子力学从头计算方法，分子力场方法的计算量要小数十倍，而且在适当的范围内，分子力场方法的计算精度与量子化学计算相差无几，因此对大分子复杂体系而言，分子力场方法是一套行之有效的方法。以分子力场为基础的分子力学计算方法在分子动力学、蒙特卡罗方法、分子对接等分子模拟方法中有着广泛的应用。 1. 力场主要讨论的对象及其
31、构成：以下参数构成一套力场函数体系需要有一套联系分子能量和构型的函数，还需要给出各种不同原子在不同成键状况下的物理参数，比如正常的键长、键角、二面角等，这些力场参数多来自实验或者量子化学计算。 v 键伸缩能：构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化；v 键角弯曲能：键角变化引起的分子能量变化；v 二面角扭曲能：单键旋转引起分子骨架扭曲所产生的能量变化； v 非键相互作用：包括范德华力、静电相互作用等与能量有关的非键相互作用； v 交叉能量项：上述作用之间耦合引起的能量变化。2. 常用力场函数和分类不同的分子力场会选取不同的函数形式来描述上述能量与体系构型之间的关系。到目前，不
32、同的研究小组设计了很多适用于不同体系的力场函数，根据他们选择的函数和力场参数，可以分为以下几类：力场CHARMM力场CVFF力场MMX力场AMBER力场COMPASS力场MMF94力场CFF力场UFF力场Dreiding力场ESFF力场第一代力场第二代力场第三代力场图6-9 分子力场的发展和分类v 传统力场(第一代力场)a) AMBER力场：由Kollman课题组开发的力场，是目前使用比较广泛的一种力场，适合处理生物大分子；b) CHARMM力场：由Karplus课题组开发，对小分子体系到溶剂化的大分子体系都有很好的拟合；c) CVFF力场：CVFF力场是一个可以用于无机体系计算的力场；d)
33、MMX力场：MMX力场包括MM2和MM3，是目前应用最为广泛的一种力场，主要针对有机小分子。v 第二代力场 第二代的势能函数形式比传统力场要更加复杂，涉及的力场参数更多，计算量也更大，当然也相应地更加准确；a) CFF力场CFF力场是一个力场家族，包括了CFF91、PCFF、CFF95等很多力场，可以进行从有机小分子、生物大分子到分子筛等诸多体系的计算；b) COMPASS力场由MSI公司开发的力场，擅长进行高分子体系的计算；c) MMF94力场Hagler开发的力场，是目前最准确的力场之一。v 通用力场(第一代力场)通用力场也叫基于规则的力场，它所应用的力场参数是基于原子性质计算所得，用户可
34、以通过自主设定一系列分子作为训练集来生成合用的力场参数；a) ESFF力场MSI公司开发的力场，可以进行有机、无机分子的计算；b) UFF力场可以计算周期表上所有元素的参数；c) Dreiding力场适用于有机小分子、大分子、主族元素的计算。3. 分子力场的计算对于分子之间的相互作用通常需要将计算体系的所有粒子纳入考虑范畴之内，那么分子之间的相互作用可以作为长程作用来处理。在长程相互作用中又分为周期性长程相互作用力和非周期性相互作用力。周期性的长程作用力通常可以利用无网格算法，包含Ewald求和算法等；而对于非周期性的相互作用通常使用快速多极子方法计算。下面针对这两个最为常见的算法做简单的叙述
35、。v Ewald求和算法Ewald求和算法是Ewald于1921年提出，最初是以计算离子晶体的能量。此方法选定一个计算系统的盒子，盒子中的粒子除了和盒子内部的粒子相互作用之外，还会与其镜像系统中的粒子发生作用，并且镜像系统和盒子内部的系统完全相同。图显示的是计算的相互作用的系统，中央的黑色的区域表示的是计算的区域，周围的镜像随着系统的距离逐渐变淡，那么这样的一个极限作用区域就像一个球，这个球称为Ewald球。假设作用的系统的边长是，含有个相互作用的粒子，其镜像的作用位置可以表示为()；其中都是正整数。v 快速多极子法(The Fast Multipole Method)在一个开放的系统中，不是
36、对所谓的晶格进行求和，而是对要求计算的整个体系的粒子对进行求和。那么这样的一个系统不同位置粒子的静电能可以简单的表述为(6.29)因此这个体系的静电作用和前文所讲述的Lennard-Jones作用项具有相同的形式。快速多极子法将需要计算的体系分为多个体积相同的格胞，由格胞内部的所有的原子计算每一个格胞的多极(电荷，偶极，四极)。而格胞和格胞之间的相互作用是利用中央多极展开的方法进行计算。因此中央多极展开的方法的条件是格胞之间的距离必须是大于一定距离的，因为大于一个格胞以上的以多极展开的方式计算相互作用的力，而格胞内部的原子对是以成对原子的作用力来计算的。【练习与思考】6-5. 查找文献，推导L
37、ennard-Jones势函数的公式。6-6. 推导EAM势函数的基本公式，并且实现简单的算法，编写基本设计程序。6-7. 以对势中的Lennard-Jones势函数为例，写出对势中势函数的基本算法，并且进行基本的程序设计。6.5 数值方法与编程技巧6.5.1 边界条件随着材料的低维化和纳米化，材料的表面和边界对于其性质的影响愈来愈大。在处理物质表面、界面，以及块体物质时，对其周期边界条件问题的考虑就变得非常必要了。对于周期性边界条件可以分为一维、二维及三维的情况。介绍薄膜(可使用二维周期边界条件)和扩大到半无限的表面情况，然后讨论液体与晶体、晶体与晶体之间形成的界面等，最后学习块状物质的三维
38、周期边界条件。1. 薄膜的边界条件(二维周期边界条件)对于薄膜情况，如图6-10所示赋予其周期性边界条件。可以认为薄膜在x-y平面内无限扩展，存在周期边界条件，而在 z方向受到限制(不赋予周期边界条件)。事实上，二维周期性边界条件也适用于在z方向不具有真正意义上的自由度的二维物质。但就物质科学的研究范畴而言，多数情况下是指在某种衬底上生成的原子或分子的薄膜。此处讨论的边界条件，除针对薄膜以外，也可近似地用于下面将要讲到的在 z 的负方向半无限地扩展的表面问题，例如表面重构等。但是，当在z方向存在异质原子层时，原来块体物质在异质原子层附近将有近于表面的性质，这种由块体向表面的过渡特征取决于原子层
39、的厚度，一般要求原子层厚度在力的中断距离以上(作用力程以上)。在处理实际问题时，只要尝试进行几次改变原子层厚度的模拟和计算，并检验计算结果是否变化就行了。图6-10 薄膜周期性边界条件示意图2. 表面边界条件(半无限的情况)对于表面存在表面重构的情况，是与薄膜的情况不同的，其在 z 轴方向扩展至半无限。处理此类系统之边界条件时，只对x-y平面和 z 的负方向赋予相应的周期性边界条件，其示意图如图6-11所示。但是，在将基本单元原样地复制中，必须考虑到这样的一种情况，那就是在距薄膜结构的基本单元最下层的相邻底层被配置了基本单元表面层的复制品，从而导致了在本来应具有块体物质性质的基本单元的最下层的
40、结构中出现了表面层的强烈影响。所以将基本单元对z轴反转(即镜像映射)进行复制后使用。然而，即使在这种情况下，考虑到基本单元在z方向的厚度对计算结果有影响，所以应当进行与前面所述同样的检验。图6-11 采用映射方法半无限表面的周期性边界条件图6-12 固定层法半无限表面的周期性边界条件在利用分子动力学方法对分子束外延的晶体生长模拟时，要充分认识到分子束外延 (MBE)是一个原子数不断增加的非平衡系。所谓MBE技术就由束源产生的原子蒸发到衬底上，从而形成晶体。因为由束源出来的原子具有动量，所以碰撞衬底，衬底的重心就变得下移，为了克服这种现象的影响，对分子束外延即在z轴方向不赋予周期性边界条件，而采
41、用通过调整配置使原子位置坐标不变的所谓固定原子层方法，然后使用边界条件，如图6-12所示。固定(原子)层还可理解为具有无限大质量的原子层。来自固定层中原子作用力的贡献可作为外力处理。为了基本单元最下层具有块体物质的性质，固定层中的原子层数目必须足够多。采用镜像映射的方法可自然地处理体系伴随温度变化而产生的膨胀和收缩。但就使用固定层的方法而言，基本单元发生变化的同时，必须将固定层的晶格常数作为外部参量变化来考虑。【练习与思考】6-8. 通过对比映射方法和固定层法的不同，理解在分子动力学中的二维边界条件的设置方法，思考在线形的一维约束条件的边界的处理方法。6-9. 根据界面的边界条件的设置方法，分
42、析讨论在分子动力学中如何构造线缺陷模型，做要的图示和语言表述，说明界面的条件的选取方法。6-10. 根据上述的边界的处理方法，试思考块体材料的边界条件的实现方法，并且给出简单的实例。6.5.2 时间积分处理计算机模拟方法的基点是利用现代计算机高速和精确的优点，对几百个以至上千分子的运动方程进行数值积分。有许多不同的积分方法，它们的效率和方便程度各异。基本问题就是用有限差分法来对二阶常微分方程进行积分。而对于分子动力学而言其基本的工作也就是对运动方程进行时间上的积分。这样需要将运动方程离散化为有限的差分方程，即是将时间离散化为有限的格点。主要的方法是利用泰勒展开采取不同的截断半径。而相邻的格点之
43、间的距离就是所谓的时间步长，这样在知道了原始的位置坐标之后，在知道某些物理量对于时间t的微分量之后就可以很容易确定下一个时间体系的坐标。这里很容易出现所谓的两类误差，即是截断误差和舍入误差。为此我们需要在实际的模拟中，应该根据不同的问题来确定所需要的精度，为此发展了几种时间积分的方法。1. Verlet方法Verlet法经常用于NEV系综的分子动力学微分方程的数值积分求解。若将在和的位置坐标分别用时刻t的位置坐标进行泰勒展开，则有(6.30)(6.31)为简单起见，以后均省略有关原子标识符号I，V(t)为速度，为了消去速度项，将式(6.30)和式(6.31)相加，并将r()移至右边，同时利用，
44、则有(6.32)上式(6.32)是差分公式，其中()的奇次项全被抵消，同时忽略及其他高次项。差分精度(误差)为0()。若由式(6-36)减去式(6-37)则可得到速度表达式(6.33)因为是消去了关于的偶数方项，所以速度的误差为的量级。上述式(6.32)右边第一、二项均是的量级项，而第三项是量级项。这样若将大的绝对值项和小的绝对项混合，则随着重复计算次数的增多，其误差发生积累。为防止误差积累，Hockney提出了所谓蛙跃(Leap-frog)法。2. 蛙跃法首先，利用速度将二阶常微分方程转换为二个一阶常微分方程。亦即(6.34)速度的微分由时刻，和的差分给出(6.35)根据此式，在时刻的速度给
45、出如下(6.36)另一方面，利用时刻和时刻的差分可给出位置矢量的微分(6.37)根据上式可给出时刻的位置矢量表达式(6.38)上式(6.38)所表示的意义是速度和位置矢量的变化依赖步长。利用时刻的位置坐标可计算出力，利用此力的结果求出时刻()的速度，进而，利用此速度求出在()时刻的位置坐标，这样就完成了第一步的计算。正如以上所述，如果使用Verlet法或蛙跃法，则可求解关于NEV系综的线性常微分方程。但是，在此之外的NTV，NHP，NTP等系综给出的是非线性常微分方程。所以用Verlet法或蛙跃法不能进行求解。关于这些非线性常微分方程的求解要采用下面介绍的所谓Gear方法。Hockey提出的L
46、eap-frog算法是Verlet算法的变化，这种方法涉及半时间间隔的速度。这种算法与Verlet算法相比有两个优点：(1)包括显速度项；(2)收敛速度快，计算量小。这种算法明显的缺陷是位置和速度不同步。3. Gear方法因为原子轨迹是时间的连续函数，所以可将()时刻的位置、速度、加速度等对时刻泰勒展开(6.39)若使用给出的新的位置矢量，则可计算()时刻的力F()和加速度()。根据加速度可评价预测的加速度的误差，亦即(6.40)若将此差值加于预测项上，则可得到修正项(6.41)(6.42)4. Velocity-Verlet算法这种算法可以同时给出位置、速度和加速度，并且不牺牲精度，给出了显
47、速度项，计算量适中，目前应用比较广泛。(6.43)5. Gear的预测-校正算法这种算法分三步来完成。首先，根据泰勒展开，预测新的位置、速度和加速度。然后，根据新的计算的力计算加速度。最后这个加速度再由与泰勒级数展开式中的加速度进行比较，两者之差在校正步里用来校正位置和速度项。这种方法的缺点就是占有计算机的内存大。【练习与思考】6-11. 查找文献，总结分子动力学中现有的时间算法，列表说明各种不同的时间算法优势。6-12. 推导出Verlet算法的基本方程，并编程实现。6-13. Hockey提出的Leap-frog算法是Verlet算法的变化，这种方法涉及半时间间隔的速度，查找原始文献，推导出Leap-frog算法的基本方程，结合练习6-12中的问题，修改程序实现Leap-frog算法。6-14. 根据文中的叙述，详细的推导出Gear算法的基本方程，并编程实现。6-15. 查阅文献，分析讨论时间可逆算法刘维公式在分子动力学中应用。6.6 计算结果的解析方法总体来说分子动力学(MD)目前能够讨论的问题已经可以涉及到比较广阔的范畴，从基本的模拟体系的位形，到现在的各种材料的光谱学
