范俊宇鄭朝陽蘇艷?趙紀(jì)軍

1)(大連理工大學(xué),三束材料改性教育部重點(diǎn)實驗室,大連 116024)2)(中國工程物理研究院流體物理研究所,沖擊波與爆轟物理國家實驗室,綿陽 621900)(2016年10月9日收到;2016年11月21日收到修改稿)

專題:高壓下物質(zhì)的新結(jié)構(gòu)與新性質(zhì)研究進(jìn)展

單軸壓縮下固態(tài)硝基苯的第一性原理研究?

范俊宇1)鄭朝陽2)蘇艷1)?趙紀(jì)軍1)

1)(大連理工大學(xué),三束材料改性教育部重點(diǎn)實驗室,大連 116024)2)(中國工程物理研究院流體物理研究所,沖擊波與爆轟物理國家實驗室,綿陽 621900)(2016年10月9日收到;2016年11月21日收到修改稿)

采用第一性原理密度泛函理論結(jié)合經(jīng)典色散修正方法,對固態(tài)硝基苯在單軸壓縮下的基本結(jié)構(gòu)關(guān)系進(jìn)行了計算.靜水壓縮和單軸壓縮都壓縮到初始平衡體積的70%.將靜水壓下優(yōu)化后的晶胞體積、晶格參數(shù)以及平衡條件下的晶格能與實驗值進(jìn)行了比較,均符合較好.同時,為了充分地表征固態(tài)硝基苯的各向異性,將硝基苯沿著三個晶格矢量的方向進(jìn)行單軸壓縮,把每個方向的應(yīng)力張量、能帶帶隙、每個原子能量的改變分別作為體積壓縮比的函數(shù)進(jìn)行了比較和分析.其中,最顯著的各向異性效應(yīng)是在體積壓縮比為0.76時,沿X軸壓縮導(dǎo)致硝基苯能帶帶隙閉合,體系呈金屬化;而靜水壓縮或沿Y軸和Z軸壓縮時體系始終呈半導(dǎo)體狀態(tài),帶隙均大于1.59 eV.為了充分理解這一各向異性特性,我們計算了硝基苯晶體的局域態(tài)密度和電荷密度分布,并對金屬化現(xiàn)象做出了合理的分析和解釋.在不同的壓力加載條件下,通過對不同物理量的計算,發(fā)現(xiàn)X軸方向是硝基苯晶體內(nèi)部最敏感的方向.這些各向異性效應(yīng)的研究將有助于人們在原子尺度上深入理解固態(tài)硝基苯的物理化學(xué)性質(zhì).

硝基苯,單軸壓縮,靜水壓,含能材料

1 引 言

含能材料是一種運(yùn)用于國防工業(yè)領(lǐng)域的重要材料,其自身存儲著較高的化學(xué)能量,受到外界擾動可以快速地釋放熱量或氣態(tài)產(chǎn)物.一般來說,含能材料分為炸藥、推進(jìn)劑和煙火劑等.在過去的數(shù)十年里,許多研究團(tuán)隊對含能材料的研究做了非常多的努力和貢獻(xiàn)[1?4].充分理解材料的宏觀特性與微觀結(jié)構(gòu)變化之間的關(guān)系是含能材料的研究熱點(diǎn),尤其是爆轟過程以及爆轟的初始階段對于含能材料的研究至關(guān)重要.明確材料的爆轟行為與分子結(jié)構(gòu)的關(guān)系,可以使我們更加深入地了解含能材料的性能,從而有效地設(shè)計出新型高效的炸藥[5].然而,由于含能爆炸物實驗研究的危險性以及實驗數(shù)據(jù)的缺乏,理論模擬變得尤為重要,并且可以在一定程度上指導(dǎo)和輔助實驗的進(jìn)行.典型的含能材料通常是有機(jī)分子晶體,例如硝基甲烷、黑索金、FOX-7等[6?9].人們運(yùn)用密度泛函理論,在原子和分子尺度上研究含能有機(jī)分子晶體在高壓環(huán)境下的物理和化學(xué)性質(zhì),這對理解其高壓行為(包括靜水壓縮和單軸壓縮)有很大的幫助.

固態(tài)硝基苯(其化學(xué)式為C6H5NO2)通常作為含能分子晶體的原型體系而受到關(guān)注和研究,并在實驗和理論方面均取得了一定的進(jìn)展[10?17].實驗上,通過X-射線衍射[10,11]、微波[12]、電子衍射[13]等技術(shù)測定硝基苯的晶體結(jié)構(gòu),如Boese等[10]和Trotter[11]先后在不同的溫度條件下測定了常壓下硝基苯的晶體結(jié)構(gòu).如圖1所示,硝基苯的晶體結(jié)構(gòu)屬于P 21/c空間群、單斜晶系,且每個晶胞內(nèi)包含4個硝基苯分子.在103 K時,硝基苯的實驗晶格參數(shù)是:a=3.8014?,b=11.6153?,c=12.9843?,β=94.98?.由于硝基苯分子間存在弱相互作用,鄰近的氫原子和氧原子能夠形成氫鍵,使得硝基苯晶體形成與三氨基三硝基苯(TATB)晶體類似的層狀結(jié)構(gòu)[10,18,19].高壓下,Kozu等[16]通過爆炸平面波生成器和金剛石對頂砧裝置測得硝基苯的雨貢紐曲線和靜水壓縮的PV關(guān)系.Kobayashi和Sekine[17]利用單脈沖激光拉曼系統(tǒng)測定了硝基苯晶體的實時拉曼光譜,認(rèn)為分子間作用力使NO2伸縮振動在沖擊壓縮下有較小的藍(lán)移.

