基于分子动力学的常用力场、算法及结果分析

1、基于分子动力学的常用力场、算法及结果分析摘要：分子动力学模拟是随着计算机技术发展而兴起的一项新技术，结合了物理，数学和化学的模拟方法。分子动力学模拟具有沟通宏观特性与微观结构的作用，对于许多在理论分析和实验都难以了解和观察到的现象能给出一定的微观解释，可以分析其中的物理现象和本质。本文主要对分子动力学中常用力场，算法及结果分析方法进行综述。最后基于计算机硬件的发展与模拟算法的优化，对该领域未来发展进行展望。关键词：分子动力学；力场；算法；结果分析；引言早在20世纪50年代末，人们就已经开始进行计算模拟实验。一开始，这些计算实验的对象只是一些小的模型体系，比如几个硬球的集合。但早在那时，Alde

2、r和Wainwright1就已经预见到, 这样的计算机模拟实验将会成为联结宏观实验现象和其微观本质的重要桥梁。尽管解析理论进行模拟也可以部分解释实验，但是往往附带了很多近似，而计算机模拟能够提供体系的详尽信息，具有更大的优势。由于宏观和微观在尺度上是截然不同的，为了能够在两者之间建立切实的关联，即微观模拟体系能够反映相应的宏观实验现象，需要通过周期性边界条件对模拟对象体系进行周期性的复制，以避免在实际中并不存在的边缘效应。分子动力学模拟主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和

3、其他宏观性质。通过分子模拟技术，人们可以在分子、原子尺度上对所关心材料的物理、化学性质以及与之相对应的各种现象进行模拟，摒弃传统材料科学研究中被称之为“爱迪生方式”的炒菜式研究，利用计算机对材料的组成、性质、制备工艺、加工工艺等环节进行虚拟实验，获得大量数据，再与少量验证性实验进行对照，分析物理现象和本质，从而大大的降低材料的研发成本和时间。这种新材料的设计模式已经很大程度上取代了传统的新材料设计方式，成为新材料设计的主要方式之一。2011年，美国启动一项价值超过5亿美元的“先进制造业伙伴关系”计划，而以材料计算与模拟为主的“材料基因组计划”则是其中重要内容之一2。目前，对分子动力学模拟的力场

4、、算法及结果分析方法很多，但基本上是分散在英文文献中，由于翻译和理解上的差异，影响了国内从事分子模拟科研工作者的理解交流。本文整理了目前分子动力学中主流的力场、算法及结果分析，并在最后对该领域进行了展望。1 力场1.1 力场简介分子动力学模拟是计算庞大复杂系统的有效方法，它以力场为依据，力场的完备与否决定计算的可靠程度。分子的总能量是动能与势能之和，分子的势能通常表示为简单的几何坐标的函数。一般势能中包括：（1）范德华力，与能量有关的非键相互作用交叉能量项，（2）构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化，（3）键角变化引起的分子能量变化，（4）单轴旋转引起分子骨架扭曲所产生的

5、能量变化，（5）离平面振动项，共平面原子的中心原子离平面小幅振动的势能，（6）库伦作用项，带电荷粒子间存在的静电吸引或排斥作用的势能。力场可以看作是势能面的经验表达式，它是分子动力学模拟的基础、力场是通过原子位置计算体系能量的，与之前的量子力学方法相比，大大节约了计算时间，可用于计算包含上万粒子数目的体系。势能函数在大多数情况下将描述分子几何形变最大程度地简化为仅仅使用简谐项和三角函数来实现，而非键原子之间的相互作用，则只采用库伦相互作用和兰纳-琼斯势相结合来描述。势能函数的可靠性主要取决于力场参数准确性，而力场参数通常通过拟合实验观测数据和量子力学从头计算得到的数据。目前在生物大分子体系模拟

