(完整版)第6章分子动力学方法汇总.docx

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

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

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

4、nwright 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.2分子动力学的计算框架6.2.1基本思想分子动力学方法是沿用牛顿运动方程或拉格朗日方程、哈密顿方程， 通过考察粒子的运动来研究多粒子系统的物理性质。在处理孤立粒子或原子团簇

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

8、给出模拟体系的统计信息，主要包括复杂的光学，电学和一些其他的物理信息。其基本思想，总结如图6-1 所示。势函数，相原子的坐标热力学性质运动方程位置互作用（牛顿方程，哈密顿方程，三维结构系综，外界拉格朗日方动力学性质条件程）原子的坐标和速度原子运动光学性质初始化条件动力学方程一次信息二次信息图 6-1分子动力学的一般思想6.2.2计算流程总体来说，分子动力学的基本计算有如下四个步骤：1.模型的设定，也就是势函数的选取势函数的研究和物理系统中对物质的描述研究息息相关。分子动力学的模拟最早使用是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是Lennard-Jones 势函数，还有嵌入原

9、子EAM 势函数，不同的物质状态描述用不同的势函数。模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。关于势函数选取将在6.4 节进一步讨论。Computational Materials Science ： From Basic Principle to Practical Design Methodology第 6 章 分子动力学方法-129-2. 给定初始条件运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如Verlet 算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。一般意义上讲系统的初始条件不可能知道

10、，实际上也不需要精确选择代求系统的初始条件，因为当模拟时间足够长时，系统就会忘掉初始条件(对于无记忆的体系而言)。当然，合理的初始条件可以加快系统趋于平衡的时间和步程，获得好的精度。常用的初始条件可以选择为，令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零； 令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。设定坐标和相互作用t+t计算作用在原子上的力Vi (tt / 2)Vi (tt / 2)tFi (t ) / miRi (tt)Ri ( t)t Vi (tt / 2)计算物理量并

11、对其结果进行统计tt maxNY结束图 6-2 分子动力学计算流程3. 趋于平衡计算在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。 但这样计算出的系统是不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。为使得系统平衡， 模拟中设计一个趋衡过程，即在这个过程中，增加或者从系统中移出能量，直到持续给出确定的能量值，称此时系统已经达到平衡，达到平衡的时间称为弛豫时间。分子动力学中， 积分步程的大小选择十分重要，这决定了模拟所需要的时间。为了减小误差，步长要小，但过小系统模拟的弛豫时间就长。因此根据经验选择适当的步长。如，对一个具有几百个氩气分子的体系，采用Lennard

12、-Jones 势函数，发现取 h 为 0.01量级，可以得到很好的结果。这里选择的h 是没有量纲的，实际上这样选择的h 对应的时间在 10-14s 的量级。如果模拟 1000 步，系统达到平衡，弛豫时间只有10-11s。4. 宏观物理量的计算Computational Materials Science ： From Basic Principle to Practical Design Methodology-130-第 6 章 分子动力学方法实际计算宏观的物理量往往是在模拟的最后阶段进行的，它是沿相空间轨迹求平均来计算得到的， 根据各态历经假说，时间平均代替系综平均。关于各态历经假说将在第

13、8 章的蒙特卡洛方法中深入讨论。6.2.3经典分子动力学中近似处理实际的计算体系由于其外在条件的复杂性，往往需要对计算的系统进行一定的近似处理，在经典分子动力学计算中，引进以下近似处理：? 经典粒子相互作用，不考虑电子相互作用量子效应。? 力的作用形式，由参数可调的相互作用势函数决定，并经过实验测定来验证。? 模拟体系与实际体系相差较大，一般需要采用周期边界来扩展计算体系。时间平均是在有限时间内完成。【练习与思考】6-1.查找文献，根据上述的分子动力学的处理流程图，编写实现分子动力学的简单程序，可参考 Daan F 和 Berend S 编著的分子模拟一书。6.3分子动力学的系综系综 (ens

14、emble)是指在一定的宏观条件下(约束条件 )，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统集合，或称为统计系综。系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念； 系综是统计理论的一种表述方式， 系综理论使统计物理成为普遍的微观统计理论；系综并不是实际的物体，构成系综的系统才是实际物体。约束条件是由一组外加宏观参量来表示。在平衡统计力学范畴下，可以用来处理稳定系综。6.3.1常用系综分类1. 正则系综 (canonical ensemble)全称应为 “宏观正则系综”，简写为 NVT ，即表示具有确定的粒子数(N) 、体积 (V) 、温度(T) 。正则系综是分

