方鎂石高壓結(jié)構(gòu)預(yù)測和高溫結(jié)構(gòu)穩(wěn)定性研究*

宋婷 孫小偉 魏小平 歐陽玉花 張春林 郭鵬 趙煒

(蘭州交通大學(xué)數(shù)理學(xué)院,蘭州 730070)

1 引 言

方鎂石是地球下地幔礦物—鎂方鐵礦(Mg,Fe)O的主要組分,化學(xué)組成為氧化鎂(MgO).在下地幔中,方鎂石的含量僅次于鈣鈦礦(MgSiO3)[1],同時方鎂石也是耐火陶瓷工業(yè)的基本材料,對其高溫高壓性質(zhì)的研究具有重要的地球物理意義和重要的工業(yè)應(yīng)用價值.由于MgO在凝聚態(tài)物理和地球科學(xué)研究中的重要性,其高溫高壓熱力學(xué)性質(zhì)已經(jīng)通過靜高壓、動高壓實驗及計算機模擬方法進行了廣泛的研究[2-5].盡管MgO材料的基礎(chǔ)物理性質(zhì)目前已經(jīng)得到了較好的了解,但理論預(yù)測的熔化線和高壓下實驗測量結(jié)果之間存在巨大的分歧[3,6],為澄清歧見人們展開了對MgO高壓結(jié)構(gòu)的進一步研究[7,8];MgO高壓相變研究以及基于各種可能新相結(jié)構(gòu)的熱力學(xué)特性的可靠性預(yù)測正逐漸成為材料高壓科學(xué)和凝聚態(tài)物理領(lǐng)域的一個新課題.

有關(guān)MgO晶體的高壓物性研究,早期主要集中在相變部分.MgO具有簡單的立方巖鹽結(jié)構(gòu)(即NaCl結(jié)構(gòu),用B1表示),高壓下會發(fā)生從巖鹽結(jié)構(gòu)到立方氯化銫結(jié)構(gòu)(即CsCl結(jié)構(gòu),用B2表示)的相轉(zhuǎn)變.1995年,Duffy等[9]的實驗研究發(fā)現(xiàn)壓力高達227 GPa時,巖鹽結(jié)構(gòu)的MgO依然穩(wěn)定存在,和其后諸多理論計算的預(yù)測結(jié)果相一致[10-12].促使人們對MgO高壓結(jié)構(gòu)進行進一步深入研究的一個重要原因,是其理論預(yù)測和實驗測量的熔化曲線之間存在很大的分歧,而不同理論方法的差別和實驗系統(tǒng)誤差不足以說明得到的熔化數(shù)據(jù)存在巨大分歧的原因.1994年Zerr和Boehler[13]給出了MgO壓力僅上升到31.5 GPa的靜高壓實驗熔化溫度數(shù)據(jù),測量得到的熔化線斜率為36 K/GPa,而2008年Zhang和Fei[3]的沖擊波實驗測量值為221 K/GPa,兩者之間相差6倍多.分析誤差產(chǎn)生的原因,除樣品遭受了非靜水壓或熱壓力的影響外,還有可能是靜高壓裝置中判斷樣品發(fā)生熔化的依據(jù)有問題,即對樣品表面觀測到的對流也許僅僅意味著樣品表面而非樣品內(nèi)部發(fā)生了熔化;在沖擊壓實驗中,Anderso和Duba[14]也提出過同樣的問題,他們指出由Hugoniot聲速測量結(jié)果計算的熔化溫度中存在“過熱”熔化現(xiàn)象,即晶體在高于自身熔點的溫度下發(fā)生熔化的現(xiàn)象,但并沒有對這個問題做進一步的定量分析;本課題組曾利用單相分子動力學(xué)模擬方法將0.1 MPa,31.5 GPa,47 GPa和135 GPa壓力下的分子動力學(xué)熱不穩(wěn)定性數(shù)據(jù)和實驗數(shù)據(jù)比較得到的過熱率經(jīng)過檢驗,用到了高壓情況,外推得到了MgO達150 GPa的高壓熔化曲線[15],結(jié)果遠高于Zerr和Boehler[13]靜高壓實驗給出的數(shù)據(jù).

