分子动力学方法.docx

《分子动力学方法.docx》由会员分享，可在线阅读，更多相关《分子动力学方法.docx（10页珍藏版）》请在课桌文档上搜索。

1、分子动力学方法一、引言计算机模拟中的另一类确定性模拟方法，即统计物理中的所谓合于动力学方法(MOleCUIarDynamicsMethod)o这种方法是按该体系内部的内禀动力学规律来计算并确定位形的转变。它首先需要建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性，从而得到系统的宏观性质。在这样的处理过程中我们可以看出：MD方法中不存在任何随机因素。在MD方法处理过程中方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中,每个分子都各自服从经典的牛

2、顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现Boltzman的统计力学途径。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，讨算量大，占内存也多、本节将介绍分子动力学方法及其应用。原则上，MD方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其他的微观粒子。实际上，MD模拟方法和随机模拟方法一样都面临着两个基本限制：一个是有限观测时间的限制；另一个是

3、有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下(即粒子数日趋于无穷时)的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响体系的某些性质。对于MD方法，向然的系综是微正则系综，这时能量是运动常量。然而，当我们想要研究温度和(或)压力是运动常量的系统时，系统不再是封闭的。例如当温度为常量的系统可以认为系统是放置在一个热俗中。当然，在MD方法中我们只是在想像中将系统放入热浴中。实际上，在模拟计算中具体所采取的做法是对一些自由度加以约束。例如在恒温体系的情况下，

4、体系的平均动能是一个不变量。这时我们可以设计一个算法，使平均动能被约束在一个给定值上。由于这个约束，我们并不是在真正处理一个正则系综，而实际上仅仅是复制了这个系综的位形部分。只要这一约束不破坏从一个状态到另一个状态的马尔科夫特性，这种做法就是正确的。不过其动力学性质可能会受到这一约束的影响。自五十年代中期开始，MD方法得到了广泛的应用。它与蒙特卡洛方法一起已经成为计算机模拟的重要方法。应用MD方法取得了许多重要成果，例如气体或液体的状态方程、相变问题、吸附问题等，以及非平衡过程的研究。其应用已从化学反应、生物学的蛋白质到重离子碰撞等广泛的学科研究领域。二、分子运动方程的数值求解采用MD方法时，

5、必须对一组分于运动微分方程做数值求解。从计算数学的角度来看，这个求解是一个初值问题。实际上计算数学为了求解这种问题己经发展了许多的算法，但并不是所有的这些算法都可以用来解决物理问题。下面我们先以一个一维谐振子为例，来看一下如何用计算机数值计算方法求解初值问题。一维谐振子的经典哈密顿量为/1这里的哈密顿量(即能量)为守恒量。假定初始条件为x(p)、P(O),则它的哈密顿方程是对时间的一阶微分方程更=蛆，电=-”=七dIaPmdzax(2.2)现在我们要用数值积分方法计算在相空间中的运动轨迹(X(t)、p(t)o我们采用有限差分法，将微分方程变为有限差分方程，以便在计算机上做数值求解，并得到空间坐

6、标和动量随时间的演化关系。首先，我们取差分计算的时间步长为h,采用我们有限差分法中的一阶微分形式的向前差商表示，即直接运用展开到h的一阶泰勒展开公式/(r+)=(r)+O(2)即y()-y(r)山h(2.3)则微分方程(2.2)可以被改写为差分形dxx(r+h)-x(r)M)dr=、m(2.4)电“凶+A)-K)=*r)也h(2.5)将上面两个公式整理后，我们得到解微分方程(2.2)的欧拉(EUler)算法x(+)-x(r)+步区m(2.6)p(r+h)=pQ)-M)(27)这是x(t)、p(t)的一组递推公式、有了初始条件x(0)、P(O),就可以一步一步地使用前一时刻的坐标、动量值确定下一

7、时刻的坐标、动量值。这个方法是一步法的典型例子。由于在实际数值计算时h的大小是有限的，因而在上述算法中微分被离散化为差分形式来计算时总是有误差的。可以证明一步法的局部离散化误差与总体误差是相等的，都为O(h2)的量级。在实际应用中，适当地选择h的大小是十分重要的。h取得太大，得到的结果偏离也大，甚至于连能量都不守恒：h取得太小，有可能结果仍然不够好。这就要求我们改进计算方法，进一步考虑二步法。实际上泰勒,开户的一般形式/(r)w(/)+o(*)E(2.8)其中O(hn+)表示误差的数量级。前面叙述的欧拉算法就是取n=l。现在考虑公式(2.8)中直到含h的二次项的展开(即取n=2),则得到/(，

8、+M=刖+g+(空df2d?(2.9)诉)(2.10)/(f-)f(t)fg+驾+dr2dx2将上面两式相加、减得到含二阶和一阶导数的公式%=+)-2(r)+/(r-)J(2.11)苣=f4+A)-f(，一)也2h(2.12)令f(l)=x(t),利用牛顿第二定律公式产chF公式(ZU)写为坐标的递推公式Mr十九)=-Mr-力)+2x(f)+h2上m(2.13)公式(2.12)写为计算动量的公式得到p(f)=mi()=nv(f)-Z+)-x-)J2(2.14)这样我们就推导出了一个比(2.6)和(2.7)更精确的递推公式。这是二步法的一种，称为Verle方法。还有其他一些二步法，如龙格一库塔(

9、RUnge-KUtta)方法等。当然我们还可以建立更高阶的多步算法，然而大部分更高阶的方法所需要的内存比一步法和二步法所需要的大得多，并且有些更高阶的方法还需要用迭代来解出隐式给定的变量，内存的需求量就更大。并且当今的计算机都仅仅只有有限的内存，因而并不是所有的高阶算法都适用于物理系统的计算机计算。在实际数值计算中，我们必须特别注意舍入误差和稳定性问题。为了减少舍入误差，我们可以采用高精度计算,并且要避免相近大小的数相消，以及数量级相差很大的两个数相加和注意运算顺序。三、分子动力学模拟的基本步骤在计算机上对分子系统的MD模拟的实际步骤可以划分为四步；首先是设定模拟所采用的模型：第二，给定初始条

