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

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

(中南大學(xué) 化學(xué)化工學(xué)院,湖南 長(zhǎng)沙 410083)

引 言

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

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

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

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

1 計(jì)算方法

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

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

Δ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原子的標(biāo)準(zhǔn)生成熱;E0是電子能與零點(diǎn)能之和;Hcorr是分子從0K到298K的熱力學(xué)矯正值,孤立原子C、H、O、N的ΔfHm,0K實(shí)驗(yàn)值可由NIST網(wǎng)站[21]查得。

原子當(dāng)量反應(yīng)法(Atom Equivalent)[16]與原子化能相似,都是利用量子化學(xué)方法計(jì)算由分子形成原子時(shí)的能量變化,并且結(jié)合已知原子的相關(guān)信息,從而預(yù)測(cè)分子的生成熱。原子當(dāng)量方法計(jì)算公式如下:

ΔHi=Ei-∑jnjεj

(5)

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

等鍵反應(yīng)法(Isodesmic Reaction)[17]是計(jì)算化合物生成熱的一種較為精確的方法。在設(shè)計(jì)的等鍵反應(yīng)體系中,化學(xué)鍵的類型和數(shù)目相同,產(chǎn)物分子和反應(yīng)物分子電子環(huán)境相似,故電子相關(guān)造成的誤差可以相互抵消,使得計(jì)算生成熱的誤差大大降低。等鍵反應(yīng)方法基于Hess定律[22],通過設(shè)計(jì)等鍵反應(yīng),計(jì)算參與等鍵反應(yīng)的分子的電子能量(E0)、零點(diǎn)振動(dòng)能(ZPE)、熱校正值(HT),利用下列等式計(jì)算等鍵反應(yīng)的焓變?chǔ)298,由NIST網(wǎng)站[21]查的參考物質(zhì)在標(biāo)準(zhǔn)狀態(tài)下的生成熱值,通過式(6)～(10)計(jì)算得到目標(biāo)分子的氣態(tài)生成熱。

(6)

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

(7)

ΔZPE=∑ZPEP-∑ZPER

(8)

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

(9)

ΔH298=ΔE0+ΔZPE+ΔHT

(10)

2 結(jié)果與討論

本研究以1,2,3-三唑和1,2,4-三唑?yàn)榛A(chǔ),唑環(huán)N1引入羥基,由H、—N═N—、—NH—NH—、—NH—、—CH2—、四嗪、吡唑和呋咱連接成雙環(huán)三唑,設(shè)計(jì)16種橋連三唑-N-氧化物分子骨架類型,通過—CH3、—NH2、—OH、—CN、—NO2、—NHNO2、—NHNH2、—C(NO2)2、—C(NO2)3等基團(tuán)調(diào)節(jié)能量和穩(wěn)定性,獲得160種新型含能分子結(jié)構(gòu)式如圖1所示。

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

2.1 幾何優(yōu)化

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

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

圖2 —H取代基系列的結(jié)構(gòu)優(yōu)化構(gòu)型圖Fig.2 Optimized configurations of H series

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

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

2.2 不同方法計(jì)算的生成熱對(duì)比

用G4MP2-6X方法計(jì)算了常用炸藥TNT、NTO、FOX-7、TATB與LLM-105的固態(tài)生成熱,結(jié)果見表1。與NIST網(wǎng)站查閱的實(shí)驗(yàn)值對(duì)比,誤差不超過40kJ /mol,證明G4MP2-6X方法預(yù)測(cè)的標(biāo)準(zhǔn)生成熱結(jié)果可信。

表1 C, H, O, N原子的計(jì)算焓值和實(shí)驗(yà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計(jì)算方法得到傳統(tǒng)含能化合物的HOF(s, 298K)值及與NIST實(shí)驗(yàn)值對(duì)比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)理論水平下,分別用原子化能法、原子當(dāng)量法與等鍵反應(yīng)法對(duì)160種設(shè)計(jì)的含能分子的標(biāo)準(zhǔn)摩爾生成熱進(jìn)行計(jì)算,并與G4MP2-6X方法直接讀取的標(biāo)準(zhǔn)摩爾生成熱值進(jìn)行對(duì)比,結(jié)果如表2所示。

表2 C、H、O、N原子當(dāng)量Table 2 Atom equivalents of C, H, N and O

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

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

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

圖3 不同方法所計(jì)算的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為按照基團(tuán)分類的160個(gè)化合物的HOF值(原子化能方法計(jì)算)。

圖4 原子化能方法計(jì)算的160個(gè)化合物的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系列?！狢(NO2)2的氮含量與CN鍵數(shù)量都不如—C(NO2)3,但HOF值更高,這可能因?yàn)橘啥趸鶊F(tuán)呈平面構(gòu)型,共軛的結(jié)構(gòu)使得生成熱提高。同一橋基連接的聯(lián)三唑同分異構(gòu)體(聯(lián)1,2,3-三唑,聯(lián)1,2,4-三唑),橋連雙1,2,3-三唑的HOF更大,這可能因?yàn)?,2,3-三唑的生成焓(272kJ mol-1)比1,2,4-三唑(109kJ mol-1)高。對(duì)于10種不同取代基,四嗪聯(lián)三唑分子骨架都具有最高的HOF值,可能因?yàn)樗泥号c三唑環(huán)共軛,分子生成熱提高。

3 結(jié) 論

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

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

圖S1 —CH3取代基系列分子的幾何優(yōu)化構(gòu)型圖Fig.S1 Optimized configurations of —CH3 series

圖S2 —OH取代基系列分子的幾何優(yōu)化構(gòu)型圖Fig.S2 Optimized configurations of —OH series

圖S3 —NH2取代基系列分子的幾何優(yōu)化構(gòu)型圖Fig.S3 Optimized configurations of —NH2 series

圖S4 —CN取代基系列分子的幾何優(yōu)化構(gòu)型圖Fig.S4 Optimized configurations of —CN series

圖S5 —NO2取代基系列分子的幾何優(yōu)化構(gòu)型圖Fig.S5 Optimized configurations of —NO2 series

(相關(guān)支撐材料見文后)

支撐材料

1 計(jì)算方法

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

原子當(dāng)量反應(yīng)法中使用如表2所示Rice[1]擬合的原子當(dāng)量。

構(gòu)建等鍵反應(yīng)用到的分子包括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)水平進(jìn)行幾何優(yōu)化、振動(dòng)頻率與單點(diǎn)能計(jì)算,由NIST網(wǎng)站[2]查得標(biāo)況下的生成焓值。

圖S6 —NHNO2取代基系列分子的幾何優(yōu)化構(gòu)型圖Fig.S6 Optimized configurations of —NHNO2 series

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

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

2 幾何優(yōu)化構(gòu)型圖