圖1 硝基苯分子結(jié)構(gòu)(左)和晶體結(jié)構(gòu)(右)圖Fig.1.The molecu lar structu re(left)and crystal structu re(right)of nitrobenzene.

在理論研究方面,Domenicano等[14]通過分子軌道方法研究了硝基苯分子的非平面構(gòu)型.而Clarkson和Smith[15]采用B3 LYP/6-311+G**計算方案研究了硝基苯晶體的拉曼和紅外振動光譜.Wang等[20]運(yùn)用色散修正的密度泛函理論(DFT)計算了硝基苯在靜水壓0—10 GPa范圍內(nèi)分子結(jié)構(gòu)的變化和拉曼振動模式.雖然有關(guān)硝基苯晶體的理論工作較少,但是為硝基苯的研究工作打下了基礎(chǔ).

固態(tài)硝基苯可以用作生產(chǎn)炸藥的原料,而液態(tài)硝基苯往往是工業(yè)污水的成分之一,因此這兩種狀態(tài)的硝基苯的分解機(jī)理和過程被廣泛關(guān)注[21?24].Pruitt和Goebbert[21]結(jié)合質(zhì)譜分析方法和密度泛函理論計算,發(fā)現(xiàn)固態(tài)硝基苯等芳香族化合物π軌道上額外電子的作用,導(dǎo)致陰離子的C—N解離能比中性的C—N解離能要低.Dong和Sang[24]利用連續(xù)式全混流型反應(yīng)系統(tǒng)觀測到硝基苯在超臨界水環(huán)境下能夠首先分解為苯和硝酸鹽,而硝基和超臨界水的氧化作用使得反應(yīng)進(jìn)一步生成碳的一氧化物和二氧化物.

盡管這些理論和實驗工作對硝基苯晶體的研究做了很大的努力,但在高壓下固態(tài)硝基苯的結(jié)構(gòu)轉(zhuǎn)變以及應(yīng)力分布等理論研究卻很少.含能分子晶體單軸壓縮[25?27]的研究可以得到晶體沿著不同方向的物理化學(xué)特性,所以對固態(tài)硝基苯晶體單軸壓縮行為的研究是很有必要的.由于固態(tài)硝基苯可以看作是準(zhǔn)正交的晶體結(jié)構(gòu),在本文中我們運(yùn)用色散能量修正下的密度泛函理論[28,29]對單軸壓縮下的硝基苯晶體進(jìn)行了系統(tǒng)的計算研究.

2 計算方法

本文中的DFT計算[30]采用CASTEP(Cambridge Sequential Total Energy Package)程序[31]來完成. 離子和電子相互作用采用平面波基組下的模守恒贗勢來近似描述,廣義梯度近似下的Perdew-Burke-Ernzerhof(PBE)泛函[32]用于處理體系的交換關(guān)聯(lián)勢能.對于分子間的非共價弱相互作用的處理,采用Grimme[28]提出的經(jīng)驗的能量色散修正方案,即DFT-D2.運(yùn)用Broyden-Fletcher-Gold farb-Shanno算法對晶體結(jié)構(gòu)充分弛豫,其中設(shè)定每個原子的能量收斂標(biāo)準(zhǔn)是5.0×10?6eV,最大的力收斂標(biāo)準(zhǔn)是0.01 eV/?,最大的晶胞應(yīng)力收斂標(biāo)準(zhǔn)是0.02 GPa.體系能量截斷能和布里淵區(qū)k點(diǎn)取樣是決定控制DFT計算的關(guān)鍵參數(shù)[33],所以我們進(jìn)行了一系列的收斂性測試來保證計算的數(shù)值精度,發(fā)現(xiàn)截斷能取700 eV和k點(diǎn)平均間隔取0.05??1時,體系總能保持穩(wěn)定,其收斂精度為1.6meV/atom.

初始的晶胞結(jié)構(gòu)根據(jù)Boese等[10]的實驗數(shù)據(jù)搭建,在對晶胞形狀、晶格參數(shù)和原子坐標(biāo)沒有約束的條件下,對晶體結(jié)構(gòu)進(jìn)行了充分地弛豫,從而計算出0 K時平衡的晶體結(jié)構(gòu).其中,靜水壓縮時,外界對晶體結(jié)構(gòu)從0 GPa開始持續(xù)加壓直到20 GPa,且在壓縮過程中,晶胞的體積、形狀、晶胞中原子的坐標(biāo)可以充分地弛豫.對于單軸壓縮,我們保持單斜晶系和β角不變,分別只沿著三個坐標(biāo)軸的方向之一進(jìn)行單軸壓縮(X軸對應(yīng)晶格矢量a,Y軸對應(yīng)晶格矢量b,Z軸對應(yīng)晶格矢量c),在整個壓縮過程中,每一步的壓縮都保持2%的比例,直至達(dá)到其平衡體積的70%.單軸壓縮時,固定晶胞形狀,僅弛豫原子的坐標(biāo).

