分子动力学模拟方法

第四章分子动力学模拟方法•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）分子动力学简史粒子的运动取决于经典力学(牛顿定律（F=ma）课程讲解内容：经典分子动力学(ClassicalMolecularDynamics)分子动力学方法基础：原理：计算一组分子的相空间轨道，其中每个分子各自服从牛顿运动定律：11112)(21NiNijijNiiirUmpHiiiiimdtdmvrp111111)()(NiNijijNiNijijirUrdtdijrFp)0(0itirr)0(0itidtvdr初始条件：分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法，可以预测纳米尺度上的材料动力学特性。通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与原子运动路径相关的基本过程。在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。分子动力学方法是确定性方法，一旦初始构型和速度确定了，分子随时间所产生的运动轨迹也就确定了。分子动力学方法特征：分子动力学的算法：有限差分方法一、Verlet算法粒子位置的Taylor展开式：粒子位置：粒子速度：t(t)t)-(t(t)2t)(t2iarrriiit2t)-(tt)(t(t)iiirrv粒子加速度：3i2iit(t)61t(t)21t(t)(t)t)(tbavrrii3i2iit(t)61t(t)21t(t)(t)t)(tbavrrii(t)(t)imiiFa开始运动时需要r(t-Δt)：t(0)(0)t)(ivrr+缺点：Verlet算法处理速度非常笨拙Verlet算法的表述：算法启动①规定初始位置②规定初始速度③扰动初始位置：④计算第n步的力⑤计算第n+1步的位置：计算第n步的速度：重复④至⑥t(0)(0)t)(ivrrt(t)t)-(t(t)2t)(t2iarrriiit2t)-(tt)(t(t)iiirrvVerlet算法程序：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=(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优点：1、精确，误差O(Δ4)2、每次积分只计算一次力3、时间可逆缺点：1、速度有较大误差O(Δ2)2、轨迹与速度无关，无法与热浴耦联Verlet算法的优缺点：二、蛙跳(Leap-frog)算法：半步算法t(t)t)21-(tt)21(tiavviitt)21(t(t)t)(tivrrii2t)21-(tt)21(t(t)iiivvv1.首先利用当前时刻的加速度，计算半个时间步长后的速度：2.计算下一步长时刻的位置：3.计算当前时刻的速度：t-Δt/2tt+Δt/2t+Δtt+3Δt/2t+2Δtvrv开始运动时需要v(-Δt/2)：t/2(0)(0)t/2)(iavvLeap-frog算法的表述：算法启动①规定初始位置②规定初始速度③扰动初始速度：④计算第n步的力⑤计算第n+1/2步的速度：⑥计算第n+1步的位置：⑦计算第n步的速度：⑧重复④至⑦t/2(0)(0)t/2)(iavvt(t)t)21-(tt)21(tiavviitt)21(t(t)t)(tivrrii2t)21-(tt)21(t(t)iiivvvLeap-frog算法的优缺点：优点：1、提高精确度2、轨迹与速度有关，可与热浴耦联缺点：1、速度近似2、比Verlet算子多花时间三、VelocityVerlet算法：2iit(t)21t(t)(t)t)(tavrriitt)](t(t)[21(t)t)(tiiaavviit(t)21(t)t)21(tiavviitt)(t21t)21(tt)(tiavviitt)21(t(t)t)(tivrrii等价于优点：速度计算更加准确VelocityVerlet算法的表述：算法启动①规定初始位置②规定初始速度③计算第n+1步的位置：④计算第n+1步的力⑤计算第n+1步的速度：⑥重复③至⑤2iit(t)21t(t)(t)t)(tavrriitt)](t(t)[21(t)t)(tiiaavviiVerlet三种形式算法的比较：VerletLeap-frogVelocityVerlet四、预测－校正(Predictor-Corrector)格式算法：1.预测(Predictor)阶段：其基本思想是Taylor展开，)()()()()()(21)()()()(61)(21)()()(232ttttttttttttttttttttttttppppbbbaabavvbavrr根据新的原子位置rp，可以计算获得校正后的ac(t+t),定义预测误差:)()()(ttttttpcaaa)()()()()()()()()()()()(3210ttctttttctttttctttttctttpcpcpcpcabbaaaavvarr利用此预测误差，对预测出的位置、速度、加速度等量进行校正：2.校正(Corrector)阶段：预测阶段运动方程的变换：定义一组矢量:333322221061d)(d21d)(dd)(d)(ttttttttttrrrrrrrr)()()()(1000310032101111)()()()(32103210ttttttttttttpppprrrrrrrr)()()()()()(21)()()()(61)(21)()()(232ttttttttttttttttttttttttppppbbbaabavvbavrr校正阶段运动方程的变换：rrrrrrrrr)()()()()()()()(321032103210ccccttttttttttttttttppppccccr的形式：C0，C1，C2，C3的值以及C0，取决于运动方程的阶数。①一阶运动方程：)(rrf)()(11ttttpcrrrValuesc0c1c2c3c4c535/1211/243/813/41/65251/720111/121/31/24695/288125/2435/725/481/120②二阶运动方程之一：)(rrf)()(22ttttpcrrrValuesc0c1c2c3c4c5301141/65/611/3519/1203/411/21/1263/20251/360111/181/61/60②二阶运动方程之二：),(rrrf)()(22ttttpcrrrValuesc0c1c2c3c4c5301141/65/611/3519/903/411/21/1263/16251/360111/181/61/60五、积分时间步长t的选择：室温下，∆t≈1fs(femtosecond10-15s)，温度越高，∆t应该减小太长的时间步长会造成分子间的激烈碰撞，体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力微正则系综分子动力学(NVEMD)它是分子动力学方法的最基本系综具有确定的粒子数N，能量E和体积V算法：①规定初始位置和初始速度②对运动方程积分若干步③计算势能和动能④若能量不等于所需要的值，对速度进行标度⑤重复②至④，直到系统平衡微正则系综(NVE)MD模拟算法的流程图：给定每个分子的初始位置ri(0)和速度vi(0)计算每个分子的受力Fi和加速度ai解运动方程并求出每个分子运动一个时间步长后到达的位置所具有的速度统计系统的热力学性质及其它物理量统计性质不变？打印结果，结束YesNo移动所有分子到新的位置并具有当前时刻的速度微正则系综MD模拟程序F3讲解（LJ,NVE）：无因次量：/*rr3*/*EE/*ff/*3PP/*kTT2/12*)/(mtt2/1*)/(vvmMD模拟中几个热力学量的计算：对于由N个单原子组成的系统：①动能和温度：TNkmEBNiK23v2112i采用对比量：/*kTT2/1*)/(vvm/*EE*12*i*23v21NTENiK2111π4)ru2)(r(drNrUUUUcrNiNijijclrcc对于LJ流体：39331επ38cclrcrrNU②势能：4)(612rrrUc/*rr3*/*EE采用对比量：114)(6*12***rrrUc3*9***1131π38ccrrNUlrc③内能：KEUE内能由势能和动能组成：采用对比量：***KEUE④压力：393232επ316cclrcrrPlrccBPVWTkP/*rr3*/*EE采用对比量：/*3PP/*ff111)(31NiNijijccrWWijijijijijijcrUrWrrfr)(31)(******lrccPVWTP3*9*2**1132π316cclrcrrP111***)(31NiNijijccrWW********)()(jjjjjjcrUrWiiiiiirrfr练习：推导LJ流体分子间力的表达式(fx,fy,fz及其对比量)：4)(612rrrU势能函数形式：力：r)()(rUr