15、子动力学方法模拟处理的典型代表。假定N 个粒子处在体积为V 的盒子内，将其埋入温度恒为T 的热浴中。此时，总能量(E)和系统压强 (P)可能在某一平均值附近起伏变化。 平衡体系代表封闭系统，与大热源热接触平衡的恒温系统。正则系综的特征函数是亥姆霍兹自由能F(N,V,T) 。2. 微正则系综 (micro-canonical ensemble)简写为 NVE ，即表示具有确定的粒子数 (N) 、体积 (V) 、总能量 (E)。微正则系综广泛被应用在分子动力学模拟中。 假定 N 个粒子处在体积为 V 的盒子内， 并固定总能量 (E) 。此时，系综的温度 (T) 和系统压强 (P)可能在某一平均值附

16、近起伏变化。平衡体系为孤立系统，与外界即无能量交换，也无粒子交换。微正则系综的特征函数是熵S(N,V,E)。3. 等温等压 (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

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

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

19、(t ) ，为求出该粒子 (i) 在时刻 tt 的位置，经常使用的方法是Verlet 算法 (对于该算法的具体细节在下文将有详细的推导)。根据系统的动能及统计物理学中的能量均分定理，则系统的平衡温度可表达为21pi2(6.1)T2 imikB N f根据此方法， 在晶体中粒子的配置与排布、晶形变化等结构相变也可以用计算机模拟进行研究。根据物理化学的基本气体的基本方程可以知道，体系的NEV 系综下的计算体系主要变化的量有两个压力和温度。下文就从恒温和恒压两个条件下对NEV 系综的基本方法做简单的推导。1. 恒温方法Nse S 和 Hoover W G 引入与热源相关的参数，尝试构造恒温状态， 亦

20、即具有确定N ，V ， T 值的系统，对此可设想为与大热源接触而达到平衡的系统。由于热源很大，交换能量不会改变热源的温度。在两者建立平衡后，系统将与热源具有相同的温度，系统与热源合起来构成一个复合系统，如图6-3 所示。这样，系统中微观粒子的动力学方程：dqipi(6.2)dtmidpipi(6.3)dtdqidpi23 N kB Tex 2(6.4)dti 2mi2QComputational Materials Science ： From Basic Principle to Practical Design Methodology-132-第 6 章分子动力学方法式中为系统的势函数，

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

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

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

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

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

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

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

28、6.13)(r1 ,r2 ,L L , rN )2Wtr ( G )i 12mi2式中，i 是与格矢 si相应的原子共扼动量，GhT h 为数值行列式，中是多粒子体系的势函数，是 h 中 与 体 积 V 相 应 的 共 扼 动 量 ，是基本单元的体积，且由(h) a(b c) 给出， W 为相应的质量， Pex 为外部压强，代表各向异性应力张量。2. 温度恒定体系Computational Materials Science ： From Basic Principle to Practical Design Methodology第 6 章 分子动力学方法-135-接下来讨论满足温度恒定的基

29、本方法。恒温体系比定压体系要抽象， 简单地给出视觉化图像是困难的。势能在原子坐标和动量之外引入了附加自由度f ，并由下式用f 把现实系统与具有能量确定的 (微正则系统 )假想系统联系起来。dt(6.14)dtf式中， t为现实系统的时间，t 为假想系统的时间。该假想(力学 ) 系统的哈密顿量H 可以写成N( pi )22H( r1 , r2 ,L , rN )PfTex (lnf )(6.15)2m f2g kBi 12Mi式中， Pi是在假想系统中原子的动量，Pf 是 f 对应的共扼动量，g 是由公式 gN 1 给出的假想系统的自由度数目，Tex 二是热浴的温度，kB 为玻尔兹曼常数。显然，

30、式错误 !未找到引用源。 中右边第一项代表粒子的动能，第二项为粒子系统的势能，第三项是热浴的假想动能， 第四项是热浴的势能；因为在假想系统中，力学系综是微正则分布的，所以其配分函数 NEW 变成具有总能量 E 的状态总数L(HE)dr1 dr2 LdrN dp1dp2 LdpN dfdPf(6.16)NEWLdr1 dr2 L drN dp1dp2 L dpN dfdPfL(HE)dr1 dr2 LdrN dp1dp2 LdpN dfdPf(6.17)NEWLdr1 dr2 L drN dp1dp2 L dpN dfdPf若将此配分函数积分变换到现实系统中，则可得到在恒定温度T 下的正则分布玻

