含能材料三種生成熱計算方法的比較

劉 玥,易平平,林晨晨,劉 超,易小藝,何 飄

(中南大學 化學化工學院,湖南 長沙 410083)

引 言

含能材料的生成熱是評價爆轟性能與穩定性的重要參數,其中,固相生成熱是計算爆速(D)和爆壓(P)必不可少的數據[1-3]。對于一般化合物,可通過查閱官方網站(如NIST網站)及熱力學手冊獲得生成熱數據,或是先測定燃燒熱再進行計算。但是對于含能材料而言,由于其高能量特性,它們的生成熱往往難以實測。因此,采用理論方法高效精確地計算其生成熱,是研究含能材料性質的重要方法,也是影響預測爆轟性能準確性的一大關鍵因素,備受實驗和理論科學家的關注。

隨著計算機技術的發展,應用理論計算準確地預測小分子的生成熱已經成為一種方便易行的方法[4-8]。計算化學主要的理論方法有:量子力學、分子力學、統計動力學和分子力學以及以上方法的組合[9]。量子化學是用量子力學的原理研究原子、分子和晶體的電子層結構、化學鍵理論、分子間作用力、化學反應理論、各種光譜、波譜和電子能譜的理論以及無機和有機化合物、生物大分子和各種功能材料的結構和性能關系的科學[10]。量子化學理論方法可分為從頭算方法( ab initio )[11]、半經驗分子軌道法[12]和密度泛函理論( DFT )方法[13]。從頭算法一般比較耗時,需要大量的內存和磁盤空間[11];半經驗分子軌道法比第一性原理計算快很多,可是如果計算的分子與參數化時使用的分子結構不相近,半經驗方法可能給出完全錯誤的結果[12]。而DFT方法則具有明顯的優勢,它可以直接確定精確的基態能量和電子密度[14],用電子密度取代波函數作為研究的基本量,大大簡化了電子結構的計算,并且可以較為準確地預測分子的幾何構型及熱化學性質[13]。

用量子化學方法計算生成熱的一大難點是計算結果不直接提供生成熱,只給出分子的總能量與熱力學焓值等數據,需借助輔助反應才能求得化合物的生成熱。常用的輔助反應包括原子化能反應法[15]、原子當量反應法[16]、等鍵反應法[17]。運用等鍵反應進行誤差抵消,雖然可以大大減少計算量,但是由于方法構造的不確定,將導致同一等級計算結果的不確定。運用原子化能法可以直接計算分子的生成熱,對于給定的計算方法,計算結果相對固定;但是為了得到準確的數值,往往需要非常高等級的計算方法和計算資源。自1989年John Pople[18]等發表Gaussian-1(簡稱G1)以來,已經有G2、G3和G4等四代方法。高斯理論方法(Gaussian theory)是迄今為止最為廣泛使用的計算化合物能量的方法之一,在缺少實驗值的情況下,其計算結果也被當作標準數值使用,然而高精度方法如G4對于多分子體系而言耗時太過昂貴。2011年,Bun Chan 等[19]提出了G4(MP2)-6X熱力學組合方法,該方法相對于G4(MP2)最主要變化,是把CCSD(T)的CCSD和(T)部分的相關能進行系數校正,并把MP2換成了SCS-MP2,而且給出了可直接輸出H(0)和H(T)數值的腳本,耗時與G4(MP2)相當,而精度接近G4。

本研究設計了160種新型含能分子,并在B3LYP/6-311G(d, p)//MP2/6-311++G(d, p)下用3種輔助反應(原子化能方法、等鍵反應方法與原子當量方法)計算標準摩爾生成熱值,與G4(MP2)-6X方法結合原子化能方法計算的標準摩爾生成熱值做相關性分析,并對4種方法進行了比較研究。

1 計算方法

運用Gaussion 09[20]軟件,在B3YLP/6-311G( d, p)理論水平下,對設計的160個含能分子進行幾何優化,并進行振動頻率計算,獲得零點振動能ZPE,所得結構均證實為勢能面上的極小點(無虛頻);在MP2/6-311++G(d, p) 理論水平下獲得更精確的電子能E(MP2),在此基礎上分別借助原子化能反應、原子當量反應、等鍵反應求得化合物的生成熱。同時,將G4(MP2)-6X方法中的幾何優化與振動頻率計算的泛函與基組采用B3LYP/6-311G( d, p)以節約計算資源,并在Bun Chan[19]給出的腳本文件中做出了修改以直接輸出298K下的焓值H(298K),結合原子化能方法計算得到化合物的氣相生成熱。

