第6章分子动力学方法.pdf

《第6章分子动力学方法.pdf》由会员分享，可在线阅读，更多相关《第6章分子动力学方法.pdf（28页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、第 6 章分子动力学方法- 127 -第第6 6章章分子动力学方法分子动力学方法经典分子动力学方法无疑是材料， 尤其是大分子体系和大体系模拟有效的方法之一。 分子动力学可以用于 NPT ，NVE ，NVT等不同系综的计算，是一种基于牛顿力学确定论的热力学计算方法。 与蒙特卡罗法相比在宏观性质计算上具有更高的准确度和有效性， 可以广泛应用于物理，化学，生物，材料，医学等各个领域。本章在介绍分子动力学的基本概念的基础上，简单介绍了分子动力学的基本思想， 势函数分类和基本方程。 然后介绍了分子动力学的常用系综和典型的 NPT ，NVE ，NVT系综基本方程。结合材料建模中的基本简化方法和技巧，阐述了

2、边界条件和时间积分的数值处理技巧。 最后，利用统计力学的基本概念给出分子动力学的计算信息的解析方式。并且结合Materials Explore软件计算分析了 CNT 的几何结构稳定性。6.16.1 引言引言分子动力学方法(MolecularDynamics, MD) 方法是一种按该体系内部的内禀动力学规律来计算并确定位形的变化的确定性模拟方法。 首先需要在给定的外界条件下建立一组粒子的运动方程， 然后通过直接对系统中的一个个粒子运动方程进行数值求解， 得到每个时刻各个分子的坐标与动量， 即在相空间的运动轨迹， 再利用统计力学方法得到多体系统的静态和动态特性，从而获得系统的宏观性质。 可以看出，

3、分子动力学方法中不存在任何随机因素，这个也是分子动力学方法和后文要提到的蒙特卡洛方法的区别之一。 在分子动力学方法的处理过程中， 方程组的建立是通过对物理体系的微观数学描述给出的。 在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学定律(或者是拉格朗日方程)。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日函数来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现玻尔兹曼的统计力学途径。 这种方法可以处理与时间有关的过程， 因而可以处理非平衡态问题。 但是分子动力学方法的计算机程序相对蒙特卡罗较复杂， 其计算成本较高。分子动力学方法发展历史改革经历了近60 年，分子动力学方

4、法是 20 世纪 50 年代后期由 Alder B J和 Wainwright T E 创造发展的。 Alder和 Wainwright在 1957 年利用分子动力学模拟，发现了“ 刚性球组成的集合系统会发生由其液相到结晶相的相转变” ，后来人们称这种相变为 Alder相变。 其结果表明， 不具有引力的粒子系统也具有凝聚态。 到 20 世纪 70 年代，产生了刚性体系的动力学方法，成功地被应用于水和氮等分子性溶液体系的处理； 1972 年，Less A W 和 Edwards S F等发展了该方法并扩展到了存在速度梯度，即处于非平衡状态的系统。之后，此方法被 Ginan M J等推广到了具有温

5、度梯度的非平衡系统， 从而构造并形成了所谓的非平衡分子动力学方法体系。 进人 20 世纪 80 年代之后， 出现了在分子内部对一部分自由度施加约束条件的新的分子动力学方法，从而使 分子动力学方法可适用于类似蛋白质等生物大分子的解析与设计。分子动力学方法真正作为材料科学领域的一个重要研究方法，开始于由 Andersen,Parrinello和 Rahman 等创立恒压分子动力学方法和N se等完成恒Computational Materials Science：From Basic Principle to Practical Design Methodology- 128 -第 6 章分子动力

6、学方法温分子动力学方法的建立及在应用方面的成功。 后来， 针对势函数模型化比较困难的半导体和金属等，1985 年人们又提出了将基于密度泛函理论的电子论和分子动力学方法有机统一起来的所谓 Car-Parrinello方法，亦即第一性原理分子动力学方法。这样，分子动力学的方法进一步得到发展和完善， 它不仅可以处理半导体和金属的间题， 同时还可应用于处理有机物和化学反应。关于 Car-Parrinello分子动力学方法将在第 7 章重点学习讨论。本章重点讨论经典分子动力学方法的基本原理和计算方法。6.26.2 分子动力学的计算框架分子动力学的计算框架6.2.16.2.1 基本思想基本思想分子动力学方