模擬地球深部的高溫高壓條件的實驗手段主要有靜高壓和動高壓實驗技術(shù).靜高壓技術(shù)由于受到金剛石壓砧所能實現(xiàn)的壓強上限的限制,測量能達到的壓力范圍有限,目前金剛石單晶壓砧技術(shù)可以達到的最高壓強約為350 GPa[16],但實驗上100 GPa以上高壓的實現(xiàn)和控制非常困難;采用沖擊壓縮技術(shù)雖然可以在比較高的壓力和溫度范圍內(nèi)研究材料的高壓物性,但是由于沖擊波傳播速度很快,在樣品中持續(xù)的時間很短,會給實際測量帶來很多困難.Zhang和Fei[3]研究報道了單晶MgO在114和192 GPa沖擊壓縮下的雨貢紐數(shù)據(jù),結(jié)合前人的沖擊波數(shù)據(jù)分析了沿MgO的P-V雨貢紐線在(170 ± 10)GPa存在體積不連續(xù)的原因,認為是MgO發(fā)生從NaCl立方體結(jié)構(gòu)到NiAs六角密堆積結(jié)構(gòu)(用B81表示)的相變所引起的;Aguado和Madden[6]的模擬計算顯示MgO從纖鋅礦相(用B4表示)發(fā)生熔化的溫度比從巖鹽相發(fā)生熔化的溫度低很多,可以解釋高壓下理論計算的熔化線和實驗結(jié)果相差很大的矛盾.最近,盡管針對MgO的B81和B4結(jié)構(gòu),國際上開展了一些新的預(yù)測性研究工作[7,8],但MgO存在B81和B4新結(jié)構(gòu)的合理性還需要進一步的論證.

眾所周知,材料不同的結(jié)構(gòu)必然導(dǎo)致不同的性質(zhì),MgO晶體各種熱力學(xué)特性的可靠預(yù)測研究,其高壓結(jié)構(gòu)的細致研究是前提.本文利用基于密度泛函理論(DFT)的第一性原理計算方法,對MgO在0—800 GPa壓力范圍內(nèi)的各種可能結(jié)構(gòu)及利用基于粒子群優(yōu)化算法的卡利普索(CALYPSO)晶體結(jié)構(gòu)預(yù)測方法得到的潛在結(jié)構(gòu)進行了系統(tǒng)研究,發(fā)現(xiàn)MgO在0—580 GPa的壓力范圍內(nèi)一直以穩(wěn)定巖鹽結(jié)構(gòu)存在,隨著壓力的增加,會發(fā)生從B1到B2結(jié)構(gòu)的相轉(zhuǎn)變;盡管B81和B4結(jié)構(gòu)能合理解釋沖擊壓縮實驗中沿MgO的P-V雨貢紐線在(170 ± 10)GPa存在體積不連續(xù)和高壓下理論計算的熔化線與實驗結(jié)果相差很大的原因,但在壓力高達800 GPa下都不能夠穩(wěn)定存在.在MgO高壓結(jié)構(gòu)預(yù)測基礎(chǔ)上,本文分別利用殼層和呼吸殼層模型分子動力學(xué)方法,結(jié)合檢驗的有效經(jīng)驗勢參數(shù),在500—10000 K的溫度范圍和0—150 GPa的壓力范圍內(nèi)對其穩(wěn)定巖鹽結(jié)構(gòu)的熔化特性進行了模擬.

2 計算方法

2.1 MgO高壓結(jié)構(gòu)預(yù)測及優(yōu)化

本文使用由吉林大學(xué)馬琰銘小組開發(fā)的擁有自主知識產(chǎn)權(quán)的基于粒子群優(yōu)化算法的CALYPSO軟件包[17],預(yù)測MgO的1—4倍分子式組成的晶胞分別在0,135,300和500 GPa壓力下的晶體結(jié)構(gòu).該軟件只需給定化學(xué)配比,就能在特定的壓力和溫度條件下預(yù)測材料的基態(tài)和亞穩(wěn)態(tài)結(jié)構(gòu).這種預(yù)測結(jié)構(gòu)的方法已經(jīng)成功應(yīng)用在大量實驗合成的單質(zhì)和化合物的常壓及高壓晶體結(jié)構(gòu)研究中[18].例如,2014年,Li等[19]利用該方法對H2S在10—200 GPa壓力區(qū)間的結(jié)構(gòu)進行預(yù)測,提出了兩個能量穩(wěn)定且具有金屬性的新高壓相P-1和Cmca,并首次預(yù)言高壓下H2S的高溫超導(dǎo)電性;在這一預(yù)測工作的啟發(fā)下,Drozdov等[20]開展了S-H體系的高壓超導(dǎo)實驗研究,發(fā)現(xiàn)S-H體系在高壓下呈現(xiàn)兩個超導(dǎo)態(tài),其中低溫高壓下測得的Tc與Li等[19]的計算數(shù)據(jù)基本吻合,而室溫退火后測得的Tc達到驚人的203 K;2017年,Peng等[21]利用該方法預(yù)言了更多富氫包合物結(jié)構(gòu)高壓高溫超導(dǎo)體的存在.本文對預(yù)測得到的結(jié)構(gòu)采用第一性原理計算的VASP軟件包進行幾何優(yōu)化[22],電子和電子之間的交換關(guān)聯(lián)勢采用廣義梯度近似(GGA)下的Perdew-Burke-Ernzerhof泛函進行處理[23],計算過程中Mg的2p6,3s2和O的2s2,2p4電子均被看作價電子處理,價電子和芯電子之間的相互作用采用投影綴加平面波描述[24],倒易空間布里淵區(qū)k點采用Monkhorst-Pack方法選取[25],精確幾何優(yōu)化中的最大分割間隔為2π×0.08 ?—1.

