分子动力学模拟与实践（分子动力学模拟的发展历史-1）

1.1.1分子动力学模拟的概念

分子模拟(molecular modeling 或 molecular simulation)是一类通过计算机 模拟来研究分子或分子体系结构与性质的重要研究方法，包括分子力学 (molecular mechanics, MM)、Monte Carlo (MC)模拟、分子动力学(molecular dynamics, MD)模拟等。这些方法均以分子或分子体系的经典力学模型为基 础，或通过优化单个分子总能量的方法得到分子的稳定构型(MM)；或通过反 复采样分子体系位形空间并计算其总能量的方法，得到体系的最可几构型与热力 学平衡性质(MC)；或通过数值求解分子体系经典力学运动方程的方法得到体 系的相轨迹，并统计体系的结构特征与性质(MD)。目前，得益于分子模拟理 论、方法及计算机技术的发展，分子模拟已经成为继实验与理论手段之后，从分子水平了解和认识世界的第三种手段。

1.1.2 MD模拟的早期历史

最早的MD模拟在1957年就已实现。当时，Alder和Wainwright通过计算机模拟的方法，研究了从32个到500个刚性小球分子系统的运动。模拟开始时, 这些小球分子被置于有序分布的格点上，具有大小相同的速度，但速度方向随机 分布。除相互间的完全弹性碰撞外，刚性小球分子之间没有任何相互作用，小球 分子在碰撞间隙做匀速直线运动。在经过一段时间的模拟，系统中的刚性小球分 子速度达到Maxwell-Boltzmann分布后，他们分另U根据维力定理和径向分布函数 计算了系统的压力，发现两种方法得到的结果一致1959年，他们提出可以 把MD模拟方法推广到更复杂的具有方阱势的分子体系，模拟研究分子体系的 结构和性质図。

1964年，Rahman模拟研究了具有Lennard-Jones势函数的864个Ar原子 体系，得到了与状态方程有关的性质、径向分布函数、速度自相关函数、均方位移等。此后，分子模拟工作者广泛模拟研究了具有不同势函数参数的Len- nard-Jones模型分子体系，得到了体系的结构及其各种热力学性质，探讨了 Lennard-Jones势函数参数对体系结构与性质的影响，建立了 Lennard-Jones势 函数参数与模型分子体系结构及性质之间的关系。

1.1.3分子体系运动方程的数值解

在刚性小球分子体系中，除碰撞瞬间外分子间没有任何相互作用，分子的运 动轨迹由一系列的折线组成。因此，对刚性小球分子体系的MD模拟，算法的 核心是计算刚性小球分子间的碰撞时间及其碰撞前后的运动方向和速度的变化， 而不是直接数值求解刚性小球分子体系的运动方程。相反，对于Lennard-Jones 分子，分子间一直存在相互作用，分子的运动轨迹复杂，必须通过求解体系的运 动方程来确定，这就促进了运动方程数值方法的发展。

单原子分子没有内部结构，计算量小，容易实现MD模拟。相反，多原子 分子具有复杂的内部结构，运动方程更加复杂，计算量大，难以实现MD模拟。 因此，直到20世纪70年代，才实现水分子体系、正烷桂分子体系等多原 子分子体系的MD模拟。

对于多原子分子体系，如果包括化学键伸缩势、键角弯曲势等所有分子内的 相互作用势都由势函数描述，不存在所谓的约束，则体系运动方程的数值解与单 原子分子相同，没有本质区别。常用的MD数值积分方法包括Verlet算法、 预测-校正算法及其由Verlet算法衍生的蛙跳法凹、速度Verlet法〔回等。但是， 如果分子中的化学键等被约束在固定长度或整个分子被约束成刚体，则体系的运 动方程将完全不同。