3 結(jié)果與討論

3.1 靜水壓縮下硝基苯的結(jié)構(gòu)變化

表1對計算得到的固態(tài)硝基苯的晶格參數(shù)和晶格能與實驗數(shù)據(jù)做了比較.通過晶格參數(shù)的測定,不僅可以了解分子晶體的幾何結(jié)構(gòu),也可以由此來判定所采用的理論方法對分子間相互作用描述的可靠程度.初始的硝基苯晶胞運(yùn)用PBE泛函加色散修正方法得到充分的優(yōu)化,與實驗數(shù)據(jù)[10]進(jìn)行比較,其中晶格參數(shù)a,b,c的誤差分別是?2.98%,0.01%,?1.67%,而晶胞體積的誤差是?4.35%.注意到實驗是在103 K溫度下,而我們的計算溫度是0 K,因此推測溫度效應(yīng)是導(dǎo)致常壓下晶格參數(shù)低估的原因之一.從計算結(jié)果來看,DFT-D方案低估了常壓下的晶格參數(shù).

體系的晶格能可以用來表征分子晶體中分子間相互作用的強(qiáng)度,其定義為氣態(tài)分子與固態(tài)晶體之間的能量差.在本文中,將一個孤立的硝基苯分子放置在15?×15?×15?的立方超元胞內(nèi)模擬氣態(tài)分子,計算得到的硝基苯晶格能是0.7243 eV,這與Caillet和Claverie[34]實驗測得的硝基苯晶體晶格能在0.6635—0.6852 eV之間相比,高估了5.71%.結(jié)合對晶格參數(shù)的描述,我們認(rèn)為DFT-D在常壓下對分子間相互作用的描述并不充分,導(dǎo)致計算的晶格參數(shù)和晶格能與實驗值比較有不同程度的誤差.

為了考察硝基苯晶體在高壓下結(jié)構(gòu)的變化并檢驗色散修正方案在高壓下描述晶體結(jié)構(gòu)的可靠性,圖2給出了理論計算和實驗[16]的晶胞體積隨壓力的變化.從圖2可以看出,隨著壓力的增大,晶胞體積單調(diào)遞減.整體上,本文的計算很好地再現(xiàn)了體積隨壓力變化的趨勢.然而,在給定的壓力下,計算的晶胞體積比實驗數(shù)據(jù)略小,推測這是由于DFT-D框架下對分子間作用力的描述不充分和上面提到的溫度效應(yīng)所導(dǎo)致的.但是,我們注意到隨著壓力的增大,相同壓力下理論值與實驗值更加接近.例如在0 GPa時,體積的偏差是?4.35%,而在接近4 GPa時,偏差約為?1.93%,這也意味著高壓下DFT-D對分子間相互作用力的描述比較可靠.

圖2 晶胞體積隨壓力的變化Fig.2.The relationshipbetween the systempressure and unit cell volume.

為了檢驗高壓下晶格參數(shù)的變化是否具有可靠性,圖3給出了理論和實驗的晶格參數(shù)隨壓力的變化,發(fā)現(xiàn)晶格參數(shù)a,b,c隨壓力的變化趨勢與實驗[16]保持了較好的一致.同時,本文的結(jié)果也同Wang等[20]在DFT-D框架下對硝基苯晶體在靜水壓縮下的計算結(jié)果基本一致.在5 GPa壓力下,晶格參數(shù)a,b,c的體積壓縮比分別為10.01%,3.30%,6.74%,可見a方向的體積壓縮比明顯大于b和c方向的體積壓縮比.這主要源于硝基苯晶體中分子間氫鍵網(wǎng)絡(luò)主要沿著晶格參數(shù)b和c所在平面內(nèi),導(dǎo)致氫鍵網(wǎng)絡(luò)的層內(nèi)b,c方向不易壓縮,而層間方向a易壓縮.

有機(jī)分子晶體往往存在分子間的氫鍵相互作用,雖然這種相互作用的強(qiáng)度比離子鍵、共價鍵要低,但卻在凝聚態(tài)物質(zhì)如層狀固體、分子晶體中起到主導(dǎo)作用.圖1是硝基苯分子的結(jié)構(gòu)圖,根據(jù)第一性原理計算,硝基苯晶體中C—H...O氫鍵的存在會對晶體的分子構(gòu)型產(chǎn)生一定的影響.我們注意到以下三點(diǎn):

表1 零壓下晶格參數(shù)和晶格能與實驗數(shù)據(jù)的比較Tab le 1.Comparison of experimental data with lattice parameter and lattice energy under the cond ition of zeropressure.

