過(guò)渡金屬鉬高壓物性的第一原理研究*

趙 凱，蔡靈倉(cāng)，張修路，羅 雰

(1.中國(guó)工程物理研究院流體物理研究所沖擊波物理與爆轟物理實(shí)驗(yàn)室，四川綿陽(yáng) 621999；2.西南科技大學(xué)極端條件物質(zhì)特性實(shí)驗(yàn)室，四川綿陽(yáng) 621010)

1 引 言

過(guò)渡金屬M(fèi)o是高壓物理中常用的基準(zhǔn)材料之一，它的物態(tài)方程(Equation of State,EOS)被用于金剛石壓砧(Diamond Anvil Cell,DAC)實(shí)驗(yàn)中紅寶石熒光的標(biāo)定。自20世紀(jì)90年代以來(lái)，實(shí)驗(yàn)和理論計(jì)算對(duì)于Mo的高壓相結(jié)構(gòu)及其穩(wěn)定性一直存在爭(zhēng)議[1-6]。在室溫下，Mo以體心立方(bcc)結(jié)構(gòu)的形式存在。然而，由于過(guò)渡金屬d電子軌道未完全填充的特殊性，高溫高壓下Mo的相轉(zhuǎn)變路徑仍未得到清楚的闡釋。因此，研究高壓下Mo的相行為，有助于形成對(duì)d電子過(guò)渡金屬高壓物性的系統(tǒng)了解。

Hixson等人[1]實(shí)驗(yàn)測(cè)量了Mo的Hugoniot聲速，發(fā)現(xiàn)其在210 GPa壓強(qiáng)附近出現(xiàn)拐折，故認(rèn)為在沖擊熔化之前Mo先發(fā)生了固-固相變，但隨后的靜高壓實(shí)驗(yàn)結(jié)果[2-3]并不支持這種論斷。與此同時(shí)，Moriarty[4]和Boettger等人[5]以及Belonoshko等人[6]基于第一原理，計(jì)算研究了Mo的體心立方(bcc)、面心立方(fcc)及密排六方(hcp)結(jié)構(gòu)的相穩(wěn)定性，他們的結(jié)果也不支持在210 GPa附近發(fā)生相變。Cazorla等人[7]考慮非諧效應(yīng)之后進(jìn)行計(jì)算，結(jié)果表明，在0～600 GPa范圍內(nèi)，Mo的fcc和hcp高溫相不會(huì)比bcc相更穩(wěn)定。Wang等人[8]根據(jù)粒子群優(yōu)化算法對(duì)鉬的高壓穩(wěn)定結(jié)構(gòu)進(jìn)行了探索，預(yù)測(cè)了零溫下在660 GPa處的bcc→dhcp相變。以上研究多集中于討論Mo高壓相的熱力學(xué)穩(wěn)定性，而決定其動(dòng)力學(xué)穩(wěn)定性的晶格動(dòng)力學(xué)卻鮮有涉及。

本研究利用第一原理密度泛函方法，研究高壓下Mo的bcc和dhcp相的總能和彈性常數(shù)，并利用線性響應(yīng)理論計(jì)算聲子色散關(guān)系；再基于準(zhǔn)諧德拜(Debye)模型研究其熱力學(xué)性質(zhì)，如德拜溫度、熱膨脹系數(shù)和熱容量等。

2 計(jì)算方法

2.1 總能和聲子譜計(jì)算

采用CASTEP程序包[9]計(jì)算Mo在零溫下的總能、彈性常數(shù)和聲子譜。交換相關(guān)泛函選用廣義梯度近似(Generalized Gradient Approximation，GGA)的Wu-Cohen(WC)方法[10]。采用Vanderbilt超軟贗勢(shì)[11]描述離子實(shí)與價(jià)電子之間的相互作用，所使用的電子組態(tài)是4s24p64d55s1。布里淵區(qū)的k點(diǎn)積分采用Monkhorst-Pack(MP)方法[12]。經(jīng)過(guò)仔細(xì)測(cè)試，平面波截?cái)嗄苓x為1 keV，布里淵區(qū)的k點(diǎn)采樣分別選為28×28×28(bcc)和25×25×8(dhcp)。為了消除Pulay應(yīng)力的影響，計(jì)算中使用了有限基組修正。自洽場(chǎng)計(jì)算的收斂閾值設(shè)置為1.0 μeV/atom。幾何優(yōu)化采用BFGS算法[13]，自洽場(chǎng)的收斂條件設(shè)置為：原子能量偏差小于0.2 μeV/atom，最大海爾曼-費(fèi)曼力偏差小于0.5 eV/nm，最大應(yīng)力偏差小于0.1 GPa，最大位移偏差小于0.000 2 nm。

