分子动力学模拟方法优秀课件.ppt

《分子动力学模拟方法优秀课件.ppt》由会员分享，可在线阅读，更多相关《分子动力学模拟方法优秀课件.ppt（72页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、分子分子分子分子动动力学模力学模力学模力学模拟拟方法方法方法方法第1页，本讲稿共72页1957年年：基于刚球势的分子動力学法（：基于刚球势的分子動力学法（Alder and Wainwright）1964年年：利用：利用Lennard-Jone势函数法对液态势函数法对液态氩氩性质的模拟（性质的模拟（Rahman）1971年年：模拟具有分子：模拟具有分子团簇团簇行为的行为的水水的性质（的性质（Rahman and Stillinger）1977年年：约束动力学方法（：约束动力学方法（Rychaert,Ciccotti&Berendsen;van Gunsteren）1980年年：恒压条件下的动力

2、学方法（：恒压条件下的动力学方法（Andersen法、法、Parrinello-Rahman法）法）1983年年：非平衡态动力学方法（：非平衡态动力学方法（Gillan and Dixon）1984年年:恒温条件下的动力学方法恒温条件下的动力学方法(Berendsen et al.)1984年年：恒温条件下的动力学方法（：恒温条件下的动力学方法（Nos-Hoover法）法）1985年年：第一原理分子動力学法（：第一原理分子動力学法（Car-Parrinello法）法）1991年年：巨正则系综巨正则系综的分子动力学方法（的分子动力学方法（Cagin and Pettit）分子动力学简史第2页，本

3、讲稿共72页粒子的运动取决于经典力学(牛顿定律（F=ma）课程讲解内容：经典分子动力学(ClassicalMolecularDynamics)第3页，本讲稿共72页分子动力学方法基础：分子动力学方法基础：分子动力学方法基础：分子动力学方法基础：原理：计算一组分子的相空间轨道，其中每个分子各自服从牛顿运动定律：初始条件：第4页，本讲稿共72页n n分分分分子子子子动动动动力力力力学学学学是是是是在在在在原原原原子子子子、分分分分子子子子水水水水平平平平上上上上求求求求解解解解多多多多体体体体问问问问题题题题的的的的重重重重要要要要的的的的计计计计算算算算机机机机模模模模拟拟拟拟方方方方法法法法，

4、可以预测纳米尺度上的材料动力学特性。可以预测纳米尺度上的材料动力学特性。可以预测纳米尺度上的材料动力学特性。可以预测纳米尺度上的材料动力学特性。n n通通通通过过过过求求求求解解解解所所所所有有有有粒粒粒粒子子子子的的的的运运运运动动动动方方方方程程程程，分分分分子子子子动动动动力力力力学学学学方方方方法法法法可可可可以以以以用用用用于于于于模模模模拟拟拟拟与与与与原子运动路径相关的基本过程。原子运动路径相关的基本过程。原子运动路径相关的基本过程。原子运动路径相关的基本过程。n n在在在在分分分分子子子子动动动动力力力力学学学学中中中中，粒粒粒粒子子子子的的的的运运运运动动动动行行行行为为为为

5、是是是是通通通通过过过过经经经经典典典典的的的的NewtonNewtonNewtonNewton运运运运动动动动方方方方程程程程所所所所描描描描述。述。述。述。n n分分分分子子子子动动动动力力力力学学学学方方方方法法法法是是是是确确确确定定定定性性性性方方方方法法法法，一一一一旦旦旦旦初初初初始始始始构构构构型型型型和和和和速速速速度度度度确确确确定定定定了了了了，分分分分子子子子随随随随时间所产生的运动轨迹也就确定了。时间所产生的运动轨迹也就确定了。时间所产生的运动轨迹也就确定了。时间所产生的运动轨迹也就确定了。分子动力学方法特征：分子动力学方法特征：分子动力学方法特征：分子动力学方法特征

6、：第5页，本讲稿共72页分子动力学的算法：有限差分方法分子动力学的算法：有限差分方法分子动力学的算法：有限差分方法分子动力学的算法：有限差分方法一、一、一、一、VerletVerlet算法算法算法算法粒子位置的Taylor展开式：粒子位置：粒子位置：粒子速度：粒子速度：粒子加速度：粒子加速度：开始运动时需要开始运动时需要r(t-t)：+缺点：缺点：缺点：缺点：VerletVerlet算法处理速度非常笨拙算法处理速度非常笨拙算法处理速度非常笨拙算法处理速度非常笨拙第6页，本讲稿共72页VerletVerlet算法的表述：算法的表述：算法的表述：算法的表述：n算法启动规定初始位置规定初始位置规定初

7、始位置规定初始位置规定初始速度规定初始速度规定初始速度规定初始速度扰动初始位置：扰动初始位置：扰动初始位置：扰动初始位置：计算第计算第计算第计算第n n步的力步的力步的力步的力计算第计算第计算第计算第n+1n+1步的位置：步的位置：步的位置：步的位置：n n计算第计算第计算第计算第n n步的速度：步的速度：步的速度：步的速度：n n重复重复重复重复至至至至第7页，本讲稿共72页VerletVerlet算法程序：算法程序：算法程序：算法程序：Do 100 I=1,N RXNEWI=2.0*RX(I)RXOLD(I)+DTSQ*AX(I)RYNEWI=2.0*RY(I)RYOLD(I)+DTSQ*

8、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)=RZNEWI100 CONTINUE第8页，本讲稿共72页n n优点：优点：优点：优点：1 1、精确，误差、精确，误差、精确，误差、精确，误差O(4)O(4)2 2、每次积分只计算一次力、每次积分只计算一次力、每次积分只计

