分子动力学模拟方法

,Click to edit Master title style,Click to edit Master text styles,Second level,Third level,Fourth level,Fifth level,2021/3/13,*,*,单击此处编辑母版标题样式,2021/3/13,*,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,2021/3/13,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,2021/3/13,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,2021/3/13,*,分子动力学模拟方法,1957年：基于刚球势的分子動力学法Alder and Wainwright,1964年：利用Lennard-Jone势函数法对液态氩性质的模拟Rahman,1971年：模拟具有分子团簇行为的水的性质Rahman and Stillinger,1977年：约束动力学方法Rychaert, Ciccotti van Gunsteren,1980年：恒压条件下的动力学方法Andersen法、Parrinello-Rahman法,1983年：非平衡态动力学方法Gillan and Dixon,1984年: 恒温条件下的动力学方法(Berendsen et al.),1984年：恒温条件下的动力学方法Nos-Hoover法,1985年：第一原理分子動力学法Car-Parrinello法,1991年：巨正那么系综的分子动力学方法Cagin and Pettit,分子动力学简史,2021/3/13,2,粒子的运动取决于经典力学(牛顿定律F=ma,课程讲解内容：经典分子动力学,(Classical Molecular Dynamics),2021/3/13,3,分子动力学方法根底：,原理：,计算一组分子的相空间轨道，其中每个分子各自服从牛顿运动定律：,初始条件：,2021/3/13,4,分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法，可以预测纳米尺度上的材料动力学特性。,通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与原子运动路径相关的根本过程。,在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。,分子动力学方法是确定性方法，一旦初始构型和速度确定了，分子随时间所产生的运动轨迹也就确定了。,分子动力学方法特征：,2021/3/13,5,分子动力学的算法：有限差分方法,一、,Verlet,算法,粒子位置的,Taylor,展开式：,粒子位置：,粒子速度：,粒子加速度：,开场运动时需要r(t-t)：,+,缺点：,Verlet,算法处理速度非常笨拙,2021/3/13,6,Verlet,算法的表述：,算法启动,规定初始位置,规定初始速度,扰动初始位置：,计算第,n,步的力,计算第,n+1,步的位置：,计算第,n,步的速度：,重复至,2021/3/13,7,Verlet,算法程序：,Do 100 I = 1, N,RXNEWI = 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 = ( RXNEWI RXOLD(I) ) / DT2,VYI = ( RYNEWI RYOLD(I) ) / DT2,VZI = ( RZNEWI RZOLD(I) ) / DT2,RXOLD(I) = RX(I),RYOLD(I) = RY(I),RZOLD(I) = RZ(I),RX(I) = RXNEWI,RY(I) = RYNEWI,RZ(I) = RZNEWI,100 CONTINUE,2021/3/13,8,优点：,1、准确，误差O(4),2、每次积分只计算一次力,3、时间可逆,缺点：,1、速度有较大误差O(2),2、轨迹与速度无关，无法与热浴耦联,Verlet,算法的优缺点：,2021/3/13,9,二、蛙跳,(Leap-frog),算法：半步算法,1.,首先利用当前时刻的加速度，计算半个时间步长后的速度：,2.,计算下一步长时刻的位置：,3.,计算当前时刻的速度：,t-,t/2,t,t+,t/2,t+,t,t+3,t/2,t+2,t,v,r,v,开场运动时需要v(-t/2)：,2021/3/13,10,Leap-frog,算法的表述：,算法启动,规定初始位置,规定初始速度,扰动初始速度：,计算第,n,步的力,计算第,n+1/2,步的速度：,计算第,n+1,步的位置：,计算第,n,步的速度：,重复至,2021/3/13,11,Leap-frog,算法的优缺点：,优点：,1、提高准确度,2、轨迹与速度有关，可与热浴耦联,缺点：,1、速度近似,2、比Verlet算子多花时间,2021/3/13,12,三、,Velocity Verlet,算法：,等价于,优点：速度计算更加准确,2021/3/13,13,Velocity Verlet,算法的表述：,算法启动,规定初始位置,规定初始速度,计算第,n+1,步的位置：,计算第,n+1,步的力,计算第,n+1,步的速度：,重复,至,2021/3/13,14,Verlet三种形式算法的比较：,Verlet,Leap-frog,Velocity Verlet,2021/3/13,15,四、预测校正,(Predictor-Corrector),格式算法：,预测(Predictor)阶段：其根本思想是Taylor展开，,2021/3/13,16,根据新的原子位置,r,p,，可以计算获得校正后的,a,c,(,t,+,t,),定义预测误差,:,利用此预测误差，对预测出的位置、速度、加速度等量进展校正：,校正,(Corrector),阶段：,2021/3/13,17,预测阶段运动方程的变换：,定义一组矢量,:,2021/3/13,18,校正阶段运动方程的变换：,的形式：,C,0,，,C,1,，,C,2,，,C,3,的值以及,C0,，,取决于运动方程的阶数。,2021/3/13,19,一阶运动方程：,Values,c,0,c,1,c,2,c,3,c,4,c,5,3,5/12,1,1/2,4,3/8,1,3/4,1/6,5,251/720,1,11/12,1/3,1/24,6,95/288,1,25/24,35/72,5/48,1/120,2021/3/13,20,二阶运动方程之一：,Values,c,0,c,1,c,2,c,3,c,4,c,5,3,0,1,1,4,1/6,5/6,1,1/3,5,19/120,3/4,1,1/2,1/12,6,3/20,251/360,1,11/18,1/6,1/60,2021/3/13,21,二阶运动方程之二：,Values,c,0,c,1,c,2,c,3,c,4,c,5,3,0,1,1,4,1/6,5/6,1,1/3,5,19/90,3/4,1,1/2,1/12,6,3/16,251/360,1,11/18,1/6,1/60,2021/3/13,22,五、积分时间步长,t,的选择：,室温下，,t, 1 fs (femtosecond 10,-15,s),，温度越高，,t,应该减小,太长的时间步长会造成分子间的剧烈碰撞，体系数据溢出;,太短的时间步长会降低模拟过程搜索相空间的能力,2021/3/13,23,微正那么系综分子动力学(NVE MD),它是分子动力学方法的最根本系综,具有确定的粒子数,N,，能量,E,和体积,V,算法：,规定初始位置和初始速度,对运动方程积分假设干步,计算势能和动能,假设能量不等于所需要的值，对速度进展标度,重复至，直到系统平衡,2021/3/13,24,微正那么系综(NVE)MD模拟算法的流程图：,给定每个分子的初始位置,r,i,(0),和速度,v,i,(0),计算每个分子的受力,F,i,和加速度,a,i,解运动方程并求出每个分子运动一个时间步长后到达的位置所具有的速度,统计系统的热力学性质及其它物理量,统计性质不变？,打印结果，完毕,Yes,No,移动所有分子到新的位置并具有当前时刻的速度,2021/3/13,25,微正那么系综MD模拟程序F3讲解LJ, NVE：,无因次量：,2021/3/13,26,MD,模拟中几个热力学量的计算：,对于由,N,个单原子组成的系统：,动能和温度：,采用比照量：,2021/3/13,27,对于,LJ,流体：,势能：,采用比照量：,2021/3/13,28,内能：,内能由势能和动能组成：,采用比照量：,2021/3/13,29,压力：,采用比照量：,2021/3/13,30,2021/3/13,31,练习：,推导LJ流体分子间力的表达式(fx, fy, fz及其比照量)：,势能函数形式：,力：,采用比照量：,=x, y, z,2021/3/13,32,LJ,分子间的维里项：,=x, y, z,2021/3/13,33,采用比照量的运动方程形式：(以蛙跳(Leap-frog)算法为例),采用比照量：,2021/3/13,34,最终得到：,同理得到：,2021/3/13,35,速度的标度,(Velocity Scaling),：,根据能量均分原理，可知：,标度因子,:,比照量,速度标度,:,或,2021/3/13,36,微正那么系综MD模拟程序F3讲解LJ, NVE：,初始化：,READ (*,(A) TITLE ! 运行作业题目,READ (*,*) NSTEP ! 运行步数,READ (*,*) IPRINT ! 打印步数,READ (*,(A) CNFILE ! 位型文件,READ (*,*) DENS ! 比照密度,READ (*,*) RTEMP ! 比照温度,READ (*,*) RCUT ! 比照截断半径,READ (*,*) DT ! 比照时间步长,CALL READCN ( CNFILE ),2021/3/13,37,初始位型：,面心立方,(,face-centered cubic, FCC),：,每面中心有一格点,2021/3/13,38,体心立方,(body,-centered cubic, BCC),：,简单立方,(simple,cubic, SC),：,2021/3/13,39,XL,初始位型：,面心立方,(,FCC) (,程序,F23),NC=(REAL(N)/4.0)*(1.0/3.0),XL = 1.0 / REAL ( NC ),Y = 0.5 * XL,R(1)=(0, 0, 0) R(2)=(0, Y, Y),R(3)=(Y, 0, Y) R(4)=(Y, Y, 0),M = 0,DO 10 I = 1, NC,DO 10 J = 1, NC,DO 10 K = 1, NC,DO 11 IJ = 1, 4,RX(IJ+M)=RX(IJ) + XL*(K-1),RY(IJ+M)=RY(IJ) + XL*(J -1),RZ(IJ+M)=RZ(IJ) + XL*(I -1),11 CONTINUE,M = M + 4,10 CONTINUE,2021/3/13,40,DO 100 I = 1, N,100 CONTINUE,将模拟盒子的中心移到原点：,2021/3/13,41,初始速度：,简单的选择：,V,random(-0.5, 0.5),=x, y, z,标度因子,:,速度标度,:,2021/3/13,42,FACTOR = SQRT ( RTEMP ),DO 100 I = 1, N,VX(I) = FACTOR * ( RANF(DUMMY) - 0.5 ),VY(I) = FACTOR * ( RANF(DUMMY) - 0.5 ),VZ(I) = FACTOR * ( RANF(DUMMY) - 0.5 ),CONTINUE,随机安排初始速度：,2021/3/13,43,标度初始速度：,DO 200 I = 1, N,SUMKX = SUMKX + VX(I)*2,SUMKY = SUMKY + VY(I)*2,SUMKZ = SUMKZ + VZ(I)*2,200 CONTINUE,BEITAX =SQRT(RTEMP/SUMKX),BEITAY =SQRT(RTEMP/SUMKY),BEITAZ =SQRT(RTEMP/SUMKZ),DO 300 I = 1, N,VX(I) = VX(I) * BEITAX,VY(I) = VY(I) * BEITAY,VZ(I) = VZ(I) * BEITAZ,300 CONTINUE,标度因子,:,2021/3/13,44,DO 200 I = 1, N,SUMX = SUMX + VX(I),SUMY = SUMY + VY(I),SUMZ = SUMZ + VZ(I),CONTINUE,SUMX = SUMX / REAL ( N ),SUMY = SUMY / REAL ( N ),SUMZ = SUMZ / REAL ( N ),DO 300 I = 1, N,VX(I) = VX(I) - SUMX,VY(I) = VY(I) - SUMY,VZ(I) = VZ(I) - SUMZ,300 CONTINUE,控制体系的总动量为零：,2021/3/13,45,从,Maxwell,分布中抽样：,x,x,d,x,高斯,(Gauss),分布,:,对于等几率随机试验,(Bernoulli,试验,),，大量的试验结果满足高斯分布,2021/3/13,46,麦克斯韦速度分布定律,:,由于：,=x, y, z,单位体积的分子再每个分量上的速度分布实际上就是高斯分布。,2021/3/13,47,从,Maxwell,分布中抽样：,高斯,(Gauss),分布的随机数生成方法,:,生成随机数：,i,，,i=1, 2, , 12,2021/3/13,48,DO 10 I = 1, 12,SUM = SUM + RANF ( DUMMY ),10 CONTINUE,R2 = R * R,GAUSS = ( A9 * R2 + A7 ) * R2 + A5 ) * R2 + A3 ) * R2 +A1 ) * R,高斯,(Gauss),分布的随机数生成,(,程序,F24),2021/3/13,49,FACTOR = SQRT ( RTEMP ),DO 100 I = 1, N,VX(I) = FACTOR * GAUSS(DUMMY),VY(I) = FACTOR * GAUSS(DUMMY),VZ(I) = FACTOR * GAUSS(DUMMY),CONTINUE,控制总动量为零：同前面一样处理。,从,Maxwell,分布中抽样分布中随机安排初始速度：,2021/3/13,50,微正那么系综MD模拟程序F3讲解(LJ, NVE) ：,量纲变换：,SIGMA = ( DENS / REAL ( N ) ) * ( 1.0 / 3.0 ),RCUT = RCUT * SIGMA,DT,DT * SIGMA,DENS = DENS / ( SIGMA * 3 ),模拟盒子的边长为,1,2021/3/13,51,长程校正：,微正那么系综MD模拟程序F3讲解(LJ, NVE) ：,SR3 = ( SIGMA / RCUT ) * 3,SR9 = SR3 * 3,SIGCUB = SIGMA * 3,VLRC = ( 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 ),2021/3/13,52,算法：算法启动,微正那么系综MD模拟程序F3讲解(LJ, NVE) ：,CALL FORCE ( -DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ),CALL MOVE ( -DT ),CALL FORCE ( -DT, SIGMA, RCUT, V, VC, W ),CALL FORCE ( DT, SIGMA, RCUT, V, VC, W ),CALL KINET ( OLDK ),CALL MOVE ( DT ),CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ),CALL KINET ( NEWK ),2021/3/13,53,算法：差分格式：,SR2 = SIGSQ / RIJSQ,VIJ = 4.0 * ( SR12 - SR6 ),WIJ = 24.0 * ( 2.0 * SR12 - SR6 ),VELIJ = WIJ * DT / RIJSQ,DVX = VELIJ * RXIJ,.,VXI = VXI + DVX,.,VX(J) = VX(J) DVX,V = V + VIJ,W = W + WIJ,CALL MOVE(DT),2021/3/13,54,DO 1000 I = 1, N,RX(I) = RX(I) + VX(I) * DT,RY(I) = RY(I) + VY(I) * DT,RZ(I) = RZ(I) + VZ(I) * DT,1000,CONTINUE,MOVE(DT):,2021/3/13,55,速度的标定只用于平衡阶段,DO 200 I = 1, N,SUMK = SUMKX + VX(I)*2+VY(I)*2 + VZ(I)*2,200 CONTINUE,BEITA =SQRT( 3.0 * RTEMP / SUMK ),DO 300 I = 1, N,VX(I) = VX(I) * BEITA,VY(I) = VY(I) * BEITA,VZ(I) = VZ(I) * BEITA,300 CONTINUE,2021/3/13,56,正那么系综分子动力学(NVT MD),具有确定的粒子数,N,，温度,T,和体积,V,速度的直接标度,热浴方法,(Andersen Thermostat),约束方法,(,阻尼力方法,),系统扩展方法,(Extended Systems Method),问题的关键：温度的约束,Nose-Hoover,方法,2021/3/13,57,一、热浴方法,(Andersen Thermostat),引入一个与虚拟粒子碰撞的随机力,想象系统浸在热浴当中,系统和热浴间的相互作用强度由随机碰撞的频率决定,碰撞的几率等于,Nu,dt,如果一个粒子经历碰撞，它的速度将从约束温度下的,Maxwell,分布中随机抽取,总能量和总动量均不守恒,2021/3/13,58,二、约束方法,是等动能,(,Iso-Kinetics,),分子动力学方法,系统的运动方程为：,引入阻尼系数,以保证将温度约束在恒定值,根据高斯最小约束原理：,2021/3/13,59,三、,Nose-Hoover,扩展方法,根本思想：,设想原系统与一个耦合系统共同组成一个扩展系统，允许热流在原系统和耦合系统之间交换。,Q,：等效质量,S:,扩展坐标变量,： 热力学阻尼系数,L,： 扩展系统的自由度,2021/3/13,60,Predictor-corrector algorithm is straightforward,Verlet algorithm is feasible, but tricky to implement,积分方案：,update of,x,depends on p,update of p depends on,x,2021/3/13,61,Nos-Hoover,方法正确地描述了,NVT,系综中的动量和位型，而等动能方法只正确地描述了后者。,扩展系统的哈密顿量,Hamiltonian,守恒：,2021/3/13,62,1.,复杂分子体系的势能函数形式：,Potential Energy = Stretching Energy +,Bending Energy +,Torsion Energy +,Non-Bonded Interaction Energy,这些方程与描述原子或键各种不同行为的参数就构,成了力场，,force-field.,UFF, OPLS, Amber, CVFF, Compass,分子模拟方法补充介绍：,2021/3/13,63,Bond Stretching Energy,k,b,is the spring constant of the bond.,r,0,is the bond length at equilibrium.,Unique k,b,and r,0,assigned for each bond pair, i.e. C-C, O-H,2021/3/13,64,Bending Energy,k,is the spring constant of the bend.,0,is the bond length at equilibrium.,Unique parameters for angle bending are assigned to each bonded triplet of atoms based on their types (e.g. C-C-C, C-O-C, C-C-H, etc.),2021/3/13,65,The “Hookeian potential,k,b,and,k,broaden or steepen the slope of the parabola.,The larger the value of,k, the more energy is required to deform an angle (or bond) from its equilibrium value.,2021/3/13,66,Torsion Energy,A,controls the amplitude of the curve,n,controls its periodicity,shifts the entire curve along the rotation angle axis (,).,The parameters are determined from curve fitting.,Unique parameters for torsional rotation are assigned to each bonded quartet of atoms based on their types (e.g. C-C-C-C, C-O-C-N, H-C-C-H, etc.),2021/3/13,67,Non-bonded Energy,A,determines the degree the attractiveness,B,determines the degree of repulsion,q,is the charge,A,determines the degree the attractiveness,B,determines the degree of repulsion,q,is the charge,2021/3/13,68,2.,长程静电力的处理方式：,Ewald 加和Ewald Summation,2021/3/13,69,有关,Ewald,加和的说明：,Basic form requires an O(N,2,) calculation,efficiency can be introduced to reduce to O(N,3/2,),good value of a is 5L, but should check for given application,can be extended to sum point dipoles,Other methods are in common use,reaction field,particle-particle/particle mesh (better for larger systems),fast multipole (better for larger systems),2021/3/13,70,3. Verlet,近邻表,(Verlet Neighbor lists),：,Keep an expanded neighbor list,r,l,r,c,so that neighbor pairs need not be calculated every timestep,r,l:,chosen such that update lists,every 10-20 timesteps,Good for,N, 1000,2021/3/13,72,谢谢大家！,