圖3 壓力依賴下的晶格參數(shù)的變化關(guān)系以及實驗數(shù)據(jù)與理論值的比較Fig.3.Pressure dependence of lattice constants for Nitrobenzene.Experimental resu lts compared with present theoretical resu lts.

1)在氣態(tài)硝基苯分子中,C—C鍵的長度關(guān)于C—N鍵對稱分布,C1—C2,C2—C3,C5—C6和C6—C1鍵長均為1.396 ?,C3—C4和C4—C5鍵長均為1.400?;當(dāng)孤立的硝基苯分子形成晶體時,6個C—C鍵的鍵長發(fā)生了對稱性的伸縮,C1—C2,C2—C3,C5—C6鍵長均為1.399?,C6—C1鍵長為1.398?,C3—C4和C4—C5鍵長均為1.393?;

2)孤立硝基苯分子結(jié)合成晶體時,5個C—H鍵都一致地縮短了0.001—0.002?,而C—N鍵的縮短最為明顯,達(dá)到0.018?;與此相反,兩個N—O鍵卻均伸長了0.006?;

3)O=N鍵和O—N鍵向內(nèi)收縮使得O=N—O角由124.85?減小至123.35?.

因此,從硝基苯形成晶體的分子幾何構(gòu)型變化可以看出,分子間相互作用起到了關(guān)鍵作用,不能忽視.

3.2 單軸壓縮下硝基苯的結(jié)構(gòu)性質(zhì)

圖4 應(yīng)力張量本征值σxx,σyy,σzz(分別對應(yīng)a,b,c晶格矢量)隨體積壓縮比的變化Fig.4.Relationshipbetween stress tensor eigenvalues σxx,σyy,σzz(corresponding tolattice vectors of a,b,c,respectively)and compression ratioV/V0.

為了模擬單軸壓縮對硝基苯結(jié)構(gòu)的影響,初始的晶格參數(shù)按照固定的比例進(jìn)行了壓縮.根據(jù)文獻(xiàn)中研究單斜晶系應(yīng)力張量的方案[25?27],我們對每個方向和每一步的單軸壓縮時的主應(yīng)力,即應(yīng)力張量的本征值進(jìn)行了測定.將應(yīng)力張量的三個本征值σxx,σyy,σzz依次排序,并且標(biāo)記為σ1,σ2,σ3.圖4給出了每個單軸壓縮方向上應(yīng)力張量的三個主應(yīng)力σ1,σ2,σ3隨體積壓縮比V/V0(取值范圍在1—0.7)的變化.同時,靜水壓縮下的壓力與體積壓縮比的關(guān)系作為對比也呈現(xiàn)在圖4中.對于不同的壓縮方向,應(yīng)力張量的三個本征值隨體積壓縮比表現(xiàn)出不同的變化趨勢.三個主應(yīng)力都隨著體積壓縮比的減小而增大,但是不同的變化幅度和趨勢表現(xiàn)出硝基苯晶體各向異性的特征.其中,沿著Y軸和Z軸方向,主應(yīng)力值σ1,σ2,σ3在不同方向之間的差值保持在1.2—6.4 GPa區(qū)間內(nèi).沿X軸方向的變化呈現(xiàn)出更加明顯的可壓縮性,在體積壓縮比V/V0大于0.94時,主應(yīng)力σ1沿X軸方向的增長率明顯高于沿Y軸和Z軸方向的增長率,并且在體積壓縮比V/V0等于0.7時,沿X軸方向壓縮得到的主應(yīng)力值σ1高達(dá)40.95 GPa,相當(dāng)于沿Y軸(16.42 GPa)和Z軸(17.62 GPa)所得的主應(yīng)力值的兩倍之多.而對于主應(yīng)力值σ2和σ3,沿X 軸壓縮所引發(fā)的應(yīng)力變化趨勢卻大相徑庭,在體積壓縮比V/V0等于0.7時,主應(yīng)力值σ2沿X軸壓縮的值是4.40 GPa,遠(yuǎn)遠(yuǎn)低于沿Y軸和Z軸對應(yīng)的應(yīng)力值.與σ2類似,主應(yīng)力值σ3在體積壓縮比V/V0達(dá)到0.7時僅為5.52 GPa,對于不同的壓縮強(qiáng)度,σ3值的變化并不明顯.綜上所述,我們認(rèn)為硝基苯晶體中,沿氫鍵網(wǎng)絡(luò)的層間方向,即X軸方向,相對于Y軸和Z軸方向?qū)τ趹?yīng)力的變化更加敏感;且靜水壓縮作為單軸壓縮的對比,靜水壓力都略低于沿Y軸和Z軸壓縮的壓力本征值.

