保福成，彭汝芳，張朝陽

（1.西南科技大學(xué)環(huán)境友好能源材料國家重點實驗室，四川 綿陽 621010；2.中國工程物理研究院化工材料研究所，四川 綿陽 621999）

1 引言

含能材料是儲存大量化學(xué)能的亞穩(wěn)態(tài)物質(zhì)，包括炸藥、推進劑和煙火藥［1-5］，廣泛應(yīng)用于軍事和民用領(lǐng)域。作為軍事武器的能量來源，含能材料優(yōu)異的爆轟性能和較低的感度決定了軍事裝備的先進性。過去幾十年里，諸多致力于追求高性能含能材料分子設(shè)計的策略已被提出，如設(shè)計具有高張力鍵的籠形硝基化合物［6-7］和分子中含有大量的N—N 鍵和N—C 鍵的高氮化合物［8-10］。此外，為了緩解能量與感度間的矛盾，人們提出了在分子中引入氫鍵的策略，改善分子的堆積模式以在晶體水平上緩解能量-安全性間矛盾［11-12］。隨著計算機技術(shù)的快速發(fā)展和各種預(yù)測模型的建立，高通量分子設(shè)計正逐步成為材料設(shè)計的主流［13-19］，含能材料也將不會例外［20-23］。

爆轟性能包括爆熱、爆速、爆壓、爆溫和爆容，是含能材料的基本性能，也是含能分子高通量設(shè)計要首先考慮的問題之一。生成焓（HOF）是化合物的基本熱力學(xué)性質(zhì)與爆轟性能正相關(guān)［24-25］，是預(yù)測爆轟性能的必要參數(shù)。目前，等鍵反應(yīng)結(jié)合密度泛函理論（DFT）計算的方案［26-27］已經(jīng)能夠給出化學(xué)精度的HOF；但是，基于優(yōu)化的等鍵反應(yīng)方案的HOF 計算程序開發(fā)難度較高，目前還只能手動進行。原子化方案［28-31］通過計算分子內(nèi)能和組成原子的氣態(tài)標準HOF 獲得分子的HOF，能夠很方便地實現(xiàn)程序化；因此，高通量HOF 計算通常采用原子化方案。然而，原子化方案須同昂貴的高水平量子化學(xué)方法相結(jié)合才能準確計算分子和原子的內(nèi)能。常用的量子化學(xué)方法有半經(jīng)驗方法，Gong等［32］用PM6 法和PM3 法對56 種高能材料的生成焓進行了預(yù)測，預(yù)測偏差較大，其均方根偏差分別為42.09 kJ·mol-1和58.83 kJ·mol-1；高精度從頭算的完全機組方法CBS，Montgomery 等［33］基于CBS-Q 方法對中小分子進行熱化學(xué)計算，精度可達4.18～8.36 kJ·mol-1，該方法計算精度較高，但所耗機時較多，不適用于大分子體系。目前密度泛函理論方法得到廣泛認可和使用，Rice等［34］在B3LYP/6-31G（d）水平上計算了35 種CHNO分子的HOF，預(yù)測氣相分子HOF 與實驗值的均方根偏差 為12.97 kJ·mol-1，最 大 偏 差 為25.52 kJ·mol-1。然而，在高通量分子設(shè)計篩選中，高精度方法會急劇增加計算成本，而半經(jīng)驗方法又難以滿足計算精度需要。

因此，為給高通量篩選的理論方法選擇提供依據(jù)，一方面須加快基于等鍵反應(yīng)方案的HOF 計算的程序化進程；另一方面，須評估不同理論水平基于原子化方案計算HOF 的偏差對爆轟性能的影響程度?；诖?，本研究對3 種不同水平理論方法，半經(jīng)驗PM6 方法、密度泛函理論方法B3LYP/6-31G（d，p）和高精度的完全基組CBS-4M 方法，在面向高通量含能分子設(shè)計篩選的爆轟性能預(yù)測中的適宜性進行了評估，以滿足高通量篩選所要求的程序化、預(yù)測精度與計算效率。