獲得MgO預(yù)測結(jié)構(gòu)后,應(yīng)用CASTEP軟件包對預(yù)測結(jié)構(gòu)、實驗結(jié)構(gòu)及其他可能存在的結(jié)構(gòu)進行結(jié)構(gòu)優(yōu)化并計算相關(guān)性質(zhì)[26],這里電子與電子間的交換關(guān)聯(lián)能采用2008年由Perdew等[27]修訂的Perdew-Burke-Ernzerhof形式的GGA方法,該形式的交換相關(guān)泛函能夠提高密堆固體的平衡特性;芯電子與價電子間的相互作用勢由超軟贗勢[28]實現(xiàn).為了確定計算收斂性,文中對平面波截斷能和k點采樣進行了收斂性測試.針對MgO的B1,B2,B4等7種不同結(jié)構(gòu),截斷能均取600 eV,就能保證高壓下的收斂精度.零壓和不同外壓條件下的MgO晶體結(jié)構(gòu)和原子位置優(yōu)化采用Broyden,Fletcher,Goldfarb和Shanno提出并發(fā)展起來的幾何算法[29].聲子譜的計算采用線性響應(yīng)方法[30].

2.2 分子動力學(xué)模擬技術(shù)

論文針對MgO巖鹽結(jié)構(gòu),研究了高壓下溫度對該結(jié)構(gòu)穩(wěn)定性的影響.巖鹽結(jié)構(gòu)的MgO模擬體系由5×5×5個單位原胞組成,通過對多個單位原胞組成的超胞施加三維周期性邊界條件,從而使離子的運動空間成為贗無限.粒子的初始速度根據(jù)模擬體系設(shè)定的溫度賦值,壓力由維里定理得到,模擬中長程靜電力的處理采用Ewald求和技術(shù)[31],分別在實空間和倒空間中進行計算.模擬選用等溫等壓系綜,即NPT系綜.模擬時間步長設(shè)為1 fs,每個平衡態(tài)計算20000步(20 ps),前10000步(10 ps)是趨向平衡階段,以使體系平衡到所設(shè)定的溫度和壓力,然后再計算10000步(10 ps),以測量指定壓力下MgO模擬體系的溫度和體積.

對MgO晶體進行經(jīng)典分子動力學(xué)模擬研究,離子間相互作用勢函數(shù)和可靠作用勢參數(shù)的選取最為關(guān)鍵.本文中,MgO晶體各離子間的總相互作用能Vij由長程庫侖能和短程非庫侖作用能組成,短程相互作用取Buckingham形式：

(1)式右邊各項分別表示庫侖能、排斥能項和范德瓦耳斯項(偶極-偶極項),其中Z為有效電荷;r為原子之間的相互作用距離;A和ρ是排斥勢參數(shù);C為范德瓦耳斯常數(shù).為了很好地描述離子的極化特性,模擬中考慮了廣為使用的殼層模型[32],即認為每一個點離子由核和殼兩部分組成,若核所帶電荷為X,殼為Y,則這個點離子的總電荷為X+Y.相同離子的核和殼之間不存在靜電相互作用,認為它們之間由一彈簧相連,這個彈簧的彈性常數(shù)為K;該離子在周圍其他離子的電場中得以極化,設(shè)極化率為η,則η與該離子的殼電荷Y及彈性常數(shù)K有如下的關(guān)系式：

引入殼層模型能很好地描述離子的極化特性,模擬選用了兩套非常相似的殼層模型勢參數(shù)：Lewis和Catlow[33]給出的勢參數(shù)(SM-LC)及Stoneham和Sangster[34]給出的勢參數(shù)(SM-SS),來模擬MgO體系的熔化特性.

在殼層模型基礎(chǔ)上考慮壓縮效應(yīng)就是呼吸殼模型,該模型允許氧離子的排斥半徑在晶體中其他離子的作用下各向同性地變形,每個離子的核和呼吸殼層通過諧波勢連接,因此,系統(tǒng)的總勢能包括呼吸能量項Ubre：

