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

《第6章分子动力学方法汇总.pdf》由会员分享，可在线阅读，更多相关《第6章分子动力学方法汇总.pdf（28页珍藏版）》请在三一文库上搜索。

1、第 6 章分子动力学方法 Computational Materials Science ：From Basic Principle to Practical Design Methodology - 127 - 第6章分子动力学方法 经典分子动力学方法无疑是材料，尤其是大分子体系和大体系模拟有效的方法之一。分 子动力学可以用于NPT，NVE ， NVT 等不同系综的计算，是一种基于牛顿力学确定论的热 力学计算方法。 与蒙特卡罗法相比在宏观性质计算上具有更高的准确度和有效性，可以广泛 应用于物理，化学，生物，材料，医学等各个领域。本章在介绍分子动力学的基本概念的基 础上， 简单介绍了分子动力学

2、的基本思想，势函数分类和基本方程。然后介绍了分子动力学 的常用系综和典型的NPT，NVE ，NVT 系综基本方程。结合材料建模中的基本简化方法和 技巧， 阐述了边界条件和时间积分的数值处理技巧。最后， 利用统计力学的基本概念给出分 子动力学的计算信息的解析方式。并且结合Materials Explore 软件计算分析了CNT 的几何 结构稳定性。 6.1引言 分子动力学方法(Molecular Dynamics, MD) 方法是一种按该体系内部的内禀动力学规律 来计算并确定位形的变化的确定性模拟方法。首先需要在给定的外界条件下建立一组粒子的 运动方程， 然后通过直接对系统中的一个个粒子运动方程

3、进行数值求解，得到每个时刻各个 分子的坐标与动量，即在相空间的运动轨迹，再利用统计力学方法得到多体系统的静态和动 态特性，从而获得系统的宏观性质。可以看出，分子动力学方法中不存在任何随机因素，这 个也是分子动力学方法和后文要提到的蒙特卡洛方法的区别之一。在分子动力学方法的处理 过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中， 每个分子都各自服从经典的牛顿力学定律(或者是拉格朗日方程)。每个分子运动的内禀动力 学是用理论力学上的哈密顿量或者拉格朗日函数来描述，也可以直接用牛顿运动方程来描 述。确定性方法是实现玻尔兹曼的统计力学途径。这种方法可以处理与时间有关的过

4、程，因 而可以处理非平衡态问题。但是分子动力学方法的计算机程序相对蒙特卡罗较复杂，其计算 成本较高。 分子动力学方法发展历史改革经历了近60 年，分子动力学方法是20 世纪 50 年代后期 由 Alder B J 和 Wainwright T E创造发展的。 Alder 和 Wainwright 在 1957 年利用分子动力学 模拟， 发现了 “ 刚性球组成的集合系统会发生由其液相到结晶相的相转变” ，后来人们称这种 相变为 Alder 相变。其结果表明， 不具有引力的粒子系统也具有凝聚态。到 20 世纪 70 年代， 产生了刚性体系的动力学方法，成功地被应用于水和氮等分子性溶液体系的处理；1

5、972 年， Less A W 和 Edwards S F 等发展了该方法并扩展到了存在速度梯度，即处于非平衡状态的 系统。 之后， 此方法被Ginan M J 等推广到了具有温度梯度的非平衡系统，从而构造并形成 了所谓的非平衡分子动力学方法体系。进人 20 世纪 80 年代之后， 出现了在分子内部对一部 分自由度施加约束条件的新的分子动力学方法，从而使分子动力学方法可适用于类似蛋白 质等生物大分子的解析与设计。分子动力学方法真正作为材料科学领域的一个重要研究方 法，开始于由Andersen, Parrinello 和 Rahman 等创立恒压分子动力学方法和N se 等完成恒 第 6 章分子

6、动力学方法 Computational Materials Science ：From Basic Principle to Practical Design Methodology - 128 - 温分子动力学方法的建立及在应用方面的成功。后来，针对势函数模型化比较困难的半导体 和金属等， 1985 年人们又提出了将基于密度泛函理论的电子论和分子动力学方法有机统一 起来的所谓Car-Parrinello 方法，亦即第一性原理分子动力学方法。这样，分子动力学的方 法进一步得到发展和完善，它不仅可以处理半导体和金属的间题，同时还可应用于处理有机 物和化学反应。关于Car-Parrinello 分

