分子动力学在非晶中研究.doc

分子动力学在非晶中的应用
摘要：为非晶材料以特殊的结构和优异的性能而受到广泛关注，本文简要介绍了非晶材料和分子动力学的发展以及分子动力学的计算分析方法，重点说明分子动力学在非晶材料研究中的应用，概述了运用分子动力学对非晶材料的结构、晶化行为、凝固过程以及原子团簇的研究，对分子动力学在非晶中的运用和未来发展做了展望。
关键词：非晶材料 分子动力学 应用

Application of Molecular Dynamics on Amorphous
Abstract: Amorphous materials attract wide interest because of their unique structures as well as super properties.The development of amorphous materials and molecular dynamics and the simulation andanalysis methods are briefly introduced in this article. The applications of molecular dynamics on amorphous were mainly described. The structure, crystallization, solidification and atomic clusters of amorphous investigated by molecular dynamics are introduced in detail. The further development and the application of molecular dynamics on amorphous were prospected.
Keywords: Amorphous , Molecular Dynamics, Application

1 引言
非晶材料作为新型材料领域中的一种，它起源于20世纪60年代。1960年，美国加州理工学院的杜威等[1]采用极速冷却法制备出了第一片厚度只有20μm的Au75Si25非晶合金（金属玻璃）薄片，这标志着非晶合金这种新型材料的诞生。它与普通的晶体材料相比，具有结构简单、无取向等特点，同时它比晶体材料具有更加优异的性能，如高强度、高硬度、高耐磨性和高耐蚀性能等，是潜力极大的工程结构材料。而这些优异的力学、物理和化学性能，主要由其特殊的微观结构决定。非晶合金的性能受到成分和微观结构的制约。不同成分的非晶合金具有不同的团簇结构和团簇堆垛方式。因此非晶态材料微观结构、非晶形成过程及其内在机制是现在的研究热点，
但是由于实验方法和实验条件的局限性，仅仅通过实验方法不能彻底了解非晶态材料的微观结构、非晶形成过程及其内在机制。当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；虽然计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是材料形成过程中许多与原子有关的微观细节，在实验中无法获得，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在材料科学研究中显得非常有吸引力[2]。
分子动力学是一种用牛顿力学规律来模拟分子运动的方法，可以计算体系的热力学量和其他宏观性质，根本问题是要确定一群有相互作用的粒子在时空中的演化规律，即确定各个粒子在什么时刻什么位置是怎么样的运动状态。其步骤为：首先建立数学模型，指定粒子运动时遵循的规律以及粒子间相作用的方式，然后用计算机计算这些粒子的运动轨道，分析出系统的各种性质。1957年Alder 等首先在硬球模型下采用分子动力学研究气体和液体的状态方程,开创了用分子动力学模拟方法研究物质宏观性质的先例。1972年Less等发展了该方法并扩展了存在速度梯度的非平衡系统。1980年Andersen等[3]创造了恒压分子动力学方法。1983年Gillan等将该方法推广到具有温度梯度的非平衡系统,从而形成了非平衡系统分子动力学方法体系。1984年Nose等完成了恒温分子动力学方法的创建。1985年针对势函数模型化比较困难的半导体和金属等，Car等[4]提出了将电子论与分子动力学方法有机统一起来的第一性原理分子动力学方法。1991年Cagin等[5]进一步提出了应用于处理吸附问题的巨正则系综分子动力学方法。
2 计算与分析方法
2.1 计算方法
2.1.1 第一性原理分子动力学
第一性原理[6]计算是指仅仅利用普朗克常数、电子质量和电量三个基本物理常数以及元素的原子序数，不借助于任何经验参数，计算体系全部电子的分子积分，求解薛定愕方程，得到体系基态的电子结构和总能量。然而由于多粒子体系的复杂性,薛定愕方程无法精确求解，必须引入一些基本近似，这就是基于绝热近似和单电子近似的密度泛函理论[7]。
2.1.2 经典分子动力学
经典分子动力学(Classical MD)通过实验结果或经验模型确定原子间作用势，计算量较小，可以解决较大规模的问题，但是可移植性(Transferability)差。针对不同的问题,可能需要确定不同的经验参数。
在20世纪80年代以前，分子动力学模拟一般都采用对势模型(Pair potential) ，该模型仅考虑近邻原子间的库仑作用力和短程相互作用，并认为系统能量为各粒子能量总和。对势可以比较好地描述除金属和半导体以外的几乎所有无机化合物。比较常用的对势有硬球势、Lennard-Jones(LJ) 势、Morse 势、Johnson势等，它们在特定的问题中均有各自的优越性。实际上，在多原子体系中，一个原子的位置不同将影响空间一定范围内的电子云分布，从而影响其他原子之间的有效相互作用，因此，人们开始考虑粒子间的多体作用(Many-body effects)，构造出多体势结构。
多体势于20世纪80年代初期开始出现，Daw等[8]在1984年首次提出了嵌入原子法( Embedded-atommethod, EAM)。EAM势很好地描述了金属原子之间的相互作用，是描述金属体系最常用的一种势函数。对于由共价键结合的有机分子以及半导体材料并不适用。为更好描述各种含有共价键作用的物质，人们考虑了电子云的非球形对称，将EAM势推广到共价健材料。为此，Baskes等提出了修正嵌入原子核法(MEAM)。从某种意义上说这个模型是半经验的，因为它从局域电子密度观点出发解决全部问题，使用的参数有从实验中获得的数据(如晶格常数、转变能、体积弹性模量、弹性系数等)。
2.1.3 现代分子动力学
为了克服传统分子动力学可移植性差这一缺陷，人们考虑直接从量子力学(Quantum mechanics, QM)轨道理论出发获取原子间作用势。基于QM的分子动力学称之为现代分子动力学(Modern MD)，也称之为从头分子动力学(Ab initio molecular dynamics, AIMD)。密度泛函分子动力学(DFMD)和第一原理分子动力学(FPMD)是比较常用的。
DFMD是基于量子力学密度泛函理论(Densityfunctionaltheory) ,直接从量子力学基本原理考虑电子云结构,模拟更为准确,可移植性更好,但计算量大。密度泛函理论是在量子理论基础上建立起来的,从波函数出发定义电子的密度,赋予波函数确切的物理意义,通过求解Schrodinger 方程,确定电子的密度,再根据能量与密度的关系给出系统的能量。
2.2 结构分析方法
2.2.1 双体分布函数
双体分布函数(pair distribution function)又称偶分布函数，是研究非晶态体系的一种重要的结构参数。双体分布函数分为Faber-Ziman[9]和Bhatia-Thomton[10]两种形式，我们一般所说的是FZ形式，计算公式为：

