1、第四章第四章 分子动力学模拟方法分子动力学模拟方法1957年年：基于刚球势的分子動力学法（：基于刚球势的分子動力学法（Alder and Wainwright）1964年年：利用：利用Lennard-Jone势函数法对液态势函数法对液态氩氩性质的模拟（性质的模拟（Rahman）1971年年：模拟具有分子：模拟具有分子团簇团簇行为的行为的水水的性质（的性质（Rahman and Stillinger）1977年年：约束动力学方法（：约束动力学方法（Rychaert,Ciccotti&Berendsen;van Gunsteren）1980年年：恒压条件下的动力学方法（：恒压条件下的动力学方法（A
2、ndersen法、法、Parrinello-Rahman法）法）1983年年：非平衡态动力学方法（：非平衡态动力学方法（Gillan and Dixon）1984年年:恒温条件下的动力学方法恒温条件下的动力学方法(Berendsen et al.)1984年年：恒温条件下的动力学方法（：恒温条件下的动力学方法（Nos-Hoover法）法）1985年年：第一原理分子動力学法（：第一原理分子動力学法（Car-Parrinello法）法）1991年年：巨正则系综巨正则系综的分子动力学方法（的分子动力学方法（Cagin and Pettit）分子动力学简史分子动力学简史粒子的运动取决于经典力学粒子的运
3、动取决于经典力学(牛顿定律（牛顿定律（F=ma）课程讲解内容：经典分子动力学课程讲解内容：经典分子动力学(Classical Molecular Dynamics)原理：原理：计算一组分子的相空间轨道，其中每个分子各自服从计算一组分子的相空间轨道，其中每个分子各自服从牛顿运动定律：牛顿运动定律：11112)(21NiNijijNiiirUmpHiiiiimdtdmvrp 111111)()(NiNijijNiNijijirUrdtdijrFp)0(0itirr)0(0itidtvdr初始条件：初始条件：粒子位置的粒子位置的Taylor展开式：展开式：粒子位置粒子位置：粒子速度粒子速度：t(t)
4、t)-(t(t)2 t)(t2iarrriii t2t)-(tt)(t(t)iiirrv粒子加速度：粒子加速度：3i2iit(t)61 t(t)21 t(t)(t)t)(tbavrrii3i2iit(t)61 t(t)21 t(t)(t)t)(tbavrrii (t)(t)imiiFa开始运动时需要开始运动时需要r(t-t)：t(0)(0)t)(ivrr+n 算法启动算法启动 规定初始位置规定初始位置 规定初始速度规定初始速度 扰动初始位置：扰动初始位置：计算第计算第n步的力步的力 计算第计算第n+1步的位置：步的位置：u 计算第计算第n步的速度：步的速度：u 重复重复至至t(0)(0)t)(
5、ivrr t(t)t)-(t(t)2 t)(t2iarrriii t2t)-(tt)(t(t)iiirrvDo 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
6、 RY(I)=RYNEWI RZ(I)=RZNEWI100 CONTINUEt(t)t)21-(t t)21(tiavviit t)21(t(t)t)(tivrrii2t)21-(tt)21(t(t)iiivvv1.首先利用当前时刻的加速度，计算半个时间步长后的速度：首先利用当前时刻的加速度，计算半个时间步长后的速度：2.计算下一步长时刻的位置：计算下一步长时刻的位置：3.计算当前时刻的速度：计算当前时刻的速度：t-t/2tt+t/2t+tt+3t/2t+2tvrv开始运动时需要开始运动时需要v(-t/2)：t/2(0)(0)t/2)(iavvn 算法启动算法启动 规定初始位置规定初始位置 规
7、定初始速度规定初始速度 扰动初始速度：扰动初始速度：计算第计算第n步的力步的力 计算第计算第n+1/2步的速度：步的速度：计算第计算第n+1步的位置：步的位置：计算第计算第n步的速度：步的速度：重复重复至至t/2(0)(0)t/2)(iavvt(t)t)21-(t t)21(tiavviit t)21(t(t)t)(tivrrii2t)21-(tt)21(t(t)iiivvv2iit(t)21 t(t)(t)t)(tavrriitt)(t(t)21(t)t)(tiiaavviit(t)21(t)t)21(tiavviitt)(t 21 t)21(t t)(tiavviit t)21(t(t)t
8、)(tivrrii等价于等价于n 算法启动算法启动 规定初始位置规定初始位置 规定初始速度规定初始速度 计算第计算第n+1步的位置：步的位置：计算第计算第n+1步的力步的力 计算第计算第n+1步的速度：步的速度：重复重复至至2iit(t)21 t(t)(t)t)(tavrriitt)(t(t)21(t)t)(tiiaavviiVerletLeap-frogVelocity Verlet1.预测预测(Predictor)阶段：其基本思想是阶段：其基本思想是Taylor展开，展开，)()()()()()(21)()()()(61)(21)()()(232tttttttttttttttttttttt
9、ttppppbbbaabavvbavrr根据新的原子位置根据新的原子位置r rp p，可以计算获得校正后的，可以计算获得校正后的a ac c(t t+t t),),定义预测误差定义预测误差:)()()(ttttttpcaaa)()()()()()()()()()()()(3210ttctttttctttttctttttctttpcpcpcpcabbaaaavvarr利用此预测误差，对预测出的位置、速度、加速度等量进行校正：利用此预测误差，对预测出的位置、速度、加速度等量进行校正：2.校正校正(Corrector)阶段：阶段：u预测阶段运动方程的变换：预测阶段运动方程的变换：定义一组矢量定义一组
10、矢量:333322221061d)(d21d)(dd)(d)(ttttttttttrrrrrrrr)()()()(1000310032101111)()()()(32103210ttttttttttttpppprrrrrrrr)()()()()()(21)()()()(61)(21)()()(232ttttttttttttttttttttttttppppbbbaabavvbavrru校正阶段运动方程的变换：校正阶段运动方程的变换：rrrrrrrrr)()()()()()()()(321032103210ccccttttttttttttttttppppccccr的形式：的形式：C0，C1，C2，
11、C3的值以及的值以及C0，取决于运动方程的阶数。取决于运动方程的阶数。一阶运动方程：一阶运动方程：)(rrf)()(11ttttpcrrrValuesc0c1c2c3c4c535/1211/2 43/813/41/65251/720111/121/31/24 695/288125/2435/725/481/120 二阶运动方程之一：二阶运动方程之一：)(rrf)()(22ttttpcrrrValuesc0c1c2c3c4c5301141/65/611/3519/1203/411/21/1263/20251/360111/181/61/60 二阶运动方程之二：二阶运动方程之二：),(rrr f)
12、()(22ttttpcrrrValuesc0c1c2c3c4c5301141/65/611/3519/903/411/21/1263/16251/360111/181/61/60t室温下，室温下，t 1 fs(femtosecond 10-15s)，温，温度越高，度越高，t 应该减小应该减小n 太长的时间步长会造成分子间的激烈碰撞，体系数据溢出太长的时间步长会造成分子间的激烈碰撞，体系数据溢出;n 太短的时间步长会降低模拟过程搜索相空间的能力太短的时间步长会降低模拟过程搜索相空间的能力 u 它是分子动力学方法的最基本系综它是分子动力学方法的最基本系综u 具有确定的粒子数具有确定的粒子数N，能量
13、，能量E和体积和体积Vu 算法：算法：规定初始位置和初始速度规定初始位置和初始速度 对运动方程积分若干步对运动方程积分若干步 计算势能和动能计算势能和动能 若能量不等于所需要的值，对速度进行标度若能量不等于所需要的值，对速度进行标度 重复重复至，直到系统平衡至，直到系统平衡微正则系综微正则系综(NVE)MD模拟算法的流程图：模拟算法的流程图：给定每个分子的初始位置给定每个分子的初始位置ri(0)和速度和速度vi(0)计算每个分子的受力计算每个分子的受力Fi和加速度和加速度ai解运动方程并求出每个分子运动一个时间步解运动方程并求出每个分子运动一个时间步长后到达的位置所具有的速度长后到达的位置所具
14、有的速度统计系统的热力学性质及其它物理量统计系统的热力学性质及其它物理量统计性质不变？统计性质不变？打印结果，结束打印结果，结束YesNo移动所有分子到新的位置并具有当前时刻的移动所有分子到新的位置并具有当前时刻的速度速度微正则系综微正则系综MD模拟程序模拟程序F3讲解（讲解（LJ,NVE）：）：无因次量：无因次量：/*rr 3*/*EE/*ff/*3PP/*kTT 2/12*)/(mtt 2/1*)/(vvmMD模拟中几个热力学量的计算：模拟中几个热力学量的计算：对于由对于由N个单原子组成的系统：个单原子组成的系统：动能和温度：动能和温度：TNkmEBNiK23v2112i采用对比量：采用对
15、比量：/*kTT 2/1*)/(vvm/*EE*12*i*23v21NTENiK21114)ru 2)(r(drNrUUUUcrNiNijijclrcc 对于对于LJ流体：流体：3933138 cclrcrrNU势能：势能：4)(612rrrUc/*rr 3*/*EE 采用对比量：采用对比量：114)(6*12*rrrUc3*9*113138 ccrrNUlrc内能：内能：KEUE内能由势能和动能组成：内能由势能和动能组成：采用对比量：采用对比量：*KEUE压力：压力：393232316 cclrcrrPlrccBPVWTkP/*rr 3*/*EE 采用对比量：采用对比量：/*3PP/*ff
16、111)(31NiNijijccrWWijijijijijijcrUrWrrfr)(31)(*lrccPVWTP3*9*2*1132316 cclrcrrP 111*)(31NiNijijccrWW*)()(jjjjjjcrUrWiiiiiirrfr练习：练习：推导推导LJ流体分子间力的表达式流体分子间力的表达式(fx,fy,fz及其对比量及其对比量)：4)(612rrrU势能函数形式：势能函数形式：力：力：r)()(rUrf2612 2184 rrrrrrrUrUf/*rr 采用对比量：采用对比量：/*ff=x,y,z2*6*12*121184 rrrrf*)(fr rWcLJ分子间的维里项
17、：分子间的维里项：2*6*12*121184 rrrrf=x,y,z 121184 )(6*12*rrrWc2*)(rrrWfct(t)t)21-(t t)21(timFvviit t)21(t(t)t)(tivrrii2t)21-(tt)21(t(t)iiivvv(t)ia/*ff2/12*)/(mtt 2/1*)/(vvm采用对比量：采用对比量：/*rr 2/12*2/1*2/1*t )(t )t21-(t )t21(tmmmm*i*i*iFvv2/12/12*t )(t)t21-(t )t21(tmmm*i*i*iFvv*t)(t)t21-(t )t21(t*i*i*iFvv最终得到：最
18、终得到：*t)t21(t)(t )t(t*i*i*ivrr2)t21-(t)t21(t(t)*i*i*ivvv同理得到：同理得到：mTkvBx2根据能量均分原理，可知：根据能量均分原理，可知：2/12*2/122/122/1333vTvmTkvNmTNkTTreqrefBrefBref标度因子标度因子:对比量对比量速度标度速度标度:oldnewvv*old*newvv 或或微正则系综微正则系综MD模拟程序模拟程序F3讲解（讲解（LJ,NVE）：）：初始化：初始化：READ(*,(A)TITLE !运行作业题目运行作业题目READ(*,*)NSTEP !运行步数运行步数READ(*,*)IPRI
19、NT !打印步数打印步数READ(*,(A)CNFILE !位型文件位型文件READ(*,*)DENS !对比密度对比密度READ(*,*)RTEMP !对比温度对比温度 READ(*,*)RCUT !对比截断半径对比截断半径READ(*,*)DT !对比时间步长对比时间步长CALL READCN(CNFILE)初始位型：初始位型：u面心立方面心立方(face-centered cubic,FCC)：每面中心有一格点每面中心有一格点 u体心立方体心立方(body-centered cubic,BCC)：u简单立方简单立方(simple cubic,SC)：XL初始位型：初始位型：面心立方面心立
20、方(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=0DO 10 I=1,NCDO 10 J=1,NCDO 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+410 CONTINUEDO 100 I=1,N RX(I)=RX(I)-0.5 RY(I)=
21、RY(I)-0.5 RZ(I)=RZ(I)-0.5100 CONTINUE将模拟盒子的中心移到原点：将模拟盒子的中心移到原点：初始速度：初始速度：n 简单的选择：简单的选择：V random(-0.5,0.5)=x,y,zmTkvB22/1*2/12*2/122/1TTvTvmTkTTreqreqrefBref标度因子标度因子:*old*newvv 速度标度速度标度:*2*TvFACTOR=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(
22、DUMMY)-0.5)100 CONTINUE*2*Tv随机安排初始速度：随机安排初始速度：标度初始速度：标度初始速度：SUMKX=0.0SUMKY=0.0SUMKZ=0.0DO 200 I=1,N SUMKX=SUMKX+VX(I)*2 SUMKY=SUMKY+VY(I)*2 SUMKZ=SUMKZ+VZ(I)*2200 CONTINUEBEITAX=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)=
23、VZ(I)*BEITAZ300 CONTINUE2/1*2/12*TTvTreqreq标度因子标度因子:*old*newvv SUMX=0.0SUMY=0.0SUMZ=0.0DO 200 I=1,N SUMX=SUMX+VX(I)SUMY=SUMY+VY(I)SUMZ=SUMZ+VZ(I)200 CONTINUESUMX=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)-SUMZ300 CONTINUE控制体系的总动量为零：控制体系的总
24、动量为零：n 从从Maxwell分布中抽样：分布中抽样：x)(xP0 x xdx高斯高斯(Gauss)分布分布:对于等几率随机试验对于等几率随机试验(Bernoulli试验试验)，大量的试验结果满足高斯分布大量的试验结果满足高斯分布 222)(2exp)(21xxxxxP麦克斯韦速度分布定律麦克斯韦速度分布定律:dvekTmNdvvfTkmvB22122mTkvB2由于：由于：222)(2exp)(21xxxxxP=x,y,z单位体积的分子再每个分量上的速度分布实际上就是单位体积的分子再每个分量上的速度分布实际上就是高斯分布。高斯分布。n 从从Maxwell分布中抽样：分布中抽样：高斯高斯(G
25、auss)分布的随机数生成方法分布的随机数生成方法:222)(2exp)(21xxxxxP生成随机数：生成随机数：i，i=1,2,1246121iiR)R a R)a R)a )R a Ra (123252729SUM=0.0DO 10 I=1,12 SUM=SUM+RANF(DUMMY)10 CONTINUER =(SUM-6.0)/4.0R2=R*RGAUSS=(A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R高斯高斯(Gauss)分布的随机数生成分布的随机数生成(程序程序F24)FACTOR=SQRT(RTEMP)DO 100 I=1,N VX(I)=FACTOR*GAU
26、SS(DUMMY)VY(I)=FACTOR*GAUSS(DUMMY)VZ(I)=FACTOR*GAUSS(DUMMY)100 CONTINUE控制总动量为零：同前面一样处理。控制总动量为零：同前面一样处理。从从Maxwell分布中抽样分布中随机安排初始速度：分布中抽样分布中随机安排初始速度：微正则系综微正则系综MD模拟程序模拟程序F3讲解讲解(LJ,NVE)：量纲变换：量纲变换：SIGMA =(DENS/REAL(N)*(1.0/3.0)RCUT =RCUT*SIGMADT DT*SIGMADENS =DENS/(SIGMA*3)333NVN*模拟盒子的边长为模拟盒子的边长为1VN3*L2/1
27、*)/(mtt 长程校正：长程校正：微正则系综微正则系综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)393*398cclrcrrNU393*32916 cclrclrcrrNVPW算法：算法启动算法：算法启动微正则系综微正则系综MD模拟程序模拟程序F3讲解讲解(LJ,NVE)：CALL F
28、ORCE(-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)算法：差分格式：算法：差分格式：SR2 =SIGSQ/RIJSQVIJ =4.0*(SR12-SR6)WIJ =24.0*(2.0*SR12-SR6)VELIJ=WIJ*DT/RIJSQDVX
29、=VELIJ*RXIJ.VXI =VXI+DVX.VX(J)=VX(J)DVXV =V+VIJW =W+WIJCALL MOVE(DT)114)(6*12*crrrU 121184 )(6*12*rrrWc2*)(rrrWfc*t)(t)t21-(t )t21(t*i*i*iFvvDO 1000 I=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):*t)t21(t)(t )t(t*i*i*ivrr 速度的标定（只用于平衡阶段）速度的标定（只用于平衡阶段）SUMK=0.
30、0DO 200 I=1,N SUMK=SUMKX+VX(I)*2+VY(I)*2+VZ(I)*2200 CONTINUEBEITA=SQRT(3.0*RTEMP/SUMK)DO 300 I=1,N VX(I)=VX(I)*BEITA VY(I)=VY(I)*BEITA VZ(I)=VZ(I)*BEITA300 CONTINUE2/1*2/12*3TTvTreqreq*old*newvv u 具有确定的粒子数具有确定的粒子数N，温度，温度T和体积和体积V速度的直接标度速度的直接标度热浴方法热浴方法(Andersen Thermostat)约束方法约束方法(阻尼力方法阻尼力方法)系统扩展方法系统扩
31、展方法(Extended Systems Method)u 问题的关键：温度的约束问题的关键：温度的约束Nose-Hoover方法方法一、热浴方法一、热浴方法(Andersen Thermostat)u 引入一个与虚拟粒子碰撞的随机力引入一个与虚拟粒子碰撞的随机力u 想象系统浸在热浴当中想象系统浸在热浴当中系统和热浴间的相互作用强度由随机碰撞的频率决定系统和热浴间的相互作用强度由随机碰撞的频率决定碰撞的几率等于碰撞的几率等于Nudt如果一个粒子经历碰撞，它的速度将从约束温度下的如果一个粒子经历碰撞，它的速度将从约束温度下的Maxwell分布中随机抽取分布中随机抽取u 总能量和总动量均不守恒总能
32、量和总动量均不守恒二、约束方法二、约束方法u 是等动能是等动能(Iso-Kinetics)分子动力学方法分子动力学方法u 系统的运动方程为：系统的运动方程为：m/pr pprfp),(u 引入阻尼系数引入阻尼系数 以保证将温度约束在恒定值以保证将温度约束在恒定值u 根据高斯最小约束原理：根据高斯最小约束原理：iiiii2pfp三、三、Nose-Hoover扩展方法扩展方法基本思想：基本思想：设想原系统与一个耦合系统共同组成一个扩展系统，设想原系统与一个耦合系统共同组成一个扩展系统，允许热流在原系统和耦合系统之间交换允许热流在原系统和耦合系统之间交换。Q：等效质量：等效质量S:扩展坐标变量扩展坐
33、标变量：热力学阻尼系数热力学阻尼系数L：扩展系统的自由度扩展系统的自由度u Predictor-corrector algorithm is straightforwardu Verlet algorithm is feasible,but tricky to implement积分方案：积分方案：update of depends on pupdate of p depends on Nos-Hoover方法正确地描述了方法正确地描述了NVT系综中的动系综中的动量和位型，而等动能方法只正确地描述了后者。量和位型，而等动能方法只正确地描述了后者。扩展系统的哈密顿量扩展系统的哈密顿量Hamilt
34、onian守恒：守恒：1.复杂分子体系的势能函数形式：复杂分子体系的势能函数形式：这些方程与描述原子或键各种不同行为的参数就构这些方程与描述原子或键各种不同行为的参数就构成了力场，成了力场，UFF,OPLS,Amber,CVFF,Compass分子模拟方法补充介绍：分子模拟方法补充介绍：kb is the spring constant of the bond.r0 is the bond length at equilibrium.Unique kb and r0 assigned for each bond pair,i.e.C-C,O-Hk is the spring constant
35、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.)kb and k broaden or steepen the slope of the parabola.The larger the value of k,the more energy is required to deform an an
36、gle(or bond)from its equilibrium value.A controls the amplitude of the curven controls its periodicityf 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
37、 on their types(e.g.C-C-C-C,C-O-C-N,H-C-C-H,etc.)A determines the degree the attractivenessB determines the degree of repulsionq is the chargeA determines the degree the attractivenessB determines the degree of repulsionq is the charge2.长程静电力的处理方式：长程静电力的处理方式：Ewald 加和（加和（Ewald Summation）ijjirqqU 有关有关
38、Ewald加和的说明：加和的说明：n Basic form requires an O(N2)calculation efficiency can be introduced to reduce to O(N3/2)good value of a is 5L,but should check for given application can be extended to sum point dipolesn Other methods are in common use reaction field particle-particle/particle mesh(better for larger systems)fast multipole(better for larger systems)3.Verlet近邻表近邻表(Verlet Neighbor lists)：Keep an expanded neighbor list rl rcso that neighbor pairs need not be calculated every timestep rl:chosen such that update lists every 10-20 timesteps Good for N 1000