采用密度泛函微擾理論[14]計(jì)算聲子色散曲線。其中，布里淵區(qū)k點(diǎn)積分采用14×14×14的MP網(wǎng)格，使用模守恒贗勢(shì)[15]描述芯態(tài)電子與價(jià)電子的相互作用，電子組態(tài)為4d55s1，動(dòng)力學(xué)矩陣采用11×11×11的q點(diǎn)網(wǎng)格進(jìn)行計(jì)算。

2.2 準(zhǔn)諧德拜近似

采用準(zhǔn)諧德拜模型[16]研究Mo在高溫高壓下的熱力學(xué)性質(zhì)，非平衡態(tài)Gibbs自由能G*可以寫作

G*(V；p,T)=E(V)+pV+Avib[Θ(V)；T]

(1)

式中：V、p、T分別為體積、壓強(qiáng)和溫度，pV表示恒定的靜水壓條件；E(V)是晶體的基態(tài)總能；Θ是德拜溫度；Avib是Helmholtz振動(dòng)自由能，可以寫作

(2)

式中：n是元胞中的原子數(shù);kB是玻爾茲曼常數(shù)；德拜積分D定義為

(3)

對(duì)于各向同性固體，德拜特征溫度Θ可以表示為

(4)

式中：?是約化普朗克常量，M是每個(gè)元胞的分子質(zhì)量，BS是絕熱體積模量，μ是泊松比，f(μ)由下式給出

(5)

在給定的壓力和溫度下，對(duì)非平衡Gibbs函數(shù)G*(V;P,T)求V的極小值，即

[?G*(V；p,T)/?Vp,T=0

(6)

求解(6)式方程，可以獲得熱力學(xué)物態(tài)方程，從而求出等溫體積模量BT、熱容量(定容熱容CV和定壓熱容Cp)，以及熱膨脹系數(shù)α

(7)

(8)

Cp=CV(1+αγT)

(9)

(10)

式中：γ是Grüneisen系數(shù)，定義為

(11)

3 結(jié)果與討論

3.1 結(jié)構(gòu)相的穩(wěn)定性

為了研究Mo的bcc→dhcp相變的臨界壓力，我們計(jì)算了0～1 000 GPa范圍內(nèi)bcc和dhcp相Mo的零溫焓值Hbcc和Hdhcp，結(jié)果如圖1所示。可以看出,Mo在約620 GPa處發(fā)生了bcc→dhcp相變，這與前人的計(jì)算結(jié)果[8]存在約40 GPa的差異，主要原因是采用的波函數(shù)基組和交換相關(guān)泛函不同?；跍?zhǔn)諧德拜模型，研究了不同溫度下Mo的等溫壓縮線。如圖2所示，我們的結(jié)果與Hixson等人[17]的經(jīng)驗(yàn)計(jì)算結(jié)果符合得很好，但是293 K溫度下的體積卻明顯小于Wang等人[18]基于經(jīng)典平均場(chǎng)理論計(jì)算的結(jié)果，這是因?yàn)閃ang等人采用全勢(shì)線性綴加平面波的方法求解基態(tài)能量，而本研究使用的是贗勢(shì)平面波方法?？梢钥闯觯?20 GPa相變點(diǎn)附近，Mo原子體積有輕微的變化。

圖1 bcc和dhcp鉬的零溫焓值之差ΔH隨外界靜水壓的變化關(guān)系Fig.1 Enthalpy difference between bcc and dhcp Mo vs. pressure

圖2 不同溫度下計(jì)算得到的bcc和dhcp鉬的等溫壓縮線Fig.2 Calculated isotherms of bcc and dhcp Mo at given temperatures

3.2 晶格動(dòng)力學(xué)性質(zhì)

圖3、圖4分別給出了零壓和17.2 GPa下bcc-Mo的聲子譜，并與Farber等人[19]的非彈性X射線散射實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比?？梢钥闯?，當(dāng)外界壓力為17.2 GPa時(shí)，沿著Γ-N方向，TA[110]?001?橫聲學(xué)模的實(shí)驗(yàn)值明顯高于計(jì)算值；而沿著Γ-H方向，TA模在q=0.8處出現(xiàn)明顯的過(guò)度彎曲，這與Farber等人[19]計(jì)算結(jié)果的趨勢(shì)一致。零壓時(shí)，我們的計(jì)算結(jié)果與實(shí)驗(yàn)值在大體趨勢(shì)上符合很好，但TA[110]?001?橫聲學(xué)模在布里淵區(qū)邊界(即N點(diǎn)附近)與實(shí)驗(yàn)值差距仍較為明顯。

