量子化学软件仿真(四):计算中性分子的生成焓
dracula14292018/08/26化学 IP:江苏

<title></title>

前三篇文章分论述了含能材料生成焓与密度的仿真计算,其中分子密度仿真介绍的较为透彻,可以达到主流学术论文中的精度水平。但是生成焓采用了PM7半经验法,误差较大,检索文献结果显示应用也较少。所以补齐较高精度的生成焓仿真方法成为了含能材料设计工作的最后一环。 

本文仍然沿用前文中的ChemBioOffice 2014Multiwfn软件,但更新了GaussianG09W A.02版本。 

一、仿真方法 

查阅大量文献可知目前的仿真方法有:1、半经验法AM1PM3MNDO等。2、完全基组法CBS-QB3CBS-4M等。3、密度泛函B3LYP/6-31G(d)4Hartree-Fock法。5Gaussian1~4系列。6W1~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和更小的分子可以稳定输出结果,PETNHMX及更大的分子则会报错。受此局限性,本文不再继续讨论该方法。 

Gaussian1~4系列与W1~4系列虽然精度较高,但耗时过长。如@天泽测试W1U精度时,C2H6分子花费了5个小时才输出结果,可以想象这种效率是无法胜任日常含能材料设计工作的。Gaussian1~4系列从耗时最长的G2到计算较快的G3(MP2)二者相差5倍有余,但同完全基组法一样,在win7平台上的G09W A.02无法得到PETN这样的大分子体系运行结果。因此这两种方法也不在本文后续讨论范围内。 

最为广泛讨论和应用的密度泛函方法,准确的说是B3LYP泛函搭配不同的基组。该方法同完全基组法一样,在计算耗时与准确性方面达到了良好的平衡。更令人满意的是安装G09W A.02win7平台可以稳定运行如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中,得到气态生成焓。

QQ截图20180826003111.jpg

<title></title>

1 高斯输出文件中的scf

QQ截图20180826003018.jpg

图2 四个“yes”表示几何优化收敛无虚频

<title></title>

ΔH(gas,298 K) = Ei - ∑ njej                              (1)

表2 拟合参数表

<title></title>

<title></title>

Atom Equivalents 

e (hartrees) 

-38.121621 

-0.592039 

-54.774096 

-75.161771 

C' 

-38.121380 

N' 

-54.765886 

O' 

-75.157348 

<title></title>