7、法是沿用牛顿运动方程或拉格朗日方程、 哈密顿方程， 通过考察粒子的运动来研究多粒子系统的物理性质。在处理孤立粒子或原子团簇(不考虑与其他粒子的相互作用)时，单纯地对牛顿运动方程进行积分求解就行了；而在处理凝聚系统时，可以认为将所考虑的 N 个粒子放入具有一定体积的容器中，在这样组成的封闭系中，建立直角坐标(x,y,z) ，每个粒子的位置就由三个坐标分量决定，通过求解3N 个联立方程组就可以得到保守系的总能量。分子动力学的基本思想在于通过设定原子之间的相互作用(势函数)和相关的系综(亦即作用对象和条件)，确定其基本的模拟范畴。在给定的牛顿方程、拉格朗日方程或者是哈密顿方程进行时间的迭代， 在达到

8、指定的力学收敛条件之后得到一个最终的位置坐标， 然后通过相关的位形信息和运动的速度和加速度等信息通过热力统计学的方法， 给出模拟体系的统计信息，主要包括复杂的光学，电学和一些其他的物理信息。其基本思想，总结如图6-1所示。势函数，相互作用系综，外界条件运动方程（牛顿方程，哈密顿方程，拉格朗日方程）原子的坐标位置三维结构动力学性质原子的坐标和速度原子运动光学性质热力学性质初始化条件动力学方程一次信息图 6-1分子动力学的一般思想二次信息6.2.26.2.2 计算流程计算流程总体来说，分子动力学的基本计算有如下四个步骤：1.1. 模型的设定，也就是势函数的选取模型的设定，也就是势函数的选取势函数的

9、研究和物理系统中对物质的描述研究息息相关。 分子动力学的模拟最早使用是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是 Lennard-Jones 势函数，还有嵌入原子 EAM势函数，不同的物质状态描述用不同的势函数。模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。关于势函数选取将在6.4节进一步讨论。Computational Materials Science：From Basic Principle to Practical Design Methodology第 6 章分子动力学方法- 129 -2.2. 给定初始条件给定初始条件运动方程的求解需要知道粒子的初始

10、位置和速度， 不同的算法要求不同的初始条件。 如Verlet 算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。一般意义上讲系统的初始条件不可能知道， 实际上也不需要精确选择代求系统的初始条件，因为当模拟时间足够长时，系统就会忘掉初始条件(对于无记忆的体系而言)。当然，合理的初始条件可以加快系统趋于平衡的时间和步程，获得好的精度。常用的初始条件可以选择为， 令初始位置在差分划分网格的格子上， 初始速度则从玻尔兹曼分布随机抽样得到； 令初始位置随机的偏离差分划分网格的格子上， 初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔

11、兹曼分布随机抽样得到。设定坐标和相互作用t+t计算作用在原子上的力Vi(tt/2)Vi(tt/2)t Fi(t)/miRi(tt)Ri(t)t Vi(tt/2)计算物理量并对其结果进行统计NttmaxY图6-23.3. 趋于平衡计算趋于平衡计算结束分子动力学计算流程在边界条件和初始条件给定后就可以解运动方程， 进行分子动力学模拟。 但这样计算出的系统是不会具有所要求的系统能量， 并且这个状态本身也还不是一个平衡态。 为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，增加或者从系统中移出能量， 直到持续给出确定的能量值，称此时系统已经达到平衡，达到平衡的时间称为弛豫时间。分子动力学中，积分

12、步程的大小选择十分重要， 这决定了模拟所需要的时间。 为了减小误差，步长要小，但过小系统模拟的弛豫时间就长。因此根据经验选择适当的步长。如，对一个具有几百个氩气分子的体系，采用Lennard-Jones 势函数，发现取h 为 0.01量级，可以得到很好的结果。这里选择的h 是没有量纲的，实际上这样选择的h 对应的时间在 10-14s的量级。如果模拟 1000 步，系统达到平衡，弛豫时间只有10-11s。4.4. 宏观物理量的计算宏观物理量的计算Computational Materials Science：From Basic Principle to Practical Design Met

