<title></title>
此楼可删
83156
%7B%22isLastPage%22%3Atrue%2C%22notes%22%3A%5B%5D%2C%22pid%22%3A%22850049%22%2C%22tid%22%3A%2283156%22%2C%22mainForumsId%22%3A%5B%2283%22%5D%2C%22categoriesId%22%3A%5B%22426%22%5D%2C%22tcId%22%3A%5B%5D%7D
%7B%22isEditMode%22%3Afalse%7D
量子化学软件仿真（四）：计算中性分子的生成焓
<title></title>
前三篇文章分论述了含能材料生成焓与密度的仿真计算，其中分子密度仿真介绍的较为透彻，可以达到主流学术论文中的精度水平。但是生成焓采用了PM7半经验法，误差较大，检索文献结果显示应用也较少。所以补齐较高精度的生成焓仿真方法成为了含能材料设计工作的最后一环。
本文仍然沿用前文中的ChemBioOffice 2014、Multiwfn软件，但更新了Gaussian到G09W A.02版本。
一、仿真方法
查阅大量文献可知目前的仿真方法有：1、半经验法AM1、PM3、MNDO等。2、完全基组法CBS-QB3、CBS-4M等。3、密度泛函B3LYP/6-31G(d)。4、Hartree-Fock法。5、Gaussian1~4系列。6、W1~4系列。
表1 仿真方法对比表
| <title></title>方法 | 耗时 | 精度 | 适用分子体系 |
|---|---|---|---|
| 半经验法<title></title> | 少 | 低 | 大 |
| 完全基组法 | 中等 | 良好 | 中 |
| 密度泛函 | 中等 | 良好 | 大 |
| Gaussian1~4系列<title></title> | 较长 | 高 | 小 |
| W1~4系列 | 极长 | 极高 | 极 小 |
<title></title>
和PM3类似的PM7在之前有过介绍，该方法耗时极短即使计算ONC这样复杂的分子时也仅需运行几秒，但只能得到气相生成焓且精度不高。该类半经验法仅能提供有限的参考，或在对比大批量的分子时提供一个宏观的指导方向。
德国的T.M.Klapotke[1]研究组偏爱使用完全基组法中的CBS-4M，在他署名的论文中几乎全部使用了CBS-4M来得到材料的生成焓。该方法在计算耗时与准确性方面取得了一个较好的平衡值，以NG为算例，耗时40 min，误差18.7 kj/mol。但是该版本高斯在win7平台上无法计算大分子体系，RDX和更小的分子可以稳定输出结果，PETN与HMX及更大的分子则会报错。受此局限性，本文不再继续讨论该方法。
Gaussian1~4系列与W1~4系列虽然精度较高，但耗时过长。如@天泽测试W1U精度时，C2H6分子花费了5个小时才输出结果，可以想象这种效率是无法胜任日常含能材料设计工作的。Gaussian1~4系列从耗时最长的G2到计算较快的G3(MP2)二者相差5倍有余，但同完全基组法一样，在win7平台上的G09W A.02无法得到PETN这样的大分子体系运行结果。因此这两种方法也不在本文后续讨论范围内。
最为广泛讨论和应用的密度泛函方法，准确的说是B3LYP泛函搭配不同的基组。该方法同完全基组法一样，在计算耗时与准确性方面达到了良好的平衡。更令人满意的是安装G09W A.02的win7平台可以稳定运行如CL-20这样的大分子，即该方法可以满足含能材料设计的生成焓仿真需求。
<title></title>
二、基组选择
在参考Rice[2][3]与Politzer[4]的论文后，在综合精度与耗时后选取6-31(d)基组。该基组只针对重原子增加了极化函数，且不引入弥散函数，这点使得6-31(d)基组特别廉价高效，适合个人电脑平台使用。
三、仿真计算
本文选取了28个算例样本，按分子结构分为脂肪烃硝基类、芳烃硝基类、杂环硝胺类、笼形硝胺类、硝酸酯类、叠氮化合物类等；按物态又分为液态、固态。选取的样本原则上要求能查找到生成焓的实验值，经过不同研究组可重复实验值的为佳。这就导致近年来的新型高氮材料几乎都无法纳入比较，本文只能以5-氨基-1-H-四唑、1-H-四氮唑及叠氮氰作为高氮类含能材料的替代样本。
3.1 软件设置与数据处理
同之前一样，本文依然依托ChemBio 3D平台绘制分子结构，然后利用ChemBio 3D内的高斯接口设置高斯运行参数与关键词。
在完成分析后，用记事本打开后缀名为out的高斯输出文件，查找关键词“yes”确保优化后的构型无虚频。更改查找关键词为“scf”，读取文件末尾出的scf参数。将该数据代入公式1中，得到气态生成焓。
<title></title>
图1 高斯输出文件中的scf项
图2 四个“yes”表示几何优化收敛无虚频
<title></title>
ΔH(gas,298 K) = Ei - ∑ njej (1)
表2 拟合参数表
<title></title>
<title></title>
Atom Equivalents | e (hartrees) |
C | -38.121621 |
H | -0.592039 |
N | -54.774096 |
O | -75.161771 |
C' | -38.121380 |
N' | -54.765886 |
O' | -75.157348 |
<title></title>
表2中的C、H、N、O指的是仅参与形成单键的原子，C'、N'、O'指的是参与形成双键或三键的原子。
例：ΔH(硝基苯,gas,298 K)=627.51*（-436.750584875-5*-0.592039 -6*-38.121380-1*-54.765886-2*-75.157348）=11.59 kcal
使用Multiwfn打开ChemBio 3D输出的fchk文件，运行完毕后，读取分子表面积A、电荷平衡度ν、静电势总偏差σtot2(详细操作见：XXXXXXXXXXXXXXXXXXXXXXXX/t/80709)。根据物态的不同分别代入公式2或公式3得到蒸发焓或升华焓。
ΔHvap = aA0.5 + b(νσtot2)0.5 + c (2)
ΔHsub = aA2 + b(νσtot2)0.5 + c (3)
表3 拟合参数表
<title></title>
| a(kcal/mol-Å-1) | b(kcal/mol) | c(kcal/mol) |
Vaporization | 1.818689 | 1.3321583 | 216.142460 |
Sublimation | 4.234303×10-4 | 2.5793785 | 26.7335407 |
<title></title>将以上数据代入公式4或5，即可得到物质的标准生成焓。
ΔH(Liquid,298 K) = ΔH(gas,298 K) - ΔHvap (4)
ΔH(Solid,298 K) = ΔH(gas,298 K) - ΔHsub (5)
<title></title>
3.2 误差分析
<title></title>
表4 算例样本表
<title></title>
| Compound | C | H | O | N | C' | O' | N' | SCF | Hf(298K,g)(kcal) | SA | ν | σtot2 | △H(kcal) | Hf(298K,s)(kcal) | Reference(kcal) | Deviation(kcal) |
1 | Nitroethane | 2 | 5 | 0 | 0 | 0 | 2 | 1 | -284.3280928 | -27.66 | 110.19336 | 0.2299714 | 102.66089 | 9.421714419 | -37.07841627 | -34.4[2] <title></title> | -2.68 <title></title> |
2 | Methyl nitrate | 1 | 3 | 1 | 0 | 0 | 2 | 1 | -320.1894413 | -30.97 | 100.3202 | 0.2491877 | 80.02642 | 8.022411861 | -38.99018515 | -37.23[5] | -1.76 <title></title> |
3 | Nitromethane | 1 | 3 | 0 | 0 | 0 | 2 | 1 | -245.0093332 | -19.46 | 89.5873 | 0.2464283 | 106.93485 | 7.910032638 | -27.37110713 | -26.9[2] | -0.47 <title></title> |
4 | Nitrofuran | 0 | 3 | 1 | 0 | 4 | 2 | 1 | -434.5158721 | -7.46 | 131.93155 | 0.2497634 | 153.89428 | 13.00634359 | -20.46248865 | -24.85[6] | 4.39 |
5 | 1-Nitrobutane | 4 | 9 | 0 | 0 | 0 | 2 | 1 | -362.9551514 | -37.48 | 153.72657 | 0.1824654 | 89.42271 | 11.78791813 | -49.27180591 | -46.03[2] | -3.24 <title></title> |
6 | Nitrobenzene | 0 | 5 | 0 | 0 | 6 | 2 | 1 | -436.7505849 | 11.59 | 149.95022 | 0.2151663 | 160.52117 | 13.95719466 | -2.36576074 | 2.98[2] <title></title> | -5.35 <title></title> |
7 | Nitromethylbenzene | 1 | 7 | 0 | 0 | 6 | 2 | 1 | -476.0616771 | 8.21 | 172.24593 | 0.1993000 | 129.21776 | 14.4868408 | -6.279696758 | -5.46[2] <title></title> | <title></title><title></title><title></title>-0.82 |
8 | NG | 3 | 5 | 3 | 0 | 0 | 6 | 3 | -958.1671991 | -72.22 | 219.57709 | 0.1206177 | 125.32423 | 15.98651057 | -88.20163415 | -88.4[2] | 0.20 <title></title> |
9 | 3,4-Furazandimethanol | 2 | 4 | 3 | 0 | 2 | 4 | 4 | -900.018683 | 8.61 | 216.82424 | 0.1238876 | 211.34416 | 17.45420511 | -8.842279795 | -11.4[2] <title></title> | 2.56 <title></title> |
10 | PETN | 5 | 8 | 4 | 0 | 0 | 8 | 4 | -1316.468542 | -97.08 | 290.2866 | 0.1059363 | 116.06916 | 37.99210431 | -135.0757682 | -128.7[2] <title></title> | -6.38 <title></title> |
11 | Cyanogen azide | 0 | 0 | 0 | 0 | 1 | 0 | 4 | -257.0029926 | 114.16 |
|
|
|
|
| 108[2] <title></title> | 6.16 <title></title> |
12 | Nitrosobenzene | 0 | 5 | 0 | 0 | 6 | 1 | 1 | -361.5397793 | 45.14 |
|
|
|
|
| 48[2] <title></title> | -2.86 <title></title> |
13 | 1,3,5-Trinitrobenzene | 0 | 3 | 0 | 0 | 6 | 6 | 3 | -845.736688 | 5.93 | 201.96146 | 0.1910062 | 109.43278 | 22.33019759 | -16.39709401 | -8.9[2] <title></title> | -7.50 <title></title> |
14 | TATB | 0 | 6 | 0 | 3 | 6 | 6 | 3 | -1011.833117 | 7.17 | 221.09914 | 0.2495522 | 116.37801 | 27.86630248 | -20.69340454 | -33.4[5] <title></title> | 12.71 [2] <title></title> |
15 | 4-Nitroaniline | 0 | 6 | 0 | 1 | 6 | 2 | 1 | -492.1087075 | 16.62 | 163.74721 | 0.2493327 | 336.49932 | 28.24631941 | -11.62701895 | -9.2[2] <title></title> | -2.43 |
16 | RDX | 3 | 6 | 0 | 3 | 0 | 6 | 3 | -897.4083381 | 45.68 | 199.25744 | 0.1898157 | 145.4652 | 23.63192403 | 22.04631414 | 18.9[2] <title></title> | 3.15 |
17 | HMX | 4 | 8 | 0 | 4 | 0 | 8 | 4 | -1196.545772 | 60.08 | 237.30493 | 0.1936042 | 122.65332 | 29.68068742 | 30.39486955 | 24.5[2] <title></title> | 5.89 |
18 | 1-Nitronaphthalene | 0 | 7 | 0 | 0 | 10 | 2 | 1 | -590.3876978 | 31.98 | 193.6462 | 0.2025959 | 159.65969 | 23.8145386 | 8.161610455 | 10.93[2] <title></title> | -2.77 |
19 | 2-Methyl-4,6-dinitrophenol | 1 | 6 | 1 | 0 | 6 | 4 | 2 | -755.7940233 | -43.27 | 200.85197 | 0.2361404 | 93.17833 | 22.44750726 | -65.71638882 | -66.7[2] <title></title> | 0.98 |
20 | α-CL-20 | 6 | 6 | 0 | 6 | 0 | 12 | 6 | -1791.181111 | 143.65 | 318.77201 | 0.0641395 | 211.24676 | 45.78809818 | 97.85955051 | 90.2[2] <title></title> | 7.66 |
21 | Hexanitroethane | 2 | 0 | 0 | 0 | 0 | 12 | 6 | -1306.656829 | 43.87 | 212.25543 | 0.0648325 | 121.01005 | 19.56773388 | 24.29818594 | 25.9[2] <title></title> | -1.60 |
22 | FOX-7 | 0 | 4 | 0 | 2 | 2 | 4 | 2 | -598.3079141 | 7.75 | 148.68958 | 0.1889045 | 462.40696 | 26.73519888 | -18.98048354 | -31.99[5] | 13.01 |
23 | 5-Amino-1-H-tetrazole | 0 | 3 | 0 | 2 | 1 | 0 | 3 | -313.6067212 | 85.73 | 108.42971 | 0.2144450 | 512.63398 | 25.28909108 | 60.44487565 | 49.64[6] | 10.80 |
24 | 1-H-Tetrazole | 0 | 2 | 0 | 1 | 1 | 0 | 3 | -258.2508986 | 79.26 | 93.09171 | 0.1999160 | 401.92447 | 20.05717316 | 59.20566963 | 56.38[6] | 2.83 |
25 | 2,4,6-Trinitroaniline | 0 | 4 | 0 | 1 | 6 | 6 | 3 | -901.1059873 | 3.95 | 208.75047 | 0.2231595 | 122.18479 | 25.18705483 | -21.23957391 | -17.4[2] | -3.84 |
26 | Tetryl | 1 | 5 | 0 | 1 | 6 | 8 | 4 | -1144.856665 | 31.28 | 250.57862 | 0.1794349 | 142.22035 | 32.88364606 | -1.599323085 | 9.8[2] | -11.40 |
27 | NQ | 0 | 4 | 0 | 2 | 1 | 2 | 2 | -409.8537789 | 19.09 | 123.45334 | 0.2150268 | 513.04602 | 26.81174739 | -7.724718134 | -22.2[5] | 14.48 |
28 | NTO | 0 | 2 | 0 | 2 | 2 | 3 | 2 | -521.9752369 | 2.26 | 135.63393 | 0.1721624 | 256.09303 | 18.18319141 | -15.91847449 | -24.1[5] | 8.18 |
29 | ENDA | 2 | 6 | 0 | 2 | 0 | 4 | 2 | -599.4869417 | 11.23 | 168.57101 | 0.1698469 | 231.51106 | 21.4731916 | -10.24682764 | -24.8[5] | 14.55 |
30 | 3-Amino-5-nitro-l,2,4-triazol | 0 | 3 | 0 | 2 | 2 | 2 | 3 | -502.0950324 | 52.96 | 141.51939 | 0.1783582 | 534.33416 | 26.92752481 | 26.02837401 | 14.6[7] | 11.43 |
对9个液态样本，2个气态样本，19个固态样本进行仿真计算后结果见图3。30个材料的平均绝对误差为5.74 kcal/mol，均方根误差（rms）为7.21 kcal/mol，最大偏差14.55 kcal/mol。分析数据可以发现，含-NH2基团的材料仿真计算值偏高，且-NH2基团在分子量中占比大者更容易数据偏高。
一般含能材料分子量约为150~500，爆热840~1550 kcal/kg。假设某新型含能材料分子量为200，爆热1000 kcal。将最大偏差14.55 kcal/mol代入，则误差为7.3%。这样的误差虽然稍高，但该方法用于指导含能材料设计是足以胜任的。
四、后续研究
Rice[3]为了将精度再提高一个层次，使用了6-311++G（2df，2p）基组。通过对35个材料的数据处理，重新拟合了C、H、N、O、C'、N'、O'及a、b、c参数。该方法虽然精度有所提高，但引入了弥散函数后耗时却大幅增加十几倍。在个人电脑平台上，该方法不具有可行性。在经过多种方法对比后，我认为本文的密度泛函B3LYP/6-31G(d)方法对硬件的要求低，能运算大分子体系却耗时短，精度足够高，是最为适应爱好者条件的方法。
Reference
[1]Chemistry of High Energy Materials 2th
[2]Predicting Heats of Formation of Energetic Materials Using Quantum Mechanical Calculations
[3]Improved Prediction of Heats of Formation of Energetic Materials Using Quantum Chemical Calculations
[4]Calculation of heats of sublimation and solid phase heats of formation
[5]explosives 6th
[6]XXXXXXXXXXXXXXXXXXXXXXXX
[7]Synthesis, Properties and Performance of the High Explosive ANTA
致谢：天泽同学，没有与你的共同讨论，这个方法可能就被我搁置了。在发现CBS-4M无法计算大分子体系后，我本意是暂时不想继续了，因为日复一日的查资料、跑程序、验证结果、结果出错、再查论文是一件很费精力的事。在你孜孜不倦的找我讨论下，我只好一遍遍的打开高斯，验证一个个的想法，最后才有了文章里的内容。感谢你的督促与帮助！
[修改于 6年3个月前 - 2018/08/26 15:15:46]
居然有用6-311++G（2df，2p）的，加这么多弥散和极化倒不如提高基组而且据说pople基组构...
准确的说也是用6-31G(d）优化了构型再用6-311++G（2df，2p）处理得到分子表面积A、电荷平衡度ν、静电势总偏差σtot2的，不过即使是这样也会耗费了相当长的时间。
我是从大量论文里挑选方法，自己跑测试对比的。确实有pople基组不如def2的说法，但是奇怪的是pople基组仍然是学术论文中的主流方法，最后我也就沿用了下来。个人看来，Rice和Politzer通过人工拟合参数的方式，在理论基础上通过统计学消除了部分误差，这种方法虽然粗糙，但实际中却是行之有效的。既缩小了误差，又不增加运行时间，这可能是该方法得以流行的一大原因吧。
引用
加载评论中，请稍候...
200字以内，仅用于支线交流，主线讨论请采用回复功能。
PM系列的优势就是快。。。在非超算的基础上还是很有意义的
PM7数秒就能出结果，速度优势太强了。但是Calculating Heat of Formation Values of Energetic Compounds A Comparative Study一文也进行了45个含氮含能材料的计算对比，从文中可以看出PM7的误差还是要较DTF大许多的。如果时间允许的话，自己设计几个含能材料分子在自己电脑上跑一下不算什么困难。若是要大量对比分子结构的不同有何差异趋势时，倒是可以用PM7处理大量材料。
引用
加载评论中，请稍候...
200字以内，仅用于支线交流，主线讨论请采用回复功能。
引用
加载评论中，请稍候...
200字以内，仅用于支线交流，主线讨论请采用回复功能。