分子動力學(xué)中勢函數(shù)研究

1、󰀁第20卷第2期2004年6月天󰀁津󰀁理󰀁工󰀁學(xué)󰀁院󰀁學(xué)󰀁報󰀁Jun.2004󰀁󰀁󰀁󰀁文章編號:1004󰀁2261(2004)02󰀁0101󰀁05分子動力學(xué)中勢函數(shù)研究󰀁陳󰀁強1,曹紅紅1,黃海波1,2(1.北京航空航天大學(xué)理學(xué)院,北京100083;2.天津理工學(xué)院基礎(chǔ)教育學(xué)院,天津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反映了原子的大小.根據(jù)量子力學(xué)二次微擾論的偶極子󰀁偶極子相互作用可導(dǎo)出n=12,這一項描述了范德瓦耳斯力.后一項是排斥力,其來源之一是原子核之間的庫侖斥力,來源之二是電子之間由于泡利不相容原理產(chǎn)生的交疊能.󰀁󰀁在實際應(yīng)用之中Lennard󰀁Jones勢常取󰀁󰀁󰀁V(r)=V0(r0/r)12-2(r0/r)6(2)󰀁󰀁雖然上式是Lennard󰀁Jones勢最為常用的形式,但也經(jīng)常采用其他

5、形式的排斥項,特別是指數(shù)形式󰀁exp(-R/󰀁).用指數(shù)排斥項和倒六次方吸引項構(gòu)成的勢,通常稱為Buckingham勢.󰀁󰀁(2)式形式的Lennard󰀁Jones勢常用來模擬惰性氣體的凝聚態(tài).對于分子晶體,例如CO2等可以采用原子-原子之間勢能的加和來表示分子之間的相互作1󰀁對勢模型󰀁收稿日期:2003󰀁10󰀁17󰀁󰀁陳󰀁󰀁󰀂󰀂),󰀁102

6、󰀁天󰀁津󰀁理󰀁工󰀁學(xué)󰀁院󰀁學(xué)󰀁報󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁第20卷󰀁第2期󰀁󰀁󰀁1929年,Morse注意到雙原子分子的振動譜的量子力學(xué)問題可用指數(shù)形式的勢函數(shù)解析地解決,并發(fā)現(xiàn)計算

7、結(jié)果與實驗一致函數(shù):1󰀁󰀁Daw和Baskes在準(zhǔn)原子近似和有效介質(zhì)理論的基礎(chǔ)上根據(jù)密度泛函理論認(rèn)為:某原子的原子核除了受到周圍其他原子核的排斥作用外,還受到該原子的核外電子及其周圍其他原子產(chǎn)生的背景電子的靜電作用.于是,在EAM模型中,由N個原子組成的系統(tǒng)總能量可表示為:󰀁󰀁󰀁Etot=󰀁ij(rij)+󰀁Fi(󰀁(6)i)󰀁1j=󰀁2i=1󰀁ii=1其中󰀁ij是原子芯之間的兩體作用對勢,Fi是把原子i

8、埋入電子密度為󰀁i的位置處所需要的能量.其中󰀁i為系統(tǒng)中所有其他原子在i原子處產(chǎn)生的局域背景電子密度,可見這一項描述了勢能的多體性質(zhì).󰀁󰀁EAM勢具體表達形式是通過擬和完整純金屬體系的點陣常數(shù)、升華能、彈性常數(shù)、空位形成能、合金潛熱等基本物理量來確定.早期勢函數(shù)的擬和過程比較復(fù)雜,后來Foiles等人建立了一種比較簡捷的擬和方法,即利用Rose普適函數(shù):󰀁󰀁󰀁E(r)=-Ec1+b(-1)exp-b(-1)rere(7)󰀁eB)、󰀁e為平衡原子體積、E