13、hodology- 130 -第 6 章分子动力学方法实际计算宏观的物理量往往是在模拟的最后阶段进行的， 它是沿相空间轨迹求平均来计算得到的，根据各态历经假说，时间平均代替系综平均。 关于各态历经假说将在第8 章的蒙特卡洛方法中深入讨论。6.2.36.2.3 经典分子动力学中近似处理经典分子动力学中近似处理实际的计算体系由于其外在条件的复杂性，往往需要对计算的系统进行一定的近似处理，在经典分子动力学计算中，引进以下近似处理：经典粒子相互作用，不考虑电子相互作用量子效应。力的作用形式，由参数可调的相互作用势函数决定，并经过实验测定来验证。模拟体系与实际体系相差较大，一般需要采用周期边界来扩展计算

14、体系。时间平均是在有限时间内完成。【练习与思考】【练习与思考】6-1. 查找文献，根据上述的分子动力学的处理流程图，编写实现分子动力学的简单程序，可参考 Daan F和 Berend S编著的分子模拟一书。6.36.3 分子动力学的系综分子动力学的系综系综(ensemble) 是指在一定的宏观条件下(约束条件)，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统集合，或称为统计系综。系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念； 系综是统计理论的一种表述方式， 系综理论使统计物理成为普遍的微观统计理论；系综并不是实际的物体，构成系综的系统才是实际物体。约束条件是由一

15、组外加宏观参量来表示。在平衡统计力学范畴下，可以用来处理稳定系综。6.3.16.3.1 常用系综分类常用系综分类1.1. 正则系综正则系综(canonical ensemble)(canonical ensemble)全称应为“ 宏观正则系综” ，简写为NVT ，即表示具有确定的粒子数(N)、体积(V)、温度(T) 。正则系综是分子动力学方法模拟处理的典型代表。假定N 个粒子处在体积为 V 的盒子内，将其埋入温度恒为T 的热浴中。此时，总能量(E) 和系统压强(P) 可能在某一平均值附近起伏变化。平衡体系代表封闭系统， 与大热源热接触平衡的恒温系统。 正则系综的特征函数是亥姆霍兹自由能F (N

16、 , V ,T )。2.2. 微正则系综微正则系综(micro-canonical ensemble)(micro-canonical ensemble)简写为 NVE ，即表示具有确定的粒子数(N)、体积(V)、总能量(E) 。微正则系综广泛被应用在分子动力学模拟中。 假定 N 个粒子处在体积为 V 的盒子内， 并固定总能量(E) 。 此时，系综的温度(T) 和系统压强(P) 可能在某一平均值附近起伏变化。平衡体系为孤立系统，与外界即无能量交换，也无粒子交换。微正则系综的特征函数是熵S(N , V ,E)。3.3. 等温等压等温等压(constant-pressure, constant-t

17、emperature)(constant-pressure, constant-temperature)简写为 NPT ，即表示具有确定的粒子数(N)、压强(P) 、温度(T) 。其总能量(E) 和系统体积(V)可能存在起伏。体系是可移动系统壁情况下的恒温热浴。特征函数是吉布斯自由能G (N ,P,T )。Computational Materials Science：From Basic Principle to Practical Design Methodology第 6 章分子动力学方法- 131 -4.4. 等压等焓等压等焓(constant-pressure, constant-e

18、nthalpy)(constant-pressure, constant-enthalpy)简写为 NPH ，即表示具有确定的粒子数(N)、压强(P) 、焓(H)。由于HEPV，故在该系综下进行模拟时要保持压力与焓值为固定， 其调节技术的实现也有一定的难度， 这种系综在实际的分子动力学模拟中已经很少遇到了。5.5. 巨正则系综巨正则系综(grand canonical ensemble)(grand canonical ensemble)简写为 VT，即表示具有确定的粒体积(V)、温度(T) 和化学势( )。巨正则系综通常是蒙特卡罗模拟的对象和手段。此时、系统能量 (E) 、压强(P) 和粒子

19、数(N)会在某一平均值附近有一个起伏。体系是一个开放系统，与大热源大粒子源热接触平衡而具有恒定的温度(T) 。特征函数是马休(Massieu) 函数 J( ，V ，T)。对于不同的系综， 由于其选择的运动方程的求解方法亦有不同。 针对不同的计算体系具体的细节比较复杂， 下文将以 NEV 和 NTP 这两类简单的系综运用两类不同的条件方法下求解的不同方程做简要的介绍。6.3.26.3.2 NEVNEV系综基本方程系综基本方程NEV系综是一类较为常见的系综。对于该系综，经典分子动力学方法一般采用差分近似法求解牛顿运动方程，并追踪系统的时间变化。考虑由 N 个粒子构成的系统，假设第 i个粒子在时刻x