在不同的壓力加載條件下,對每個原子的能量變化也進(jìn)行了計算.圖5給出了沿著不同方向單軸壓縮以及靜水壓縮時,平均每個原子能量隨體積壓縮比的變化.在圖5中,同樣可以觀察到能量變化在三個方向上的各向異性行為.隨著體積壓縮比V/V0的增大,每個方向上原子的能量變化幅度也同時增大,并且由于單軸壓縮方向晶格的束縛,單軸壓縮時原子能量的變化總是大于靜水壓縮時能量的變化幅度.在體積壓縮比V/V0等于0.86時,沿X軸方向原子的能量變化開始明顯高于沿Y軸和Z軸方向,在V/V0為0.7時,能量變化率達(dá)到0.24 eV/atom.沿Y軸方向的能量變化曲線上,體積壓縮比在0.8—0.78之間,能量出現(xiàn)微小的波動,即在體積壓縮比0.8—0.78之間,每個原子的能量變化率增大.我們注意到,在應(yīng)力隨體積壓縮比的變化中同樣可以觀測到類似的變化.

與此同時,我們還考察了硝基苯分子結(jié)構(gòu)的變化,如圖6所示.可以發(fā)現(xiàn)在體積壓縮比V/V0為0.78—0.8時,硝基苯分子C—N鍵鍵長的變化幅度明顯大于其他壓縮區(qū)間,我們認(rèn)為硝基苯分子C—N鍵鍵長的變化與每個原子能量的變化以及應(yīng)力的變化有密切的關(guān)系,C—N鍵的鍵長幅度變化增大導(dǎo)致每個原子能量和應(yīng)力及其變化率增大.在過去對含能分子晶體的理論研究中,比如TATB[18]和硝基甲烷[35],C—N鍵對壓力的變化也比較敏感.因此我們推測在外界加載條件下,C—N鍵是影響硝基苯晶體物理化學(xué)性質(zhì)變化的關(guān)鍵因素.

圖5 單軸壓縮與靜水壓縮下,每個原子的能量變化Fig.5.Change in energy per atomupon uniaxial and hyd rostatic compressions.

圖6 C—N鍵長隨體積壓縮比的變化Fig.6.Change of C—N bond length under d iff erent kinds of compression.

能帶帶隙是描述晶體電子結(jié)構(gòu)的重要物理量之一.圖7給出了不同方向的單軸壓縮與靜水壓縮下,硝基苯的能帶帶隙隨體積壓縮比的變化.相比于靜水壓縮和沿X軸方向壓縮,沿著Y軸和Z軸壓縮時體系的能帶帶隙緩慢地減小,而在體積壓縮比V/V0大于0.8時,帶隙的減小速率明顯加快,并且在體積壓縮比V/V0達(dá)到0.7時,帶隙寬度趨于一致.對于沿X軸方向壓縮,帶隙的變化趨勢表現(xiàn)出特有的規(guī)律,體積壓縮比V/V0在1—0.76區(qū)間時,帶隙寬度單調(diào)遞減,減小率遠(yuǎn)遠(yuǎn)大于沿Y軸、Z軸壓縮和靜水壓縮.在整個壓縮區(qū)間內(nèi),沿X軸帶隙寬度變化了2.25 eV,而沿Y軸、Z軸和靜水壓縮僅僅變化了0.77,0.54和0.86 eV.值得注意的是,沿X軸壓縮并且體積壓縮比V/V0達(dá)到0.76,對應(yīng)的體系平均壓力為9.27 GPa時,體系的能帶帶隙閉合,硝基苯晶體由半導(dǎo)體轉(zhuǎn)變?yōu)榻饘?這是由于硝基苯晶體中X軸方向是氫鍵網(wǎng)絡(luò)的層間方向,因此在X軸方向容易壓縮使體系金屬化.

為了進(jìn)一步討論金屬化的晶體結(jié)構(gòu),圖8給出了發(fā)生金屬化轉(zhuǎn)變時硝基苯晶體的總態(tài)密度和分態(tài)密度.從總態(tài)密度圖可以看到:在費(fèi)米能級附近,價帶頂?shù)姆鍋碓从趐態(tài)電子的貢獻(xiàn),而價帶頂其他的峰值來源于s態(tài)和p態(tài)電子的貢獻(xiàn),同時導(dǎo)帶底的峰值也主要來自p態(tài)電子的貢獻(xiàn).Cui等[36]在β-HMX(奧克托今)晶體電子結(jié)構(gòu)的研究中也得到了類似的結(jié)論.在分態(tài)密度圖中,我們可以看出在價帶頂和導(dǎo)帶底的峰值主要來自C,O,N原子p態(tài)電子的作用.值得注意的是,C原子s態(tài)和p態(tài)電子對所有峰值的貢獻(xiàn)都大于其他原子,特別是費(fèi)米能級上的電子態(tài)主要來自于C原子的p態(tài).我們注意到,在零壓平衡條件(即V/V0=1)下,硝基苯層間距離(即沿著a方向的近鄰苯環(huán)之間的距離)為3.326?,當(dāng)V/V0等于0.76,即體系金屬化時,苯環(huán)間的距離縮短為2.680?.所以,我們構(gòu)造了包含層狀結(jié)構(gòu)的超胞來計算V/V0=1和0.76時體系的電荷密度圖,圖9(a)和圖9(b)分別表示硝基苯在常壓下和金屬化時的電荷密度.從圖9可以看出,在結(jié)構(gòu)金屬化時,苯環(huán)上的部分電荷發(fā)生重疊,而硝基上的電荷沒有發(fā)生重疊,也就是說近鄰苯環(huán)上的π電子發(fā)生了一定的交疊.因此我們推斷體系金屬化是由于C原子上離域π電子導(dǎo)致的.

