verlet算法.doc

从计算数学的观点来看，MD方法是一个初值问题，主要用有限差分方法来求解。能够求解Newton运动方程的有限差分算法很多，但分子动力学计算原子数目较多，所以对时间和空间要求都比较高，需要研究适合于分子动力学的时间积分算法。
时间积分算法基于有限差分法，对时间进行离散，为时间步长，由时刻的变量及其时间导数推演得到+时刻的物理量，逐步计算即得到整个历史过程的物理量。这个过程是近似的，不可避免地引入误差，包括截断误差和舍入误差。前者主要是有限差分算法采用Taylor展开截断某些高阶项引入的，后者是计算机本身数值精度引入的。两种误差都能通过减小时间步长而减小。对于比较大的，截断误差起主导作用；但随着下降截断误差迅速降低，舍入误差随下降变化不大，并在较小的时起主导作用。分子动力学时间步长一般随模拟的材料、温度、物理过程而变化，一般为fs (秒)量级。
进行数值求解的算法有多种，其共同的出发点是，将粒子的位置和其它力学量如速度、加速度等展开为Taylor级数。最常用的数值求解算法有Verlet算法，leap-frog算法，Beeman算法及Gear所提出的校正预测法。
1.Verlet算法
Verlet提出的Verlet算法在分子动力学中运用最为广泛，也是最简单的。它运用原子在时刻的位置和加速度及时刻的位置，计算出时刻的位置。将粒子的位置以Taylor公式展开，即

将式中的换为-，得

两式相减得速度式

两式相加得

因，故依据上式可由及的位置预测时的位置。
Verlet算法执行简明，需要的内存小，但其缺点在于速度式中含有项，而实际计算中通常选取很小的值，且位置要通过小项与非常大的两项与的差的相加得到，容易造成精度损失。另外，它不是一个自启动算法，新位置必须由时刻与前一时刻的位置得到。在t=O时刻，只有一组位置，所以必须通过其他方法得到的位置。获得时刻的位置的方法之一时应用近似式。

2.跳蛙法（leap-frog method）
鉴于Verlet算法的一些缺点，Hockney提出了跳蛙法。跳蛙法是从Verlet法推导出来的。它在半个积分时间步得到速度，并利用这一速度计算新的位置。位置和速度表达式为

计算时假设已知与，由时的位置计算质点所受的力与加速度，再预测时间为时的速度，以此类推。
根据与，可得时间为时的速度为

利用跳蛙法计算仅需储存与两类信息，既节约储存空间，而且准确性及稳定性较高。跳蛙法不需要计算下一步位置就可以得出速度，但需要注意的是速度并未与位置在同一时间定义，结果是动能和势能也未同时定义。所以不能直接计算总能量。
3.Beeman法
Beeman法是除Verlet跳蛙法外，另一种较为常见的方法。其积分公式如下

此方法需储存、与，储存量大于Verlet的跳蛙法。但优点在于可以引用较长的积分间隔。Beeman方法所引用的积分步长可为Verlet方法的倍，而具有相同的准确性。
4.校正预测法（predictor-corrector method）
Gaer提出了基于预测-校正积分方法的校正预测法。因为经典运动粒子的轨迹为连续的，故于时间时的位置、速度等可由时间的Taylor展开式预测得

由于这些物理量来自Taylor展开式，并非由解Newton运动方程式而来，所以式（2.17）中所产生的速度、加速度等并非完全正确。为了解决这个问题，用所预测的位置计算所受的力及正确的加速度。设正确的加速度与预测的加速度之间的误差为

得各量之校正式为

式中，、、、均为常数。这仅为Gear的一次预测校正法，也可以推展至更高次的校正法。对于大多数MD应用来说，Verlet型的算法就可以满足基本的要求了。但有时采用高阶算法更方便。一般来说Gear算法中的Taylor展开阶数越高，精度越好。但占用内存比较大，这对于大量分子系统的模拟不如Verlet算法方便。当减小时，Gear的高次预测校正法误差较小，但增大时，Verlet方法较佳。

