分子动力学模拟方法

分子动力学模拟方法第一页，共七十二页，2022年，8月28日1957年：基于刚球势的分子動力学法（AlderandWainwright）1964年：利用Lennard-Jone势函数法对液态氩性质的模拟（Rahman）1971年：模拟具有分子团簇行为的水的性质（RahmanandStillinger）1977年：约束动力学方法（Rychaert,Ciccotti&Berendsen;vanGunsteren）1980年：恒压条件下的动力学方法（Andersen法、Parrinello-Rahman法）1983年：非平衡态动力学方法（GillanandDixon）1984年:恒温条件下的动力学方法(Berendsenetal.)1984年：恒温条件下的动力学方法（Nosé-Hoover法）1985年：第一原理分子動力学法（→Car-Parrinello法）1991年：巨正则系综的分子动力学方法（CaginandPettit）分子动力学简史第二页，共七十二页，2022年，8月28日粒子的运动取决于经典力学(牛顿定律（F=ma）课程讲解内容：经典分子动力学

(ClassicalMolecularDynamics)第三页，共七十二页，2022年，8月28日分子动力学方法基础：原理：计算一组分子的相空间轨道，其中每个分子各自服从牛顿运动定律：初始条件：第四页，共七十二页，2022年，8月28日分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法，可以预测纳米尺度上的材料动力学特性。通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与原子运动路径相关的基本过程。在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。分子动力学方法是确定性方法，一旦初始构型和速度确定了，分子随时间所产生的运动轨迹也就确定了。分子动力学方法特征：第五页，共七十二页，2022年，8月28日分子动力学的算法：有限差分方法一、Verlet算法粒子位置的Taylor展开式：粒子位置：粒子速度：粒子加速度：开始运动时需要r(t-Δt)：+缺点：Verlet算法处理速度非常笨拙第六页，共七十二页，2022年，8月28日Verlet算法的表述：算法启动规定初始位置规定初始速度扰动初始位置：计算第n步的力计算第n+1步的位置：计算第n步的速度：重复④至⑥第七页，共七十二页，2022年，8月28日Verlet算法程序：Do100I=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))/DT2VYI=(RYNEWI–RYOLD(I))/DT2VZI=(RZNEWI–RZOLD(I))/DT2RXOLD(I)=RX(I)RYOLD(I)=RY(I)RZOLD(I)=RZ(I)RX(I)=RXNEWIRY(I)=RYNEWIRZ(I)=RZNEWI100CONTINUE第八页，共七十二页，2022年，8月28日优点：1、精确，误差O(Δ4)2、每次积分只计算一次力3、时间可逆缺点：1、速度有较大误差O(Δ2)2、轨迹与速度无关，无法与热浴耦联Verlet算法的优缺点：第九页，共七十二页，2022年，8月28日二、蛙跳(Leap-frog)算法：半步算法1.首先利用当前时刻的加速度，计算半个时间步长后的速度：2.计算下一步长时刻的位置：3.计算当前时刻的速度：t-Δt/2tt+Δt/2t+Δtt+3Δt/2t+2Δtvrv开始运动时需要v(-Δt/2)：第十页，共七十二页，2022年，8月28日Leap-frog算法的表述：算法启动规定初始位置规定初始速度扰动初始速度：计算第n步的力计算第n+1/2步的速度：计算第n+1步的位置：计算第n步的速度：重复④至⑦第十一页，共七十二页，2022年，8月28日Leap-frog算法的优缺点：优点：1、提高精确度2、轨迹与速度有关，可与热浴耦联缺点：1、速度近似2、比Verlet算子多花时间第十二页，共七十二页，2022年，8月28日三、VelocityVerlet算法：等价于优点：速度计算更加准确第十三页，共七十二页，2022年，8月28日VelocityVerlet算法的表述：算法启动规定初始位置规定初始速度计算第n+1步的位置：计算第n+1步的力计算第n+1步的速度：重复③至⑤第十四页，共七十二页，2022年，8月28日Verlet三种形式算法的比较：VerletLeap-frogVelocityVerlet第十五页，共七十二页，2022年，8月28日四、预测－校正(Predictor-Corrector)格式算法：预测(Predictor)阶段：其基本思想是Taylor展开，第十六页，共七十二页，2022年，8月28日根据新的原子位置rp，可以计算获得校正后的ac(t+t),定义预测误差:利用此预测误差，对预测出的位置、速度、加速度等量进行校正：校正(Corrector)阶段：第十七页，共七十二页，2022年，8月28日预测阶段运动方程的变换：定义一组矢量:第十八页，共七十二页，2022年，8月28日校正阶段运动方程的变换：的形式：C0，C1，C2，C3的值以及C0，取决于运动方程的阶数。第十九页，共七十二页，2022年，8月28日一阶运动方程：Valuesc0c1c2c3c4c535/1211/2

43/813/41/65251/720111/121/31/24