式中,ki是呼吸離子i的彈性常數(shù);A0,i是自由離子i的排斥半徑;Ai是呼吸殼層模型中允許離子i的排斥半徑在晶體中其他離子作用下各向同性變形的排斥半徑[35].這些短程對勢參數(shù)、殼層電荷Y、描述殼層模型極化特性的彈性常數(shù)K和呼吸殼層模型中描述呼吸離子的彈性常數(shù)k通過擬合MgO晶體的結(jié)構(gòu)參數(shù)得到.最近,Ball和Grimes擬合了一套新的MgO經(jīng)驗勢參數(shù)(記為BG)[36],其中Mg和O離子僅帶部分電荷(± 1.7e);正如Henkelman等[36]指出的：除非離子之間的距離非常短,否則Ball和Grimes擬合的MgO經(jīng)驗勢參數(shù)[36]與Lewis和Catlow擬合的勢參數(shù)[33]相比,差別很小.這些擬合勢參數(shù)一并由表1和表2列出[33-36],以做全面比較.

表1 巖鹽結(jié)構(gòu)的MgO晶體特性模擬的短程對勢參數(shù)Table 1.Parameters of the short-range pair potentials for simulation of MgO crystal properties with NaCl-type structure.

表2 巖鹽結(jié)構(gòu)的MgO晶體特性模擬的殼層及呼吸殼層參數(shù)Table 2.Shell and breathing shell parameters of MgO crystal characteristic simulation with NaCltype structure.

2.3 計算方法可靠性檢驗

圖1 模擬得到的零溫下MgO巖鹽結(jié)構(gòu)的體積比率隨壓力的變化及和實驗結(jié)果的比較Fig.1.Volume ratios of MgO with NaCl-type structure under pressures at zero temperature,in comparison with the experimental data.

為了檢驗計算方法及選取參數(shù)的可靠性,本文利用殼層和呼吸殼層模型分子動力學(xué)方法及基于DFT的GGA及局域密度近似(LDA)方法,分別計算了零溫下MgO巖鹽結(jié)構(gòu)從零到地球核幔邊界壓力(135 GPa)范圍內(nèi)的體積比率隨壓力的變化,所得結(jié)果與靜壓和非靜壓X射線衍射實驗結(jié)果[2]及金剛石壓砧(DAC)靜高壓裝置得到的結(jié)果[37]進行了比較,如圖1所示.在準(zhǔn)靜水壓條件下,Fei[2]用氖做傳壓介質(zhì)和氯化鈉做壓標(biāo)測得了溫度為300 K、壓強上升到23 GPa時MgO的DAC靜態(tài)壓縮實驗值,呼吸殼層模型分子動力學(xué)模擬的壓力隨體積比率變化的關(guān)系和Fei[2]上升到23 GPa時的靜態(tài)數(shù)據(jù)吻合得很好,結(jié)合SM-LC勢參數(shù)的殼層模型分子動力學(xué)模擬結(jié)果和Duffy等[9]整理的上升到100 GPa時的沖擊壓縮數(shù)據(jù)一致.同時,呼吸殼層模型分子動力學(xué)模擬結(jié)果和第一性原理GGA計算結(jié)果在高壓下符合得很好,由于GGA和LDA方法上的差別,使得基于LDA的第一性原理計算往往低估晶格常數(shù)從而低估模擬體積,但這種低估在低壓下不明顯.和結(jié)合SM-SS勢參數(shù)的殼層模型分子動力學(xué)模擬結(jié)果相比,呼吸殼層模型在描述MgO高壓物態(tài)方程性質(zhì)時,壓縮效應(yīng)體現(xiàn)的尤為明顯.

3 分析和討論

3.1 高壓結(jié)構(gòu)預(yù)測

常壓下MgO以B1結(jié)構(gòu)存在,高壓下會發(fā)生從B1到B2結(jié)構(gòu)的相變[9-12].Zhang和Fei[3]研究報道了單晶MgO在114和192 GPa沖擊壓縮下的Hugnonit數(shù)據(jù),結(jié)合前人的沖擊波數(shù)據(jù)分析了沿MgO的P-VHugnonit線在(170 ± 10)GPa存在體積不連續(xù)的原因,認為是MgO發(fā)生從B1立方體結(jié)構(gòu)到B81六角密堆積結(jié)構(gòu)的相變所引起的[3];Aguado和Madden[6]的模擬計算顯示MgO從B4纖鋅礦結(jié)構(gòu)發(fā)生熔化的溫度比從B1結(jié)構(gòu)發(fā)生熔化的溫度低很多,可以解釋高壓下理論計算的熔化線和實驗結(jié)果相差很大的矛盾.