边界条件
执行分子动力计算通常选取一定数目N的分子，将其置于一立方体的盒子中。设盒子的边长为，则其体积为=。若分子的质量为，则系统的密度为：

计算系统的密度应等于实验所测定的密度。
介质的宏观性质是由大量粒子所为，所以在模拟中，粒子的数量只有足够多，才能真实地再现介质的宏观性质。然而由于计算机的计算能力的限制，我们在模拟时粒子数一般应少于10000个。同时粒子数量越少，表面效应越明显，因为此时大量的粒子是处于表面的。这样少的粒子组成的集团和真实体系相比，其表面原子与体内原子数之比显然过大，这必然造成表面效应。无论模拟的区域是否在容器中，表面的粒子受力仍不同于实际情况，从而很难正确的描述其宏观性质。为了减小有限尺寸的影响，在分子动力学模拟中采用了各种有效的边界条件。通常用于分子动力学模拟中的边界条件有两种:周期性边界条件、非周期性边界条件。
（1）模拟计算中，为使计算中系统的密度维持恒定，通常采用周期性边界条件[2]，以二维的计算系统为例，图2.3显示了二维盒中系统粒子的排列及移动方向。

图2.3 二维周期性系统的例子排列与移动 图2.4 粒子的最近映像
图中位于中央的盒子表示所模拟的系统，其周围盒子与模拟系统具有相同的排列及运动。当计算系统中任一粒子移出盒外，则必有一粒子由相对的方向移入。如图2.3中的第2个粒子。这样的限制条件使得系统中的粒子数维持恒定，密度不变，符合实际的要求。
计算系统中分子间的作用力时，采取最近镜像方法。如图2.4(粒子的最近映像)所示，计算分子2与5的作用力，是取与分子2和其最近的距离镜像分子5。因为在计算中利用最近镜像的观念，因此需采用截断半径的方法计算非键结的远程作用力，否则会因重复计算粒子所受的力而导致不正确的结果。图2.5为一般的van der Waals势能图。

图2.5 van der Waals势能图
图中显示，当>时，势能的值已非常趋近于零，表示分子间的van der Waals作用力已可以忽略不计。在很多分子动力学模拟过程中，若分子间的距离大于截断半径，则将其作用视为零。截断半径最大不能超过盒长的一半，即。一般的原子所选取的截断半径约为。
（2）分子模拟并不都是周期边界，在许多情况下，如对于液滴体系，溶液中沉淀的分子团簇，并不需要周期边界。如图2.7所示，图中的外壳，即阴影部分是体系的边界区。边界区内的分子可以设计为位置固定的。体系中部的分子是可以移动的。

图2.7 非周期边界情形
时间步长的选取
时间步长与采用的数值积分方法、研究的对象有关。在选择了数值算法后，时间步长的选取是影响计算结果的一个关键因素。这里存在一对矛盾，一方面应尽量选择大的时间步以模拟更长的时间，从而在相空间抽样的比例更大；另一方面，过大的时间步将导致结果精度降低、甚至分子动力学的崩溃。
时间步长与所研究现象的弛豫时间密切相关。一般地，时间步长可安全地取为弛豫时间的。对于金属单原子晶体结构，时间步长可大致取为10sf。但模拟系统可能存在几个不同的时间尺度，如对于分子系统，分子间和分子内的弛豫时间是不同的。因此存在一个时间步长能很好地描述一种现象、但却无法描述同一系统另一种现象的可能性。所以时间步长应根据研究对象的特征以及模拟过程中能量的变化情况来选取，不存在一个通用的时间步长。

简化单位
分子动力学方法计算体系中原子或分子的运动，若采用cgs制，原子质量以g为单位，则原子质量的量纲为(g)，位置以cm为单位则量纲为(cm)，时间步长以秒为单位则量纲约为(s)。这些量的值均很小，可能在模拟中导致计算的误差。
通常选择能量，长度及质量为转换简化单位的标准，则