10、件；第三，趋于平衡的计算过程；最后是宏观物理量的计算。下面就这四个步骤分别做简单介绍。1、模拟模型的设定设定模型是分子动力学模拟的第一步工作。例如在一个分子系统中，假定两个分子间的相互作用势为硬球势，其势函数表示为K)Jf如果，。，V(r)4e(其中是行势的最小宿(实除上,1更常用的她的49.l)Lennard-JoneS型势。它的势函数表示为(3.1)可以确定能量的单位)，这个最小值出现在距离I等于21版的地方(。可以确定为长度的单位)模型确定后，根据经典物理学的规律我们就可以知道在系综模拟中的守恒量。例如对在微正则系综的模拟中能量、动量和角动量均为守恒量。在此系综中它们分别表示为(3.2)

11、P=EPIi（33）m=%xPi（3.4）其中p=mro由于我们只限于研究大块物质在给定密度卜的性质，所以必须引进一个叫做分子动力学元胞的体积元：以维持一个恒定的密度。对气体和液体，如果所占体积足够大，并且系统处于热平衡状态的情况下，那么这个体积的形状是无关紧要的。对于晶态的系统，元胞的形状是有影响的。为了计算简便，对于气体和液体，我们取一个立方形的体积为MD元胞。设MD元胞的线度大小为L,则其体积为L3。由于引进这样的立方体箱子，将产生六个我们不希望出现的表面。模拟中碰撞这些箱的表面的粒子应当被反射回到元胞内部，特别是对粒子数目很少的系统。然而这些表面的存在对系统的任何一种性质都会有重大的影

12、响。为了减小引入的表面效应，我们采用周期性边界条件。采用这种边界条件，我们就可以消除引入的表面效应，构造出一个准无穷大的体积来更精确地代表宏观系统。实际上，这里我们做了一个假定，即让这个小体积元胞镶嵌在一个无穷大的大块物质之中。周期性边界条件的数学表示形式为其中A为任意的可观测曾3、n；n3对任意整数）聂小龙算条件就是命令基本MD元胞完全等同地重复无穷多次。具体在实现该边界条件时是这样操作的：当有一个粒子穿过基本MD元胞的六方体表面时，就让这个粒子以相同的速度穿过此表面对面的表面重新进入该MD元胞内。在分子动力学模拟中考虑粒子间的相互作用时，通常采用最小像力约定。这个约定是在由无穷重复的MD基

13、本元胞中，一个粒子只同它所在的基本元胞内的另外N-I个（设在此元胞内有N个粒子）中的每个粒子或其最邻近影像粒子发生相互作用。如果r处的粒子i同r处的粒子j之间的距离为一4）F（对一切的Q36）实际工感不药定就磷过满足不等式条件rck9=(3.9)在MD模拟过程中，温度是需要加以监测的物理量，特别是在模拟的起始阶段。根据能量均分定理，我们可以从平均动能值计算得到温度值：Z(3.10)其中d为每个粒子的自由度，如果不考虑系统所受的约束，则d=3。系统内部的位形能量的轨道平均值为4r)M(3.11)假定位势在处被截断，那么上式计算出的势能以及由此得到的总能量就包含有误差。为了对此偏差作出修正，我们采

14、用对关联函数来表示位能UIN=2pfu(r)g(r)r2drJ(3.12)式中的g(r)就是对关联函数，它是描述与时间无关的粒子间关联性的量度。g(r)的物理意义是当原点r=0处有一个粒子时，在空间位置r的点周围的体积元中单位体积内发现另一个粒子的几率。若n(r)为距离原点r到r+Ar之间的平均粒子数，则g=77而后(313)在MD模拟过程中，所有的距离已经在力的计算中得到，因而很容易计算对关联函数的值。图3.2为由计算机模拟得到的两组不同参数下的对关联函数的例子。由于位势的截断，对关联函数仅对hVL/2以下的距离有意义。在公式(3.11)中，所有的位能都加到截断距离为止，尾部修正可以取为Uc