圖2 MgO晶體的可能結(jié)構(gòu)(其中大球代表Mg原子,小球代表O原子)(a)B1;(b)B2;(c)B3;(d)B4;(e)B81;(f)B82;(g)P3m1Fig.2.Probable crystal structures of MgO (the large and small spheres represent magnesium and oxygen atoms,respectively)：(a)B1;(b)B2;(c)B3;(d)B4;(e)B81;(f)B82;(g) P3m1.

另外,還考慮了MgO的閃鋅礦結(jié)構(gòu)(用B3表示),這可以通過與其所在的堿土金屬氧化物進行對比來說明.MgO是排在BeO之后CaO之前的堿土金屬氧化物.本課題組曾對CaO的高壓物性進行過研究[38,39],選擇CaO的原因是它的B1→B2相變壓力已由實驗確定為60 GPa左右[40],可作為檢驗?zāi)M計算準(zhǔn)確與否的標(biāo)準(zhǔn),而且CaO常壓下B1結(jié)構(gòu)到高壓下B2結(jié)構(gòu)的相變并無疑議.BeO不同于其他的堿土金屬氧化物,結(jié)晶時形成的是穩(wěn)定纖鋅礦結(jié)構(gòu),而其他的堿土金屬氧化物結(jié)晶時形成的是立方巖鹽結(jié)構(gòu);根據(jù)Phillips電離度理論(電離度 ≥ 0.785的二元化合物在結(jié)晶時形成的是B1結(jié)構(gòu),而電離度 ＜ 0.785的二元化合物形成的是B4結(jié)構(gòu)或B3結(jié)構(gòu))[41],B4結(jié)構(gòu)的BeO在高壓下會轉(zhuǎn)變?yōu)锽1結(jié)構(gòu),而近年的理論研究表明,B4結(jié)構(gòu)會先轉(zhuǎn)變?yōu)锽3結(jié)構(gòu),然后再轉(zhuǎn)變?yōu)锽1結(jié)構(gòu),van Camp和van Doren[42]預(yù)測了BeO晶體B4 → B3 → B1的相變壓力分別為74和137 GPa,計算中采用了軟的非局域贗勢.直到最近,Cai等[43]采用第一性原理贗勢和GGA研究了這兩個相變的勢壘,認為B4 → B3這一相變不可能發(fā)生,僅B4 → B1這一相變可以在高壓下發(fā)生;Luo等[44]的最新計算表明隨著壓力增加到100 GPa,BeO直接由B4結(jié)構(gòu)轉(zhuǎn)變?yōu)锽1結(jié)構(gòu).事實上,由于計算誤差以及B3和B4結(jié)構(gòu)的能量差十分接近(在零溫零壓時的能量差都在平均每個原子10 meV左右),簡單的能量計算結(jié)果不能作為判斷B3結(jié)構(gòu)是否存在的有力證據(jù).同時兩種結(jié)構(gòu)非常相似(B4結(jié)構(gòu)是六角密堆形式,B3結(jié)構(gòu)是立方密堆形式),導(dǎo)致這兩種結(jié)構(gòu)的物性也非常相似,其中之一就是非常接近的內(nèi)能,這就解釋了為什么在不同的計算中有時存在B3結(jié)構(gòu),有時不存在B3結(jié)構(gòu).最近Aguado和Madden[6]的模擬計算顯示MgO從B4相發(fā)生熔化的溫度比從B1相發(fā)生熔化的溫度低很多,用來解釋高壓下熔化線的理論計算和實驗測量之間存在巨大差別的原因.如果MgO和BeO一樣,存在穩(wěn)定的B4相,也應(yīng)該存在一個和B4相非常相似的B3相.以上分析的MgO各種可能存在的結(jié)構(gòu)B1,B2,B3,B4,B81及基于粒子群優(yōu)化算法預(yù)測的晶體結(jié)構(gòu)B82和P3m1的結(jié)構(gòu)示意圖由圖2給出.

