分子动力学模拟方法.ppt
现在学习的是第1页，共69页1957年：基于刚球势的分子動力学法（Alder and Wainwright） 1964年：利用Lennard-Jone势函数法对液态氩性质的模拟（Rahman） 1971年：模拟具有分子团簇行为的水的性质（Rahman and Stillinger） 1977年：约束动力学方法（Rychaert, Ciccotti & Berendsen; van Gunsteren） 1980年：恒压条件下的动力学方法（Andersen法、Parrinello-Rahman法） 1983年：非平衡态动力学方法（Gillan and Dixon） 1984年: 恒温条件下的动力学方法(Berendsen et al.) 1984年：恒温条件下的动力学方法（Nos-Hoover法） 1985年：第一原理分子動力学法（Car-Parrinello法） 1991年：巨正则系综的分子动力学方法（Cagin and Pettit）分子动力学简史现在学习的是第2页，共69页粒子的运动取决于经典力学(牛顿定律（F=ma）课程讲解内容：经典分子动力学(Classical Molecular Dynamics)现在学习的是第3页，共69页原理：计算一组分子的相空间轨道，其中每个分子各自服从牛顿运动定律： 11112)(21NiNijijNiiirUmpHiiiiimdtdmvrp 111111)()(NiNijijNiNijijirUrdtdijrFp)0(0itirr)0(0itidtvdr初始条件：现在学习的是第4页，共69页现在学习的是第5页，共69页粒子位置的Taylor展开式：粒子位置：粒子速度： t (t) 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+现在学习的是第6页，共69页n 算法启动t (0) (0) t)(ivrr t (t) t)-(t(t)2 t)(t2iarrriii t2t)-(tt)(t (t)iiirrv现在学习的是第7页，共69页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) = RZNEWI100 CONTINUE现在学习的是第8页，共69页现在学习的是第9页，共69页t (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)(iavv现在学习的是第10页，共69页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)iiivvv现在学习的是第11页，共69页现在学习的是第12页，共69页2iit (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)(tivrrii等价于现在学习的是第13页，共69页n 算法启动2iit (t) 21 t (t) (t) t)(tavrriitt)(t (t) 21 (t) t)(tiiaavvii现在学习的是第14页，共69页VerletLeap-frogVelocity Verlet现在学习的是第15页，共69页)()()()()()(21)()()()(61)(21)()()(232ttttttttttttttttttttttttppppbbbaabavvbavrr现在学习的是第16页，共69页根据新的原子位置rp，可以计算获得校正后的ac(t+t),定义预测误差:)()()(ttttttpcaaa)()()()()()()()()()()()(3210ttctttttctttttctttttctttpcpcpcpcabbaaaavvarr利用此预测误差，对预测出的位置、速度、加速度等量进行校正：现在学习的是第17页，共69页定义一组矢量:333322221061d)(d21d)(dd)(d)(ttttttttttrrrrrrrr)()()()(1000310032101111)()()()(32103210ttttttttttttpppprrrrrrrr)()()()()()(21)()()()(61)(21)()()(232ttttttttttttttttttttttttppppbbbaabavvbavrr现在学习的是第18页，共69页rrrrrrrrr )()()()()()()()(321032103210ccccttttttttttttttttppppccccr的形式：C0， C1， C2， C3的值以及C0，取决于运动方程的阶数。现在学习的是第19页，共69页)(rrf)()(11ttttpcrrrValuesc0c1c2c3c4c535/1211/2 43/813/41/65251/720111/121/31/24 695/288125/2435/725/481/120现在学习的是第20页，共69页)(rrf )()(22ttttpcrrrValuesc0c1c2c3c4c5301141/65/611/3519/1203/411/21/1263/20251/360111/181/61/60现在学习的是第21页，共69页),(rrr f)()(22ttttpcrrrValuesc0c1c2c3c4c5301141/65/611/3519/903/411/21/1263/16251/360111/181/61/60现在学习的是第22页，共69页t室温下， t 1 fs (femtosecond 10-15s)，温度越高，t 应该减小n 太长的时间步长会造成分子间的激烈碰撞，体系数据溢出; n 太短的时间步长会降低模拟过程搜索相空间的能力 现在学习的是第23页，共69页u 它是分子动力学方法的最基本系综u 具有确定的粒子数N，能量E和体积Vu 算法：现在学习的是第24页，共69页微正则系综(NVE)MD模拟算法的流程图：给定每个分子的初始位置ri(0)和速度vi(0)计算每个分子的受力Fi和加速度ai解运动方程并求出每个分子运动一个时间步长后到达的位置所具有的速度统计系统的热力学性质及其它物理量统计性质不变？ 打印结果，结束YesNo移动所有分子到新的位置并具有当前时刻的速度现在学习的是第25页，共69页微正则系综MD模拟程序F3讲解（LJ, NVE）：无因次量：/*rr 3*/*EE /*ff/*3PP /*kTT 2/12*)/(mtt 2/1*)/(vvm现在学习的是第26页，共69页MD模拟中几个热力学量的计算：TNkmEBNiK23v2112i/*kTT 2/1*)/(vvm/*EE *12*i*23v21NTENiK现在学习的是第27页，共69页21114) ru 2)(r(drNrUUUUcrNiNijijclrcc 对于LJ流体：3933138 cclrcrrNU 4)(612rrrUc/*rr 3*/*EE 114)(6*12*rrrUc3*9*113138 ccrrNUlrc现在学习的是第28页，共69页KEUE*KEUE现在学习的是第29页，共69页393232316 cclrcrrPlrccBPVWTkP/*rr 3*/*EE /*3PP /*ff111)(31NiNijijccrWWijijijijijijcrUrWrrfr)(31 )(现在学习的是第30页，共69页*lrccPVWTP3*9*2*1132316 cclrcrrP 111* )(31NiNijijccrWW*)( )(jjjjjjcrUrWiiiiiirrfr现在学习的是第31页，共69页练习：推导LJ流体分子间力的表达式(fx, fy, fz及其对比量)： 4)(612rrrUr)()(rUrf2612 2184 rrrrrrrUrUf/*rr /*ff2*6*12* 121184 rrrrf现在学习的是第32页，共69页* )(fr rWc2*6*12* 121184 rrrrf 121184 )(6*12*rrrWc2* )( rrrWfc现在学习的是第33页，共69页t (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 现在学习的是第34页，共69页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*t )t21(t )(t )t(t*i*i*ivrr2)t21-(t)t21(t (t)*i*i*ivvv现在学习的是第35页，共69页mTkvBx22/12*2/122/122/1333vTvmTkvNmTNkTTreqrefBrefBrefoldnewvv *old*newvv 现在学习的是第36页，共69页微正则系综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 )现在学习的是第37页，共69页每面中心有一格点 现在学习的是第38页，共69页现在学习的是第39页，共69页XLNC=(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 CONTINUE现在学习的是第40页，共69页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页，共69页mTkvB22/1*2/12*2/122/1TTvTvmTkTTreqreqrefBref*old*newvv *2*Tv现在学习的是第42页，共69页*2*Tv随机安排初始速度：现在学习的是第43页，共69页标度初始速度：2/1*2/12*TTvTreqreq*old*newvv 现在学习的是第44页，共69页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) - SUMZ300 CONTINUE控制体系的总动量为零：现在学习的是第45页，共69页x)(xP0 xxdx对于等几率随机试验(Bernoulli试验)，大量的试验结果满足高斯分布 222)(2exp)(21xxxxxP现在学习的是第46页，共69页麦克斯韦速度分布定律:dvekTmNdvvfTkmvB22122mTkvB2 222)(2exp)(21xxxxxP单位体积的分子再每个分量上的速度分布实际上就是高斯分布。现在学习的是第47页，共69页 222)(2exp)(21xxxxxP46121iiR)R a R ) a R ) a )R a Ra ( 123252729现在学习的是第48页，共69页SUM = 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)现在学习的是第49页，共69页从Maxwell分布中抽样分布中随机安排初始速度：现在学习的是第50页，共69页微正则系综MD模拟程序F3讲解(LJ, NVE) ：SIGMA = ( DENS / REAL ( N ) ) * ( 1.0 / 3.0 )RCUT = RCUT * SIGMADT DT * SIGMADENS = DENS / ( SIGMA * 3 )333NVN*模拟盒子的边长为1VN3*L2/1*)/(mtt 现在学习的是第51页，共69页微正则系综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现在学习的是第52页，共69页微正则系综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 )现在学习的是第53页，共69页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) 114)(6*12*crrrU 121184 )(6*12*rrrWc2* )( rrrWfc*t )(t )t21-(t )t21(t*i*i*iFvv现在学习的是第54页，共69页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) * DT1000 CONTINUEMOVE(DT):*t )t21(t )(t )t(t*i*i*ivrr现在学习的是第55页，共69页2/1*2/12*3TTvTreqreq*old*newvv 现在学习的是第56页，共69页u 具有确定的粒子数N，温度T和体积Vu 问题的关键：温度的约束Nose-Hoover方法现在学习的是第57页，共69页一、热浴方法 (Andersen Thermostat)u 引入一个与虚拟粒子碰撞的随机力u 想象系统浸在热浴当中u 总能量和总动量均不守恒现在学习的是第58页，共69页二、约束方法u 是等动能(Iso-Kinetics)分子动力学方法u 系统的运动方程为：m/pr pprfp),(u 引入阻尼系数以保证将温度约束在恒定值u 根据高斯最小约束原理：iiiii2pfp现在学习的是第59页，共69页三、Nose-Hoover扩展方法基本思想：设想原系统与一个耦合系统共同组成一个扩展系统，允许热流在原系统和耦合系统之间交换。Q：等效质量S: 扩展坐标变量： 热力学阻尼系数L： 扩展系统的自由度现在学习的是第60页，共69页u Predictor-corrector algorithm is straightforwardu Verlet algorithm is feasible, but tricky to implement积分方案：update of depends on pupdate of p depends on 现在学习的是第61页，共69页Nos-Hoover方法正确地描述了NVT系综中的动量和位型，而等动能方法只正确地描述了后者。扩展系统的哈密顿量Hamiltonian守恒：现在学习的是第62页，共69页1. 复杂分子体系的势能函数形式：UFF, OPLS, Amber, CVFF, Compass分子模拟方法补充介绍：现在学习的是第63页，共69页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-H现在学习的是第64页，共69页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.)现在学习的是第65页，共69页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 angle (or bond) from its equilibrium value. 现在学习的是第66页，共69页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 on their types (e.g. C-C-C-C, C-O-C-N, H-C-C-H, etc.)现在学习的是第67页，共69页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 charge现在学习的是第68页，共69页2. 长程静电力的处理方式：Ewald 加和（Ewald Summation）ijjirqqU 现在学习的是第69页，共69页
