面向高通量含能分子設(shè)計(jì)篩選的三種生成焓計(jì)算方法評(píng)估

保福成，彭汝芳，張朝陽(yáng)

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

1 引言

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

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

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

2 研究方法

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

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

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

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

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

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

（3）計(jì)算原子化反應(yīng)能（∑D0(M))，即，產(chǎn)物與反應(yīng)物的內(nèi)能差，kJ·mol-1。再通過(guò)式（1）計(jì)算獲得0 K條件下的標(biāo)準(zhǔn)摩爾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）通過(guò)式（2）進(jìn)行溫度從0 K 到298.15 K 的HOF 校正。

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

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

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

式中，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 計(jì)算機(jī)時(shí)對(duì)比

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

圖2 3 種方法的計(jì)算機(jī)時(shí)對(duì)比Fig.2 Comparison of the machine-time by three calculation methods

3.2 生成焓的計(jì)算偏差

表1 列出了20 種含能分子的實(shí)驗(yàn)密度（ρ）、升華焓（Δ(sub)）。研究通過(guò)3 種理論方法獲得298.15 K下的氣相HOF 和固相HOF，結(jié)果也見表1。由表1 可以看出，3 種不同精度的方法計(jì)算結(jié)果差異較大，而B3LYP 方法與CBS 方法計(jì)算結(jié)果較為接近。以高精度的CBS 方法獲得的固相HOF 為標(biāo)準(zhǔn)，基于PM6 和B3LYP 法獲得的固相HOF 的絕對(duì)偏差如圖3a 所示，從圖3a 中可以看出，不同精度的方法計(jì)算獲得的HOF 絕對(duì)差值較大。總體上PM6 計(jì)算結(jié)果的偏差最大，如圖3a中紅色顯示，除分子較大的CL-20 外，B3LYP方法獲得的HOF 偏差較小，如圖3a 中藍(lán)色顯示。3 種方法計(jì)算結(jié)果對(duì)比中，PM6 方法相對(duì)于CBS 方法的平均絕對(duì)偏差為65.1 kJ·mol-1，絕對(duì)偏差最大的分子是DNTAT，為151.5 kJ·mol-1，絕 對(duì) 偏 差 最 小 的 分 子 是TNAZ，為4.2 kJ·mol-1。而B3LYP 方法相對(duì)于CBS 方法的平均絕對(duì)偏差為34.2 kJ·mol-1，絕對(duì)偏差最大的分子是CL-20，為132.6 kJ·mol-1，絕對(duì)偏差最小的分子是RDX，為-1.3 kJ·mol-1。此外，計(jì)算了分子的HOF 在單位質(zhì)量的絕對(duì)偏差和單位體積的絕對(duì)偏差如圖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的計(jì)算結(jié)果中占比較小。

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

表1 20 種CHNO 含能分子的實(shí)驗(yàn)密度、升華焓及3 種方法計(jì)算所得氣相生成焓Δ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 含能分子的實(shí)驗(yàn)密度、升華焓及3 種方法計(jì)算所得氣相生成焓Δ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 生成焓計(jì)算偏差對(duì)爆轟性能的影響評(píng)估

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

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