量子化学软件计算:以三肼基均三嗪的密度、生成焓为例

文 号

825022

5 回复 / 1396 浏览


dracula14291 年前 -2016-08-29 23:03825022
在当前的含能材料研究中,通过Gaussian等量子化学软件对目标分子的结构进行优化并计算密度、生成焓等参数是设计、筛选含能材料的必要步骤。本文将使用ChemBioOffice 2014、MOPAC、Gaussian 03w、Multiwfn软件,以三肼基均三嗪(THT)这一物质为例进行简单的计算。不过本人也只是刚刚涉足量子化学软件领域,许多地方做的不够好甚至错误百出,希望有朋友能指出问题,大家多多探讨共同进步。
一、ChemDraw2014
相对于Chem3D我更偏爱ChemDraw一些,因为ChemDraw更简单易懂。Chem3D则主要用来进行量子化学计算,使用ChemDraw来完成分子结构图的绘制。
1、选择左侧工具栏中的Cyclohexane Ring(环己烷),绘制出一个简单的环己烷结构。
2、再点击Multiple Bonds(多键),把环己烷修改成苯环。
3、使用Text(文本)功能,单击需要更改部位,将C原子变为N原子,这样就完成了均三嗪骨架的绘制。
4、点击Solid Bond(单键),绘制出三乙基均三嗪的结构。
5、依照3将乙基上的C原子更改为N原子,使乙基变为肼基。
6、将文件输出为THT.cdx。
269159

图1.1 经ChemDraw绘制出的THT结构图
7、以Chem3D打开THT.cdx,从Claculations里选择Mopac Interface模块,不进行修改使用默认的PM7半经验算法计算THT在标准状态下的气相生成焓,为99.01 Kcal/mol。
269160

图1.2 计算
269161

图1.3 默认的参数
8、将文件再次另存,保存为THT.gjf,为下一步导入Gaussian做准备。
二、Gaussian 03w
使用高斯是为了将分子结构优化,并输出为fch文件,使之与Multiwfn进行计算并得到最终结果。
1、用Gview打开上一步得到的THT.gjf文件,点击Claculate并设置B3LYP 方法,在6 –31G** 基组上对结构进行优化。并在Link 0选项里设置输出chk文件。
269162

图2.1 设置基组
269164

图2.2 设置输出chk文件
2、点击submit,启动Gaussian 03w,运算开始。运算结束后得到THT2.LOG和THT.chk两个文件。
3、关闭Gaussian 03w程序,并再次打开,使用Utilities下的Formchk功能,将chk文件转换成Multiwfn可以读取的fch文件(如图2.3)。
269165

图2.3 将chk转为fch
三、Multiwfn
以上如此繁杂的步骤只是为了得到Multiwfn可以运行fch文件,使用Multiwfn计算出分子表面积(\(A\))、电荷平衡度(\(V\))和静电势总偏差(\(σ_\text{tot}\)),然后用Politzer公式(1)计算出蒸发焓。
$$ΔH_\text{sub} = \alpha A^{2+}\beta(vσ_\text{tot}^2)+\gamma \tag{1}$$再根据式(2)即可计算得固态标准摩尔生成焓。
$$ΔH_\text{f}\text{(solid,298 K)} = ΔH_\text{f}\text{(gas,298 K)} - ΔH_\text{sub}\text{(298 K)} \tag{2}$$1、右键点击THT.fch,将其默认打开方式设置成Multiwfn。
2、双击THT.fch运行Multiwfn,出现图3.1所示界面。
269166

图3.1 Multiwfn运行界面
3、输入12并回车,再输入0运行计算,结果见图3.2。
269167

图3.2 Multiwfn运行结果
其中\(A = 204.42256\ \mathrm{Angstrom^2}\),\(vσ_\text{tot}^2 = 33.7538833\ \mathrm{(kcal/mol)^2}\)
得\(ΔH_\text{sub}\text{(298 K)}= 23.71\ \mathrm{kcal/mol}\)
$$ΔH_\text{f}\text{(solid,298 K)}  = 99.01-23.71 = 75.3\ \mathrm{kcal/mol}$$4、密度则根据一个THT分子的体积是\((1×10^{-8})^3×A\ \mathrm{cm}^3\),质量为\(M/(6.02214179 ×10^{23})\ \mathrm{g}\),根据式(3)得:
$$ρ = m/V \tag{3}$$$$ρ_\text{THT} = 171/(6.02214179 × 10^{23})/(1×10^{-24}×A) = 1.389\ \mathrm{g/cm^3}$$四、验证
将5-氨基四唑的文献值与计算值进行对比,以验证本文方法在计算氮杂环化合物时的准确性。                            
\(
\begin{array}{|c|c|c|c|}
\hline
  & ΔH_\text{sub}\text{(298 K)}\mathrm{(kcal/mol)} & ΔH_\text{f}\text{(solid,298 K)}\mathrm{(kcal/mol)} & ρ\mathrm{(g/cm^3)} \\
\hline
\text{Cal} & 20.01 & 59.96 & 1.246 \\
\hline
\text{Exp} & & 49.93 & 1.711 \\
\hline
\end{array}
\)
                                              表4-1 5-氨基四唑的计算与实测对比
