分子模拟的方法及其应用

第40卷第5期当代化工Vol.40，No.52011年5月ContemporaryChemicalIndustryMay，2011基金项目：国家自然科学基金（批准号：50730009）。作者简介：李春艳，女，硕士研究生，主要从事分子力学模拟研究。E-mail:[e-mail]。分子模拟的方法及其应用李春艳，刘华，刘波涛（中北大学物理系，山西太原030001）摘要：综述了分子模拟的方法及其在材料科学中的应用，介绍了分子模拟的基本原理、牛顿运动方程及其有限差分算法、势函数的发展、平衡态系综及其温度和压力的控制，最后还指出了分子模拟的应用及其进一步的研究方向。关键词：分子模拟;有限差分法;原子间相互作用势;材料科学中图分类号：O642文献标识码：A文章编号：1671-0460（2011）05-0494-03MethodofMolecularSimulationandItsApplicationLIChun-yan,LIUHua,LIUBo-tao（NorthUniversityofChina,ShanxiTaiyuan030001,China）Abstract:Methodofmolecularsimulationanditsapplicationwerereviewed.Theprincipleofmolecularsimulation,Newtonianequationofmotionwereintroducedaswellasitsfinitedifferencetechnique,progressofpotentialfunction,realizationofequilibriumensemblesandthetemperatureandpressurecontrol.Directionsofresearchesandapplicationsofmolecularsimulationwereputforward.Keywords:Molecularsimulation;Finitedifferencetechnique;Interatomicpotential;Materialscience在高分子领域中计算机模拟的应用从20世纪60年代至今已发展到了一个崭新的阶段[1-2]。这个新阶段的特点是，计算机模拟不仅能提供定性的描述，而且能模拟出高分子材料结构与性能的定量的结果；模拟分子体系算法的发展，主要是从描述简单的非真实分子体系的算法到复杂的真实分子体系的算法。该发展克服了用蒙特卡洛法仅能够描述不同温度下分子结构的特征，但不能描述随温度变化而引起体系宏观物理性质的变化中分子结构的变化过程。1分子模拟方法1.1蒙特卡洛方法(简称MC)的基本原理MC法早在20世纪50年代就被应用到高分子科学的研究之中。先驱性的工作是Wall在为研究高分子链的排除体积问题时所做的蒙特卡洛模拟[3-5]。基本模拟过程是在一定系统条件下，将系统内粒子进行随机的位移、转动，或粒子在两相间转移位置。根据给定的分子势能函数，进行粒子间内能的加和。采用Metropolis取样方法[6-7]，生成一系列体系的微观粒子随机构型，从而逐渐趋近于平衡时的Boltzmann分布。目前蒙特卡洛的应用几乎涉及到高分子科学中的各个领域，为现代高分子科学理论基础的建立和发展起到了十分重要的推动作用。1.2分子动力学方法(简称MD)的基本原理分子动力学方法(MD)是另一种主要的计算机模拟方法，在材料科学、物理、化学等学科的各个领域得到广泛应用。模拟过程是在一定体系及分子势能函数条件下，从计算分子间作用力着手，求解牛顿运动方程，得到体系中各分子微观状态随时间的变化，再将分子的位置和动量组成的微观状态对时间平均，即可求出体系的压力、能量、粘度等宏观性质以及组成粒子的空间分布等微观结构[8]。2分子模拟的计算2.1分子动力学的运动方程分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理,对一个由N个粒子构成的孤立体系,粒子的运动由牛顿运动方程决定,也就是2d=-(,,)122drimVrrriNt∇（1）式中：mi，ri—第i个原子的质量和位置；▽i=-d/dri,V(r1,r2,rN)—体系所处的势。2.2分子动力学运动方程的求解在分子动力学中，为了得到原子的运动，一般采用以下几种有限差分法来求解运动方程。常用的有：第40卷第5期李春燕，等：分子模拟的方法及其应用495Verlet算法[9]是在60年代后期出现的,对扩散分子的质心运动的积分是最稳定的也是最常用的数值方法。Swope提出的Velocity–Verlet算法[10]可以同时给出位置、速度与加速度，并且不牺牲精度。这种算法的优点是给出了显速度项，并且计算量适中，应用比较广泛。Hockney提出的“蛙跳”Leap–forg算法[11]是Verlet算法的变化，这种算法与Verlet算法相比有两个优点:（1）包括显速度项;（2）收敛速度快,计算量小。这种算法明显的缺陷是位置和速度不同步。Beeman提出的Beeman算法[12]运用了更精确的速度表达式。它的表示式比Verlet算法复杂得多，因此计算量较大。Gear的预测——校正算法[13]分为三步来完成:首先,根据Taylor展开,预测新的位置、速度和加速度。然后,根据新的计算的力计算加速度。这个加速度再由与Taylor级数展开式中的加速度进行比较。这种方法的缺点就是占有计算机的内存大。3分子间的作用势3.1对势在分子动力学模拟的初期，人们经常采用的是对势模型有以下几种：（1）在1957年Alder和Wainwright使用间断对势[14]:当r≤r1时,E=∞;当r＞r1时,E=0。（2）连续对势Lennard-Jones(L-J)势是最老的原子间作用势之一，是在1924年由J.E.Lennard-Jone提出来的[15]，L-J势的基本形式为：()1264LJRRRσσφε⎡⎤⎛⎞⎛⎞=−⎢⎥⎜⎟⎜⎟⎝⎠⎝⎠⎢⎥⎣⎦（2）式中：σ、ε—经验常数。（3）Born-Lande(B-J)势的形式如下[16]：()24πijijijmijijZZebrrrφε=+（3）（4）Morse势[17]：根据双原子分子的振动谱[18]，提出了指数形式的相互作用势：()()()rrrAeBeαβφ−−=−（4）它有4个参数A,B,α和β，与L-J势的普适形式相类似。（5）Johnson势的形式如下[19]：()()3ijijnijnnijnrArBCrDφ=−−+−（5）3.2多体势多体势于20世纪80年代初期开始出现。多体势于20世纪80年代初期开始出现。常用的多体势函数模型主要有EAM势、FS势、TB势等。3.2.1EAM势Daw和Baskes在1984年首次提出了嵌入原子法（EAM）[20]。EAM势的基本思想是把晶体的总势能分成两部分：一部分是晶格点阵上的原子核之间的相互作用对势，另一部分是原子核镶嵌在电子云中的嵌入能，它代表多体相互作用。在嵌入原子法中，系统的总势能表示为：()()tot,12iiijijiijEFrρφ=+∑∑（6）式中：F是嵌入能：第2项是对势项，根据需要可以取不同的形式。ρi是除第i个原子以外的所有其它原子的核外电子在第i个原子处产生的电子云密度之和，可以表示为：()iiijjirρρ≠=∑，其中，ρi（rij）是第j个原子的核外电子在第i个原子处贡献的电荷密度，rij是第i个原子与第j个原子之间的距离。对于不同的金属，嵌入能函数和对势函数需要通过拟合金属的宏观参数来确定。3.2.2FS势Finnis等根据金属能带的紧束缚理论，形成了一种在数学上等同于EAM的势函数[21]，并给出了多体势的函数形式，把嵌入能函数设为平方根形式，即FS势。其中的多体项及对势()iiijjirρρ≠=∑项分别为：()()()231ijkkijKijkrARrHRrρ==−−∑()()()631ijijkkijkijkrarrHrrσ==−−∑（7）式中：当x0时，H(x)=0；当x0时，H(x)=1。Ak、RK、aK、rK为常数，且有R1＞R2,r1＞r2＞…＞r6。它们的值随金属物质的不同而有所不同。3.2.3TB势基于EAM势的势函数还有很多种[22-23]。为了把将EAM势推广到共价键材料中去，需考虑到电子云的非球形对称。于是Baskes等[24-26]提出了修正型嵌入原子核法(MEAM)。总的电子密度为下式：()()()()3220hhiiihtρρ==∑（8）给出。其中t(h)i是常系数。在MEAM里，镶嵌能项采取了如下形式：F(ρ)=AEρlnρ,其中A，E是系数，因物质不同而不同。4分子动力学模拟中平衡系综的控制方法平衡态分子动力学模拟是在一定的系综下进行的,经常用的平衡系综是NVT或NPT系综。在这两种系综中,牵涉到控制温度和压力的几种技术,分别介496当代化工2011年5月绍如下：4.1控温方法在平衡系综中，经常要把温度调整到期望值。而体系的有效温度是由动能决定的：k=∑miv2i/2=fkBT/2式中：f为体系的自由度，KB为波尔兹曼常数。控温常用的方法有:（1）速度标定法[27]系统的温度与动能存在如下的关系：()2B/23/2NkiiCiEmvNNKT//==−∑（9）式中：N—原子数；NC—约束数；KB—波尔兹曼常数；vi—原子i的速度。这种方法的基本思想是：如果t时刻的温度是T（t）,速度乘以因子λ后,温度的变化为：ΔT=（λ2-1）T（t）T为所控制体系的温度。（2）Berendsen热浴法法[28]这种也称恒温槽方法，是假设所模拟的体系与一个恒温槽连在一起,那么两者之间就可以通过热交换而使所模拟的体系达到恒温目的,方法如下:()C1/tTTλτ/=+∆−,式中,参数Г表征系统与恒温槽之间的热交换速率,Δt为分子动力学模拟的时间步长,通过vnew=λvold校正即可保持体系的温度在T0附近振动,而参数Г可用于控制这个振动幅度。（3）Nose-Hoover方法法[29-30]Nose-Hoover热浴法是通过改变模拟体系的哈密顿量来实现控温的,其基本思想就是在哈密顿量里加入一个假想的项来代表一个恒温源,具体做法如下：H=∑miv2i/2+V(r)+Qζ2/2+gKTlnS（10）式中：S和ζ分别是假想项的坐标和动量。这样体系的微分方程就变为：d/diivrr/=，()=-d/d+/iiiiaVrmvzm//，()2Bd/d=-/iiztmvgKTQ//∑（11）式中：g—体系的自由度；Q—一个可调参量,表征着假想项的质量,T为温度。4.2压力控制分子动力学模拟通常要求模拟体系的压力保持不变，而体积可以随温度的变化而改变，即进行等压模拟。这可以通过改变模拟原胞的3个方向或一个方向的尺寸来实现体积的变化。通常有：4.2.1Berendsen方法[31]这种控压方法的基本思路与基本控制方程均与他的控温方法类似,如下:μ=[1+△t(P-P0)γ/Г]1/3式中γ是一个可调参数,P0为所控制的压力。然后,在每一次的分子动力学模拟迭代中,粒子坐标x,y,z均用μ相乘,得到新的坐标,即可实现对压力P的控制。4.2.2Andersen方法[32]这种方法是通过改变体系的Lagrangian来实现的,新的Lagrangian定义如下：()'r11/2/2NNiiijijLmsGhVrWThhp=//=−+−Ω∑∑∑iiiii（12）其中：h=﹛a,b,c﹜—模拟单胞的基矢；G=h’h,原子坐标ri为ri=hs；V—原子之间的相互作用势；Ω—单胞体积；Tr—矩阵的迹；W—一个可调参量这种方法的优点是模拟单胞不仅大小可变,形状也可变。5分子模拟在材料科学中的应用分子动力学模拟计算液态金属的结构及热力学性质在国外得到了广泛地应用。其中波兰的Janusz等[33]利用多体相互作用，借助分子动力学计算了Ag、Au、Cu、Ni等面心立方金属的热力学性质。Holzman等[34]采用EAM势计算了面心立方金属液-气界面的特殊性