分子动力学中势函数研究

1、󰀁第20卷第2期2004年6月天󰀁津󰀁理󰀁工󰀁学󰀁院󰀁学󰀁报󰀁Jun.2004󰀁󰀁󰀁󰀁文章编号:1004󰀁2261(2004)02󰀁0101󰀁05分子动力学中势函数研究󰀁陈󰀁强1,曹红红1,黄海波1,2(1.北京航空航天大学理学院,北京100083;2.天津理工学院基础教育学院,天津300400)Aresear

2、chontheinteratomicpotentialinmoleculardynamics(MD)CHENQiang1,CAOHong󰀁hong1,HUANGHai󰀁bo1,2(1.InstituteofScience,BeihangUniversity,Beijing100083,China;2.FoundamentalEducationCollage,TianjinInstituteofTechnology,Tianjin300400,China)Abstract:Selectingappropriateinteratomicpotentialtodesc

3、ribetheinteratomicinteractionisthebaseofmolec󰀁ulardynamics(MD)computersimulation.Inthispaper,manyrepresentativeinteratomicpotentialshavebeenreviewed.Wegiveareviewofthemabouttheoreticbackground,modality,specialtyandapplicationarea.Keywords:moleculardynamics;computersimulation;interatomicpoten

4、tial(1)其中,V0反映了相互作用的强度;r0反映了原子的大小.根据量子力学二次微扰论的偶极子󰀁偶极子相互作用可导出n=12,这一项描述了范德瓦耳斯力.后一项是排斥力,其来源之一是原子核之间的库仑斥力,来源之二是电子之间由于泡利不相容原理产生的交叠能.󰀁󰀁在实际应用之中Lennard󰀁Jones势常取󰀁󰀁󰀁V(r)=V0(r0/r)12-2(r0/r)6(2)󰀁󰀁虽然上式是Lennard󰀁Jones势最为常用的形式,但也经常采用其他

5、形式的排斥项,特别是指数形式󰀁exp(-R/󰀁).用指数排斥项和倒六次方吸引项构成的势,通常称为Buckingham势.󰀁󰀁(2)式形式的Lennard󰀁Jones势常用来模拟惰性气体的凝聚态.对于分子晶体,例如CO2等可以采用原子-原子之间势能的加和来表示分子之间的相互作1󰀁对势模型󰀁收稿日期:2003󰀁10󰀁17󰀁󰀁陈󰀁󰀁󰀂󰀂),󰀁102

6、󰀁天󰀁津󰀁理󰀁工󰀁学󰀁院󰀁学󰀁报󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁第20卷󰀁第2期󰀁󰀁󰀁1929年,Morse注意到双原子分子的振动谱的量子力学问题可用指数形式的势函数解析地解决,并发现计算

7、结果与实验一致函数:1󰀁󰀁Daw和Baskes在准原子近似和有效介质理论的基础上根据密度泛函理论认为:某原子的原子核除了受到周围其他原子核的排斥作用外,还受到该原子的核外电子及其周围其他原子产生的背景电子的静电作用.于是,在EAM模型中,由N个原子组成的系统总能量可表示为:󰀁󰀁󰀁Etot=󰀁ij(rij)+󰀁Fi(󰀁(6)i)󰀁1j=󰀁2i=1󰀁ii=1其中󰀁ij是原子芯之间的两体作用对势,Fi是把原子i

8、埋入电子密度为󰀁i的位置处所需要的能量.其中󰀁i为系统中所有其他原子在i原子处产生的局域背景电子密度,可见这一项描述了势能的多体性质.󰀁󰀁EAM势具体表达形式是通过拟和完整纯金属体系的点阵常数、升华能、弹性常数、空位形成能、合金潜热等基本物理量来确定.早期势函数的拟和过程比较复杂,后来Foiles等人建立了一种比较简捷的拟和方法,即利用Rose普适函数:󰀁󰀁󰀁E(r)=-Ec1+b(-1)exp-b(-1)rere(7)󰀁eB)、󰀁e为平衡原子体积、E