9、算一次力、每次积分只计算一次力3 3、时间可逆、时间可逆、时间可逆、时间可逆n n缺点：缺点：缺点：缺点：1 1、速度有较大误差、速度有较大误差、速度有较大误差、速度有较大误差O(2)O(2)2 2、轨迹与速度无关，无法与热浴耦联、轨迹与速度无关，无法与热浴耦联、轨迹与速度无关，无法与热浴耦联、轨迹与速度无关，无法与热浴耦联VerletVerlet算法的优缺点：算法的优缺点：算法的优缺点：算法的优缺点：第9页，本讲稿共72页二、蛙跳二、蛙跳二、蛙跳二、蛙跳(Leap-frog)(Leap-frog)算法：半步算法算法：半步算法算法：半步算法算法：半步算法1.首先利用当前时刻的加速度，计算半个时

10、间步长后的速度：2.计算下一步长时刻的位置：3.计算当前时刻的速度：t-t/2tt+t/2t+tt+3t/2t+2tvrv开始运动时需要开始运动时需要v(-t/2)：第10页，本讲稿共72页Leap-frogLeap-frog算法的表述：算法的表述：算法的表述：算法的表述：n算法启动规定初始位置规定初始位置规定初始位置规定初始位置规定初始速度规定初始速度规定初始速度规定初始速度扰动初始速度：扰动初始速度：扰动初始速度：扰动初始速度：计算第计算第计算第计算第n n步的力步的力步的力步的力计算第计算第计算第计算第n+1/2n+1/2步的速度：步的速度：步的速度：步的速度：计算第计算第计算第计算第n

11、+1n+1步的位置：步的位置：步的位置：步的位置：计算第计算第计算第计算第n n步的速度：步的速度：步的速度：步的速度：重复重复重复重复至至至至第11页，本讲稿共72页Leap-frogLeap-frog算法的优缺点：算法的优缺点：算法的优缺点：算法的优缺点：n n优点：优点：优点：优点：1 1、提高精确度、提高精确度、提高精确度、提高精确度2 2、轨迹与速度有关，可与热浴耦联、轨迹与速度有关，可与热浴耦联、轨迹与速度有关，可与热浴耦联、轨迹与速度有关，可与热浴耦联n n缺点：缺点：缺点：缺点：1 1、速度近似、速度近似、速度近似、速度近似2 2、比、比、比、比VerletVerlet算子多花