695/288125/2435/725/481/120第二十页，共七十二页，2022年，8月28日二阶运动方程之一：Valuesc0c1c2c3c4c5301141/65/611/3519/1203/411/21/1263/20251/360111/181/61/60第二十一页，共七十二页，2022年，8月28日二阶运动方程之二：Valuesc0c1c2c3c4c5301141/65/611/3519/903/411/21/1263/16251/360111/181/61/60第二十二页，共七十二页，2022年，8月28日五、积分时间步长t的选择：室温下，∆t≈1fs(femtosecond10-15s)，温度越高，∆t

应该减小太长的时间步长会造成分子间的激烈碰撞，体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力第二十三页，共七十二页，2022年，8月28日微正则系综分子动力学(NVEMD)它是分子动力学方法的最基本系综具有确定的粒子数N，能量E和体积V算法：规定初始位置和初始速度对运动方程积分若干步计算势能和动能若能量不等于所需要的值，对速度进行标度重复②至④，直到系统平衡第二十四页，共七十二页，2022年，8月28日微正则系综(NVE)MD模拟算法的流程图：给定每个分子的初始位置ri(0)和速度vi(0)计算每个分子的受力Fi和加速度ai解运动方程并求出每个分子运动一个时间步长后到达的位置所具有的速度统计系统的热力学性质及其它物理量统计性质不变？

打印结果，结束YesNo移动所有分子到新的位置并具有当前时刻的速度第二十五页，共七十二页，2022年，8月28日微正则系综MD模拟程序F3讲解（LJ,NVE）：无因次量：第二十六页，共七十二页，2022年，8月28日MD模拟中几个热力学量的计算：对于由N个单原子组成的系统：动能和温度：采用对比量：第二十七页，共七十二页，2022年，8月28日对于LJ流体：势能：采用对比量：第二十八页，共七十二页，2022年，8月28日内能：内能由势能和动能组成：采用对比量：第二十九页，共七十二页，2022年，8月28日压力：采用对比量：第三十页，共七十二页，2022年，8月28日第三十一页，共七十二页，2022年，8月28日练习：推导LJ流体分子间力的表达式(fx,fy,fz及其对比量)：势能函数形式：力：采用对比量：=x,y,z第三十二页，共七十二页，2022年，8月28日LJ分子间的维里项：=x,y,z第三十三页，共七十二页，2022年，8月28日采用对比量的运动方程形式：

(以蛙跳(Leap-frog)算法为例)采用对比量：第三十四页，共七十二页，2022年，8月28日最终得到：同理得到：第三十五页，共七十二页，2022年，8月28日速度的标度(VelocityScaling)：根据能量均分原理，可知：标度因子:对比量速度标度:或第三十六页，共七十二页，2022年，8月28日微正则系综MD模拟程序F3讲解（LJ,NVE）：初始化：READ(*,‘(A)’)TITLE!运行作业题目READ(*,*)NSTEP!运行步数READ(*,*)IPRINT!打印步数READ(*,‘(A)’)CNFILE!位型文件READ(*,*)DENS!对比密度READ(*,*)RTEMP

!对比温度

READ(*,*)RCUT!对比截断半径READ(*,*)DT!对比时间步长CALLREADCN(CNFILE)第三十七页，共七十二页，2022年，8月28日初始位型：面心立方(face-centeredcubic,FCC)：每面中心有一格点

第三十八页，共七十二页，2022年，8月28日体心立方(body-centeredcubic,BCC)：简单立方(simplecubic,SC)：第三十九页，共七十二页，2022年，8月28日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第四十页，共七十二页，2022年，8月28日DO100I=1,NRX(I)=RX(I)-0.5RY(I)=RY(I)-0.5RZ(I)=RZ(I)-0.5100CONTINUE将模拟盒子的中心移到原点：第四十一页，共七十二页，2022年，8月28日初始速度：简单的选择：V＝random(-0.5,0.5)=x,y,z标度因子:速度标度:第四十二页，共七十二页，2022年，8月28日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随机安排初始速度：第四十三页，共七十二页，2022年，8月28日标度初始速度：SUMKX=0.0SUMKY=0.0SUMKZ=0.0DO200I=1,NSUMKX=SUMKX+VX(I)**2SUMKY=SUMKY+VY(I)**2SUMKZ=SUMKZ+VZ(I)**2200CONTINUEBEITAX=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标度因子:第四十四页，共七十二页，2022年，8月28日SUMX=0.0SUMY=0.0SUMZ=0.0DO200I=1,NSUMX=SUMX+VX(I)SUMY=SUMY+VY(I)SUMZ=SUMZ+VZ(I)CONTINUESUMX=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控制体系的总动量为零：第四十五页，共七十二页，2022年，8月28日从Maxwell分布中抽样：xxdx高斯(Gauss)分布:对于等几率随机试验(Bernoulli试验)，大量的试验结果满足高斯分布第四十六页，共七十二页，2022年，8月28日麦克斯韦速度分布定律:由于：=x,y,z单位体积的分子再每个分量上的速度分布实际上就是高斯分布。第四十七页，共七十二页，2022年，8月28日从Maxwell分布中抽样：高斯(Gauss)分布的随机数生成方法:生成随机数：i，i=1,2,…,12第四十八页，共七十二页，2022年，8月28日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)第四十九页，共七十二页，2022年，8月28日FACTOR=SQRT(RTEMP)DO100I=1,NVX(I)=FACTOR