原子化能方法(Atomization Energies)[15],首先分別計算同一等級下分子與組成原子的焓值,并利用NIST網站[21]查得孤立原子的生成熱,通過計算原子化能反應能來預測分子的氣相生成熱。計算公式如下:

ΔfHm(298K,M)=∑νiΔfHm(Atom,298K)-
[∑νiH(Atom,298K)-H(M,298K)]

(1)

ΔfHm(298K,M)=ΔfHm(0K,M)+
[Hcorr(M)-∑viHcorr(Atom)]

(2)

ΔfHm(0K,M)=∑νiΔfHm(0K,Atom)-ΔrHm(0K)

(3)

ΔrHm(0K)=∑νiE0(Atom)-E0(M)

(4)

式中:H(M,298K)為298K分子的焓值;∑νiH(Atom,298K)為298K的焓值;∑νiΔfHm(Atom,298K)為298K原子的標準生成熱;E0是電子能與零點能之和;Hcorr是分子從0K到298K的熱力學矯正值,孤立原子C、H、O、N的ΔfHm,0K實驗值可由NIST網站[21]查得。

原子當量反應法(Atom Equivalent)[16]與原子化能相似,都是利用量子化學方法計算由分子形成原子時的能量變化,并且結合已知原子的相關信息,從而預測分子的生成熱。原子當量方法計算公式如下:

ΔHi=Ei-∑jnjεj

(5)

式中:Ei是分子i的能量;εj記作一個“原子當量”,定義為εj=(Ej-xj),其中Ej是原子j的能量;nj是分子i中原子j的數目;xj是對原子j的理論水平校正。

等鍵反應法(Isodesmic Reaction)[17]是計算化合物生成熱的一種較為精確的方法。在設計的等鍵反應體系中,化學鍵的類型和數目相同,產物分子和反應物分子電子環境相似,故電子相關造成的誤差可以相互抵消,使得計算生成熱的誤差大大降低。等鍵反應方法基于Hess定律[22],通過設計等鍵反應,計算參與等鍵反應的分子的電子能量(E0)、零點振動能(ZPE)、熱校正值(HT),利用下列等式計算等鍵反應的焓變ΔH298,由NIST網站[21]查的參考物質在標準狀態下的生成熱值,通過式(6)～(10)計算得到目標分子的氣態生成熱。

(6)

ΔE0=∑E0,P-∑E0,R

(7)

ΔZPE=∑ZPEP-∑ZPER

(8)

ΔHT=∑HT,P-∑HT,R

(9)

ΔH298=ΔE0+ΔZPE+ΔHT

(10)

2 結果與討論

本研究以1,2,3-三唑和1,2,4-三唑為基礎,唑環N1引入羥基,由H、—N═N—、—NH—NH—、—NH—、—CH2—、四嗪、吡唑和呋咱連接成雙環三唑,設計16種橋連三唑-N-氧化物分子骨架類型,通過—CH3、—NH2、—OH、—CN、—NO2、—NHNO2、—NHNH2、—C(NO2)2、—C(NO2)3等基團調節能量和穩定性,獲得160種新型含能分子結構式如圖1所示。

圖1 設計的16種橋連三唑分子骨架類型Fig.1 Designed 16 kinds of molecular framework types of pontotriazole

2.1 幾何優化

為了選取合適的計算模型,本研究中160個分子的初始結構都設置為能量更低的反式構象(反式1,2,3-三唑、反式1,2,4-三唑)。幾何優化收斂后,橋基為H(1)、—N═N—(2)、四嗪(6)、吡唑(7)、呋咱(8)的化合物,環骨架仍保持和初始結構相似的平面構型,而其他橋連化合物(橋基為—NH—、—NHNH—、—CH2—),在橋基處有“折角”,使得兩個唑環無法共平面。