根据经典力学，没有内部自由度的多原子分子可以被近似为 刚体，运动状态可以分解为质心的平动和刚体的转动两种独立运动模式。其中， 质心的平动与质点力学没有任何区别，不需要另外讨论。刚体的转动通常由Euler角随时间的演化描述，可以通过数值求解相应的运动方程得到。但是，以 Euler角为变量描述刚体运动时，牵涉奇点问题，算法不稳定。为此，Evans提出了以四元数描述刚体运动状态的方法，很好地克服了奇点问题，成为常用的 MD模拟算法小•⑵。

在室温附近，大多数原子、分子处于基态，分子中的化学键长和键角均在平 衡位置附近以很小的幅度振动。因此，假设化学键长和键角在MD模拟过程中 固定不变，不会对MD模拟结果产生显著影响。但是，除键长和键角以外的其 他内部运动自由度，如二面角的旋转.仍被允许，分子不能被近似为质点或刚 体。为了利用MD模拟研究具有固定键长的约束体系，可以采用SHAKE算 法⑷、RATTLE算法等。约束体系的MD模拟，是MD模拟不可缺少的重要 方法，至今仍是非常热门的研究课题。

目前，MD模拟技术已经成熟，不但可以模拟简单的不具内部自由度的单原 子分子和刚性多原子分子，也可以模拟具有内部自由度的多原子分子。即使对于 蛋白质、DNA这样复杂的生物大分子的模拟，也已经没有任何算法上的限制， 只有计算机计算能力的限制。

1. 1.4温度与压力的调控与统计系综的实现

早期MD模拟体系被限制在一定的空间之内，具有固定的体积。此外，根 据能量守恒定律，模拟体系的总能量也被固定，不允许任何波动。从统计力学角 度分析，这样的MD模拟体系是微正则系综或NVE系综。但许多实际的化学、 生物等体系，常与一恒温热源热接触，具有固定温度，但不具有固定的总能量， 这就是正则系综或NVT系综。此外，大多数实际或人工体系，除具有固定温度 外，也具有固定的压力，但具有可变的总能量和体积，这是NPT系综。

最简单的调控体系温度的方法是变标度(scaling)恒温法。当体系的温度， 即总动能偏离设定值时，体系中所有原子或分子的速度被乘以称为标度因子的系 数，使体系的总动能回归设定动能之后，Berendsen改进了变标度恒温法， 使温度逐渐被调节到设定温度，减小了温度的波动幅度氏］。此外，Andersen还 提出了热浴法，在模拟过程中让体系中的一个或若干个原子或分子与恒温热源中 的分子发生随机碰撞，达到调整体系总体温度并使之保持恒定的目的〔闵。但是， 由这样的方法模拟得到的相轨迹不连续，不符合实际情况。

上述温度调控算法粗糙，没有严格的理论基础。具有严格理论基础的温度调 控算法是恒温扩展法(extended system),通过在模拟体系广义坐标和广义动量 外引入额外的自由度与热浴耦合的方法，达到调控温度的目的。理论上，恒温扩 展法可以通过改变或扩展实际模拟体系的哈密顿函数或拉格朗日函数实现，如 Nose恒温扩展法⑵如及其经Hoover修正的Nose-Hoover恒温算法等购。

除了调控温度技术外，压力的调控也是MD模拟的基本方法。与变标度恒 温法相似，最直接的压力调控算法是通过调整模拟体系的体积，达到间接调控体 系压力的目的。但是，如果在各原子的位置坐标直接乘以同一系数，将改变体系 中分子内原子间距或键长，得到错误的结果。因此，调控压力的算法与变标度恒 温法的实现方法不同，需通过复杂的坐标变换，在不改变分子内原子间相对距离 的条件下改变分子的相对位置，实现改变模拟体系体积的目的，这就是标度变换 恒压法。同样，压力的调控也可通过扩展系统自由度，采用恒压扩展法实现，如 Andersen恒压恒熔系综〔闵等。

扩展法的提出是MD理论和方法研究的一个里程碑。目前，通过扩展模拟 体系自由度，在体系的哈密顿函数或拉格朗日函数上添加额外的项，达到调控体 系的温度与压力的方法，仍然是一活跃的研究领域。