20、i(t)位置，速度为vi(t)受到的作用力为Fi(t)，为求出该粒子(i) 在时刻tt的位置，经常使用的方法是Verlet 算法(对于该算法的具体细节在下文将有详细的推导)。根据系统的动能及统计物理学中的能量均分定理，则系统的平衡温度可表达为2TkNfB1pi22mii(6.1)根据此方法， 在晶体中粒子的配置与排布、 晶形变化等结构相变也可以用计算机模拟进行研究。根据物理化学的基本气体的基本方程可以知道，体系的NEV 系综下的计算体系主要变化的量有两个压力和温度。下文就从恒温和恒压两个条件下对NEV 系综的基本方法做简单的推导。1.1. 恒温方法恒温方法N se S和 Hoover W G

21、引入与热源相关的参数 ， 尝试构造恒温状态， 亦即具有确定 N ，V ，T 值的系统，对此可设想为与大热源接触而达到平衡的系统。由于热源很大，交换能量不会改变热源的温度。在两者建立平衡后，系统将与热源具有相同的温度， 系统与热源合起来构成一个复合系统，如图6-3所示。这样，系统中微观粒子的动力学方程：dqipidtmidpipidtdqi(6.2)(6.3)(6.4)pi232dN kBTexdt2i2miQComputational Materials Science：From Basic Principle to Practical Design Methodology- 132 -第 6

22、 章分子动力学方法式中为系统的势函数，Q为有效质量， 表示热力学摩尔系数。式(6.2) 至式(6.4) 同往常的分子动力学方法的区别体现在式(6.3) 中， 即增加了与热源的相互作用相关的并与力的量纲相同的一项(pi)。 与热源相关的变化参数 的运动方程表明， 当系统的总能量大于NkBT时，是增加的，从而显示出使粒子速度减小那样的作用；反之则显示出使粒子速度增大；Q是表示与温度控制有关的常数。绝热壁足够大的热源所研究的粒子体系温度Tex32热交换图6-3采取扩展系方法恒温分子动力学方法原理示意图该方法的特征是系统处于微观状态的分布服从正则分布，若在推导式 (6.2) 、式(6.3) 、式(6.

23、4) 中，引人由公式d lns定义的变量，则系统的哈密顿量(即能量)可以表示为dtQ3HH0( )2NkBT lns22(6.5)式中第二项和第三项对应于与热源相联系的部分。 这说明，由于系统与热源存在热接触， 二者可以交换能量，因此粒子的能量发生变化(H0H)，同时也说明系统可能的微观状态可具有不同的能量值。2.2. 恒压方法恒压方法为了进行压力的控制，考虑如图6-4所示的方法，即利用活塞调控系统的体积变化。对于质量一定的理想气体系统，体积变大，其压力就变小，反之亦然。利用这样的原理，可以通过改变体积使压力达到某一个预定值。就实际的模拟而言，可控制体积大小均匀 (三维)地变化，将粒子的坐标、

24、动量表示成如下形式1qiV3qi1pV3pii(6.6)从而将粒子的坐标、动量与对分子动力学单元(边长LV1/3)归一化的参变量qi和pi联系起来。应用归一化的参变量并根据正则方程可得到粒子运动方程。Computational Materials Science：From Basic Principle to Practical Design Methodology第 6 章分子动力学方法- 133 -所研究粒子体系活塞外压内压图6-4恒压分子动力学方法原理示意图(实际上，活塞为三维运动)粒子系和活塞组成复合系统，其哈密顿量可由下式给出pv2HH0PexV2WH0Vi23(6.7)(6.8)p

25、i22mi1q3( V)式中，Pex是外压；pv是体积V对应的共轭动量；W是决定于体积变化速度的因子。若将由上述哈密顿量导出的运动方程用粒子的实际坐标直接写出来，则有dqipidVqidtmidt 3V(6.9)(6.10)dpidVpi()dtdt3Vqipi2qi()2imiqid ViW(2)Pex3Vdt(6.11)式(6.9) 中右边第一项是根据维里(Virial)定理求出的系统的瞬间压P。在此处理中，特别重要的是式(6.9) 右边第一项与统计力学压力表达式相同。6.3.36.3.3 NTPNTP系综质点系的基本方程系综质点系的基本方程前面已讨论了关于 NEV 系综的质点系方程。鉴于

