1、1计算材料实验理论基础原理及一般的模拟应用-晶体 si 的模拟计算摘要：分子动力学(Molecular Dynamics)模拟是一种用来计算一个经典多体体系的平衡和传递性质的方法，即组成系统的粒子的核心运动遵守经典力学定律.将组成系统的微观粒子(分子/原子)视为经典粒子，将所研究的粒子系统视为经典多 体系统。关键词：势函数 分子动力学 si 晶体 方均根一 分子动力学模拟的原理与方法 分子动力学(Molecular Dynamics)模拟是一种用来计算一个经典多体体系的平衡和传递性质的方法，即组成系统的粒子的核心运动遵守经典力学定律.将组成系统的微观粒子(分子/原子视为经典粒子，将所研究的粒子
2、系统视为经典多体系统。选定描述系统中粒子之间相互作用规律的势能函数及系统外加约束条件，积分所有粒子的 Newton 运动方程( 或广义坐标下的 Lagangian 和 Hamiltonian 方程) 1。 纪录在各时刻系统中粒子的位置与动量，从而模拟系统随时间演进的微观过程，进一步通过有关统计方法计算系统的各种参数和性质。这种方法是按该体系内部的动力学规律来计算并确定位形的转变。它首先需要建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性，从而得到系统的宏观性质。在这
3、 样的处理过程中我们可以看出，MD 方法中不存在任何随机因素。在 MD 方法处理过程中方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内察动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现Boltzman 的统计力学途径。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。 对于 MD 方法，向然的系综是微正则系综，这时能量是运动常量。然而，当我们想要研究温度和 (或)压力是运动常量的系统时，系统不再是封闭的。例如当温度为常量的系统可以认为系统是放置在一个热俗中。
4、当然，在 MD 方法中我们只是在想象中将系统放入热浴中。实际上，在模拟计算中具体所采取的做法是对一些自由度加以约束。例如在恒温体系的情况下，体系的平均动能是一个不变量。这时我们可以设计一个算法，使平均动能被约束在一个给定值上。由于这个约束，我们并不是在真正处理一个正则系综，而实际上仅仅是复制了这个系综的位形部分。 1分子动力学模拟的基本步骤 分子动力学的基本过程是对所有的原子或分子解牛顿运动方程并决定每一时间步的粒子的位置和动力要素。在计算机上对分子系统的 MD 模拟的实际步骤可以划分为三步;第一，设定模拟所采用的模型及初始条件;第二，趋于平衡的计算过程; 第三，宏观物理量的计算。1.1 模拟
5、模型的设定与条件设置通常模拟模型的设定与条件设置包括以下几个方面: 1.1.1 给定初始状态 2初始状态包括原子与分子的初始立体构型和所有原子的速率，这相当于赋予解牛顿运动微分方程的初始条件。启动分子动力学模拟系统的初始值包括粒子初始位置和初始速度。一般设定粒子的初始位置为 OK 时晶格格点上。除了初始立体构型 (坐标)，所有原子和分子的初始速度必须被指定以决定依据运动方程的每一原子和分子的动力学特征，设定一个原子在初始构型中的初始速度就等于决定了其动能。模拟系统的热特性过程中，温度是一个不可缺少的物理量。MD 模拟把系统内粒子视为经典粒子，经典的 Boltznann 统计联系了粒子速度与温度
6、的关系对于体系中每一原子的动能值的总和也决定了平衡状态下体系的温度，如果体 系被假定处于平衡态，那么处于体系中的原子的初始速度与设定的体系温度有关。 当系统处于平衡状态时，粒子速度服从 Maxwell 分布。即使粒子初始速度不服从该分布，经过长时间趋于平衡的过程，系统内粒子速度会自发调整到服从 Maxwell 分布。但是初始条件的合理选择将可以加快系统趋于平衡. 1.1.2 元胞、截断半径与边界条件 原则上，MD 方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统;既可以是点粒子体系，也可以是具有内部结构的体系; 处理的微观客体既可以是分子，也可以是其他
7、的微观粒子。 实际上，MD 模拟方法和随机模拟方法一样都面临着两个基本限制:一个是有限观测时间的限制;另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下的宏观性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入动力学元胞并设定周期性边界条件 2 MD 模拟通过在无限大的体系中取出一块小晶体作为模拟体系，即分子动力学元胞的体积元。对气体和液体，如果所占体积足够大，并且系统处于热平衡状态的情况下，那么这个体积的形状是无关紧要的。对于晶态的系统，元胞的形状是有影响的。为了计算简便，对于气体和液体，我们取一个立方形的体积为
8、 MD 元胞.设 MD 元胞的线度大小为 L，则其体积为 L3。由于引进这样的立方体箱子，将产生六个我们不希望出现的表面。模拟中碰撞这些箱的表面的粒子应当被反射回到元胞内部，特别是对粒子数目很少的系统。然而这些表面的存在对系统的任何一种性质都会有重大的影响。为了减小引入的表面效应，我们采用周期性边界条件。采用这种边界条件，我们就可以消除引入的表面效应，构造出一个准无穷大的体积来更精确地代表宏观系统。 目前，根据具体的研究问题不同，模拟体系的边界通常有以下几种处理方法: 1 固定边界条件 对于正则系综(NVT)的分子动力学模拟，有时采用固定边界条件，比如对于粒子与固体表面相互作用体系，在模拟中将
9、底部几层原子固定，以此模拟衬底。 2 自由边界条件 自由边界条件即对边界原子不加任何限制。对于有限的原子或分子团簇，边界条件是不必考虑的，因为团簇被认为是在真空中漂浮的，这种情况可用作自由边界条件。如果模拟原胞包含的粒子对于所研究的问题足够大，或者为了减小周期性边界条件或固定边界条件对所研究问题的影响(冲击波的反射等 )，MD 模拟中也常采用自由边界条件。 3 周期性边界条件这个边界条件就是命令基本 MD 元胞完全等同地重复无穷多次。使用周期性边界条件时，必须保证以下条件: (1)粒子从模拟元胞的一边逃出的同界条件 通过对模拟元胞无限复制，由有限的粒子来表示无限大系统，或者在某个方向上无限大。
10、常用于处理连续性介质，如大体积材料，薄膜的非法向方向等。周期性边界条件的数学表示形式为 3A(x)=A(x+nL), n=(n1,n2,n3)其中 A 为任意的可观测量， n1,n2,n3 为任意整数。当某一粒子从元胞的一端出去时，从另一端进入该元胞。如图(1). 图(1)元胞与相胞中的粒子示意图如果 ri 处的粒子 i 同 rj 处的粒子 j 之间的距离为: rij=min(r i+rj-nL), (对一切的 n)实际上这个约定就是通过满足不等式条件 rcL/2 来截断位势(r c 为截止距离)。由于对于含 N 个原子的体系，存在 N(N-1)这样的作用。所以为减小计算量，当两原子距离超过R
11、c(如 8A)时，往往忽略不计这种作用，尤其对于范氏华力这种近程力。但如果体系含静电引力，它是一种长程力，此时截断距离的数值应当选得很大，使得距离大于 L/2 的粒子的相互作用可以忽略，以避免有限尺寸效应。 图(2) 近邻原子截断示意图(2)分子动力学模拟的势能截断半径小于模拟元胞最小边长的一半。这样，一个粒子就不会和本身在其他元胞的影子有相互作用，也不会出现和另外一个粒子及其影像同时存在的情况。 1.1.3 原子间相互作用势 3 从前面介绍的 MD 模拟方法可以看出，势能函数的选取是材料的计算机模拟的核心问4题。势能函数选取的恰当与否直接关系到模拟结果的准确度。势能函数总是表示成粒子坐标(r
12、)的数学函数形式，它给出体系总能量与粒子坐标的关系。如果这个函数对体系的描述越准确，那么用它来计算所得到的物理体系的性质就与真实情况越接近。原子间相互作用势的构造对于 MD 方法的广泛应用至关重要.一个有 N 个原子的体系中任意两个原子间的相互作用势可以根据一些非经验方法得到，不过这样得到的原子间相互作用势形式较为复杂，不便用于计算机模拟，对较大的体系来说，这样处理比较耗费计算机机时。从二十世纪五十年代到 80 年代，大多数分子动力学模拟都是采用经验性径向对势描述粒子之间的相互作用。人们根据经验规律构造原子间相互作用势，比如最简单的两体势模型只考虑两体作用，忽略了多体效应，即原子间相互作用势只
13、与这两个原子间的距离有关，而与其它原子无关。由于模型简单，及其在计算中模拟运算量小，计算速度快，而得到广泛采用。但是两体势模型在直接给出总能量的同时，又要求应用体积依赖的能量项，使其具有致命的弱点。而且该方法未把原子间与角度有关的作用包括进来，因而不能模拟硅锗等共价晶体。因为硅、锗等半导体材料的原子键合强度依赖于周围原子的配置。Stoncham 等统计，有四十多种势函数形式被开发出用于硅材料的模拟研究，但比较广泛应用的只有几种.最常用的是 Stillinger-Weber 势与 Tersoff 势1.2 趋于平衡计算过程 按照上面给出的运动方程、边界条件和初始条件，就可以进行分子动力学模拟计算
14、。但是，这样计算出的系统不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。为了使系统达到平衡，模拟中需要一个趋衡过程。在这个过程中，增加或从系统中移出能量，直到系统具有所要求的能量。然后，再对运动方程中的时间向前积分若于步，若系统能量值基本不变。我们称这时系统己经达到平衡态。这段达到平衡所需的时间称为弛豫时间。系统处在由不平衡到平衡的趋衡阶段时，速度会自动进行调整到要求的Maxwell 分布。而且，如果模拟是从完美的原子晶体点阵开始的，系统处干最稳态的情况，势能处于最小值，在系 统的趋衡过程中，初始动能中的部分将转化为势能，如图 3-2 所示。所以，对于要求系统维持在常温的情况，在
15、趋衡的过程中需要不停标定速度，使系统维持在要求的温度。 图(3)平衡过程中的能量变化 二晶体 si 的实验模拟说明(一).Materials Explorer 是一个多功能分子动力学软件。立足于 Windows 平台的多功能分子动力学软件。拥有强大的分子动力学计算及 Monte Carlo 软件包，是结合应用领域来研究材料工程的有力工具。Materials Explorer 可以用来研5究有机物、高聚物、生物大分子、金属、陶瓷材料、半导体等晶体、非晶体、溶液，流体，液体和气体相变、膨胀、压缩系数、抗张强度、缺陷等。Materials Explorer 软件中包含2Body，3Body，EAM，
16、AMBER 等 63 个力场可供用户选择。 Materials Explorer 软件拥有完美的图形界面，方便使用者操作。功能： MD Cell Builder 创建非晶和有序体系 Add MD Cell 创建分层体系，如晶界和相界等 Molecule Generator 在吸附，CVD ， spattering 模拟中插入新的分子到 MD Cell 附加工具： Molecule Builder 提供所有类型的分子的创建 Crystal Builder 使用空间点群和不对称单元创建晶体结构，用于无机固体及分子晶体系统的研究 Polymer Module 模拟不定型聚合物 pdb bdl 支持输
17、入 pdb 格式的文件 各种灵活的模拟功能： NEV, NTV, NPH 和 NTP 全部使用 Parrinello-Rahman-Nose 方法 两种积分方法: Gear (predictor-corrector) 和 Hernandez (leap-frog) 使用 SHAKE 算法进行键约束计算 可以采用周期边界条件 研究异相系统（气-固界面，固体颗粒边界等） 多种分子模型：potential，rigid body，bond constraint，united atom model 使用 Parrinello-Rahman 和 Nose 方法进行温度和压力控制分子动力学模拟 初始驰豫阻止
18、在液体和非晶模拟中的爆炸 外场：静电场，磁场，重力场，包括球和弹性能力校正 电荷确定工具确定分子的原子电荷 轻松建立随机多组分液相或气相系统； 轻松进行分子模拟，从而研究晶体或外延生长，表面吸附和表面破坏； 利用复制，剪切和粘贴功能建立带有缺陷和杂质的模块系统 分析： Monitoring Module 显示温度，压力，内能以及其它热力学性质与时间的 2D 曲线图 3D Atomic Configuration Module 显示系统的快照，轨迹和动画 MSD Module 从输出数据计算均方位移（ MSD） 。显示 MSD 的 2D 图像，以及每个分子的自扩散系数。 PCF Module 计
19、算配对相关性函数和自由基分布函数，并描绘出相应的二维图 Interference Function Module 根据一对相关函数计算 X-ray 和中子衍射。 Voronoi Module 计算 Voronoi 多面体的数目和这个多面体的面数。利用 Voronoi分析对无定形固态进行表征 Internal Coordinate Module 计算键长，键角，二面角，或指定分子类型的里面角，在分布图中显示分布状态。 Velocity Auto-Correlation Function fc 是截断过渡函数 ,这使两原子间的作用在距离很远时光滑的趋向于零;bij 被称为键序函数。由于 bij 不
20、是一个常数,而是一个依赖 i,j 粒子的位置,并与 i 粒子周围其它的近邻原子8有关的三体函数项,因而反映了原子间的三体相互作用。Si 原子的相互作用的 Stillinger-Weber(S-W)是 1985 年由 Stillinger 和 Weber 提出的一个将二体相互作用项与三体相互作用项集合在一起的多体势。其表达是如下:SW (Stillinger-Weber)势己广泛用来研究硅的晶格常数、晶格动力学特性、点缺陷、位错理论研究、表面熔化、熔点和生长等许多其它方面。共价键的键能不仅依赖于键长而且与键角和配位数密切相关，这意味着原子间键的强度与原子所处的局域环境密切相关.同时，Stilii
21、inget-Webet 势比较倾向于描述四面体配位的结构，键能对键角和键长的依赖关系可以近似的用两体和三体相互作用较好地描述，但键能对配位数的依赖关系却很复杂。因此，一些有其他原子的体系，如 C. Si-Ge 原子体系，原子之间不仅仅只有 sp3 成键，还有sP1， Sp2 成键，这时 SW 势就不适合用于模拟计算了。步骤：Tools 选 Potential Parameter Editor 出现以下窗口后，点 File 选 New 之后下一窗口选 Tersoff 后选 From 点 OK，生成 Tersoff 势;在选 Stillinger-Waber 3-Body Part,生成Still
22、inger-Waber 势。4在 Materials Explorer 窗口中选 Simulation 选 Potential 后弹出的窗口中选 Assign 选Tersoff-Type 选 Tersoff 点 OK 确定，同理加 Stillinger Wber 势。95.摸拟条件设置在 1 atm 下温度从常温变到 2000的模拟条件。在 Materials Explorer 窗口中选 Conditions 对 16 项进行设置，Length 设为3000，Ensemble 项中选 NTP，Temperature 项中温度改为从 298K 到 2500K，Pressure 项中选 Varia
23、ble 并设为 1atm。106. 均方根值也称作为效值，它的计算方法是先平方、再平均、然后开方。 有时也称方均根。英语写为:Root Mean Square（RMS）. 美国传统词典的定义为： The square of the average of sqoares of a set of numbers. 即：将 N 个项的平方和除以 N 后开平方的结果，即方均根的结果。计算不同时刻下的方均根位移（MSD） 、对关联函数、扩散系数。117列出最前、中间某时刻及最后时刻的原子构型。1000k 时的原子8列出 v-t，T-t，H-t ，a 、b、c-t，、-t 的关系函数。 在 7 中运行完毕
24、后 MD status 窗口单击 OK，在 Materials Explorer 窗口中单击Results 选 Monitoring 后得出 P、V、T 关于时间 t 的关系函数；在 Monitoring窗口中选 Graph 点 2D-Graph 后选 H、abc、 点 Apply。后的图129得出两种势下的熔点。由关 T-t 系图得：Tersoff 势下的熔点：Stillingwaber 势下的熔点：10重复上述模拟，设置点、线缺陷，比较分析缺陷对熔点的影响。参考文献与1陈舜麟.计算材料科学，第一版.北京:化学工业出版社,2005.60-672罗熙淳.基于分子动力学的微纳米加工表面形成机理研究哈尔滨:哈尔滨工业大 学.2002.19-20 3MoriarytJA.Phys.Rev.lett.,995,55:1502-1505 4陈舜麟.计算材料科学，第一版.北京:化学工业出版社，2005.152-154.5 http:/