分子动力计算的初始设定
均匀相的液体系统的初始位置通常定位其晶体的结构，或均匀分布的面心立方结构。若分子体积很大，则以分子间不产生重叠为原则，近似系统的平衡结构愈佳。若复杂系统，则由分子力学能量最小化的结果为模拟起点。
初速度有两种常用的设置方法:一种是取平均速度，方向余弦取
之间的随机数，相对简单;另一种是按靠近Maxwell分布来取样，这种方法相比之下，更容易达到热力学平衡，即

计算前检查粒子的速度分布以使各方向的总动量为0。在进行用分子动力学方法模拟平衡问题的时候，初始一段时间的模拟结果往往丢弃掉，初始条件一般不会影响最终结果，但会影响模拟的收敛速度。选择合适的初始条件可节省所需的模拟时间。

边界处理
由于在之前跳蛙法中分子的速度对位置产生影响使分子已不在原来的坐标上，所以要进行边界处理以保证每个分子都在边界内。这里的处理方式由于X方向和Y方向的边界情况的不同而不同。
对于X方向仍然考虑其为周期边界，即一旦有分子越出晶胞，则会有另一个粒子从相邻的晶胞中进入这个晶胞，就好像这个分子本身还在晶胞中一样。这样无论分子在X方向上的位置怎样变化都不会影响到流体密度和状态。
对于Y方向，因为其有固定壁面，不能对其应用周期边界。而是要删除所有晶胞外墙壁处可能的搜索伙伴，即相邻晶胞，所以这个边界处理稍显复杂。一种方法可以用一层假设固定的分子来模拟墙壁的作用，将分子连接到一个恒温器上，墙就可以把流体内增加的热量转移出去。另一种方法是“随机”墙，即当一个分子试图穿过墙它被反射会内部。墙的粗糙程度体现在将反射分子速度的方向随机化，并调整其速度大小以匹配墙的温度。这个方法很简单，但是如果使用不当，可能会妨碍运动方程的整体性，尤其是当墙处于高密度、高温度的区域时。这里采用“随机”墙方法来处理Y方向的边界，并显示出了不错的效果。

3.6 分子碰撞计算
对于分子数为N的模拟体系，原则上任何两个分子间都存在相互作用，那么计算时须进行次运算，非常耗时。为提高效率，在实际模拟过程中用势能截断的方法，这里选取LJ势。可以预先选定一个截断半径，截断半径一般取或(为分子的直径)，只计算以截断半径为球体以内的分子之间的作用力，而当分子之间的距离超过截断半径时，分子间的相互作用就可以忽略。注意，周期性边界条件时，应小于即晶胞直径的一半。此时需要考虑的是给定的粒子与任一其它粒子的最近的周期映像相互作用。
为了计算特定分子都与哪些分子有力的作用，需要构造一个如图3.3的邻接列表。

图3.3 分子的邻接列表示意图
在划分晶胞的基础上，在可能发生相互作用的粒子对之间建立一个列表，命名为邻接列表。考虑到截断半径内的分子可能会运动到球的外部，设想一个具有厚度为的球壳，在时间内，原来位于半径为的球内的粒子都仍然在半径为=+的球内。计算1粒子受力时，只计算2、3、4三个粒子对1粒子的作用力之和，其它粒子的作用忽略不计。模拟过程中计算以内的粒子与1粒子间距的变化。比如，当5、6粒子与1粒子距离小于或7与1粒子距离小于时更新列表。
对邻接列表的处理实质上是要对粒子间相互作用进行描述。粒子碰撞后加速度也随即改变。这里要做如下假设：
(1)分子是球形和准刚性的，在碰撞前后分子的形状不发生改变。
(2)分子之间的碰撞是瞬时的，并且是完全弹性碰撞。
(3)碰撞时忽略所受的其它外力。
(4)碰撞的接触为点接触，不考虑分子的旋转运动。
(5)每个分子在一个时间步之内最多只能有一次碰撞。
(6)两分子碰撞过程中，不与其它分子发生碰撞。
在建立邻接列表之后，依据邻接列表所表述的可能有分子作用力的分子对，对它们的加速度依此进行改变。通过势能函数（L-J势能）

与位置积分

以及受力，由可推导出

从而可得加速度是当前时刻每个在其截断半径内的分子对其作用的矢量叠加。