圖7 單軸壓縮與靜水壓縮下能帶帶隙隨體積壓縮比V/V0的變化Fig.7.Changes in band gapupon uniaxial and hyd rostatic compressions as a function of V/V0.

圖8 硝基苯晶體金屬化時的分態(tài)密度圖,虛線表示費(fèi)米能級Fig.8.Partial density of states(PDOS)of nitrobenzene at metallization cond ition.The Fermi energy is shown as a dashed vertical line.

圖9 (a)常壓下體系電荷密度圖;(b)金屬化時體系電荷密度圖Fig.9.(a)Charge density of systemunder ambient cond ition;(b)charge density of systemat metallization.

3.3 硝基苯晶體單軸壓縮下的剪切應(yīng)力

剪切應(yīng)力是引起材料塑性形變的重要驅(qū)動力,也是表征材料或者晶體各向異性的主要物理量之一.對剪切應(yīng)力的計算有利于我們進(jìn)一步地理解含能分子晶體中由沖擊波引起的各向異性特征,可以通過對晶體特定的晶向?qū)嵤﹩屋S壓縮實現(xiàn)沖擊波加載.因此,我們根據(jù)應(yīng)力張量的三個本征值,通過方程τ=(σi?σj)/2來計算體系內(nèi)部不同的剪切應(yīng)力值,其中,τ表示剪切應(yīng)力;i/j,兩者取值范圍都是1,2,3.剪切應(yīng)力τxy可通過i=1和j=2計算得到,同理在i=1,j=3和i=2,j=3時可計算得到τxz和τyz.根據(jù)計算結(jié)果,我們發(fā)現(xiàn)剪切應(yīng)力τxy和τxz在不同的壓縮方向上遠(yuǎn)遠(yuǎn)大于剪切應(yīng)力τyz,使得硝基苯晶體的各向異性特征更加明顯,因此這里我們重點(diǎn)討論τxy和τxz在不同壓縮方向的行為.圖10(a)和圖10(b)中給出了τxy和τxz在不同壓縮方向上隨體積壓縮比V/V0的變化.體積壓縮比在1—0.92范圍內(nèi),剪切力τxy沿X,Y,Z軸三個單軸壓縮方向基本沒有差異并且保持著較低的壓力值,但體積壓縮比V/V0=0.92后,沿X軸方向的剪切力呈指數(shù)型增大,在體積壓縮比V/V0=0.7時,該方向的剪應(yīng)力值高達(dá)18.28 GPa,而沿著Y軸和Z軸壓縮的剪應(yīng)力只有2.08 GPa和0.05 GPa.而對于剪應(yīng)力τxz,其變化趨勢與τxy非常相似,在體積壓縮比V/V0等于0.9以后,沿三個方向壓縮產(chǎn)生的剪應(yīng)力表現(xiàn)出明顯不同的變化趨勢,沿Y軸和Z軸增長幅度依然保持在0—4 GPa內(nèi),體積壓縮比V/V0=0.7時,剪應(yīng)力τxz沿X軸方向壓縮也達(dá)到了17.72 GPa.因此,通過對剪切應(yīng)力的分析,我們進(jìn)一步理解了硝基苯晶體的各向異性,同時可以再次預(yù)言X軸方向是硝基苯晶體中最敏感的方向.

圖10 單軸壓縮下切應(yīng)力隨體積壓縮比的變化 (a)τxy;(b)τxzFig.10.Change of shear stresswith compression ratioupon uniaxial compression:(a) τxy;(b) τxz.

4 結(jié) 論

本文通過密度泛函理論結(jié)合能量修正方案DFT-D2,對硝基苯晶體在靜水壓縮和單軸壓縮下的結(jié)構(gòu)和各向異性特征進(jìn)行了系統(tǒng)的計算研究.在靜水壓縮下,DFT-D方案低估了常壓下的晶胞體積和晶格參數(shù)a,b,c.而在高壓下,計算的晶胞體積與實驗值很接近,比如在4 GPa時誤差在1%以內(nèi),計算結(jié)果比較可靠.同時計算得到的硝基苯晶體的晶格能與實驗得出的晶格能相比,能量范圍略有高估,誤差為5.7%.綜合考慮這些計算得到的物理量及其與實驗的誤差,我們可以認(rèn)定PBE泛函結(jié)合能量修正(DFT-D2)的方法對于描述硝基苯晶體的結(jié)構(gòu)以及隨后的單軸壓縮計算是可靠的.

