第6章-分子动力学方法
第6章分子动力学方法
经典分子动力学方法无疑是材料，尤其是大分子体系和大体系模拟有效的方法之一。分子动力学可以用于NPT，NVE，NVT等不同系综的计算，是一种基于牛顿力学确定论的热力学计算方法。与蒙特卡罗法相比在宏观性质计算上具有更高的准确度和有效性，可以广泛应用于物理，化学，生物，材料，医学等各个领域。本章在介绍分子动力学的基本概念的基础上，简单介绍了分子动力学的基本思想，势函数分类和基本方程。然后介绍了分子动力学的常用系综和典型的NPT，NVE，NVT系综基本方程。结合材料建模中的基本简化方法和技巧，阐述了边界条件和时间积分的数值处理技巧。最后，利用统计力学的基本概念给出分子动力学的计算信息的解析方式。并且结合Materials Explore软件计算分析了CNT的几何结构稳定性。
6.1引言
分子动力学方法(Molecular Dynamics, MD)方法是一种按该体系部的禀动力学规律来计算并确定位形的变化的确定性模拟方法。首先需要在给定的外界条件下建立一组粒子的运动方程，然后通过直接对系统中的一个个粒子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计力学方法得到多体系统的静态和动态特性，从而获得系统的宏观性质。可以看出，分子动力学方法中不存在任何随机因素，这个也是分子动力学方法和后文要提到的蒙特卡洛方法的区别之一。在分子动力学方法的处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学定律(或者是拉格朗日方程)。每个分子运动的禀动力学是用理论力学上的哈密顿量或者拉格朗日函数来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现玻尔兹曼的统计力学途径。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是分子动力学方法的计算机程序相对蒙特卡罗较复杂，其计算成本较高。
分子动力学方法发展历史改革经历了近60年，分子动力学方法是20世纪50年代后期由Alder B J 和Wainwright T E创造发展的。Alder和Wainwright在1957年利用分子动力学模拟，发现了“刚性球组成的集合系统会发生由其液相到结晶相的相转变”，后来人们称这种相变为Alder相变。其结果表明，不具有引力的粒子系统也具有凝聚态。到20世纪70年代，产生了刚性体系的动力学方法，成功地被应用于水和氮等分子性溶液体系的处理；1972年，Less A W和 Edwards S F 等发展了该方法并扩展到了存在速度梯度，即处于非平衡状态的系统。之后，此方法被 Ginan M J等推广到了具有温度梯度的非平衡系统，从而构造并形成了所谓的非平衡分子动力学方法体系。进人20世纪80年代之后，出现了在分子部对一部分自由度施加约束条件的新的分子动力学方法，从而使分子动力学方法可适用于类似蛋白质等生物大分子的解析与设计。分子动力学方法真正作为材料科学领域的一个重要研究方法，开始于由Andersen, Parrinello和Rahman等创立恒压分子动力学方法和Nǒse 等完成恒温分子动力学方法的建立及在应用方面的成功。后来，针对势函数模型化比较困难的半导体和金属等，1985 年人们又提出了将基于密度泛函理论的电子论和分子动力学方法有机统一起来的所谓Car-Parrinello 方法，亦即第一性原理分子动力学方法。这样，分子动力学的方法进一步得到发展和完善，它不仅可以处理半导体和金属的间题，同时还可应用
于处理有机物和化学反应。关于Car-Parrinello 分子动力学方法将在第7章重点学习讨论。本章重点讨论经典分子动力学方法的基本原理和计算方法。
6.2 分子动力学的计算框架
6.2.1 基本思想
分子动力学方法是沿用牛顿运动方程或拉格朗日方程、哈密顿方程，通过考察粒子的运动来研究多粒子系统的物理性质。在处理孤立粒子或原子团簇(不考虑与其他粒子的相互作用)时，单纯地对牛顿运动方程进行积分求解就行了；而在处理凝聚系统时，可以认为将所考虑的N 个粒子放入具有一定体积的容器中，在这样组成的封闭系中，建立直角坐标(x, y, z)，每个粒子的位置就由三个坐标分量决定，通过求解3N 个联立方程组就可以得到保守系的总能量。分子动力学的基本思想在于通过设定原子之间的相互作用(势函数)和相关的系综(亦即作用对象和条件)，确定其基本的模拟畴。在给定的牛顿方程、拉格朗日方程或者是哈密顿方程进行时间的迭代，在达到指定的力学收敛条件之后得到一个最终的位置坐标，然后通过相关的位形信息和运动的速度和加速度等信息通过热力统计学的方法，给出模拟体系的统计信息，主要包括复杂的光学，电学和一些其他的物理信息。其基本思想，总结如图6-1所示。
一次信息二次信息动力学方程初始化条件
图 6-1 分子动力学的一般思想
6.2.2 计算流程
总体来说，分子动力学的基本计算有如下四个步骤：
1. 模型的设定，也就是势函数的选取
势函数的研究和物理系统中对物质的描述研究息息相关。分子动力学的模拟最早使用是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是Lennard-Jones 势函数，还有嵌入原子EAM 势函数，不同的物质状态描述用不同的势函数。模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。关于势函数选取将在6.4节进一步讨论。
2. 给定初始条件
运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如Verlet 算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。
一般意义上讲系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条件，因为当模拟时间足够长时，系统就会忘掉初始条件(对于无记忆的体系而言)。当然，合理的初始条件可以加快系统趋于平衡的时间和步程，获得好的精度。
常用的初始条件可以选择为，令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。
图6-2 分子动力学计算流程
3.趋于平衡计算
在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，增加或者从系统中移出能量，直到持续给出确定的能量值，称此时系统已经达到平衡，达到平衡的时间称为弛豫时间。
分子动力学中，积分步程的大小选择十分重要，这决定了模拟所需要的时间。为了减小误差，步长要小，但过小系统模拟的弛豫时间就长。因此根据经验选择适当的步长。如，对一个具有几百个氩气分子的体系，采用Lennard-Jones势函数，发现取h为0.01量级，可以得到很好的结果。这里选择的h是没有量纲的，实际上这样选择的h对应的时间在10-14s 的量级。如果模拟1000步，系统达到平衡，弛豫时间只有10-11s。
4.宏观物理量的计算
实际计算宏观的物理量往往是在模拟的最后阶段进行的，它是沿相空间轨迹求平均来计算得到的，根据各态历经假说，时间平均代替系综平均。关于各态历经假说将在第8章的蒙特卡洛方法中深入讨论。
6.2.3经典分子动力学中近似处理
实际的计算体系由于其外在条件的复杂性，往往需要对计算的系统进行一定的近似处理，在经典分子动力学计算中，引进以下近似处理：
?经典粒子相互作用，不考虑电子相互作用量子效应。
?力的作用形式，由参数可调的相互作用势函数决定，并经过实验测定来验证。
?模拟体系与实际体系相差较大，一般需要采用周期边界来扩展计算体系。
时间平均是在有限时间完成。
【练习与思考】
6-1.查找文献，根据上述的分子动力学的处理流程图，编写实现分子动力学的简单程序，可参考Daan F和Berend S编著的《分子模拟》一书。
6.3分子动力学的系综
系综(ensemble)是指在一定的宏观条件下(约束条件)，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统集合，或称为统计系综。系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念；系综是统计理论的一种表述方式，系综理论使统计物理成为普遍的微观统计理论；系综并不是实际的物体，构成系综的系统才是实际物体。约束条件是由一组外加宏观参量来表示。在平衡统计力学畴下，可以用来处理稳定系综。
6.3.1常用系综分类
1.正则系综(canonical ensemble)
全称应为“宏观正则系综”，简写为NVT，即表示具有确定的粒子数(N)、体积(V)、温度(T)。正则系综是分子动力学方法模拟处理的典型代表。假定N个粒子处在体积为V的盒子，将其埋入温度恒为T的热浴中。此时，总能量(E)和系统压强(P)可能在某一平均值附近起伏变化。平衡体系代表封闭系统，与大热源热接触平衡的恒温系统。正则系综的特征函数是亥姆霍兹自由能(,,)
F N V T。
2.微正则系综(micro-canonical ensemble)
简写为NVE，即表示具有确定的粒子数(N)、体积(V)、总能量(E)。微正则系综广泛被应用在分子动力学模拟中。假定N个粒子处在体积为V的盒子，并固定总能量(E)。此时，系综的温度(T)和系统压强(P)可能在某一平均值附近起伏变化。平衡体系为孤立系统，与外界即无能量交换，也无粒子交换。微正则系综的特征函数是熵(,,)
S N V E。
3.等温等压(constant-pressure, constant-temperature)
简写为NPT，即表示具有确定的粒子数(N)、压强(P)、温度(T)。其总能量(E)和系统体积(V)可能存在起伏。体系是可移动系统壁情况下的恒温热浴。特征函数是吉布斯自由能G N P T。
(,,)
4.等压等焓(constant-pressure, constant-enthalpy)
简写为NPH ，即表示具有确定的粒子数(N)、压强(P)、焓(H)。由于H E PV =+，故在该系综下进行模拟时要保持压力与焓值为固定，其调节技术的实现也有一定的难度，这种系综在实际的分子动力学模拟中已经很少遇到了。
5. 巨正则系综(grand canonical ensemble)
简写为μVT，即表示具有确定的粒体积(V)、温度(T)和化学势(μ)。巨正则系综通常是蒙特卡罗模拟的对象和手段。此时、系统能量(E)、压强(P)和粒子数(N)会在某一平均值附近有一个起伏。体系是一个开放系统，与大热源大粒子源热接触平衡而具有恒定的温度(T)。特征函数是马休(Massieu)函数J(μ，V ，T)。
对于不同的系综，由于其选择的运动方程的求解方法亦有不同。针对不同的计算体系具体的细节比较复杂，下文将以NEV 和NTP 这两类简单的系综运用两类不同的条件方法下求解的不同方程做简要的介绍。
6.3.2 NEV 系综基本方程
NEV 系综是一类较为常见的系综。对于该系综，经典分子动力学方法一般采用差分近似法求解牛顿运动方程，并追踪系统的时间变化。考虑由N 个粒子构成的系统，假设第i 个粒子在时刻()i x t 位置，速度为()i v t 受到的作用力为()i F t ，为求出该粒子(i)在时刻t t +?的位置，经常使用的方法是Verlet 算法(对于该算法的具体细节在下文将有详细的推导)。
根据系统的动能及统计物理学中的能量均分定理，则系统的平衡温度可表达为 2212i i B f i p T k N m ??=? ? ????∑ (6.1)
根据此方法，在晶体中粒子的配置与排布、晶形变化等结构相变也可以用计算机模拟进行研究。
根据物理化学的基本气体的基本方程可以知道，体系的NEV 系综下的计算体系主要变化的量有两个压力和温度。下文就从恒温和恒压两个条件下对NEV 系综的基本方法做简单的推导。
1. 恒温方法
N óse S 和Hoover W G 引入与热源相关的参数?，尝试构造恒温状态，亦即具有确定N ，V ，T 值的系统，对此可设想为与大热源接触而达到平衡的系统。由于热源很大，交换能量不会改变热源的温度。在两者建立平衡后，系统将与热源具有相同的温度，系统与热源合起来构成一个复合系统，如图6-3所示。
这样，系统中微观粒子的动力学方程：
i i i dq p dt m = (6.2)
i i i dp p dt dq ????Φ=-- ??? (6.3) 232[]22i B ex i i p d N k T dt m Q ???=-??? ???
∑ (6.4) 式中Φ为系统的势函数，Q 为有效质量，?表示热力学摩尔系数。式(6.2)至式(6.4)
同往常的分子动力学方法的区别体现在式(6.3)中，即增加了与热源的相互作用相关的并与力的量纲相同的一项(i p ?)。与热源相关的变化参数?的运动方程表明，当系统的总能量大于32
B Nk T 时，?是增加的，从而显示出使粒子速度减小那样的作用；反之则显示出使粒子
速度增大；Q 是表示与温度控制有关的常数。
图6-3 采取扩展系方法恒温分子动力学方法原理示意图
该方法的特征是系统处于微观状态的分布服从正则分布，若在推导式(6.2)、式(6.3)、式(6.4)中，引人由公式
ln d s dt ?=定义的变量，则系统的哈密顿量(即能量)可以表示为 203()ln 22B Q H H Nk T s ?=+?+ (6.5) 式中第二项和第三项对应于与热源相联系的部分。这说明，由于系统与热源存在热接触，二者可以交换能量，因此粒子的能量发生变化(H 0→H)，同时也说明系统可能的微观状态可具有不同的能量值。
2. 恒压方法
为了进行压力的控制，考虑如图6-4所示的方法，即利用活塞调控系统的体积变化。对于质量一定的理想气体系统，体积变大，其压力就变小，反之亦然。利用这样的原理，可以通过改变体积使压力达到某一个预定值。就实际的模拟而言，可控制体积大小均匀(三维)地变化，将粒子的坐标、动量表示成如下形式
1'31'3i i i i q V q p V p -?=????=?? (6.6)
从而将粒子的坐标、动量与对分子动力学单元(边长1/3L V =)归一化的参变量i q 和i p 联系起来。应用归一化的参变量并根据正则方程可得到粒子运动方程。
图6-4 恒压分子动力学方法原理示意图(实际上，活塞为三维运动)
粒子系和活塞组成复合系统，其哈密顿量可由下式给出
202v ex p H H P V W =++ (6.7) 212330()2q i i i p H V
V m -??=+Φ ???
∑ (6.8) 式中，ex P 是外压；v p 是体积V 对应的共轭动量；W 是决定于体积变化速度的因子。若将由上述哈密顿量导出的运动方程用粒子的实际坐标直接写出来，则有
''3i i i i dq p q dV dt m dt V =+? (6.9) '''()3i i i dp p dV dt dt V q ?Φ=--?? (6.10)
'2''22[()]()3i i i i i i ex p q m q d V W P V dt ?Φ-?∑??=-∑ (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. 压力恒定体系
为简单起见，以由圆柱容器充满气体和调节活塞构成的系统为例，如图6-5所示。设活塞的重量为W ，则在平衡状态下，由活塞的重量产生的外部压力(ex P )与由气体产生的压平
衡。如果外压ex P 变大，为了保持平衡，则气体部分的体积将收缩，从而压变大。
图6-5 活塞与气体容器构成的示意图
因此，把体积作为系统的动力学参数建立其运动方程，就好像维持压力恒定。其方程可写为
22()ex d V W P V dt ?=Λ-? (6.12)
式中，V 为体积，Λ是由维里(Virial)定理求出的压，W 是对应于V 的假想质量。图6-5为活塞与气体容器构成的示意图。将上述想法推广到一般情况下的晶格系统，如图6-6所示，对一般的晶体，其基本单元可由三个平移矢量 (,,)x y z a a a a =,(,,)x y z b b b b =, (,,)x y z c c c c =决定其形状和大小。将,
,a b c 的三个分量写成矩阵h ，该矩阵规定了晶体基本单元的形状和大小，这样h 也成为原子的实际坐标i r 和晶格矢量i s 之间的变换矩阵，亦即
x y z x y z x y z a a a h b b b c c c ????=?????? i i r hs =
c
a
图6-6 基于单元形状示意图
在式(6.15)中，体积V 作为动力学参数，而在一般晶体系统中，h 应该作为动力学参数。若给出此动力学系统的哈密顿量，则有
1121()1(,,,)()222t T N i i N ex i i G tr H r r r P tr G m W ππ-=∏∏=+Φ++Ω+Γ∑ (6.13) 式中，i π是与格矢i s 相应的原子共扼动量，T G h h =为数值行列式，Φ中是多粒子体系
的势函数，∏是h 中与体积V 相应的共扼动量，Ω是基本单元的体积，且由()()
h a b c Ω=?=??给出，W 为相应的质量，ex P 为外部压强，Γ代表各向异性应力量。
2. 温度恒定体系
接下来讨论满足温度恒定的基本方法。恒温体系比定压体系要抽象，简单地给出视觉化图像是困难的。势能在原子坐标和动量之外引入了附加自由度f ，并由下式用f 把现实系统与具有能量确定的(微正则系统)假想系统联系起来。 dt dt f '= (6.14) 式中，t 为现实系统的时间，'t 为假想系统的时间。该假想(力学)系统的哈密顿量'H 可以写成 221221()(,,,)(ln )22N f i N B ex i i P p H r r r g k T f M m f ='''''=+Φ++???∑ (6.15)
式中，i P 是在假想系统中原子的动量，f P 是f 对应的共扼动量，g 是由公式1g N =+给出的假想系统的自由度数目，ex T 二是热浴的温度，B k 为玻尔兹曼常数。显然，式(6.18)中
右边第一项代表粒子的动能，第二项为粒子系统的势能，第三项是热浴的假想动能，第四项是热浴的势能；因为在假想系统中，力学系综是微正则分布的，所以其配分函数NEW 'Z 变成具有总能量 E 的状态总数
12121212()N N f NEW
N N f H E dr dr dr dp dp dp dfdP dr dr dr dp dp dp dfdP δ''''''-'Z =''''''?????? (6.16) 12121212()N N f NEW N N f H E dr dr dr dp dp dp dfdP dr dr dr dp dp dp dfdP δ''''''-'Z =''''''?????? (6.17) 若将此配分函数积分变换到现实系统中，则可得到在恒定温度T 下的正则分布(也称为玻尔兹曼分布)的配分函数 NTV Z
121212
12exp()N N B NTV N N H dr dr dr dp dp dp k T dr dr dr dp dp dp -Z =??????? (6.18) 式中，H 是由下式给出的现实系统的哈密顿量 2121(,,,)2N i N i i p H r r r m ==+Φ∑ (6.19)
'H 是假想系统的守恒量，而H 对真实系统来讲，由于存在热交换将不是守恒量。为了判断计算结果是否合理，只要简单地检验'H 是否保持恒定(守恒)就行了。
3. NTP 系综
到此，已经导出了压力恒定系综和温度恒定系综的哈密顿量，将两个系综的哈密顿量组合起来就可以得到NTP 系综的哈密顿量，其模型如图6-7所示。
NTP 系综的哈密顿量为
1122212()(,,
,)221()(ln )22t T N i i N i i f ex B ex G tr H r r r m f W f P P tr G gk T f M ππ-=∏∏=+Φ+?+?Ω+Γ++?∑ (6.20)
基本单元
图6-7 NTP 系统示意图
若根据此哈密顿量导出正则方程式，并用速度置换动量，则可得到各变数的运动方程
2121()N i i ij ij i i j d s f m s m G G s f dt χ-==-+∑ (6.21)
22()ex d h Wfh W P h f dt σσ=Λ--Γ- (6.22) 2221
[()()()]()N T T i i i B ex i d f f M m hs hs Wtr h h gk T f M f dt ==+-?+∑ (6.23) 式中，变量上方的圆点表示对时间的微分，而,,x σΛ分别由下列各式给出
1ij ij ij ij d r dr χψ=- (6.24)
111[()()()]N N N i i i ij ij ij i i j i m hs hs x r r ==>Λ=?+?Ω∑∑∑ (6.25)
1h σ-=Ω (6.26) 4. NPT 系综分子动力学求解
若在三维周期边界条件下求解上述方程，则可模拟固体的结构相变、溶解现象和晶化过程。各变量的运动方程式(6.21)，式(6.22)和式(6.23)是非线性常微分方程，不能用贝鲁勒(Berrele)法进行求解。因此，对这些方程通常用牛顿迭代方法求解。而问题的关键在于确定式(6.22) 和式(6.23)中的假想质量。所考察基本单元的粒子体系存在着由粒子作用势以及粒子数与粒子质量共同决定的固有压力和温度等的起伏变化。若想再现这些波动周期那样给出M 和W ，则由体积变化或温度变化引起的粒子系统对平衡态的偏离将迅速回复而变成稳定状态。因此，若设此周期分别为P t 和T t ，则M 和W 可以表示为
26()2T B t M Nk T π= (6.27) 23()()2P t L W k π
= (6.28) 式中，T 为温度，B k 为玻尔兹曼常数，N 为原子数，L 为基本单元的线尺度，k 为压
缩率，P t 和T t 分别根据对NEV 系综的计算结果而定。
【练习与思考】
6-2.推演NEV系综的牛顿方程，哈密顿方程，说明各个物理量的意义，在那些实际的计算体系中得到应用。
6-3.根据NPT系综的推导方法，对NPH系综的哈密顿方程做简单的推导。
6-4.详细的分析不同的系综的使用畴，给出各个系综的使用条件，列表进行简单的比较。
6.4原子势函数和分子力场构造
分子动力学中存在着两个基本的概念，一是在前面提到的系综的问题，它所确定的是模拟系统的大小，规格和相关的模拟条件等信息；二是本节需要讲述的势函数。简单的讲，势函数就是用来描述原子和原子之间相互作用的。随着势函数的不断发展，现有的势函数主要分成了两个大的部分。即是经典的势函数和电子理论的势函数。如图6-8所示，经典分子动力学的基本的势函数的一个简单分类。
图6-8经典的分子动力学相互作用势分类
作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、Morse势等二体势模型，对于金属计算，主要采用Morse势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是相对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大的困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的Lennard-Jones，Morse等势模型的应用仍然非常广泛。6.4.1经典理论的原子势函数
在分子动力学的框架之下，早期的势函数是主要以描述原子和原子之间的相互作用而产生和发展起来的。早在1903年的时候，Mie G就研究了两个粒子之间的相互作用势，他给出的势函数是由两项构成的，一项代表原子之间的相互吸引，另一项则代表着原子之间的相
互排斥。这就是非常著名的Lennard-Jones势函数的原始思考。对于其他有关势函数的发展在这里就不做过多的赘述了。
从图6-9可以清楚的知道原子之间的相互作用可以分为四类：对势、对泛函势、组合势、组合泛函势。下面就对相关的势函数的基本的发展做一个简单的总结。
表6-1 原子中相互作用势函数
6.4.2分子力场
分子力场根据量子力学的波恩—奥本海默近似，一个分子的能量可以近似看作构成分子的各个原子的空间坐标的函数，简单地讲就是分子的能量随分子构型的变化而变化，而描述这种分子能量和分子结构之间关系的就是分子力场函数。分子力场函数为来自实验结果的经验公式，可以讲对分子能量的模拟比较粗糙，但是相比于精确的量子力学从头计算方法，分子力场方法的计算量要小数十倍，而且在适当的围，分子力场方法的计算精度与量子化学计算相差无几，因此对大分子复杂体系而言，分子力场方法是一套行之有效的方法。以分子力场为基础的分子力学计算方法在分子动力学、蒙特卡罗方法、分子对接等分子模拟方法中有着广泛的应用。
1.力场主要讨论的对象及其构成：
以下参数构成一套力场函数体系需要有一套联系分子能量和构型的函数，还需要给出各种不同原子在不同成键状况下的物理参数，比如正常的键长、键角、二面角等，这些力场参数多来自实验或者量子化学计算。
?键伸缩能：构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化；
?键角弯曲能：键角变化引起的分子能量变化；
?二面角扭曲能：单键旋转引起分子骨架扭曲所产生的能量变化；
?非键相互作用：包括德华力、静电相互作用等与能量有关的非键相互作用；
?交叉能量项：上述作用之间耦合引起的能量变化。
2.常用力场函数和分类
不同的分子力场会选取不同的函数形式来描述上述能量与体系构型之间的关系。到目前，不同的研究小组设计了很多适用于不同体系的力场函数，根据他们选择的函数和力场参数，可以分为以下几类：
图6-9 分子力场的发展和分类
? 传统力场(
第一代力场)
a) AMBER 力场：由Kollman 课题组开发的力场，是目前使用比较广泛的一种力场，适合处
理生物大分子；
b) CHARMM 力场：由Karplus 课题组开发，对小分子体系到溶剂化的大分子体系都有很好
的拟合；
c) CVFF 力场：CVFF 力场是一个可以用于无机体系计算的力场；
d) MMX 力场：MMX 力场包括MM2和MM3，是目前应用最为广泛的一种力场，主要针对有机
小分子。
? 第二代力场
第二代的势能函数形式比传统力场要更加复杂，涉及的力场参数更多，计算量也更大，当然也相应地更加准确；
a) CFF 力场CFF 力场是一个力场家族，包括了CFF91、PCFF 、CFF95等很多力场，可以进
行从有机小分子、生物大分子到分子筛等诸多体系的计算；
b) COMPASS 力场由MSI 公司开发的力场，擅长进行高分子体系的计算；
c) MMF94力场Hagler 开发的力场，是目前最准确的力场之一。
? 通用力场(第一代力场)
通用力场也叫基于规则的力场，它所应用的力场参数是基于原子性质计算所得，用户可以通过自主设定一系列分子作为训练集来生成合用的力场参数；
a) ESFF 力场MSI 公司开发的力场，可以进行有机、无机分子的计算；
b) UFF 力场可以计算周期表上所有元素的参数；
c) Dreiding 力场适用于有机小分子、大分子、主族元素的计算。
3. 分子力场的计算
对于分子之间的相互作用通常需要将计算体系的所有粒子纳入考虑畴之，那么分子之间的相互作用可以作为长程作用来处理。在长程相互作用中又分为周期性长程相互作用力和非周期性相互作用力。周期性的长程作用力通常可以利用无网格算法，包含Ewald 求和算法等；而对于非周期性的相互作用通常使用快速多极子方法计算。下面针对这两个最为常见的算法做简单的叙述。
? Ewald 求和算法
Ewald 求和算法是Ewald 于1921年提出，最初是以计算离子晶体的能量。此方法选定第一代力场 第二代力场 第三代力场
一个计算系统的盒子，盒子中的粒子除了和盒子部的粒子相互作用之外，还会与其镜像系统中的粒子发生作用，并且镜像系统和盒子部的系统完全相同。图显示的是计算的相互作用的系统，中央的黑色的区域表示的是计算的区域，周围的镜像随着系统的距离逐渐变淡，那么这样的一个极限作用区域就像一个球，这个球称为Ewald 球。
假设作用的系统的边长是L ，含有n 个相互作用的粒子，其镜像的作用位置可以表示为(,,iL jL kL ±±±)；其中,,i j k 都是正整数。
? 快速多极子法(The Fast Multipole Method)
在一个开放的系统中，不是对所谓的晶格进行求和，而是对要求计算的整个体系的粒子对进行求和。那么这样的一个系统不同位置粒子的静电能可以简单的表述为 1()N j i j i j q r r r φ==-∑ (6.29)
因此这个体系的静电作用和前文所讲述的Lennard-Jones 作用项具有相同的形式。快速多极子法将需要计算的体系分为多个体积相同的格胞，由格胞部的所有的原子计算每一个格胞的多极(电荷，偶极，四极…)。而格胞和格胞之间的相互作用是利用中央多极展开的方法进行计算。因此中央多极展开的方法的条件是格胞之间的距离必须是大于一定距离的，因为大于一个格胞以上的以多极展开的方式计算相互作用的力，而格胞部的原子对是以成对原子的作用力来计算的。
【练习与思考】
6-5. 查找文献，推导Lennard-Jones 势函数的公式。
6-6. 推导EAM 势函数的基本公式，并且实现简单的算法，编写基本设计程序。
6-7. 以对势中的Lennard-Jones 势函数为例，写出对势中势函数的基本算法，并且进行基
本的程序设计。
6.5 数值方法与编程技巧
6.5.1 边界条件
随着材料的低维化和纳米化，材料的表面和边界对于其性质的影响愈来愈大。在处理物质表面、界面，以及块体物质时，对其周期边界条件问题的考虑就变得非常必要了。对于周期性边界条件可以分为一维、二维及三维的情况。介绍薄膜(可使用二维周期边界条件)和扩大到半无限的表面情况，然后讨论液体与晶体、晶体与晶体之间形成的界面等，最后学习块状物质的三维周期边界条件。
1. 薄膜的边界条件(二维周期边界条件)
对于薄膜情况，如图6-10所示赋予其周期性边界条件。可以认为薄膜在x-y 平面无限扩展，存在周期边界条件，而在 z 方向受到限制(不赋予周期边界条件)。事实上，二维周期性边界条件也适用于在z 方向不具有真正意义上的自由度的二维物质。但就物质科学的研究畴而言，多数情况下是指在某种衬底上生成的原子或分子的薄膜。
此处讨论的边界条件，除针对薄膜以外，也可近似地用于下面将要讲到的在 z 的负方向半无限地扩展的表面问题，例如表面重构等。但是，当在z 方向存在异质原子层时，原来块体物质在异质原子层附近将有近于表面的性质，这种由块体向表面的过渡特征取决于原子
层的厚度，一般要求原子层厚度在力的中断距离以上(作用力程以上)。在处理实际问题时，只要尝试进行几次改变原子层厚度的模拟和计算，并检验计算结果是否变化就行了。
图6-10 薄膜周期性边界条件示意图
2. 表面边界条件(半无限的情况)
对于表面存在表面重构的情况，是与薄膜的情况不同的，其在 z 轴方向扩展至半无限。处理此类系统之边界条件时，只对x-y 平面和 z 的负方向赋予相应的周期性边界条件，其示意图如图6-11所示。但是，在将基本单元原样地复制中，必须考虑到这样的一种情况，那就是在距薄膜结构的基本单元最下层的相邻底层被配置了基本单元表面层的复制品，从而导致了在本来应具有块体物质性质的基本单元的最下层的结构中出现了表面层的强烈影响。所以将基本单元对z 轴反转(即镜像映射)进行复制后使用。然而，即使在这种情况下，考虑到基本单元在z 方向的厚度对计算结果有影响，所以应当进行与前面所述同样的检验。
基本单元复制单元
复制单元
z 表层
图6-11 采用映射方法半无限表面的周期性边界条件
基本单元复制单元
复制单元
图6-12 固定层法半无限表面的周期性边界条件
在利用分子动力学方法对分子束外延的晶体生长模拟时，要充分认识到分子束外延 (MBE)是一个原子数不断增加的非平衡系。所谓MBE 技术就由束源产生的原子蒸发到衬底上，从而形成晶体。因为由束源出来的原子具有动量，所以碰撞衬底，衬底的重心就变得下移，为了克服这种现象的影响，对分子束外延即在z 轴方向不赋予周期性边界条件，而采用通过调整配置使原子位置坐标不变的所谓固定原子层方法，然后使用边界条件，如图6-12所示。
固定(原子)层还可理解为具有无限大质量的原子层。来自固定层中原子作用力的贡献可作为外力处理。为了基本单元最下层具有块体物质的性质，固定层中的原子层数目必须足够多。采用镜像映射的方法可自然地处理体系伴随温度变化而产生的膨胀和收缩。但就使用固定层的方法而言，基本单元发生变化的同时，必须将固定层的晶格常数作为外部参量变化来考虑。
【练习与思考】
6-8. 通过对比映射方法和固定层法的不同，理解在分子动力学中的二维边界条件的设置方
法，思考在线形的一维约束条件的边界的处理方法。
6-9. 根据界面的边界条件的设置方法，分析讨论在分子动力学中如何构造线缺陷模型，做
要的图示和语言表述，说明界面的条件的选取方法。
6-10. 根据上述的边界的处理方法，试思考块体材料的边界条件的实现方法，并且给出简单
的实例。
6.5.2 时间积分处理
计算机模拟方法的基点是利用现代计算机高速和精确的优点，对几百个以至上千分子的运动方程进行数值积分。有许多不同的积分方法，它们的效率和方便程度各异。基本问题就是用有限差分法来对二阶常微分方程进行积分。而对于分子动力学而言其基本的工作也就是对运动方程进行时间上的积分。这样需要将运动方程离散化为有限的差分方程，即是将时间离散化为有限的格点。主要的方法是利用泰勒展开采取不同的截断半径。而相邻的格点之间的距离就是所谓的时间步长，这样在知道了原始的位置坐标之后，在知道某些物理量对于时间t 的微分量之后就可以很容易确定下一个时间体系的坐标。这里很容易出现所谓的两类误差，即是截断误差和舍入误差。为此我们需要在实际的模拟中，应该根据不同的问题来确定所需要的精度，为此发展了几种时间积分的方法。
1. Verlet 方法
Verlet 法经常用于NEV 系综的分子动力学微分方程的数值积分求解。若将在t t +?和t t -?的位置坐标分别用时刻t 的位置坐标进行泰勒展开，则有 2()()()()()2
a t r t t r t t V t t +?=+??+??+ (6.30) 2()()()()()2a t r t t r t t V t t -?=-??+??- (6.31)
为简单起见，以后均省略有关原子标识符号I ，V(t)为速度，为了消去速度项，将式(6.30)
和式(6.31)相加，并将r(t t -?)移至右边，同时利用()()F t a t m =
，则有 2()()2()()()F t r t t r t r t t t m +?≈--?+?? (6.32) 上式(6.32)是差分公式，其中(t ?)的奇次项全被抵消，同时忽略4()t ?及其他高次项。差分精度(误差)为0(4()t ?)。若由式(6-36)减去式(6-37)则可得到速度表达式 ()()()2r t t r t t V t t +?--?≈? (6.33) 因为是消去了关于t ?的偶数方项，所以速度的误差为2()t ?的量级。上述式(6.32)右边第一、二项均是0()t ?的量级项，而第三项是2()t ?量级项。这样若将大的绝对值项和小的绝对项混合，则随着重复计算次数的增多，其误差发生积累。为防止误差积累，Hockney 提出了所谓蛙跃(Leap-frog)法。
2. 蛙跃法
首先，利用速度将二阶常微分方程转换为二个一阶常微分方程。亦即 dV m F dt dr V dt
?=????=?? (6.34) 速度的微分由时刻，/2t t +?和/2t t -?的差分给出 ()()()22t t V t V t F t t m ??+--=? (6.35)
根据此式，在时刻/2t t +?的速度给出如下 ()()()22t t t V t V t F t m ???+=-+? (6.36) 另一方面，利用t t +?时刻和t 时刻的差分可给出位置矢量的微分 ()()()2
V t t r t t V t t +?-?=+? (6.37) 根据上式可给出时刻t 的位置矢量表达式 ()()()2t V t t r t t V t ?+?=+??+ (6.38)
上式(6.38)所表示的意义是速度和位置矢量的变化依赖步长/2t ?。利用时刻的位置坐标可计算出力，利用此力的结果求出时刻(/2t t +?)的速度，进而，利用此速度求出在(/2t t +?)时刻的位置坐标，这样就完成了第一步的计算。正如以上所述，如果使用Verlet 法或蛙跃法，则可求解关于NEV 系综的线性常微分方程。但是，在此之外的NTV ，NHP ，NTP 等系综给出的是非线性常微分方程。所以用Verlet 法或蛙跃法不能进行求解。关于这些非线性常微分方程的求解要采用下面介绍的所谓Gear 方法。Hockey 提出的Leap-frog 算法是Verlet 算法的变化，这种方法涉及半时间间隔的速度。这种算法与Verlet 算法相比有两个优点：
(1)包括显速度项；(2)收敛速度快，计算量小。这种算法明显的缺陷是位置和速度不同步。
3. Gear 方法
因为原子轨迹是时间的连续函数，所以可将(/2t t +?)时刻的位置、速度、加速度等对
时刻t 泰勒展开 232()()()()()()()26()()()()() 2()()() ()() P P P P a t b t r t t r t tV t t t t b t r t t V t ta t a t t a t t b t b t t b t +?=+?+??+?+?+?=+?+++?=+??++?=+ ????????? (6.39)
若使用给出的新的位置矢量()p r t t +?，则可计算(t t +?)时刻的力F(t t +?)和加速度c a ／(t t +?)。根据加速度c a 可评价预测的加速度()p r t t +?的误差，亦即
()()()C P a t t a t t a t t ?+?=+?-+? (6.40) 若将此差值加于预测项上，则可得到修正项
0123()()()()()()()()()()()()C P C P C P C P r t t r t t C a t t V t t V t t C a t t a t t a t t C a t t b t t b t t C a t t ??+?=+?+?+???+?=+?+?+?????+?=+?+?+????+?=+?+?+???? (6.41)
1()()()2
11()()()22
111()[()()]222r t t r t t V t t v t t v t t t a t v t v t t v t t +?=+??+?+?=-?+??=+?+-? (6.42)
4. Velocity-Verlet 算法
这种算法可以同时给出位置、速度和加速度，并且不牺牲精度，给出了显速度项，计算量适中，目前应用比较广泛。 22()
()()()()2
1()()[()()]21()()()()22
1()()()()22a t r t t r t t V t t r t t v t t a t a t t a t r t t v t t a t t r t t v t t t +?=+??+??+?=+?++?+?=+??+?+?=+?+?? (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)目前能够讨论的问题已经可以涉及到比较广阔的畴，从基本的模拟体系的位形，到现在的各种材料的光谱学性能的计算。但是总体来说MD 的计算结果仅仅分了两个层次。第一个是计算模拟直接的输出，基本上包含了在平衡态的体系位形信息和体系运动的基本动力学信息；第二个是基于现代统计热力学的平均统计得到的一些宏观体系的综合信息，包含了光学，电学，磁学和光谱学等等。但是，在材料科学中，不仅要很好地知道原子水平上的三维结构，而且还必须计算出表征宏观物性(热力学性质、光学性质等)的物理量值。为了灵活应用分子动力学方法，下面对与这些物性相关的物理参数值的计算方法作介绍。
6.6.1 宏观参量的统计平均
由分子动力学基本方程出发进行计算的输出结果是关于原子的位置坐标和速度。物性参数值可根据位置坐标和速度通过统计处理求出。在统计物理中，可利用系综微观量的统计平均值来计算物性参数值，亦即： 1(,)(,)!ens A dpdrA p r p p r N <>=??
(6.44) 式中，ens 是系综(Ensemble)的英文缩写，(,)p p r 是分布函数(几率密度)。几率密度因所考察的系综不同而不同。例如，对于微正则系综(NEV 系)，(,)p p r 可表达为
(,)((,))NEW p p r H p r E δ=- (6.45) 而对于正则系综(NTV 系)，等温、等压系综(NPT 系)和巨正则系综，分别由下列各式给出：
(,)(,)exp[]NTV B H p r p p r k T -=
(,)(,)exp[]NPT B H p r PV p p r k T -+= (6.46) (,)(,)exp[]TV B H p r N p p r k T
μμ--= 系综之间的演变可由相应的公式联系起来。首先，配分函数Z 和热力学势lg B k T Z ψ=-，分别由下列式子给出
0exp()f F Z dF Ff Z ∞=-??
f F B k TFf ψ=ψ+ (6.47) 式中，在(f ,F )中，设定1B k T
β=，则可有以下组合：
(,)E β——NTV 和NEV 系综之间的转变
(,)P V β——NTP 和NTV 系综之间的转变
(,)N βμ-——μVT 和NTV 系综之间的转变
在上面最后一个演变例子中，关于N 的积分应理解为对N ＝0，1，…的求和。在NEV 系综情况下，F Z 表示所有状态数W ，热力学势F Ψ为(-TS)；对于NTP 和μVT 系综，其热力学势分别为μN 和(-PV)。此时，物理量在系综中的平均值的变换可由下式给出
exp()exp[()]f f F F A dF Ff A ββ<>=Φ-Φ+<>? (6.48)
在热力学极限情况下，上述各系综的平均是相等的
f F A A <>=<> (6.49)
到此，容易发现，有了上述这些分布函数，代入式(6.50)进行积分运算，理论上讲，可以求出各种各样的物理量。然而，在通常情况下，式(6.50)不能给出真正的解析解。在理想的刚体球势情况下，可以进行解析求解，而对于连续函数势，只能是数值求解。不难看出，这些解始终都是近似解。从历史的角度来说，正是为了超越上述积分方程带来的限制，人们引人了所谓蒙特卡罗方法和分子动力学方法。
在分子动力学方法中，人们使用了所谓的时间平均等于系综平均的各态历经假设，亦即
ens time A A <>=<> (6.49) 在热力学统计物理学中，各态历经假设没有严格地证明，它的正确性是由它的推论被实验结果所证实而被肯定。
式(6.59)中time A <>由下式给出
1
1lim (,)(,)time obs obs A A p r dt A p r t τ<>==∑? (6.50)
在气体动力论中，原子的速度分布为Maxwell 分布。在统计力学中，瞬间的热力学性质可以就其平均值做Gaussian 分布描述。所以，对任意一个状态函数X(t)，其瞬间的值是作无规地起伏变化，而此状态值X(t)则依据它的平均值定义了其分布函数。但是要使得性质分布函数X(t)的取样点为可幸赖的，则模拟的向空间长度要多长？在此，先来分辨三种形式的统计平均值：轨迹平均(trajectory average)、取样平均(sample average)、取样平均的期望值(the expected value of the sample average)。
轨迹平均：考虑一个独立系统的平衡相空间轨迹，并且想象在轨迹上取M ∞个不同的
取样点，这些取样点是均匀分布在轨迹在线，相隔t ?时间的相邻点。由于轨迹是代表所有的平衡分布，而且t ?很小，所以M ∞是一个非常大的值。对一个与时间有关的性质物理量
()x t ，时间平均值可取为1()M k X x k t M ∞∞??=??∑轨迹平均值X ??是我们想要知道的量。
对于在任何时刻的()x t 相对于平均值X ??的变异数为 ()()()221
1M k V X x t S X x k t M ∞∞?-?≡=?-????-∑ (6.51)
取样平均：虽然我们想要知道系统性质的平均值(也就是轨迹平均值)X ??，但是却又不希望对所有轨迹上的M ∞做计算。所以找了一个取代的方法：在轨迹上取一小段，并在此
段上取M 点做计算(M M ∞)。也就是在一小段轨迹上取连续相邻的M 个点，使成一个集合，
此集合称为「样本」。这个样本有它自己的平均值