贗勢平面波方法的一大優(yōu)點就是它能夠自動根據(jù)原子的受力情況來調(diào)整原子的位置和晶胞參數(shù),直到所有原子的受力都為零,從而使整個體系的總能達到最低,找到給定條件下的最穩(wěn)定結(jié)構(gòu).圖3給出了利用平面波贗勢結(jié)合GGA,DFT方法得到的MgO的7種可能結(jié)構(gòu)的每個分子式的焓差隨壓力的變化,其中以B1結(jié)構(gòu)作為參考.由于實驗一般都是在一定溫度T、壓力P下進行的(由于熱脹冷縮,使體積V固定的實驗很難進行),所以嚴(yán)格來說此時應(yīng)以吉布斯函數(shù)G判定相的穩(wěn)定性.當(dāng)溫度為零時,G在數(shù)值上等于焓H,這時可通過熱力學(xué)函數(shù)焓來判斷熱力學(xué)相的穩(wěn)定性.本工作中所加外壓均為流體靜壓力.在給定外壓下,各相的相對穩(wěn)定性可由焓H=E+PV來確定,焓最小的結(jié)構(gòu)最穩(wěn)定.可以看出零溫下MgO最穩(wěn)定的結(jié)構(gòu)為B1結(jié)構(gòu),隨著壓力增加到580 GPa,MgO會發(fā)生從B1到B2結(jié)構(gòu)的相轉(zhuǎn)變,這符合壓致相變中壓力下的配位原則,即高壓相的配位數(shù)(B2結(jié)構(gòu)的配位數(shù)為8)大于等于低壓相的配位數(shù)(B1結(jié)構(gòu)的配位數(shù)為6);B2相是MgO的高壓相,屬體心立方結(jié)構(gòu),空間群為無實驗晶格參數(shù),模擬得到的晶格參數(shù)為2.660 ?;其他各相在0—800 GPa的壓力范圍內(nèi)都不能夠穩(wěn)定存在.

圖3 計算得到的MgO晶體在零溫下的B1,B2,B3,B4,B81,B82和P3m1各可能結(jié)構(gòu)的相對焓隨壓力的變化Fig.3.Calculated relative enthalpies of MgO with B1,B2,B3,B4,B81,B82andP3m1 phases as a function of pressure at zero temperature.

聲子譜,即格波的角頻率與波矢量的關(guān)系曲線,通過對晶體材料聲子譜研究可以明確材料是否具有動力學(xué)穩(wěn)定性特征;同時,根據(jù)聲子譜中譜線重疊以及各原子聲子態(tài)密度峰的位置可反映出原子間是否存在相似或相同的振動狀態(tài),從而能夠推斷出原子間是否存在相互作用以及作用力的強弱.為了檢驗本文提出的MgO晶體7種可能結(jié)構(gòu)的動力學(xué)穩(wěn)定性,在DFT框架下,采用線性響應(yīng)方法[30]補充計算了零壓下B1,B2,B3,B4,B81,B82和P3m1結(jié)構(gòu)的聲子譜和B2結(jié)構(gòu)在相變壓力580 GPa下的聲子譜,如圖4所示.可以看出,零壓下MgO晶體B2結(jié)構(gòu)的聲子譜有虛頻存在,說明該結(jié)構(gòu)在這種狀態(tài)下是不穩(wěn)定的,當(dāng)壓力達到580 GPa時,B2結(jié)構(gòu)成為穩(wěn)定結(jié)構(gòu),符合前面熱力學(xué)穩(wěn)定性計算結(jié)果;零壓下B3,B4,B81,B82和P3m1結(jié)構(gòu)的聲子譜在整個布里淵區(qū)均沒有虛頻出現(xiàn),這意味著它們在零壓下是動力學(xué)穩(wěn)定的,是MgO的亞穩(wěn)結(jié)構(gòu).

3.2 高溫結(jié)構(gòu)穩(wěn)定性

固體發(fā)生熔化時,固、液兩相的吉布斯自由能相等,而體積發(fā)生突變,且伴隨著熵增,因此,屬于一級相變.在熔化過程中,其抗剪切能力消失,在結(jié)構(gòu)上由長程有序結(jié)構(gòu)變?yōu)闊o序結(jié)構(gòu),固、液兩相之間在晶體學(xué)上沒有任何明確的位向關(guān)系,因此屬于重構(gòu)性相變.物質(zhì)在極端壓力條件下的熔化行為是一種非常復(fù)雜的物理過程,在寬廣的溫度和壓強范圍內(nèi),要對物質(zhì)熔化等熱力學(xué)性質(zhì)做出合理描述,將涉及復(fù)雜的離子間的相互作用問題,選擇考慮極化效應(yīng)的殼層模型和在考慮極化效應(yīng)基礎(chǔ)上考慮壓縮效應(yīng)的呼吸殼層模型,模擬了零壓下巖鹽結(jié)構(gòu)的MgO在500—5000 K溫度范圍內(nèi)的體系的摩爾體積和總能隨溫度的變化,如圖5和圖6所示.