12、时间算子多花时间算子多花时间算子多花时间第12页，本讲稿共72页三、三、三、三、Velocity VerletVelocity Verlet算法：算法：算法：算法：等价于优点：速度计算更加准确优点：速度计算更加准确优点：速度计算更加准确优点：速度计算更加准确第13页，本讲稿共72页Velocity VerletVelocity Verlet算法的表述：算法的表述：算法的表述：算法的表述：n算法启动规定初始位置规定初始位置规定初始位置规定初始位置规定初始速度规定初始速度规定初始速度规定初始速度计算第计算第计算第计算第n+1n+1步的位置：步的位置：步的位置：步的位置：计算第计算第计算第计算第n+

13、1n+1步的力步的力步的力步的力计算第计算第计算第计算第n+1n+1步的速度：步的速度：步的速度：步的速度：重复重复重复重复至至至至第14页，本讲稿共72页VerletVerlet三种形式算法的比较：三种形式算法的比较：三种形式算法的比较：三种形式算法的比较：VerletLeap-frogVelocity Verlet第15页，本讲稿共72页四、预测校正四、预测校正四、预测校正四、预测校正(Predictor-Corrector)(Predictor-Corrector)格式算法：格式算法：格式算法：格式算法：1.1.预测预测预测预测(Predictor)(Predictor)阶段：其基本思想

14、是阶段：其基本思想是阶段：其基本思想是阶段：其基本思想是TaylorTaylor展开，展开，展开，展开，第16页，本讲稿共72页根据新的原子位置rp，可以计算获得校正后的ac(t+t),定义预测误差:利用此预测误差，对预测出的位置、速度、加速度等量进行校正：2.2.校正校正校正校正(Corrector)(Corrector)阶段：阶段：阶段：阶段：第17页，本讲稿共72页n n预测阶段运动方程的变换：预测阶段运动方程的变换：预测阶段运动方程的变换：预测阶段运动方程的变换：定义一组矢量:第18页，本讲稿共72页n n校正阶段运动方程的变换：校正阶段运动方程的变换：校正阶段运动方程的变换：校正阶段

15、运动方程的变换：的形式：C0，C1，C2，C3的值以及C0，取决于运动方程的阶数。第19页，本讲稿共72页 一阶运动方程：一阶运动方程：一阶运动方程：一阶运动方程：Valuesc0c1c2c3c4c535/1211/243/813/41/65251/720111/121/31/24695/288125/2435/725/481/120第20页，本讲稿共72页 二阶运动方程之一：二阶运动方程之一：二阶运动方程之一：二阶运动方程之一：Valuesc0c1c2c3c4c5301141/65/611/3519/1203/411/21/1263/20251/360111/181/61/60第21页，本讲

16、稿共72页 二阶运动方程之二：二阶运动方程之二：二阶运动方程之二：二阶运动方程之二：Valuesc0c1c2c3c4c5301141/65/611/3519/903/411/21/1263/16251/360111/181/61/60第22页，本讲稿共72页五、积分时间步长五、积分时间步长五、积分时间步长五、积分时间步长t的选择：的选择：的选择：的选择：室温下，室温下，t 1 fs(femtosecond 10-15s)，温度越高，温度越高，t 应该减小应该减小n 太长的时间步长会造成分子间的激烈碰撞，体系数据溢出;n 太短的时间步长会降低模拟过程搜索相空间的能力 第23页，本讲稿共72页微正

17、则系综分子动力学微正则系综分子动力学微正则系综分子动力学微正则系综分子动力学(NVE MD)(NVE MD)u它是分子动力学方法的最基本系综u具有确定的粒子数N，能量E和体积Vu算法：规定初始位置和初始速度规定初始位置和初始速度规定初始位置和初始速度规定初始位置和初始速度对运动方程积分若干步对运动方程积分若干步对运动方程积分若干步对运动方程积分若干步计算势能和动能计算势能和动能计算势能和动能计算势能和动能若能量不等于所需要的值，对速度进行标度若能量不等于所需要的值，对速度进行标度若能量不等于所需要的值，对速度进行标度若能量不等于所需要的值，对速度进行标度重复重复重复重复至至至至，直到系统平衡，