15、=2pu(r)g(r)r2dr(3.14)压强可以通过计算在面积元dA的法线方向上净动量转移的时间平均值来得到，也可以利用含对关联函数的维里状态方程计算。该维里状态方程可以写为(3.15)P=伙BT至于势能的计算，我们可以把积分划分为两项，一项是由相互作用力程之内的贡献引起的，一项是对位势截断的改正项：P=MT-公病-PC(3.16)其中长程改正项为:户C=詹。(噂向7)下面我们将讨论具体如何进行MD模拟：图6.2由计算机模拟得到的两组不同参数下的对关联函数4、平衡态分子动力学模拟在经典MD模拟方法的应用当中，存在着对两种系统状态的MD模拟。一种是对平衡态的MD模拟，另一类是对非平衡态的MD模

16、拟。对平衡态系综MD模拟又可以分为如下类型：微正则系综的MD(NVE)模拟，正则系综的MD(NVT)模拟，等温等压系综MD(NFT)模拟和等熔等压系综MD(NPH)模拟等。下面我们仅对平衡态的MD方法中前两类模拟做简单的介绍。1)、微正则系综的MD模拟在进行对微正则系综的MD模拟时，首先我们要确定所采用的相互作用模型。我们假定一个孤立的多粒子体系，其粒子间的相互作用位势是球对称的，则其哈密顿量可以写为H=挈+歙V)(4.(1)其中丐是第i个粒子与第j个粒子之间的距离。在这个微正则系综中，由于这个系统的哈密顿量中不显式地出现时间关联，因而系统的能量是个守恒量。系统的体积和粒子数也是不变的。此外；

17、由于整个系统并未运动，所以整个系统的总动量P恒等于零。这就是系统受到的四个约束。由该系统的哈密顿量可以推导出牛顿方程形式的运动方程组=G=L2N)也my(4.(2)要用数值求解的方法解出(4.2)微分方程组，类似于本章第二节中介绍的Veriet方法方程组(4.2)的求解变成求解方程组：ri(t+Zo=2rf(O-riG-)+Fi(t)2fm,0“,N)该方程组反映出：从前面I和l-h时刻这两步的空间坐标位置及I时刻的作用力，就可以算出下一步t+h时刻的坐标位置。下面为了将(4.3)式写成更简洁的形式，令t.=nh，r严=也)，Fp)=Fi(J)(4.4)则从(4.3)式可以得到如下差分方程组的

18、形式产D=才)-rFF+FtWm,/=U.)f4s.如果己知一组初始空间位置,，则通过求解方程组(4.5)一步步得到L2),z,。fiii由空间坐标又可以算出粒子的运动速度为S)=卜产DI)/2(4.6)这里在第n+1步算出的速度是前一时刻；即第n步的速度、因而动能的计算比势能的计算落后一步。根据上述原理我们可以将粒子数恒定、体积恒定、能量恒定的微正则系综(NVE)的MD模拟步骤设计如下：(1)给定初始空间位置ro),r),(i=l,2,3,N)ii(2)计算在第n步时粒子所受的力凡”)：F(H)=F(Z)oiiin利用公式八用)=2r(n)-r(n-i)+F(n)h2m计算在第n+1步时所有

19、粒子所处的空间位置(用)iiiii(4)计算第n步的速度：V(n)=(r(n)-(n-I)2iii(5)返回到步骤(2),开始下一步的模拟计算。如前所述，用上述形式的Verlet算法，动能的计算比势能的计算落后一步。此外，这种算法不是自启动的。要真正求出微分方程组(4.2)的解，除了需要给出初始空间位置厂(。)外，还要求另外给出一i组空间位置上。实际，有时候采用改进后的计算方法可能更方便；即把N个粒于的初始位置放在网I格的格点上，然后加以扰动。如果初始条件是空间位置和速度，则采用下面的公式来计算空间位置Ki力=6)+肋程+坪0%2/6a7)然后再按上述模拟步骤进行计算。Verlel算法的速度变

20、型形式将会使其数值计算的稳定性得到加强。下面我们就此做简单介绍。令则公式(4.5)写为rn)=r-0+HT)ZW=ZfT)+mlF(4.9)上式在数学上与(4.5)式是等价的，并称为相加形式。由此VerIet算法的速度形式的模拟步骤可以表述为(1)给定初始空间位置-,(i=l,2,3,N)(2)给定初始速度W(3)利用公式：r(/Hi)=r()+/?v(n)+F()/j2/2/w,计算在第n+1步时所有粒子所处的空间位置r(i)o(4)计算在第n+1步时所有粒子的速度也叫：v(+i)=v(n)+(F(/+i)+F2,p*=。.沿134UE2.5279.139.57-1421.9820.15-1142.923.6275.11土9.72-1496.4521.61-1221.38OTV初以上百分比2.50.729710.0250.19650.00347.083.60.71920.0250.19490.00346.42图4.1动能淌犯程糊（%*1网）600MD步数图4.2位解湖叱过稗糊（TW.53MD步数T*=2.53P0.6362S5025据W齿用壬*9图4总颍演魂场对呼=2哪MD步数0EJkW