由圖5和圖6可以看出,巖鹽結(jié)構(gòu)的MgO在加熱到一定溫度時,對應(yīng)的摩爾體積和總能發(fā)生了明顯的躍變：利用呼吸殼層模型模擬得到的躍變溫度為4490 K,利用殼層模型分子動力學(xué)且結(jié)合SM-SS和SM-LC勢參數(shù)模擬得到的躍變溫度分別為4495和3894 K,利用BG作用勢參數(shù)模擬得到的躍變溫度為3796 K,其中利用SM-LC和BG勢參數(shù)模擬得到的結(jié)果非常接近.熔化是一個動力學(xué)過程,根據(jù)現(xiàn)代熔化理論,晶體的熔化溫度可根據(jù)一定的熔化機制來進行修正[46],即過熱率可按下式來修正：

圖4 計算得到的MgO晶體B1 (a),B3 (b),B4 (c),B81(d),B82(e),P3m1 (f)和B2 (g)結(jié)構(gòu)在零壓下的聲子譜和B2結(jié)構(gòu)在相變壓力為580 GPa下的聲子譜(h)Fig.4.Calculated phonon spectra of MgO with B1 (a),B3 (b),B4 (c),B81(d),B82(e),P3m1 (f)and B2 (g)phases at 0 GPa and with B2 phase at 580 GPa (h).

式中,T為模擬觀測到的溫度;Tm為材料實際熔化溫度.我們曾利用單相分子動力學(xué)模擬方法將0.1 MPa,31.5 GPa,47 GPa和135 GPa壓力下的分子動力學(xué)熱不穩(wěn)定性溫度和實驗數(shù)據(jù)比較來確定過熱率,得到的不同壓力下的過熱率基本一致[15].事實上,過熱熔化跟升溫速率有很大的關(guān)系,而壓力對其影響不大,可據(jù)此來修正晶體過熱熔化溫度,本文模擬的加溫間隔為100 K.和實驗熔化溫度3083 K進行比較[47],得到呼吸殼層模型模擬的過熱率為0.456,殼層模型分子動力學(xué)結(jié)合SMSS和SM-LC勢參數(shù)模擬得到的過熱率分別為0.458和0.263,BG作用勢參數(shù)模擬的過熱率為0.231.根據(jù)晶體熔化均勻形核理論[46]和動力學(xué)災(zāi)變理論[48]定義的過熱極限來看,過熱發(fā)生在一個很大的范圍之內(nèi),即：θ在0.1—2.0范圍,本文得到的過熱率符合該過熱極限.從這個意義上來講,本文選取的殼層和呼吸殼層模型均可用來研究MgO的高壓熔化特性和高溫結(jié)構(gòu)穩(wěn)定性.

圖5 分子動力學(xué)模擬得到的零壓下的MgO巖鹽結(jié)構(gòu)的體積隨溫度的變化Fig.5.Molecular dynamics calculated volume of MgO with NaCl-type structure as a function of temperature at zero pressure.

圖7給出了單相分子動力學(xué)模擬得到的巖鹽結(jié)構(gòu)的MgO壓力達150 GPa的熔化相圖,該相圖對模擬得到的熱不穩(wěn)定性溫度結(jié)合過熱率進行了修正,同時和Zerr和Boehler[13]利用DAC技術(shù)得到的實驗數(shù)據(jù)以及Belonoshko和Dubrovinsky[49]的兩相分子動力學(xué)模擬、Wang經(jīng)驗?zāi)Ｐ蚚50]與Lindemann熔化方程[51]得到的結(jié)果做了比較.為了消除“加熱模擬體系直至熔化”方法產(chǎn)生的過熱,Luo等[52]提出了“過熱-過冷循環(huán)”方法,單相方法模擬固體熔化會引起過熱,如果用單相方法來模擬液體結(jié)晶同樣會導(dǎo)致過冷,過熱與過冷的程度相當(dāng),因此如果采用單相模擬的方法,那么過熱與過冷從體積變化圖上將形成一個封閉的環(huán),對過熱與過冷進行適當(dāng)?shù)摹捌骄本涂梢远ǔ鋈刍瘻囟?“過熱-過冷循環(huán)”方法即補充逆向液體結(jié)晶過冷模擬以消除過熱的方法,實際上是對過熱與過冷的一種平均,但其本質(zhì)上也是單相模擬.兩相模擬的主要思想是先構(gòu)建一個固-液共存的系統(tǒng),固體和液體之間存在交界面,然后固定壓力把系統(tǒng)的溫度升到接近熔化溫度,如果此時的溫度低于實際的熔化溫度,那么液體部分就要逐漸結(jié)晶,固-液交界面就要向液體一邊移動直至最終完全結(jié)晶,如果此時的溫度高于實際的熔化溫度,那么固體部分就要逐漸熔化,固-液界面就要向固體一側(cè)移動直至完全熔化.和兩相模擬相比,單相模擬方法本質(zhì)上是一種“加熱模擬體系直至熔化”的模擬方法,但是這種方法由于加熱比較快,所以容易引起過熱現(xiàn)象,即晶體在高于自身熔點的溫度下發(fā)生熔化的現(xiàn)象.這里針對呼吸殼層模型、SM-SS和SM-LC殼層模型、BG作用勢參數(shù)模型計算的熱不穩(wěn)定性溫度采用的修正過熱率分別為0.456,0.458,0.263和0.231.由圖7可以看出,本文呼吸殼層模型模擬結(jié)果與Belonoshko和Dubrovinsky[49]的兩相分子動力學(xué)模擬結(jié)果符合,殼層模型和BG作用勢參數(shù)模型模擬結(jié)果和Wang經(jīng)驗?zāi)Ｐ蚚50]預(yù)測熔化溫度符合程度達到140 GPa.Zerr和Boehler[13]的DAC實驗結(jié)果明顯偏低于所有理論結(jié)果,這是由于表面效應(yīng)的原因造成的,與大塊晶體相比,細微顆粒的熔點要低得多,原因是細微顆粒具有大的比表面積和表面能,在溫度很低的情況下就能在表面發(fā)生熔化,而在DAC實驗中觀察到的熔化溫度剛好是發(fā)生表面熔化時的溫度.

