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

分子动力学模拟及其在材料中的研究进展汇总文档视界分子动力学模拟及其在材料中的研究进展汇总分子动力学模拟及其在材料中的研究进展汇总分子动力学模拟及其在材料中的研究进展摘要：本文综述了分子动力学模拟技术的发展，介绍了分子动力学的分类、运动方程的求解、初始条件和边界条件的选取、平衡系综及其控制、感兴趣量的提取以及分子动力学模拟在材料中的研究进展。关键词：分子动力学模拟平衡态系综金属材料感兴趣量径向分布函数引言科学工作者在长期的科学研究实践中发现，当实验研究方法不能知足研究工作的需求时，用计算机模拟却能够提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但能够用来指导实验，并验证某些理论假设，进而促进理论和实验的发展。十分是在材料构成经过中很多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即能够方便地得到。这种优点使分子动力学模拟在金属材料研究中显得非常有吸引力。分子动力学MD(MolecularDynamics)模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。MD模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的构造和性质。该模拟技术主要涉及粒子运动的动力学问题,与蒙特卡罗模拟方法(简称MC)相比,分子动力学是一种“确定性方法,它所计算的是时间平均,而MC进行的是系综平均。然而根据统计力学各态历经假设,时间平均等价于系综平均。因而,两种方法严格的比拟计算能给出几乎一样的结果。经典的分子动力学方法是Alder等于1957年提出并首先在“硬球液体模型下应用，发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统会发生有液相到结晶相的转变。后来人们称这种相变为Alder相变。Rahman于1963年采用连续势模型研究了液体的分子动力学模拟。1972年Less等发展了该方法并扩展了存在速度梯度的非平衡系统。1980年Andersen等创造了恒压分子动力学方法。1983年Gillan等将该方法推广到具有温度梯度的非平衡系统，进而构成了非平衡系统分子动力学方法体系。1984年Nose等完成了恒温分子动力学方法的创立。1985年针对势函数模型化比拟困难的半导体和金属等，Car等提出了将电子论与分子动力学方法有机统一起来的第一性原理分子动力学方法。1991年Cagin等1进一步提出了应用于处理吸附问题的巨正则系综分子动力学方法。20世纪80年代后期，计算机技术飞速发展，加上多体势函数的提出与发展，使分子动力学模拟技术有了进一步的发展。1.分子动力学分类分子动力学的目的是研究体系中与时间和温度等有关的性质而不只是静力学模拟中研究的构型方面。分子动力学假定原子的运动是由牛顿运动方程决定的，这意味着原子的运动是与特定的轨道联络在一起的。分子动力学模拟的关键问题是原子间作用势确实定，主要是求解下述牛顿运动方程组。其中Ma为原子质量，Ra为原子空间位置，t表示时间，F为原子间作用力。N)。确定原子间的互相作用力F，也就是确定原子间作用势E(R确定原子间作用势，必须知道相应的电子基态。电子基态的计算是一个非常复杂的量子多体问题，即解多体薛定愕(Schrodinger)方程(式(2):式中：Etot表示系统的总能量，ri表示第i电子的空间坐标，(ri，R)是系统波函数。系统的哈密顿算子(Hamiltonian)H可表示为:其中P、p分别表示核和电子的动量算子，M、m分别表示核和电子的质量，、表示原子核的序号，Z表示电荷数。式(3)右端第一、二项分别表示核和电子的动能，第三项表示电子间的互相作用势，第四项表示核和电子的互相作用势，第五项表示核间的互相作用势。根据Born-Oppenheimer近似(电子云构造受核运动的影响极小)，系统的薛定愕方程可分离为原子核薛定愕方程和电子薛定愕方程。而电子薛定愕方程进一步可写为:其中(ri，R)为电子的波函数，E(R)的物理意义是核静止时系统的基态能量，是核坐标R。的函数，能够理解为原子间作用势。当F确定时，就能够通过求解牛顿运动方程分析系统的力学行为。事实上，求解薛定愕方程是非常困难的，因而通常是通过试验拟合或半经历解法得到原子间作用势，然后求得系统能量。也就是讲，分子动力学模拟通常是经历或半经历的。根据对原子间作用势不同的简化处理方法，分子动力学可划分为经典分子动力学和当代分子动力学。1经典分子动力学经典分子动力学(ClassicalMD)通过实验结果或经历模型确定原子间作用势，计算量较小，能够解决较大规模的问题，但是可移植性(Transferability)差。针对不同的问题，可能需要确定不同的经历参数。在20世纪80年代以前，分子动力学模拟一般都采用对势模型(Pairpotential)，该模型仅考虑近邻原子间的库仑作用力和短程互相作用，并以为系统能量为各粒子能量总和。对势能够比拟好地描绘除金属和半导体以外的几乎所有无机化合物。比拟常用的对势有硬球势、Lennard-Jones(LJ)势、Morse势、Johnson势等，它们在特定的问题中均有各自的优越性。实际上，在多原子体系中一个原子的位置不同将影响空间一定范围内的电子云分布，进而影响其他原子之间的有效互相作用，因而，人们开场考虑粒子间的多体作用(Many-bodyeffects)，构造出多体势构造。多体势于20世纪80年代初期开场出现，Daw等在1984年初次提出了嵌人原子法(Embedded-atommethod,EAM),EAM势很好地描绘了金属原子之间的互相作用，是描绘金属体系最常用的一种势函数。对于由共价键结合的有机分子以及半导体材料并不适用。为更好描绘各种含有共价键作用的物质，人们考虑了电子云的非球形对称，将EAM势推广到共价健材料。为此，Baskes等提出了修正嵌人原子核法(MEAM)。从某种意义上讲这个模型是半经历的，由于它从局域电子密度观点出发解决全部问题，使用的参数有从实验中获得的数据(如晶格常数、转变能、体积弹性模量、弹性系数等)。2当代分子动力学为了克制传统分子动力学可移植性差这一缺陷，人们考虑直接从量子力学(Quantummechanics,QM)轨道理论出发获取原子间作用势。基于QM的分子动力学称之为当代分子动力学(ModernMD)，也称之为从头分子动力学(Abinitiomolec-ulardynamics,AIMD)。密度泛函分子动力学(DFMD)和第一原理分子动力学(FPMD)是比拟常用的。DFMD是基于量子力学密度泛函理论(Densityfunctionaltheory)，直接从量子力学基本原理考虑电子云构造，模拟更为准确，可移植性更好，但计算量大。密度泛函理论是在量子理论基础上建立起来的，从波函数出发定义电子的密度，赋予波函数确切的物理意义，通过求解Schrodinger方程，确定电子的密度，再根据能量与密度的关系给出系统的能量。第一原理分子动力学(First-principlesmoleculardynamics,FPMD)是利用第一原理法对电子构造进行计算，解决材料中各元素间的成键、结合和相稳定性，材料的力学行为与电子构造和成键性质、电荷分布的主要方向等。Smargiassi等给出了一种FPMD算法，它把系统总能量进一步细分为8个部分，并利用各部分现成的显式结果，只要凡由LDA得到。这种方法最关键的一点是计算中不需要波函数的显式构造，使其计算量大大减小。对于不同的应用对象原子间势函数的势参数互不一样，势参数确实定一般有3种方法:通过实验值(如晶格常数、弹性常数、内聚能和空位构成能等)拟合势参数；通过蒙特卡罗方法确定势参数；通过基于量子力学得到的各种微观信息来确定势参数。2.运动方程的求解在分子动力学中，系统中原子的一系列位形是通过牛顿运动方程积分得到的。为了得到原子的运动，一般采用各种有限差分法来求解运动方程。常用的几种算法如下:(1)Verlet算法是一种目前应用最广泛的数值积分求解运动方程组的算法。由系统的哈密顿量能够推导出牛顿方程形式的运动方程组:令则能够得到如下差分方程组的形式:由空间坐标能够算出粒子的运动速度为:令则(8)式能够写为:利用上式则能够在计算中得到同一时间步上的空间位置和速度，并且提高了数值计算的稳定性。(2)Swope提出的Velocity-Verlet算法能够同时给出位置、速度与加速度，并且不牺牲精度。这种算法的优点是给出了显式速度项，并且计算量适中，目前应用比拟广泛。3.分子动力学模拟的初始条件和边界条件系统的初始位形和初始速度可通过实验数据、理论模型或两者的结合来确定。假如模拟系统无固定晶格构造，则每个原子的位置可用舍选法或Petropolis等方法从初始密度分布n(r)得到；每个原子的初速度可从初始温度分布Tc(r)下的Maxwell-Boltzmann分布来随机选取。对于流体分子系统，粒子的初始速度按高斯分布来选取比拟合理。物体的宏观性质是大量粒子的统计行为，模拟系统的粒子数必须非常大才能准确地再现系统的行为。分子动力学模拟采用周期性边界条件的假定。取模拟系统为中心元胞，一般为立方体，体系中含有几百至上万个粒子，使中心元胞在三维空间上重复排列，于是系统粒子的象粒子将在三维空间周期性出现，体系势能包括粒子与象粒子间的互相作用在内。应用周期性边界条件后，整个体系就变得粒子数赝无限了。周期性边界条件是一个近似，它给体系强加一个实际并不存在的长程序，假如是带电粒子系统还会人为造成极化。显然元胞所含粒子数N不能太小，否则将因边界条件的影响而极大地偏离实际情况。4.分子动力学模拟的系综在分子动力学模拟方法的应用中，存在着对两种系统状态的MD模拟。一种是对平衡态的MD模拟，另一种是对非平衡态的MD模拟，两者的差异在于平衡系统的感兴趣量及边界条件与时间无关。为了准确反映实际的物理经过，分子动力学模拟总是在一定的系综下进行，详细有微正则系综，即体系的原子数(N)、体积(V)和能量(E)都保持不变；正则系综，即体系的原子数(N)、体积(V)和温度(T)都保持不变，并且总动量为零(P=0)；等温等压系综，即系统的原子数(N)、压力(P)和温度(T)都保持不变；等压等烩系综，即保持系统的原子数(N)、压力(P)和熔值(H)都不变。5.感兴趣量的提取求出模拟系统中原子在相空间的轨迹后，应用统计物理原理即可得到系统的宏观特性和构造特点。在此讨论几种常用物理参量的统计形式。(1)径向分布函数径向分布函数(Radialdistributionfunction,RDF)或偶关联函数(Paircorrelationfunction,PCF)表征着构造的无序化程度，被广泛地应用于研究液态和非晶态的构造。径向分布函数为:其中:g(r)是范围内找到1个原子的几率，P是系统的平均数密度，R是原子位置,是Dirac符号，N为原子数.(2)静态构造因子静态构造因子(Staticstructurefactor,SF)也是判定构造无序化程度的物理量,其表达式为:其中:N代表总原子数，K为倒格矢，ri为原子的位置矢量。对理想晶体而言，其静态构造因子为1，而对理想流体，则为0。 (3)局部晶序分析局部晶序分析(Localcrystallineorder,LCO)又称应用共近邻分析或对分析技术，径向分布函数和静态构造因子都只能从总体上对构造作无序与有序的判定，不能分析原子短程排布的几何特点。为了克制这个缺点，Honeycutt-Andersen便提出了HA键型指数法(对分析技术方法)，通过HA指数法计算不同温度下原子间键合类型及指数。HA键型指数法是目前对液态、非晶态金属等无序系统的微观构造进行分析研究的一种重要方法。根据这种方法，用4位数ijkl描绘原子所属的状态:i代表两个原子的成键关系,i=1为成键,i=2为未成键；j为成键两原子的共用近期邻原子数；k代表共用近期邻原子之间的成键数；对前三位数一样而构造不同的用不同的l值来区分。(4)原子团簇类型指数法对于由特定键型种类和数目所组成的不同构型的原子多面体构造，采用“原子成团类型指数法(Cluster-typeindexmethod,CTIM)来描绘。DWQi等提出用4个数码来表示每一种相关原子团簇，第一个数码表示(与中心原子)组成原子团的原子数目(又称配位数)；第二、第三、第四个数码分别表示成键原子与中心原子组成1441,1551和1661键型的数目。用这种方法，能够很方便地表示二十面体、Frank-Kasper多面体、Bernal系多面体和其它与此相关的多面体。(5)均方位移和扩散系数粒子位移平方的平均值称为均方位移(Meansquaredisplacement,MSD)，其定义式为:其中ri(0)为原子i在零时刻的位移；ri(t)为原子i在t时刻的位移。根据爱因斯坦扩散定律，均方位移随时间的变化表征了液态金属原子的扩散行为。它与扩散系数(D)存在如下关系:式中:c为常数。(6)配位数配位数(Coordinationnumber,CN)，即原子的第一近期邻原子的个数。配位数的多少显示了该原子周围原子分布密度的大小，其常作为构造分析的一种辅助手段。(7)键取向序参数采用键取向序参数(Bond-orientationparameters)来确定金属及合金的液态结构中局域取向对称特征。在该方法中，1个原子与周围原子的成键连接与1套球谐函数相关联。这里的成键不指化学键，而是与近邻原子某种方面的规定。6.分子动力学模拟在金属材料中的应用分子动力学可划分为经典分子动力学和当代分子动力学。经典分子动力学计算量较小，能够解决较大规模的问题，但针对不同的问题，可能需要确定不同的经历参数。而当代分子动力学直接从量子力学轨道理论出发获取原子间作用势，不需要经历参数，准确性高，但计算量比拟大，一般用来解决较小规模的问题。目前这两种方法在金属材料研究中都得到了大量的应用。在经典分子动力学方面，湖南大学对含有10万个Al原子的液态金属系统在凝固经过中团簇构造的构成特性进行了模拟研究，并采用原子团类型指数法(CTIM)来描绘各种类型的团簇构造组态。同时对液态金属Na在4种不同冷速下的快速凝固经过进行了模拟跟踪研究。采用双体分布函数g(r)曲线、HoneycuttAndersen键型指数法和原子团类型指数法对凝固经过中微观构造的变化进行了分析2,3。山东大学采用EAM势研究了Au纳米团簇的熔化及经过，发现金属纳米团簇熔化首先从外表开场，逐步向中心区域推进。并计算了中介区域团簇的尺寸、熔化温度、外表能、嫡、熔等热力学量以及均方根位移(RMSD)等动力学量，为研究纳米团簇提供了定量数据。另外还发现，液态团簇在2种不同的冷速下冷却得到2种不同的固态组织:晶体团簇和非晶团簇。近来，吴爱玲等采用FS互相作用势，通过双体分布函数、键对分析技术、键取向序等多种方法，对液银快速冷凝经过的微观构造转变特性作了分析，给出了连续快速冷凝经过中液银原子间依靠互相作用力构成的独特的微观构造图像，并考察了冷却经过中体系能量和元胞体积随温度的变化。模拟结果表明，在快速冷凝经过中液态Ag没有构成bcc构造的倾向4,5。中国科学技术大学采用镶嵌原子模型研究了液态Cu在不同冷却条件下的晶化与玻璃化行为。模拟结果表明，在较慢的冷却条件下，液态Cu构成晶态；在较快的冷却条件下，液态Cu构成非晶态。与晶态Cu比拟，金属玻璃Cu有较高的能量和较大的体积，其内部晶格畸变导致了本征应力的存在6。重庆大学利用分子动力学模拟方法，研究了不同原子数Cu纳米团簇在热化和冷凝经过中构造和热力学性质的变化，模型采用的是Johnson的EAM作用势。结果表明:铜团簇的熔点和凝固点随其尺寸线性增加，并逐步向大块晶体靠拢；所有团簇的凝固点都低于熔点，出现凝固经过中的滞后现象；在熔点和凝固点附近，团簇都具有负热容特性，这是由相变前后团簇内部构造突变引起的7.SnkanthSastry等利用分子动力学模拟方法，借助L-J势研究了不同冷却条件下液态向玻璃态转变经过中的能量平均平方位移、密度自相关函数、态密度等随温度变化的规律，揭示了不同冷却条件下非晶构成时能量变化的差异及原子运动的不同。S.Ozgen利用分子动力学模拟方法，采用Sutton-Chen的EAM作用势研究的结果表明:慢速冷却条件下构成晶体，快速冷却条件下构成非晶，而且在连续慢速冷却经过中观察到在玻璃转变温度之上的某一温度时，系统的团簇处于复杂的排列状态8。俄罗斯的BlagoveshchenskiiNM9选用对势模型对Pb-Kalloys进行分子动力学模拟，同时也进行实验验证，发现实验与模拟结果一致，讲明该模拟方法的有效性和准确性。上海中科院的ZhouG10选用EAM势研究了液态金属金的凝固行为，发现冷却速率对金属的最终构造有很大的影响，不同的冷却速率产生不同的构造。土耳其的CelikFA11采用EAM势对含不同原子数的纯铝进行熔化、冷却研究，用键取向序参数、径向分布函数、配位数进行构造分析，发现除了以2048个原子为研究对象的纯铝最终相为非晶外，其它的均为混合晶体，而且该晶体的类型取决于研究系统中原子的个数。综上可知，分子动力学模拟在深人研究液体构造，揭示金属熔体的构造演变、非晶倾向及热力学性质计算等方面具有很大的发展和应用前景。金属熔体构造是一个重要的研究领域，采用分子动力学方法从原子层次上描绘熔体系统的原子组态及其凝固经过中的演变经过，进行了微观和宏观的良好结合。进一步扩大分子动力学在该领域的应用显然是凝聚态物理的一个热门发展方向。在当代分子动力学的应用方面，山东大学运用密度泛函B3LYP方法对AlmCun(m,n,用液体构造近似描绘非晶态性质。贵州大学化学系15采用分子动力学模拟方法研究了不同温度、不同压力下液态Ale03体系的构造，讨论了温度和压力对体系构造的影响，模拟计算了体系的双体相关函数、近期邻配位数、键长值以及键角分布，并给出了原子水平的瞬时构造图形，计算结果与已报道的实验结果相符。Y.Pouillon等采用第一性原理分子动力学模拟小团簇7CuOm的构造和性能。中国科学院金属研究所戚力等16通过分子动力学法对CuAg合金熔化及凝固经过进行了模拟，发现Cu40Ag60合金构成非晶态需要极大的冷却速率。RajeshC等17利用从头计算分子动力学方法模拟了小团簇Pb的几何、电子构造。边秀房等18利用从头算分子动力学方法模拟了在一定压力下液态GaSb的构造变化。7.结束语经典分子动力学和当代分子动力学都得到了广泛的应用，而原子间的作用势是分子动力学的重要因素，因而结合这两种分子动力学方法，能够考虑将试验拟合和量子力学方法结合起来获取原子间作用势函数，以愈加方便和准确地进行微观研究。近年来，随着计算机和计算数学算法的迅速发展，分子动力学模拟得到了长足的发展，将成为金属材料研究必不可少的工具。