由吡唑(7)、呋咱(8)橋接的雙環三唑化合物,均翻轉為順式構型,說明由這兩種橋基連接的化合物更趨于順式構象;而由H(1)、—N═N—(2)、—NH—(3)、—NHNH—(4)、—CH2—(5)、四嗪(6)連接的雙環三唑化合物趨于反式構象。本研究的160個分子的幾何優化結構圖見支撐材料中圖S1～S9,圖2為—H取代基系列的幾何優化構型圖。

圖2 —H取代基系列的結構優化構型圖Fig.2 Optimized configurations of H series

以H(1)、—N═N—(2)、四嗪(6)、吡唑(7)、呋咱(8)為橋基的橋連雙環三唑環,環骨架仍保持和初始結構相似的平面構型,表明當唑環與橋基共軛共平面時,分子的穩定性更高;此外,—N═N—、四嗪、吡唑、呋咱這類橋基,由于體系中引入了N—C、N—N鍵,將有利于提高含能化合物的爆炸性質[23]。而—NH—(3)、—NHNH—(4)、—CH2—(5) 橋連化合物中,兩個唑環面存在“折角”,兩個唑環和橋基之間由于電子排斥作用而呈異面構型。這可能是因為—NH—與—NHNH—的孤對電子在sp3雜化軌道,無法與唑環的游離電子共軛,而—CH2—沒有多余的孤對電子或空軌道;這三種橋基的H原子處在sp3雜化軌道,受空間排斥的影響,兩個唑環與橋基均呈異面構型。

由吡唑(7)、呋咱(8)橋接的雙環三唑化合物均翻轉為了順式構型,說明由這兩種橋基連接的化合物更趨于順式構象;這可能因為反式構型下,吡唑上的H與唑環羥基上的O的距離為1.55605?,相比于順式構型下吡唑H與羥基O的距離為2.47471?,原子間空間斥力更弱;呋咱系列化合物在順式構型下取代基對稱向兩邊伸展,而反式構型下取代基受呋咱的空間排斥。

2.2 不同方法計算的生成熱對比

用G4MP2-6X方法計算了常用炸藥TNT、NTO、FOX-7、TATB與LLM-105的固態生成熱,結果見表1。與NIST網站查閱的實驗值對比,誤差不超過40kJ /mol,證明G4MP2-6X方法預測的標準生成熱結果可信。

表1 C, H, O, N原子的計算焓值和實驗焓值Table 1 Enthalpies for atoms C, H, N and O and their literature values for atomic ΔH°f,298/(kJ/mol)

表1 G4MP2-6X計算方法得到傳統含能化合物的HOF(s, 298K)值及與NIST實驗值對比Table 1 Comparison of the HOF(s, 298K) values calculated by G4MP2-6X and experimental values by NIST

在B3LYP/6-311G(d,p)//MP2/6-311++G(d,p)理論水平下,分別用原子化能法、原子當量法與等鍵反應法對160種設計的含能分子的標準摩爾生成熱進行計算,并與G4MP2-6X方法直接讀取的標準摩爾生成熱值進行對比,結果如表2所示。

表2 C、H、O、N原子當量Table 2 Atom equivalents of C, H, N and O

表2 不同方法計算的HOF(g, 298K)值Table 2 The HOF(g, 298K) values calculated by different methods

表3 G4(MP2)-6X方法結合原子化能方法的計算數據Table 3 Calculation data of G4(MP2)-6X method combing the Atomization energies method

為了對比3種輔助反應(原子化能方法、等鍵反應方法與原子當量方法)求得的標準摩爾生成熱與G4(MP2)-6X方法結合原子化能方法輸出的標準摩爾生成熱的相關性,分別繪制了3種方法與G4(MP2)-6X方法的線性相關圖,并擬合了輔助方法與G4(MP2)-6X方法的生成熱的一元回歸方程,如圖3所示。其中,原子化能方法的相關性最高,相關系數為0.88513;等鍵反應方法次之,相關系數為0.43492;原子當量方法的相關性最低,僅有0.15977。對于本研究的160個含能分子,等鍵反應的相關性較低,可能是因為以—H(1)、—N═N—(2)、四嗪(6)、吡唑(7)、呋咱(8)為連接基團的分子,它們的共軛環骨架在等鍵拆分時被破壞,并且無法利用已有氣相實驗生成熱數據的小分子建立等鍵方程以還原共軛的電子環境,從而帶來誤差。