6、中使用最为广泛的分子力场是CHARMM力场3和AMBER力场4，也是早期研究生物大分子的分子力场，其现有的力场参数仍在不断优化，并且涵盖的分子类型也在扩大。粗粒化模型在计算生物物理研究中越来越引起人们的关注5, 6由于该模型中定义了粗粒化粒子，对应于全原子模型中的若干原子或原子基团甚至分子，减少了体系中的粒子数和自由度，使得模拟的时间和空间尺度得以大幅度提高，虽然会丢失一些原子细节信息，但是这种模型是应用于研究缓慢的生物现象或依赖于大组装体的生物现象7，如生物膜的波动，对它的模拟需要巨大的膜片。1.2常见力场分子动力计算体系由最初的单原子分子系统延伸至多原子分子、聚合物分子、生化分子系统，力场

7、也随着系统复杂度的增加而增加其复杂性。针对特定的目的，力场分为许多不同的形式，具有不同的使用范围与局限性，在执行分子动力计算时，选择适合的力场极为重要，往往决定了计算成果的优劣与可靠性。下面就一些常用力场做简单介绍。传统力场：AMBER力场4：由Kollman课题组开发的力场，是目前使用比较广泛的一种力场，适合处理较小的蛋白质、核酸、多糖等生化分子，此力场参数全部来自计算结果与实验值的比对。CHARMM力场3：由Karplus课题组开发，CHARM力场参数除了来自计算结果与实验值的比对外，还引用了大量的量子计算结果为依据。此力场可应用于研究许多分子系统，包括小的有机分子，溶液，聚合物，生化分子

8、等。几乎除了有机金属分子外，通常皆可得到与实验值相近的结果。MMX力场8, 9：此力场为Allinger等人所发展，依其发展的先后顺序分别为MM2，MM3，MM4，MM+等。MM力场将一些常见的原子细分，如将碳原子分为sp3，sp2，sp，酮基碳，环丙烷碳，碳自由基，碳阳离子等。这些不同形态的碳原子具有不同的力场参数。 此力场适用于各种有机化合物，自由基，离子。在MM形式的力场中仔细考虑了许多交叉作用项，其结果往往优于其他形式的力场。相对的，其力场形式较为复杂，比较不易程序化，计算也较费时。CVFF力场：为Dauber Osguthope等所发展。此力场最初以生化分子为主，其力场参数适用于胺基

9、酸，水，及各种官能基。其后，经过不断强化，CVFF力场可适用于各种多肽，蛋白质和大量有机分子。第二代力场：第二代的势能函数形式比传统力场要更加复杂，涉及的力场参数更多，计算量也更大，当然也相应地更加准确。CFF力场：CFF力场是一个力场家族，包括了CFF91、PCFF、CFF95等很多力场，可以进行从有机小分子、生物大分子到分子筛等诸多体系的计算。COMPASS力场10：由MSI公司开发的力场，擅长进行高分子体系的计算。MMF94力场：Hagler开发的力场，是目前最准确的力场之一。SFF力场：SI公司开发的力场，可以进行有机、无机分子的计算。UFF力场：可以计算周期表上所有元素的参数。Dre

10、iding力场11：用于有机小分子、大分子、主族元素的计算。2 分子动力学关键算法分子动力学发展到现在成为一种稳健的理论工具，并因此成为许多实验方法必不可少的补充，离不开前人在热力学平衡模拟和非平衡模拟方面所做的大量工作。有两个算法非常值得我们回顾，那就是Nosé-Hoover恒温算法12, 13和Nosé-Andersen恒压算法14。这两个算法分别实现了正则体系的模拟以及恒压恒焓体系的模拟。二者的结合，实现了目前常用的恒温恒压体系的模拟，而恒温恒压体系最接近大多数实验条件。在这些算法中，格点加和算法15旨在处理长程静电相互作用。而与其相关的，是一种使用不同时间步长16,

11、 17模拟运动演化的数值方法，即利用Trotter因子分解，将力场中短程力和长程力的贡献分开处理，对于运动较快的和较慢的自由度也采用分别处理的方法。这样一来，对于运动方程进行积分时就可以采用不同的时间步长。与多时间步长方法紧密相关的另一种方法是利用完整约束18，去除体系中一些较“硬”（具有较高振动频率）的自由度，从而可增加时间步长。例如，可以通过Shake和Rattle算法18将含有氢原子的化学键保持在其平衡值不变，以消除其振动自由度，这样就可以使用较大的积分步长，而又不致影响到系统总能量的守恒。这一方法的根本目的是在不增加计算消耗的情况下能够获得更多的采样。但是，由于增加了额外的计算求解约束