其中，L为模拟立方盒子的边长；NN分别为α原子、β原子的个数；n指以第i个原子为中心，在半径r到r+△r的球壳范围内β原子的个数。的物理意义是以原子为中心，在半径r到r+△r的球壳范围内β原子的几率。
BT形式的双体分布函数为gNN(r)，gCC(r)和gNC(r)。对于二元合金，BT形式的双体分布函数与FZ形式的双体分布函数存在数量上的关系：

其中g11(r), g12(r), g22(r)为FZ形式的双体分布函数，ci为i种类原子的原子浓度。从BT形式的双体分布函数与FZ形式的双体分布函数之间的关系不难得出，gNN(r)的物理意义为忽略体系的原子种类，以任意一个原子为中心在距离为r的距离上发现一个原子的几率。gCC(r)的物理意义为以一个原子为中心在距离为r的距离上发现与中心原子相同种类的原子以及相异种类原子的几率之差。gNC(r)的物理意义为以一个原子为中心在距离为r的距离上发现1种类的原子与2种类的原子的几率之差。由此看以看出，出了gNN(r)表示一种几率，gCC(r)和gNC(r)表示的均为一种几率差，在某种程度上显示的是从不同的角度表征化学序随距离的变化。三元甚至更多元合金的gNN(r)和gCC(r)可以通过以下关系得到：