2中的CHNO指的是仅参与形成单键的原子,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(详细操作见:https://www.kechuang.org/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>将以上数据代入公式45,即可得到物质的标准生成焓。 

Δ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' 

O' 

N' 

SCF 

Hf(298K,g)(kcal) 

SA 

ν 

σtot2 

△H(kcal) 

Hf(298K,s)(kcal) 

Reference(kcal) 

Deviation(kcal) 

Nitroethane 

-284.3280928 

-27.66 

110.19336 

0.2299714 

102.66089 

9.421714419 

-37.07841627 

-34.4[2]

<title></title>

-2.68 

<title></title>

Methyl nitrate 

-320.1894413 

-30.97 

100.3202 

0.2491877 

80.02642 

8.022411861 

-38.99018515 

-37.23[5]

-1.76

<title></title>

Nitromethane 

-245.0093332 

-19.46 

89.5873 

0.2464283 

106.93485 

7.910032638 

-27.37110713 

-26.9[2]

-0.47

<title></title>

Nitrofuran 

-434.5158721 

-7.46 

131.93155 

0.2497634 

153.89428 

13.00634359 

-20.46248865 

-24.85[6]

4.39

1-Nitrobutane 

-362.9551514 

-37.48 

153.72657 

0.1824654 

89.42271 

11.78791813 

-49.27180591 

-46.03[2]

-3.24

<title></title>

Nitrobenzene 

-436.7505849 

11.59 

149.95022 

0.2151663 

160.52117 

13.95719466 

-2.36576074 

2.98[2]

<title></title>

-5.35

<title></title>

Nitromethylbenzene 

-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

NG 

-958.1671991 

-72.22 

219.57709 

0.1206177 

125.32423 

15.98651057 

-88.20163415 

-88.4[2]

0.20

<title></title>

3,4-Furazandimethanol
dinitrate 

-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 

-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 

-257.0029926 

114.16 

 

 

 

 

 

108[2]

<title></title>

6.16

<title></title>

12 

Nitrosobenzene 

-361.5397793 

45.14 

 

 

 

 

 

48[2]

<title></title>

-2.86

<title></title>

13 

1,3,5-Trinitrobenzene 

-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 

-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 

-492.1087075 

16.62 

163.74721 

0.2493327 

336.49932 

28.24631941 

-11.62701895 

-9.2[2]

<title></title>

-2.43 

16 

RDX 

-897.4083381 

45.68 

199.25744 

0.1898157 

145.4652 

23.63192403 

22.04631414 

18.9[2]

<title></title>

3.15 

17 

HMX 

-1196.545772 

60.08 

237.30493 

0.1936042 

122.65332 

29.68068742 

30.39486955 

24.5[2]

<title></title>

5.89 

18 

1-Nitronaphthalene 

10 

-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 

-755.7940233 

-43.27 

200.85197 

0.2361404 

93.17833 

22.44750726 

-65.71638882 

-66.7[2]

<title></title>

0.98 

20 

α-CL-20 

12 

-1791.181111 

143.65 

318.77201 

0.0641395 

211.24676 

45.78809818 

97.85955051 

90.2[2]

<title></title>

7.66 

21 

Hexanitroethane 

12 

-1306.656829 

43.87 

212.25543 

0.0648325 

121.01005 

19.56773388 

24.29818594 

25.9[2]

<title></title>

-1.60 

22 

FOX-7 

-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 

-313.6067212 

85.73 

108.42971 

0.2144450 

512.63398 

25.28909108 

60.44487565 

49.64[6]

10.80 

24 

1-H-Tetrazole 

-258.2508986 

79.26 

93.09171 

0.1999160 

401.92447 

20.05717316 

59.20566963 

56.38[6]

2.83 

25 

2,4,6-Trinitroaniline 

-901.1059873 

3.95 

208.75047 

0.2231595 

122.18479 

25.18705483 

-21.23957391 

-17.4[2]

-3.84 

26 

Tetryl 

-1144.856665 

31.28 

250.57862 

0.1794349 

142.22035 

32.88364606 

-1.599323085 

9.8[2]

-11.40 

27 

NQ 

-409.8537789 

19.09 

123.45334 

0.2150268 

513.04602 

26.81174739 

-7.724718134 

-22.2[5]

14.48 

28 

NTO 

-521.9752369 

2.26 

135.63393 

0.1721624 

256.09303 

18.18319141 

-15.91847449 

-24.1[5]

8.18 

29 

ENDA 

-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 

-502.0950324 

52.96 

141.51939 

0.1783582 

534.33416 

26.92752481 

26.02837401 

14.6[7]

11.43 

9个液态样本,2个气态样本,19个固态样本进行仿真计算后结果见图330个材料的平均绝对误差为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++G2df2p)基组。通过对35个材料的数据处理,重新拟合了CHNOC'、N'、O'及abc参数。该方法虽然精度有所提高,但引入了弥散函数后耗时却大幅增加十几倍。在个人电脑平台上,该方法不具有可行性。在经过多种方法对比后,我认为本文的密度泛函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无法计算大分子体系后,我本意是暂时不想继续了,因为日复一日的查资料、跑程序、验证结果、结果出错、再查论文是一件很费精力的事。在你孜孜不倦的找我讨论下,我只好一遍遍的打开高斯,验证一个个的想法,最后才有了文章里的内容。感谢你的督促与帮助!

[修改于 5年8个月前 - 2018/08/26 15:15:46]

+1  学术分    博丽灵梦    2018/08/27 理论高度指导实践
来自:数理化 / 化学
13
1
已屏蔽 原因:{{ notice.reason }}已屏蔽
{{notice.noticeContent}}
~~空空如也
dracula1429 作者
5年8个月前 修改于 5年8个月前 IP:江苏
850050

<title></title>

此楼可删

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
dracula1429作者
5年8个月前 IP:江苏
850052

最后加个几种方法的直观对比

图片1_副本.png

 

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
哈哈哈哈哈士奇
5年8个月前 IP:辽宁
850053
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
ne555
5年8个月前 IP:山东
850080

居然有用6-311++G2df2p)的,加这么多弥散和极化倒不如提高基组

而且据说pople基组构建的不太好,不如用def2

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
St_Maxwell
5年8个月前 IP:江西
850085
引用ne555发表于4楼的内容
居然有用6-311++G(2df,2p)的,加这么多弥散和极化倒不如提高基组而且据说pople基组构...

很多文献是这样的,肯定不能盲目学。

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
博丽灵梦
5年8个月前 IP:北京
850112

PM系列的优势就是快。。。在非超算的基础上还是很有意义的

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
天泽
5年8个月前 修改于 5年8个月前 IP:美国
850115

顶一下  KC.png
 顺便说一下,含氧化氮(N→O) 组的化合物,比如在 TTTO,LAX-112,LLM-105 里,在计算气态生成焓的时候,使用普通 N 和 O原子参数。原因是 N→O 里的 N 不一定是饱和键。个人使用这个方法的时候,使用N' 和 O' 或者 N' 和 O 的参数的误差比 N 和 O 高好几倍。