12、方程，总体上所提高的计算效率并不理想。但是，在人们孜孜以求不断增大时间和空间尺度的过程中，真正革命性的成就莫过于并行架构和空间分解算法19的出现，使得计算时间可以随着处理器数目的增加而线性地减少。 然而，这一伟大的成就并没有让人们追求最长模拟时间或最大模拟体系的欲望得到彻底的满足，反而变得更加强烈，同时，在一定程度上也遮住了在改进原子间相互作用力的表达方面所取得的辉煌成就的光芒，特别是极化率20和分布式多极矩21的引入，当然这也相应地增加了计算的开销，但这是为了提高计算精度而必须付出的代价，另一方面也多少忽视了在表征运动方程数值积分的误差方面所付出的巨大努力。分子动力学模拟要求按一定的时间步长

13、对经典运动方程进行不连续的积分，而时间步长受制于分子体系中运动最快的自由度22，在对生物体系的全原子描述中，由于要描述包括氢键在内的化学键的振动，为了确保体系能量守恒，时间步长只能在飞秒数量级。因此，需要进行数百万甚至数十亿次积分才能在与生物体系相应的时间尺度上描述体系的行为。如前所述，尽可能减少计算能量时的开销是非常重要的。另一方面，在对势加和的方法中，理论上的计算开销正比于N2，即体系中粒子数的平方。在实际计算中15，对静电相互作用中的短程力部分可采用球形截断策略，而求解Poisson方程时则可采用离散化的处理方案，这样，就可以在倒易空间里处理长程相互作用，从而使实际的计算开销减少到N l

14、ogN。在利用大型中央处理器CPU阵列运行大规模并行分子动力学程序的做法出现之后大约二十年，数值模拟领域又成为了另一场革新的舞台，这一革新是由通用图形处理器GPU的出现而引发的。如今，通用图形处理器已经变成一种并不昂贵的多核通用处理器，它能够以并行的方式处理浮点数据运算，而制造成本的降低得益于其在消费类电子产品，尤其是在视频游戏控制平台方面的广泛应用。曾经风靡一时的程序，如NAMD23和GROMACS24等都很快作出了调整以适应这种新型的异构化并行架构25，同时，新的专门针对图形处理器而设计的分子动力学代码如HOOMD-blue26等也应运而生。3 结果分析方法3.1 常用宏观统计分析方法对分

15、布函数、径向分布函数和静态结构因子都是表征体系结构的宏观统计参量，三者之间可以进行数值的相互转化。由于静态结构因子也可以通过X衍射或中子衍射在实验中得到，因此对分布函数、径向分布函数和静态结构因子这3种方法就可以作为连接计算机模拟与实验分析的桥梁，成为分子动力学模拟中最常使用的结构分析方法。对分布函数（PDF）反映的是从一个任意指定的“中心”粒子出发，到半径为的位置上出现其它粒子的几率（单位体积内的粒子数目）27。如图1所示，取一原子为观测中心，以1 Å为单位向外画出一层层的同心球壳，然后统计每一壳层里的原子数密度（粒子数n(r)/体积V）与平均数密度（0总粒子数/体积V）比值，即为

16、对分布函数，常用g(r)表示。图1 对分布函数与原子结构关系示意图对于g(r)的分析，主要是基于以下几点：峰的分布规律。根据对分布函数的物理含义，对于晶体，由于原子的周期性排列，每层原子会对应出现分布峰值，而在层与层的原子之间出现原子的几率为0。对于液体结构，由于短程有序，在第一近邻处会出现较为尖锐的峰值，在中远程的范围里，液体结构的对分布函数会出现漫高峰，并且趋近于1。这说明在液态结构中，相对于一个中心原子，在无穷远处会总会有一个原子的存在。第一峰的高度和形状。第一峰代表第一近邻原子之间的结合强度，如果它的外形比较尖锐，表明在此半径范围内的原子数密度比平均密度要高很多，中心原子与最近邻原子的

