高壓下鉬的聲子色散的改進(jìn)分析型嵌入原子法研究*

張曉軍，王安祥 ，高 賓，陳長(zhǎng)樂(lè)

(1.西安工程大學(xué)理學(xué)院，陜西西安 710048; 2.西北工業(yè)大學(xué)陜西省凝聚態(tài)結(jié)構(gòu)與性質(zhì)重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072)

1 引 言

聲子色散關(guān)系是研究固體材料物理性質(zhì)如比熱、熱導(dǎo)、超導(dǎo)、聲速、聲-電相互作用的重要信息[1]。此外，高壓下聲子譜中的負(fù)頻模式是評(píng)判固體材料結(jié)構(gòu)相變的重要依據(jù)[2]。因此固體聲子色散關(guān)系的實(shí)驗(yàn)和理論研究一直吸引著眾多科研工作者[3-4]。對(duì)于金屬或合金在常壓下的聲子色散關(guān)系已存在大量的實(shí)驗(yàn)和理論研究結(jié)果。近年來(lái)，人們趨于對(duì)特殊條件下(如高溫和高壓)聲子色散的研究。實(shí)驗(yàn)方面，把金剛石壓砧技術(shù)和非彈性X射線散射技術(shù)或非彈性中子散射技術(shù)相結(jié)合來(lái)測(cè)量高壓下固體的聲子色散關(guān)系[5-6]。另一方面，試圖從理論上重現(xiàn)這些實(shí)驗(yàn)結(jié)果或預(yù)言高壓下的聲子色散關(guān)系。Fang等[7]利用密度泛函理論在準(zhǔn)諧和近似下計(jì)算了壓強(qiáng)分別為23和55 GPa時(shí)金屬銥沿對(duì)稱方向上的聲子譜。Hu等[8]利用第一原理計(jì)算了不同壓強(qiáng)下金屬鈦Ti沿4個(gè)對(duì)稱方向的聲子色散曲線。這些方法物理意義深刻、明確且能提供原子間相互作用的精確信息，在計(jì)算金屬和合金的聲子色散關(guān)系方面也非常成功。但是這些方法要求計(jì)算機(jī)有相當(dāng)大的計(jì)算能力，而且需要花費(fèi)相當(dāng)多的計(jì)算時(shí)間。同時(shí)，由于計(jì)算局限性，該方法最多只能處理上百個(gè)原子的系統(tǒng)，往往不能高效地描述實(shí)際的材料結(jié)構(gòu),尤其對(duì)于復(fù)雜的體系難以普遍應(yīng)用。由于這些原因，研究人員開始著手利用經(jīng)驗(yàn)多體勢(shì)模型來(lái)研究高壓下的聲子色散關(guān)系。與第一原理計(jì)算相比，經(jīng)驗(yàn)多體勢(shì)模型能夠快速提供原子或分子間的相互作用信息，且能處理幾千甚至上萬(wàn)個(gè)原子的龐雜系統(tǒng)。Kazanc等[9]曾利用經(jīng)驗(yàn)多體勢(shì)即Sutton-Chen式嵌入原子法(SCEAM)計(jì)算了不同壓強(qiáng)下鎳鈀無(wú)序合金沿[100]、[110]和[111] 3個(gè)高對(duì)稱方向上的聲子色散曲線。我們也曾用張邦維等發(fā)展的改進(jìn)分析型嵌入原子法(MAEAM)模擬了高壓下金屬銅沿3個(gè)高對(duì)稱方向和4個(gè)低對(duì)稱方向上的聲子色散曲線[10]。然而，這些結(jié)果都沒(méi)有相應(yīng)的實(shí)驗(yàn)數(shù)據(jù)相佐證，難于說(shuō)明經(jīng)驗(yàn)多體勢(shì)是否能合理描述高壓下原子間的相互作用。本研究把張邦維等發(fā)展的改進(jìn)分析型嵌入原子法模型與晶格動(dòng)力學(xué)理論相結(jié)合，模擬高壓下金屬鉬(Mo)沿3個(gè)高對(duì)稱方向上的聲子色散曲線，并和對(duì)應(yīng)的實(shí)驗(yàn)結(jié)果進(jìn)行比較。

2 改進(jìn)分析型嵌入原子法模型