2 研究方法

本研究旨在獲得高通量的計算方法，而不是要進行高通量的分子設(shè)計。研究的對象為20 種CHNO 類含能分子，涵蓋了硝基化合物炸藥、硝胺炸藥、硝酸酯炸藥、呋咱類、嗪類及疊氮類炸藥等，它們的分子結(jié)構(gòu)如圖1 所示。選用的CHNO 含能分子，一方面有實驗HOF 值，便于比較，能確定計算方法的準確性；另一方面，結(jié)構(gòu)上具有代表性和多樣性，如鏈狀、環(huán)狀與籠形的分子結(jié)構(gòu)，同時具有C─NO2、N─NO2和O─NO2官能團，待進行高通量設(shè)計的分子大多具有這一類似的結(jié)構(gòu)。

圖1 20 種含能分子的結(jié)構(gòu)Fig.1 Structures of twenty energetic molecules

由于計算方法對CHNO 體系的可行性沒有結(jié)構(gòu)上的限制，因此3 種不同水平的方法可計算得到以上20 種含能分子標準狀態(tài)下的氣相 HOF（Δf(M，g，298.15 K)）。對于Δf(M，g，298.15 K)，采用PM6 方法（VAMP 程序包）計算便可直接獲得；而采 用 DFT-B3LYP/6-31G（d，p）和 CBS-4M 方 法（Gaussian 03 程序包）計算時，結(jié)合了原子化方案?；谠踊桨赣嬎悝(M，g，298.15 K)方法簡要介紹如下［28-31］：

（1）參考了文獻［31］構(gòu)造M 原子化反應(yīng)。

式中，x1，x2，x3，x4表示原子個數(shù)。

（2）采用上述量子化學(xué)方法對分子M 進行結(jié)構(gòu)優(yōu)化，頻率計算無虛頻確認勢能面上穩(wěn)定的結(jié)構(gòu)。

（3）計算原子化反應(yīng)能（∑D0(M))，即，產(chǎn)物與反應(yīng)物的內(nèi)能差，kJ·mol-1。再通過式（1）計算獲得0 K條件下的標準摩爾HOF（Δf（M，g，0 K）。

式 中，Δf(X，g，0 K)表 示 原 子X 在0 K 下 的HOF，kJ·mol-1；C、H、N 和O 原子的Δf(X，g，0 K)分別為711.19、216.02、470.82 kJ·mol-1和246.81 kJ·mol-1。

（4）通過式（2）進行溫度從0 K 到298.15 K 的HOF 校正。