26、NTV 和 NHP 系综的基本方程包含于NTP 系综的基本方程之中，在此只对NTP 系综的基本方程作一些说明。NTP 系综是在宏观上具有确定的粒子数 N 、温度 T 和压力 P 的系综。在此系综中，T 和 P 可作为外部参数给出。那么，怎样才能满足N ，T，P 恒定的条件呢？显然，粒子数N 保持一定是最简单的，因为只要使所考察的基本单元内的粒子数恒定就行了。1.1. 压力恒定体系压力恒定体系为简单起见，以由圆柱容器内充满气体和调节活塞构成的系统为例，如图6-5所示。设活塞的重量为 W ，则在平衡状态下，由活塞的重量产生的外部压力(Pex)与由气体产生的内压平衡。如果外压Pex变大，为了保持平衡

27、，则气体部分的体积将收缩，从而内压变大。Computational Materials Science：From Basic Principle to Practical Design Methodology- 134 -第 6 章分子动力学方法Pex图6-5活塞与气体容器构成的示意图因此，把体积作为系统的动力学参数建立其运动方程， 就好像维持压力恒定。其方程可写为d2VW2(Pex) Vdt(6.12)式中，V为体积，是由维里(Virial)定理求出的内压，W是对应于V的假想质量。图 6-5为活塞与气体容器构成的示意图。 将上述想法推广到一般情况下的晶格系统， 如图 6-6所示，对一般的晶体

28、，其基本单元可由三个平移矢量a(ax,ay,az),b(bx,by,bz),c(cx,cy,cz)决定其形状和大小。将a,b,c的三个分量写成矩阵h，该矩阵规定了晶体基本单元的形状和大小，这样h也成为原子的实际坐标ri和晶格矢量si之间的变换矩阵，亦即axhbxcxaybycyazbzcz错误！未找到引用源。错误！未找到引用源。错误！未找到引用源。错误！未找到引用源。rihsicba图6-6基于单元形状示意图在式(0.2) 中，体积V作为动力学参数，而在一般晶体系统中，h应该作为动力学参数。若给出此动力学系统的哈密顿量，则有Hi 1NitG1i2mi(r1,r2,tr (T)1,rN)Pext

29、r ( G )2W2(6.13)式中，i是与格矢si相应的原子共扼动量，GhTh为数值行列式，中是多粒子体系的势函数，是h中与体积 V相应的共扼动量，是基本单元的体积，且由(h)a (bc)给出，W为相应的质量，Pex为外部压强，代表各向异性应力张量。2.2. 温度恒定体系温度恒定体系Computational Materials Science：From Basic Principle to Practical Design Methodology第 6 章分子动力学方法- 135 -接下来讨论满足温度恒定的基本方法。 恒温体系比定压体系要抽象， 简单地给出视觉化图像是困难的。势能在原子坐标

30、和动量之外引入了附加自由度f，并由下式用f把现实系统与具有能量确定的(微正则系统)假想系统联系起来。dtdtf(6.14)式中，t为现实系统的时间，t为假想系统的时间。该假想(力学)系统的哈密顿量H 可以写成(pi)2H(r1,r2,2i 12mifN),rNPf22Mg kBTex(lnf)(6.15)式中，Pi是在假想系统中原子的动量，Pf是f对应的共扼动量，g是由公式gN1给出的假想系统的自由度数目，Tex二是热浴的温度，kB为玻尔兹曼常数。显然，式错误！未错误！未找到引用源。找到引用源。中右边第一项代表粒子的动能， 第二项为粒子系统的势能， 第三项是热浴的假想动能，第四项是热浴的势能；

31、 因为在假想系统中，力学系综是微正则分布的，所以其配分函数变成具有总能量 E 的状态总数NEWNEWNEW(HE)drdrdr dp dpdp dfdPdrdrdr dp dpdp dfdP(HE)drdrdr dp dpdp dfdPdrdrdr dp dpdp dfdP12N12N12N12Nf12N12N12N12Nff(6.16)f(6.17)若将此配分函数积分变换到现实系统中，则可得到在恒定温度T 下的正则分布(也称为玻尔兹曼分布)的配分函数ZNTVNTVHBexp(k T)drdr12drNdp1dp2dpNdpNdr1dr2NdrNdp1dp2(6.18)式中，H是由下式给出的现