7、子动力学方法将在第7 章重点学习讨论。本章重点讨 论经典分子动力学方法的基本原理和计算方法。 6.2分子动力学的计算框架 6.2.1基本思想 分子动力学方法是沿用牛顿运动方程或拉格朗日方程、哈密顿方程， 通过考察粒子的运 动来研究多粒子系统的物理性质。在处理孤立粒子或原子团簇(不考虑与其他粒子的相互作 用)时，单纯地对牛顿运动方程进行积分求解就行了；而在处理凝聚系统时，可以认为将所 考虑的 N 个粒子放入具有一定体积的容器中，在这样组成的封闭系中，建立直角坐标(x, y, z)，每个粒子的位置就由三个坐标分量决定，通过求解3N 个联立方程组就可以得到保守系 的总能量。 分子动力学的基本思想在于

8、通过设定原子之间的相互作用(势函数 )和相关的系综 (亦即作用对象和条件)，确定其基本的模拟范畴。在给定的牛顿方程、拉格朗日方程或者是 哈密顿方程进行时间的迭代，在达到指定的力学收敛条件之后得到一个最终的位置坐标，然 后通过相关的位形信息和运动的速度和加速度等信息通过热力统计学的方法，给出模拟体系 的统计信息，主要包括复杂的光学，电学和一些其他的物理信息。其基本思想，总结如图 6-1 所示。 势函数，相 互作用 系综，外界 条件 运动方程 （牛顿方程， 哈密顿方程， 拉格朗日方 程） 原子的坐标 位置 原子的坐标 和速度 热力学性质 动力学性质 光学性质 原子运动 三维结构 一次信息二次信息

9、动力学方程初始化条件 图 6-1 分子动力学的一般思想 6.2.2计算流程 总体来说，分子动力学的基本计算有如下四个步骤： 1.模型的设定，也就是势函数的选取 势函数的研究和物理系统中对物质的描述研究息息相关。分子动力学的模拟最早使用是 硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是Lennard-Jones 势函数， 还有嵌入原子EAM 势函数，不同的物质状态描述用不同的势函数。模型势函数一旦确定， 就可以根据物理学规律求得模拟中的守恒量。关于势函数选取将在6.4 节进一步讨论。 第 6 章分子动力学方法 Computational Materials Science ：Fro

10、m Basic Principle to Practical Design Methodology - 129 - 2.给定初始条件 运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如 Verlet 算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或 者一组零时刻的速度值。 一般意义上讲系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条 件，因为当模拟时间足够长时，系统就会忘掉初始条件(对于无记忆的体系而言)。当然，合 理的初始条件可以加快系统趋于平衡的时间和步程，获得好的精度。 常用的初始条件可以选择为，令初始位置在差分划分网格

11、的格子上，初始速度则从玻尔 兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零； 令 初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。 设定坐标和 相互作用 t+ t 计算作用在原子上的力 计算物理量并对其结果 进行统计 结束 tt max Y N (/ 2)(/ 2)( ) / ()( )(/ 2) iiii iii V ttVtttF tm R ttRtt V tt 图6-2 分子动力学计算流程 3.趋于平衡计算 在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出 的系统是不会具有所要求的系统能量，并且这个

12、状态本身也还不是一个平衡态。为使得系统 平衡， 模拟中设计一个趋衡过程，即在这个过程中，增加或者从系统中移出能量，直到持续 给出确定的能量值，称此时系统已经达到平衡，达到平衡的时间称为弛豫时间。 分子动力学中， 积分步程的大小选择十分重要，这决定了模拟所需要的时间。为了减小 误差，步长要小，但过小系统模拟的弛豫时间就长。因此根据经验选择适当的步长。如，对 一个具有几百个氩气分子的体系，采用Lennard-Jones 势函数，发现取h 为 0.01 量级，可以 得到很好的结果。这里选择的h 是没有量纲的，实际上这样选择的h 对应的时间在10-14s 的 量级。如果模拟1000 步，系统达到平衡，

13、弛豫时间只有10-11s。 4.宏观物理量的计算 第 6 章分子动力学方法 Computational Materials Science ：From Basic Principle to Practical Design Methodology - 130 - 实际计算宏观的物理量往往是在模拟的最后阶段进行的，它是沿相空间轨迹求平均来计 算得到的， 根据各态历经假说，时间平均代替系综平均。关于各态历经假说将在第8 章的蒙 特卡洛方法中深入讨论。 6.2.3经典分子动力学中近似处理 实际的计算体系由于其外在条件的复杂性，往往需要对计算的系统进行一定的近似处 理，在经典分子动力学计算中，引进以下

