分子动力学方法

），选择某种MD算法求解具有初值问题的运动方程，便得到相空间轨迹{r(N)(t)}，{p(N)(t)}对轨迹平均的宏观物理量A的表示为(3.7)如果宏观物理量为动能，它的平均为(3.8)由于在模拟过程中计算出的动能值是在不连续的路径上的值，因此公式（3.8）可以表示为在时间的各个间断点μ上计算动能的平均值（3.9）在MD模拟过程中，温度是需要加以监测的物理量，特别是在模拟的起始阶段。根据能量均分定理，我们可以从平均动能值计算得到温度值：（3.10）其中d为每个粒子的自由度，如果不考虑系统所受的约束，则d=3。系统内部的位形能量的轨道平均值为（3.11）假定位势在处被截断，那么上式计算出的势能以及由此得到的总能量就包含有误差。为了对此偏差rC作出修正，我们采用对关联函数来表示位能（3.12）式中的g(r)就是对关联函数，它是描述与时间无关的粒子间关联性的量度。g(r)的物理意义是当原点r=0处有一个粒子时，在空间位置r的点周围的体积元中单位体积内发现另一个粒子的几率。若n(r)为距离原点r到r+Δr之间的平均粒子数，则（3.13）在MD模拟过程中，所有的距离已经在力的计算中得到，因而很容易计算对关联函数的值。图3.2为由计算机模拟得到的两组不同参数下的对关联函数的例子。由于位势的截断，对关联函数仅对＜以rCL/2下的距离有意义。在公式（3.11）中，所有的位能都加到截断距离为止，尾部修正可以取为（3.14）压强可以通过计算在面积元dA的法线方向上净动量转移的时间平均值来得到，也可以利用含对关联函数的维里状态方程计算。该维里状态方程可以写为5（3.15）至于势能的计算，我们可以把积分划分为两项，一项是由相互作用力程之内的贡献引起的，一项是对位势截断的改正项：（3.16）其中长程改正项为：（3.17）下面我们将讨论具体如何进行MD模拟：图6.2由计算机模拟得到的两组不同参数下的对关联函数4、平衡态分子动力学模拟在经典MD模拟方法的应用当中，存在着对两种系统状态的MD模拟。一种是对平衡态的MD模拟，另一类是对非平衡态的MD模拟。对平衡态系综MD模拟又可以分为如下类型：微正则系综的MD（NVE）模拟，正则系综的MD（NVT）模拟，等温等压系综MD（NPT）模拟和等焓等压系综MD（NPH）模拟等。下面我们仅对平衡态的MD方法中前两类模拟做简单的介绍。1）、微正则系综的MD模拟在进行对微正则系综的MD模拟时，首先我们要确定所采用的相互作用模型。我们假定一个孤立的多粒子体系，其粒子间的相互作用位势是球对称的，则其哈密顿量可以写为（4.1）其中rij是第i个粒子与第j个粒子之间的距离。在这个微正则系综中，由于这个系统的哈密顿量中不显式地出现时间关联，因而系统的能量是个守恒量。系统的体积和粒子数也是不变的。此外；由于整个系统并未运动，所以整个系统的总动量P恒等于零。这就是系统受到的四个约束。由该系统的哈密顿量可以推导出牛顿方程形式的运动方程组(4.2)要用数值求解的方法解出（4.2）微分方程组，类似于本章第二节中介绍的Veriet方法方程组（4.2）的求解变成求解方程组：(4.3)6该方程组反映出：从前面t和t-h时刻这两步的空间坐标位置及t时刻的作用力，就可以算出下一步t+h时刻的坐标位置。下面为了将（4.3）式写成更简洁的形式，令(4.4)则从（4.3）式可以得到如下差分方程组的形式(4.5)如果已知一组初始空间位置{{r(0)},{r(1)}，则通过求解方程组（4.5）一步步得到{r(2)},{r(3)},……。iiii由空间坐标又可以算出粒子的运动速度为(4.6)这里在第n＋l步算出的速度是前一时刻；即第n步的速度、因而动能的计算比势能的计算落后一步。根据上述原理我们可以将粒子数恒定、体积恒定、能量恒定的微正则系综（NVE）的MD模拟步骤设计如下：（1）给定初始空间位置{{r(0)},{r(1)}，(i=1,2,3,……N)ii（2）计算在第n步时粒子所受的力{F(n)}：F(n)F(t)。iiin（3）利用公式r(n1)2r(n)r(n1)F(n)h2/m计算在第n+l步时所有粒子所处的空间位置{r(n1)}iiiii（4）计算第n步的速度：v(n)(r(n1)r(n1))/2hiii（5）返回到步骤（2），开始下一步的模拟计算。如前所述，用上述形式的Verlet算法，动能的计算比势能的计算落后一步。此外，这种算法不是自启动的。要真正求出微分方程组（4.2）的解，除了需要给出初始空间位置{r(0)}外，还要求另外给出一i组空间位置{r(1)}。实际，有时候采用改进后的计算方法可能更方便；即把N个粒于的初始位置放在网i格的格点上，然后加以扰动。如果初始条件是空间位置和速度，则采用下面的公式来计算空间位置{r(1)}i（4.7）然后再按上述模拟步骤进行计算。verlet算法的速度变型形式将会使其数值计算的稳定性得到加强。下面我们就此做简单介绍。令（4.8）则公式（4.5）写为（4.9）上式在数学上与（4.5）式是等价的，并称为相加形式。由此Verlet算法的速度形式的模拟步骤可以表述为（l）给定初始空间位置{r(1)}，(i=1,2,3,……N)i7(2）给定初始速度{v(1)}i（3）利用公式：r(n1)r(n)hv(n)F(n)h2/2m，计算在第n＋1步时所有粒子所处的空间位置iiii{r(n1)}。i（4）计算在第n+1步时所有粒子的速度{v(n1)}：v(n1)v(n)h(F(n1)F(n))/2miiiii（5）返回到步骤（3），开始第n+2步的模拟计算。Verlet速度形式的算法比前一种算法好些。它不仅可以在计算中得到同一时间步上的空间位置和速度，井且数值计算的稳定性也提高了。一般情况下，对于给定能量的系统不可能给出精确的初始条件。这时需要先给出一个合理的初始条件，然后在模拟过程中逐渐调节系统能量达到给定值。其步骤为：首先将运动方程组解出若干步的结果；然后计算出动能和位能；假如总能量不等于给定恒定值，则通过对速度的调整来实现能量守恒。也就是将速度乘以一个标度（scaling）回子，该因子一般取为（4.10）然后再回到第一步，对下一时刻的运动方程求解。反复进行上面的过程，直到系统达到平衡。这样的模拟过程也称为平衡化阶段。采用对速度标度的办法，可以使速度发生很大变化。为了消除可能带来的效应，必须要有足够的时间让系统再次建立平衡。在到达趋衡阶段以后，必须检验粒子的速度分布是否符合麦克斯韦-波尔兹曼分布。2）、正则系综的MD模拟在统计物理中的正则系综模拟是针对一个粒子数N、体积V、温度T和总动量（pp0）ii为守恒量的系综（NVT）。这种情况就如同一个系统置于热浴之中，此时系统的能量可能有涨落，但系统温度则已经保持恒定。在正则系综的MD模拟中施加的约束与微正则系综中的不一样。正则系综MD方法是在运动方程组上加上动能恒定（即温度恒定）的约束，而不是像做正则系综的MD模拟中对运动方程加上能量恒定的约束。在正则系综MD的平衡化过程中，速度标度因子一般选下面的形式较为合适（4.11）我们可将正则系综MD的Verlet算法的速度形式的模拟具体步骤列在下面：（l）给定初始空间位置{t(1)}:(i1,2,,N)i（2）给定初始速度{v(1)}i（3）利用公式r(n1)r(n)hv(n)F(n)h2/2m计算在第n+1步时所有粒子所处的空间位置iiii{r(n1)}i8（4）计算在第n＋1步时所有粒于的速度{v(n1)v(n)h(F(n1)F(n))/2m}，动能和速度标度因iiii子（5）计算将速度{vn1}标度因于β值，并让该值作为下一次计算时，第n+1步粒子的速度i{vn1}{vn1}：ii（6）返回到步骤（3），开始第n＋2步的模拟计算。按照上面的步骤，对时间进行一步步的循环。待系统达到平衡后，则退出循环。这就是正则系综的MD模拟过程。下面我们举一个微正则系综的MD模拟的应用示例来看看模拟的结果。例：对一个总能量确定的单原于（氩）粒子系统的MD模拟计算。我们具体选取256个原子的模拟。粒于间的相互作用位势为Lennard-Jones势（4.12）其中-ε为位势的极小值（取ε为能量单位），其位置在r=21/6σ处、该体系的粒于限制在一个立方体的箱子，边界上采用最小象力约定。我们采用自然单位制。长度和时间的标度单位分别为σ和(mσ2/48ε)1/2（对氩原于该时间单位为310-12秒），这样就使得运动方程为无量纲形式。模拟时我们考虑两个相图上的点(T*，ρ*)=(2.53，0.636)，（0.722，0.83134），它们分别具有两种立方体的尺寸，即L＝7.83和L＝6.75。初始条件假定为：各个原子处于一个面心立方格子的格点上，而速度按相应温度下的波尔兹曼分布抽样取值。位势的截断取两个值和，用以比较其rc=2.5rc=3.6对模拟结果的影响。在执行平衡化过程中，调节粒子速度的标度因于为(4.13)反复上面的速度调节，直到系统能量达到给定值、在这个例子中，平衡化过程用了1000步MD模拟。模拟结果列于表4.1中。表中的误差为误差。系统总动能的模拟演化过程由图4.1给出。实际上，图中显示出在数百步后动能就达到平衡了。图4.2则显示出位能的平衡化过程。系统总能量的平衡化过程则由图4.3表示，其平衡化是通过对粒子速度的调节跳跃式地达到的。图4.4为动能的分布图，模拟得到的平均速度为v＝0.3654，而理论上该值应当是v1.13T*/240.3668。这个结果已经是相当不错了，因为我们只对256个粒子的系统进行了模拟。而且速度大于平均速度的粒子数所占百分比与期望值46.7％也一致。表中的数据表明模拟结果与所选择的截断距离值变化并不灵敏。9表4.1对256个粒子的氩原子系统进行1000步微正则系综MD模拟的结果趋衡到T*=2.53，ρ*=0.636趋衡到T*=0.722，ρ*=0.83134图4.1动能演化过程图（T*＝2.53）图4.2位能演化过程图（T*＝2.53）图4.1总能量演化过程图（T*＝2.53）图4.2动能的分布图（T*＝2.53）10