其中N为元素种类。gNC(r)的表达式是不对称的，因此我们认为三元及以上的合金并不存在gNC(r)。这是因为尽管可以比较在某个距离上发现不同元素的几率，但是几率差并不是唯一的。gCC(r)是发现同类和异类的几率差，无论体系中存在多少种元素，这个几率差我唯一的
2.2.2 键角分布函数
键角分布函数表征的是一种三体分布关系。首先确定中心原子的所有最近邻原子，再以中心原子为顶角以任意两个配位原子为两端计算角度的大小，这些角度的分布即为键角分布函数。具体计算步骤为:
假设i原子有Ni个最近邻原子，j和k原子为i原子的最近邻原子，他们与i的距离分布为rij和rik，j原子和k原子的距离为rjk，那么rij和rik之间的键角大小为：

键角分布函数为：

键角分布函数反映的虽然是一种三体关系但是θ的大小与rjk的大小存在一定的关系，一般情况下，θ越大，rjk也越大。所以键角分布函数在某种程度上反映的是一种双体关系。
2.2.3 键对分析
双体分布函数是三维空间结构在一维尺度上的统计平均信息，仅用它们还无法全面准确地描述液态或非晶态结构特征。Honeycutt和Andersen[11]发展的键对分析技术(bond pair analysis)进一步揭示了局域原子构型的几何特征。该方法利用两个原子及其与周围共有原子之间的成键关系来表征这对原子，它能够有效地描述液态、非晶态、晶态及其相互转变过程中几何结构演化的特点。这里所谓的成键一般是指两个原子之间的距离小于或等于由径向分布函数曲线上第一个峰谷位置所确定的最近邻距离。在键对分析技术中，采用四个指数(i j k l)来描述原子的局域环境：以A、B两个原子为例，当A和B原子的间距小于其成键距离时，认为这两个原子成键，则i=l，反之i=2；j是同时与A和B原子成键的最近邻原子数；k为这j个原子之间的成键数目；l作为附加指数，区别当前三个指数一致时不同的成键情况。
图2-3示意性地给出了几种典型的存在于液态、晶态或非晶态中的原子键对：fcc晶体结构以1441键对为特征键对；hcp以1421和1422键对为特征键对；bcc晶体结构以1441和1661键对为特征键对；二十面体结构以1541和1551键对为特征键对，其与1431键对是液态和非晶态中大量存在的典型键对。然而，必须看到，键对反映的是仅仅是团族结构的部分信息。正二十面体由12个1551组成，正十二面体数目多表明1551键对数目必然很多。但是1551数目很多并不能表明正二十面体数目很多。键对本质上表明的是团簇的对称性，1551数目的多少表明的是团族五次对称的状况而非二十面体团族数目的多少。