在重現(xiàn)或預(yù)測(cè)固體材料聲子色散關(guān)系的計(jì)算中，原子之間的相互作用勢(shì)是必不可少的。經(jīng)常采用的兩體勢(shì)有Morse勢(shì)[11]和Murrell Sorbie勢(shì)[12]，多體勢(shì)模型有密度泛函理論(DFT)[13]、嵌入原子勢(shì)模型(EAM)[14]和緊束縛理論(TBT)[15]等。其中，由Zhang等[16-17]發(fā)展的改進(jìn)分析型嵌入原子法模型已成功地計(jì)算了面心立方金屬的界面能[18]、晶界能[19]以及體心立方金屬的空位遷移能[17]等。改進(jìn)分析型嵌入原子法的基本公式為[16-17]

式中:Et是系統(tǒng)的總能量，F(xiàn)(ρi)是在除第i個(gè)原子外的其他原子組成的基體中再嵌入第i個(gè)原子的嵌入能，僅是其他原子在第i個(gè)原子所在處產(chǎn)生的背景電子密度ρi的函數(shù)，其中f(rij)為單個(gè)孤立原子的球型電子密度分布函數(shù)，rij是第i個(gè)原子和第j個(gè)原子間的距離；φ(rij)是第i個(gè)原子和第j個(gè)原子間的相互作用能；M(Pi)是修正項(xiàng)，表示原子電子密度非球型對(duì)稱分布所引起的系統(tǒng)總能量的變化。嵌入函數(shù)F(ρi)、原子間的相互作用勢(shì)φ(rij)、修正項(xiàng)M(Pi)和電子密度函數(shù)f(rij)采用如下形式[16-17]

式中:下標(biāo)“e”表示平衡狀態(tài)，re表示在平衡狀態(tài)下純?cè)鼐w中原子的最近鄰距離，其他模型參數(shù)如F0、n、α、ks(s=-1～4)見表1。

表1 金屬鉬的模型參數(shù)Table 1 The model parameters of metal Mo

3 晶格振動(dòng)模型

圖1 體心立方晶體結(jié)構(gòu) Fig.1 Bcc crystal structure

常溫常壓下，金屬鉬屬體心立方結(jié)構(gòu)，其晶體結(jié)構(gòu)如圖1所示，晶格常數(shù)為a。在晶體結(jié)構(gòu)中設(shè)坐標(biāo)原點(diǎn)處的原子為參考原子n，根據(jù)晶格動(dòng)力學(xué)理論[20]，參考原子n的振動(dòng)方程為

(8)

(9)

求解(8)式，可得

(10)

(11)

顯然，(11)式是A(n)的3個(gè)線性齊次方程，因?yàn)榫Ц裾駝?dòng)總是存在的，故A(n)有非零解的條件是

|Dα β(n,m)-δα βω2(q)|=0

(12)

通過(guò)聯(lián)立(1)式、(9)和(12)式，就可以得到體系的振動(dòng)頻率ω(q)和波矢之間的關(guān)系，即聲子色散關(guān)系。

4 模擬方法和結(jié)果分析

在晶格動(dòng)力學(xué)理論和改進(jìn)分析型嵌入原子法的框架下，采用數(shù)值計(jì)算的方法模擬了高壓下金屬鉬的聲子色散曲線。在模擬過(guò)程中，選取8a×8a×8a的晶塊作為模擬單元，其中包括坐標(biāo)原點(diǎn)在內(nèi)的4a×4a×4a晶塊作為內(nèi)部單元，其余原子為外部單元，外部單元是為了確保內(nèi)部單元邊界處和近邊界處的原子周圍有足夠多的近鄰原子與之相互作用，同時(shí)假設(shè)所有的原子滿足三維周期邊界條件。選取對(duì)勢(shì)函數(shù)的截止距離為0.626 nm，應(yīng)用Matlab語(yǔ)言編輯計(jì)算機(jī)程序計(jì)算參考原子與其他原子之間的原子力常數(shù)和動(dòng)力學(xué)矩陣元，把這些值代入久期方程((12)式)獲得不同波矢所對(duì)應(yīng)的聲子頻率值，具體的計(jì)算細(xì)節(jié)可參考文獻(xiàn)[21]和文獻(xiàn)[22]。在計(jì)算不同壓強(qiáng)下的聲子色散關(guān)系時(shí)，采用Rose-Vinet[23-24]方程來(lái)控制加載在體系上的靜壓力。根據(jù)Rose-Vinet方程，金屬鉬的壓強(qiáng)和壓縮的關(guān)系可表示為

(13)

圖2 金屬鉬的壓強(qiáng)和體積間的關(guān)系 Fig.2 Pressure versus volume for Mo

