分子动力学方法PPT课件.ppt

分子动力学方法 1 3 1基本原理 计算机模拟分类 1 随机模拟方法 优点 随机模拟方法计算的程序简单 占内存少 但是该方法难于处理非平衡态的问题 2 分子动力学方法 MolecularDynamics或简称MD 可以处理非平衡态问题 但是使用该方法的程序较复杂 计算量大 占内存也多 分子动力学方法利用牛顿古典力学来计算许多分子在相空间中的轨迹 2 3 1基本原理 分子动力学 MD 方法的出发点是物理系统的确定的微观描述 哈密顿描述方程 拉格朗日方程或牛顿运动方程 因此 分子动力学方法是用运动方程来计算多体或少体系统的性质 结果得到的既有系统的静态特性 也有动态特性 蒙特卡罗方法只能得到系统的位形特性 MD方法的具体做法是在计算机上求运动力程的数值解 为此 通过适当的格式对方程进行近似 使之适于在计算机上求数值解 其实质是计算一组分子的相空间轨道 其中每个分子各自都服从经典运动定律 这里的系统不仅是点粒子系统 也包括具有内部结构的粒子组成的系统 3 3 1基本原理 分子动力学是在原子 分子水平上求解多体问题的重要的计算机模拟方法 通过求解所有粒子的运动方程 分子动力学方法可以用于模拟与粒子运动路径相关的基本过程 在分子动力学中 粒子的运动行为是通过经典的Newton运动方程所描述 4 3 1 1粒子运动方程的数值求解 粒子体系的运动方程 Lagrangian方程Lagrangian函数的定义为则运动的Lagrangian方程为 5 3 1 1粒子运动方程的数值求解 粒子的运动方程 Lagrangian方程 6 3 1 1粒子运动方程的数值求解 单原子的牛顿运动方程 7 3 1 1粒子运动方程的数值求解 Verlet算法r t t r t v t t 1 2 a t t2 1 r t t r t v t t 1 2 a t t2 2 将上面两式相加 得到 r t t 2r t r t t a t t2 3 v t t v t a t t 1 2 b t t2 4 a t t a t b t t 5 将 5 式的b t 代入 4 得到 v t t v t 1 2 a t a t t t 6 8 3 1 1粒子运动方程的数值求解 其他求解算法 Leap frogalgorithmr t t r t v t 1 2 t tv t 1 2 t v t 1 2 t a t tBeeman salgorithmr t t r t v t t 2 3 a t t2 1 6 a t t t2v t t v t v t t 1 3 a t t 5 6 a t t 1 6 a t t t 9 3 1 2热力学量的计算 在物理系统的计算机模拟中 系综平均必须用时间平均代替 在通常的模拟中 粒子数N和体积V是固定的 给定初始位置rN 0 和初始动量pN 0 后 一个MD算法将从运动方程生成轨道 rN t pN t 轨道平均的定义为假定能量守恒 并且轨道在一切具有同一能量的相同休积内经历相同的时间 则轨道平均等于微正则系综平均 10 3 1 2热力学量的计算 孤立系统的总能量是一个守恒量 沿着分子动力学模拟生成的任何一条轨道 能量应保持不变 即 孤立系统的动能和势能不是守恒量 它们的大小沿着生成的轨道逐点变化 11 3 1 2热力学量的计算 生成的动能路径是不连续的 必须在各个时间间断点上计算动能的值以求平均 其中从平均动能可以计算系统的温度 温度是一个重要的物理量 需要加以监测 特别是在模拟的起始阶段 12 3 1 2能量均分定理 在热力学极限下 一切系综都是等同的 并且可以应用能量均分定理 热力学极限是指粒子数 或体积 趋向无穷大时的极限 一般宏观物体包含了1023个粒子 可以认为是满足热力学极限的 相变只有在热力学极限下才会发生 系统的哈密顿量为则 13 3 1 2能量均分定理 由于系统的每个粒子有三个自由度 因此假定位势在处被截断 系统内部的位形能量的轨道平均值为其中 14 3 1 2能量均分定理 由于位势被截断 总能量和势能含有误差 为了估计必须作出改正 所有的内部位形能都加到截止距离为止 尾部改正可以取其中是g r 对关联函数 它是在原点r 0处有一个粒子时 在r周围的体积元dr内找到一个粒子的概率 令n r 为离一个给定粒子的距离在r和r dr之间的平均粒子数 于是 15 3 1 2能量均分定理 其他的量也需要尾部改正 以压强的计算为例 这时维里 virial 状态方程成立 至于势能的计算 可以把积分分成两项 一项是由相互作用力程之内的贡献引起的 一项是对位势截断的改正项 长程改正项为 16 3 2分子动力学模拟的基本步骤 分子动力学模拟的实际步骤可以划分为四步 设定模拟所采用的模型 给定初始条件 趋于平衡的计算过程 宏观物理量的计算 17 3 2分子动力学模拟的基本步骤 1 模拟模型的设定硬球势Lennard Jones型势 18 根据经典物理学的规律 可以知道在系综模拟中的守恒量 微正则系综的模拟中能量 动量和角动量均为守恒量 在此系综中他们分别表示为 元胞周期性边界条件 3 2分子动力学模拟的基本步骤 19 分子动力学模拟的最小像力约定示意图 最小像力约定 20 2 给定初始条件给定粒子的初始位置和速度的数值 1 令初始位置在差分网格格子上 初始速度从玻尔兹曼分布随机抽样得到 2 令初始位置随机地偏离差分网格格子 初始速度为零 3 令初始位置随机地偏离差分网格格子 初始速度从玻尔兹曼分布随机抽样得到 3 趋于平衡使系统达到平衡 模拟中需要一个趋衡过程 在这个过程中 我们增加或从系统中移出能量 直到系统具有所要求的能量 3 2分子动力学模拟的基本步骤 21 3 3平衡态分子动力学模拟 1 微正则系综的分子动力学模拟粒子数恒定 体积恒定 能量恒定 整个系统的总动量恒等于零 分子动力学模拟步骤如下 Verlet算法 1 给定初始空间位置 2 计算在第n步时粒子所受的力 3 计算在第n 1步时所有粒子所处的空间位置 4 计算第n步的速度 5 返回到步骤 2 开始下一步的模拟计算 22 Verlet算法的速度形式 1 给定初始空间位置 2 给定初始速度 3 利用公式 计算在第n 1步时所有粒子所处的空间位置 4 计算在第n 1步时所有粒子的速度 5 返回到步骤 3 开始第n 2步的模拟计算 3 3平衡态分子动力学模拟 23 3 3平衡态分子动力学模拟 一般对于能量确定的系统不可能给出精确的初始条件 可以将初始位置固定在格子的格点上 而初始速度则由波尔兹曼分布得出 以上述方法建立的系统一般不会具有所需要的能量 需要在模拟过程中逐渐调节系统能量达到给定值 步骤为 1 对运动方程求解若干步 2 计算动能和势能 3 若能量不等于所需的值 对速度进行标度 4 从第一步开始重复 直至系统达到平衡为止 24 速度标度因子 2 正则系综的分子动力学模拟速度标度因子 3 3平衡态分子动力学模拟 25 2020 2 5 26 3 3平衡态分子动力学模拟 正则系综分子动力学的模拟具体步骤 Verlet算法的速度形式 27 1 给定初始空间位置 2 给定初始速度 3 利用公式 计算在第n 1步时所有粒子所处的空间位置 4 计算在第n 1步时所有粒子的速度 动能和速度标度因子 5 计算将速度乘以标度因子的值 并让该值作为下一次计算时 第n 1步粒子的速度 6 返回到步骤 3 开始第n 2步的模拟计算 28 3 4MD在材料科学中的应用 3 4 1高分子链动力学模拟计明娟等用高温淬火分子动力学模拟方法研究了甲硫氨酸一脑啡肽 Met enkephalin 在真空中的构象性质 经聚类分析和能量优化得到了十三个低能构象 并与吗啡进行了空间拟合 29 3 4MD在材料科学中的应用 按图l的序列构造出Met enkephalin分子并进行初始的优化 然后以每0 1ps温度增加l0K进行动力学运算 加热到800K 在此温度下进行高温淬火分子动力学模拟计算 平衡时间为10ps 模拟时间为30ps 时间步长为1 0fs 取样间隔为0 2ps 在正则 NTV 系综下运行 对平衡后的构象作进一步分析 根据它们骨架 构象分布的分析 以及所有结构中每一个残基Co C骨架片段之间均方根偏差不同性的计算产生五个等级 将均方根偏差最小间隔 0 1nm rms 0 2nm 的结构分成十三个簇 从每个簇中选取最低能量构象作为此簇代表 最后将十三个所选结构中的每一个结构利用最陡下降法和共轭梯度法分别运行1000步进行完全的能量优化 直至rms能量偏差 0 004kJ mol而收敛 然后将优化后的十三个构象与吗啡进行空间拟合 30 图3 2经QMD模拟所得的十三个低能构象 31 3 4MD在材料科学中的应用 模拟结果表明 Met enkephalin是一种构象柔性大的多肽分子 可具有多种构象 主要为Gly2 Gly3和Gly3 Phe4之间的 折叠型式 大多数的结构具有由氢键连接两个末端残基所形成的假大环 侧链和芳香环定位在假大环的同侧 酪氨酸残基中的氨基N相当于吗啡中哌啶环上的季胺N 酪氨酸残基中的酚环相当于吗啡中的酚环 苯丙氨酸残基中的芳香环相当于吗啡中环己烯环 此肽与吗啡在空间结构 药效基团 上非常相似 都满足阿片受体的要求 所以 它们可作用于同一受体 32 3 4MD在材料科学中的应用 3 4 2脆性断裂模拟单德彬等为了研究在微观尺度上单晶铜的弯曲裂纹萌生和扩展机理 建立了单晶铜弯曲变形的分子动力学模型 应用速度标定法控制温度 采用Morse势进行了单晶铜的弯曲变形分子动力学模拟 33 3 4MD在材料科学中的应用 二维分子动力学模型如图3 3 图3 3弯曲变形仿真实验模型 无凹槽模型 图3 4弯曲变形仿真实验模型 无凹槽模型 34 3 4MD在材料科学中的应用 以上模拟晶面是 100 晶面 原子位置按理想点阵排列 原子运动按牛顿原子处理 牛顿原子的作用由Morse势函数计算 仿真实施后分别在模型的左右两端3列原子上施加向下的位移 同时施加弯矩 M 计算采用速度 Verlet算法 原子的初始位置设在fcc晶体的晶格上 原子的初始速度可利用Maxwell分布确定 也可直接置为零 温度控制采用速度标定法控制在绝对零度 以避免原子热激活的复杂影响 35 图3 5单晶Cu 100 晶面无凹槽模型弯曲变形仿真实验 36 图3 6单晶Cu 100 晶面有凹槽模型弯曲变形仿真实验 37 3 4MD在材料科学中的应用 通过模拟试验结果图3 5发现 在位移和弯矩的作用下 个别原子间间隙的不断增大和应变能的不断积累最终使晶体内部出现空位 在128ps已经很明显的看到空位 这些空位位于模型中部 空位的不断合并萌生裂纹 微观裂纹的顶端又形成空位 随着时间步的推移 塑性变形加剧 裂纹逐渐萌生和扩展 发展成为贯穿模型的裂纹 形成裂纹的扩展 微观裂纹在裂纹顶端形核而造成主裂纹的连续扩展 裂尖向下萌生 后续微观裂纹的扩展类似于宏观裂纹扩展 与无凹槽模型相比 图3 6有凹槽模型其裂纹萌生和扩展的时间明显缩短 表明裂纹缺陷对断裂过程起促进作用 38 3 4MD在材料科学中的应用 研究结果表明 应变能的不断积累使晶体内部产生空位 材料的裂纹萌生于空位 空位的合并形成纳米级裂纹 后续微观裂纹的扩展类似于宏观裂纹 裂纹缺陷促进了裂纹的萌生和扩展 而且对无凹槽模型和有凹槽模型裂纹萌生和扩展的时间步长的对比研究表明 裂纹缺陷对断裂过程起促进作用 39 3 4MD在材料科学中的应用 铁中裂纹扩展的结构演化吴映飞用分子动力学模拟的方法研究了bcc Fe中I型裂纹在应力及温度场下裂尖区的结构演化问题 图3 7裂纹的几何结构示意图 40 3 4MD在材料科学中的应用 裂纹前沿沿 101 方向 裂纹面垂直于 0l0 方向 裂纹扩展沿 101 加载通过增加应力强度因子 K 实现 K1由临界应力强度因子KIC度量 基于各向同性线弹性力学的平面应变条件得到初始裂纹 进行了两组模拟 1 第一组模拟选取初始温度T 5K 应力强度因子K 1 0KIC 设在所对应的外载下存在初始裂纹 随后温度增加到l00 300 和500K 外载增加到K 2 8KIC 在这组模拟中 显示出堆垛层错 裂纹形状尖锐 且扩展速度较快 2 第二组模拟选取初始温度T 100K 应力场条件同于第一组 温度加到300和500K 外载加到K 2 8KIC 在这组模拟中 显示出了堆垛层错 位错发射 裂尖钝化和分枝 以及孪晶带 该组裂纹的扩展速度慢于第一组 41 图3 9KI 1 8KIC后 裂尖区出现了孪晶带 T 100K KI 2 5KIC 图3 8低温下裂尖区原子的ABAB型堆垛层错 T 5K KI 1 8KIC 42 3 4MD在材料科学中的应用 研究结果表明 在裂纹扩展过程中 裂尖区结构演化与初始裂纹的工作条件相关 并且依赖于温度与外载的协同作用 在相同条件下裂纹周围局域原子的能量也受初始裂纹的工作条件影响 43 3 4MD在材料科学中的应用 3 4 3表面效应对于晶界附近原子迁移的影响K Saitoh等采用分子动力学模拟方法 对微电子器件中铝芯丝材内晶界附近原子的迁移进行了分析 该模拟是基于EMT理论建立的模型 这一理论结合大量电子浓度的体效应可以验证出更加准确的晶界或表面处原子结构的研究结果 主要工作分为两部分 即晶界附近扩散性能的模拟和晶界附近原子的重新排列 它的算法是对MD时间段 2 0fs 和运动方程的积分 在该模型中 将 5 100 晶界和自由表面分别设定为xz面和yz面 44 图3 10扩散性能研究中的原子模型 45 3 4MD在材料科学中的应用 其中扩散系数D可以利用图3 10所示的模型来进行模拟 并且可以用