9、c是结合能、EcB是弹性模量.令两式相等,则当电子密度函数󰀁i已其中:b=(知时,如果知道对势形式,可确定埋入势;反之知道埋入势,也可唯一确定对势.在对势部分,Daw和Foiles在其工作中采用库仑力形式,Johnson采用Born󰀁Mayer势,Vetor等人采用Morse势.早期,电子密度函数󰀁󰀁Fock理论计算的自由原子电子i一般用Hartree密度表示,近来则常常用负指数形式.Johnson等人在处理纯铁的󰀁-󰀁相界面时,则采用了分段函数形式描述电子密度.󰀁󰀁

10、EAM势很好的描述了金属原子之间的相互作用,是描述金属体系最常用的一种势函数.2.2󰀁Finnis󰀁Sinclair势󰀁󰀁F󰀁S势和EAM势的理论基础不同,它是建立在紧束缚(TB)近似的基础上的.F󰀁S势函数可以写成如下形式:󰀁󰀁󰀁E=i=1NNN314.于是他提出如下形式的势󰀁󰀁󰀁V(r)=V0exp-2󰀁(r/r0-1)-2exp-󰀁(r/r0-1)󰀁&

11、#983041;Born和Mayer1估计碱金属离子之间的排斥项可用指数形式表示,于是提出如下形式的势函数:ZiZje2CijDij󰀁i+󰀁j-r󰀁󰀁󰀁󰀁ij=+Aijbexp-6-8r󰀁rr(4)Born󰀁Mayer排斥项等价于Morse势的第一项.其中第󰀁󰀁利用对势模型人们能很好地拟合fcc金属的弹性常数,但用于bcc金属却比较失败.此时,C12=C44,而实际上Cauchy压力P=(C12-C44)/2却往往不等于零.后来,Joh

12、nson在对势项的基础上添加一个体积相关项2:󰀁󰀁󰀁Etot=NUV(V󰀁)+󰀁󰀁2i=1󰀁jNj=1󰀁i(3)󰀁󰀁Morse势和Lennard󰀁Jones势的曲线形式非常相󰀁N󰀁V(rij)(5)󰀁󰀁其中第二项是对势项,UV是与原子体积有关的能量项,V󰀁是平均原子体积.这一想法与用赝势对简单s-p键金属离子芯的分析比较符合.然而,添加体积依

13、赖项同样会导致一些矛盾.例如,除非体积依赖项是线性的,否则采用长波近似(连续介质)模型算出的体弹性模量和采用均匀形变方法计算的结果将会不同.22󰀁无方向性多体对泛函势󰀁󰀁以上对势模型在分子晶体和离子型化合物的模拟计算之中取得了比较大的成功,但是对于过渡金属,由于它们在金属键中含有一定的共价键,所以遇到了许20󰀁N(j󰀁i󰀁Vij(Rij)-F(󰀁i)(8)󰀁2004年6月󰀁󰀁󰀁󰀁󰀁b

14、3041;󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁陈󰀁强,等:分子动力学中势函数研究N󰀁103󰀁󰀁󰀁󰀁󰀁i=j=1󰀁󰀁ij(Rij)(9)󰀁󰀁MEAM势的角度因子是通过电子密度函数引进的,在位置i处由近邻原子产生的电子密度󰀁i为:2t()(12)󰀁h=1󰀁(h)其中

15、83041;(h=0,1,2,3)是分电子密度,分别对应核外󰀁󰀁󰀁󰀁i=󰀁exp2(0)(h)󰀁󰀁V和󰀁均为经验拟合的对势函数,公式8中第一项表示核-核之间的相互作用,第二项不再是嵌入能,而表示多体相互作用的关联能,󰀁可以理解为二分量紧束缚近似下的积分平方求和.Finnis󰀁Sinclair模型最初用来描述具有bcc结构的纯金属的原子间互作用,Ackland等人将它推广到二元合金体系5.对于二元合金体系,可以用下面6个对势函数表示VAA

16、、VBB、VAB、󰀁AA、󰀁BB、󰀁AB.Ackland等假定󰀁AA、󰀁BB、VAA、VBB与浓度无关,因此可以认为VAA、VBB、󰀁AA、󰀁BB和纯金属的参数相同;函数󰀁AB取󰀁AA和󰀁BB的几何平均,这和紧束缚近似下的积分求和的解释是一致的.因此在已知纯金属势参数的情况下,仅有对势项VAB需要通过拟合而得到5.󰀁󰀁Ackland给出上述模型的函数:VAA(Rij)=󰀁AA(Rij)=V

17、AB(Rij)=󰀁AB(Rij)=K=1(h)s,p,d,f轨道上的电子,t是分电子密度的权重因子.从理论上讲,s,p,d,f电子密度分别是晶体的体积、极化、切变和缺乏反演对称的量度.分电子密度由下式给出:(0)󰀁󰀁󰀁󰀁=(󰀁)=(2)2(󰀁)=(1)2(h)󰀁fuu,v(0)(ru)(ru)f(1)󰀁󰀁fi(1)(2)jxuxv(rv)rurvii󰀁󰀁fi,ju,v(ru)f(2)(rv)uuvv2-3

18、(rurv)(󰀁)=(3)2iji󰀁f(2)(ru)f(2)(rv)u,v󰀁66AAAAA3aAKH(rK-Rij)(rK-Rij)AAAAAA󰀁f(3)(ru)f(3)ijkxiuxjuxkuxvxvxv(rv)(rurv)3(13)󰀁󰀁󰀁K=1󰀁󰀁󰀁󰀁上面各式中,u,v是近邻原子的序号,xiu,xju,xku是u原子相对于计算原子的坐标分量.其中:-1)(14)re󰀁󰀁式中re是平衡原

19、子间距,该函数描述径向电子密󰀁󰀁󰀁f(h)(r)=exp-󰀁(h)(度.从上面可以看出各个原子对计算原子位置处的电子密度的贡献是互相影响的.󰀁󰀁MEAM势更精确地描述了原子间的相互作用,但是拟合复杂.由于这个势函数以及得出的力函数非常复杂,所以计算量比较大,根据运行时的比较,同样的原子数目,MEAM势的程序计算量大概是F󰀁S势的5倍左右.同时由于势函数很精细,所以程序不容易收敛.3.2󰀁键级势(BondOrderPotentials)󰀁󰀁

20、;1985年Abell11从化学赝势理论和定域轨道基集出发的每个原子结合能的解析表达式.Tersoff根据Abell理论给出了一个描述共价键相互作用的原子间互作用势模型,如下所示:󰀁󰀁󰀁E=12AKH(RK-Rij)(RK-Rij)ABAB3aABKH(rK-Rij)(rK-Rij)K=1󰀁AA(Rij)󰀁BB(Rij)(10)󰀁󰀁H(x)是阶跃函数.与EAM相比,F󰀁S多体势形式较为简单,势参数易于拟合.此外计算量比较小,程序容易收敛,所以应用很广泛.3♦

21、41;考虑角度效应的多体势󰀁󰀁EAM势还是一个中心力场势,对非立方系的金属间化合物,应用效果很差.Baskes等人为此发展了包含角度因子的MEAM势8,9,10,MEAM势是EAM势的扩展,其计算原子体系的能量公式和式(6)一样:󰀁󰀁󰀁Etot=2i=1j=1󰀁i󰀁Ei=iVij=fc(rij)󰀁2i,󰀁2i,j󰀁ij󰀁iAijexp(-󰀁1rij)-Bijexp(-󰀁2rij)(15)&

22、#983041;󰀁其中第一项是排斥项,被归因于正交化;第二项代表键能.其中Bij由(16)式给出:󰀁󰀁NN󰀁ij(rij)+i=1󰀁NFi(󰀁i)(11)󰀁104󰀁Bij=B0exp(-zij/b)󰀁󰀁󰀁zij=k󰀁i,j天󰀁津󰀁理󰀁工󰀁学󰀁院󰀁学󰀁报󰀁󰀁b

23、3041;󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁第20卷󰀁第2期󰀁非调和项.一般来说,引入的各种作用势成分越多,势n󰀁w(rik)/w(rij)󰀁-1(16)越精确.󰀁󰀁󰀁󰀁=󰀁󰀁󰀁c+exp(-dcos󰀁ijk)w(r)=fc(r)exp(-&#

24、983041;2r)bonds󰀁kb(b-bo)2+angles󰀁k󰀁(󰀁-2󰀁0)+󰀁󰀁Abell󰀁Tersoff势可以用于模拟半导体,有机物等13许多共价键化合物.后来BrennerDW对Abell󰀁Tersoff势进行了一定的修改,将它应用于炭氢系分子。3.3󰀁Chelikowsky󰀁Phillips势󰀁󰀁ChelikowskyJR和PhillipsJC根据第一性原理求出了具有共价键

25、结合的原子体系的势,其函数形式类似于前述的Abell󰀁Tersoff势.它不仅可以处理体相,而且还能处理表面,团簇状态那样含有键和强度由悬挂键转变成反键之类的量子效应.其函数形式可以表示为14+󰀁)+󰀁󰀁nbondsZiZj1264(r)-2(r)+rk󰀁1+cos(23)󰀁󰀁非键和作用势由Lennard󰀁Jones势和库仑势构成.分子间相互作用仅有非键和作用,包括库仑势、Lennard󰀁Jones势,但有时候也需要引入氢键。3.6󰀁紧

26、束缚(TB)势模型󰀁󰀁紧束缚(Tight󰀁Binding)方法是能够把量子力学引入原子间作用势计算的最简单可行的方法16,14,它将势函数建立在可靠的理论基础上,目前关于原子间势函数的工作很多都是以它为基础的,可以参看相关文献14,17.紧束缚方法是采用弱重叠的原子基函数的线性组合做为基函数或采用平面波基集,这是紧束缚方法的基本思想.紧束缚方法能够描述非定域态,所以可以用来模拟能带结构.此外紧束缚方法还同时考虑了键的方向性。󰀁󰀁采用原子轨道的线性组合,则N个原子组成的系统的波函数可以写成:󰀁定谔方程

27、得:󰀁󰀁󰀁Hi=1N:i<j󰀁󰀁󰀁󰀁=󰀁Aexp(-󰀁1r2gijexp(-󰀁2r2ij)ij)-2rijrij(17)󰀁󰀁gij是表示含有多体效应之键和强度的项,由下面给出:󰀁󰀁󰀁gij=g0(1+󰀁g0ij)+g1(1-󰀁g1ij)sijsji󰀁󰀁󰀁sij=1+b

28、3041;cos(3󰀁ijk)󰀁3.4󰀁Stillinger󰀁Waber势(S󰀁W势)󰀁󰀁对于Si和Ge等半导体,其键和强度依赖周围原子的配置,必须引入三体力及以上的所谓多体力,S󰀁W势是由StillingerFH和WaberTA在1985年引入的13,其表达式如下:󰀁󰀁󰀁󰀁=i<j(n)=Ni=1󰀁Nc(in)󰀁i.将其带入薛(18)(19)󰀁c(i

29、n)󰀁i=E(n)i=1󰀁n)c(i󰀁i(24)其中H为系统的总哈密顿量,由正交归一化条件可得N个本征值方程:󰀁󰀁󰀁󰀁HijcjjN(n)=E(n)c(in)(25)󰀁󰀁2(rij)+󰀁i<j<k󰀁3(ri,rj,rk)(20)(21)󰀁󰀁根据泡利不相容性原理,如果系统含有2z个电子,能带的能量可以写成:󰀁󰀁󰀁Eband=2n=1

30、rij󰀁󰀁󰀁󰀁2(rij)=󰀁f()󰀁zE(n)=zrirjrk󰀁󰀁󰀁󰀁3(ri,rj,rk)=󰀁g(,)(22)󰀁󰀁式中󰀁,󰀁分别表示与键和强度和原子大小相关联的参数,f和g两个函数的具体形式可以参看文献15.S󰀁W势是保证正四面体结构最稳定的势函数,一般只用与描述Si和Ge等半导体材料.3.5󰀁有机分子中的作用力问题

31、83041;󰀁有机分子内部的相互作用有键和作用和非键和作用,如式(23)所示4,14.其中,键和作用势由键长势、键角势和扭力势三部分组成.此外还应含有平面外角󰀁2󰀁ijn=1n)(n)c(cjHji=i󰀁󰀁󰀁ijHjiij(26)其中,󰀁ij是密度矩阵元,这些矩阵元表示了电荷密度的系数.其非对角矩阵元被称为键级,由此可以得出成键态和反键态上的电子数目:󰀁󰀁󰀁󰀁ij+󰀁ji=nbond-nantibond=

32、2󰀁(󰀁Eband=kn)(n)c(cii󰀁zn=1+n)(n)c(cj)i󰀁(27)根据Hellman󰀁Feynman定理:󰀁󰀁󰀁󰀁ij(󰀁ij󰀁Hji󰀁Hij+󰀁)jikk(28)󰀁2004年6月󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁𘀀

33、1;󰀁󰀁󰀁󰀁陈󰀁强,等:分子动力学中势函数研究󰀁105󰀁󰀁󰀁由此可以计算作用在每个原子上的合力.关于紧束缚方法可以参看文献4。󰀁󰀁除此以外,人们还发展了许多其它含角度的原子间作用势:例如对于EAM势模型,Pasianot和Farkas1也提出了另一种处理角度依赖关系的改进方法。还有Pearson,Biswas,Hamann提出了一个描述三体作用的比较一般的势函数,不过参数较多.总的来看目前所提出的考虑角度效应的多体势模型尚不很完善,而且往往比较复杂不便于运算.参󰀁考♦