分子动力学模拟方法.ppt

第四章分子动力学模拟方法,1957年：基于刚球势的分子動力学法（AlderandWainwright）1964年：利用Lennard-Jone势函数法对液态氩性质的模拟（Rahman）1971年：模拟具有分子团簇行为的水的性质（RahmanandStillinger）1977年：约束动力学方法（Rychaert,CiccottivanGunsteren）1980年：恒压条件下的动力学方法（Andersen法、Parrinello-Rahman法）1983年：非平衡态动力学方法（GillanandDixon）1984年:恒温条件下的动力学方法(Berendsenetal.)1984年：恒温条件下的动力学方法（Nos-Hoover法）1985年：第一原理分子動力学法（Car-Parrinello法）1991年：巨正则系综的分子动力学方法（CaginandPettit）,分子动力学简史,粒子的运动取决于经典力学(牛顿定律（F=ma）,课程讲解内容：经典分子动力学(ClassicalMolecularDynamics),分子动力学方法基础：,原理：计算一组分子的相空间轨道，其中每个分子各自服从牛顿运动定律：,初始条件：,分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法，可以预测纳米尺度上的材料动力学特性。通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与原子运动路径相关的基本过程。在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。分子动力学方法是确定性方法，一旦初始构型和速度确定了，分子随时间所产生的运动轨迹也就确定了。,分子动力学方法特征：,分子动力学的算法：有限差分方法,一、Verlet算法,粒子位置的Taylor展开式：,粒子位置：,粒子速度：,粒子加速度：,开始运动时需要r(t-t)：,+,缺点：Verlet算法处理速度非常笨拙,Verlet算法的表述：,算法启动,规定初始位置规定初始速度扰动初始位置：计算第n步的力计算第n+1步的位置：计算第n步的速度：重复至,Verlet算法程序：,Do100I=1,NRXNEWI=2.0*RX(I)RXOLD(I)+DTSQ*AX(I)RYNEWI=2.0*RY(I)RYOLD(I)+DTSQ*AY(I)RZNEWI=2.0*RZ(I)RZOLD(I)+DTSQ*AZ(I)VXI=(RXNEWIRXOLD(I)/DT2VYI=(RYNEWIRYOLD(I)/DT2VZI=(RZNEWIRZOLD(I)/DT2RXOLD(I)=RX(I)RYOLD(I)=RY(I)RZOLD(I)=RZ(I)RX(I)=RXNEWIRY(I)=RYNEWIRZ(I)=RZNEWI100CONTINUE,优点：1、精确，误差O(4)2、每次积分只计算一次力3、时间可逆缺点：1、速度有较大误差O(2)2、轨迹与速度无关，无法与热浴耦联,Verlet算法的优缺点：,二、蛙跳(Leap-frog)算法：半步算法,1.首先利用当前时刻的加速度，计算半个时间步长后的速度：,2.计算下一步长时刻的位置：,3.计算当前时刻的速度：,t-t/2,t,t+t/2,t+t,t+3t/2,t+2t,v,r,v,开始运动时需要v(-t/2)：,Leap-frog算法的表述：,算法启动,规定初始位置规定初始速度扰动初始速度：计算第n步的力计算第n+1/2步的速度：计算第n+1步的位置：计算第n步的速度：重复至,Leap-frog算法的优缺点：,优点：1、提高精确度2、轨迹与速度有关，可与热浴耦联缺点：1、速度近似2、比Verlet算子多花时间,三、VelocityVerlet算法：,等价于,优点：速度计算更加准确,VelocityVerlet算法的表述：,算法启动,规定初始位置规定初始速度计算第n+1步的位置：计算第n+1步的力计算第n+1步的速度：重复至,Verlet三种形式算法的比较：,Verlet,Leap-frog,VelocityVerlet,四、预测校正(Predictor-Corrector)格式算法：,预测(Predictor)阶段：其基本思想是Taylor展开，,根据新的原子位置rp，可以计算获得校正后的ac(t+t),定义预测误差:,利用此预测误差，对预测出的位置、速度、加速度等量进行校正：,校正(Corrector)阶段：,预测阶段运动方程的变换：,定义一组矢量:,校正阶段运动方程的变换：,的形式：,C0，C1，C2，C3的值以及,C0，,取决于运动方程的阶数。,一阶运动方程：,二阶运动方程之一：,二阶运动方程之二：,五、积分时间步长t的选择：,室温下，t1fs(femtosecond10-15s)，温度越高，t应该减小,太长的时间步长会造成分子间的激烈碰撞，体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力,微正则系综分子动力学(NVEMD),它是分子动力学方法的最基本系综,具有确定的粒子数N，能量E和体积V,算法：,规定初始位置和初始速度对运动方程积分若干步计算势能和动能若能量不等于所需要的值，对速度进行标度重复至，直到系统平衡,微正则系综(NVE)MD模拟算法的流程图：,给定每个分子的初始位置ri(0)和速度vi(0),计算每个分子的受力Fi和加速度ai,解运动方程并求出每个分子运动一个时间步长后到达的位置所具有的速度,统计系统的热力学性质及其它物理量,统计性质不变？,打印结果，结束,Yes,No,移动所有分子到新的位置并具有当前时刻的速度,微正则系综MD模拟程序F3讲解（LJ,NVE）：,无因次量：,MD模拟中几个热力学量的计算：,对于由N个单原子组成的系统：,动能和温度：,采用对比量：,对于LJ流体：,势能：,采用对比量：,内能：,内能由势能和动能组成：,采用对比量：,压力：,采用对比量：,练习：推导LJ流体分子间力的表达式(fx,fy,fz及其对比量)：,势能函数形式：,力：,采用对比量：,=x,y,z,LJ分子间的维里项：,=x,y,z,采用对比量的运动方程形式：(以蛙跳(Leap-frog)算法为例),采用对比量：,最终得到：,同理得到：,速度的标度(VelocityScaling)：,根据能量均分原理，可知：,标度因子:,对比量,速度标度:,或,微正则系综MD模拟程序F3讲解（LJ,NVE）：,初始化：,READ(*,(A)TITLE!运行作业题目READ(*,*)NSTEP!运行步数READ(*,*)IPRINT!打印步数READ(*,(A)CNFILE!位型文件READ(*,*)DENS!对比密度READ(*,*)RTEMP!对比温度READ(*,*)RCUT!对比截断半径READ(*,*)DT!对比时间步长,CALLREADCN(CNFILE),初始位型：,面心立方(face-centeredcubic,FCC)：,每面中心有一格点,体心立方(body-centeredcubic,BCC)：,简单立方(simplecubic,SC)：,XL,初始位型：面心立方(FCC)(程序F23),NC=(REAL(N)/4.0)*(1.0/3.0)XL=1.0/REAL(NC)Y=0.5*XLR(1)=(0,0,0)R(2)=(0,Y,Y)R(3)=(Y,0,Y)R(4)=(Y,Y,0),M=0DO10I=1,NCDO10J=1,NCDO10K=1,NCDO11IJ=1,4RX(IJ+M)=RX(IJ)+XL*(K-1)RY(IJ+M)=RY(IJ)+XL*(J-1)RZ(IJ+M)=RZ(IJ)+XL*(I-1)11CONTINUEM=M+410CONTINUE,DO100I=1,NRX(I)=RX(I)-0.5RY(I)=RY(I)-0.5RZ(I)=RZ(I)-0.5100CONTINUE,将模拟盒子的中心移到原点：,初始速度：,简单的选择：,Vrandom(-0.5,0.5),=x,y,z,标度因子:,速度标度:,FACTOR=SQRT(RTEMP)DO100I=1,NVX(I)=FACTOR*(RANF(DUMMY)-0.5)VY(I)=FACTOR*(RANF(DUMMY)-0.5)VZ(I)=FACTOR*(RANF(DUMMY)-0.5)CONTINUE,随机安排初始速度：,标度初始速度：,SUMKX=0.0SUMKY=0.0SUMKZ=0.0DO200I=1,NSUMKX=SUMKX+VX(I)*2SUMKY=SUMKY+VY(I)*2SUMKZ=SUMKZ+VZ(I)*2200CONTINUE,BEITAX=SQRT(RTEMP/SUMKX)BEITAY=SQRT(RTEMP/SUMKY)BEITAZ=SQRT(RTEMP/SUMKZ)DO300I=1,NVX(I)=VX(I)*BEITAXVY(I)=VY(I)*BEITAYVZ(I)=VZ(I)*BEITAZ300CONTINUE,标度因子:,SUMX=0.0SUMY=0.0SUMZ=0.0DO200I=1,NSUMX=SUMX+VX(I)SUMY=SUMY+VY(I)SUMZ=SUMZ+VZ(I)CONTINUE,SUMX=SUMX/REAL(N)SUMY=SUMY/REAL(N)SUMZ=SUMZ/REAL(N)DO300I=1,NVX(I)=VX(I)-SUMXVY(I)=VY(I)-SUMYVZ(I)=VZ(I)-SUMZ300CONTINUE,控制体系的总动量为零：,从Maxwell分布中抽样：,x,x,dx,高斯(Gauss)分布:,对于等几率随机试验(Bernoulli试验)，大量的试验结果满足高斯分布,麦克斯韦速度分布定律:,由于：,=x,y,z,单位体积的分子再每个分量上的速度分布实际上就是高斯分布。,从Maxwell分布中抽样：,高斯(Gauss)分布的随机数生成方法:,生成随机数：i，i=1,2,12,SUM=0.0DO10I=1,12SUM=SUM+RANF(DUMMY)10CONTINUER=(SUM-6.0)/4.0R2=R*RGAUSS=(A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R,高斯(Gauss)分布的随机数生成(程序F24),FACTOR=SQRT(RTEMP)DO100I=1,NVX(I)=FACTOR*GAUSS(DUMMY)VY(I)=FACTOR*GAUSS(DUMMY)VZ(I)=FACTOR*GAUSS(DUMMY)CONTINUE控制总动量为零：同前面一样处理。,从Maxwell分布中抽样分布中随机安排初始速度：,微正则系综MD模拟程序F3讲解(LJ,NVE)：,量纲变换：,SIGMA=(DENS/REAL(N)*(1.0/3.0)RCUT=RCUT*SIGMADTDT*SIGMADENS=DENS/(SIGMA*3),模拟盒子的边长为1,长程校正：,微正则系综MD模拟程序F3讲解(LJ,NVE)：,SR3=(SIGMA/RCUT)*3SR9=SR3*3SIGCUB=SIGMA*3VLRC=(8.0/9.0)*PI*DENS*SIGCUB*REAL(N):*(SR9-3.0*SR3)WLRC=(16.0/9.0)*PI*DENS*SIGCUB*REAL(N):*(2.0*SR9-3.0*SR3),算法：算法启动,微正则系综MD模拟程序F3讲解(LJ,NVE)：,CALLFORCE(-DT,SIGMA,RCUT,NEWV,NEWVC,NEWW)CALLMOVE(-DT)CALLFORCE(-DT,SIGMA,RCUT,V,VC,W)CALLFORCE(DT,SIGMA,RCUT,V,VC,W)CALLKINET(OLDK)CALLMOVE(DT)CALLFORCE(DT,SIGMA,RCUT,NEWV,NEWVC,NEWW)CALLKINET(NEWK),算法：差分格式：,SR2=SIGSQ/RIJSQVIJ=4.0*(SR12-SR6)WIJ=24.0*(2.0*SR12-SR6)VELIJ=WIJ*DT/RIJSQDVX=VELIJ*RXIJ.VXI=VXI+DVX.VX(J)=VX(J)DVXV=V+VIJW=W+WIJCALLMOVE(DT),DO1000I=1,NRX(I)=RX(I)+VX(I)*DTRY(I)=RY(I)+VY(I)*DTRZ(I)=RZ(I)+VZ(I)*DT1000CONTINUE,MOVE(DT):,速度的标定（只用于平衡阶段）,SUMK=0.0DO200I=1,NSUMK=SUMKX+VX(I)*2+VY(I)*2+VZ(I)*2200CONTINUE,BEITA=SQRT(3.0*RTEMP/SUMK)DO300I=1,NVX(I)=VX(I)*BEITAVY(I)=VY(I)*BEITAVZ(I)=VZ(I)*BEITA300CONTINUE,正则系综分子动力学(NVTMD),具有确定的粒子数N，温度T和体积V,速度的直接标度热浴方法(AndersenThermostat)约束方法(阻尼力方法)系统扩展方法(ExtendedSystemsMethod),问题的关键：温度的约束,Nose-Hoover方法,一、热浴方法(AndersenThermostat),引入一个与虚拟粒子碰撞的随机力,想象系统浸在热浴当中,系统和热浴间的相互作用强度由随机碰撞的频率决定碰撞的几率等于Nudt如果一个粒子经历碰撞，它的速度将从约束温度下的Maxwell分布中随机抽取,总能量和总动量均不守恒,二、约束方法,是等动能(Iso-Kinetics)分子动力学方法,系统的运动方程为：,引入阻尼系数以保证将温度约束在恒定值,根据高斯最小约束原理：,三、Nose-Hoover扩展方法,基本思想：设想原系统与一个耦合系统共同组成一个扩展系统，允许热流在原系统和耦合系统之间交换。,Q：等效质量S:扩展坐标变量：热力学阻尼系数L：扩展系统的自由度,Predictor-correctoralgorithmisstraightforwardVerletalgorithmisfeasible,buttrickytoimplement,积分方案：,updateofxdependsonpupdateofpdependsonx,Nos-Hoover方法正确地描述了NVT系综中的动量和位型，而等动能方法只正确地描述了后者。,扩展系统的哈密顿量Hamiltonian守恒：,1.复杂分子体系的势能函数形式：,PotentialEnergy=StretchingEnergy+BendingEnergy+TorsionEnergy+Non-BondedInteractionEnergy这些方程与描述原子或键各种不同行为的参数就构成了力场，force-field.,UFF,OPLS,Amber,CVFF,Compass,分子模拟方法补充介绍：,BondStretchingEnergy,kbisthespringconstantofthebond.r0isthebondlengthatequilibrium.,Uniquekbandr0assignedforeachbondpair,i.e.C-C,O-H,BendingEnergy,kisthespringconstantofthebend.0isthebondlengthatequilibrium.,Uniqueparametersforanglebendingareassignedtoeachbondedtripletofatomsbasedontheirtypes(e.g.C-C-C,C-O-C,C-C-H,etc.),The“Hookeian”potential,kbandkbroadenorsteepentheslopeoftheparabola.Thelargerthevalueofk,themoreenergyisrequiredtodeformanangle(orbond)fromitsequilibriumvalue.,TorsionEnergy,Acontrolstheamplitudeofthecurvencontrolsitsperiodicityshiftstheentirecurvealongtherotationangleaxis().,Theparametersaredeterminedfromcurvefitting.Uniqueparametersfortorsionalrotationareassignedtoeachbond