由MAEAM模型所計(jì)算的壓強(qiáng)分別為0.1 MPa、17 GPa和37 GPa時(shí)金屬鉬沿[00ζ]、[0ζζ]和[ζζζ] 3個(gè)高對(duì)稱方向的聲子色散曲線如圖3、圖4、圖5所示，其中實(shí)線為模擬值，實(shí)心圓點(diǎn)為實(shí)驗(yàn)值[26-27]，群論符號(hào)Γ、H、P′、P和N表示體心立方晶格布里淵區(qū)中的不同對(duì)稱點(diǎn)，符號(hào)L和T分別表示縱振動(dòng)模和橫振動(dòng)模，ζ=q/qmax是簡(jiǎn)約波矢。根據(jù)晶格動(dòng)力學(xué)理論，體心立方布拉菲晶胞會(huì)產(chǎn)生3個(gè)聲子振動(dòng)支，其中2個(gè)為橫波聲學(xué)振動(dòng)支，1個(gè)是縱波聲學(xué)振動(dòng)支，由于四重和六重軸對(duì)稱性，2個(gè)橫波聲學(xué)振動(dòng)支往往會(huì)在某些振動(dòng)方向上兼并為1個(gè)橫波振動(dòng)支，如圖3、圖4、圖5所示，這些圖像和現(xiàn)在的計(jì)算結(jié)果相符合。從圖3可以看出，在常壓下(0.1 MPa)，我們的計(jì)算結(jié)果整體上和實(shí)驗(yàn)結(jié)果相一致，包括一些布里淵邊界點(diǎn)，例如T[00ζ]和L[00ζ] 振動(dòng)支中的H點(diǎn)，T1[0ζζ]和T2[0ζζ]振動(dòng)支中的N點(diǎn)。在L[00ζ]振動(dòng)支中簡(jiǎn)約波矢ζ在0.3～0.6之間，T[ζζζ]振動(dòng)支中ζ在0.3～0.9之間，L[0ζζ]振動(dòng)支中ζ在0.3～0.5之間。計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果在數(shù)值上有些偏差，但所對(duì)應(yīng)色散曲線的形狀卻非常相似(圖3)。壓強(qiáng)p=17 GPa時(shí)，色散曲線的理論計(jì)算實(shí)線和實(shí)驗(yàn)點(diǎn)線整體一致，特別在低頻附近(圖4)。與常壓下聲子色散曲線相比，在壓強(qiáng)為17 GPa時(shí)的色散曲線中，頻率偏差主要出現(xiàn)在T[00ζ]和L[00ζ]振動(dòng)支中的H點(diǎn)、T[ζζζ]和L[ζζζ]振動(dòng)支中的P點(diǎn)以及T1[0ζζ]和L[0ζζ]振動(dòng)支中的N點(diǎn)附近，最大偏差不超過(guò)10.8%。壓強(qiáng)p=37 GPa時(shí)，L[00ζ]振動(dòng)支的色散曲線和實(shí)驗(yàn)結(jié)果基本吻合，其他方向上色散曲線的形狀和常壓以及17 GPa時(shí)的結(jié)果非常相似，只是頻率略有上移，遺憾的是只有L[00ζ] 振動(dòng)支的實(shí)驗(yàn)結(jié)果與之相比較，而缺少其他振動(dòng)支的實(shí)驗(yàn)數(shù)據(jù)。從上面的分析可看出，我們的計(jì)算結(jié)果均與實(shí)驗(yàn)結(jié)果基本吻合，只是數(shù)值上有較些差異。考慮到計(jì)算結(jié)果是在諧和近似下所得，而實(shí)驗(yàn)結(jié)果并不能忽略非諧和效應(yīng)，因此二者之間存在差異是合理的。這表明MAEAM模型可以合理描述常壓和高壓下原子之間的相互作用。

圖3 常壓下鉬的聲子色散曲線 Fig.3 Phonon dispersion curves of Mo under normal pressures

圖5 壓強(qiáng)p=37 GPa時(shí)鉬的聲子色散曲線 Fig.5 Phonon dispersion curves of W at 37 GPa

圖6 不同高壓下鉬的聲子色散曲線 Fig.6 Phonon dispersion curves of Mo under different high pressures