圖3 零壓下bcc-Mo的聲子色散關(guān)系(空心點(diǎn)是Farber等人[19]的實(shí)驗(yàn)值)Fig.3 Phonon dispersion of bcc-Mo at 0 GPa (Open symbols are the experimental data by Farber et al.[19])

圖4 17.2 GPa下bcc-Mo的聲子色散關(guān)系(空心點(diǎn)是Farber等人[19]的實(shí)驗(yàn)值)Fig.4 Phonon dispersion of bcc-Mo at 17.2 GPa (Open symbols are the experimental data by Farber et al.[19])

為研究Mo的高壓相變，我們計(jì)算了600～1 000 GPa范圍內(nèi)Mo的聲子譜，部分結(jié)果如圖5所示。當(dāng)外界壓力從600 GPa上升到800 GPa時(shí)，TA[110]?001?橫聲學(xué)模在q=0.25附近出現(xiàn)強(qiáng)烈的軟化。當(dāng)壓力達(dá)到1 000 GPa時(shí)，TA[110]?001?模在q=0.25和0.33處均出現(xiàn)虛頻，而這兩個(gè)振動(dòng)模分別對(duì)應(yīng)bcc→dhcp和bcc→9R的結(jié)構(gòu)相變[20]。盡管bcc-Mo的動(dòng)力學(xué)不穩(wěn)定性出現(xiàn)在1 000 GPa附近，但是結(jié)構(gòu)相變可能在此之前已經(jīng)發(fā)生。這與Krasilnikov等人[21]計(jì)算的1 000 GPa聲子譜有較大差異，他們的結(jié)果中，僅在q=0.25 附近出現(xiàn)虛頻。因此，高壓下dhcp和9R結(jié)構(gòu)的相對(duì)穩(wěn)定性還需要進(jìn)一步深入的研究。為進(jìn)一步精確地判定bcc-Mo發(fā)生動(dòng)力學(xué)不穩(wěn)定的壓力點(diǎn)，我們分別計(jì)算了925和950 GPa下的TA[110]?001?模，結(jié)果如圖6所示，在950 GPa處聲子譜就已經(jīng)出現(xiàn)了虛頻。圖7給出了600和900 GPa下dhcp-Mo的聲子譜，可以看出，在600 GPa時(shí)dhcp-Mo已經(jīng)是動(dòng)力學(xué)穩(wěn)定的結(jié)構(gòu)，從而保證了本研究所確定的620 GPa相變點(diǎn)的有效性。

圖5 不同靜水壓下bcc-Mo的聲子色散關(guān)系Fig.5 Phonon dispersions of bcc-Mo under different pressures

圖6 bcc-Mo在925和950 GPa下的聲子色散關(guān)系Fig.6 Phonon dispersions of bcc-Mo at 925 and 950 GPa

圖7 不同靜水壓下dhcp-Mo的聲子色散關(guān)系Fig.7 Phonon dispersions of dhcp-Mo under different pressures

3.3 力學(xué)穩(wěn)定性

基于彈性常數(shù)的第一原理計(jì)算，預(yù)測(cè)了bcc和dhcp結(jié)構(gòu)Mo的力學(xué)穩(wěn)定性，結(jié)果分別如圖8和表1所示。對(duì)于立方晶體，在外部靜水壓作用下的力學(xué)穩(wěn)定性判據(jù)如下[22]

(12)

而六角密堆積結(jié)構(gòu)的力學(xué)穩(wěn)定性判據(jù)則是

(13)

從圖8可以看出，在0～1 000 GPa范圍內(nèi)，bcc-Mo是力學(xué)穩(wěn)定的，我們所計(jì)算的彈性穩(wěn)定性與Krasilnikov等人[21]計(jì)算結(jié)果的趨勢(shì)一致，但呈現(xiàn)出較大的波動(dòng)，表明bcc-Mo的力學(xué)性能隨著壓力的增加而表現(xiàn)出復(fù)雜的變化行為。表1列出了不同壓力下dhcp-Mo的力學(xué)穩(wěn)定性判據(jù)。結(jié)果表明，在400 GPa以上，dhcp結(jié)構(gòu)都是穩(wěn)定的。結(jié)合以上熱力學(xué)和晶格動(dòng)力學(xué)計(jì)算結(jié)果可以判斷，dhcp結(jié)構(gòu)在620 GPa以上是可以穩(wěn)定存在的相。