图1 常见的几种原子键对
2.2.4 Voronoi多面体分析
键对分析技术揭示了团族内部的对称性，但是还不能具体地给出团簇的几何图象及其特征，为此需要引入多面体分析。无序结构中Voronoi多面体的作用与晶体中的Wigner-Seitz元胞类似。它的定义为中心原子与近邻原子连线的垂直平分面所围成的最小封闭凸多面体。
Voronoi多面体确定的配位原子一般均在双体分布函数的第一峰内，因此，双体分布函数与Voronoi之间的关系就是Voronoi反映的是双体分布函数第一峰内各个原子分布的形状。Voronoi多面体是由多个键对组成的，它的指数与键对的指数存在着密切的联系。Voronoi多面体中面的边数一般与键对指数中共用原子j的个数相等。但是Voronoi多面体并不能反映组成键对的饱和性,因此假如在Voronoi多面体指数中加入组成键对的总的不饱和指数可能会更好的反映团簇的结构[12]。
3 分子动力学在非晶研究中的应用
3.1 在非晶材料结构研究中的应用
非晶材料表现出来的优异性能与其结构有着密不可分的联系，因此研究非晶材料的结构有助于我们更好地开发新的优异的非晶材料。李家云[13]采用分子动力学模拟的方法，研究了CuZr二元非晶在过冷态的结构变化，发现CuZr合金中五重对称性占据了液态结构的主导地位，随着温度的变化，CuZr合金的液态结构表现出很强的结构稳定性，这样的结构稳定性是与CuZr非晶强的非晶形成能力相对应的。苏俊亮等[14]人运用从头算分子动力学方法模拟了合金Al94-xNi6Lax(x=3-9)的液态和淬冷形成体系的结构，对其结构进行了分析，发现冷却速度足够快时液态合金可以保持其长程无序结构形成非晶态，且二者结构相似，可以用液态体系来近似描述非晶态性质。结构显示，Al原子与Ni原子之间的结合最为紧密，Al原子与自身之间的结合次之，Al原子与La原子之间的结合最弱。另外La原子的增加减缓了体系中个原子的运动，导致扩散系数减小，最终使个体系都处于亚稳的非晶态。
3.2 在非晶态材料晶化过程研究中的应用
非晶态合金在结构上处于热力学亚稳状态，在一定条件下要发生晶化转变。研究非晶态结构的晶化行为，对于深入理解非平衡相变机制和评价金属与合金的玻璃形成能力及热稳定性都有很大的帮助。
王荣山等[15]人运用分子动力学模拟的方法，采用紧束缚势函数，研究了非晶Cu在升温条件下的晶化行为，分析了系统能量、体积、双体关联函数和局部结构的变化，并结合键对分析方法计算了不同弛豫时间下典型短程有序结构的分布。结果表明：在非晶Cu升温的最初阶段，原子运动未必造成短程结构的进一步规则化；结构转变初期，首先发生1431和1541键向1421键的转变，1421键型数量在400K以上则基本呈现直线上升的趋势，接近600K时达到最大值，此后1421键的数量随温度的上升而下降，熔化时其数量急剧减少。王宇等[16]人利用分子动力学方法对非晶纯镍材料压缩变形过程中纳米晶化现象进行了模拟，研究了非晶变形过程中绝热温升对非晶晶化的影响，结果表明：绝热温升不是导致非晶晶化的主要因素。从微观结构演化的角度考察了非晶晶化过程中晶粒的形核和长大，分析发现，应变导致非晶态金属从亚稳态结构逐渐向稳态结构转变，在系统内部的部分短程序原子团逐步合并形成小的晶核，随着应变的增加，晶核逐渐长大，形成一定尺寸的纳米晶粒。刘浏等[17]人用分子动力学方法模拟了Ti-Al 合金非晶的形成与晶化过程。通过对不同冷却速率淬火过程的模拟，得到了不同合金成分的非晶转变冷速，表明40-90at%Al 的合金非晶形成能力强，临界冷却速率低。构型分析表明：非晶合金主要由二十面体团簇和初晶团簇构成，非晶退火晶化的微观机理研究说明晶化的过程实际上是体系缺陷和二十面体团簇的数量降低、初晶数量和体积增大的过程。
3.3 在非晶材料凝固过程研究中的应用
对于非晶材料来说，凝固是加工过程中重要的一步。研究非晶材料凝固过程对于获得优良性能的非晶材料具有重要意义。
邱克强等[18]采用分子动力学模拟了二元合金Cu66Ti34的凝固过程。原子间作用采用 GEAM势，利用偶关联函数，均方位移(MSD)等分析方法，研究Cu66Ti34合金在4×1013K/s 冷却速度下的玻璃化转变温度、原子的扩散行为。GEAM 势函数能够较好地反映Cu66Ti34合金玻璃化转变的动力学性质，通过偶分布函数第一谷的最小值与第一峰最大值之比获得的玻璃转变温度为600 K，与相近成分Cu50Ti50的实验值接近；在模拟条件下根据均方位移数值，可以知道Cu的扩散能力要比Ti的强，且都在600 K开始呈现明显冻结状态。合金的生成焓随温度的降低而降低，比热容与温度具有二次函数关系。可以利用Cp函数图像来估计玻璃化转变的发生。
朱志雄等[19]基于TB势函数，采用Materials Explorer分子动力学软件包较好地描述二元液态合金Ni3Al和NiAl凝固过程的结构变化特征，分析了冷却过程中不同温度下的偶分布函数、能量和体积的变化。在模拟液态合金凝固过程中，冷却速率对降温过程中液体合金的结构转变快慢和最后结构有重要影响。在较快的冷却速率(4×1013K/s) 下，Ni3Al形成了非晶结构；在相对较慢的冷却速率(4×1011 K/s)下，Ni3Al在1100 K左右结构开始发生了变化，最终形成了晶体结构。在冷却速率为 4×1013 K/s 时，NiAl 形成了非晶，在冷却速率为 4×1011 K/s 时，NiAl在810 K 有部分晶体开始形成。
厉虹等[20]以过冷液体的微观结构及原子在过冷液体中的运动情况为基础，利用分子动力学模拟方法，选取GEAM势函数来描述原子之间的相互作用力，结合偶关联函数，均方位移(MSD)等方法，通过模拟计算，对金属熔体Cu66Ti34分别在4×1011 K/s、4×1013 K/s两个冷速下，从1600K到300K的凝固过程进行了研究。结果表明：在4×1013 K/s冷速下，体系形成非晶结构，玻璃化转变温度为600K，由MSD图像可以看出在凝固过程的刚开始阶段，束缚力相同，故认为界面是平面，过一段时间后，粒子所受束缚不同，所以此时的界面应表现为非平面界面。在4×1011 K/s冷速下，体系在800 K时部分发生晶化。根据MSD数值可知，在4×1011 K/s时粒子的扩散能力要比在4×1013 K/s时的大，凝固速度降低两个数量级，而MSD的数量级增加两个数量级，说明在较低的凝固速度下，凝固时间增大，粒子有足够的时间进行重排与迁移。
3.4 在非晶材料原子团簇研究中的应用
近年来，人们对非晶材料结构展开了大量的研究，普遍认为在非晶的微观结构中存在着一些被称为原子团的基本结构单元，这些基本的结构单元对非晶的形成有着重要的作用，因此对于原子团簇的研究也是现在的热点。
刘俊等[21]采用NPT分子动力学模拟方法，应用周期边界条件，模拟了Cu-50%Ni(原子分数)合金熔体在不同冷却过程中原子团簇的演变情况，给出了以1×1014 K/s冷速冷却至室温时Cu-50%Ni非晶体系中存在的各种结构单元，并研究了Cu、Ni原子在这些结构单元中的排列情况。结果表明，较高冷速下形成的非晶具有较高的能量和较高的非晶转变温度(Tg)。在冷却过程中，原子间的短程作用逐渐加强。PCF图第一峰值逐渐增大。1551键对在非晶体系中占主导地位，且受冷速影响较大。体系中除了正二十面体外，还存在着各种缺陷多面体，其中含1551键对较多的缺陷多面体其数目也较多。相对于正二十面体,缺陷二十面体结构不稳定,在冷却到一定温度以下数目有下降的趋势，FK多面体与Bernal多面体数目始终很少。大原子(Cu)易于占据多面体顶点，而小原子(Ni)则倾向于占据各多面体中心较大的空隙。冷却速度对于非晶的形成与结构有很大的影响，随着冷速的降低，各多面体数目均有不同程度的下降。
王丽等[22]采用常温常压分子动力学模拟技术，模拟了液态Ni3Al中原子团簇在快速凝固条件下的演变过程，模型采用的是TB(tight binding)作用势。用偶分布函数、键对和多面体等结构参数来描述快速凝固条件下团簇种类和数量的变化，并将团簇结构可视化。在2000 K下，液态Ni3Al中团簇数量较少，且都是由缺陷二十面体构成；在4×1013K·S-1的冷速下，团簇的数量随温度的降低不断增加，且出现完整二十面体团簇，体系最终形成了由二十面体和缺陷二十面体团簇网络所组成的非晶结构。
庞彪等[23]通过径向分布函数和键对分析，发现在快速降温形成的非晶结构存在缺陷，随着 温度的降低，时间的延长，体系中完整二十面体个数增加，非晶中缺陷二十面体转化成完整二十面体，结构变得更加稳定。

