面向高通量含能分子设计筛选的三种生成焓计算方法评估

保福成，彭汝芳，张朝阳

（1.西南科技大学环境友好能源材料国家重点实验室，四川 绵阳 621010；2.中国工程物理研究院化工材料研究所，四川 绵阳 621999）

1 引言

含能材料是储存大量化学能的亚稳态物质，包括炸药、推进剂和烟火药［1-5］，广泛应用于军事和民用领域。作为军事武器的能量来源，含能材料优异的爆轰性能和较低的感度决定了军事装备的先进性。过去几十年里，诸多致力于追求高性能含能材料分子设计的策略已被提出，如设计具有高张力键的笼形硝基化合物［6-7］和分子中含有大量的N—N 键和N—C 键的高氮化合物［8-10］。此外，为了缓解能量与感度间的矛盾，人们提出了在分子中引入氢键的策略，改善分子的堆积模式以在晶体水平上缓解能量-安全性间矛盾［11-12］。随着计算机技术的快速发展和各种预测模型的建立，高通量分子设计正逐步成为材料设计的主流［13-19］，含能材料也将不会例外［20-23］。

爆轰性能包括爆热、爆速、爆压、爆温和爆容，是含能材料的基本性能，也是含能分子高通量设计要首先考虑的问题之一。生成焓（HOF）是化合物的基本热力学性质与爆轰性能正相关［24-25］，是预测爆轰性能的必要参数。目前，等键反应结合密度泛函理论（DFT）计算的方案［26-27］已经能够给出化学精度的HOF；但是，基于优化的等键反应方案的HOF 计算程序开发难度较高，目前还只能手动进行。原子化方案［28-31］通过计算分子内能和组成原子的气态标准HOF 获得分子的HOF，能够很方便地实现程序化；因此，高通量HOF 计算通常采用原子化方案。然而，原子化方案须同昂贵的高水平量子化学方法相结合才能准确计算分子和原子的内能。常用的量子化学方法有半经验方法，Gong等［32］用PM6 法和PM3 法对56 种高能材料的生成焓进行了预测，预测偏差较大，其均方根偏差分别为42.09 kJ·mol-1和58.83 kJ·mol-1；高精度从头算的完全机组方法CBS，Montgomery 等［33］基于CBS-Q 方法对中小分子进行热化学计算，精度可达4.18～8.36 kJ·mol-1，该方法计算精度较高，但所耗机时较多，不适用于大分子体系。目前密度泛函理论方法得到广泛认可和使用，Rice等［34］在B3LYP/6-31G（d）水平上计算了35 种CHNO分子的HOF，预测气相分子HOF 与实验值的均方根偏差 为12.97 kJ·mol-1，最 大 偏 差 为25.52 kJ·mol-1。然而，在高通量分子设计筛选中，高精度方法会急剧增加计算成本，而半经验方法又难以满足计算精度需要。

因此，为给高通量筛选的理论方法选择提供依据，一方面须加快基于等键反应方案的HOF 计算的程序化进程；另一方面，须评估不同理论水平基于原子化方案计算HOF 的偏差对爆轰性能的影响程度。基于此，本研究对3 种不同水平理论方法，半经验PM6 方法、密度泛函理论方法B3LYP/6-31G（d，p）和高精度的完全基组CBS-4M 方法，在面向高通量含能分子设计筛选的爆轰性能预测中的适宜性进行了评估，以满足高通量筛选所要求的程序化、预测精度与计算效率。

2 研究方法

本研究旨在获得高通量的计算方法，而不是要进行高通量的分子设计。研究的对象为20 种CHNO 类含能分子，涵盖了硝基化合物炸药、硝胺炸药、硝酸酯炸药、呋咱类、嗪类及叠氮类炸药等，它们的分子结构如图1 所示。选用的CHNO 含能分子，一方面有实验HOF 值，便于比较，能确定计算方法的准确性；另一方面，结构上具有代表性和多样性，如链状、环状与笼形的分子结构，同时具有C─NO2、N─NO2和O─NO2官能团，待进行高通量设计的分子大多具有这一类似的结构。

图1 20 种含能分子的结构Fig.1 Structures of twenty energetic molecules

由于计算方法对CHNO 体系的可行性没有结构上的限制，因此3 种不同水平的方法可计算得到以上20 种含能分子标准状态下的气相 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 程序包）计算时，结合了原子化方案。基于原子化方案计算Δf(M，g，298.15 K)方法简要介绍如下［28-31］：

（1）参考了文献［31］构造M 原子化反应。

式中，x1，x2，x3，x4表示原子个数。

（2）采用上述量子化学方法对分子M 进行结构优化，频率计算无虚频确认势能面上稳定的结构。

