搜索
第1页共9页分子动力学模拟及其在金属材料中的研究进展摘要本文综述了分子动力学模拟技术的发展，介绍了分子动力学的分类、运动方程的求解、初始条件和边界条件的选取、平衡系综及其控制、感兴趣量的提取以及分子动力学模拟在金属材料中的研究进展。关键词：分子动力学模拟平衡态系综金属材料感兴趣量径向分布函数引言科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在金属材料研究中显得非常有吸引力。分子动力学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分子动力学分类分子动力学的目标是研究体系中与时间和温度等有关的性质而不只是静力学模拟中研究的构型方面。分子动力学假定原子的运动是由牛顿运动方程决定的，这意味着原子的运动第2页共9页是与特定的轨道联系在一起的。分子动力学模拟的关键问题是原子间作用势的确定，主要是求解下述牛顿运动方程组。其中Ma为原子质量，Ra为原子空间位置，t表示时间，Fα为原子间作用力。确定原子间的相互作用力Fα，也就是确定原子间作用势E(RN)。确定原子间作用势，必须知道相应的电子基态。电子基态的计算是一个非常复杂的量子多体问题，即解多体薛定愕(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势等，它们在特定的问题中均有各自的优越性。实际上，在多原子体系中一个原子的位置不同将影响空间一定范围内的电子云分布，从而影响其他原子之间的有效相互作用，因此，人们开始考虑粒子间的多体作用第3页共9页(Many-bodyeffects)，构造出多体势结构。多体势于20世纪80年代初期开始出现，Daw等在1984年首次提出了嵌人原子法(Embedded-atommethod,EAM),EAM势很好地描述了金属原子之间的相互作用，是描述金属体系最常用的一种势函数。对于由共价键结合的有机分子以及半导体材料并不适用。为更好描述各种含有共价键作用的物质，人们考虑了电子云的非球形对称，将EAM势推广到共价健材料。为此，Baskes等提出了修正嵌人原子核法(MEAM)。从某种意义上说这个模型是半经验的，因为它从局域电子密度观点出发解决全部问题，使用的参数有从实验中获得的数据(如晶格常数、转变能、体积弹性模量、弹性系数等)。（2）现代分子动力学为了克服传统分子动力学可移植性差这一缺陷，人们考虑直接从量子力学(Quantummechanics,QM)轨道理论出发获取原子间作用势。基于QM的分子动力学称之为现代分子动力学(ModernMD)，也称之为从头分子动力学(Abinitiomoleculardynamics,AIMD)。密度泛函分子动力学(DFMD)和第一原理分子动力学(FPMD)是比较常用的。DFMD是基于量子力学密度泛函理论(Densityfunctionaltheory)，直接从量子力学基本原理考虑电子云结构，模拟更为准确，可移植性更好，但计算量大。密度泛函理论是在量子理论基础上建立起来的，从波函数出发定义电子的密度，赋予波函数确切的物理意义，通过求解Schrodinger方程，确定电子的密度，再根据能量与密度的关系给出系统的能量。第一原理分子动力学(First-principlesmoleculardynamics,FPMD)是利用第一原理法对电子结构进行计算，解决材料中各元素间的成键、结合和相稳定性，材料的力学行为与电子结构和成键性质、电荷分布的主要方向等。Smargiassi等给出了一种FPMD算法，它把系统总能量进一步细分为8个部分，并利用各部分现成的显式结果，只有凡由LDA得到。这种方法最关键的一点是计算中不需要波函数的显式结构，使其计算量大大减小。对于不同的应用对象原子间势函数的势参数互不相同，势参数的确定一般有3种方法:①通过实验值(如晶格常数、弹性常数、内聚能和空位形成能等)拟合势参数；②通过蒙特卡罗方法确定势参数；③通过基于量子力学得到的各种微观信息来确定势参数。2运动方程的求解在分子动力学中，系统中原子的一系列位形是通过牛顿运动方程积分得到的。为了得到原子的运动，一般采用各种有限差分法来求解运动方程。常用的几种算法如下:(1)Verlet算法是一种目前应用最广泛的数值积分求解运动方程组的算法。由系统的哈密顿量可以推导出牛顿方程形式的运动方程组:令则可以得到如下差分方程组的形式:第4页共9页由空间坐标可以算出粒子的运动速度为:令则(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感兴趣量的提取求出模拟系统中原子在相空间的轨迹后，应用统计物理原理即可得到系统的宏观特性和结构特点。在此讨论几种常用物理参量的统计形式。第5页共9页(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为成键两原子的共用最近邻原子数；