请看 Politzer, Peter, Pat Lane, and Jane S. Murray. "Some interesting aspects of N-oxides." Molecular Physics 112.5-6 (2014): 719-725.

KC.png

 

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
dracula1429作者
5年8个月前 修改于 5年8个月前 IP:江苏
850235
引用ne555发表于4楼的内容
居然有用6-311++G(2df,2p)的,加这么多弥散和极化倒不如提高基组而且据说pople基组构...

准确的说也是用6-31G(d)优化了构型再用6-311++G(2df,2p)处理得到分子表面积A、电荷平衡度ν静电势总偏差σtot2的,不过即使是这样也会耗费了相当长的时间。

我是从大量论文里挑选方法,自己跑测试对比的。确实有pople基组不如def2的说法,但是奇怪的是pople基组仍然是学术论文中的主流方法,最后我也就沿用了下来。个人看来,Rice和Politzer通过人工拟合参数的方式,在理论基础上通过统计学消除了部分误差,这种方法虽然粗糙,但实际中却是行之有效的。既缩小了误差,又不增加运行时间,这可能是该方法得以流行的一大原因吧。

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
dracula1429作者
5年8个月前 IP:江苏
850236
引用博丽灵梦发表于6楼的内容
PM系列的优势就是快。。。在非超算的基础上还是很有意义的

PM7数秒就能出结果,速度优势太强了。但是Calculating Heat of Formation Values of Energetic Compounds A Comparative Study一文也进行了45个含氮含能材料的计算对比,从文中可以看出PM7的误差还是要较DTF大许多的。如果时间允许的话,自己设计几个含能材料分子在自己电脑上跑一下不算什么困难。若是要大量对比分子结构的不同有何差异趋势时,倒是可以用PM7处理大量材料。

QQ截图20180828221002.jpg

 

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
dracula1429作者
5年8个月前 IP:江苏
850238
引用天泽发表于7楼的内容
顶一下  顺便说一下,含氧化氮(N→O) 组的化合物,比如在 TTTO,LAX-112,LLM-10...

忘了算算含N→O材料了,算例不够全面😅

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
frankmlin
2年0个月前 IP:山西
902218

例:ΔH(硝基苯,gas,298 K)=627.51*(-436.750584875-5*-0.592039 -6*-38.121380-1*-54.765886-2*-75.157348)=11.59 kcal 

问下,式子里面的627.51是做什么用的呢?我怎么记得从hatree到kcal是要乘上一个2625.5呢?求解答



引用
评论(1)
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
虚心求教
1年0个月前 IP:山西
918594

老师您好,文中的两个拟合参数表(CHON的能量、蒸发升华焓公式系数)是如何得到的数据?


引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
哈哈哈哈哈士奇
1年0个月前 IP:陕西
918611
引用虚心求教发表于12楼的内容
老师您好,文中的两个拟合参数表(CHON的能量、蒸发升华焓公式系数)是如何得到的数据?

中北的?

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论

想参与大家的讨论?现在就 登录 或者 注册

所属专业
所属分类
上级专业
同级专业
dracula1429
专家 学者 机友 笔友
文章
108
回复
1069
学术分
25
2006/07/09注册,6时13分前活动
暂无简介
主体类型:个人
所属领域:无
认证方式:手机号
IP归属地:未同步
文件下载
加载中...
{{errorInfo}}
{{downloadWarning}}
你在 {{downloadTime}} 下载过当前文件。
文件名称:{{resource.defaultFile.name}}
下载次数:{{resource.hits}}
上传用户:{{uploader.username}}
所需积分:{{costScores}},{{holdScores}}下载当前附件免费{{description}}
积分不足,去充值
文件已丢失

当前账号的附件下载数量限制如下:
时段 个数
{{f.startingTime}}点 - {{f.endTime}}点 {{f.fileCount}}
视频暂不能访问,请登录试试
仅供内部学术交流或培训使用,请先保存到本地。本内容不代表科创观点,未经原作者同意,请勿转载。
音频暂不能访问,请登录试试
支持的图片格式:jpg, jpeg, png
插入公式
评论控制
加载中...
文号:{{pid}}
投诉或举报
加载中...
{{tip}}
请选择违规类型:
{{reason.type}}

空空如也

加载中...
详情
详情
推送到专栏从专栏移除
设为匿名取消匿名
查看作者
回复
只看作者
加入收藏取消收藏
收藏
取消收藏
折叠回复
置顶取消置顶
评学术分
鼓励
设为精选取消精选
管理提醒
编辑
通过审核
评论控制
退修或删除
历史版本
违规记录
投诉或举报
加入黑名单移除黑名单
查看IP
{{format('YYYY/MM/DD HH:mm:ss', toc)}}