31、尔兹曼分布 )的配分函数Z NTVL exp(H )dr dr2L drNdp dp2LdpN11NTVkBTdr1dr2 L dr Ndp1dp 2 L dp NL式中， H 是由下式给出的现实系统的哈密顿量(也称为(6.18)N2Hpi(r1, r2 ,L , rN )(6.19)2mii 1H 是假想系统的守恒量，而H 对真实系统来讲，由于存在热交换将不是守恒量。为了判断计算结果是否合理，只要简单地检验H 是否保持恒定(守恒 )就行了。3. NTP 系综到此，已经导出了压力恒定系综和温度恒定系综的哈密顿量， 将两个系综的哈密顿量组合起来就可以得到 NTP 系综的哈密顿量，其模型如图 6-

32、7 所示。NTP 系综的哈密顿量为Nt1tr (T)Hi Gi(r1 , r2,L , rN )2mif 22W2fi 11Pf2(6.20)PexgkB Tex (ln f )tr ( G)2M2Computational Materials Science ： From Basic Principle to Practical Design Methodology- 136 -第 6 章分子动力学方法基本单元改变体积恒定压力原子热源维持恒定温度图6-7NTP系统示意图若根据此哈密顿量导出正则方程式，并用速度置换动量，则可得到各变数的运动方程2Nf&1 & &d siij sijmi (mi

33、2fG G )sidtj 1d2h&(PexhWfhW2)dtfd 2 fNT&iTf& 2M&i& &gkBTex fM ( ) mi (hs )(hs )Wtr (hh)dt 2i 1f式中，变量上方的圆点表示对时间的微分，而x, 分别由下列各式给出ij1 dijrijdrij1 N&)&)NNrij)mi (hsi(hsixij (riji1i1j ih 14. NPT 系综分子动力学求解(6.21)(6.22)(6.23)(6.24)(6.25)(6.26)若在三维周期边界条件下求解上述方程，则可模拟固体的结构相变、溶解现象和晶化过程。各变量的运动方程式(6.21) ，式 (6.22)

34、 和式 (6.23) 是非线性常微分方程，不能用贝鲁勒(Berrele) 法进行求解。因此，对这些方程通常用牛顿迭代方法求解。而问题的关键在于确定式 (6.22) 和式 (6.23) 中的假想质量。所考察基本单元内的粒子体系存在着由粒子作用势以及粒子数与粒子质量共同决定的固有压力和温度等的起伏变化。 若想再现这些波动周期那样给出 M 和 W ，则由体积变化或温度变化引起的粒子系统对平衡态的偏离将迅速回复而变成稳定状态。因此，若设此周期分别为tP 和 tT ，则 M 和W 可以表示为M6Nk BT ( tT)2(6.27)2W(3L )( tP )2(6.28)k 2式中， T 为温度， kB

35、为玻尔兹曼常数， N 为原子数， L 为基本单元的线尺度，k 为压缩率， tP 和 tT 分别根据对 NEV系综的计算结果而定。【练习与思考】Computational Materials Science ： From Basic Principle to Practical Design Methodology第 6 章 分子动力学方法-137-6-2.推演 NEV 系综的牛顿方程，哈密顿方程，说明各个物理量的意义，在那些实际的计算体系中得到应用。6-3.根据 NPT 系综的推导方法，对NPH 系综的哈密顿方程做简单的推导。6-4.详细的分析不同的系综的使用范畴，给出各个系综的使用条件， 列

36、表进行简单的比较。6.4原子势函数和分子力场构造分子动力学中存在着两个基本的概念， 一是在前面提到的系综的问题， 它所确定的是模拟系统的大小， 规格和相关的模拟条件等信息；二是本节需要讲述的势函数。简单的讲，势函数就是用来描述原子和原子之间相互作用的。随着势函数的不断发展，现有的势函数主要分成了两个大的部分。即是经典的势函数和电子理论的势函数。如图6-8 所示，经典分子动力学的基本的势函数的一个简单分类。对势无机化合物原子对泛函势金属相互作组合势有机化合物用经典组合泛函势半导体理论分刚性椭球体子圆柱体模型液晶，间Gray-Beme势高分子等作空心球模型用图 6-8 经典的分子动力学相互作用势分

37、类作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状， 动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、Morse 势等二体势模型，对于金属计算，主要采用Morse 势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM 等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是相对于二体势模型，多体势往往缺乏明确的表达式，参量很

38、多，模拟收敛速度很慢，给应用带来很大的困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的Lennard-Jones， Morse 等势模型的应用仍然非常广泛。6.4.1经典理论的原子势函数在分子动力学的框架之下， 早期的势函数是主要以描述原子和原子之间的相互作用而产生和发展起来的。早在 1903 年的时候， Mie G 就研究了两个粒子之间的相互作用势，他给出的势函数是由两项构成的， 一项代表原子之间的相互吸引， 另一项则代表着原子之间的相Computational Materials Science ： From Basic Principle to Practical Design

39、 Methodology- 138 -第 6 章分子动力学方法互排斥。这就是非常著名的 Lennard-Jones 势函数的原始思考。对于其他有关势函数的发展在这里就不做过多的赘述了。从图 6-9 可以清楚的知道原子之间的相互作用可以分为四类：对势、对泛函势、 组合势、组合泛函势。下面就对相关的势函数的基本的发展做一个简单的总结。表 6-1 原子中相互作用势函数势函数作用因素典型代表L-J 势、对势仅仅决定于原子的坐标BMH 势、GK 势对泛函势多体相互作用EAM 势有效介质理论组合势共价键的相互作用SW 势组合泛函势Abell-Tersoff 势CheliKowsky-Phillips 势6

40、.4.2分子力场分子力场根据量子力学的波恩奥本海默近似，一个分子的能量可以近似看作构成分子的各个原子的空间坐标的函数，简单地讲就是分子的能量随分子构型的变化而变化，而描述这种分子能量和分子结构之间关系的就是分子力场函数。分子力场函数为来自实验结果的经验公式， 可以讲对分子能量的模拟比较粗糙，但是相比于精确的量子力学从头计算方法，分子力场方法的计算量要小数十倍，而且在适当的范围内，分子力场方法的计算精度与量子化学计算相差无几，因此对大分子复杂体系而言，分子力场方法是一套行之有效的方法。以分子力场为基础的分子力学计算方法在分子动力学、蒙特卡罗方法、 分子对接等分子模拟方法中有着广泛的应用。1.力场

41、主要讨论的对象及其构成：以下参数构成一套力场函数体系需要有一套联系分子能量和构型的函数，还需要给出各种不同原子在不同成键状况下的物理参数，比如正常的键长、键角、二面角等，这些力场参数多来自实验或者量子化学计算。? 键伸缩能：构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化；? 键角弯曲能：键角变化引起的分子能量变化；? 二面角扭曲能：单键旋转引起分子骨架扭曲所产生的能量变化；? 非键相互作用：包括范德华力、静电相互作用等与能量有关的非键相互作用；? 交叉能量项：上述作用之间耦合引起的能量变化。2. 常用力场函数和分类不同的分子力场会选取不同的函数形式来描述上述能量与体系构型之间的关系

42、。到目前，不同的研究小组设计了很多适用于不同体系的力场函数， 根据他们选择的函数和力场参数，可以分为以下几类：Computational Materials Science ： From Basic Principle to Practical Design Methodology第 6 章 分子动力学方法-139-AMBER 力场ESFF力场CFF力场CHARMM 力场力场COMPASS力场UFF力场CVFF力场MMF94 力场Dreiding力MMX 力场场第一代力场第二代力场第三代力场图 6-9 分子力场的发展和分类? 传统力场 (第一代力场 )a) AMBER 力场：由 Kollman

43、 课题组开发的力场，是目前使用比较广泛的一种力场，适合处理生物大分子；b) CHARMM 力场：由 Karplus 课题组开发，对小分子体系到溶剂化的大分子体系都有很好的拟合；c) CVFF 力场： CVFF 力场是一个可以用于无机体系计算的力场；d) MMX 力场： MMX 力场包括 MM2 和 MM3 ，是目前应用最为广泛的一种力场，主要针对有机小分子。? 第二代力场第二代的势能函数形式比传统力场要更加复杂，涉及的力场参数更多，计算量也更大，当然也相应地更加准确；a)CFF 力场 CFF 力场是一个力场家族，包括了CFF91、 PCFF、CFF95 等很多力场，可以进行从有机小分子、生物大

