張曉軍，王安祥 ，高 賓，陳長(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ì)稱(chēng)方向上的聲子譜。Hu等[8]利用第一原理計(jì)算了不同壓強(qiáng)下金屬鈦Ti沿4個(gè)對(duì)稱(chēng)方向的聲子色散曲線。這些方法物理意義深刻、明確且能提供原子間相互作用的精確信息，在計(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)用。由于這些原因，研究人員開(kāi)始著手利用經(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ì)稱(chēng)方向上的聲子色散曲線。我們也曾用張邦維等發(fā)展的改進(jìn)分析型嵌入原子法(MAEAM)模擬了高壓下金屬銅沿3個(gè)高對(duì)稱(chēng)方向和4個(gè)低對(duì)稱(chēng)方向上的聲子色散曲線[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ì)稱(chēng)方向上的聲子色散曲線，并和對(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ì)稱(chēng)分布所引起的系統(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)見(jiàn)表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)部單元邊界處和近邊界處的原子周?chē)凶銐蚨嗟慕徳优c之相互作用，同時(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ì)稱(chēng)方向的聲子色散曲線如圖3、圖4、圖5所示，其中實(shí)線為模擬值，實(shí)心圓點(diǎn)為實(shí)驗(yàn)值[26-27]，群論符號(hào)Γ、H、P′、P和N表示體心立方晶格布里淵區(qū)中的不同對(duì)稱(chēng)點(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ì)稱(chēng)性，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ì)稱(chēng)方向上聲子色散的實(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.