圖6是壓強(qiáng)分別為60、80和100 GPa時(shí)的聲子色散曲線，其中固體實(shí)線是高壓下的色散曲線，為了方便比較，在圖6中我們還添加了常壓下的色散曲線，并用點(diǎn)線表示。從圖6中可以看出，常壓下以及高壓下(p=60、80和100 GPa)的聲子色散曲線非常相似，這種相似性并不驚奇，因?yàn)榻饘巽f保持體心立方結(jié)構(gòu)直到我們所考慮的最大壓強(qiáng)范圍內(nèi)(100 GPa)。高壓下所有振動(dòng)支的振動(dòng)頻率均高于常壓下的振動(dòng)頻率，與常壓聲子色散相比，p=60 GPa時(shí)的最大振動(dòng)頻率增加了14.3%，p=80 GPa時(shí)增加了16.3%，p=100 GPa時(shí)增加了18.1%。在所有的振動(dòng)支中的振動(dòng)頻率均隨壓強(qiáng)的增大而增大。在不考慮質(zhì)量、相互作用距離和結(jié)合能所引起的差別時(shí)，這些現(xiàn)象和高壓下體心立方金屬鎢[28]和鉭[29]的計(jì)算結(jié)果相一致。

為了近一步說(shuō)明聲子頻率和壓強(qiáng)之間的關(guān)系，模擬了金屬鉬的T2[0ζζ]和L[00ζ]振動(dòng)支中簡(jiǎn)約波矢分別為0.2、0.5、0.8 和1.0時(shí)的聲子色散頻率隨壓強(qiáng)的關(guān)系，如圖7所示。從圖7可以看出，在兩個(gè)振動(dòng)支上，4個(gè)簡(jiǎn)約點(diǎn)的聲子頻率均隨壓強(qiáng)的增加而增加。在相同振動(dòng)支中各簡(jiǎn)約點(diǎn)的聲子頻率和壓強(qiáng)之間的關(guān)系基本相同。因此，如果想知道聲子色散和壓強(qiáng)之間的關(guān)系時(shí)，我們只要考慮某一振動(dòng)方向上某一簡(jiǎn)約點(diǎn)的振動(dòng)情況即可。此外，在40～100 GPa范圍內(nèi)，各個(gè)簡(jiǎn)約點(diǎn)的頻率和壓強(qiáng)基本上呈線性關(guān)系。

圖7 不同簡(jiǎn)約波矢的聲子頻率和壓強(qiáng)之間的關(guān)系 Fig.7 Phonon frequencies versus pressures at different reduced wave vectors

5 結(jié) 論

應(yīng)用改進(jìn)分析型嵌入原子法模型計(jì)算了不同高壓下金屬鉬的原子力常數(shù)和動(dòng)力學(xué)矩陣，重現(xiàn)了3種壓強(qiáng)下金屬鉬沿3個(gè)高對(duì)稱方向上聲子色散的實(shí)驗(yàn)結(jié)果，預(yù)測(cè)了鉬在壓強(qiáng)分別為60、80和100 GPa時(shí)的聲子色散曲線。得到以下結(jié)果：(1) 壓強(qiáng)分別為0.1 MPa、17 GPa和37 GPa時(shí)金屬鉬的聲子色散曲線的模擬結(jié)果和實(shí)驗(yàn)值符合較好，特別在低頻附近二者幾乎一致。(2) 在壓強(qiáng)分別為p=60、80和100 GPa時(shí)所預(yù)測(cè)的聲子色散曲線和常壓下聲子色散曲線的形狀都非常相似，且隨著壓強(qiáng)的增大各振動(dòng)支的振動(dòng)頻率均依次增大。(3) 不同振動(dòng)支上不同簡(jiǎn)約點(diǎn)的聲子頻率均隨壓強(qiáng)的增加而增加，在40～100 GPa范圍內(nèi)，所考慮的各個(gè)簡(jiǎn)約點(diǎn)的頻率和壓強(qiáng)基本上呈線性關(guān)系。(4) 上述方法也可用于其他金屬或合金在高壓條件下聲子色散關(guān)系的計(jì)算和研究。

[1] XIE Y,ZHANG J M.Atomistic simulation of phonon dispersion for body-centred cubic alkali metals [J].Can J Phys,2008,86(6):801-805.

[2] WANG Y,HECTOR L G,ZHANG H,et al.Thermodynamics of the Ceγ-αtransition:Density functional study [J].Phys Rev B,2008,78(10):104113.

[3] WONG J,KRISCH M,FARBER D L,et al.Crystal dynamics of fcc Pu-Ga alloy by high-resolution inelastic X-ray scattering [J].Phys Rev B,2005,72(6):064115.

[4] ADITYA M V.Phonon dispersion in equiatomic Li-based binary alloys [J].Chin Phys Lett,2008,25(2):654-657.