含能分子M 爆轟性能計算需要的HOF 參數(shù)為固相標準摩爾HOF（Δf(M，s，298.15 K）），可用標準狀態(tài)下的氣相HOF（Δf(M，g，298.15 K））減去升華焓（Δ(sub)）得到。對于Δ(sub)，則按照經(jīng)驗式（3）計算［35-36］。

式中，A、ν和分別表示分子表面積、分子表面正負靜電勢平衡常數(shù)與總表面靜電勢的方差，可由Bulat 等［37］提出的計算方法計算獲得。常數(shù)α=11.171×10-4kJ·mol-1，β= 6.904 kJ·mol-1，γ=12.409 kJ·mol-1。計算Δf（M，g，298.15 K）與Δ（sub）的 差 值，便 可 得 到Δf(M，s，298.15 K)。

基于上述HOF 和實驗密度，分別采用K-J經(jīng)驗方程式見式（4）～（5）［38-40］、BKW 模型［41-44］和VLW 模型［45］對10種典型CHNO 含能分子的爆轟性能進行計算。

式中，N表示每克炸藥爆炸產(chǎn)生的氣體摩爾數(shù)，mol·g-1；Mave為氣體爆轟產(chǎn)物的平均摩爾質(zhì)量，g·mol-1；ρ為炸藥的裝載密度，g·cm-3；Q為爆熱，kJ·g-1。

3 結(jié)果和討論

3.1 計算機時對比

研 究 所 用CPU 是2.6 GHz 的Intel，Gold 6142。不同精度的量子化學(xué)方法計算HOF 所耗機時不同，3種方法的不同機時如圖2 所示，由圖2 可見，PM6 半經(jīng)驗方法所需時間最短，且無明顯分子大小依賴性，計算平均機時僅約5.0×10-4核時。DFT（B3LYP）方法所耗機時較PM6 方法高約3 個數(shù)量級，平均機時0.36 核時。而CBS 方法耗時最大，20 種分子計算的平均機時為每分子1.06 核時，特別是對大分子體系計算所耗機時巨大，如大分子CL-20 所需的機時較其他小分子體系高約1 個數(shù)量級。僅從計算機時考慮，PM6 方法能夠在最低的成本下進行高效的高通量篩選，而CBS 方法耗時較大，計算成本高，不利于高通量篩選。

圖2 3 種方法的計算機時對比Fig.2 Comparison of the machine-time by three calculation methods

3.2 生成焓的計算偏差

表1 列出了20 種含能分子的實驗密度（ρ）、升華焓（Δ(sub)）。研究通過3 種理論方法獲得298.15 K下的氣相HOF 和固相HOF，結(jié)果也見表1。由表1 可以看出，3 種不同精度的方法計算結(jié)果差異較大，而B3LYP 方法與CBS 方法計算結(jié)果較為接近。以高精度的CBS 方法獲得的固相HOF 為標準，基于PM6 和B3LYP 法獲得的固相HOF 的絕對偏差如圖3a 所示，從圖3a 中可以看出，不同精度的方法計算獲得的HOF 絕對差值較大?？傮w上PM6 計算結(jié)果的偏差最大，如圖3a中紅色顯示，除分子較大的CL-20 外，B3LYP方法獲得的HOF 偏差較小，如圖3a 中藍色顯示。3 種方法計算結(jié)果對比中，PM6 方法相對于CBS 方法的平均絕對偏差為65.1 kJ·mol-1，絕對偏差最大的分子是DNTAT，為151.5 kJ·mol-1，絕 對 偏 差 最 小 的 分 子 是TNAZ，為4.2 kJ·mol-1。而B3LYP 方法相對于CBS 方法的平均絕對偏差為34.2 kJ·mol-1，絕對偏差最大的分子是CL-20，為132.6 kJ·mol-1，絕對偏差最小的分子是RDX，為-1.3 kJ·mol-1。此外，計算了分子的HOF 在單位質(zhì)量的絕對偏差和單位體積的絕對偏差如圖3b～3c所示，基于PM6 方法和B3LYP 方法獲得平均單位質(zhì)量的HOF 偏差分別為0.30 kJ·g-1和0.15 kJ·g-1，平均單位體積的HOF偏差分別為0.37 kJ·?-3和0.19 kJ·?-3，盡管這些含能分子基于不同方法獲得的HOF 具有較高的偏差值，但其偏差在后續(xù)Q的計算結(jié)果中占比較小。

圖3 相對于CBS 方法的PM6 與B3LYP 的HOF 的計算偏差Fig.3 Error of the HOF calculated by PM6 and B3LYP compared with the CBS method

表1 20 種CHNO 含能分子的實驗密度、升華焓及3 種方法計算所得氣相生成焓Δf(g)和固相生成焓Δf(s)Table 1 The experiment density，calculated enthalpy of sublimation，gas-phase heat of formation Δf(g)and solid-state heat of formation Δf(s) calculated based on three methods for twenty energetic molecules

表1 20 種CHNO 含能分子的實驗密度、升華焓及3 種方法計算所得氣相生成焓Δf(g)和固相生成焓Δf(s)Table 1 The experiment density，calculated enthalpy of sublimation，gas-phase heat of formation Δf(g)and solid-state heat of formation Δf(s) calculated based on three methods for twenty energetic molecules

Note：the ρ，enthalpy of sublimation（Δ(sub)），heat of formation in gas state（Δf(g)）and solid state（Δf(s)）are in g·cm-3，kJ·mol-1，kJ·mol-1 and kJ·mol-1，respectively.

compounds ρ ΔH ? m(sub)ΔfH ?m(g)PM6 341.8 142.3 65.7-524.3-33.1 206.3 38.1-82.0 160.2 59.4 587.9 995.8 677.4 319.7 43.5 1078.6 362.8 95.8 572.0 809.2 B3LYP 585.8 233.9 172.8-402.9-15.1 95.0-12.1-12.6 174.5 61.5 507.5 882.0 697.1 309.2 32.2 1017.1 381.2 120.5 635.5 696.2 CBS 453.1 216.7 174.1-406.3 12.1 120.9-17.6 9.2 106.3 9.2 543.9 844.3 760.2 312.1 49.8 1094.5 397.9 91.6 652.3 727.6 B3LYP 411.3 33.5 26.4-555.2-193.3-70.7-151-200.9 33.5-68.2 360.9 749.4 608.0 182.1-111.2 866.7 237.3 6.8 424.2 583.5 ΔfH ?m(s)PM6 167.3-58.1-80.7-676.6-211.3 40.6-100.8-270.3 19.2-70.3 441.3 863.2 588.3 192.6-99.9 928.2 218.9-17.9 360.7 696.5 CL-20 HMX RDX PETN FOX-7 LLM-105 NTO TATB Tetryl TNT DHDFP DNTAT BTF NTAN NQ DADYT DNFP TNAZ DHT TAFP CBS 278.6 16.3 27.7-558.6-166.1-44.8-156.5-179.1-34.7-120.5 397.3 711.7 671.1 185.0-93.6 944.1 254.0-22.1 441.0 614.9 2.044［6］1.894［46］1.806［47］1.781［48］1.893［49］1.919［50］1.918［51］1.937［52］1.731［53］1.654［54］2.008［55］1.901［56］1.932［57］1.919［58］1.752［59］1.774［60］1.829［61］1.861［62］1.729［63］1.834［64］174.5 200.4 146.4 152.3 178.2 165.7 138.9 188.3 141.0 129.7 146.6 132.6 89.1 127.1 143.4 150.4 143.9 113.7 211.3 112.7

3.2 生成焓計算偏差對爆轟性能的影響評估

圖4 基于3 種方法計算獲得的10 種含能分子的爆熱Q 及其相對偏差REFig.4 Explosive heat Q and relative error RE of ten energetic molecules calculated based on three methods

基于3 種方法獲得的固相HOF、實驗密度和K-J、BKW 與VLW 模型預(yù)測得到10 種含能分子的D和P如圖5aⅠ～5cⅡ所示。根據(jù)K-J經(jīng)驗方程知D和p與Q密切相關(guān)，因此HOF的大小在一定程度上會影響D、p的大小。對于相同的含能分子，基于不同理論方法獲得的HOF 計算的爆轟性能不同，例如CL-20，B3LYP 法獲得的HOF 最大（585.8 kJ·mol-1），基于BKW 模型計算得到最大的D（9.8 km·s-1）和p（45200 MPa），而PM6 法獲得 的HOF 最低（341.8 kJ·mol-1），其D和p相對較小，分別為9.6 km·s-1和43500 MPa 如圖5bⅠ和5bⅡ所示。因此，HOF 的計算偏差對含能材料的爆轟性能預(yù)測有一定的影響。

基于PM6 和B3LYP 法獲得的HOF 預(yù)測得到的爆轟性能與基于CBS 結(jié)果的RE 如圖5aⅢ～5cⅣ所示，盡管不同精度的計算方法獲得的HOF 絕對偏差較大，但預(yù)測得到爆轟性能的RE 較小，如基于PM6 方法獲得的HOF 和K-J 方 程、BKW 模 型、VLW 模 型 預(yù) 測D的 平均RE 分別為1.6%、1.0%和1.5%，預(yù)測p的平均RE 分別為3.2%、2.9% 和5.3%?；贐3LYP 方法獲得的HOF 和K-J 方 程、BKW 模 型、VLW 模 型 預(yù) 測D的 平 均RE 分別為0.6%、0.4%和0.6%，預(yù)測p的平均RE 分別為1.2%、1.2%和1.9%。

圖5 基于3 種不同模型計算獲得10 種含能分子的爆速、爆壓及相對偏差Fig.5 Detonation velocity，detonation pressure，and relative error of ten energetic molecules calculated by three different models

對比基于3 種模型預(yù)測的爆轟性能結(jié)果可以看出，基于PM6 方法計算得到的HOF 預(yù)測的爆轟性能偏差最大，在高通量分子設(shè)計篩選中難以滿足計算精度要求，而基于B3LYP 方法計算得到的HOF 預(yù)測的爆轟性能的偏差在可接受的范圍內(nèi)。此外，通過比較HOF 偏差和爆轟性能偏差發(fā)現(xiàn)，盡管基于3 種精度獲得10 種含能分子的HOF 偏差較大，但該偏差對Q、D及p的影響較小。究其原因，一方面CHNO 類含能分子的Q除了與分子本身HOF 大小相關(guān)外，還與爆炸產(chǎn)物的種類和數(shù)量相關(guān)，產(chǎn)物HOF 大小對Q的影響不可忽略；另一方面，根據(jù)K-J 經(jīng)驗方程，密度對D和p的影響遠大于Q對D和p的影響，即D～Q0.25，p～Q0.5，D～ρ，p～ρ2。因此，在高通量含能分子設(shè)計篩選中，選擇中等水平的B3LYP 方法即可滿足HOF 計算的精度需要，而選擇合適的爆炸反應(yīng)氣體規(guī)則和準確預(yù)測密度的方法需要我們重點關(guān)注。

4 結(jié)論

基于3 種量子化學(xué)計算方法和原子化方案計算獲得了20 種含能分子的氣相HOF，結(jié)合升華焓獲得了其固相HOF，基于實驗密度、固相HOF 和3 種方法預(yù)測了10 種常見含能分子爆轟性能，討論了HOF 的計算偏差對爆轟性能的影響程度。評估了面向高通量含能分子篩選的3 種方法的適宜性，結(jié)論如下：

（1）用于HOF 計算的3 種不同精度的方法所耗機時差別較大，PM6 方法能夠快速獲得所需的結(jié)果，B3LYP 方法的計算機時也在可接受范圍內(nèi)，而CBS 方法耗時巨大，不滿足高通量分子設(shè)計篩選中高效的要求。

（2）3 種不同精度計算的HOF 差別較大，以高精度的CBS 方法為標準，PM6 方法計算結(jié)果的平均絕對偏差為65.1 kJ·mol-1，而B3LYP 方法計算結(jié)果的平均絕對偏差為34.2 kJ·mol-1。

（3）研究表明傳統(tǒng)CHNO 類含能材料的爆轟產(chǎn)物對爆熱Q的貢獻較大，故其HOF 的計算偏差對Q的影響較小，其中基于PM6 方法獲得的HOF 計算Q的平均RE 為6.5%，而基于B3LYP 方法獲得的HOF 計算Q的平均RE 僅為2.8%。

（4）采用K-J、BKW 和VLW 模型預(yù)測了10 種含能分子的D和p，發(fā)現(xiàn)在較大HOF 的絕對偏差下，基于PM6 方法獲得的HOF 計算D和p的平均RE 分別小于1.6%和5.3%，基于B3LYP 方法獲得的HOF 計算D和p的平均RE 分別小于0.6%和1.9%。

綜上所述，在高通量含能分子篩選中，HOF 的預(yù)測采用中等精度的B3LYP 方法即可滿足高通量分子設(shè)計的篩選需要，已將此方法置于了含能材料高通量計算平臺EM Studio 1.0，作為默認的HOF 計算方法［68］。當(dāng)然，隨著完善的基于等鍵反應(yīng)方案的HOF 程序的出現(xiàn)，準確與效率兼顧的HOF 計算方法將更有利于高通量含能分子篩選。