18、直到系统平衡，直到系统平衡，直到系统平衡第24页，本讲稿共72页微正则系综(NVE)MD模拟算法的流程图：给定每个分子的初始位置ri(0)和速度vi(0)计算每个分子的受力Fi和加速度ai解运动方程并求出每个分子运动一个时间步长后到达的位置所具有的速度统计系统的热力学性质及其它物理量统计性质不变？打印结果，结束YesNo移动所有分子到新的位置并具有当前时刻的速度第25页，本讲稿共72页微正则系综MD模拟程序F3讲解（LJ,NVE）：无因次量：第26页，本讲稿共72页MD模拟中几个热力学量的计算：对于由对于由对于由对于由N N个单原子组成的系统：个单原子组成的系统：个单原子组成的系统：个单原子组

19、成的系统：动能和温度：动能和温度：动能和温度：动能和温度：采用对比量：采用对比量：采用对比量：采用对比量：第27页，本讲稿共72页对于LJ流体：势能：势能：势能：势能：采用对比量：采用对比量：采用对比量：采用对比量：第28页，本讲稿共72页内能：内能：内能：内能：内能由势能和动能组成：内能由势能和动能组成：内能由势能和动能组成：内能由势能和动能组成：采用对比量：采用对比量：采用对比量：采用对比量：第29页，本讲稿共72页压力：压力：压力：压力：采用对比量：采用对比量：采用对比量：采用对比量：第30页，本讲稿共72页第31页，本讲稿共72页练习：推导LJ流体分子间力的表达式(fx,fy,fz及其

20、对比量)：势能函数形式：势能函数形式：势能函数形式：势能函数形式：力：力：力：力：采用对比量：采用对比量：采用对比量：采用对比量：=x,y,z=x,y,z第32页，本讲稿共72页LJLJ分子间的维里项：分子间的维里项：分子间的维里项：分子间的维里项：=x,y,z=x,y,z第33页，本讲稿共72页采用对比量的运动方程形式：采用对比量的运动方程形式：采用对比量的运动方程形式：采用对比量的运动方程形式：(以蛙跳以蛙跳以蛙跳以蛙跳(Leap-frog)(Leap-frog)算法为例算法为例算法为例算法为例)采用对比量：采用对比量：采用对比量：采用对比量：第34页，本讲稿共72页最终得到：最终得到：最

21、终得到：最终得到：同理得到：同理得到：同理得到：同理得到：第35页，本讲稿共72页速度的标度速度的标度速度的标度速度的标度(Velocity Scaling)(Velocity Scaling)：根据能量均分原理，可知：根据能量均分原理，可知：根据能量均分原理，可知：根据能量均分原理，可知：标度因子标度因子标度因子标度因子:对比量对比量对比量对比量速度标度速度标度速度标度速度标度:或或或或第36页，本讲稿共72页微正则系综MD模拟程序F3讲解（LJ,NVE）：初始化：初始化：初始化：初始化：READ(*,(A)TITLE !运行作业题目运行作业题目READ(*,*)NSTEP !运行步数运行步

22、数READ(*,*)IPRINT !打印步数打印步数READ(*,(A)CNFILE !位型文件位型文件READ(*,*)DENS !对比密度对比密度READ(*,*)RTEMP !对比温度对比温度 READ(*,*)RCUT !对比截断半径对比截断半径READ(*,*)DT !对比时间步长对比时间步长CALL READCN(CNFILE)第37页，本讲稿共72页初始位型：初始位型：初始位型：初始位型：n n面心立方面心立方面心立方面心立方 (face-centered cubic,FCC)face-centered cubic,FCC)：每面中心有一格点 第38页，本讲稿共72页n n体心立

23、方体心立方体心立方体心立方 (body-centered cubic,BCC)body-centered cubic,BCC)：n n简单立方简单立方简单立方简单立方 (simple cubic,SC)simple cubic,SC)：第39页，本讲稿共72页XL初始位型：初始位型：初始位型：初始位型：面心立方面心立方面心立方面心立方 (FCC)(FCC)(程序程序程序程序F23)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 1