经过对比可以发现生成焓计算较为符合实测值误差仅\(10.03\ \mathrm{kcal/mol}\),但是密度的误差较大达到了\(0.465\ \mathrm{g/cm^3}\)。但是本文仅选用了最为简单的方法,若是使用更为复杂的式(4),则计算得\(1.405\ \mathrm{g/cm^3}\)。
$$ρ = αm/V + β(vσ_\text{tot}^2) + γ \tag{4}$$$$α = 0.9183$$$$β = 0.0028$$$$γ = 0.0443$$五、结语
本文只是想把整个量子化学软件使用的流程展现出来,其中肯定会有着诸多的错误、疏忽,得到的计算结果也偏离实测很多。毕竟我并不了解量子化学,也不懂得使用最恰当的基组和相关参数,只是凭着多看了几篇文献,照猫画虎的采用了文献上采用的方法。再一个软件相互对接的过程中也出现过许多的错误,比如用Chem3D画的好好的图,输入到了Gview中后总是会出现单键变双键甚至化学键消失变成了孤立的原子或者基团的情况。这些问题不解决,就谈不上精确性,更谈不上指导实际工作。以我一人之力是需要很长时间才能完善这个方法的,我想求助坛友们,能共同探索,使该方法成熟并能应用于实际,使论坛除了在合成测试之外开辟出计算机仿真的路径来,将自身能力、将论坛建设推向更高的台阶。

[修改于 1 年前 - 2016-08-31 15:15:39]


718281828kc1 年前 -2016-08-29 23:07825023
过程很详尽,对于这种基组和算法,结果偏差已经算是在合理范围了。。
至于化学键消失的问题,其实不用怎么在意,键连是通过计算机判断键长来确定的,所以可能会因此判断为无化学键

dracula14291 年前 -2016-09-01 18:56825103
引用 718281828kc:
过程很详尽,对于这种基组和算法,结果偏差已经算是在合理范围了。。
至于化学键消失的问题,其实不用怎么在意,键连是通过计算机判断键长来确定的,所以可能会因此判断为无化学键
原来是这个原因,我一直以为是软件对接出错了。。。下一步就是看看把误差再缩小一点达到能指导应用的地步。

ne5553 个月前 -2017-04-15 20:41833140
再次看帖子发现一点问题,密度计算误差如此大是否是基组的问题,因为三个常数是在B3PW91/6-31G**下拟合出的

last-order3 个月前 -2017-04-15 23:40833150

我也测试了一下。分别使用B3LYP/6-31G**和M06-2X/def2TZVP进行几何优化并计算5-氨基四唑的密度。结果如下:

B3LYP/6-31G**

Volume: 630.00144 Bohr^3 ( 93.35660 Angstrom^3)
Estimated density according to mass and volume: 1.5131 g/cm^3
Minimal value: -36.60195 kcal/mol Maximal value: 54.97392 kcal/mol
Overall surface area: 387.10628 Bohr^2 ( 108.40082 Angstrom^2)
Positive surface area: 173.08492 Bohr^2 ( 48.46872 Angstrom^2)
Negative surface area: 214.02137 Bohr^2 ( 59.93210 Angstrom^2)
Overall average value: 0.00029148 a.u. ( 0.18291 kcal/mol)
Positive average value: 0.03376935 a.u. ( 21.19061 kcal/mol)
Negative average value: -0.02678298 a.u. ( -16.80659 kcal/mol)
Overall variance (sigma^2_tot): 0.00075443 a.u.^2 ( 297.06915 (kcal/mol)^2)
Positive variance: 0.00052564 a.u.^2 ( 206.98141 (kcal/mol)^2)
Negative variance: 0.00022878 a.u.^2 ( 90.08774 (kcal/mol)^2)
Balance of charges (miu): 0.21129146
Product of sigma^2_tot and miu: 0.00015940 a.u.^2 ( 62.76817 (kcal/mol)^2)
Internal charge separation (Pi): 0.02993833 a.u. ( 18.78660 kcal/mol)

0.9138*1.5131+0.0028*62.76817+0.0443=1.603

M06-2X/TZVP

Volume: 639.89024 Bohr^3 ( 94.82197 Angstrom^3)
Estimated density according to mass and volume: 1.4897 g/cm^3
Minimal value: -36.45772 kcal/mol Maximal value: 58.45159 kcal/mol
Overall surface area: 392.17047 Bohr^2 ( 109.81893 Angstrom^2)
Positive surface area: 184.59652 Bohr^2 ( 51.69230 Angstrom^2)
Negative surface area: 207.57395 Bohr^2 ( 58.12663 Angstrom^2)
Overall average value: 0.00158584 a.u. ( 0.99513 kcal/mol)
Positive average value: 0.03332935 a.u. ( 20.91450 kcal/mol)
Negative average value: -0.02664382 a.u. ( -16.71926 kcal/mol)
Overall variance (sigma^2_tot): 0.00078132 a.u.^2 ( 307.66095 (kcal/mol)^2)
Positive variance: 0.00056427 a.u.^2 ( 222.19344 (kcal/mol)^2)
Negative variance: 0.00021705 a.u.^2 ( 85.46751 (kcal/mol)^2)
Balance of charges (miu): 0.20062615
Product of sigma^2_tot and miu: 0.00015675 a.u.^2 ( 61.72483 (kcal/mol)^2)
Internal charge separation (Pi): 0.02990414 a.u. ( 18.76514 kcal/mol)

0.9138*1.4897+0.0028*61.72483+0.0443=1.578

两个结果与实验值的差别不是很大。我在想是不是楼主在计算的时候出了什么纰漏……

P.S. 英文维基上查到的5-氨基四唑的密度是1.502 g/cm^3……

[修改于 3 个月前 - 2017-04-16 12:26:51]


dracula14293 个月前 -2017-04-16 00:12833152
这篇帖子是我当时的初步摸索,后续内容见https://bbs.kechuang.org/t/80801,改正了许多错误后达到了同论文相差无几的精度。

返回 化学化工
返回 本页顶部

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


nkc Development Server https://github.com/lovetheory/nkc2

科创研究院 (c)2005-2017

蜀ICP备11004945号-2 川公网安备51010802000058号