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

(蘭州交通大學數理學院,蘭州 730070)

1 引 言

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

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

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

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

2 計算方法

2.1 MgO高壓結構預測及優(yōu)化

本文使用由吉林大學馬琰銘小組開發(fā)的擁有自主知識產權的基于粒子群優(yōu)化算法的CALYPSO軟件包[17],預測MgO的1—4倍分子式組成的晶胞分別在0,135,300和500 GPa壓力下的晶體結構.該軟件只需給定化學配比,就能在特定的壓力和溫度條件下預測材料的基態(tài)和亞穩(wěn)態(tài)結構.這種預測結構的方法已經成功應用在大量實驗合成的單質和化合物的常壓及高壓晶體結構研究中[18].例如,2014年,Li等[19]利用該方法對H2S在10—200 GPa壓力區(qū)間的結構進行預測,提出了兩個能量穩(wěn)定且具有金屬性的新高壓相P-1和Cmca,并首次預言高壓下H2S的高溫超導電性;在這一預測工作的啟發(fā)下,Drozdov等[20]開展了S-H體系的高壓超導實驗研究,發(fā)現S-H體系在高壓下呈現兩個超導態(tài),其中低溫高壓下測得的Tc與Li等[19]的計算數據基本吻合,而室溫退火后測得的Tc達到驚人的203 K;2017年,Peng等[21]利用該方法預言了更多富氫包合物結構高壓高溫超導體的存在.本文對預測得到的結構采用第一性原理計算的VASP軟件包進行幾何優(yōu)化[22],電子和電子之間的交換關聯(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預測結構后,應用CASTEP軟件包對預測結構、實驗結構及其他可能存在的結構進行結構優(yōu)化并計算相關性質[26],這里電子與電子間的交換關聯(lián)能采用2008年由Perdew等[27]修訂的Perdew-Burke-Ernzerhof形式的GGA方法,該形式的交換相關泛函能夠提高密堆固體的平衡特性;芯電子與價電子間的相互作用勢由超軟贗勢[28]實現.為了確定計算收斂性,文中對平面波截斷能和k點采樣進行了收斂性測試.針對MgO的B1,B2,B4等7種不同結構,截斷能均取600 eV,就能保證高壓下的收斂精度.零壓和不同外壓條件下的MgO晶體結構和原子位置優(yōu)化采用Broyden,Fletcher,Goldfarb和Shanno提出并發(fā)展起來的幾何算法[29].聲子譜的計算采用線性響應方法[30].

2.2 分子動力學模擬技術

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

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

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

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

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

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

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

表2 巖鹽結構的MgO晶體特性模擬的殼層及呼吸殼層參數Table 2.Shell and breathing shell parameters of MgO crystal characteristic simulation with NaCltype structure.

2.3 計算方法可靠性檢驗

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

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

3 分析和討論

3.1 高壓結構預測

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

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

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

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

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

3.2 高溫結構穩(wěn)定性

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

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

圖4 計算得到的MgO晶體B1 (a),B3 (b),B4 (c),B81(d),B82(e),P3m1 (f)和B2 (g)結構在零壓下的聲子譜和B2結構在相變壓力為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為材料實際熔化溫度.我們曾利用單相分子動力學模擬方法將0.1 MPa,31.5 GPa,47 GPa和135 GPa壓力下的分子動力學熱不穩(wěn)定性溫度和實驗數據比較來確定過熱率,得到的不同壓力下的過熱率基本一致[15].事實上,過熱熔化跟升溫速率有很大的關系,而壓力對其影響不大,可據此來修正晶體過熱熔化溫度,本文模擬的加溫間隔為100 K.和實驗熔化溫度3083 K進行比較[47],得到呼吸殼層模型模擬的過熱率為0.456,殼層模型分子動力學結合SMSS和SM-LC勢參數模擬得到的過熱率分別為0.458和0.263,BG作用勢參數模擬的過熱率為0.231.根據晶體熔化均勻形核理論[46]和動力學災變理論[48]定義的過熱極限來看,過熱發(fā)生在一個很大的范圍之內,即：θ在0.1—2.0范圍,本文得到的過熱率符合該過熱極限.從這個意義上來講,本文選取的殼層和呼吸殼層模型均可用來研究MgO的高壓熔化特性和高溫結構穩(wěn)定性.

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

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

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

圖7 分子動力學模擬得到的MgO巖鹽結構壓力達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 結 論

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