分子动力学模拟方法

第四章分子动力学模拟方法，1957: AlderandWainwright)1964基于刚性球势1964:用Lennard-Jones势函数法模拟液态氩性质(Rahman)1971:用分子团簇行为模拟水性质(RahmanandTillinger)1977:约束动力学方法(Rychaert，CiccottivanGunsteren)1980:恒压下的动力学方法(Andersen方法，Parrinelo-Rahman方法)1983:非平衡动力学方法(GillanandDixon)1984:恒温下的动力学方法(Nos-Hoover方法)1985:第一原理分子动力学方法(Car-parrinelo方法)1991:巨正则系综的CaginandPettit，分子动力学简史，粒子的运动依赖于经典力学(牛顿定律(F=ma)，课程内容：经典分子动力学，分子动力学方法的基础：原理：计算一组分子的相空间轨道，其中每个分子服从牛顿运动定律：初始条件：分子动力学是一个通过求解所有粒子的运动方程，分子动力学方法可以用来模拟与原子运动路径相关的基本过程。在分子动力学中，粒子的运动行为用经典的牛顿运动方程来描述。分子动力学方法是一种确定性方法。一旦确定了初始构型和速度，分子随时间的轨迹也就确定了。，分子动力学的特征方法：分子动力学的算法：有限差分法，第一，Verlet算法，粒子位置的泰勒展开：粒子位置：粒子速度：粒子加速度：r(t-t):缺点：Verlet算法在处理速度上非常笨拙，Verlet算法表达式：算法启动，设置初始位置设置初始速度扰动初始位置：计算步骤n中的力计算步骤n中的位置1:计算步骤n中的速度：重复步骤4到6，Verlet算法程序：Do100I=1，nrxnewi=2.0*rx(i)rxold(i)dtsq*ax(i)rynewi=2.0*ry(i)ryold(i)精度，误差0(4)2，每次积分计算一次力3，时间可逆缺点：1。速度误差0(2)2大，轨迹与速度无关，无法与热浴耦合，Verlet算法的优缺点：2。蛙跳算法：半步算法，1。先用当前时刻的加速度计算半个时间步长后的速度：2。计算下一步的位置：3。计算当前时刻的速度：需要T-T/2，T，TT/2，TT，T 3T/2，T 2T，V，R，V，V(-T/2)才能开始运动：LEAP-FROG算法表达式：算法开始。设置初始位置设置初始速度扰动初始速度：在步骤N中计算力在步骤N中计算速度1/2:在步骤N中计算位置1:在步骤N中计算速度：重复步骤4至7、LEAP-FROG算法的优缺点：1。提高精确度2。与速度相关的轨迹。与热水浴耦合的缺点：1 .速度近似值2。花费比Verlet运算符更多的时间。三、VelocityVerlet算法：等价于，优点：速度计算更精确，VelocityVerlet算法表达式：算法开始，初始位置指定，初始速度指定，步骤n 1的位置计算，步骤n 1的力计算，步骤n 1的速度计算：重复步骤3到5，Verlet算法三种形式的比较：Verlet，Leap-frog，VelocityVerlet，IV。预测-校正格式算法：预测器阶段：其基本思想是泰勒展开，根据新的原子位置rp，可以计算得到校正后的ac(t t)，并定义预测误差：使用该预测误差，等量校正预测位置、速度和加速度：校正阶段：预测阶段中运动方程的变换：一组向量的定义：校正阶段中运动方程的变换：形式：C0、C1、C2、C3和C0的值，取决于运动方程的顺序。一阶运动方程一阶运动方程二阶运动方程二阶运动方程积分时间步长的选择t在室温下t 1fs(飞秒10-15s)温度越高t应该降低时间步长太长会导致分子间强烈碰撞和系统数据溢出。时间步长过短将降低模拟过程搜索相空间的能力。微正则系综分子动力学是分子动力学方法中最基本的系综，它有一定数量的粒子数N、能量E和体积V，算法：指定初始位置和初始速度来积分运动方程，分几个步骤计算势能和动能，如果能量不等于所需值，将速度从缩放到，直到系统平衡为止， 微正则系综(NVE)分子动力学模拟算法流程图：给定每个分子的初始位置ri(0)和速度vi(0)，计算每个分子的力Fi和加速度ai，求解运动方程，找出每个分子移动一个时间步长后到达的位置的速度，热力学性质和统计系统的其他物理量，统计性质保持不变？ 打印结果，结束，是，否，将所有分子移动到一个新的位置并具有当前时间的速度，微正则系综分子动力学模拟程序F3解释(LJ，NVE):(无量纲量：分子动力学模拟中若干热力学量的计算：对于由n个单原子原子组成的系统：动能和温度：对比量：对于LJ流体：势能：对比量：内能：内能：内能由势能和动能组成：对比：压力： 练习：推导lj流体分子间力的表达式(fx，fy，fz及其对比):势能函数形式：力：利用对比：=x，y，z，LJ分子间的维里项：=x，y，z，以比较量的运动方程形式：(以Leap-frog算法为例)，采用比较量：最终结果为：同样，得出速度标度：根据能量均衡原理，可知：标度因子：比较量READ(*，*)NSTEP！运行次数READ(*，*)IPRINT！读取的打印步骤数(*、(a) cnfile！位文件READ(*，*)DENS！对比密度读数(*，*)RTEMP！对比温度READ(*，*)RCUT！对比截断半径READ(*，*)DT！对比时间步长，CALLREADCN(CNFILE)，初始位模式：面心立方，面心立方：(FCC):每个面在中心有一个网格点，体心立方，体心立方，(BCC):简单立方(sc):XL，初始位模式：面心立方(程序F23)，数控=(实数(n)/4.0)* *(1.0/3.0)XL=1.0/实数(NC) y=0.5 * xlr (1)=(0，0) r (2)=(0，0) nrx(I)=rx(I)-0.5 ry(I)=ry(I)-0.5 rz(I)=rz(I)-0.5100继续，将模拟框的中心移到原点：初始速度：简单选择：v=随机(-0.5，0.5)，=x，y，z，比例因子：速度比例：因子=SQRT(rtemep)do 100i=1，nvx (i)=因子*(等级-0) 比例因子：Sumx=0.0 sumy=0.0 sumz=0.0 do 200i=1，nsumx=Sumx VX(I)sumy=sumyvy(I)sumz=sumzvz(I)继续，Sumx=Sumx/实数(n)sumy=sumy/实数(n)sumz=sumz/实数(n) do300i=1，nvx(I)=VX(I)-sumxvy(I)=vy(I)-sum 大量的测试结果符合高斯分布，麦克斯韦速率分布定律：因为：=x，y，z，每单位体积的分子和每个组分的速率分布实际上是高斯分布。，从麦克斯韦分布采样：高斯分布的随机数生成方法：生成随机数：I，I=1，2，12，sum=0.0do10i=1，12 sum=sum rank(dummy)10 continuous=(sum-6.0)/4.0 R2=r * rguss=(a9 * R2 a7)* R2 a5)* R2 a3)* R2 a1)* r，高斯分布的随机数生成(程序f24)，因子=sqrt (rtemp) do 100i=1。nvx(I)=factor * gauss(dum miey)vy(I)=factor * gauss(dum miey)vz(I)=factor * gauss(dum miey)继续将总动量控制为零：与之前相同。从麦克斯韦分布中的采样分布中随机地安排初始速度：解释(LJ，NVE)微正则系综的多维模拟程序F3:维数变换：=(DENS/REAL(N)* *(1.0/3.0)RCUT=RCUT * SIGMA dt=dt * SIGMA DENS=DENS/(* * 3)，模拟框的边长为1、长程修正：微正则系综的多维模拟程序F3的解释(LJ，NVE):SR3=(/RCUT)* * 3SR 9=SR3 * * 3SIG NVE)，呼叫力(-DT，sigma，rcut，newv，NewVC，neww)呼叫移动(-DT)呼叫力(-DT，Sigma，rcut，V，VC，W)呼叫力(DT，Sigma，RCUT，V，VC，W)呼叫因子(oldk)呼叫移动(DT，Sigma，RCUT，NewV，newvc，NewW)呼叫因子(newk)，算法：差分格式：SR2=SIGSQ/RIJSQVIJ=4.0 *(SR12-SR6)WIJ .VXI=VXIDVX。VX(J)=VX(J)DVXNSUMK=SUMKX VX(1)* * 2 VY(1)* * 2 VZ(1)* * 2200 CONTINUE，BEITA=SQRT(3.0 * RTEMP/SUMK)DO300 I=1，NVX(1)=VX(1)* BEITAVY(1)=VY(1)* BEITAVZ(1)=VZ(1)* BEITA 300 CONTINUE，正则系综分子动力学(NVTMD)，直接标度恒温器约束方法(阻尼力方法)系统扩展方法(扩展系统方法)碰撞概率等于Nudt。如果一个粒子