圖8 bcc-Mo的力學(xué)穩(wěn)定性隨壓力的變化關(guān)系Fig.8 Mechanical stability of bcc-Mo as a function of pressure

表1 1 000 GPa范圍內(nèi)dhcp-Mo的力學(xué)穩(wěn)定性Table 1 Mechanical stability of dhcp-Mo up to 1 000 GPa

3.4 熱力學(xué)性質(zhì)

利用準(zhǔn)諧德拜模型，計(jì)算了bcc和dhcp相Mo在不同溫度下的熵S、德拜溫度Θ、熱膨脹系數(shù)α，以及Grüneisen參數(shù)γ。如圖9所示，在給定溫度下，Mo的德拜溫度隨壓力的升高而增大,而在給定壓力下，它隨著溫度的升高而減小，且在620 GPa相變點(diǎn)附近變化很小。圖10給出了熵S隨壓力和溫度的變化關(guān)系。可以看出，S隨著壓力的增加逐漸減小，隨著溫度的升高而增大，說(shuō)明外界壓力和溫度的施加分別對(duì)應(yīng)著有序和無(wú)序的競(jìng)爭(zhēng)。應(yīng)該注意到，德拜模型所能計(jì)算的熵，乃是出于對(duì)理想聲子氣體的統(tǒng)計(jì)。圖7說(shuō)明，在相變點(diǎn)附近，聲子氣體的混亂程度幾乎沒有變化。

圖9 不同溫度下Mo的德拜溫度隨壓力的變化關(guān)系Fig.9 Debye temperature of Mo vs. pressure at different temperatures

圖10 不同溫度下Mo的熵隨壓力的變化關(guān)系Fig.10 Entropy of Mo as a function of pressure at different temperatures

圖11是不同溫度下Mo的熱膨脹系數(shù)α隨壓力的變化關(guān)系。由圖11可以看出，在100 GPa以下，隨著壓力的增加,α的變化較為劇烈，且不同溫度下α的差值變化也較大；而當(dāng)壓力高于300 GPa后，隨著壓力的增加，α變化較小，且不同溫度下α的差值幾乎不變。這說(shuō)明當(dāng)壓力較小時(shí)，尚不能對(duì)溫度產(chǎn)生的膨脹效應(yīng)起到很好的抑制作用；而當(dāng)壓力足夠高時(shí)，即便非常高的溫度也難以讓原子間距明顯地增大。在620 GPa附近，bcc和dhcp結(jié)構(gòu)相的熱膨脹系數(shù)差異很小，但均隨著溫度的升高而逐漸增大。

如圖12所示，在給定溫度下,Mo的定容熱容CV隨壓力升高而減小，且溫度越低,熱容曲線的斜率越大，下降得就越快;在5 000 K以上,CV幾乎不隨壓力變化，并趨近經(jīng)典的Dulong-Petit極限3nNAkB(n是分子中的原子數(shù)，NA為阿伏伽德羅常數(shù)，kB為波爾茲曼常數(shù))，對(duì)于Mo，n=1，Dulong-Petit極限約等于24.95 J/(mol·K)。而在一定壓力下，CV隨著溫度的升高而增大。這說(shuō)明增加壓力和降低溫度對(duì)調(diào)節(jié)熱容量的變化有相同的效果。

圖11 不同溫度下Mo的熱膨脹系數(shù)隨壓力的變化關(guān)系Fig.11 Thermal expansions of Mo vs. pressure at different temperatures

圖12 不同溫度下Mo的熱容量隨壓力的變化關(guān)系Fig.12 Isochoric heat capacity of Mo as a function of pressure at different temperatures

4 結(jié) 論

基于平面波贗勢(shì)密度泛函理論，研究了過(guò)渡金屬M(fèi)o在高壓下的結(jié)構(gòu)性質(zhì)、相變、晶格動(dòng)力學(xué)以及熱力學(xué)性質(zhì)。計(jì)算得出的等溫壓縮線與前人的理論計(jì)算結(jié)果符合得很好。通過(guò)比較bcc和dhcp相的焓值之差，預(yù)測(cè)了620 GPa處發(fā)生bcc→dhcp的結(jié)構(gòu)相變。聲子譜的計(jì)算結(jié)果在較低壓力下與實(shí)驗(yàn)非常吻合，而在較高壓力下，bcc結(jié)構(gòu)可能會(huì)向dhcp或9R結(jié)構(gòu)轉(zhuǎn)變，其相變路徑還有待進(jìn)一步深入的研究。最后，系統(tǒng)地研究了高壓高溫下Mo的熱力學(xué)性質(zhì)，如德拜溫度、熵、熱膨脹系數(shù)以及熱容量等。計(jì)算結(jié)果表明，在相變臨界點(diǎn)620 GPa附近，Mo的熱力學(xué)性質(zhì)并沒有發(fā)生顯著的變化；而升高壓力或降低溫度，對(duì)于熱力學(xué)性質(zhì)的調(diào)制有相同的效果。