32、实系统的哈密顿量pi2H(r1,r2,i 12mi,rN)(6.19)H 是假想系统的守恒量，而H对真实系统来讲，由于存在热交换将不是守恒量。为了判断计算结果是否合理，只要简单地检验H 是否保持恒定(守恒)就行了。3.3. NTPNTP 系综系综到此， 已经导出了压力恒定系综和温度恒定系综的哈密顿量， 将两个系综的哈密顿量组合起来就可以得到 NTP 系综的哈密顿量，其模型如图6-7所示。NTP 系综的哈密顿量为Hi 1NitG1i2mif2(r1,r2,Pf2tr (T),rN)2W2f(6.20)1Pextr ( G )gkBTex(lnf)22MComputational Material

33、s Science：From Basic Principle to Practical Design Methodology- 136 -第 6 章分子动力学方法基本单元改变体积恒定压力原子热源维持恒定温度图6-7NTP 系统示意图若根据此哈密顿量导出正则方程式，并用速度置换动量，则可得到各变数的运动方程Nd2sifmi2ijsijmi(G1G )sifdtj 1(6.21)(6.22)(6.23)d2hWfhW2(Pexh )fdtNd2ffM2mi(hsi)T(hsi) Wtr(hTh) gkBTex fM ( )2fdti 1式中，变量上方的圆点表示对时间的微分，而x, ,分别由下列各式

34、给出4.4. NPTNPT 系综分子动力学求解系综分子动力学求解若在三维周期边界条件下求解上述方程， 则可模拟固体的结构相变、 溶解现象和晶化过程。各变量的运动方程式 (6.21) ，式(6.22) 和式(6.23) 是非线性常微分方程，不能用贝鲁勒(Berrele)法进行求解。因此，对这些方程通常用牛顿迭代方法求解。而问题的关键在于确定式(6.22)和式(6.23) 中的假想质量。所考察基本单元内的粒子体系存在着由粒子作用势以及粒子数与粒子质量共同决定的固有压力和温度等的起伏变化。 若想再现这些波动周期那样给出M和W，则由体积变化或温度变化引起的粒子系统对平衡态的偏离将迅速回复而变成稳定状态

35、。因此，若设此周期分别为tP和tT，则M和W可以表示为M6NkBT(tT2)2ij1dijrijdrij(6.24)(6.25)(6.26)NN1Nmi(hsi) (hsi)xij(rijrij)i 1i 1j ih1(6.27)(6.28)W(3LtP2)()k2式中，T为温度，kB为玻尔兹曼常数，N为原子数，L为基本单元的线尺度，k为压缩率，tP和tT分别根据对 NEV 系综的计算结果而定。【练习与思考】【练习与思考】Computational Materials Science：From Basic Principle to Practical Design Methodology第 6

36、 章分子动力学方法- 137 -6-2. 推演 NEV系综的牛顿方程，哈密顿方程，说明各个物理量的意义，在那些实际的计算体系中得到应用。6-3. 根据 NPT 系综的推导方法，对NPH 系综的哈密顿方程做简单的推导。6-4. 详细的分析不同的系综的使用范畴， 给出各个系综的使用条件， 列表进行简单的比较。6.46.4 原子势函数和分子力场构造原子势函数和分子力场构造分子动力学中存在着两个基本的概念， 一是在前面提到的系综的问题， 它所确定的是模拟系统的大小，规格和相关的模拟条件等信息；二是本节需要讲述的势函数。简单的讲，势函数就是用来描述原子和原子之间相互作用的。 随着势函数的不断发展， 现有

37、的势函数主要分成了两个大的部分。即是经典的势函数和电子理论的势函数。如图6-8所示，经典分子动力学的基本的势函数的一个简单分类。对势原子相互作用经典理论分子间作用对泛函势组合势组合泛函势无机化合物金属有机化合物半导体刚性椭球体圆柱体模型Gray-Beme势空心球模型液晶，高分子等图6-8经典的分子动力学相互作用势分类作用势的选择与动力学计算的关系极为密切， 选择不同的作用势， 体系的势能面会有不同的形状， 动力学计算所得的分子运动和分子内部运动的轨迹也会不同， 进而影响到抽样的结果和抽样结果的势能计算， 在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多

38、采用Lennard-Jones 、Morse 势等二体势模型，对于金属计算，主要采用 Morse 势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。 但是相对于二体势模型， 多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大的困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的Lennard-Jones ，Morse 等势模型的应用仍然非常广泛。6.4.16.4.1 经典理论的原子势函数经典理论的原子势函数在分子动力学的框架之下， 早期的势函数是主

