请教有关计算生成焓的一些问题

kimariyb

各位老师好，学生是刚入门 3 个月的萌新。学生在学习初级量化班的热力学量的计算这一部分时，想要进一步学习如何计算生成焓。于是学生参考了公社里有关生成焓计算的帖子，结合不同的理论方法计算了六个小分子，在 298.15 K 下的生成焓。所使用的程序以及对应的版本为：

1. Gaussian 16 C01（PS. Gaussian 16 A03 似乎在算 G4、G4(MP2) 这类 G4 开头的热力学组合方法时会有 bug，即 out 文件末尾热力学数据输出不全。）

2. Orca 5.0.4

3. Shermo 2.3.6

4. EasyShermo （PS. 学生为了偷懒，自行编写的批量执行 Shermo 的 Python 脚本）

5. Multiwfn 3.8(dev) （PS. 用于生成 Orca 输入文件 .inp）

表的上半部分是用 G16 结合 CBS-QB3、G4(MP2)、G4 三种热力学组合方法计算得到的生成焓数据。计算步骤为：

1. 首先分别用 CBS-QB3、G4(MP2)、G4 计算 N2(g)、H2(g)、Cl2(g)、F2(g)、O2(g)、C(g)，得到各单质的 H(T)。其中 C(g) 和 O2(g) 都计算的三重态。

2. 接着分别用 CBS-QB3、G4(MP2)、G4 计算 NH3(g)、HF(g)、HCl(g)、H2O(g)、CO2(g)、CH4(g)，得到生成物的 H(T)。

3. 最后按照生成焓的定义，按照反应系数，将生成物的 H(T) 减去对应单质的 H(T) 得到 ΔH(T)。其中 CO2(g) 和 CH4(g)，最后再加上石墨的升华焓 715 kJ/mol。
接着学生又使用了常规方法，即中低级别做优化振动分析，高级别做高精度单点的方法计算热力学量。计算步骤为：

1. 首先在 B3LYP/6-311G** 的级别下对所有涉及到的单质、分子进行优化和振动分析。其中 C(g) 只在 B3LYP/6-311G** 级别下做振动分析。

2. 在已经优化好的单质、分子结构下，分别在 PWPB95-D3(BJ)/def2-QZVPP、revDSD-PBEP86-D3(BJ)/def2-QZVPP 以及 CCSD(T)/cc-PVQZ 三种级别下计算高精度单点。其中 PWPB95-D3(BJ)/def2-QZVPP 级别的单点是在 Orca 下计算的，其余两种级别下的单点都是在 Gaussian 计算的。

3. 找到 B3LYP/6-311G** 级别下的 ZPE 矫正因子，并通过 Shermo 分别计算这些单质、分子的热力学量。由于涉及到的文件比较多，所以学生通过自己开发的 Python 脚本 EasyShermo 通过批量调用
Shermo xxx.out -E xxx -ilowfreq 2 -sclZPE 0.9882 命令计算每个单质、分子的热力学量，并将 Shermo 的输出结果保存在同名 txt 文件中。有关 EasyShermo 的具体信息可以浏览 https://github.com/kimariyb/easy-Shermo

4. 最后也和热力学组合方法一样，按照反应系数，将生成物的 H(T) 减去对应单质的 H(T) 得到 ΔH(T)。其中 CO2(g) 和 CH4(g)，最后再加上石墨的升华焓 715 kJ/mol。

  

最后学生将计算得到的数据绘制成了上表，对此学生有几个问题想问一下各位老师：

1. 学生计算生成焓的步骤是否合理，如有不合理之处，还请老师们指出。

2. 如果计算合理，例如 CO2(g) 的生成焓，使用热力学组合方法和常规方法计算的误差都在 4 kJ/mol 也就是 1 kcal/mol 以上，这种误差是否可以接受呢。

3. 例如 CO2(g)、CH4(g) 这种尽管用上了 CCSD(T)/cc-PVQZ 误差还是在 1 kcal/mol 以上，如果想要继续提高精度，应该使用什么方法呢。

感谢各位老师耐心看完，还请赐教。Gaussian、Orca 计算的输入输出文件以及 Shermo 的输出信息均已附在附件中。

20230723.zip (5.51 MB, 下载次数 Times of downloads: 24)

2023-7-27 02:08 上传 Uploaded

点击下载
Click to download

dantevinsky

你这个生成焓差1kcal/mol以上的都好多啊

kimariyb

dantevinsky 发表于 2023-7-27 08:47
你这个生成焓差1kcal/mol以上的都好多啊

是的，所以发帖问问各位老师计算流程是否合理

wzkchem5

误差在这个量级还是比较正常的。如果要提高精度的话，单点能可以考虑用CCSD(T) cc-pVQZ、cc-pV5Z外推到基组极限，或者用CCSD(T)-F12/cc-pVQZ-F12；ZPE和焓校正可以考虑用VPT2考虑非谐性效应

kimariyb

wzkchem5 发表于 2023-7-27 16:37
误差在这个量级还是比较正常的。如果要提高精度的话，单点能可以考虑用CCSD(T) cc-pVQZ、cc-pV5Z外推到基组 ...

好的，感谢王老师的解答

sobereva

CO2这种微小体系可以用W1热力学组合方法算焓，是Gaussian里算热力学量能做到的最精确的结果，原子化能误差平均明显小于1 kcal/mol。注意实验的生成焓、升华焓本身误差就不小（相对于极高精度量子化学计算而言），不要对相符程度期望得过高。极高精度热力学计算（亚化学精度层次）的精度往往比热力学实验还要高。至于计算细节流程，我目前没时间细看。

kimariyb

sobereva 发表于 2023-7-28 02:04
CO2这种微小体系可以用W1热力学组合方法算焓，是Gaussian里算热力学量能做到的最精确的结果，原子化能误差 ...

感谢卢老师的回答，学生了解了

kimariyb

sobereva 发表于 2023-7-28 02:04
CO2这种微小体系可以用W1热力学组合方法算焓，是Gaussian里算热力学量能做到的最精确的结果，原子化能误差 ...

学生在卢老师提醒后用 W1BD 方法算了一下 NH3 的生成焓，和实验的误差只有 0.5 kcal/mol 左右，比 G4 方法算的误差小了 1 kcal/mol 左右。学生还想问一下，如果用常规方法计算一些比较大的体系的生成焓（热力学组合方法算不动的情况），误差在 1 ~ 2 kcal/mol 是否也可以接受呢。

sobereva

kimariyb 发表于 2023-7-28 11:54
学生在卢老师提醒后用 W1BD 方法算了一下 NH3 的生成焓，和实验的误差只有 0.5 kcal/mol 左右，比 G4 方 ...

有机体系误差能在一点几kcal/mol就算不错了。涉及过渡金属的反应的误差普遍更大