14、近似处理： 经典粒子相互作用，不考虑电子相互作用量子效应。 力的作用形式，由参数可调的相互作用势函数决定，并经过实验测定来验证。 模拟体系与实际体系相差较大，一般需要采用周期边界来扩展计算体系。 时间平均是在有限时间内完成。 【练习与思考】 6-1.查找文献，根据上述的分子动力学的处理流程图，编写实现分子动力学的简单程序， 可参考 Daan F 和 Berend S 编著的分子模拟一书。 6.3分子动力学的系综 系综 (ensemble)是指在一定的宏观条件下(约束条件 )，大量性质和结构完全相同的、处 于各种运动状态的、各自独立的系统集合，或称为统计系综。系综是用统计方法描述热力学 系统的统

15、计规律性时引入的一个基本概念；系综是统计理论的一种表述方式，系综理论使统 计物理成为普遍的微观统计理论；系综并不是实际的物体，构成系综的系统才是实际物体。 约束条件是由一组外加宏观参量来表示。在平衡统计力学范畴下，可以用来处理稳定系综。 6.3.1常用系综分类 1.正则系综 (canonical ensemble) 全称应为 “ 宏观正则系综” ，简写为 NVT ，即表示具有确定的粒子数(N)、体积 (V) 、温度 (T)。正则系综是分子动力学方法模拟处理的典型代表。假定N 个粒子处在体积为V 的盒子 内，将其埋入温度恒为T 的热浴中。此时，总能量(E)和系统压强 (P)可能在某一平均值附近

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

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

18、子数(N)、压强 (P)、焓 (H) 。由于 HEPV，故在 该系综下进行模拟时要保持压力与焓值为固定，其调节技术的实现也有一定的难度，这种系 综在实际的分子动力学模拟中已经很少遇到了。 5.巨正则系综 (grand canonical ensemble) 简写为 VT ，即表示具有确定的粒体积(V) 、温度 (T)和化学势 ( )。巨正则系综通常是蒙 特卡罗模拟的对象和手段。此时、系统能量(E)、压强 (P)和粒子数 (N)会在某一平均值附近有 一个起伏。体系是一个开放系统，与大热源大粒子源热接触平衡而具有恒定的温度(T)。特 征函数是马休 (Massieu)函数 J( ，V，T)。 对于不

19、同的系综，由于其选择的运动方程的求解方法亦有不同。针对不同的计算体系具 体的细节比较复杂， 下文将以NEV 和 NTP 这两类简单的系综运用两类不同的条件方法下求 解的不同方程做简要的介绍。 6.3.2NEV 系综基本方程 NEV 系综是一类较为常见的系综。对于该系综，经典分子动力学方法一般采用差分近 似法求解牛顿运动方程，并追踪系统的时间变化。考虑由N 个粒子构成的系统，假设第i 个粒子在时刻( ) i x t位置，速度为( ) i v t 受到的作用力为( ) i F t，为求出该粒子(i)在时刻 tt 的 位置，经常使用的方法是Verlet 算法 (对于该算法的具体细节在下文将有详细的推

20、导)。 根据系统的动能及统计物理学中的能量均分定理，则系统的平衡温度可表达为 2 21 2 i i Bfi p T kNm (6.1) 根据此方法， 在晶体中粒子的配置与排布、晶形变化等结构相变也可以用计算机模拟进 行研究。 根据物理化学的基本气体的基本方程可以知道，体系的NEV 系综下的计算体系主要变 化的量有两个压力和温度。下文就从恒温和恒压两个条件下对NEV 系综的基本方法做简单 的推导。 1.恒温方法 N se S和 Hoover W G 引入与热源相关的参数，尝试构造恒温状态，亦即具有确定N， V，T 值的系统，对此可设想为与大热源接触而达到平衡的系统。由于热源很大，交换能量 不会改