39、要以描述原子和原子之间的相互作用而产生和发展起来的。早在 1903 年的时候，Mie G 就研究了两个粒子之间的相互作用势，他给出的势函数是由两项构成的， 一项代表原子之间的相互吸引， 另一项则代表着原子之间的相Computational Materials Science：From Basic Principle to Practical Design Methodology- 138 -第 6 章分子动力学方法互排斥。这就是非常著名的 Lennard-Jones 势函数的原始思考。对于其他有关势函数的发展在这里就不做过多的赘述了。从图 6-9可以清楚的知道原子之间的相互作用可以分为四类：

40、对势、 对泛函势、 组合势、组合泛函势。下面就对相关的势函数的基本的发展做一个简单的总结。表 6-1 原子中相互作用势函数势函数对势作用因素仅仅决定于原子的坐标典型代表L-J势、BMH势、GK 势对泛函势组合势组合泛函势多体相互作用共价键的相互作用EAM势有效介质理论SW 势Abell-Tersoff势CheliKowsky-Phillips势6.4.26.4.2 分子力场分子力场分子力场根据量子力学的波恩奥本海默近似， 一个分子的能量可以近似看作构成分子的各个原子的空间坐标的函数， 简单地讲就是分子的能量随分子构型的变化而变化， 而描述这种分子能量和分子结构之间关系的就是分子力场函数。 分子

41、力场函数为来自实验结果的经验公式，可以讲对分子能量的模拟比较粗糙， 但是相比于精确的量子力学从头计算方法， 分子力场方法的计算量要小数十倍， 而且在适当的范围内， 分子力场方法的计算精度与量子化学计算相差无几，因此对大分子复杂体系而言， 分子力场方法是一套行之有效的方法。 以分子力场为基础的分子力学计算方法在分子动力学、 蒙特卡罗方法、 分子对接等分子模拟方法中有着广泛的应用。1.1.力场主要讨论的对象及其构成：力场主要讨论的对象及其构成：以下参数构成一套力场函数体系需要有一套联系分子能量和构型的函数， 还需要给出各种不同原子在不同成键状况下的物理参数，比如正常的键长、 键角、二面角等，这些力

42、场参数多来自实验或者量子化学计算。键伸缩能：构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化；键角弯曲能：键角变化引起的分子能量变化；二面角扭曲能：单键旋转引起分子骨架扭曲所产生的能量变化；非键相互作用：包括范德华力、静电相互作用等与能量有关的非键相互作用；交叉能量项：上述作用之间耦合引起的能量变化。2.2.常用力场函数和分类常用力场函数和分类不同的分子力场会选取不同的函数形式来描述上述能量与体系构型之间的关系。到目前， 不同的研究小组设计了很多适用于不同体系的力场函数， 根据他们选择的函数和力场参数，可以分为以下几类：Computational Materials Science：

43、From Basic Principle to Practical Design Methodology第 6 章分子动力学方法- 139 -AMBER 力场CFF力场CHARMM力场力场CVFF 力场MMF94力场MMX力场第一代力场第二代力场Dreiding 力COMPASS力场UFF 力场ESFF力场场第三代力场图6-9分子力场的发展和分类传统力场传统力场( (第一代力场第一代力场) )a)AMBER力场：由Kollman 课题组开发的力场，是目前使用比较广泛的一种力场，适合处理生物大分子；b)CHARMM力场：由 Karplus课题组开发，对小分子体系到溶剂化的大分子体系都有很好的拟合

44、；c)CVFF 力场：CVFF 力场是一个可以用于无机体系计算的力场；d)MMX力场：MMX力场包括 MM2和 MM3 ，是目前应用最为广泛的一种力场，主要针对有机小分子。第二代力场第二代力场第二代的势能函数形式比传统力场要更加复杂，涉及的力场参数更多，计算量也更大，当然也相应地更加准确；a)CFF 力场 CFF 力场是一个力场家族，包括了CFF91 、PCFF 、CFF95 等很多力场，可以进行从有机小分子、生物大分子到分子筛等诸多体系的计算；b)COMPASS力场由 MSI 公司开发的力场，擅长进行高分子体系的计算；c)MMF94力场 Hagler开发的力场，是目前最准确的力场之一。通用力