[1] HIXSON R S,BONESS D A,SHANER J W,et al.Acoustic velocities and phase transitions in molybdenum under strong shock compression [J].Phys Rev Lett,1989,62(6):637.

[2] VOHRA Y K,RUOFF A L.Static compression of metals Mo,Pb,and Pt to 272 GPa:Comparison with shock data [J].Phys Rev B,1990,42(13):8651.

[3] SANTAMARA-PéREZ D,ROSS M,ERRANDONEA D,et al.X-ray diffraction measurements of Mo melting to 119 GPa and the high pressure phase diagram [J].J Chem Phys,2009,130(12):124509.

[4] MORIARTY J A.Ultrahigh-pressure structural phase transitions in Cr,Mo,and W [J].Phys Rev B,1992,45(5):2004.

[5] BOETTER J C.Relativistic effects on the structural phase stability of molybdenum [J].J Phys:Condens Matter,1999,11(16):3237.

[6] BELONOSHKO A B,BURAKOVSKY L,CHEN S P,et al.Molybdenum at high pressure and temperature:Melting from another solid phase [J].Phys Rev Lett,2008,100(13):135701.

[7] CAZORLA C,ALFD,GILLAN M J.Constraints on the phase diagram of molybdenum from first-principles free-energy calculations [J].Phys Rev B,2012,85(6):064113.

[8] WANG B,ZHANG G B,WANG Y X.Predicted crystal structures of molybdenum under high pressure [J].J Alloys Compd,2013,556:116-120.

[9] PAYNE M C,TETER M P,ALLAN D C,et al.Iterative minimization techniques forabinitiototal-energy calculations:Molecular dynamics and conjugate gradients [J].Rev Mod Phys,1992,64(4):1045.

[10] WU Z,COHEN R E.More accurate generalized gradient approximation for solids [J].Phys Rev B,2006,73(23):235116.

[11] VANDERBILT D.Soft self-consistent pseudopotentials in a generalized eigenvalue formalism [J].Phys Rev B,1990,41(11):7892.

[12] MONKHORST H J,PACK J D.Special points for Brillouin-zone integrations [J].Phys Rev B,1976,13(12):5188.

[13] PFROMMER B G,CTé M,LOUIE S G,et al.Relaxation of crystals with the quasi-Newton method [J].J Comput Phys,1997,131:233-240.

[14] REFSON K,TULIP P R,CLARK S J.Variational density-functional perturbation theor for dielectrics and lattice dynamics [J].Phys Rev B,2006,73(15):155114.

[15] LIN J S,QTEISH A,PAYNE M C,et al.Optimized and transferable nonlocal separableabinitiopseudopotentials [J].Phys Rev B,1993,47(8):4174.

[16] OTERO-DE-LA-ROZA A,PéREZ D A,LUAA V.GIBBS2:A new version of the quasiharmonic model code.Ⅱ.Models for solid-state thermodynamics,features and implementation [J].Comput Phys Commun,2011,182(10):2232-2248.

[17] HIXSON R S,FRITZ J N.Shock compression of tungsten and molybdenum [J].J Appl Phys,1992,71(4):1721.

[18] WANG Y,CHEN D,ZHANG X.Calculated equation of state of Al,Cu,Ta,Mo,and W to 1 000 GPa [J].Phys Rev Lett,2000,84(15):3220.

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

[20] GRIMVALL G,K?PE B M,OZOLINV,et al.Lattice instabilities in metallic elements [J].Rev Mod Phys,2012,84:945.

[21] KRASILNIKOV O M,BELOV M P,LUGOVSKOY A V,et al.Elastic properties,lattice dynamics and structural transitions in molybdenum at high pressures [J].Comput Mater Sci,2014,81:313.

[22] SIN’KO G V,SMIRNOV N A.Abinitiocalculations of elastic constants and thermodynamic properties of bcc,fcc and hcp Al crystals under pressure [J].J Phys:Condens Matter,2002,14:6989.