24、0 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 CONTINUE第40页，本讲稿共72页DO 100 I=1,N RX(I)=RX(I)-0.5 RY(I)=RY(I)-0.5 RZ(I)=RZ(I)-0.5100 CONTINUE将模拟盒子的中心移到原点：将模拟盒子的中心移到原点：将模拟盒子的中心移到原点：将模拟盒子的中心移到原点：第41页，本讲稿共72页初始速度：初

25、始速度：初始速度：初始速度：n n简单的选择：简单的选择：简单的选择：简单的选择：V V random(-0.5,0.5)random(-0.5,0.5)=x,y,z=x,y,z标度因子标度因子标度因子标度因子:速度标度速度标度速度标度速度标度:第42页，本讲稿共72页FACTOR=SQRT(RTEMP)FACTOR=SQRT(RTEMP)DO 100 I=1,NDO 100 I=1,N VX(I)=FACTOR VX(I)=FACTOR *(RANF(DUMMY)-0.5)*(RANF(DUMMY)-0.5)VY(I)=FACTOR VY(I)=FACTOR *(RANF(DUMMY)-0.

26、5)*(RANF(DUMMY)-0.5)VZ(I)=FACTOR VZ(I)=FACTOR *(RANF(DUMMY)-0.5)*(RANF(DUMMY)-0.5)100100 CONTINUE CONTINUE随机安排初始速度：第43页，本讲稿共72页标度初始速度：SUMKX=0.0SUMKX=0.0SUMKY=0.0SUMKY=0.0SUMKZ=0.0SUMKZ=0.0DO 200 I=1,NDO 200 I=1,N SUMKX=SUMKX+VX(I)*2 SUMKX=SUMKX+VX(I)*2 SUMKY=SUMKY+VY(I)*2 SUMKY=SUMKY+VY(I)*2 SUMKZ=S

27、UMKZ+VZ(I)*2 SUMKZ=SUMKZ+VZ(I)*2200 CONTINUE200 CONTINUEBEITAX=SQRT(RTEMP/SUMKX)BEITAX=SQRT(RTEMP/SUMKX)BEITAY=SQRT(RTEMP/SUMKY)BEITAY=SQRT(RTEMP/SUMKY)BEITAZBEITAZ =SQRT(RTEMP/SUMKZ)=SQRT(RTEMP/SUMKZ)DO 300 I=1,NDO 300 I=1,N VX(I)=VX(I)*BEITAX VX(I)=VX(I)*BEITAX VY(I)=VY(I)*BEITAY VY(I)=VY(I)*BEITA

28、Y VZ(I)=VZ(I)*BEITAZ VZ(I)=VZ(I)*BEITAZ300 CONTINUE300 CONTINUE标度因子标度因子标度因子标度因子:第44页，本讲稿共72页SUMX=0.0SUMX=0.0SUMY=0.0SUMY=0.0SUMZ=0.0SUMZ=0.0DO 200 I=1,NDO 200 I=1,N SUMX=SUMX+VX(I)SUMX=SUMX+VX(I)SUMY=SUMY+VY(I)SUMY=SUMY+VY(I)SUMZ=SUMZ+VZ(I)SUMZ=SUMZ+VZ(I)200200 CONTINUE CONTINUESUMX=SUMX/REAL(N)SUMY

29、=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控制体系的总动量为零：第45页，本讲稿共72页n n从从从从MaxwellMaxwell分布中抽样：分布中抽样：分布中抽样：分布中抽样：xxdx高斯高斯高斯高斯(Gauss)(Gauss)分布分布分布分布:对于等几率随机试验(Bernoulli试验)，大量的试验结果满足高斯分布第46页，本讲稿共72页麦克斯韦速度分布定律:由于：由于：由于：由于：=x,y,z=x,y,z单位体积的分子再每

30、个分量上的速度分布实际上就是高斯分布。第47页，本讲稿共72页n n从从从从MaxwellMaxwell分布中抽样：分布中抽样：分布中抽样：分布中抽样：高斯高斯高斯高斯(Gauss)(Gauss)分布的随机数生成方法分布的随机数生成方法分布的随机数生成方法分布的随机数生成方法:生成随机数：生成随机数：生成随机数：生成随机数：i i，i=1,2,12i=1,2,12第48页，本讲稿共72页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)*