45、场通用力场( (第一代力场第一代力场) )通用力场也叫基于规则的力场， 它所应用的力场参数是基于原子性质计算所得， 用户可以通过自主设定一系列分子作为训练集来生成合用的力场参数；a)ESFF 力场 MSI 公司开发的力场，可以进行有机、无机分子的计算；b)UFF 力场可以计算周期表上所有元素的参数；c)Dreiding 力场适用于有机小分子、大分子、主族元素的计算。3.3.分子力场的计算分子力场的计算对于分子之间的相互作用通常需要将计算体系的所有粒子纳入考虑范畴之内， 那么分子之间的相互作用可以作为长程作用来处理。 在长程相互作用中又分为周期性长程相互作用力和非周期性相互作用力。周期性的长程作

46、用力通常可以利用无网格算法，包含 Ewald 求和算法等； 而对于非周期性的相互作用通常使用快速多极子方法计算。 下面针对这两个最为常见的算法做简单的叙述。EwaldEwald 求和算法求和算法Ewald 求和算法是 Ewald 于 1921 年提出，最初是以计算离子晶体的能量。此方法选定Computational Materials Science：From Basic Principle to Practical Design Methodology- 140 -第 6 章分子动力学方法一个计算系统的盒子， 盒子中的粒子除了和盒子内部的粒子相互作用之外， 还会与其镜像系统中的粒子发生作用，

47、 并且镜像系统和盒子内部的系统完全相同。 图显示的是计算的相互作用的系统，中央的黑色的区域表示的是计算的区域，周围的镜像随着系统的距离逐渐变淡，那么这样的一个极限作用区域就像一个球，这个球称为Ewald 球。假设作用的系统的边长是L，含有n个相互作用的粒子，其镜像的作用位置可以表示为(iL, jL, kL)；其中i ,j,k都是正整数。快速多极子法快速多极子法(The Fast Multipole Method)(The Fast Multipole Method)在一个开放的系统中， 不是对所谓的晶格进行求和， 而是对要求计算的整个体系的粒子对进行求和。那么这样的一个系统不同位置粒子的静电能

48、可以简单的表述为(ri)j 1Nqjrirj(6.29)因此这个体系的静电作用和前文所讲述的 Lennard-Jones 作用项具有相同的形式。快速多极子法将需要计算的体系分为多个体积相同的格胞， 由格胞内部的所有的原子计算每一个格胞的多极(电荷，偶极，四极)。而格胞和格胞之间的相互作用是利用中央多极展开的方法进行计算。 因此中央多极展开的方法的条件是格胞之间的距离必须是大于一定距离的， 因为大于一个格胞以上的以多极展开的方式计算相互作用的力， 而格胞内部的原子对是以成对原子的作用力来计算的。【练习与思考】【练习与思考】6-5. 查找文献，推导 Lennard-Jones 势函数的公式。6-6

49、. 推导 EAM势函数的基本公式，并且实现简单的算法，编写基本设计程序。6-7. 以对势中的 Lennard-Jones势函数为例，写出对势中势函数的基本算法，并且进行基本的程序设计。6.56.5 数值方法与编程技巧数值方法与编程技巧6.5.16.5.1 边界条件边界条件随着材料的低维化和纳米化， 材料的表面和边界对于其性质的影响愈来愈大。 在处理物质表面、界面， 以及块体物质时，对其周期边界条件问题的考虑就变得非常必要了。对于周期性边界条件可以分为一维、二维及三维的情况。介绍薄膜 (可使用二维周期边界条件)和扩大到半无限的表面情况， 然后讨论液体与晶体、晶体与晶体之间形成的界面等， 最后学习

50、块状物质的三维周期边界条件。1.1. 薄膜的边界条件薄膜的边界条件( (二维周期边界条件二维周期边界条件) )对于薄膜情况，如图 6-10所示赋予其周期性边界条件。可以认为薄膜在 x-y平面内无限扩展，存在周期边界条件，而在 z方向受到限制(不赋予周期边界条件)。事实上，二维周期性边界条件也适用于在z方向不具有真正意义上的自由度的二维物质。 但就物质科学的研究范畴而言，多数情况下是指在某种衬底上生成的原子或分子的薄膜。此处讨论的边界条件，除针对薄膜以外，也可近似地用于下面将要讲到的在 z 的负方向半无限地扩展的表面问题，例如表面重构等。但是，当在z方向存在异质原子层时，原来块体物质在异质原子层