分子动力学模拟方法

《分子动力学模拟方法》由会员分享，可在线阅读，更多相关《分子动力学模拟方法（72页珍藏版）》请在金锄头文库上搜索。

1、第四章 分子动力学模拟方法1957年：基于刚球势的分子動力学法（Alder and Wainwright） 1964年：利用Lennard-Jone势函数法对液态氩性质的模拟（Rahman） 1971年：模拟具有分子团簇行为的水的性质（Rahman and Stillinger） 1977年：约束动力学方法（Rychaert, Ciccotti van Gunsteren） 1980年：恒压条件下的动力学方法（Andersen法、Parrinello-Rahman法） 1983年：非平衡态动力学方法（Gillan and Dixon） 1984年: 恒温条件下的动力学方法(Berendsen

2、et al.) 1984年：恒温条件下的动力学方法（Nos-Hoover法） 1985年：第一原理分子動力学法（Car-Parrinello法） 1991年：巨正则系综的分子动力学方法（Cagin and Pettit）分子动力学简史粒子的运动取决于经典力学( 牛顿定律（F=ma）课程讲解内容：经典分子动力学 (Classical Molecular Dynamics)分子动力学方法基础：原理： 计算一组分子的相空间轨道，其中每个分子各自服从 牛顿运动定律：初始条件：n 分子动力学是在原子、分子水平上求解多体问题的重要的计 算机模拟方法，可以预测纳米尺度上的材料动力学特性。n 通过求解所有粒子

3、的运动方程，分子动力学方法可以用于模 拟与原子运动路径相关的基本过程。n 在分子动力学中，粒子的运动行为是通过经典的Newton运动 方程所描述。n 分子动力学方法是确定性方法，一旦初始构型和速度确定了 ，分子随时间所产生的运动轨迹也就确定了。分子动力学方法特征：分子动力学的算法：有限差分方法一、Verlet算法粒子位置的Taylor展开式：粒子位置：粒子速度：粒子加速度：开始运动时需要r(t-t)：+缺点：Verlet算法处理速度非常笨拙Verlet算法的表述：n 算法启动 规定初始位置 规定初始速度 扰动初始位置： 计算第n步的力 计算第n+1步的位置：u 计算第n步的速度：u 重复至Ve

4、rlet算法程序：Do 100 I = 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(

5、I) = RXNEWIRY(I) = RYNEWIRZ(I) = RZNEWI 100 CONTINUEn 优点： 1、精确，误差O(4) 2、每次积分只计算一次力 3、时间可逆n 缺点： 1、速度有较大误差O(2) 2、轨迹与速度无关，无法与热浴耦联Verlet算法的优缺点：二、蛙跳(Leap-frog)算法：半步算法1. 首先利用当前时刻的加速度，计算半个时间步长后的速度：2. 计算下一步长时刻的位置：3. 计算当前时刻的速度：t-t/2tt+t/2t+tt+3t/2t+2tvrv开始运动时需要v(-t/2)：Leap-frog算法的表述：n 算法启动 规定初始位置 规定初始速度 扰动初始

6、速度： 计算第n步的力 计算第n+1/2步的速度： 计算第n+1步的位置： 计算第n步的速度： 重复至Leap-frog算法的优缺点：n 优点： 1、提高精确度 2、轨迹与速度有关，可与热浴耦联n 缺点： 1、速度近似 2、比Verlet算子多花时间三、Velocity Verlet算法：等价于优点：速度计算更加准确Velocity Verlet算法的表述：n 算法启动 规定初始位置 规定初始速度 计算第n+1步的位置： 计算第n+1步的力 计算第n+1步的速度： 重复至Verlet三种形式算法的比较：VerletLeap- frogVelocity Verlet四、预测校正(Predicto

7、r-Corrector)格式算法：1. 预测(Predictor)阶段：其基本思想是Taylor展开，根据新的原子位置rp，可以计算获得校正后的ac(t+t),定义预测误差:利用此预测误差，对预测出的位置、速度、加速度等量进行校正：2. 校正(Corrector)阶段：u预测阶段运动方程的变换：定义一组矢量:u校正阶段运动方程的变换：的形式：C0， C1， C2， C3的值以及C0，取决于运动方程的阶数。 一阶运动方程：Valuesc0c1c2c3c4c535/1211/2 43/813/41/6 5251/7 20111/121/31/24 695/28 8125/24 35/725/481

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

9、VE MD)u 它是分子动力学方法的最基本系综u 具有确定的粒子数N，能量E和体积Vu 算法： 规定初始位置和初始速度 对运动方程积分若干步 计算势能和动能 若能量不等于所需要的值，对速度进行标度 重复至，直到系统平衡微正则系综(NVE)MD模拟算法的流程图：给定每个分子的初始位置ri(0)和速度vi(0)计算每个分子的受力Fi和加速度ai解运动方程并求出每个分子运动一个时间步 长后到达的位置所具有的速度统计系统的热力学性质及其它物理量统计性质不变？ 打印结果，结束YesNo移动所有分子到新的位置并具有当前时刻的 速度微正则系综MD模拟程序F3讲解（LJ, NVE）：无因次量：MD模拟中几个热