（3）计算原子化反应能（∑D0(M))，即，产物与反应物的内能差，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 参数为固相标准摩尔HOF（Δf(M，s，298.15 K）），可用标准状态下的气相HOF（Δf(M，g，298.15 K））减去升华焓（Δ(sub)）得到。对于Δ(sub)，则按照经验式（3）计算［35-36］。

式中，A、ν和分别表示分子表面积、分子表面正负静电势平衡常数与总表面静电势的方差，可由Bulat 等［37］提出的计算方法计算获得。常数α=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经验方程式见式（4）～（5）［38-40］、BKW 模型［41-44］和VLW 模型［45］对10种典型CHNO 含能分子的爆轰性能进行计算。

式中，N表示每克炸药爆炸产生的气体摩尔数，mol·g-1；Mave为气体爆轰产物的平均摩尔质量，g·mol-1；ρ为炸药的装载密度，g·cm-3；Q为爆热，kJ·g-1。

3 结果和讨论

3.1 计算机时对比

研 究 所 用CPU 是2.6 GHz 的Intel，Gold 6142。不同精度的量子化学方法计算HOF 所耗机时不同，3种方法的不同机时如图2 所示，由图2 可见，PM6 半经验方法所需时间最短，且无明显分子大小依赖性，计算平均机时仅约5.0×10-4核时。DFT（B3LYP）方法所耗机时较PM6 方法高约3 个数量级，平均机时0.36 核时。而CBS 方法耗时最大，20 种分子计算的平均机时为每分子1.06 核时，特别是对大分子体系计算所耗机时巨大，如大分子CL-20 所需的机时较其他小分子体系高约1 个数量级。仅从计算机时考虑，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，结果也见表1。由表1 可以看出，3 种不同精度的方法计算结果差异较大，而B3LYP 方法与CBS 方法计算结果较为接近。以高精度的CBS 方法获得的固相HOF 为标准，基于PM6 和B3LYP 法获得的固相HOF 的绝对偏差如图3a 所示，从图3a 中可以看出，不同精度的方法计算获得的HOF 绝对差值较大。总体上PM6 计算结果的偏差最大，如图3a中红色显示，除分子较大的CL-20 外，B3LYP方法获得的HOF 偏差较小，如图3a 中蓝色显示。3 种方法计算结果对比中，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 在单位质量的绝对偏差和单位体积的绝对偏差如图3b～3c所示，基于PM6 方法和B3LYP 方法获得平均单位质量的HOF 偏差分别为0.30 kJ·g-1和0.15 kJ·g-1，平均单位体积的HOF偏差分别为0.37 kJ·Å-3和0.19 kJ·Å-3，尽管这些含能分子基于不同方法获得的HOF 具有较高的偏差值，但其偏差在后续Q的计算结果中占比较小。

图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 模型预测得到10 种含能分子的D和P如图5aⅠ～5cⅡ所示。根据K-J经验方程知D和p与Q密切相关，因此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 的计算偏差对含能材料的爆轰性能预测有一定的影响。

基于PM6 和B3LYP 法获得的HOF 预测得到的爆轰性能与基于CBS 结果的RE 如图5aⅢ～5cⅣ所示，尽管不同精度的计算方法获得的HOF 绝对偏差较大，但预测得到爆轰性能的RE 较小，如基于PM6 方法获得的HOF 和K-J 方 程、BKW 模 型、VLW 模 型 预 测D的 平均RE 分别为1.6%、1.0%和1.5%，预测p的平均RE 分别为3.2%、2.9% 和5.3%。基于B3LYP 方法获得的HOF 和K-J 方 程、BKW 模 型、VLW 模 型 预 测D的 平 均RE 分别为0.6%、0.4%和0.6%，预测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 种模型预测的爆轰性能结果可以看出，基于PM6 方法计算得到的HOF 预测的爆轰性能偏差最大，在高通量分子设计筛选中难以满足计算精度要求，而基于B3LYP 方法计算得到的HOF 预测的爆轰性能的偏差在可接受的范围内。此外，通过比较HOF 偏差和爆轰性能偏差发现，尽管基于3 种精度获得10 种含能分子的HOF 偏差较大，但该偏差对Q、D及p的影响较小。究其原因，一方面CHNO 类含能分子的Q除了与分子本身HOF 大小相关外，还与爆炸产物的种类和数量相关，产物HOF 大小对Q的影响不可忽略；另一方面，根据K-J 经验方程，密度对D和p的影响远大于Q对D和p的影响，即D～Q0.25，p～Q0.5，D～ρ，p～ρ2。因此，在高通量含能分子设计筛选中，选择中等水平的B3LYP 方法即可满足HOF 计算的精度需要，而选择合适的爆炸反应气体规则和准确预测密度的方法需要我们重点关注。

4 结论

基于3 种量子化学计算方法和原子化方案计算获得了20 种含能分子的气相HOF，结合升华焓获得了其固相HOF，基于实验密度、固相HOF 和3 种方法预测了10 种常见含能分子爆轰性能，讨论了HOF 的计算偏差对爆轰性能的影响程度。评估了面向高通量含能分子筛选的3 种方法的适宜性，结论如下：

（1）用于HOF 计算的3 种不同精度的方法所耗机时差别较大，PM6 方法能够快速获得所需的结果，B3LYP 方法的计算机时也在可接受范围内，而CBS 方法耗时巨大，不满足高通量分子设计筛选中高效的要求。

（2）3 种不同精度计算的HOF 差别较大，以高精度的CBS 方法为标准，PM6 方法计算结果的平均绝对偏差为65.1 kJ·mol-1，而B3LYP 方法计算结果的平均绝对偏差为34.2 kJ·mol-1。

（3）研究表明传统CHNO 类含能材料的爆轰产物对爆热Q的贡献较大，故其HOF 的计算偏差对Q的影响较小，其中基于PM6 方法获得的HOF 计算Q的平均RE 为6.5%，而基于B3LYP 方法获得的HOF 计算Q的平均RE 仅为2.8%。

（4）采用K-J、BKW 和VLW 模型预测了10 种含能分子的D和p，发现在较大HOF 的绝对偏差下，基于PM6 方法获得的HOF 计算D和p的平均RE 分别小于1.6%和5.3%，基于B3LYP 方法获得的HOF 计算D和p的平均RE 分别小于0.6%和1.9%。

综上所述，在高通量含能分子筛选中，HOF 的预测采用中等精度的B3LYP 方法即可满足高通量分子设计的筛选需要，已将此方法置于了含能材料高通量计算平台EM Studio 1.0，作为默认的HOF 计算方法［68］。当然，随着完善的基于等键反应方案的HOF 程序的出现，准确与效率兼顾的HOF 计算方法将更有利于高通量含能分子筛选。