我們對硝基苯晶體在X,Y,Z晶軸方向進(jìn)行了單軸壓縮研究,在每一個壓縮方向上,計算了相應(yīng)的應(yīng)力張量、剪切應(yīng)力、帶隙寬度、每個原子的能量.計算結(jié)果顯示了硝基苯晶體在單軸壓縮下的各向異性特征.首先,對于應(yīng)力張量的三個本征值σ1,σ2,σ3,沿X軸方向壓縮的σ1值遠(yuǎn)高于沿Y軸和Z軸的應(yīng)力值,而沿X軸方向壓縮得到的σ2和σ3值又遠(yuǎn)小于其他兩個方向.其次,對于體系的帶隙寬度,沿X軸方向壓縮時,在體積壓縮比V/V0達(dá)到0.76時,帶隙閉合,體系呈金屬化,源于近鄰分子中苯環(huán)之間離域π電子軌道的重疊.而對于每個原子的能量變化,同樣也是沿著X軸壓縮能量變化最大.最后,我們對體系中最大的兩個剪切應(yīng)力τxy和τxz進(jìn)行了分析,發(fā)現(xiàn)沿X軸壓縮得到的剪切應(yīng)力值都遠(yuǎn)遠(yuǎn)大于沿Y,Z軸方向.因此,我們推測硝基苯晶體內(nèi)部X軸方向比Y和Z軸方向更加敏感.在實驗上,沿著這三個晶軸方向進(jìn)行沖擊壓縮,將呈現(xiàn)出不同的宏觀物性變化規(guī)律.

[1]Zheng Z Y,ZhaoJ J 2016Chin.Phys.B25 076202

[2]Sikder A,Sikder N 2004J.Hazard.Mater.112 1

[3]Politzer P,Mu rray J S,SeminarioJ M,Lane P,G rice ME,Concha MC 2001J.Mol.Struc.:Theochem573 1

[4]Zheng Z Y,ZhaoJ J 2015Chin.J.High Pressure Phys.29 81(in Chinese)[鄭朝陽,趙紀(jì)軍 2015高壓物理學(xué)報29 81]

[5]Fried L E,Manaa MR,Pagoria P F,Simpson R L 2001Annu.Rev.Mater.Res.31 291

[6]Zhang L,Chen L 2013Acta Phys.Sin.62 138201(in Chinese)[張力,陳朗2013物理學(xué)報 62 138201]

[7]Cheng HP,Dan J K,Huang Z M,Peng H,Chen G H2013Acta Phys.Sin.62 163102(in Chinese)[程和平,但加坤,黃智蒙,彭輝,陳光華2013物理學(xué)報62 163102]

[8]Meng Z R,Zhang W B,Du Y,Shang L P,Deng H2015Acta Phys.Sin.64 073302(in Chinese)[孟增睿,張偉斌,杜宇,尚麗平,鄧琥2015物理學(xué)報64 073302]

[9]Zhang L,Chen L 2014Acta Phys.Sin.63 098105(in Chinese)[張力,陳朗2014物理學(xué)報 63 098105]

[10]Boese R,Bl?ser D,Nussbaumer M,Krygowski TM1992Struct.Chem.3 363

[11]Trotter J 1959Acta Crystallogr.12 884

[12]Larsen N W 2010J.Mol.Struct.963 100

[13]BorisenkoKB,Hargittai I1996J.Mol.Struct.382 171

[14]DomenicanoA,Schu ltz G,Hargittai I,ColapietroM,Portalone G,George P,Bock C W 1989Struct.Chem.1 107

[15]C larkson J,Smith W E 2003J.Mol.Struct.655 413

[16]Kozu N,AraiM,Tamu ra M,Fu jihisa H,AokiK,Yoshida M2000Jpn.J.Appl.Phys.39 4875

[17]Kobayashi T,Sekine T2000Phys.Rev.B62 5281

[18]Liu H,ZhaoJ,Du J,Gong Z,Ji G,W ei D 2007Phys.Lett.A367 383

[19]Chen F,Zhang H,ZhaoF,Li Q l,Qu J Y 2008J.Mol.Struc.:Theochem.864 89

[20]W ang W P,Liu F S,Liu Q J,Liu Z T2016Comput.Theor.Chem.1075 98

[21]Pruitt C J M,Goebbert D J 2013Chem.Phys.Lett.580 21

[22]Fayet G,Joubert L,Rotureau P,AdamoC 2008J.Phys.Chem.A112 4054

[23]Pein BC,Sun Y,D lott D D 2013J.Phys.Chem.A117 6066

[24]Dong S L,Sang D P 1996J.Hazard.Mater.51 67

[25]Con roy M,Oleynik I,Zybin S,W hite C 2008Phys.Rev.B77 094107

[26]Conroy M,Oleynik I,Zybin S,W hite C 2009J.Phys.Chem.A113 3610

[27]Margetis D,Kaxiras E,E lstner M,FrauenheimT,Manaa MR 2002J.Chem.Phys.117 788

[28]G rimme S 2011W ires.Comput.Mol.Sci.1 211

[29]G rimme S,Antony J,Ehrlich S,Krieg H2010J.Chem.Phys.132 154104

[30]Parr R G,Yang W 1995Annu.Rev.Phys.Chem.46 701

[31]Segall M,Lindan P J,Probert MA,Pickard C,HasnipP,Clark S,Payne M2002J.Phys.:Condens.Mat.14 2717

[32]PerdewJP,Burke K,ErnzerhofM1996Phys.Rev.Lett.77 3865

