MS动力学模拟

第3章 铁基块体非晶合金-纳米晶转变的动力学模拟过程

3.1 Discover模块

3.1.1 原子力场的分配

在使用Discover模块建立基于力场的计算中，涉及几个步骤。主要有：选择力场、指定原子类型、计算或指定电荷、选择non-bond cutoffs。

在这些步骤中，指定原子类型和计算电荷一般是自动执行的。然而，在某些情形下需要手动指定原子类型。原子定型使用预定义的规则对结构中的每个原子指定原子类型。在为特定的系统确定能量和力时，定型原子使工作者能使用正确的力场参数。通常，原子定型由Discover使用定型引擎的基本规则来自动执行，所以不需要手动原子定型。然而，在特殊情形下，人们不得不手动的定型原子，以确保它们被正确地设置。

图 3-1

1） 计算并显示原子类型：点击Edit→Atom Selection，如图3-1所示

图3-2

弹出对话框，如图3-2所示

从右边的…的元素周期表中选择Fe，再点Select，此时所建晶胞中所有Fe原子都将被选中，原子被红色线圈住即表示原子被选中。再编辑集合，点击Edit→Edit 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力场的选择

1）Energy，见图3-6。

图3-6

力场的选择：

力场是经典模拟计算的核心，因为它代表着结构中每种类型的原子与围绕着它的原子是如何相互作用的。对系统中的每个原子，力场类型都被指定了，它描述了原子的局部环境。力场包括描述属性的不同的信息，如平衡键长度和力场类型对之间的电子相互作用。常见力场有COMPASS、CVFF和PCFF。

Select下拉菜单中有三个选项：

①COMPASS 力场：COMPASS 力场是第一个把以往分别处理的有机分子体系的力场与无机分子体系的力场统一的分子力场。COMPASS 力场能够模拟小分子与高分子，一些金属离子、金属氧化物与金属。在处理有机与无机体系时，采用分类别处理的方式，不同的体系采用不同的模型，即使对于两类体系的混合，仍然能够采用合理的模型描述。

②CVFF力场：CVFF 力场全名为一致性价力场(consistant valence force field),最初以生化分子为主，适应于计算氨基酸、水及含各种官能团的分子体系。其后，经过不断的强化，CVFF 力场可适用于计算多肽、蛋白质与大量的有机分子。此力场以计算系统的结构与结合能最为准确，亦可提供合理的构型能与振动频率。

③PCFF力场：PCFF为一致性力场，增加一些金属元素的力参数，可以模拟含有相应原子的分子体系，其参数的确定除大量的实验数据外，还需要大量的量子力学计算结果。

3.1.3 非键的设置

打开Non-bond选项卡，见图3-7。

图3-7

非键作用力包括范德华力和库伦力。这里将两者都选上，为的是后期做minimizer优化原子位置时精确度更高，因为考虑了作用力因素多，即两者都考虑了。

Summation method（模拟方法）：

①Atom Based：atom based基于原子的总量，包括一个原子的截断距离，一个原子的缓冲宽度距离；为直接计算法，即直接计算原子对之间的非键相互作用，当原子对超出一定距离（截断半径cutoff distance）时，即认为原子对之间相互作用为零（注：cutoff distance指范德瓦尔斯作用力和库仑力的范围，比如：设定截断半径为5，则表示已分子或原子中心为圆心，以5为半径作圆，半径以外的作用力都不考虑）。此方法计算量较小，但是可能导致能量和其导数的不连续性。当原子对间距离在Cut Off半径附近变化时，由于前一步考虑了原子对之间的相互作用，而后一步不考虑，由此会导致能量发生跳跃。当然，对于较小的体系，则可以设置足够大的Cutoff半径来保证所有的相互作用都被考虑进来。见图3-8。

图3-8

②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(N^3/2),体系较大时，

3.1.6.1 系综简介

系综（ensemble）是指具有相同条件系统（system）的集合。平衡态的分子动力学模拟，总是在一定的系综下进行。系综是统计力学中非常重要的概念，系统的一切统计特性基本都是以系综为起点推导得到的。实际应用时，要注意选择适当的系综，如（N，T，P） 常用于研究材质的相变化等。

①在微正则系综（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）是守恒的，例如节流膨胀就是一恒焓过程。在模拟中较少见。

图3-20

3.1.6.2 系综控温机制

系综的控温：温度调控机制可以使系统的温度维持在给定值，也可以根据外界环境的温度使系统温度发生涨落。一个合理的温控机制能够产生正确的统计系综，即调温后各粒子位形发生的概率可以满足统计力学法则。系综控温机制主要有：Velocity Scale、Nose、Berendsen。

图3-21

图3-21的Thermostat下拉菜单有四个，见图3-22：

①Velocity Scale（直接速度标定法）：系统温度和粒子的速度直接相关，可以通过调整粒子的速度使系统温度维持在目标值。实际分子动力学模拟中，并不需要对每一步的速度都进行标定，而是每隔一定的积分步，对速度进行周期性的标定，从而使系统温度在目标值附近小幅波动。直接速度标定法的优点是原理简单，易于程序编制。缺点是模拟系统无法和任何一个统计力学的系综对应起来；突然的速度标定引起体系能量的突然改变，致使模拟系统和真实结构的平衡态相差较远。

