分子动力学模拟方法VIP

第四章分子动力学模拟方法 1957年 基于刚球势的分子動力学法 AlderandWainwright 1964年 利用Lennard Jone势函数法对液态氩性质的模拟 Rahman 1971年 模拟具有分子团簇行为的水的性质 RahmanandStillinger 1977年 约束动力学方法 Rychaert CiccottivanGunsteren 1980年 恒压条件下的动力学方法 Andersen法 Parrinello Rahman法 1983年 非平衡态动力学方法 GillanandDixon 1984年 恒温条件下的动力学方法 Berendsenetal 1984年 恒温条件下的动力学方法 Nos Hoover法 1985年 第一原理分子動力学法 Car Parrinello法 1991年 巨正则系综的分子动力学方法 CaginandPettit 分子动力学简史 粒子的运动取决于经典力学 牛顿定律 F ma 课程讲解内容 经典分子动力学 ClassicalMolecularDynamics 分子动力学方法基础 原理 计算一组分子的相空间轨道 其中每个分子各自服从牛顿运动定律 初始条件 分子动力学是在原子 分子水平上求解多体问题的重要的计算机模拟方法 可以预测纳米尺度上的材料动力学特性 通过求解所有粒子的运动方程 分子动力学方法可以用于模拟与原子运动路径相关的基本过程 在分子动力学中 粒子的运动行为是通过经典的Newton运动方程所描述 分子动力学方法是确定性方法 一旦初始构型和速度确定了 分子随时间所产生的运动轨迹也就确定了 分子动力学方法特征 分子动力学的算法 有限差分方法 一 Verlet算法 粒子位置的Taylor展开式 粒子位置 粒子速度 粒子加速度 开始运动时需要r t t 缺点 Verlet算法处理速度非常笨拙 Verlet算法的表述 算法启动 规定初始位置规定初始速度扰动初始位置 计算第n步的力计算第n 1步的位置 计算第n步的速度 重复 至 Verlet算法程序 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 算法 半步算法 1 首先利用当前时刻的加速度 计算半个时间步长后的速度 2 计算下一步长时刻的位置 3 计算当前时刻的速度 t t 2 t t t 2 t t t 3 t 2 t 2 t v r v 开始运动时需要v t 2 Leap frog算法的表述 算法启动 规定初始位置规定初始速度扰动初始速度 计算第n步的力计算第n 1 2步的速度 计算第n 1步的位置 计算第n步的速度 重复 至 Leap frog算法的优缺点 优点 1 提高精确度2 轨迹与速度有关 可与热浴耦联缺点 1 速度近似2 比Verlet算子多花时间 三 VelocityVerlet算法 等价于 优点 速度计算更加准确 VelocityVerlet算法的表述 算法启动 规定初始位置规定初始速度计算第n 1步的位置 计算第n 1步的力计算第n 1步的速度 重复 至 Verlet三种形式算法的比较 Verlet Leap frog VelocityVerlet 四 预测 校正 Predictor Corrector 格式算法 预测 Predictor 阶段 其基本思想是Taylor展开 根据新的原子位置rp 可以计算获得校正后的ac t t 定义预测误差 利用此预测误差 对预测出的位置 速度 加速度等量进行校正 校正 Corrector 阶段 预测阶段运动方程的变换 定义一组矢量 校正阶段运动方程的变换 的形式 C0 C1 C2 C3的值以及 C0 取决于运动方程的阶数 一阶运动方程 二阶运动方程之一 二阶运动方程之二 五 积分时间步长 t的选择 室温下 t 1fs femtosecond10 15s 温度越高 t应该减小 太长的时间步长会造成分子间的激烈碰撞 体系数据溢出 太短的时间步长会降低模拟过程搜索相空间的能力 微正则系综分子动力学 NVEMD 它是分子动力学方法的最基本系综 具有确定的粒子数N 能量E和体积V 算法 规定初始位置和初始速度对运动方程积分若干步计算势能和动能若能量不等于所需要的值 对速度进行标度重复 至 直到系统平衡 微正则系综 NVE MD模拟算法的流程图 给定每个分子的初始位置ri 0 和速度vi 0 计算每个分子的受力Fi和加速度ai 解运动方程并求出每个分子运动一个时间步长后到达的位置所具有的速度 统计系统的热力学性质及其它物理量 统计性质不变 打印结果 结束 Yes No 移动所有分子到新的位置并具有当前时刻的速度 微正则系综MD模拟程序F3讲解 LJ NVE 无因次量 MD模拟中几个热力学量的计算 对于由N个单原子组成的系统 动能和温度 采用对比量 对于LJ流体 势能 采用对比量 内能 内能由势能和动能组成 采用对比量 压力 采用对比量 练习 推导LJ流体分子间力的表达式 fx fy fz及其对比量 势能函数形式 力 采用对比量 x y z LJ分子间的维里项 x y z 采用对比量的运动方程形式 以蛙跳 Leap frog 算法为例 采用对比量 最终得到 同理得到 速度的标度 VelocityScaling 根据能量均分原理 可知 标度因子 对比量 速度标度 或 微正则系综MD模拟程序F3讲解 LJ NVE 初始化 READ A TITLE 运行作业题目READ NSTEP 运行步数READ IPRINT 打印步数READ A CNFILE 位型文件READ DENS 对比密度READ RTEMP 对比温度READ RCUT 对比截断半径READ DT 对比时间步长 CALLREADCN CNFILE 初始位型 面心立方 face centeredcubic FCC 每面中心有一格点 体心立方 body centeredcubic BCC 简单立方 simplecubic SC 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 0DO10I 1 NCDO10J 1 NCDO10K 1 NCDO11IJ 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 11CONTINUEM M 410CONTINUE DO100I 1 NRX I RX I 0 5RY I RY I 0 5RZ I RZ I 0 5100CONTINUE 将模拟盒子的中心移到原点 初始速度 简单的选择 V random 0 5 0 5 x y z 标度因子 速度标度 FACTOR SQRT RTEMP DO100I 1 NVX I FACTOR RANF DUMMY 0 5 VY I FACTOR RANF DUMMY 0 5 VZ I FACTOR RANF DUMMY 0 5 CONTINUE 随机安排初始速度 标度初始速度 SUMKX 0 0SUMKY 0 0SUMKZ 0 0DO200I 1 NSUMKX SUMKX VX I 2SUMKY SUMKY VY I 2SUMKZ SUMKZ VZ I 2200CONTINUE BEITAX SQRT RTEMP SUMKX BEITAY SQRT RTEMP SUMKY BEITAZ SQRT RTEMP SUMKZ DO300I 1 NVX I VX I BEITAXVY I VY I BEITAYVZ I VZ I BEITAZ300CONTINUE 标度因子 SUMX 0 0SUMY 0 0SUMZ 0 0DO200I 1 NSUMX SUMX VX I SUMY SUMY VY I SUMZ SUMZ VZ I CONTINUE SUMX SUMX REAL N SUMY SUMY REAL N SUMZ SUMZ REAL N DO300I 1 NVX I VX I SUMXVY I VY I SUMYVZ I VZ I SUMZ300CONTINUE 控制体系的总动量为零 从Maxwell分布中抽样 x x dx 高斯 Gauss 分布 对于等几率随机试验 Bernoulli试验 大量的试验结果满足高斯分布 麦克斯韦速度分布定律 由于 x y z 单位体积的分子再每个分量上的速度分布实际上就是高斯分布 从Maxwell分布中抽样 高斯 Gauss 分布的随机数生成方法 生成随机数 i i 1 2 12 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 R 高斯 Gauss 分布的随机数生成 程序F24 FACTOR SQRT RTEMP DO100I 1 NVX I FACTOR GAUSS DUMMY VY I FACTOR GAUSS DUMMY VZ I FACTOR GAUSS DUMMY CONTINUE控制总动量为零 同前面一样处理 从Maxwell分布中抽样分布中随机安排初始速度 微正则系综MD模拟程序F3讲解 LJ NVE 量纲变换 SIGMA DENS REAL N 1 0 3 0 RCUT RCUT SIGMADT DT SIGMADENS DENS SIGMA 3 模拟盒子的边长为1 长程校正 微正则系综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 算法 算法启动 微正则系综MD模拟程序F3讲解 LJ NVE CALLFORCE DT SIGMA RCUT NEWV NEWVC NEWW CALLMOVE DT CALLFORCE DT SIGMA RCUT V VC W CALLFORCE DT SIGMA RCUT V VC W CALLKINET OLDK CALLMOVE DT CALLFORCE DT SIGMA RCUT NEWV NEWVC NEWW CALLKINET NEWK 算法 差分格式 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 DVX V V VIJW W WIJCALLMOVE DT DO1000I 1 NRX I RX I VX I DTRY I RY I VY I DTRZ I RZ I VZ I DT1000CONTINUE MOVE DT 速度的标定 只用于平衡阶段 SUMK 0 0DO200I 1 NSUMK SUMKX VX I 2 VY I 2 VZ I 2200CONTINUE BEITA SQRT 3 0 RTEMP SUMK DO300I 1 NVX I VX I BEITAVY I VY I BEITAVZ I VZ I BEITA300CONTINUE 正则系综分子动力学 NVTMD 具有确定的粒子数N 温度T和体积V 速度的直接标度热浴方法 AndersenThermostat 约束方法 阻尼力方法 系统扩展方法 ExtendedSystemsMethod 问题的关键 温度的约束 Nose Hoover方法 一 热浴方法 AndersenThermostat 引入一个与虚拟粒子碰撞的随机力 想象系统浸在热浴当中 系统和热浴间的相互作用强度由随机碰撞的频率决定碰撞的几率等于Nu dt如果一个粒子经历碰撞 它的速度将从约束温度下的Maxwell分布中随机抽取 总能量和总动量均不守恒 二 约束方法 是等动能 Iso Kinetics 分子动力学方法 系统的运动方程为 引入阻尼系数 以保证将温度约束在恒定值 根据高斯最小约束原理 三 Nose Hoover扩展方法 基本思想 设想原系统与一个耦合系统共同组成一个扩展系统 允许热流在原系统和耦合系统之间交换 Q 等效质量S 扩展坐标变量 热力学阻尼系数L 扩展系统的自由度 Predictor correctoralgorithmisstraightforwardVerletalgorithmisfeasible buttrickytoimplement 积分方案 updateofxdependsonpupdateofpdependsonx Nos Hoover方法正确地描述了NVT系综中的动量和位型 而等动能方法只正确地描述了后者 扩展系统的哈密顿量Hamiltonian守恒 1 复杂分子体系的势能函数形式 PotentialEnergy StretchingEnergy BendingEnergy TorsionEnergy Non BondedInteractionEnergy这些方程与描述原子或键各种不同行为的参数就构成了力场 force field UFF OPLS Amber CVFF Compass 分子模拟方法补充介绍 BondStretchingEnergy kbisthespringconstantofthebond r0isthebondlengthatequilibrium Uniquekbandr0assignedforeachbondpair i e C C O H BendingEnergy k isthespringconstantofthebend 0isthebondlengthatequilibrium Uniqueparametersforanglebendingareassignedtoeachbondedtripletofatomsbasedontheirtypes e g C C C C O C C C H etc The Hookeian potential kbandk broadenorsteepentheslopeoftheparabola Thelargerthevalueofk themoreenergyisrequiredtodeformanangle orbond fromitsequilibriumvalue TorsionEnergy Acontrolstheamplitudeofthecurvencontrolsitsperiodicityshiftstheentirecurvealongtherotationangleaxis Theparametersaredeterminedfromcurvefitting Uniqueparametersfortorsionalrotationareassignedtoeachbondedquartetofatomsbasedontheirtypes e g C C C C C O C N H C C H etc Non bondedEnergy AdeterminesthedegreetheattractivenessBdeterminesthedegreeofrepulsionqisthecharge Adetermi