*GAUSS(DUMMY)VY(I)=FACTOR

*GAUSS(DUMMY)VZ(I)=FACTOR

*GAUSS(DUMMY)CONTINUE控制总动量为零：同前面一样处理。从Maxwell分布中抽样分布中随机安排初始速度：第五十页，共七十二页，2022年，8月28日微正则系综MD模拟程序F3讲解(LJ,NVE)

：量纲变换：SIGMA=(DENS/REAL(N))**(1.0/3.0)RCUT=RCUT*SIGMADT＝DT*

SIGMADENS=DENS/(SIGMA**3)模拟盒子的边长为1第五十一页，共七十二页，2022年，8月28日长程校正：微正则系综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)第五十二页，共七十二页，2022年，8月28日算法：算法启动微正则系综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)第五十三页，共七十二页，2022年，8月28日算法：差分格式：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)–DVX………V=V+VIJW=W+WIJCALLMOVE(DT)第五十四页，共七十二页，2022年，8月28日DO1000I=1,N

RX(I)=RX(I)+VX(I)*DT

RY(I)=RY(I)+VY(I)*DT

RZ(I)=RZ(I)+VZ(I)*DT1000

CONTINUEMOVE(DT):第五十五页，共七十二页，2022年，8月28日速度的标定（只用于平衡阶段）SUMK=0.0DO200I=1,NSUMK=SUMKX+VX(I)**2+VY(I)**2+VZ(I)**2200CONTINUEBEITA=SQRT(3.0*RTEMP/SUMK)DO300I=1,NVX(I)=VX(I)*BEITAVY(I)=VY(I)*BEITAVZ(I)=VZ(I)*BEITA300CONTINUE第五十六页，共七十二页，2022年，8月28日正则系综分子动力学(NVTMD)具有确定的粒子数N，温度T和体积V速度的直接标度热浴方法(AndersenThermostat)约束方法(阻尼力方法)系统扩展方法(ExtendedSystemsMethod)问题的关键：温度的约束Nose-Hoover方法第五十七页，共七十二页，2022年，8月28日一、热浴方法(AndersenThermostat)引入一个与虚拟粒子碰撞的随机力想象系统浸在热浴当中系统和热浴间的相互作用强度由随机碰撞的频率决定碰撞的几率等于Nu×dt如果一个粒子经历碰撞，它的速度将从约束温度下的Maxwell分布中随机抽取总能量和总动量均不守恒第五十八页，共七十二页，2022年，8月28日二、约束方法是等动能(Iso-Kinetics)分子动力学方法系统的运动方程为：引入阻尼系数以保证将温度约束在恒定值根据高斯最小约束原理：第五十九页，共七十二页，2022年，8月28日三、Nose-Hoover扩展方法基本思想：设想原系统与一个耦合系统共同组成一个扩展系统，允许热流在原系统和耦合系统之间交换。Q：等效质量S:扩展坐标变量：热力学阻尼系数L：扩展系统的自由度第六十页，共七十二页，2022年，8月28日Predictor-correctoralgorithmisstraightforwardVerletalgorithmisfeasible,buttrickytoimplement积分方案：updateofxdependsonpupdateofpdependsonx第六十一页，共七十二页，2022年，8月28日Nosé-Hoover方法正确地描述了NVT系综中的动量和位型，而等动能方法只正确地描述了后者。扩展系统的哈密顿量Hamiltonian守恒：第六十二页，共七十二页，2022年，8月28日1.复杂分子体系的势能函数形式：PotentialEnergy=StretchingEnergy+BendingEnergy+TorsionEnergy+ Non-BondedInteractionEnergy这些方程与描述原子或键各种不同行为的参数就构成了力场，

force-field.UFF,OPLS,Amber,CVFF,Compass分子模拟方法补充介绍：第六十三页，共七十二页，2022年，8月28日BondStretchingEnergykbisthespringconstantofthebond.r0isthebondlengthatequilibrium.Uniquekbandr0assignedforeachbondpair,i.e.C-C,O-H第六十四页，共七十二页，2022年，8月28日BendingEnergykisthespringconstantofthebend.0isthebondlengthatequilibrium.Uniqueparametersforanglebendingareassignedtoeachbondedtripletofatomsbasedontheirtypes(e.g.C-C-C,C-O-C,C-C-H,etc.)第六十五页，共七十二页，2022年，8月28日The“Hookeian”potentialkbandkbroadenorsteepentheslopeoftheparabola.Thelargerthevalueofk,themoreenergyisrequiredtodeformanangle(orbond)fromitsequilibriumvalue.第六十六页，共七十二页，2022年，8月28日TorsionEnergyAcontrolstheamplitudeofthecurvencontrolsitsperiodicityshiftstheentirecurvealongtherotationangleaxis().Theparametersaredeterminedfromcurvefitting.Uniqueparametersfortorsionalrotationareassignedtoeachbondedquartetofatomsbasedontheirtypes(e.