17、相互结合强度也比较大。第一峰与第一谷的关系。Abraham28等人提出一种经验性方法，根据对分布函数中第一谷值gmin与第一峰值gmax的比例，定义参 量R=gmin/gmax，用来确定玻璃转变点Tg，如图2(a)所示。第二峰的形状与分布。第二峰代表着中程序的连接强度，如果第一峰和第二峰比较明显，同时第一峰与第二峰之间有一个很深的波谷，基本上表明原子是以一种共价键的形式存在。如果第二峰出现了劈裂，则这是出现非晶结构的典型特征29, 30。图2 PDF参量、原子体积与温度关系径向分布函数（RDF）表示在半径为r处厚度r的球壳内的平均原子数，其数值等于该位置出现粒子的几率g(r)与球壳面积之积。常

18、用符号G(r)表示。由于语言转换和理解上的差异，径向分布函数与对分布函数在中外文的文献中经常被互用，有时使用的是径向分布函数G(r)概念，表达出来的却是对分布函数g(r)含义。其实二者无论是在图像上还是物理意义上都有明确的区别31, 32。如图所示，对分布函数g(r)表示的是概率，对于液体或非晶结构，分布函数曲线随着半径r的增加应该是归一的，如图3(a)；而径向分布函数G(r)则表示在球面上粒子的统计值，曲线随着半径的增加整体上升，如图3(b)。图3 对分布函数与径向分布函数的关系静态结构因子（SSF）是一个连接实验分析与模拟分析的重要参量，对于衍射分析来说，它表征的是材料对射线的散射能力，反

19、映结构的平均信息。静态结构因子(S(k)-1)与对分布函数(g(r)-1)互为傅里叶变换，在粒子的空间结构信息中，静态结构因子为倒易空间，对分布函数为实空间。峰位置的相对关系。由于结构因子是倒易空间，变量k与距中心原子的半径r成反比关系。根据第二峰位与第一峰位k2/k1的比值，可以将液体分为3种类型。第一种为大部分的液态金属，比值为1.86；第二种为Ga, Sn, Bi等液体，比值为1.96，其特征为在第一峰右边存在一个肩峰；第三种为In和Tl，比值为1.8833；预峰的特征。在第一峰的前面，如果出现小的预峰，反映的是体系中存在着中程序结构特征33；熔点温度峰值特征。通过对多种液体的统计发现，

20、液态结构因子的第一峰高度在熔点附近都会达到一个同样的高度，即S(k) 2.8，这个规律叫做Hansen-Verlet凝固判据，是可以用来作为熔化和凝固的判据之一34。3.2 常用共近邻分析方法共近邻分析（CNA）是用来描述一个原子周围近邻原子的排列特征的方法，根据计算方法的不同，可包括键对分析、键角分析、泰森多边形指数分析、配位数等。最近邻的距离一般是采用对分布函数的第一低谷处的原子位置值rmin。键对分析是研究液体和固体中团簇的一种常用的有效方法35。这种方法用4个指数i, j, l, m对原子进行分类标定。i=1表示研究的两个原子对能够成键，i=2表示两个原子对不能成键，一般我们都关注于成

21、键原子之间的关系。j表示两个被研究原子对中相同的近邻原子数。1则表示j个相同的近邻原子中成键的原子对数。m表示l中成键近邻原子的排布情况，如果成键相连，则m为1，否则m为236。键角分布就是体系中所有原子与其它最近邻原子成键所构成的键角的分布，所获取的是原子在三维空间上的排布信息，因而可以弥补对分布函数g(r)仅仅依据原子距离结构获取信息的不足。然而由于温度的作用，原子之间的振动会影响实际的角度分布情况，因此从角度的分布变化上，可以看出体系中的二十面体或四面体的扭曲和破坏情况37。泰森多边形是一种几何算法，在气象、物流、商业等多个行业和领域中已被广泛利用，也常被称为Voronoi网格。它是由两