44、分子到分子筛等诸多体系的计算；b) COMPASS 力场由 MSI 公司开发的力场，擅长进行高分子体系的计算；c) MMF94 力场 Hagler 开发的力场，是目前最准确的力场之一。? 通用力场 (第一代力场 )通用力场也叫基于规则的力场， 它所应用的力场参数是基于原子性质计算所得，用户可以通过自主设定一系列分子作为训练集来生成合用的力场参数；a) ESFF 力场 MSI 公司开发的力场，可以进行有机、无机分子的计算；b) UFF 力场可以计算周期表上所有元素的参数；c) Dreiding 力场适用于有机小分子、大分子、主族元素的计算。3. 分子力场的计算对于分子之间的相互作用通常需要将计算体系的所有粒子纳入考虑范畴之内，那么分子之间的相互作用可以作为长程作用来处理。在长程相互作用中又分为周期性长程相互作用力和非周期性相互作用力。周期性的长程作用力通常可以利用无网格算法，包含Ewald 求和算法等； 而对于非周期性的相互作用通常使用快速多极子方法计算。下面针对这两个最为常见的算法做简单的叙述。? Ewald 求和算法Ewald 求和算法是 Ewald 于 1921 年提出，最初是以计算离子晶体的能量。此方法选定Computational Materials Science ： From Basic Principle to Practical