31、R高斯(Gauss)分布的随机数生成(程序F24)第49页，本讲稿共72页FACTOR=SQRT(RTEMP)FACTOR=SQRT(RTEMP)DO 100 I=1,NDO 100 I=1,N VX(I)=FACTOR VX(I)=FACTOR *GAUSS(DUMMY)*GAUSS(DUMMY)VY(I)=FACTOR VY(I)=FACTOR *GAUSS(DUMMY)*GAUSS(DUMMY)VZ(I)=FACTOR VZ(I)=FACTOR *GAUSS(DUMMY)*GAUSS(DUMMY)100100 CONTINUE CONTINUE控制总动量为零：同前面一样处理。控制总动量为

32、零：同前面一样处理。控制总动量为零：同前面一样处理。控制总动量为零：同前面一样处理。从Maxwell分布中抽样分布中随机安排初始速度：第50页，本讲稿共72页微正则系综MD模拟程序F3讲解(LJ,NVE)：量纲变换：量纲变换：量纲变换：量纲变换：SIGMA =(DENS/REAL(N)*(1.0/3.0)RCUT =RCUT*SIGMADT DT*SIGMADENS =DENS/(SIGMA*3)模拟盒子的边长为1第51页，本讲稿共72页长程校正：长程校正：长程校正：长程校正：微正则系综MD模拟程序F3讲解(LJ,NVE)：SR3=(SIGMA/RCUT)*3SR9=SR3*3SIGCUB=S

33、IGMA*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)第52页，本讲稿共72页算法：算法启动算法：算法启动算法：算法启动算法：算法启动微正则系综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

34、)CALL KINET(OLDK)CALL MOVE (DT)CALL FORCE(DT,SIGMA,RCUT,NEWV,NEWVC,NEWW)CALL KINET(NEWK)第53页，本讲稿共72页算法：差分格式：算法：差分格式：算法：差分格式：算法：差分格式：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+WIJCALL MOVE(DT)第54页，本讲稿共72页DO1000I

35、=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):第55页，本讲稿共72页速度的标定（只用于平衡阶段）速度的标定（只用于平衡阶段）速度的标定（只用于平衡阶段）速度的标定（只用于平衡阶段）SUMK=0.0SUMK=0.0DO 200 I=1,NDO 200 I=1,N SUMK=SUMKX+VX(I)*2+VY(I)*2+VZ(I)*2 SUMK=SUMKX+VX(I)*2+VY(I)*2+VZ(I)*2200 CONTINUE200 CONTINUEBEITA=SQR

36、T(3.0*RTEMP/SUMK)BEITA=SQRT(3.0*RTEMP/SUMK)DO 300 I=1,NDO 300 I=1,N VX(I)=VX(I)*BEITA VX(I)=VX(I)*BEITA VY(I)=VY(I)*BEITA VY(I)=VY(I)*BEITA VZ(I)=VZ(I)*BEITA VZ(I)=VZ(I)*BEITA300 CONTINUE300 CONTINUE第56页，本讲稿共72页正则系综分子动力学正则系综分子动力学正则系综分子动力学正则系综分子动力学(NVT MD)(NVT MD)u具有确定的粒子数N，温度T和体积V速度的直接标度速度的直接标度速度的直接

37、标度速度的直接标度热浴方法热浴方法热浴方法热浴方法 (Andersen Thermostat)(Andersen Thermostat)约束方法约束方法约束方法约束方法(阻尼力方法阻尼力方法阻尼力方法阻尼力方法)系统扩展方法系统扩展方法系统扩展方法系统扩展方法(Extended Systems Method)(Extended Systems Method)u问题的关键：温度的约束Nose-Hoover方法第57页，本讲稿共72页一、热浴方法(Andersen Thermostat)u引入一个与虚拟粒子碰撞的随机力u想象系统浸在热浴当中系统和热浴间的相互作用强度由随机碰撞的频率决定系统和热浴间