22、个邻点连接直线的垂直平分线组成的连续多边形组成，二维泰森多边形原理示意图如图所示。这种方法可以描述在平面上的N个点，按照最邻近原则划分与它的最近邻区域的相互关系。利用Voronoi算法，将离散的原子进行三维空间上的划分，利用所得到的空间结构来表征材料的微观结构，目前也已成为结构表征的一种常用手段。泰森多边形指数常用<n3, n4, n5, n6>的形式表示，其中ni代表着Voronoi 多面体具有的i边形数38。第一近邻的配位数也是一种传统的分析微观结构的方式，此概念首先由阿尔弗雷德·维尔纳提出，原来更多地用于化学结构的表征，是配位化学的基础。在微观结构的表征中，配位数是

23、用来描述中心原子第一壳层内原子的平均数目，反映的是中心原子与其它原子的结合能力和配位关系，描述的是体系中粒子排列的紧密程度，配位数越大，粒子排列越紧密。一般来说配位数是通过对径向分布函数的第一峰进行积分来得到。对于Voronoi指数分析而言，由于泰森多边形的算法就是计算中心原子与配位原子之间的连线平分面，因此，通过对Voronoi指数的累加，也能得到中心原子的配位数。对于晶体结构，通过配位数可以判断出晶体的结构；对于液态的非晶态结构，配位数可以作为一个发生结构转变的敏感参量，为结构转变的判断提供依据。例如，对于液体中发生的液液相变现象，众多研究39都将配位数的变化作为判断的依据之一。4 展望如

24、果分子体系足够小，我们又有足够的耐心，即便只使用为数不多的处理器，完成微秒甚至是毫秒级时间尺度的模拟已不成问题，几年前就已经实现。然而，只有大规模并行架构，特别是千兆级运算速度的超级计算机的出现, 我们才得以处理数百万个原子组成的分子体系。但现阶段，毫秒级时间尺度或者百万原子空间尺度的分子模拟仍然构成技术上的巨大挑战，这恰恰又昭示着一个新时代的来临。随着第一台千兆级浮点运算速度的超级计算机投入使用，科学家已经开始期待着亿亿级计算机的到来并开始考虑这种新型架构将能够处理的分子体系的种类。假使在未来的几十年，分子动力学的发展仍然紧跟摩尔定律所预测的计算机发展趋势，van Gunsteren40乐观

25、地预测，未来20年内，我们将可以在纳秒级时间尺度上模拟一个完整的细菌，如大肠杆菌，大约20年后，可以模拟一个完整的哺乳动物细胞。这些基于目前硬件性能和软件发展水平所作出的推测固然很振奋人心，但我们应该清楚地认识到横亘在模拟时间尺度和空间尺度之间的鸿沟依然难以逾越。在未来的计算机模拟中，利用直接无偏分子动力学模拟研究生物大体系，有可能提供时间相关和动力学相关信息。而实现这一目标，我们仍然还有很长的路要走。参考文献：1 Alder B J, Wainwright T E. Phase transition for a hard sphere systemJ. The Journal of chem

26、ical physics, 1957, 27(5): 1208.2 赵继成. 材料基因组计划简介J. 自然杂志, 2014, 36(2): 89-104.3 Brooks B R, Bruccoleri R E, Olafson B D, et al. CHARMM: a program for macromolecular energy, minimization, and dynamics calculationsJ. Journal of computational chemistry, 1983, 4(2): 187-217.4 Weiner S J, Kollman P A, Cas

27、e D A, et al. A new force field for molecular mechanical simulation of nucleic acids and proteinsJ. Journal of the American Chemical Society, 1984, 106(3): 765-784.5 Izvekov S, Voth G A. A multiscale coarse-graining method for biomolecular systemsJ. The Journal of Physical Chemistry B, 2005, 109(7):

28、 2469-2473.6 Marrink S J, Risselada H J, Yefimov S, et al. The MARTINI force field: coarse grained model for biomolecular simulationsJ. The Journal of Physical Chemistry B, 2007, 111(27): 7812-7824.7 Klein M L, Shinoda W. Large-scale molecular dynamics simulations of self-assembling systemsJ. Scienc

29、e, 2008, 321(5890): 798-800.8 Allinger N L, Yuh Y H, Lii J H. J Am Chem Soc 111: 8551;(b) Lii JHJ. Allinger NL (1989) J Am Chem Soc, 1989, 111: 8566.9 Allinger N L, Chen K, Lii J H. An improved force field (MM4) for saturated hydrocarbonsJ. Journal of computational chemistry, 1996, 17(56): 642-668.1