21、变热源的温度。在两者建立平衡后，系统将与热源具有相同的温度，系统与热源合起 来构成一个复合系统，如图6-3 所示。 这样，系统中微观粒子的动力学方程： ii i dqp dtm (6.2) i i i dp p dtdq (6.3) 2 32 22 i Bex i i pd NkT dtmQ (6.4) 第 6 章分子动力学方法 Computational Materials Science ：From Basic Principle to Practical Design Methodology - 132 - 式中为系统的势函数，Q 为有效质量，表示热力学摩尔系数。式(6.2)至式 (6.

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

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

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

25、 V W (6.7) 212 33 0 () 2 q i i i p HVV m (6.8) 式中， ex P 是外压； v p 是体积V对应的共轭动量；W是决定于体积变化速度的因子。若将由 上述哈密顿量导出的运动方程用粒子的实际坐标直接写出来，则有 3 iii i dqpqdV dtmdtV (6.9) () 3 ii i dppdV dtdtVq (6.10) 2 2 2 () () 3 i i i i ii ex p q mqd V WP Vdt (6.11) 式(6.9)中右边第一项是根据维里(Virial) 定理求出的系统的瞬间压P。在此处理中，特别 重要的是式 (6.9)右边第一项

26、与统计力学压力表达式相同。 6.3.3NTP 系综质点系的基本方程 前面已讨论了关于NEV 系综的质点系方程。鉴于NTV 和 NHP 系综的基本方程包含于 NTP 系综的基本方程之中，在此只对NTP 系综的基本方程作一些说明。NTP 系综是在宏观 上具有确定的粒子数N、温度 T 和压力 P 的系综。在此系综中，T 和 P 可作为外部参数给 出。那么，怎样才能满足N，T，P 恒定的条件呢？显然，粒子数N 保持一定是最简单的， 因为只要使所考察的基本单元内的粒子数恒定就行了。 1.压力恒定体系 为简单起见，以由圆柱容器内充满气体和调节活塞构成的系统为例，如图6-5 所示。设 活塞的重量为W，则在平

27、衡状态下，由活塞的重量产生的外部压力( ex P )与由气体产生的内 压平衡。如果外压 ex P 变大，为了保持平衡，则气体部分的体积将收缩，从而内压变大。 第 6 章分子动力学方法 Computational Materials Science ：From Basic Principle to Practical Design Methodology - 134 - Pex 图6-5活塞与气体容器构成的示意图 因此， 把体积作为系统的动力学参数建立其运动方程，就好像维持压力恒定。其方程可 写为 2 2 () ex d V WPV dt (6.12) 式中，V为体积，是由维里 (Virial)

28、 定理求出的内压，W是对应于V的假想质量。图6-5 为活塞与气体容器构成的示意图。将上述想法推广到一般情况下的晶格系统，如图 6-6 所示， 对一般的晶体，其基本单元可由三个平移矢量(,) xyz aaaa,(,) xyz bbbb, (,) xyz cccc决 定其形状和大小。将, ,a b c的三个分量写成矩阵h，该矩阵规定了晶体基本单元的形状和大 小，这样h也成为原子的实际坐标 i r 和晶格矢量 i s 之间的变换矩阵，亦即 xyz xyz xyz aaa hbbb ccc 错误！未找到引用源。 ii rhs错误！未找到引用源。 c a b 图6-6基于单元形状示意图 在式 (0.2)

29、中，体积V作为动力学参数，而在一般晶体系统中，h应该作为动力学参数。 若给出此动力学系统的哈密顿量，则有 1 12 1 ()1 (,)() 222 t T N ii Nex i i Gtr Hr rrPtrG mW (6.13) 式中， i是与格矢i s 相应的原子共扼动量， T Gh h为数值行列式，中是多粒子体系 的 势 函 数 ，是h中 与 体 积V相 应 的 共 扼 动 量 ，是 基 本 单 元 的 体 积 ， 且 由 ( )()habc给出，W为相应的质量， ex P 为外部压强，代表各向异性应力张量。 2.温度恒定体系 第 6 章分子动力学方法 Computational Mate

30、rials Science ：From Basic Principle to Practical Design Methodology - 135 - 接下来讨论满足温度恒定的基本方法。恒温体系比定压体系要抽象，简单地给出视觉化 图像是困难的。势能在原子坐标和动量之外引入了附加自由度f ，并由下式用f 把现实系 统与具有能量确定的(微正则系统 )假想系统联系起来。 dt dt f (6.14) 式中，t为现实系统的时间， t为假想系统的时间。该假想(力学 )系统的哈密顿量H 可 以写成 2 2 12 2 1 () (,)(ln) 22 N f i NBex i i Pp Hr rrg kTf

31、Mm f (6.15) 式中， i P 是在假想系统中原子的动量， f P 是 f 对应的共扼动量，g是由公式1gN给 出的假想系统的自由度数目， ex T 二是热浴的温度， B k 为玻尔兹曼常数。显然，式错误！未 找到引用源。 中右边第一项代表粒子的动能，第二项为粒子系统的势能，第三项是热浴的假 想动能，第四项是热浴的势能；因为在假想系统中，力学系综是微正则分布的，所以其配分 函数 NEW变成具有总能量 E 的状态总数 1212 1212 () NNf NEW NNf HE dr drdr dpdpdp dfdP dr drdr dp dpdp dfdP (6.16) 1212 1212

32、() NNf NEW NNf HE dr drdr dpdpdp dfdP dr drdr dp dpdp dfdP (6.17) 若将此配分函数积分变换到现实系统中，则可得到在恒定温度T 下的正则分布 (也称为 玻尔兹曼分布 )的配分函数 NTV Z 1212 1212 exp() NN B NTV NN H dr drdr dp dpdp k T dr drdrdp dpdp (6.18) 式中， H 是由下式给出的现实系统的哈密顿量 2 12 1 ( ,) 2 N i N i i p Hrrr m (6.19) H 是假想系统的守恒量，而H 对真实系统来讲，由于存在热交换将不是守恒量。为

33、了 判断计算结果是否合理，只要简单地检验H 是否保持恒定 (守恒 )就行了。 3.NTP 系综 到此，已经导出了压力恒定系综和温度恒定系综的哈密顿量，将两个系综的哈密顿量组 合起来就可以得到NTP 系综的哈密顿量，其模型如图6-7 所示。 NTP 系综的哈密顿量为 1 12 22 1 2 () ( ,) 22 1 ()(ln) 22 t T N ii N i i f exBex Gtr Hr rr m fWf P PtrGgkTf M (6.20) 第 6 章分子动力学方法 Computational Materials Science ：From Basic Principle to Pra

34、ctical Design Methodology - 136 - 改变体积 恒定压力 热源维持恒 定温度 基本单元 原子 图6-7NTP系统示意图 若根据此哈密顿量导出正则方程式，并用速度置换动量，则可得到各变数的运动方程 2 1 2 1 () N i iijijii j d sf msmGG s fdt (6.21) 2 2 () ex d hWfh WPh fdt (6.22) 2 2 2 1 () ()()() N TT iiiBex i d ff Mm hshsWtr h hgk TfM fdt (6.23) 式中，变量上方的圆点表示对时间的微分，而,x分别由下列各式给出 1ij i

35、j ijij d rdr (6.24) 11 1 ()()() NNN iiiijijij iiji m hshsxrr(6.25) 1 h(6.26) 4.NPT 系综分子动力学求解 若在三维周期边界条件下求解上述方程，则可模拟固体的结构相变、溶解现象和晶化过 程。各变量的运动方程式(6.21)，式 (6.22)和式 (6.23)是非线性常微分方程，不能用贝鲁勒 (Berrele)法进行求解。因此，对这些方程通常用牛顿迭代方法求解。而问题的关键在于确定 式(6.22) 和式 (6.23)中的假想质量。所考察基本单元内的粒子体系存在着由粒子作用势以及 粒子数与粒子质量共同决定的固有压力和温度等

36、的起伏变化。若想再现这些波动周期那样给 出 M 和W，则由体积变化或温度变化引起的粒子系统对平衡态的偏离将迅速回复而变成稳 定状态。因此，若设此周期分别为 P t 和 T t ，则 M 和W可以表示为 2 6() 2 T B t MNk T(6.27) 23 ()() 2 PtL W k (6.28) 式中，T为温度， B k 为玻尔兹曼常数， N为原子数，L为基本单元的线尺度，k为压 缩率， P t 和 T t 分别根据对NEV 系综的计算结果而定。 【练习与思考】 第 6 章分子动力学方法 Computational Materials Science ：From Basic Princi

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

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

39、计算宏观体积和微观成分关系的时候主要采用刚球模型的二 体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、Morse 势等二体势模型，对于 金属计算，主要采用Morse 势，但是由于通过实验拟合的对势容易导致柯西关系，与实验 不符，因此在后来的模拟中有人提出采用EAM 等多体势模型，或者采用第一性原理计算结 果通过一定的物理方法来拟合二体势函数。但是相对于二体势模型，多体势往往缺乏明确的 表达式，参量很多，模拟收敛速度很慢，给应用带来很大的困难，因此在一般应用中，通过 第一性原理计算结果拟合势函数的Lennard-Jones，Morse 等势模型的应用仍然非常广泛。 6.4.1

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

41、多的赘述了。 从图 6-9 可以清楚的知道原子之间的相互作用可以分为四类：对势、对泛函势、 组合势、 组合泛函势。下面就对相关的势函数的基本的发展做一个简单的总结。 表 6-1 原子中相互作用势函数 势函数作用因素典型代表 对势仅仅决定于原子的坐标 L-J 势、 BMH 势、 GK 势 对泛函势多体相互作用 EAM 势 有效介质理论 组合势共价键的相互作用SW 势 组合泛函势 Abell-Tersoff 势 CheliKowsky-Phillips 势 6.4.2分子力场 分子力场根据量子力学的波恩奥本海默近似，一个分子的能量可以近似看作构成分子 的各个原子的空间坐标的函数，简单地讲就是分子的

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

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

44、分子动力学方法 Computational Materials Science ：From Basic Principle to Practical Design Methodology - 139 - 图 6-9 分子力场的发展和分类 传统力场 (第一代力场 ) a)AMBER 力场：由 Kollman 课题组开发的力场，是目前使用比较广泛的一种力场，适合 处理生物大分子； b)CHARMM力场：由Karplus 课题组开发，对小分子体系到溶剂化的大分子体系都有很 好的拟合； c)CVFF 力场： CVFF 力场是一个可以用于无机体系计算的力场； d)MMX 力场： MMX力场包括MM2 和

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