②Nose：该方法可以把任何数量的原子与一个热浴耦合起来，可以消除局域的相关运动，而且可以模拟宏观系统的温度涨落现象。

Andersen:

③Berendsen控温机制: 又称Berendsen外部热浴法。其基本思想是假设系统和一个恒温的外部热浴耦合在一起，通过热浴吸收和释放能量来调节系统的温度，使之与恒温热浴保持一致。 对速度每一步进行标定，以保持温度的变化率与热浴和系统的温差（Tbath-T(t)）成比例。

图3-22

3.1.6.3 系综控压机制

图3-23

图3-23下拉菜单有3项：

①Andersen：假定系统与外界“活塞”耦合，当外部压强不能补偿系统内部压强时，“活塞”运动引起系统均匀地膨胀或收缩，最终使得系统压强等于外部压强。Andersen方法具有重要的意义，后来的各种压力控制方法基本都是基于Andersen思想发展起来的。

②Berendsen：这种方法是假想把系统与一“压浴”相耦合。

③Parrinello：这种方法允许原胞的形状与体积同时发生变化，以达到与外压平衡。这种方法是对Anderson调压方法的一种扩展，可以实现对原胞施加拉伸剪切以及混合加载情况的模拟，因此在对材料的力学性质的分子动力学模拟中，得到了广泛地应用。

3.2 Forcite模块

3.2.1 Quench（快冷）

在工具栏上点击

按钮，选择calculation，弹出对话框，如图3-24所示。

图3-24

选择Quench（快冷，淬火），再点击More…出现如图3-25所示对话框：

图3-25

再点击Dynamics options的more…出现如图3-26、3-27所示：

图3-26

图3-27

Initial velocities：第一次由于设置速度，所以只能选择Random（随即速度），第二次以及以后运行则可选择Current（当前速度）了，此时速度为上一次结束的速度。

图3-28

注意：模拟退火的时候要加力。即Include forces要选上，如图3-28所示。 运行结束后会得到一些文件，有1）3D Atomistic.xtd，这是快冷后得到的结构，见图3-29。2）Status.txt以及3D Atomistic.txt包含了快冷过程的相关参数设置以及结果数据。3）3D Atomistic Temperature.xcd描述了温度与时间的关系，见图3-30。4）3D Atomistic Energies.xcd描述了几种能量（势能、动能、非键能以及总能量）随时间的变化关系（见图3-31）等。

图3-29

图3-30

图3-31

对图3-29所示的结果做径向分布函数分析，得到如图3-32的图像，图像表明快冷结果得到非晶合金。 图3-32

3.2.2退火（anneal）

选择退火（anneal）如图3-33所示

图3-33

点击more…出现图3-34所示对话框：

图3-34

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×10×500）

设置好数据后，点击More…，出现图3-35所示对话框。

图3-35

目标温度根据快冷得到700K的结构而设定为700K，中间最高温度（Mid-cycle temperature）分别设为900K、880K、860K、850K、840K、835K、830K、825K、820K、810K十组数据。因其晶化温度大致在825K左右【】。

Advanced选项卡中需要注意的是要选上Include forces。见图3-36。

图3-36

再对十组所得结构作X衍射分析，得到10组XRD图。

3.3 Reflex 模块

如图3-37所示，调出Powder Diffraction工具 从工具栏选择Reflex

工具，或者从菜单栏选择Modules | Reflex，然后选择

Powder Diffraction，弹出对话框，见图3-38。

图3-37

图3-38

Powder Diffraction对话框由8个不同的选项卡组成，包括你需要的所有设置。

Diffractometer－设置基本的扫描设置，例如2-theta范围和线性变化； Radiation－设置不同的衍射线类型，可以选择X射线、电子和中子射线； Profiles－设置粉末衍射图显示的峰形函数并加宽显示衍射图； Sample－设置样品尺寸；

Temperature Factors－包括控制修正原子热振动对衍射图的影响； Asymmetry－控制用于修改峰形的任何不对称性修正； Experimental Data－允许你添加实验数据进行对比；

Display－设置常规的显示属性，这对控制图形数据是很重要的。 本次设置只把2θ角的宽度范围设置改为5°～90°，其余设置不变。再点击Calculate即可得出X衍射图谱。

根据X衍射图谱，再结合谢乐公式计算晶粒尺寸，得出退火温度在在830K时，晶粒尺寸为8.436nm，835K的晶粒尺寸为15.736nm。两个温度下的所得结果分别见图3-39和3-40。相应的X衍射图见3-41和3-42。（Scherrer公式，德拜-谢乐公式，表达式为D=Kλ/Bcosθ，其中K为Scherrer常数，其值为0.89；D为晶粒尺寸（nm）；B为积分半高宽度；θ为衍射角；λ为X射线波长，为0.154056 nm）

图3-39

图3-40

图3-41

图3-42

本文来源：https://www.bwwdw.com/article/ggja.html