5 结论与展望
与其他模拟方法相比，分子动力学模拟有正确的理论基础，保证了模拟结果的精确性；可以检测与时间依存的性质和行为，用分子动力学模拟可以处理与时间依存的动态现象,因此一些与时间有关的宏观量如扩散系数的模拟必须应用分子动力学。目前其发展已比较成熟，已问世的相关软件可以解决很多问题。在非晶的结构等信息的研究中，利用分子动力学方法进行研究可以解决很多以往解决不了或者有难度的问题。
计算机技术的飞速发展和诸多科学家对各算法的改进，势必将日益推进分子动力学模拟技术的发展，期望将来能够将分子动力学模拟推广至较长时间、更多的分子系统的模拟。
参考文献
[1] W.T.Klement, R.H. Willens, P.Duwez. Non-crystalline structure in solidified gold-silicon alloys.Nature,1960,187:869.
[2] 赵素，李金富，周尧和. 分子动力学模拟及其在材料科学中的应用[J]. 材料导报，2007,21(4):5-9
[3] Anderson HC. Molecular dynamics simulations at constant pressand/or temperature[J]. J ChemPhys,1980,72:2384
[4] Car R, Parrinello M. Unified approach for molecular dynamics and density-functional theory[J]. Phys Rev Lett ,1985,55:2471
[5] 吴兴惠,项金钟.现代材料计算与设计教程[M]. 北京:电子工业出版社,2002.122
[6] 廖沐真,吴国是,刘洪霖. 量子化学从头计算方法[M].北京,清华大学出版社,1984.
[7] G. Kresse, J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev, B 47(1993)558
[8] Daw M S, Baskes M I. Embedded-atommethod:Derivation and application to impurities, surfaces, and other defects in metals[J]. Phys Rev B, 1984,29:6443
[9] T. E. Faber, J. M. Ziman. A theory of the electrical properties of liquid metal, Phil Mag.11(1965)153
[10] A. B. Bhatia and D. E. Thornton, Structural Aspects of the Electrical Resistivity of Binary
Alloys,Phys. Rev. B 2(1970)3004
[11] J. D. Honeycutt, H. C. Andersen, Molecular dynamics study of melting and freezing of small Lennard-Jones clusters J. Phys. Chem. 91(1987)4950
[12] 潘少鹏. Fe基合金液态及非晶态结构的分子动力学研究[D]. 山东大学.
[13] 李家云. 分子动力学模拟研究CuZr大块非晶在过冷态的液态结构稳定性[J]. 哈尔滨师范大学自然科学学报，2013,29(6):71-73
[14] 苏俊亮. 液态和非晶态Al94-xNi6Lax合金结构的从头计算分子动力学模拟[D].南京师范大学，2011.
[15] 王荣山, 侯怀宇, 陈国良. 非晶 Cu 在晶化过程中的分子动力学模拟[J]. 金属学报, 2009, 45(6): 692-696.
[16] 王宇，王秀喜，王海龙. 非晶材料压缩变形中纳米晶化现象的分子动力学模拟[J].金属学报，2006,42(10):1071-1074.
[17] 刘浏,贾兴旺,付勇. 分子动力学模拟Ti-Al合金非晶的形成与晶化 [J].西华大学学报,2009, 26(3):6-10
[18] 邱克强，厉虹，孙晶，等. Cu66Ti34非晶合金凝固过程的分子动力学模拟[J].中国有色金属学报,2011, 21(9):2151-2156.
[19] 朱志雄，张鸿，刘超峰，等. Ni-Al合金凝固过程的分子动力学模拟[J].中国有色金属学报,2009，19(8):1409-1416.
[20] 厉虹.非晶合金凝固的界面结构及其热力学特征的研究[D].沈阳工业大学，2010..
[21] 刘俊，赵九洲. Cu-50%Ni合金快速凝固过程中原子团簇演变的分子动力学模拟[J].金属学报,2005,41(2):219-224
[22] 王丽，衣粟，边秀房. Ni3Al合金液态与非晶中的原子团簇[J].物理化学学报,2002,18(4):297-301.
[23] 庞彪,李家云,等. 分子动力学模拟研究Cu60Zr40非晶中的二十面体结构[J].哈尔滨师范大学自然科学学报,2012,28(5):35-37.

10