Discover 模块动力学模拟.docx

第 3 章 铁基块体非晶合金-纳米晶转变的动力学模拟过程3.1 Discover 模块动力学模拟3.1.1 原子力场的分配在使用 Discover 模块建立基于力场的计算中，涉及几个步骤。主要有：选 择力场、指定原子类型、计算或指定电荷、选择 non-bond cutoffs。 在这些步骤中，指定原子类型和计算电荷一般是自动执行的。然而，在某些 情形下需要手动指定原子类型。 原子定型使用预定义的规则对结构中的每个原子 指定原子类型。 在为特定的系统确定能量和力时， 定型原子使工作者能使用正确 的力场参数。 通常， 原子定型由 Discover 使用定型引擎的基本规则来自动执行， 所以不需要手动原子定型。然而，在特殊情形下，人们不得不手动的定型原子， 以确保它们被正确地设置。图 3-1 调出选择原子窗口图 3-2 选择原子窗口计算并显示原子类型：点击 EditAtom Selection，如图 3-1 所示。弹出对话 框，如图 3-2 所示。从右边的的元素周期表中选择 Fe，再点 Select，此时所建 晶胞中所有 Fe 原子都将被选中，原子被红色线圈住即表示原子被选中。再编辑 集合，点击 EditEdit Sets，如图 3-3、3-4 所示。图 3-3 编辑集合图 3-4 设定新集合弹出对话框见图 3-4，点击 New.，给原子集合设定一个名字。这里设置为 Fe，则 3D 视图中会显示“Fe”字样，再分配力场：在工具栏上点击 Discover 按 钮 ，从下拉列表中选择 Setup，显示 Discover Setup 对话框，选择 Typing 选项 卡，见图 3-5。图3-5 给原子添加力场在 Forcefield types 里选择相应原子力场，再点 Assign（分配）按钮进行原子 力 场 分 配 。 注 意 原 子 力 场 中 的 价 态 要 与 Properties Project 里 的 原 子 价 态 （Formalcharge）一致。3.1.2 体系力场的选择点击 Energy 选项卡，见图 3-6。图 3-6 Energy 选项卡图 3-7 力场下拉菜单力场的选择： 力场是经典模拟计算的核心， 因为它代表着结构中每种类型的原子与围绕着 它的原子是如何相互作用的。对系统中的每个原子，力场类型都被指定了，它描 述了原子的局部环境。 力场包括描述属性的不同的信息， 如平衡键长度和力场类 型对之间的电子相互作用。常见力场有 COMPASS、CVFF 和 PCFF。 Select 下拉菜单中有三个选项，见图 3-7： COMPASS 力场：COMPASS 力场是第一个把以往分别处理的有机分子 体系的力场与无机分子体系的力场统一的分子力场。 COMPASS 力场能够模拟小 分子与高分子，一些金属离子、金属氧化物与金属。在处理有机与无机体系时， 采用分类别处理的方式， 不同的体系采用不同的模型， 即使对于两类体系的混合， 仍然能够采用合理的模型描述。本文采用此力场。 CVFF 力场： CVFF 力场全名为一致性价力场(consistant valence force field), 最初以生化分子为主， 适应于计算氨基酸、 水及含各种官能团的分子体系。 其后， 经过不断的强化，CVFF 力场可适用于计算多肽、蛋白质与大量的有机分子。此 力场以计算系统的结构与结合能最为准确，亦可提供合理的构型能与振动频率。 PCFF 力场：PCFF 为一致性力场，增加一些金属元素的力参数，可以模 拟含有相应原子的分子体系， 其参数的确定除大量的实验数据外，还需要大量的 量子力学计算结果16。3.1.3 非键的设置打开 Non-bond 选项卡，见图 3-8。图3-8 非键选项卡非键作用力包括范德华力和库伦力。这里将两者都选上，为的是后期做 minimizer 优化原子位置时精确度更高，因为考虑的作用力因素多，即两者都考 虑了。 Summation method（模拟方法） ： Atom Based：atom based 基于原子的总量，包括一个原子的截断距离，一 个原子的缓冲宽度距离； 为直接计算法， 即直接计算原子对之间的非键相互作用， 当原子对超出一定距离（截断半径 cutoff distance）时，即认为原子对之间相互 作用为零（注：cutoff distance 指范德瓦尔斯作用力和库仑力的范围，比如：设 定截断半径为 5，则表示已分子或原子中心为圆心，以 5 为半径作圆，半径以外 的作用力都不考虑） 。此方法计算量较小，但是可能导致能量和其导数的不连续 性。 当原子对间距离在 Cutoff 半径附近变化时， 由于前一步考虑了原子对之间的 相互作用， 而后一步不考虑， 由此会导致能量发生跳跃。 当然， 对于较小的体系， 则可以设置足够大的 Cutoff 半径来保证所有的相互作用都被考虑进来。 见图 3-8。图3-9 非键图Group Based：group based 基于电子群的，总量中包括一个原子的截断距离， 一个原子的缓冲宽度距离；大多数的分子力场都包括了每个原子之间点电荷 的库仑相互作用。甚至在电中性的物种中也存在点电荷，例如水分子。点电荷实 际上反映了分子中不同原子的电负性。 在模拟中， 点电荷一般是通过电荷平衡法 （charge equilibrium）评价或者力场定义的电荷来分配的。当评价点电荷时，一 定要小心不要在使用 Cutoff 技术时引入错误的单极项。 要了解到这一点， 可以参 看如下事实：两个单极，当只有 1e.u.电荷时，在 10A 的位置上其相互作用大约 为 33Kcal；而对于由单位单极分离 1A 所形成的两个偶极，相同距离其相互作用 能不超过 0.3Kcal/mol。 很明显，忽略单极-单极相互作用会导致错误的结果，而忽略偶极-偶极相互 作用则是适度的近似。然而，如果单极相互作用处理不清的话，仍然会出问题。 当 non-bond Cutoff 使用基于原子-原子基组时，就可能发生，会人为将偶极劈裂 为两个“假”的单极（当一个偶极原子在 Cutoff 内，另一个在其外） 。这就不是忽 略了相对较小的偶极-偶极相互作用，而是人为引入了作用较大的单极-单极相互 作用。为了避免这种人为现象，Materials Studio 引入了在 Charge Groups 之上的 Cutoff。 一个“Charge Group”是一个小的原子基团，其原子彼此接近，净电荷为 0 或 者接近于 0。 在实际应用中， Charge Group 一般是常见的化学官能团， 例如羰基、 甲基或者羧酸基团的净电荷接近于中性 Charge Group。 Charge Group 之间的距离 为一个官能团中心到另一个官能团中心的距离 R，Cutoff 设置与 Atom Based 相 类似。 Ewald Summation：Ewald 是在周期性系统内计算 Non-bond 的一种技术。 Ewald 是计算长程静电相互作用能的一种算法。Ewald 加和方法比较合适于结晶 固体。原因在于无限的晶格内，Cutoff 方法会产生较大的误差。然而，此方法放 也可以用于无定形固体和溶液体系。Ewald 计算量较大，为 o(N3/2)，体系较大 时，会占用较多的内存并花费较长的时间19。 cell multipole cell based：只能用于基于指定数量层。 一般情况下，基于 Atom 适合于孤立体系，对于周期性体系计算量较小，但 是准确性较差； 基于 Group 适合于周期性和非周期性体系， 计算的准确性好一些， 计算量最小；Ewald 适合于周期性能体系，计算最为准确，但计算量最大。图3-10 Atom Based Cutoff窗口本次模拟选择 Atom Based 模拟方法，弹出对话框，见图 3-10。 Cutoff distance（截断距离） ：指的是范德瓦尔斯作用力和库仑力的范围，见 图 3-9。 Buffer width：缓冲宽度距离。 Setup 其他选项保留默认设置即可。3.1.4 结构优化在工具栏上点击 Discover 按钮 ，然后选择 Minimizer。或者从菜单栏选择Modules | Discover | Minimizer，见图 3-11。显示 Discover Minimizer 对话框，可 以进行几何结构优化计算。优化前（Minimizer） ，先查看所有原子是否都已分配 力场，如果没有，可以手动添加，在 Properties Explorer 中双击 Forcefield type， 然后修改力场类型即可。 其次在 Min 之前， 需要把晶体结构所有原子重新固定。 minimizer 只是对结构进行优化，以达到能量最小化。在作动力学（dynamics） 之前最好执行 minimizer， 因为如果不执行 minimizer， 则计算收敛时间会比较长， 能量波动会比较大，而且计算有可能出错。图3-11 几何优化窗口优化方法 Method：最陡下降法（Steepest Descent） 、共轭梯度法（Conjugate Gradient） 、牛顿方法（Newton）和综合法（Smart Minimizer） 。 Convergence level：收敛精度水平。 Maximum iteration：最大迭代数。 Optimize cell 选中的话表示优化晶胞参数和原子位置。 MS Discover 结构优化原理分子的势能一般为键合（键长、键角、二面角、扭转角等）和非键合相互作 用 （静电作用、 范德华作用等） 能量项的加和， 总势能是各类势能之和， 如下式： 总势能 = 范德华非键结势能 + 键伸缩势能 + 键角弯曲势能 + 双面角扭 曲势能 + 离平面振动势能 + 库伦静电势能 + . 除了一些简单的分子以外， 大多数的势能是分子中一些复杂形势的势能的组 合。 势能为分子中原子坐标的函数，由原子不同的坐标所得到的势能构成势能面 （Potential Energy Surface，PES）。势能越低，构象越稳定，在系统中出现的机 率越大；反之，势能越高，构象越不稳定，在系统中出现的机率越小。通常势能 面可得到许多极小值的位置，其中对应于最低能量的点称为全局最小值（Global Energy Minimum），相当于分子最稳定的构象。由势能面求最低极小值的过程 称为能量最小化（Energy Minimum），其所对应的结构为最优化结构（Optimized Structure），能量最小化过程，亦是结构优化的过程。 通过最小化算法进行结构优化时， 应避免陷入局部最小值 （local minimum） ， 也就是避免仅得到某一构象附近的相对稳定的构象， 而力求得到全局最小值， 即 实现全局优化。 分子力学的最小化算法能较快进行能量优化， 但它的局限性在于 易陷入局部势阱， 求得的往往是局部最小值， 而要寻求全局最小值只能采用系统 搜寻法或分子动力学法。在 Materials Studio 的 Discover 模块中，能量最小化算 法有以下四种： 最陡下降法（Steepest Descent），为一经典的方法，通过迭代求导，对多 变量的非线性目标函数极小化， 按能量梯度相反的方向对坐标添加位移， 即能量 函数的负梯度方向是目标函数最陡下降的方向，所以称为最陡下降法。此法计算 简单，速度快，但在极小值附近收敛性不够好，造成移动方向正交。最陡下降法 适用于优化的最初阶段。 共轭梯度法（Conjugate Gradient），在求导时，目标函数下降方向不是仅 选取最陡下降法所采用的能量函数的负梯度方向，而是选取两个共轭梯度方向， 即前次迭代时的能量函数负梯度方向与当前迭代时的能量函数负梯度方向的线 性组合。此法收敛性较好，但对分子起始结构要求较高，因此常与最陡下降法联 合使用，先用最陡下降法优化，再用共轭梯度法优化至收敛。 牛顿方法（Newton），以二阶导数方法求得极小值。此法的收敛很迅速， 也常与最陡下降法联合使用。 综合法（Smart Minimizer），该方法可以混合最陡下降法，共轭梯度法和 牛顿法进行结构优化，在 MS 中是可选择的。Smart Minimizer 中，牛顿法可以设 定最大的原子数， 如果体系的原子数大于所设定的值， 则计算是会自动地转为前 面设定的收敛法（共轭梯度法或最陡下降法），收敛精度会改为共轭梯度法的默认收敛精度值。 点开各种方法后面的 More，见图 3-12、3-13、3-14、3-15，可设定收敛精度 （Convergence），算法（Algorithm）和一维搜索（Line search，指每一次迭代中 的精度）等。图3-12 综合法窗口图3-13 最陡下降法窗口图3-14 共轭梯度法窗口图3-15 牛顿法窗口当 Job 结束后，结果被返回到 Disco Min 目录，最小化的结构被命名为 3D Atomistic.xsd，并被保存在 “3D Atomistic Disco Min”目录。 还生成一个名为 “3D Atomistic.out”的文本文档，它包含了有关计算的所有能量信息。同时还生成 “SimulationEnergies.xcd” 它显示了能量随迭代次数的变化情况， ， 由于优化是使 结构更加稳定，所以能量也随之降低，最终趋于某一值，如图 3-16 所示。本次 模拟结果前后对比见图 3-17，图3-16 能量-迭代次数函数图（a） 几何优化前（二维）（b） 几何优化后（二维）（c） 几何优化前（三维）（d） 几何优化后（三维）图3-17 结构优化前后对比结构优化后，原子会有所重排，使结构更加稳定。3.1.5 高温弛豫打开 discover 下的 Dynamics，如图 3-18 所示。图3-18 Dynamics窗口Ensemble（系综） NVE、NVT、NPT、NPH。 ： Temperature：目标温度。 Pressure：给系统所施加的压力。 Number of steps：整个动力学所运行的总步数。 Time step：每一动力学步骤所花费的时间（单步长时间） 。 Dynamics time：Number of steps Time step（总模拟时间） 。 Trajectory Save：Coordinates 表示保存坐标；Final Structure 表示只保存最终 结构；Full 表示保存所有。 Frame output every： 若输入 5000， 则表示每 5000 步输出一帧 （即晶体结构） 。 运行结束后，可以通过调用 Animation 观看三维动画，见图 3-19、3-20。图3-19 调出Animation方式图3-20 Animation操作图动画工具条可以控制三维窗口中动画文件的显示，见图 3-20。 它包含以下 命令： Play Backwards：倒映动画文件。 Step Backwards：每次向后放一帧。 Stop：停止放映。 Step Forwards：每次一帧加速放映。 Play：放映动画。 Pause：暂停放映，再按一次后继续放映。 Animation Mode：显示动画模式下拉菜单。 3.1.6.1 系综简介 系综（ensemble）是指具有相同条件系统（system）的集合。平衡态的分子 动力学模拟，总是在一定的系综下进行。系综是统计力学中非常重要的概念，系 统的一切统计特性基本都是以系综为起点推导得到的。实际应用时， 要注意选择 适当的系综，如（N，T，P） 常用于研究材质的相变化等。 NVE（微正则系综） 在微正则系综（micrononical ensemble）中，模型体系的粒子数 N、体积 V 及内能（热力学能）E（在热力学通常用 U 表示内能）保持不变。是一个孤立、 保守的系统。 值得注意的是： 体系总能量， 即势能和动能的总和， 是保持守恒的，常被用来判断积分的精度固定不变。 它对应于绝热过程， 即体系与环境没有热交 换，不存在温度 T 和压力 p 的控制因素。由于体系的能量 E 是守恒的，体系的 动能和势能之间互转化。一般说，给定能量的精确初始条件是无法得到的。能量 的调整通过对速度的标度进行， 这种标度可能使系统失去平衡， 迭代弛豫达到平 衡。 NVT 系综（正则系综） 正则系综（canonical ensemble）中，体系的粒子数 N、体积 V 及温度 T 保 持不变，且总动量保持不变。因此正则系综动力学有时也被称为恒温动力学。为 了控制体系的温度，就需要设置一个“虚拟”的热浴环境，与体系进行能量交换。 常用的热浴 （bath） 包括： Nose-Hoover， Berendsen， Andersen 以及 “velocity scaling （速度标定） ”方法等。 NPT 系综（恒温恒压系综） 恒温恒压系综中，体系的粒子数 N、压力 P、温度 T 都是恒定不变的。恒温 恒压系综允许体系的“体积”发生变化。这里的体积的变化有两种方式，一种是 只变化尺寸而保持形状（比如对于晶体来说，晶格类型维持不变，但是晶胞参数 中的 a，b，c 可以变化） ，另一种是同时变化形状和尺寸（即晶格类型和晶胞参 数都可以变化） 。压强 P 与体积共轭，控压可以通过标度系统的体积来实现。目 前有许多调压的方法都是采用的这个原理。 NPH 系综（恒焓恒压系综） NPH 系综中体系的粒子数 N、压力 P 及体系的焓 H（H=E+pV）是守恒的， 例如节流膨胀就是一恒焓过程。在模拟中较少见17。 点击 More出现如图 3-21 所示对话框，Energy deviation 表示每一步所允许 的最大能量偏离为 5000kcal/mol。图3-21 能量偏离设定NPT 系综最符合实验条件， 但是在 NPT 系综下跑动力学过程中晶胞密度变化 太大，不符合实际情况，因此本文采用 NVT 系综。 3.1.6.2 系综控温机制 系综的控温：温度调控机制可以使系统的温度维持在给定值， 也可以根据外 界环境的温度使系统温度发生涨落。 一个合理的温控机制能够产生正确的统计系 综， 即调温后各粒子位形发生的概率可以满足统计力学法则。 系综控温机制主要有：Velocity Scale、Nose、Berendsen。图 3-22 控温机制的设定图 3-23 几种控温机制图 3-22 的 Thermostat 下拉菜单有四个，见图 3-23： Velocity Scale（直接速度标定法） ：系统温度和粒子的速度直接相关，可 以通过调整粒子的速度使系统温度维持在目标值。 实际分子动力学模拟中，并不 需要对每一步的速度都进行标定，而是每隔一定的积分步，对速度进行周期性的 标定，从而使系统温度在目标值附近小幅波动。 直接速度标定法的优点是原理简 单，易于程序编制。缺点是模拟系统无法和任何一个统计力学的系综对应起来； 突然的速度标定引起体系能量的突然改变， 致使模拟系统和真实结构的平衡态相 差较远。 Nose： 该方法可以把任何数量的原子与一个热浴耦合起来， 可以消除局域 的相关运动，而且可以模拟宏观系统的温度涨落现象。 Andersen：体系与一强加了指定温度的热浴相耦合。 Berendsen 控温机制: 又称 Berendsen 外部热浴法。 其基本思想是假设系统 和一个恒温的外部热浴耦合在一起， 通过热浴吸收和释放能量来调节系统的温度， 使之与恒温热浴保持一致。 对速度每一步进行标定，以保持温度的变化率与热 浴和系统的温差（Tbath-T(t)）成比例。 当系综选定 NPT 时，控温机制应用 Nose，因此本文采用 Nose 控温机制。 3.1.6.3 系综控压机制图3-24 几种控压机制图 3-24 下拉菜单有 3 项： Andersen：假定系统与外界“活塞”耦合，当外部压强不能补偿系统内部 压强时， “活塞”运动引起系统均匀地膨胀或收缩，最终使得系统压强等于外部压强。Andersen 方法具有重要的意义，后来的各种压力控制方法基本都是基于 Andersen 思想发展起来的。 Berendsen：这种方法是假想把系统与一“压浴”相耦合。 Parrinello：这种方法允许原胞的形状与体积同时发生变化，以达到与外压 平衡。这种方法是对 Anderson 调压方法的一种扩展，可以实现对原胞施加拉伸 剪切以及混合加载情况的模拟，因此在对材料的力学性质的分子动力学模拟中， 得到了广泛地应用。 由于本文采用 NPT 系综，压力一定，所以将会看到控压机制不可选。 Dynamics 运行结束后，所得结构如图 3-25 所示。图3-25 高温弛豫后的3D图由于高温（2000K）弛豫，高于熔化温度，所以此时体系处于液态状态，因此原 子处于远程无序状态。 Project Explorer 中还生成了其他一些文件， 如能量-模拟时间曲线图， 见图 3-26。 温度-模拟时间曲线图，见图 3-27。以及一些输入输出文件。图 3-26图 3-27高温弛豫中原子通过迁移、运动或者扩散，逐步降低原来的高内能态，向稳 定的低内能态转变。因此能量随时间的推移将降低并趋于某一值，而温度逐渐稳 定在设定的 2000K 上下做微小变动。3.2 Forcite 模块动力学模拟3.2.1 Quench（快冷）在工具栏上点击 按钮， 选择 calculation， 弹出对话框， 如图 3-28 所示。图 3-28 快冷窗口图 3-29 快冷参数设置选择 Quench（快冷，淬火） ，再点击 More出现如图 3-29 所示对话框： 再点击 Dynamics options 的 more出现如图 3-30、3-31 所示：图 3-30 快冷动力学窗口图 3-31 速度下拉菜单Initial velocities： 第一次由于设置速度， 所以只能选择 Random （随即速度） ， 第二次以及以后运行则可选择 Current（当前速度）了，此时速度为上一次结束 的速度。图3-32注意：模拟退火的时候要加力。即 Include forces 要选上，如图 3-32 所示。 运行结束后会得到一些文件，有 1）3D Atomistic.xtd，这是快冷后得到的结 构，见图 3-33（a）和（b） ；2）Status.txt 以及 3D Atomistic.txt 包含了快冷过程 的相关参数设置以及结果数据； 3D Atomistic Temperature.xcd 描述了温度与时 3） 间的关系，见图 3-34；4）3D Atomistic Energies.xcd 描述了几种能量（势能、动 能、非键能以及总能量）随时间的变化关系（见图 3-35）等。（a）二维 图 3-33 快冷所得结构（b）三维图 3-34 快冷后的温度-时间曲线图3-35 快冷后的能量-时间曲线对图 3-33 所示的结果做径向分布函数分析，得到如图 3-36 的图像，表明快冷结 果得到非晶合金。图3-36 径向分布函数3.2.2 anneal（退火）选择退火（anneal）如图 3-37 所示。图 3-37 退火窗口图 3-38 退火动力学窗口点击 more出现图 3-38 所示对话框： Annealing cycles：运行一次退火所作的退火循环次数。 Initial temperature：一次退火循环的起始温度也是退火循环的终止温度。 Mid-cycle temperature：一次退火循环包括升温过程和降温过程中的最高温 度。 Heating ramps per cycle：一次循环中加热过程的温度梯度步数，冷却过程的 温度下降梯度（cooling ramps per cycle）步数与加热过程的温度梯度步数相等。 Dynamics steps per ramp：每一温度梯度的动力学步数。 Total number of steps：Annealing cycles （Heating ramps per cycle + cooling Heating ramps per cycle） Dynamics steps per ramp （即上图中的总步数=5 500） 10 设置好数据后，点击 More，出现图 3-39 所示对话框。图 3-39 退火动力学选项图 3-40 Advanced 选项目 标 温 度 根 据 快 冷 得 到 700K 的 结 构 而 设 定 为 700K， 中 间 最 高 温 度 （Mid-cycle temperature）分别设为 900K、880K、860K、850K、840K、835K、 830K、825K、820K、810K 十组数据。 。Advanced 选项卡中需要注意的是要选上 Include forces。 见图 3-40。 3-41 和 3-42 图 分别是在 830K 和 840K 温度下退火所得结构。图3-41 830K退火结构图3-42 840K退火结构图3-43 830K径向分布函数曲线图3-44 835K径向分布函数曲线再对十组所得结构作 X 衍射分析，得到 10 组 XRD 图。对图 3-41 和 3-42 所示结构分别做径向分布函