分子动力学模拟及其在材料中的研究进展

1、材料计算设计基础材料计算设计基础学 号：流水号：姓 名：完成日期： 分子动力学模拟及其在材料中的研究进展摘要：本文综述了分子动力学模拟技术的发展，介绍了分子动力学的分类、运动方程的求解、初始条件和边界条件的选取、平衡系综及其控制、感兴趣量的提取以及分子动力学模拟在材料中的研究进展。关键词：分子动力学模拟 平衡态系综 金属材料 感兴趣量 径向分布函数 引言科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是在

2、材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在金属材料研究中显得非常有吸引力。分子动力学MD (Molecular Dynamics)模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。MD模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。该模拟技术主要涉及粒子运动的动力

3、学问题,与蒙特卡罗模拟方法(简称MC)相比,分子动力学是一种“确定性方法”,它所计算的是时间平均,而MC进行的是系综平均。然而按照统计力学各态历经假设,时间平均等价于系综平均。因此,两种方法严格的比较计算能给出几乎相同的结果。经典的分子动力学方法是Alder等于1957年提出并首先在“硬球”液体模型下应用，发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统会发生有液相到结晶相的转变”。后来人们称这种相变为Alder相变。Rahman于1963年采用连续势模型研究了液体的分子动力学模拟。1972年Less等发展了该方法并扩展了存在速度梯度的非平衡系统。1980年Ande

4、rsen等创造了恒压分子动力学方法。1983年Gillan等将该方法推广到具有温度梯度的非平衡系统，从而形成了非平衡系统分子动力学方法体系。1984年Nose等完成了恒温分子动力学方法的创建。1985年针对势函数模型化比较困难的半导体和金属等，Car等提出了将电子论与分子动力学方法有机统一起来的第一性原理分子动力学方法。1991年Cagin等1进一步提出了应用于处理吸附问题的巨正则系综分子动力学方法。20世纪80年代后期，计算机技术飞速发展，加上多体势函数的提出与发展，使分子动力学模拟技术有了进一步的发展。1. 分子动力学分类分子动力学的目标是研究体系中与时间和温度等有关的性质而不只是静力学模

