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

基于分子动力学的常用力场、算法及结果分析摘要：分子动力学模拟是随着计算机技术发展而兴起的一项新技术，结合了物理，数学和化学的模拟方法。分子动力学模拟具有沟通宏观特性与微观结构的作用，对于许多在理论分析和实验都难以了解和观察到的现象能给出一定的微观解释，可以分析其中的物理现象和本质。本文主要对分子动力学中常用力场，算法及结果分析方法进行综述。最后基于计算机硬件的发展与模拟算法的优化，对该领域未来发展进行展望。关键词：分子动力学；力场；算法；结果分析；引言早在20世纪50年代末，人们就已经开始进行计算模拟实验。一开始，这些计算实验的对象只是一些小的模型体系，比如几个硬球的集合。但早在那时，Alder和Wainwright[1]就已经预见到,这样的计算机模拟实验将会成为联结宏观实验现象和其微观本质的重要桥梁。尽管解析理论进行模拟也可以部分解释实验，但是往往附带了很多近似，而计算机模拟能够提供体系的详尽信息，具有更大的优势。由于宏观和微观在尺度上是截然不同的，为了能够在两者之间建立切实的关联，即微观模拟体系能够反映相应的宏观实验现象，需要通过周期性边界条件对模拟对象体系进行周期性的复制，以避免在实际中并不存在的边缘效应。分子动力学模拟主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。通过分子模拟技术，人们可以在分子、原子尺度上对所关心材料的物理、化学性质以及与之相对应的各种现象进行模拟，摒弃传统材料科学研究中被称之为“爱迪生方式”的炒菜式研究，利用计算机对材料的组成、性质、制备工艺、加工工艺等环节进行虚拟实验，获得大量数据，再与少量验证性实验进行对照，分析物理现象和本质，从而大大的降低材料的研发成本和时间。这种新材料的设计模式已经很大程度上取代了传统的新材料设计方式，成为新材料设计的主要方式之一。2011年，美国启动一项价值超过5亿美元的“先进制造业伙伴关系”计划，而以材料计算与模拟为主的“材料基因组计划”则是其中重要内容之一[2]。目前，对分子动力学模拟的力场、算法及结果分析方法很多，但基本上是分散在英文文献中，由于翻译和理解上的差异，影响了国内从事分子模拟科研工作者的理解交流。本文整理了目前分子动力学中主流的力场、算法及结果分析，并在最后对该领域进行了展望。1力场1.1力场简介分子动力学模拟是计算庞大复杂系统的有效方法，它以力场为依据，力场的完备与否决定计算的可靠程度。分子的总能量是动能与势能之和，分子的势能通常表示为简单的几何坐标的函数。一般势能中包括：（1）范德华力，与能量有关的非键相互作用交叉能量项，（2）构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化，（3）键角变化引起的分子能量变化，（4）单轴旋转引起分子骨架扭曲所产生的能量变化，（5）离平面振动项，共平面原子的中心原子离平面小幅振动的势能，（6）库伦作用项，带电荷粒子间存在的静电吸引或排斥作用的势能。力场可以看作是势能面的经验表达式，它是分子动力学模拟的基础、力场是通过原子位置计算体系能量的，与之前的量子力学方法相比，大大节约了计算时间，可用于计算包含上万粒子数目的体系。势能函数在大多数情况下将描述分子几何形变最大程度地简化为仅仅使用简谐项和三角函数来实现，而非键原子之间的相互作用，则只采用库伦相互作用和兰纳-琼斯势相结合来描述。势能函数的可靠性主要取决于力场参数准确性，而力场参数通常通过拟合实验观测数据和量子力学从头计算得到的数据。目前在生物大分子体系模拟中使用最为广泛的分子力场是CHARMM力场[3]和AMBER力场[4]，也是早期研究生物大分子的分子力场，其现有的力场参数仍在不断优化，并且涵盖的分子类型也在扩大。粗粒化模型在计算生物物理研究中越来越引起人们的关注[5,6]由于该模型中定义了粗粒化粒子，对应于全原子模型中的若干原子或原子基团甚至分子，减少了体系中的粒子数和自由度，使得模拟的时间和空间尺度得以大幅度提高，虽然会丢失一些原子细节信息，但是这种模型是应用于研究缓慢的生物现象或依赖于大组装体的生物现象[7]，如生物膜的波动，对它的模拟需要巨大的膜片。1.2常见力场分子动力计算体系由最初的单原子分子系统延伸至多原子分子、聚合物分子、生化分子系统，力场也随着系统复杂度的增加而增加其复杂性。针对特定的目的，力场分为许多不同的形式，具有不同的使用范围与局限性，在执行分子动力计算时，选择适合的力场极为重要，往往决定了计算成果的优劣与可靠性。下面就一些常用力场做简单介绍。传统力场：AMBER力场[4]：由Kollman课题组开发的力场，是目前使用比较广泛的一种力场，适合处理较小的蛋白质、核酸、多糖等生化分子，此力场参数全部来自计算结果与实验值的比对。CHARMM力场[3]：由Karplus课题组开发，CHARM力场参数除了来自计算结果与实验值的比对外，还引用了大量的量子计算结果为依据。此力场可应用于研究许多分子系统，包括小的有机分子，溶液，聚合物，生化分子等。几乎除了有机金属分子外，通常皆可得到与实验值相近的结果。MMX力场[8,9]：此力场为Allinger等人所发展，依其发展的先后顺序分别为MM2，MM3，MM4，MM+等。MM力场将一些常见的原子细分，如将碳原子分为sp3，sp2，sp，酮基碳，环丙烷碳，碳自由基，碳阳离子等。这些不同形态的碳原子具有不同的力场参数。此力场适用于各种有机化合物，自由基，离子。在MM形式的力场中仔细考虑了许多交叉作用项，其结果往往优于其他形式的力场。相对的，其力场形式较为复杂，比较不易程序化，计算也较费时。CVFF力场：为DauberOsguthope等所发展。此力场最初以生化分子为主，其力场参数适用于胺基酸，水，及各种官能基。其后，经过不断强化，CVFF力场可适用于各种多肽，蛋白质和大量有机分子。第二代力场：第二代的势能函数形式比传统力场要更加复杂，涉及的力场参数更多，计算量也更大，当然也相应地更加准确。CFF力场：CFF力场是一个力场家族，包括了CFF91、PCFF、CFF95等很多力场，可以进行从有机小分子、生物大分子到分子筛等诸多体系的计算。COMPASS力场[10]：由MSI公司开发的力场，擅长进行高分子体系的计算。MMF94力场：Hagler开发的力场，是目前最准确的力场之一。SFF力场：SI公司开发的力场，可以进行有机、无机分子的计算。UFF力场：可以计算周期表上所有元素的参数。Dreiding力场[11]：用于有机小分子、大分子、主族元素的计算。2分子动力学关键算法分子动力学发展到现在成为一种稳健的理论工具，并因此成为许多实验方法必不可少的补充，离不开前人在热力学平衡模拟和非平衡模拟方面所做的大量工作。有两个算法非常值得我们回顾，那就是Nosé-Hoover恒温算法[12,13]和Nosé-Andersen恒压算法[14]。这两个算法分别实现了正则体系的模拟以及恒压恒焓体系的模拟。二者的结合，实现了目前常用的恒温恒压体系的模拟，而恒温恒压体系最接近大多数实验条件。在这些算法中，格点加和算法[15]旨在处理长程静电相互作用。而与其相关的，是一种使用不同时间步长[16,17]模拟运动演化的数值方法，即利用Trotter因子分解，将力场中短程力和长程力的贡献分开处理，对于运动较快的和较慢的自由度也采用分别处理的方法。这样一来，对于运动方程进行积分时就可以采用不同的时间步长。与多时间步长方法紧密相关的另一种方法是利用完整约束[18]，去除体系中一些较“硬”（具有较高振动频率）的自由度，从而可增加时间步长。例如，可以通过Shake和Rattle算法[18]将含有氢原子的化学键保持在其平衡值不变，以消除其振动自由度，这样就可以使用较大的积分步长，而又不致影响到系统总能量的守恒。这一方法的根本目的是在不增加计算消耗的情况下能够获得更多的采样。但是，由于增加了额外的计算求解约束方程，总体上所提高的计算效率并不理想。但是，在人们孜孜以求不断增大时间和空间尺度的过程中，真正革命性的成就莫过于并行架构和空间分解算法[19]的出现，使得计算时间可以随着处理器数目的增加而线性地减少。然而，这一伟大的成就并没有让人们追求最长模拟时间或最大模拟体系的欲望得到彻底的满足，反而变得更加强烈，同时，在一定程度上也遮住了在改进原子间相互作用力的表达方面所取得的辉煌成就的光芒，特别是极化率[20]和分布式多极矩[21]的引入，当然这也相应地增加了计算的开销，但这是为了提高计算精度而必须付出的代价，另一方面也多少忽视了在表征运动方程数值积分的误差方面所付出的巨大努力。分子动力学模拟要求按一定的时间步长对经典运动方程进行不连续的积分，而时间步长受制于分子体系中运动最快的自由度[22]，在对生物体系的全原子描述中，由于要描述包括氢键在内的化学键的振动，为了确保体系能量守恒，时间步长只能在飞秒数量级。因此，需要进行数百万甚至数十亿次积分才能在与生物体系相应的时间尺度上描述体系的行为。如前所述，尽可能减少计算能量时的开销是非常重要的。另一方面，在对势加和的方法中，理论上的计算开销正比于N2，即体系中粒子数的平方。在实际计算中[15]，对静电相互作用中的短程力部分可采用球形截断策略，而求解Poisson方程时则可采用离散化的处理方案，这样，就可以在倒易空间里处理长程相互作用，从而使实际的计算开销减少到NlogN。在利用大型中央处理器CPU阵列运行大规模并行分子动力学程序的做法出现之后大约二十年，数值模拟领域又成为了另一场革新的舞台，这一革新是由通用图形处理器GPU的出现而引发的。如今，通用图形处理器已经变成一种并不昂贵的多核通用处理器，它能够以并行的方式处理浮点数据运算，而制造成本的降低得益于其在消费类电子产品，尤其是在视频游戏控制平台方面的广泛应用。曾经风靡一时的程序，如NAMD[23]和GROMACS[24]等都很快作出了调整以适应这种新型的异构化并行架构[25]，同时，新的专门针对图形处理器而设计的分子动力学代码如HOOMD-blue[26]等也应运而生。3结果分析方法3.1常用宏观统计分析方法对分布函数、径向分布函数和静态结构因子都是表征体系结构的宏观统计参量，三者之间可以进行数值的相互转化。由于静态结构因子也可以通过X衍射或中子衍射在实验中得到，因此对分布函数、径向分布函数和静态结构因子这3种方法就可以作为连接计算机模拟与实验分析的桥梁，成为分子动力学模拟中最常使用的结构分析方法。对分布函数（PDF）反映的是从一个任意指定的“中心”粒子出发，到半径为ｒ的位置上出现其它粒子的几率（单位体积内的粒子数目）[27]。如图1所示，取一原子为观测中心，以1Å为单位向外画出一层层的同心球壳，然后统计每一壳层里的原子数密度（粒子数n(r)/体积V）与平均数密度（ρ0＝总粒子数/体积V）比值，即为对分布函数，常用g(r)表示。图1对分布函数与原子结构关系示意图对于g(r)的分析，主要是基于以下几点：峰的分布规律。根据对分布函数的物理含义，对于晶体，由于原子的周期性排列，每层原子会对应出现分布峰值，而在层与层的原子之间出现原子的几率为0。对于液体结构，由于短程有序，在第一近邻处会出现较为尖锐的峰值，在中远程的范围里，液体结构的对分布函数会出现漫高峰，并且趋近于1。这说明在液态结构中，相对于一个中心原子，在无穷远处会总会有一个原子的存在。第一峰的高度和形状。第一峰代表第一近邻原子之间的结合强度，如果它的外形比较尖锐，表明在此半径范围内的原子数密度比平均密度要高很多，中心原子与最近邻原子的相互结合强度也比较大。第一峰与第一谷的关系。Abraham[28]等人提出一种经验性方法，根据对分布函数中第一谷值gmin与第一峰值gmax的比例，定义参量R=gmin/gmax，用来确定玻璃转变点Tg，如图2(a)所示。第二峰的形状与分布。第二峰代表着中程序的连接强度，如果第一峰和第二峰比较明显，同时第一峰与第二峰之间有一个很深的波谷，基本上表