46、以通过自主设定一系列分子作为训练集来生成合用的力场参数； a)ESFF 力场 MSI 公司开发的力场，可以进行有机、无机分子的计算； b)UFF 力场可以计算周期表上所有元素的参数； c)Dreiding 力场适用于有机小分子、大分子、主族元素的计算。 3.分子力场的计算 对于分子之间的相互作用通常需要将计算体系的所有粒子纳入考虑范畴之内，那么分子 之间的相互作用可以作为长程作用来处理。在长程相互作用中又分为周期性长程相互作用力 和非周期性相互作用力。周期性的长程作用力通常可以利用无网格算法，包含Ewald 求和 算法等； 而对于非周期性的相互作用通常使用快速多极子方法计算。下面针对这两个最为

47、常 见的算法做简单的叙述。 Ewald求和算法 Ewald 求和算法是Ewald 于 1921 年提出，最初是以计算离子晶体的能量。此方法选定 力场 CHARMM力场 CVFF力场 MMX 力场 AMBER力场 COMPASS 力场 MMF94 力场 CFF力场 UFF力场 Dreiding力 场 ESFF 力场 第一代力场第二代力场第三代力场 第 6 章分子动力学方法 Computational Materials Science ：From Basic Principle to Practical Design Methodology - 140 - 一个计算系统的盒子，盒子中的粒子除了和

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