[33]Monkhorst HJ,Pack J D 1976Phys.Rev.B13 5188

[34]Caillet J T,C laverie P 1975Acta Crystallogr.Sec.A31 448

[35]Liu H,ZhaoJ J,W ei D Q,Gong Z Z 2006J.Chem.Phys.124 124501

[36]Cui HL,Ji G F,ZhaoJ J,ZhaoF,Chen X R,Zhang Q M,W ei D Q 2010Mol.Simulat.36 670

PACS:61.50.Ks,61.50.Lt,71.20.–bDOI:10.7498/aps.66.036101

First-principle simu lation of solid n itrobenzene under un iaxial compression?

Fan Jun-Yu1)Zheng Zhao-Yang2)Su Yan1)?ZhaoJi-Jun1)

1)(Key Laboratory ofMaterials Mod ifi cation by Laser,Ion and E lectron Beams,Dalian University of Technology,Ministry of Education,Dalian 116024,China)2)(National Key Laboratory of Shock Wave and Detonation Physics,Institute of Fluid Physics,China Academy of Engineering Physics,Mianyang 621900,China)(Received 9 October 2016;revised manuscript received 21 November 2016)

Energeticmaterials(EMs)including explosives,propellants and pyrotechnics have been widely used for themilitary and many other purposes.Solid nitrobenzene(an organic molecular crystal)cou ld be considered as a prototype of energetic material.Uptonow,numerous studies have been devoted tocrystal structures,spectrumproperties and decomposition mechanisms for solid nitrobenzene experimentally and theoretically.However there has been a lack of the comprehensive understanding of the anisotropic characteristics under diff erent loading conditions.Thus we investigate the hydrostatic and uniaxial compressions along three diff erent lattice directions todetermine this anisotropic eff ect.In this work,the density functional theory calculations are performed based on Cambridge Sequential Total Energy Package(CASTEP)code using normconserving pseudopotentials and a kinetic energy cutoff of 700 eV.The generalized gradient approximation with the Perdew-Burke-Ernzerhof parameterization is used.Monkhorst-Pack k-point meshes with a density of 0.05 ??1are used for Brillouin-zone integration.The empirical dispersion correction by G rimme is taken toaccount for week intermolecular interactions.The hyd rostatic compressions are applied from0 GPa to20 GPa.Cell volume,lattice shape and coordinates of the atoms cou ld be fully relaxed.while uniaxial compression is applied upto70%of the equilibriumcell volume in steps of 2%along their lattice directions respectively.At each compression step,on ly atomic coordinates are allowed torelax,with the lattice fixed.The equilibriumlattice structures under hydrostatic compressions are obtained by full relaxation at 0 Ktemperature.In ambient condition,the calcu lated volume and parameter of the unit cellare underestimated compared with the experimentaldata,and corresponding errors are?2.98%,0.01%,?4.39%,5.71%respectively.In contrast,the calculated lattice energy is overestimated compared with the range of experimental resultswith 5.71%of the error.In high pressure condition,the volume and cell parameter of the unit cell as a function of compression ratioare plotted and compared with the experimental data.The theoretical and experimental values are close with the increase of the pressure,for instant,the error decreases from?4.39%at 0 GPa to?1.93%at 4 GPa.On the other hand,the uniaxial compression is applied along the directions of three lattice vectors.The changes of stress tensor,band gap,energy per atomas a function of compression ratioare alsoplotted and discussed,which can characterize the anisotropic eff ect of solid nitrobenzene.Themost noticeable eff ect of anisotropy in solid nitrobenzene is themetallization at V/V0=0.76 compressed along the X axis,while the solid nitrobenzene under hyd rostatic pressure or other uniaxial compressions uptoV/V0=0.76 remains semiconductor with band gaplarger than 1.591 eV.By analyzing the local density of states and charge density distribution of nitrobenzene crystal,we con fi rmthat the metallization is caused by the overlapof theπelectron frombenzene ring.Through calculating diff erent physical parameters,we find that X axis is themost sensitive direction of nitrobenzene crystal.The studies of anisotropic eff ects are expected toshed light on the physicaland chemical propertiesof solid nitrobenzene on an atomistic scale and provide several insights for experiments.

nitrobenzene,uniaxial compression,hydrostatic pressure,energeticmaterials

10.7498/aps.66.036101

?國防基礎(chǔ)科研核基礎(chǔ)科學(xué)挑戰(zhàn)計劃(批準(zhǔn)號:JCKY 2016212A501)、國家自然科學(xué)基金(批準(zhǔn)號:11674046)、中國博士后科學(xué)基金(批準(zhǔn)號:2016M592704)和大連理工大學(xué)超算中心資助的課題.

?通信作者.E-mail:su.yan@d lu t.edu.cn

*Project supported by the Science Challenging Programof the National Defense Basic Scientifi c Research of China(G rant No.JCKY 2016212A501),the National Natural Science Foundation of China(G rant No.11674046),the China Postdoctoral Science Foundation(G rant No.2016M592704),and the Supercompu ting Center of Dalian University of Technology,China.

?Corresponding au thor.E-mail:su.yan@d lut.edu.cn