30、0 Sun H. COMPASS: an ab initio force-field optimized for condensed-phase applications overview with details on alkane and benzene compoundsJ. The Journal of Physical Chemistry B, 1998, 102(38): 7338-7364.11 Rappé A K, Casewit C J, Colwell K S, et al. UFF, a full periodic table force field for m

31、olecular mechanics and molecular dynamics simulationsJ. Journal of the American chemical society, 1992, 114(25): 10024-10035.12 Nosé S. A molecular dynamics method for simulations in the canonical ensembleJ. Molecular physics, 1984, 52(2): 255-268.13 Hoover W G. Canonical dynamics: equilibrium

32、phase-space distributionsJ. Physical review A, 1985, 31(3): 1695.14 Andersen H C. Molecular dynamics simulations at constant pressure and/or temperatureJ. The Journal of chemical physics, 1980, 72(4): 2384-2393.15 Darden T, York D, Pedersen L. Particle mesh Ewald: An N log (N) method for Ewald sums

33、in large systemsJ. The Journal of chemical physics, 1993, 98(12): 10089-10092.16 Tuckerman M, Berne B J, Martyna G J. Reversible multiple time scale molecular dynamicsJ. The Journal of chemical physics, 1992, 97(3): 1990-2001.17 Izaguirre J A, Reich S, Skeel R D. Longer time steps for molecular dyna

34、micsJ. The Journal of chemical physics, 1999, 110(20): 9853-9864.18 Miyamoto S, Kollman P A. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water modelsJ. Journal of computational chemistry, 1992, 13(8): 952-962.19 Brown D, Clarke J H R, Okuda M, et al. A domain decomposit

35、ion parallelization strategy for molecular dynamics simulations on distributed memory machinesJ. Computer Physics Communications, 1993, 74(1): 67-80.20 Chipot C, Ángyán J G. Continuing challenges in the parametrization of intermolecular force fields. Towards an accurate description of elec

36、trostatic and induction termsJ. New journal of chemistry, 2005, 29(3): 411-420.21 Sagui C, Pedersen L G, Darden T A. Towards an accurate representation of electrostatics in classical force fields: Efficient implementation of multipolar interactions in biomolecular simulationsJ. The Journal of chemic

37、al physics, 2004, 120(1): 73-87.22 Nanoscience: Nanobiotechnology and nanobiologyM. Springer Science & Business Media, 2009.23 Plimpton S. Fast parallel algorithms for short-range molecular dynamicsJ. Journal of computational physics, 1995, 117(1): 1-19.24 Hess B, Kutzner C, Van Der Spoel D, et

38、al. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulationJ. Journal of chemical theory and computation, 2008, 4(3): 435-447.25 Stone J E, Phillips J C, Freddolino P L, et al. Accelerating molecular modeling applications with graphics processorsJ. Journal of com

39、putational chemistry, 2007, 28(16): 2618-2640.26 Anderson J A, Lorenz C D, Travesset A. General purpose molecular dynamics simulations fully implemented on graphics processing unitsJ. Journal of Computational Physics, 2008, 227(10): 5342-5359.27 March N H, Tosi M P, March N H. Introduction to liquid

40、 state physicsM. Singapore: World Scientific, 2002.28 Abraham F F. An isothermalisobaric computer simulation of the supercooledliquid/glass transition region: Is the shortrange order in the amorphous solid fcc?J. The journal of chemical physics, 1980, 72(1): 359-365.29 Cheng Y Q, Ma E. Atomic-level

41、structure and structureproperty relationship in metallic glassesJ. Progress in Materials Science, 2011, 56(4): 379-473.30 Qi L, Dong L F, Zhang S L, et al. Glass formation and local structure evolution in rapidly cooled Pd Ni alloy melt under high pressureJ. Physics Letters A, 2008, 372(5): 708-711.31 Kim T H, Kelton K F. Structural study of supercooled liquid transition metalsJ. The Journal of chemical physics, 2007, 126(5): 054513.32 Lai