38、的相互作用强度由随机碰撞的频率决定系统和热浴间的相互作用强度由随机碰撞的频率决定系统和热浴间的相互作用强度由随机碰撞的频率决定碰撞的几率等于碰撞的几率等于碰撞的几率等于碰撞的几率等于NuNu dtdt如如如如果果果果一一一一个个个个粒粒粒粒子子子子经经经经历历历历碰碰碰碰撞撞撞撞，它它它它的的的的速速速速度度度度将将将将从从从从约约约约束束束束温温温温度度度度下下下下的的的的MaxwellMaxwell分布中随机抽取分布中随机抽取分布中随机抽取分布中随机抽取u总能量和总动量均不守恒第58页，本讲稿共72页二、约束方法u是等动能(Iso-Kinetics)分子动力学方法u 系统的运动方程为：u引

39、入阻尼系数以保证将温度约束在恒定值u根据高斯最小约束原理：第59页，本讲稿共72页三、Nose-Hoover扩展方法基本思想：设想原系统与一个耦合系统共同组成一个扩展系统，允许热流在原系统和耦合系统之间交换。Q：等效质量S:扩展坐标变量：热力学阻尼系数L：扩展系统的自由度第60页，本讲稿共72页u Predictor-corrector algorithm is straightforwardu Verlet algorithm is feasible,but tricky to implement积分方案：update of x depends on pupdate of p depends

40、 on x第61页，本讲稿共72页Nos-Hoover方方法法正正确确地地描描述述了了NVT系系综综中中的的动动量量和和位型，而等动能方法只正确地描述了后者。位型，而等动能方法只正确地描述了后者。扩展系统的哈密顿量Hamiltonian守恒：第62页，本讲稿共72页1.复杂分子体系的势能函数形式：Potential Energy=Stretching Energy+Potential Energy=Stretching Energy+Bending Energy+Bending Energy+Torsion Energy+Torsion Energy+Non-Bonded Interactio

41、n Energy Non-Bonded Interaction Energy这些方程与描述原子或键各种不同行为的参数就构这些方程与描述原子或键各种不同行为的参数就构这些方程与描述原子或键各种不同行为的参数就构这些方程与描述原子或键各种不同行为的参数就构成了力场，成了力场，成了力场，成了力场，force-fieldforce-field.UFF,OPLS,Amber,CVFF,Compass分子模拟方法补充介绍：第63页，本讲稿共72页Bond Stretching EnergyBond Stretching Energykb is the spring constant of the bond

42、.r0 is the bond length at equilibrium.Unique kb and r0 assigned for each bond pair,i.e.C-C,O-H第64页，本讲稿共72页Bending EnergyBending Energyk 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 the

43、ir types(e.g.C-C-C,C-O-C,C-C-H,etc.)第65页，本讲稿共72页The“Hookeian”potentialThe“Hookeian”potentialkb 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.第66页，本讲稿共72页Torsion EnergyTorsion EnergyA cont

44、rols 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 on their types(e.g.C-C-C-C,C-O-C-N,H-C-C-H,e

45、tc.)第67页，本讲稿共72页Non-bonded EnergyNon-bonded EnergyA 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 charge第68页，本讲稿共72页2.长程静电力的处理方式：Ewald 加和（Ewald Summation）第69页，本讲稿共72页有关

46、Ewald加和的说明：n nBasic form requires an O(NBasic form requires an O(N2 2)calculation)calculation efficiency can be introduced to reduce to O(Nefficiency can be introduced to reduce to O(N3/23/2)good value of a is 5L,but should check for given applicationgood value of a is 5L,but should check for given

47、application can be extended to sum point dipolescan be extended to sum point dipolesn nOther methods are in common useOther methods are in common use reaction fieldreaction field particle-particle/particle mesh(better for larger systems)particle-particle/particle mesh(better for larger systems)fast

48、multipole(better for larger systems)fast multipole(better for larger systems)第70页，本讲稿共72页3.Verlet近邻表(VerletNeighborlists)：n nKeep an expanded neighbor list Keep an expanded neighbor list r rl l r rc cso that neighbor pairs need not be calculated every timestep rl:chosen such that update lists every 10-20 timesteps Good for N 1000第72页，本讲稿共72页