5、拟中研究的构型方面。分子动力学假定原子的运动是由牛顿运动方程决定的，这意味着原子的运动是与特定的轨道联系在一起的。分子动力学模拟的关键问题是原子间作用势的确定，主要是求解下述牛顿运动方程组。 其中Ma为原子质量，Ra为原子空间位置，t表示时间，F为原子间作用力。确定原子间的相互作用力F，也就是确定原子间作用势E(RN)。确定原子间作用势，必须知道相应的电子基态。电子基态的计算是一个非常复杂的量子多体问题，即解多体薛定愕(Schrodinger)方程(式(2):式中：Etot表示系统的总能量，ri表示第i电子的空间坐标，(ri，R)是系统波函数。系统的哈密顿算子(Hamiltonian)H可表示

6、为:其中P、p分别表示核和电子的动量算子，M、m分别表示核和电子的质量，、表示原子核的序号，Z表示电荷数。式(3)右端第一、二项分别表示核和电子的动能，第三项表示电子间的相互作用势，第四项表示核和电子的相互作用势，第五项表示核间的相互作用势。根据Born-Oppenheimer近似(电子云结构受核运动的影响极小)，系统的薛定愕方程可分离为原子核薛定愕方程和电子薛定愕方程。而电子薛定愕方程进一步可写为:其中(ri，R)为电子的波函数，E(R)的物理意义是核静止时系统的基态能量，是核坐标R。的函数，可以理解为原子间作用势。当F确定时，就可以通过求解牛顿运动方程分析系统的力学行为。事实上，求解薛定愕

7、方程是非常困难的，因此通常是通过试验拟合或半经验解法得到原子间作用势，然后求得系统能量。也就是说，分子动力学模拟通常是经验或半经验的。根据对原子间作用势不同的简化处理方法，分子动力学可划分为经典分子动力学和现代分子动力学。（1）经典分子动力学经典分子动力学(ClassicalMD)通过实验结果或经验模型确定原子间作用势，计算量较小，可以解决较大规模的问题，但是可移植性(Transferability)差。针对不同的问题，可能需要确定不同的经验参数。在20世纪80年代以前，分子动力学模拟一般都采用对势模型(Pairpotential)，该模型仅考虑近邻原子间的库仑作用力和短程相互作用，并认为系统

8、能量为各粒子能量总和。对势可以比较好地描述除金属和半导体以外的几乎所有无机化合物。比较常用的对势有硬球势、Lennard-Jones(LJ)势、Morse势、Johnson势等，它们在特定的问题中均有各自的优越性。实际上，在多原子体系中一个原子的位置不同将影响空间一定范围内的电子云分布，从而影响其他原子之间的有效相互作用，因此，人们开始考虑粒子间的多体作用(Many-bodyeffects)，构造出多体势结构。多体势于20世纪80年代初期开始出现，Daw等在1984年首次提出了嵌人原子法(Embedded-atommethod,EAM),EAM势很好地描述了金属原子之间的相互作用，是描述金属体

9、系最常用的一种势函数。对于由共价键结合的有机分子以及半导体材料并不适用。为更好描述各种含有共价键作用的物质，人们考虑了电子云的非球形对称，将EAM势推广到共价健材料。为此，Baskes等提出了修正嵌人原子核法(MEAM)。从某种意义上说这个模型是半经验的，因为它从局域电子密度观点出发解决全部问题，使用的参数有从实验中获得的数据(如晶格常数、转变能、体积弹性模量、弹性系数等)。（2）现代分子动力学为了克服传统分子动力学可移植性差这一缺陷，人们考虑直接从量子力学(Quantummechanics,QM)轨道理论出发获取原子间作用势。基于QM的分子动力学称之为现代分子动力学(ModernMD)，也称

10、之为从头分子动力学(Abinitiomolec- ulardynamics,AIMD)。密度泛函分子动力学(DFMD)和第一原理分子动力学(FPMD)是比较常用的。DFMD是基于量子力学密度泛函理论(Densityfunctionaltheory)，直接从量子力学基本原理考虑电子云结构，模拟更为准确，可移植性更好，但计算量大。密度泛函理论是在量子理论基础上建立起来的，从波函数出发定义电子的密度，赋予波函数确切的物理意义，通过求解Schrodinger方程，确定电子的密度，再根据能量与密度的关系给出系统的能量。第一原理分子动力学(First-principlesmoleculardynamics

11、,FPMD)是利用第一原理法对电子结构进行计算，解决材料中各元素间的成键、结合和相稳定性，材料的力学行为与电子结构和成键性质、电荷分布的主要方向等。Smargiassi等给出了一种FPMD算法，它把系统总能量进一步细分为8个部分，并利用各部分现成的显式结果，只有凡由LDA得到。这种方法最关键的一点是计算中不需要波函数的显式结构，使其计算量大大减小。对于不同的应用对象原子间势函数的势参数互不相同，势参数的确定一般有3种方法:通过实验值(如晶格常数、弹性常数、内聚能和空位形成能等)拟合势参数；通过蒙特卡罗方法确定势参数；通过基于量子力学得到的各种微观信息来确定势参数。2. 运动方程的求解在分子动力

12、学中，系统中原子的一系列位形是通过牛顿运动方程积分得到的。为了得到原子的运动，一般采用各种有限差分法来求解运动方程。常用的几种算法如下:(1) Verlet算法是一种目前应用最广泛的数值积分求解运动方程组的算法。由系统的哈密顿量可以推导出牛顿方程形式的运动方程组:令 则可以得到如下差分方程组的形式:由空间坐标可以算出粒子的运动速度为:令则(8)式可以写为: 利用上式则可以在计算中得到同一时间步上的空间位置和速度，并且提高了数值计算的稳定性。(2) Swope提出的Velocity-Verlet算法可以同时给出位置、速度与加速度，并且不牺牲精度。这种算法的优点是给出了显式速度项，并且计算量适中，

13、目前应用比较广泛。3. 分子动力学模拟的初始条件和边界条件系统的初始位形和初始速度可通过实验数据、理论模型或两者的结合来确定。如果模拟系统无固定晶格结构，则每个原子的位置可用舍选法或Petropolis等方法从初始密度分布n(r)得到；每个原子的初速度可从初始温度分布Tc(r)下的Maxwell-Boltzmann分布来随机选取。对于流体分子系统，粒子的初始速度按高斯分布来选取比较合理。物体的宏观性质是大量粒子的统计行为，模拟系统的粒子数必须非常大才能准确地再现系统的行为。分子动力学模拟采用周期性边界条件的假定。取模拟系统为中心元胞，一般为立方体，体系中含有几百至上万个粒子，使中心元胞在三维空

14、间上重复排列，于是系统粒子的象粒子将在三维空间周期性出现，体系势能包括粒子与象粒子间的相互作用在内。应用周期性边界条件后，整个体系就变得粒子数赝无限了。周期性边界条件是一个近似，它给体系强加一个实际并不存在的长程序，如果是带电粒子系统还会人为造成极化。显然元胞所含粒子数N不能太小，否则将因边界条件的影响而极大地偏离实际情况。4. 分子动力学模拟的系综在分子动力学模拟方法的应用中，存在着对两种系统状态的MD模拟。一种是对平衡态的MD模拟，另一种是对非平衡态的MD模拟，两者的差别在于平衡系统的感兴趣量及边界条件与时间无关。为了准确反映实际的物理过程，分子动力学模拟总是在一定的系综下进行，具体有微正

15、则系综，即体系的原子数(N)、体积(V)和能量(E)都保持不变；正则系综，即体系的原子数(N)、体积(V)和温度(T)都保持不变，并且总动量为零(P=0)；等温等压系综，即系统的原子数(N)、压力(P)和温度(T)都保持不变；等压等烩系综，即保持系统的原子数(N)、压力(P)和熔值(H)都不变。5. 感兴趣量的提取求出模拟系统中原子在相空间的轨迹后，应用统计物理原理即可得到系统的宏观特性和结构特点。在此讨论几种常用物理参量的统计形式。(1)径向分布函数径向分布函数(Radialdistributionfunction,RDF)或偶关联函数(Paircorrelationfunction,PCF

16、)表征着结构的无序化程度，被广泛地应用于研究液态和非晶态的结构。径向分布函数为:其中:g(r)是范围内找到1个原子的几率，P是系统的平均数密度，R是原子位置,是Dirac符号，N为原子数.(2) 静态结构因子静态结构因子(Static structure factor, SF)也是判断结构无序化程度的物理量,其表达式为:其中:N代表总原子数，K为倒格矢，ri为原子的位置矢量。对理想晶体而言，其静态结构因子为1，而对理想流体，则为0。(3) 局部晶序分析局部晶序分析(Local crystalline order, LCO)又称应用共近邻分析或对分析技术，径向分布函数和静态结构因子都只能从总体上

17、对结构作无序与有序的判断，不能分析原子短程排布的几何特点。为了克服这个缺点，Honeycutt-Andersen便提出了HA键型指数法(对分析技术方法)，通过HA指数法计算不同温度下原子间键合类型及指数。HA键型指数法是目前对液态、非晶态金属等无序系统的微观结构进行分析研究的一种重要方法。根据这种方法，用4位数ijkl描述原子所属的状态:i代表两个原子的成键关系,i=1为成键,i=2为未成键；j为成键两原子的共用最近邻原子数；k代表共用最近邻原子之间的成键数；对前三位数相同而结构不同的用不同的l值来区分。(4)原子团簇类型指数法对于由特定键型种类和数目所组成的不同构型的原子多面体结构，采用“原

18、子成团类型指数法(Cluster-typeindexmethod,CTIM)”来描述。DWQi等提出用4个数码来表示每一种相关原子团簇，第一个数码表示(与中心原子)组成原子团的原子数目(又称配位数)；第二、第三、第四个数码分别表示成键原子与中心原子组成1441,1551和1661键型的数目。用这种方法，可以很方便地表示二十面体、Frank-Kasper多面体、Bernal系多面体和其它与此相关的多面体。(5)均方位移和扩散系数粒子位移平方的平均值称为均方位移(Meansquaredisplacement,MSD)，其定义式为: 其中ri(0)为原子i在零时刻的位移；ri(t )为原子i在t时刻

19、的位移。根据爱因斯坦扩散定律，均方位移随时间的变化表征了液态金属原子的扩散行为。它与扩散系数(D)存在如下关系:式中:c为常数。 (6)配位数配位数(Coordinationnumber,CN)，即原子的第一最近邻原子的个数。配位数的多少显示了该原子周围原子分布密度的大小，其常作为结构分析的一种辅助手段。(7)键取向序参数采用键取向序参数(Bond-orientationparameters)来确定金属及合金的液态结构中局域取向对称特征。在该方法中，1个原子与周围原子的成键连接与1套球谐函数相关联。这里的成键不指化学键，而是与近邻原子某种方面的规定。6. 分子动力学模拟在金属材料中的应用分子动

20、力学可划分为经典分子动力学和现代分子动力学。经典分子动力学计算量较小，可以解决较大规模的问题，但针对不同的问题，可能需要确定不同的经验参数。而现代分子动力学直接从量子力学轨道理论出发获取原子间作用势，不需要经验参数，准确性高，但计算量比较大，一般用来解决较小规模的问题。目前这两种方法在金属材料研究中都得到了大量的应用。在经典分子动力学方面，湖南大学对含有10万个Al原子的液态金属系统在凝固过程中团簇结构的形成特性进行了模拟研究，并采用原子团类型指数法(CTIM)来描述各种类型的团簇结构组态。同时对液态金属Na在4种不同冷速下的快速凝固过程进行了模拟跟踪研究。采用双体分布函数g(r)曲线、Hon

21、eycuttAndersen键型指数法和原子团类型指数法对凝固过程中微观结构的变化进行了分析2,3。山东大学采用EAM势研究了Au纳米团簇的熔化及过程，发现金属纳米团簇熔化首先从表面开始，逐步向中心区域推进。并计算了中介区域团簇的尺寸、熔化温度、表面能、嫡、熔等热力学量以及均方根位移(RMSD)等动力学量，为研究纳米团簇提供了定量数据。另外还发现，液态团簇在2种不同的冷速下冷却得到2种不同的固态组织:晶体团簇和非晶团簇。近来，吴爱玲等采用FS相互作用势，通过双体分布函数、键对分析技术、键取向序等多种方法，对液银快速冷凝过程的微观结构转变特性作了分析，给出了连续快速冷凝过程中液银原子间依靠相互作

22、用力形成的独特的微观结构图像，并考察了冷却过程中体系能量和元胞体积随温度的变化。模拟结果表明，在快速冷凝过程中液态Ag没有形成bcc结构的倾向4,5。中国科学技术大学采用镶嵌原子模型研究了液态Cu在不同冷却条件下的晶化与玻璃化行为。模拟结果表明，在较慢的冷却条件下，液态Cu形成晶态；在较快的冷却条件下，液态Cu形成非晶态。与晶态Cu比较，金属玻璃Cu有较高的能量和较大的体积，其内部晶格畸变导致了本征应力的存在6。重庆大学利用分子动力学模拟方法，研究了不同原子数Cu纳米团簇在热化和冷凝过程中结构和热力学性质的变化，模型采用的是Johnson的EAM作用势。结果表明:铜团簇的熔点和凝固点随其尺寸线

23、性增加，并逐渐向大块晶体靠拢；所有团簇的凝固点都低于熔点，出现凝固过程中的滞后现象；在熔点和凝固点附近，团簇都具有负热容特性，这是由相变前后团簇内部结构突变引起的7.SnkanthSastry等利用分子动力学模拟方法，借助L-J势研究了不同冷却条件下液态向玻璃态转变过程中的能量平均平方位移、密度自相关函数、态密度等随温度变化的规律，揭示了不同冷却条件下非晶形成时能量变化的差异及原子运动的不同。S.Ozgen利用分子动力学模拟方法，采用Sutton-Chen的EAM作用势研究的结果表明:慢速冷却条件下形成晶体，快速冷却条件下形成非晶，而且在连续慢速冷却过程中观察到在玻璃转变温度之上的某一温度时，

24、系统的团簇处于复杂的排列状态8。俄罗斯的BlagoveshchenskiiNM9选用对势模型对Pb-Kalloys进行分子动力学模拟，同时也进行实验验证，发现实验与模拟结果一致，说明该模拟方法的有效性和准确性。上海中科院的Zhou G10选用EAM势研究了液态金属金的凝固行为，发现冷却速率对金属的最终结构有很大的影响，不同的冷却速率产生不同的结构。土耳其的CelikFA11采用EAM势对含不同原子数的纯铝进行熔化、冷却研究，用键取向序参数、径向分布函数、配位数进行结构分析，发现除了以2048个原子为研究对象的纯铝最终相为非晶外，其它的均为混合晶体，而且该晶体的类型取决于研究系统中原子的个数。综上可知，分子动力学模拟在深人研究液体结构，揭示金属熔体的结构演变、非晶倾向及热力学性质计算等方面具有很大的发展和应用前景。金属熔体结构是一个重要的研究领域，采用分子动力学方法从原子层次上描述熔体系统的原子组态及其凝固过程中的演变过程，进行了微观和宏观的良好结合。进一步扩大分子动力学在该领域的应用显然是凝聚态物理的一个热门发展方向。在现代分子动力学的应用方面，山东大学运用密度泛函B3LYP方法对AlmCun(m