10、力学量的计算：对于由N个单原子组成的系统：动能和温度：采用对比量：对于LJ流体：势能：采用对比量：内能：内能由势能和动能组成：采用对比量：压力：采用对比量：练习：推导LJ流体分子间力的表达式(fx, fy, fz及其对比量)：势能函数形式：力：采用对比量：=x, y, zLJ分子间的维里项：=x, y, z采用对比量的运动方程形式： (以蛙跳(Leap-frog)算法为例)采用对比量：最终得到：同理得到：速度的标度(Velocity Scaling)：根据能量均分原理，可知：标度因子:对比量速度标度:或微正则系综MD模拟程序F3讲解（LJ, NVE）： 初始化：READ (*,(A) TITL

11、E ! 运行作业题目READ (*,*) NSTEP ! 运行步数READ (*,*) IPRINT ! 打印步数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)

12、：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 = 0 DO 10 I = 1, NC DO 10 J = 1, NC DO 10 K = 1, NCDO 11 IJ = 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) 11 CON

13、TINUEM = M + 4 10 CONTINUEDO 100 I = 1, NRX(I) = RX(I) - 0.5RY(I) = RY(I) - 0.5RZ(I) = RZ(I) - 0.5100 CONTINUE将模拟盒子的中心移到原点：初始速度：n 简单的选择： V random(-0.5, 0.5)=x, y, z标度因子:速度标度:FACTOR = SQRT ( RTEMP )DO 100 I = 1, NVX(I) = FACTOR * ( RANF(DUMMY) - 0.5 )VY(I) = FACTOR * ( RANF(DUMMY) - 0.5 )VZ(I) = FACT

14、OR * ( RANF(DUMMY) - 0.5 )100 CONTINUE随机安排初始速度：标度初始速度：SUMKX = 0.0SUMKY = 0.0SUMKZ = 0.0DO 200 I = 1, NSUMKX = SUMKX + VX(I)*2SUMKY = SUMKY + VY(I)*2SUMKZ = SUMKZ + VZ(I)*2200 CONTINUEBEITAX =SQRT(RTEMP/SUMKX)BEITAY =SQRT(RTEMP/SUMKY) BEITAZ =SQRT(RTEMP/SUMKZ)DO 300 I = 1, NVX(I) = VX(I) * BEITAXVY(I

15、) = VY(I) * BEITAYVZ(I) = VZ(I) * BEITAZ300 CONTINUE标度因子:SUMX = 0.0SUMY = 0.0SUMZ = 0.0DO 200 I = 1, NSUMX = 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, NVX(I) = VX(I) - SUMXVY(I) = VY(I) - SUMYVZ(I) = VZ(I) - SUMZ300 CONTINUE控制体系的总动量为零：n 从Maxwell分布中抽样：xxdx高斯(Gauss)分布:对于等几率随机试验(Bernoulli试验)， 大量的试验结果满足高斯分布麦克斯韦速度分布定律:由于：=x, y, z单位体积的分子再每个分量上的速度分布实际上就是 高斯分布。n 从Maxwell分布中抽样：高斯(Gauss)分布的随机数生成方法:生成随机数：i，i=1, 2, , 12SUM = 0.0DO 10 I = 1, 12SUM = SUM + RANF ( D