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

宋婷 孫小偉 魏小平 歐陽(yáng)玉花 張春林 郭鵬 趙煒

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

1 引 言

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

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

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

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

2 計(jì)算方法

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

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

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

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

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

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

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

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

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

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

表1 巖鹽結(jié)構(gòu)的MgO晶體特性模擬的短程對(duì)勢(shì)參數(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 計(jì)算方法可靠性檢驗(yàn)

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

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

3 分析和討論

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

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

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

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

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

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

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

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

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

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