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

趙 凱，蔡靈倉，張修路，羅 雰

(1.中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理實驗室，四川綿陽 621999；2.西南科技大學極端條件物質特性實驗室，四川綿陽 621010)

1 引 言

過渡金屬Mo是高壓物理中常用的基準材料之一，它的物態方程(Equation of State,EOS)被用于金剛石壓砧(Diamond Anvil Cell,DAC)實驗中紅寶石熒光的標定。自20世紀90年代以來，實驗和理論計算對于Mo的高壓相結構及其穩定性一直存在爭議[1-6]。在室溫下，Mo以體心立方(bcc)結構的形式存在。然而，由于過渡金屬d電子軌道未完全填充的特殊性，高溫高壓下Mo的相轉變路徑仍未得到清楚的闡釋。因此，研究高壓下Mo的相行為，有助于形成對d電子過渡金屬高壓物性的系統了解。

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

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

2 計算方法

2.1 總能和聲子譜計算

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

采用密度泛函微擾理論[14]計算聲子色散曲線。其中，布里淵區k點積分采用14×14×14的MP網格，使用模守恒贗勢[15]描述芯態電子與價電子的相互作用，電子組態為4d55s1，動力學矩陣采用11×11×11的q點網格進行計算。

2.2 準諧德拜近似

采用準諧德拜模型[16]研究Mo在高溫高壓下的熱力學性質，非平衡態Gibbs自由能G*可以寫作

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

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

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

(6)

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

(7)

(8)

Cp=CV(1+αγT)

(9)

(10)

式中：γ是Grüneisen系數，定義為

(11)

3 結果與討論

3.1 結構相的穩定性

為了研究Mo的bcc→dhcp相變的臨界壓力，我們計算了0～1 000 GPa范圍內bcc和dhcp相Mo的零溫焓值Hbcc和Hdhcp，結果如圖1所示。可以看出,Mo在約620 GPa處發生了bcc→dhcp相變，這與前人的計算結果[8]存在約40 GPa的差異，主要原因是采用的波函數基組和交換相關泛函不同。基于準諧德拜模型，研究了不同溫度下Mo的等溫壓縮線。如圖2所示，我們的結果與Hixson等人[17]的經驗計算結果符合得很好，但是293 K溫度下的體積卻明顯小于Wang等人[18]基于經典平均場理論計算的結果，這是因為Wang等人采用全勢線性綴加平面波的方法求解基態能量，而本研究使用的是贗勢平面波方法。可以看出，在620 GPa相變點附近，Mo原子體積有輕微的變化。

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

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

3.2 晶格動力學性質

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

圖3 零壓下bcc-Mo的聲子色散關系(空心點是Farber等人[19]的實驗值)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的聲子色散關系(空心點是Farber等人[19]的實驗值)Fig.4 Phonon dispersion of bcc-Mo at 17.2 GPa (Open symbols are the experimental data by Farber et al.[19])

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

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

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

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

3.3 力學穩定性

基于彈性常數的第一原理計算，預測了bcc和dhcp結構Mo的力學穩定性，結果分別如圖8和表1所示。對于立方晶體，在外部靜水壓作用下的力學穩定性判據如下[22]

(12)

而六角密堆積結構的力學穩定性判據則是

(13)

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

圖8 bcc-Mo的力學穩定性隨壓力的變化關系Fig.8 Mechanical stability of bcc-Mo as a function of pressure

表1 1 000 GPa范圍內dhcp-Mo的力學穩定性Table 1 Mechanical stability of dhcp-Mo up to 1 000 GPa

3.4 熱力學性質

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

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

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

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

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

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

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

4 結 論

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

[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.