9、c是結(jié)合能、EcB是彈性模量.令兩式相等,則當(dāng)電子密度函數(shù)󰀁i已其中:b=(知時,如果知道對勢形式,可確定埋入勢;反之知道埋入勢,也可唯一確定對勢.在對勢部分,Daw和Foiles在其工作中采用庫侖力形式,Johnson采用Born󰀁Mayer勢,Vetor等人采用Morse勢.早期,電子密度函數(shù)󰀁󰀁Fock理論計算的自由原子電子i一般用Hartree密度表示,近來則常常用負(fù)指數(shù)形式.Johnson等人在處理純鐵的󰀁-󰀁相界面時,則采用了分段函數(shù)形式描述電子密度.󰀁󰀁

10、EAM勢很好的描述了金屬原子之間的相互作用,是描述金屬體系最常用的一種勢函數(shù).2.2󰀁Finnis󰀁Sinclair勢󰀁󰀁F󰀁S勢和EAM勢的理論基礎(chǔ)不同,它是建立在緊束縛(TB)近似的基礎(chǔ)上的.F󰀁S勢函數(shù)可以寫成如下形式:󰀁󰀁󰀁E=i=1NNN314.于是他提出如下形式的勢󰀁󰀁󰀁V(r)=V0exp-2󰀁(r/r0-1)-2exp-󰀁(r/r0-1)󰀁&

11、#983041;Born和Mayer1估計堿金屬離子之間的排斥項可用指數(shù)形式表示,于是提出如下形式的勢函數(shù):ZiZje2CijDij󰀁i+󰀁j-r󰀁󰀁󰀁󰀁ij=+Aijbexp-6-8r󰀁rr(4)Born󰀁Mayer排斥項等價于Morse勢的第一項.其中第󰀁󰀁利用對勢模型人們能很好地擬合fcc金屬的彈性常數(shù),但用于bcc金屬卻比較失敗.此時,C12=C44,而實際上Cauchy壓力P=(C12-C44)/2卻往往不等于零.后來,Joh

12、nson在對勢項的基礎(chǔ)上添加一個體積相關(guān)項2:󰀁󰀁󰀁Etot=NUV(V󰀁)+󰀁󰀁2i=1󰀁jNj=1󰀁i(3)󰀁󰀁Morse勢和Lennard󰀁Jones勢的曲線形式非常相󰀁N󰀁V(rij)(5)󰀁󰀁其中第二項是對勢項,UV是與原子體積有關(guān)的能量項,V󰀁是平均原子體積.這一想法與用贗勢對簡單s-p鍵金屬離子芯的分析比較符合.然而,添加體積依

13、賴項同樣會導(dǎo)致一些矛盾.例如,除非體積依賴項是線性的,否則采用長波近似(連續(xù)介質(zhì))模型算出的體彈性模量和采用均勻形變方法計算的結(jié)果將會不同.22󰀁無方向性多體對泛函勢󰀁󰀁以上對勢模型在分子晶體和離子型化合物的模擬計算之中取得了比較大的成功,但是對于過渡金屬,由于它們在金屬鍵中含有一定的共價鍵,所以遇到了許20󰀁N(j󰀁i󰀁Vij(Rij)-F(󰀁i)(8)󰀁2004年6月󰀁󰀁󰀁󰀁󰀁b

14、3041;󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁陳󰀁強,等:分子動力學(xué)中勢函數(shù)研究N󰀁103󰀁󰀁󰀁󰀁󰀁i=j=1󰀁󰀁ij(Rij)(9)󰀁󰀁MEAM勢的角度因子是通過電子密度函數(shù)引進的,在位置i處由近鄰原子產(chǎn)生的電子密度󰀁i為:2t()(12)󰀁h=1󰀁(h)其中

15、83041;(h=0,1,2,3)是分電子密度,分別對應(yīng)核外󰀁󰀁󰀁󰀁i=󰀁exp2(0)(h)󰀁󰀁V和󰀁均為經(jīng)驗擬合的對勢函數(shù),公式8中第一項表示核-核之間的相互作用,第二項不再是嵌入能,而表示多體相互作用的關(guān)聯(lián)能,󰀁可以理解為二分量緊束縛近似下的積分平方求和.Finnis󰀁Sinclair模型最初用來描述具有bcc結(jié)構(gòu)的純金屬的原子間互作用,Ackland等人將它推廣到二元合金體系5.對于二元合金體系,可以用下面6個對勢函數(shù)表示VAA

16、、VBB、VAB、󰀁AA、󰀁BB、󰀁AB.Ackland等假定󰀁AA、󰀁BB、VAA、VBB與濃度無關(guān),因此可以認(rèn)為VAA、VBB、󰀁AA、󰀁BB和純金屬的參數(shù)相同;函數(shù)󰀁AB取󰀁AA和󰀁BB的幾何平均,這和緊束縛近似下的積分求和的解釋是一致的.因此在已知純金屬勢參數(shù)的情況下,僅有對勢項VAB需要通過擬合而得到5.󰀁󰀁Ackland給出上述模型的函數(shù):VAA(Rij)=󰀁AA(Rij)=V

17、AB(Rij)=󰀁AB(Rij)=K=1(h)s,p,d,f軌道上的電子,t是分電子密度的權(quán)重因子.從理論上講,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原子相對于計算原子的坐標(biāo)分量.其中:-1)(14)re󰀁󰀁式中re是平衡原

19、子間距,該函數(shù)描述徑向電子密󰀁󰀁󰀁f(h)(r)=exp-󰀁(h)(度.從上面可以看出各個原子對計算原子位置處的電子密度的貢獻是互相影響的.󰀁󰀁MEAM勢更精確地描述了原子間的相互作用,但是擬合復(fù)雜.由于這個勢函數(shù)以及得出的力函數(shù)非常復(fù)雜,所以計算量比較大,根據(jù)運行時的比較,同樣的原子數(shù)目,MEAM勢的程序計算量大概是F󰀁S勢的5倍左右.同時由于勢函數(shù)很精細(xì),所以程序不容易收斂.3.2󰀁鍵級勢(BondOrderPotentials)󰀁󰀁

20、;1985年Abell11從化學(xué)贗勢理論和定域軌道基集出發(fā)的每個原子結(jié)合能的解析表達式.Tersoff根據(jù)Abell理論給出了一個描述共價鍵相互作用的原子間互作用勢模型,如下所示:󰀁󰀁󰀁E=12AKH(RK-Rij)(RK-Rij)ABAB3aABKH(rK-Rij)(rK-Rij)K=1󰀁AA(Rij)󰀁BB(Rij)(10)󰀁󰀁H(x)是階躍函數(shù).與EAM相比,F󰀁S多體勢形式較為簡單,勢參數(shù)易于擬合.此外計算量比較小,程序容易收斂,所以應(yīng)用很廣泛.3♦

21、41;考慮角度效應(yīng)的多體勢󰀁󰀁EAM勢還是一個中心力場勢,對非立方系的金屬間化合物,應(yīng)用效果很差.Baskes等人為此發(fā)展了包含角度因子的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天󰀁津󰀁理󰀁工󰀁學(xué)󰀁院󰀁學(xué)󰀁報󰀁󰀁b

23、3041;󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁第20卷󰀁第2期󰀁非調(diào)和項.一般來說,引入的各種作用勢成分越多,勢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勢可以用于模擬半導(dǎo)體,有機物等13許多共價鍵化合物.后來BrennerDW對Abell󰀁Tersoff勢進行了一定的修改,將它應(yīng)用于炭氫系分子。3.3󰀁Chelikowsky󰀁Phillips勢󰀁󰀁ChelikowskyJR和PhillipsJC根據(jù)第一性原理求出了具有共價鍵

25、結(jié)合的原子體系的勢,其函數(shù)形式類似于前述的Abell󰀁Tersoff勢.它不僅可以處理體相,而且還能處理表面,團簇狀態(tài)那樣含有鍵和強度由懸掛鍵轉(zhuǎn)變成反鍵之類的量子效應(yīng).其函數(shù)形式可以表示為14+󰀁)+󰀁󰀁nbondsZiZj1264(r)-2(r)+rk󰀁1+cos(23)󰀁󰀁非鍵和作用勢由Lennard󰀁Jones勢和庫侖勢構(gòu)成.分子間相互作用僅有非鍵和作用,包括庫侖勢、Lennard󰀁Jones勢,但有時候也需要引入氫鍵。3.6󰀁緊

26、束縛(TB)勢模型󰀁󰀁緊束縛(Tight󰀁Binding)方法是能夠把量子力學(xué)引入原子間作用勢計算的最簡單可行的方法16,14,它將勢函數(shù)建立在可靠的理論基礎(chǔ)上,目前關(guān)于原子間勢函數(shù)的工作很多都是以它為基礎(chǔ)的,可以參看相關(guān)文獻14,17.緊束縛方法是采用弱重疊的原子基函數(shù)的線性組合做為基函數(shù)或采用平面波基集,這是緊束縛方法的基本思想.緊束縛方法能夠描述非定域態(tài),所以可以用來模擬能帶結(jié)構(gòu).此外緊束縛方法還同時考慮了鍵的方向性。󰀁󰀁采用原子軌道的線性組合,則N個原子組成的系統(tǒng)的波函數(shù)可以寫成:󰀁定諤方程

27、得:󰀁󰀁󰀁Hi=1N:i<j󰀁󰀁󰀁󰀁=󰀁Aexp(-󰀁1r2gijexp(-󰀁2r2ij)ij)-2rijrij(17)󰀁󰀁gij是表示含有多體效應(yīng)之鍵和強度的項,由下面給出:󰀁󰀁󰀁gij=g0(1+󰀁g0ij)+g1(1-󰀁g1ij)sijsji󰀁󰀁󰀁sij=1+b

28、3041;cos(3󰀁ijk)󰀁3.4󰀁Stillinger󰀁Waber勢(S󰀁W勢)󰀁󰀁對于Si和Ge等半導(dǎo)體,其鍵和強度依賴周圍原子的配置,必須引入三體力及以上的所謂多體力,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為系統(tǒng)的總哈密頓量,由正交歸一化條件可得N個本征值方程:󰀁󰀁󰀁󰀁HijcjjN(n)=E(n)c(in)(25)󰀁󰀁2(rij)+󰀁i<j<k󰀁3(ri,rj,rk)(20)(21)󰀁󰀁根據(jù)泡利不相容性原理,如果系統(tǒng)含有2z個電子,能帶的能量可以寫成:󰀁󰀁󰀁Eband=2n=1

30、rij󰀁󰀁󰀁󰀁2(rij)=󰀁f()󰀁zE(n)=zrirjrk󰀁󰀁󰀁󰀁3(ri,rj,rk)=󰀁g(,)(22)󰀁󰀁式中󰀁,󰀁分別表示與鍵和強度和原子大小相關(guān)聯(lián)的參數(shù),f和g兩個函數(shù)的具體形式可以參看文獻15.S󰀁W勢是保證正四面體結(jié)構(gòu)最穩(wěn)定的勢函數(shù),一般只用與描述Si和Ge等半導(dǎo)體材料.3.5󰀁有機分子中的作用力問題

31、83041;󰀁有機分子內(nèi)部的相互作用有鍵和作用和非鍵和作用,如式(23)所示4,14.其中,鍵和作用勢由鍵長勢、鍵角勢和扭力勢三部分組成.此外還應(yīng)含有平面外角󰀁2󰀁ijn=1n)(n)c(cjHji=i󰀁󰀁󰀁ijHjiij(26)其中,󰀁ij是密度矩陣元,這些矩陣元表示了電荷密度的系數(shù).其非對角矩陣元被稱為鍵級,由此可以得出成鍵態(tài)和反鍵態(tài)上的電子數(shù)目:󰀁󰀁󰀁󰀁ij+󰀁ji=nbond-nantibond=

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

33、1;󰀁󰀁󰀁󰀁陳󰀁強,等:分子動力學(xué)中勢函數(shù)研究󰀁105󰀁󰀁󰀁由此可以計算作用在每個原子上的合力.關(guān)于緊束縛方法可以參看文獻4。󰀁󰀁除此以外,人們還發(fā)展了許多其它含角度的原子間作用勢:例如對于EAM勢模型,Pasianot和Farkas1也提出了另一種處理角度依賴關(guān)系的改進方法。還有Pearson,Biswas,Hamann提出了一個描述三體作用的比較一般的勢函數(shù),不過參數(shù)較多.總的來看目前所提出的考慮角度效應(yīng)的多體勢模型尚不很完善,而且往往比較復(fù)雜不便于運算.參󰀁考♦