圖3 不同方法所計算的HOF值與G4(MP2)-6X方法輸出的HOF值的擬合曲線圖Fig.3 Fitting curves of HOF values calculated by different methods with the values by G4(MP2)-6X methods

圖4為按照基團分類的160個化合物的HOF值(原子化能方法計算)。

圖4 原子化能方法計算的160個化合物的HOF(g,298K)值Fig.4 HOF(g,298K) values of 160 kinds of compounds calculated by atomization energy

由圖4可知,取代基—C(NO2)2、—C(NO2)3、—NHNO2和—CN顯著提高了分子的生成熱,而—OH與—NH2取代基系列的生成熱低于—H系列。—C(NO2)2的氮含量與CN鍵數量都不如—C(NO2)3,但HOF值更高,這可能因為偕二硝基基團呈平面構型,共軛的結構使得生成熱提高。同一橋基連接的聯三唑同分異構體(聯1,2,3-三唑,聯1,2,4-三唑),橋連雙1,2,3-三唑的HOF更大,這可能因為1,2,3-三唑的生成焓(272kJ mol-1)比1,2,4-三唑(109kJ mol-1)高。對于10種不同取代基,四嗪聯三唑分子骨架都具有最高的HOF值,可能因為四嗪與三唑環共軛,分子生成熱提高。

3 結 論

(1)探究了目前常用3種輔助方法計算的HOF值與G4(MP2)-6X方法結合原子化能方法輸出HOF值的相關性。首先在B3LYP/6-311G(d, p)理論水平下對160種新型含能分子進行結構優化、頻率計算,并進一步在更高級別的MP2/6-311++G(d, p)理論水平下對分子進行了單點能計算。分別讀取熱力學數值與電子能量,利用3種常用的輔助反應(原子化能方法、等鍵反應方法和原子當量方法)計算HOF值,對比G4(MP2)-6X方法輸出的HOF值,并擬合了與G4(MP2)-6X方法的相關性。

(2)3種方法的相關性為:原子化能方法(0.88513)>等鍵反應方法( 0.43492 )>原子當量方法(0.15977)。建立了線性回歸方程,可采用該方程結合B3LYP/6-311G(d, p)//MP2/6-311++G(d, p)理論水平獲得G4(MP2)-6X方法的計算精度。本研究計算的160個物質在不同理論水平和生成焓計算方法下的生成焓,可供廣大研究人員參考使用。

圖S1 —CH3取代基系列分子的幾何優化構型圖Fig.S1 Optimized configurations of —CH3 series

圖S2 —OH取代基系列分子的幾何優化構型圖Fig.S2 Optimized configurations of —OH series

圖S3 —NH2取代基系列分子的幾何優化構型圖Fig.S3 Optimized configurations of —NH2 series

圖S4 —CN取代基系列分子的幾何優化構型圖Fig.S4 Optimized configurations of —CN series

圖S5 —NO2取代基系列分子的幾何優化構型圖Fig.S5 Optimized configurations of —NO2 series

(相關支撐材料見文后)

支撐材料

1 計算方法

原子化能方法使用的原子焓值如表1所示。

原子當量反應法中使用如表2所示Rice[1]擬合的原子當量。

構建等鍵反應用到的分子包括1,2,4-三唑、1,2,3-三唑、NH3、N2H4、CH3NH2、N2H2、H2、C2H6、C2H4、CH3NHNH2、CH4、CH3OH、C2H5NOHC2H5、C2H5NHC2H5、CH3CH2CH3、C3H4N2、CH3OCH3、CH3CN、CH3NO2、CH3NHCH3、C2H5OC2H5, 在B3LYP/6-311G(d, p)//MP2/6-311++G(d, p)水平進行幾何優化、振動頻率與單點能計算,由NIST網站[2]查得標況下的生成焓值。

圖S6 —NHNO2取代基系列分子的幾何優化構型圖Fig.S6 Optimized configurations of —NHNO2 series

圖S8 —C(NO2)2取代基系列分子的幾何優化構型圖Fig.S8 Optimized configurations of —C(NO2)2 series

圖S9 —C(NO2)3取代基系列分子的幾何優化構型圖Fig.S9 Optimized configurations of —C(NO2)3 series

2 幾何優化構型圖