[5] DANIEL L F,KRISCH M,ANTONANGELI D,et al.Lattice dynamics of molybdenum at high pressure [J].Phys Rev Lett,2006,96(4):115502.

[6] MITTAL R,CHAPLOT S L,CHOUDHURY N,et al.Inelastic neutron scattering,lattice dynamics and high-pressure phase stability in LuPO4and YbPO4[J].J Phys-Condens Mat,2007,19(44):6202.

[7] FANG H,LIU B,GU M.High-pressure lattice dynamic and thermodynamic properties of Ir by first-principles calculation [J].Physica B,2010,405(2):732-737.

[8] HU C E,ZENG Z Y,ZHANG L.Theoretical investigation of the high pressure structure,lattice dynamics,phase transition,and thermal equation of state of titanium metal [J].J Appl Phys,2010,107(9):093509.

[9] KAZANC S.The effects on the lattice dynamical properties of the temperature and pressure in random NiPd alloy [J].Can J Phys,2013,91(10):833-838.

[10] ZHANG X J,CHEN C L,FENG F L.High-pressure phonon dispersion of copper by using the modified analytic embedded atom method [J].Chin Phys B,2013,22(9):096301.

[11] MORSE P M.Diatomic molecules according to the wave mechanics.Ⅱ.Vibrational levels [J].Phys Rev,1929,34(1):57-64.

[12] MURRELL J N,SORBIE K S.New analytic form for potential energy curves of stable diatomic states [J].J Chem Soc Faraday Trans,1974,70(2):1552-1556.

[13] SEBETCI A.A density functional study of bare and hydrogenated platinum clusters [J].Chem Phys,2006,331(1):9-18.

[14] DAW M S,BASKES M I.Embedded atom method:A review derivation and application to impurities,surface and other defects in metals [J].Phys Rev B,1984,29(12):6443-6447.

[15] BOYKIN T B,VOGL P.Dielectric response of molecules in empirical tight-binding theory [J].Phys Rev B,2002,65(3):035202.

[16] OUYANG Y F,ZHANG B W.Analytic embedded atom potentials for bcc metals:application to calculating the thermodynamic data of bcc alloys [J].Phys Lett A,1994,192(1):79-86.

[17] HU W Y,SHU X L,ZHANG B W.Point-defect properties in body-centered cubic transition metals with analytic EAM interatomic potentials [J].Comput Mater Sci,2002,23(1):175-189.

[18] MA F,ZHANG J M,XU K W.Theoretical analysis of interface energy for unrelaxed Ag (001)/Ni (001) twist interface boundaries with MAEAM [J].Surf Interface Anal,2004,36(4):355-359.

[19] ZHANG J M,WEI X M,XIN H,et al.Atomic scale calculation of energies of Cu (001) twist boundaries [J].Chin Phys,2004,14(5):1015-1020.

[20] DOVE M T.Introducton to lattice dynamics [M].New York:Cambridge University Press,1993.

[21] KAZANC S,OZGEN S.Pressure effect on phonon frequencies in some transition metals:A molecular dynamics study [J].Physica B,2005,365(1/2/3/4):185-192.

[22] KAZANC S,OZGEN S,ADIGUZEL O.Pressure effects on martensitic transformation under quenching process in a molecular dynamics model of NiAl alloy [J].Physica B,2003,334(3/4):375-381.

[23] ROSE J H,SMITH J R,GUINEA F,et al.Universal features of the equation of state of metals [J].Phys Rev B,1984,29(6):2963-2969.

[24] VINET P,ROSE J H,FERRANTE J,et al.Universal features of the equation of state of solids [J].J Phys-Condens Mat,1989,1(11):1941-1963.

[25] MAO H K,BELL P M,SHANER J W,et al.Specific volume measurements of Cu,Mo,Pd,and Ag and calibration of the rubyR1fluorescence pressure gauge from 0.06 to 1 Mbar [J].J Appl Phys,1978,49(6):3276-3283.

[26] WOODS A D B,CHEN S H.Lattice dynamics of molybdenum [J].Solid State Commun,1964,2(8):233-242.

[27] FARBER D L,KRISCH M,ANTONANGELI D,et al.Lattice dynamics of molybdenum at high pressure [J].Phys Rev Lett,2006,96(11):115502.

[28] ZHANG X J,CHEN C L.Pressure dependence of phonon dispersion in bcc tungsten [J].Chin J Phys,2013,51(2):359-367.

[29] TAIOLI S,CAZORLA C,GILLAN M J,et al.Melting curve of tantalum from first principles [J].Phys Rev B,2007,75(21):214103.