圖6 分子動力學(xué)模擬得到的零壓下的MgO巖鹽結(jié)構(gòu)的總能隨溫度的變化Fig.6.Molecular dynamics calculated total energy of MgO with NaCl-type structure as a function of temperature at zero pressure.

圖7 分子動力學(xué)模擬得到的MgO巖鹽結(jié)構(gòu)壓力達150 GPa的熔化相圖Fig.7.Molecular dynamics calculated melting phase diagram of MgO with NaCl-type structure as a function of pressure up to 150 GPa.

4 結(jié) 論

本文利用基于DFT的第一性原理計算方法,對MgO各種可能存在的結(jié)構(gòu)B1,B2,B3,B4,B81及基于粒子群優(yōu)化算法預(yù)測的晶體結(jié)構(gòu)B82和P3m1的穩(wěn)定性在0—800 GPa的壓力范圍內(nèi)進行了系統(tǒng)研究,發(fā)現(xiàn)MgO在0—580 GPa的壓力范圍內(nèi)一直以穩(wěn)定巖鹽結(jié)構(gòu)存在,隨著壓力增加到580 GPa,MgO會發(fā)生從B1到B2結(jié)構(gòu)的相轉(zhuǎn)變.盡管B81結(jié)構(gòu)能合理解釋沖擊壓縮實驗中沿MgO的P-V雨貢紐線在(170 ± 10)GPa存在體積不連續(xù)的原因,B4結(jié)構(gòu)亦能合理解釋高壓下MgO理論計算的熔化線和實驗結(jié)果相差很大的原因,但聲子譜計算表明MgO晶體推測的B81和B4結(jié)構(gòu)及其他結(jié)構(gòu)在零壓下僅為MgO晶體的亞穩(wěn)結(jié)構(gòu).在確定MgO固-固相變物理圖景的基礎(chǔ)上,采用經(jīng)典分子動力學(xué)方法,分別引入能很好描述離子極化特性的殼層模型和離子壓縮效應(yīng)的呼吸殼層模型,對MgO巖鹽結(jié)構(gòu)的高溫穩(wěn)定性進行了模擬研究,發(fā)現(xiàn)利用呼吸殼層模型模擬得到的熱不穩(wěn)定性溫度為4490 K,利用殼層模型分子動力學(xué)且結(jié)合Stoneham和Sangster[34]給出的SMSS勢參數(shù)及Lewis和Catlow[33]給出的SM-LC勢參數(shù)模擬得到的熱不穩(wěn)定性溫度分別為4495和3894 K.和實驗熔化溫度相比,呼吸殼層模型模擬的過熱率為0.456,殼層模型分子動力學(xué)結(jié)合SMSS和SM-LC勢參數(shù)模擬得到的過熱率分別為0.458和0.263,符合根據(jù)晶體熔化均勻形核理論和動力學(xué)災(zāi)變理論定義的0.1—2.0的過熱極限.MgO晶體巖鹽結(jié)構(gòu)壓力達150 GPa的熔化相圖計算表明利用呼吸殼層模型、殼層模型模擬得到的結(jié)果和文獻給出的兩相分子動力學(xué)模擬結(jié)果符合,DAC實驗結(jié)果明顯偏低于理論結(jié)果的原因是由于表面效應(yīng)造成,即在實驗中觀察到的熔化溫度是發(fā)生表面熔化時的溫度而不是整個塊體材料的真正熔化溫度.