Verlet算法

Verlet算法从计算数学的观点来看，MD方法是一个初值问题，主要用有限差分方法来求解。能够求解Newton运动方程的有限差分算法很多，但分子动力学计算原子数目较多，所以对时间和空间要求都比较高，需要研究适合于分子动力学的时间积分算法。时间积分算法基于有限差分法，对时间进行离散，t为时间步长，由t时刻的变量及其时间导数推演得到t+t时刻的物理量，逐步计算即得到整个历史过程的物理量。这个过程是近似的，不可避免地引入误差，包括截断误差和舍入误差。前者主要是有限差分算法采用Taylor展开截断某些高阶项引入的，后者是计算机本身数值精度引入的。两种误差都能通过减小时间步长t而减小。对于比较大的t，截断误差起主导作用；但随着t下降截断误差迅速降低，舍入误差随t下降变化不大，并在较小的t时起主导作用。分子动力学时间步长一般随模拟的材料、温度、物理过程而变化，一般为fs(1510秒)量级。进行数值求解的算法有多种，其共同的出发点是，将粒子的位置和其它力学量如速度、加速度等展开为Taylor级数。最常用的数值求解算法有Verlet算法，leap-frog算法，Beeman算法及Gear所提出的校正预测法。1.Verlet算法Verlet提出的Verlet算法在分子动力学中运用最为广泛，也是最简单的。它运用原子在t时刻的位置)(tr和加速度)(ta及tt时刻的位置，计算出tt时刻的位置。将粒子的位置以Taylor公式展开，即...))((!21)()()(222ttrdtdttrdtdtrttr将式中的t换为-t，得...))((!21)()()(222ttrdtdttrdtdtrttr两式相减得速度式)()(21)(ttrttrtdtdrtv两式相加得222))(()(2)()(ttrdtdtrttrttr因)()(22tatrdtd，故依据上式可由t及tt的位置预测tt时的位置。Verlet算法执行简明，需要的内存小，但其缺点在于速度式中含有t1项，而实际计算中通常选取很小的t值，且位置)(ttr要通过小项与非常大的两项)(2tr与)(ttr的差的相加得到，容易造成精度损失。另外，它不是一个自启动算法，新位置必须由t时刻与前一时刻tt的位置得到。在t=O时刻，只有一组位置，所以必须通过其他方法得到tt的位置。获得tt时刻的位置的方法之一时应用近似式)0(*)0()(vtrtr。2.跳蛙法（leap-frogmethod）鉴于Verlet算法的一些缺点，Hockney提出了跳蛙法。跳蛙法是从Verlet法推导出来的。它在半个积分时间步得到速度，并利用这一速度计算新的位置。位置和速度表达式为)21()()()()21()21(ttvtrttrttattvttviiiiii计算时假设已知)21(ttvi与)(tri，由t时的位置)(tri计算质点所受的力与加速度)(tai，再预测时间为tt21时的速度)21(ttvi，以此类推。根据)21(ttvi与)21(ttvi，可得时间为t时的速度为)21()21(21)(ttvttvtviii利用跳蛙法计算仅需储存)21(ttvi与)(tri两类信息，既节约储存空间，而且准确性及稳定性较高。跳蛙法不需要计算下一步位置就可以得出速度，但需要注意的是速度并未与位置在同一时间定义，结果是动能和势能也未同时定义。所以不能直接计算总能量。3.Beeman法Beeman法是除Verlet跳蛙法外，另一种较为常见的方法。其积分公式如下tttatattatvttvtttatattvtrttriiiiiiiiii)]()(5)(2[61)()()]()(4[61)()()(2此方法需储存)(tri、)(tvi与)(ttai，储存量大于Verlet的跳蛙法。但优点在于可以引用较长的积分间隔t。Beeman方法所引用的积分步长t可为Verlet方法的4~3倍，而具有相同的准确性。4.校正预测法（predictor-correctormethod）Gaer提出了基于预测-校正积分方法的校正预测法。因为经典运动粒子的轨迹为连续的，故于时间tt时的位置、速度等可由时间t的Taylor展开式预测得...)()(...)()()(...)(21)()()(...)(61)(21)()()(232tbttbttbtattattbttatvttvttbttattvtrttrpppp由于这些物理量来自Taylor展开式，并非由解Newton运动方程式而来，所以式（2.17）中所产生的速度、加速度等并非完全正确。为了解决这个问题，用所预测的位置)(ttrp计算所受的力及正确的加速度)(ttap。设正确的加速度与预测的加速度之间的误差为)()()(ttattcttapc得各量之校正式为)()()()()()()()()()()()(3210ttacttbttbttacttattattacttvttvttacttrttrpcpcpcpc式中，0c、1c、2c、3c均为常数。这仅为Gear的一次预测校正法，也可以推展至更高次的校正法。对于大多数MD应用来说，Verlet型的算法就可以满足基本的要求了。但有时采用高阶算法更方便。一般来说Gear算法中的Taylor展开阶数越高，精度越好。但占用内存比较大，这对于大量分子系统的模拟不如Verlet算法方便。当t减小时，Gear的高次预测校正法误差较小，但t增大时，Verlet方法较佳。边界条件执行分子动力计算通常选取一定数目N的分子，将其置于一立方体的盒子中。设盒子的边长为L，则其体积为V=3L。若分子的质量为m，则系统的密度为：3LNm计算系统的密度应等于实验所测定的密度。介质的宏观性质是由大量粒子所为，所以在模拟中，粒子的数量只有足够多，才能真实地再现介质的宏观性质。然而由于计算机的计算能力的限制，我们在模拟时粒子数一般应少于10000个。同时粒子数量越少，表面效应越明显，因为此时大量的粒子是处于表面的。这样少的粒子组成的集团和真实体系相比，其表面原子与体内原子数之比显然过大，这必然造成表面效应。无论模拟的区域是否在容器中，表面的粒子受力仍不同于实际情况，从而很难正确的描述其宏观性质。为了减小有限尺寸的影响，在分子动力学模拟中采用了各种有效的边界条件。通常用于分子动力学模拟中的边界条件有两种:周期性边界条件、非周期性边界条件。（1）模拟计算中，为使计算中系统的密度维持恒定，通常采用周期性边界条件[2]，以二维的计算系统为例，图2.3显示了二维盒中系统粒子的排列及移动方向。图2.3二维周期性系统的例子排列与移动图2.4粒子的最近映像图中位于中央的盒子表示所模拟的系统，其周围盒子与模拟系统具有相同的排列及运动。当计算系统中任一粒子移出盒外，则必有一粒子由相对的方向移入。如图2.3中的第2个粒子。这样的限制条件使得系统中的粒子数维持恒定，密度不变，符合实际的要求。计算系统中分子间的作用力时，采取最近镜像方法]2[。如图2.4(粒子的最近映像)所示，计算分子2与5的作用力，是取与分子2和其最近的距离镜像分子5。因为在计算中利用最近镜像的观念，因此需采用截断半径cr的方法计算非键结的远程作用力，否则会因重复计算粒子所受的力而导致不正确的结果。图2.5为一般的vanderWaals势能图。图2.5vanderWaals势能图图中显示，当rcr时，势能的值已非常趋近于零，表示分子间的vanderWaals作用力已可以忽略不计。在很多分子动力学模拟过程中，若分子间的距离大于截断半径，则将其作用视为零。截断半径最大不能超过盒长的一半，即2Lrc。一般的原子所选取的截断半径约为A15~10。（2）分子模拟并不都是周期边界，在许多情况下，如对于液滴体系，溶液中沉淀的分子团簇，并不需要周期边界。如图2.7所示，图中的外壳，即阴影部分是体系的边界区。边界区内的分子可以设计为位置固定的。体系中部的分子是可以移动的。图2.7非周期边界情形时间步长的选取时间步长与采用的数值积分方法、研究的对象有关。在选择了数值算法后，时间步长的选取是影响计算结果的一个关键因素。这里存在一对矛盾，一方面应尽量选择大的时间步以模拟更长的时间，从而在相空间抽样的比例更大；另一方面，过大的时间步将导致结果精度降低、甚至分子动力学的崩溃。时间步长与所研究现象的弛豫时间密切相关。一般地，时间步长可安全地取为弛豫时间的51~101。对于金属单原子晶体结构，时间步长可大致取为10sf。但模拟系统可能存在几个不同的时间尺度，如对于分子系统，分子间和分子内的弛豫时间是不同的。因此存在一个时间步长能很好地描述一种现象、但却无法描述同一系统另一种现象的可能性。所以时间步长应根据研究对象的特征以及模拟过程中能量的变化情况来选取，不存在一个通用的时间步长。简化单位分子动力学方法计算体系中原子或分子的运动，若采用cgs制，原子质量以g为单位，则原子质量的量纲为2310(g)，位置以cm为单位则量纲为810(cm)，时间步长以秒为单位则量纲约为161310~10(s)。这些量的值均很小，可能在模拟中导致计算的误差。通常选择能量，长度及质量m为转换简化单位的标准，则TkTb*EE*tmt212*FF*LL*maa/*21*/muu分子动力计算的初始设定均匀相的液体系统的初始位置通常定位其晶体的结构，或均匀分布的面心立方结构。若分子体积很大，则以分子间不产生重叠为原则，近似系统的平衡结构愈佳。若复杂系统，则由分子力学能量最小化的结果为模拟起点。初速度有两种常用的设置方法:一种是取平均速度mkT3，方向余弦取1~1之间的随机数，相对简单;另一种是按靠近Maxwell分布来取样，这种方法相比之下，更容易达到热力学平衡，即zyxTkmvkmTvpBB,,,21exp)2()(22/1计算前检查粒子的速度分布以使各方向的总动量为0。在进行用分子动力学方法模拟平衡问题的时候，初始一段时间的模拟结果往往丢弃掉，初始条件一般不会影响最终结果，但会影响模拟的收敛速度。选择合适的初始条件可节省所需的模拟时间。边界处理由于在之前跳蛙法中分子的速度对位置产生影响使分子已不在原来的坐标上，所以要进行边界处理以保证每个分子都在边界内。这里的处理方式由于X方向和Y方向的边界情况的不同而不同。对于X方向仍然考虑其为周期边界，即一旦有分子越出晶胞，则会有另一个粒子从相邻的晶胞中进入这个晶胞，就好像这个分子本身还在晶胞中一样。这样无论分子在X方向上的位置怎样变化都不会影响到流体密度和状态。对于Y方向，因为其有固定壁面，不能对其应用周期边界。而是要删除所有晶胞外墙壁处可能的搜索伙伴，即相邻晶胞，所以这个边界处理稍显复杂。一种方法可以用一层假设固定的分子来模拟墙壁的作用，将分子连接到一个恒温器上，墙就可以把流体内增加的热量转移出去。另一种方法是“随机”墙，即当一个分子试图穿过墙它被反射会内部。墙的粗糙程度体现在将反射分子速度的方向随机化，并调整其速度大小以匹配墙的温度。这个方法很简单，但是如果使用不当，可能会妨碍运动方程的整体性，尤其是当墙处于高密度、高温度的区域时。这里采用“随机”墙方法来处理Y方向的边界，并显示出了不错的效果。3.6分子碰撞计算对于分子数为N的模拟体系，原则上任何两个分子间都存在相互作用