49、能可以简单的表述为 1 ( ) N j i j ij q r rr (6.29) 因此这个体系的静电作用和前文所讲述的Lennard-Jones 作用项具有相同的形式。快速多极 子法将需要计算的体系分为多个体积相同的格胞，由格胞内部的所有的原子计算每一个格胞 的多极 (电荷，偶极，四极) 。而格胞和格胞之间的相互作用是利用中央多极展开的方法进 行计算。 因此中央多极展开的方法的条件是格胞之间的距离必须是大于一定距离的，因为大 于一个格胞以上的以多极展开的方式计算相互作用的力，而格胞内部的原子对是以成对原子 的作用力来计算的。 【练习与思考】 6-5.查找文献，推导Lennard-Jones 势函数的公式。 6-6.推导 EAM 势函数的基本公式，并且实现简单的算法，编写基本设计程序。 6-7.以对势中的Lennard-Jones 势函数为例，写出对势中势函数的基本算法，并且进行基 本的程序设计。 6.5数值方法与编程技巧 6.5.1边界条件 随着材料的低维化和纳米化，材料的表面和边界对于其性质的影响愈来愈大。在处理物 质表面、界面， 以及块体物质时，对其周期边界条件问题的考虑就变得非常必要了。对于周 期性边界条件可以分为一维、二维及三维的情况。介绍薄膜(可使用二维周期边