笔者同时在计算化学公社更新了帖子《EasyShermo: 萌新的全自动用 Shermo 做批处理的脚本》一文，现已被计算化学公社社长卢天收录为计算化学公社精华帖以及 2023 年 7 月计算化学网文 hightlight。
地址为：http://bbs.keinsci.com/forum.php?mod=viewthread&tid=38871&fromuid=47478
前言
热力学是物理化学中非常重要的一门学科，它研究的是物质的能量转化和热力学性质，如热力学平衡、热力学势等。在化学领域中，热力学的应用非常广泛，例如用来解释化学反应的热力学特性、计算热力学稳定性、预测化学反应的平衡常数等。
在热力学中，热力学量极其重要，热力学量是热力学中用来描述物质热力学性质的物理量，它们反映了物质在不同状态下的热力学特性。基于量子化学的热力学计算方法可以通过计算这些热力学量，从而预测化学反应的热力学性质，如反应热、熵、自由能等。本文只介绍气相下热力学量的计算方法并基于理想气体假设。不介绍如何计算溶液环境下的热力学量。
在量子化学中，常见的热力学量的计算方法主要有两种，一种是常规方法，另外一种是热力学组合方法。本文将介绍如何使用这两种方法计算热力学量。本文中所有使用的程序以及版本如下所示，所使用的计算环境为 Rocky Linux 8.8 64bit。
Gaussian 16 A.03 Linuxorca 5.0.4 LinuxShermo 2.3.6EasyShermo 1.1.0
注！Gaussian 16 A.03 在使用 G4、G4MP2 计算热力学量时会有 Bug，后改为 Gaussian 16 C.01
常规方法
气相下计算热力学量的常规合理的方法就是，先用中或中低级别优化结构并做振动分析得到热力学矫正量，再用高级别在此结构下算一次高精度单点，最后两者相加。
笔者一般使用 Gaussian 做几何优化、振动分析得到热力学矫正量，再使用 Orca 做一次高精度单点。对于小体系（1~8 原子以内），可以在 6-311G** 或 def-TZVP 基组下结合对当前体系合适的泛函做几何优化和振动分析（实在觉得计算资源有多，用 def2-TZVP 也不是不行）。接着用 Orca 在 Approximated CCSD(T)/CBS with help of MP2 (cc-pVTZ->QZ extrapolation) 级别下做单点任务。如果用 Gaussian 算高精度单点的话，可以直接用 CCSD(T)/cc-pVQZ，或者是差一点的 revDSD-PBEP86-D3(BJ)/def2-QZVPP，在 Gaussian 16 里可以通过以下关键词使用这个级别。
1 | # DSDPBEP86/def2QZVPP IOp(3/125=0079905785,3/78=0429604296,3/76=0310006900,3/74=1004) |
如果是稍微大一点的体系，如果照用以上的级别就不太好，计算时间就是一个问题。笔者比较倾向用 6-311G* 或 6-311G** 基组结合一个合适的泛函做几何优化和振动分析。接着用 Orca 在 RI-PWPB95-D3(BJ)/def2-QZVPP 或者 DLPNO-CCSD(T)/cc-pVTZ with normalPNO and RIJK 级别算单点。如果用 Gaussian 在 CCSD(T)/cc-pVTZ 级别下算稍微大一点体系的单点，很难算得动。但是在 Orca 中，则可以使用上述的黑科技级别计算单点，相比于 CCSD(T)/cc-pVTZ 级别便宜许多，也可以试试 CCSD(T)-F12/cc-pVDZ-F12 with RI 这个级别，不过这些级别除了 CCSD(T)/cc-pVTZ 都只能在 Orca 里用，这也是笔者比较倾向用 Orca 算单点的原因。
如果觉得 RI-PWPB95-D3(BJ)/def2-QZVPP 的结果不太满意，可以试试用 RI-revDSD-PBEP86-D4/def2-TZVPP 级别算。如果只用 Gaussian 算大一点体系的单点，revDSD-PBEP86-D3(BJ)/def2-TZVPP 是很好的选择了。如果计算资源实在不忍直视，那就只能降低到用 6-31G* 基组做几何优化和振动分析，用 M06-2X/def2-TZVP 算单点。
以上的泛函和基组的选择并不一定满足所有情况，如果想要判断自己选择的基组或泛函是否合理，可以浏览 Sobereva 的博文《谈谈量子化学中基组的选择》和《简谈量子化学计算中 DFT 泛函的选择》。
热力学组合方法
热力学组合方法是量子化学计算热力学量的另外一种方法，之所以成为热力学组合，是因为这是一系列方法的组合。整套过程自动完成，不需要提前做几何优化、振动分析。在 Gaussian 中使用热力学组合方法来算热力学量，只用写热力学组合方法的名字就可以了，不需要写其他关键词。常见的热力学组合方法名称、精度和耗时在 Sobereva 的帖子《常见热力学组合方法精度、耗时的关系》里进行了说明。
根据 Sob 老师的帖子，性价比比较高的热力学组合方法分别为 CBS-QB3、G4(MP2) 以及 G4。除此之外，G4(MP2)-6X 也很值得使用，不过要在 Gaussian 中使用 G4(MP2)-6X 必须使用一个模板：
1 | %chk=mol.chk |
任务结束后，再将此 Perl 脚本：G4MP2_6x.pl 和 Gaussian 输出文件文件（假设为 G4MP2_6x.out）都放在一个文件夹下，运行 G4MP2_6x.pl G4MP2_6x.out 命令，Perl 脚本就会自动把相关数据从 Gaussian 输出文件中提取出来并进行处理，默认是在标况下算的。如果是 G4、CBS-QB3 这类的热力学组合方法计算出来的焓、自由能都可以直接在 Gaussian 输出文件中得到。
1 | Temperature= 298.150000 Pressure= 1.000000 |
EasyShermo 批处理
在介绍 EasyShermo 前，首先介绍一下 Shermo。Shermo 是 Sobereva 开发的一个免费的可以独立运行的计算分子热力学数据的程序，Shermo 在热力学数据计算方面既强大又灵活方便，计算结果也比其他量化程序可靠很多。如果想要了解 Shermo 的更多信息，可以浏览：
- Shermo 的官方网站：Shermo
- Shermo 的中文教程：使用Shermo结合量子化学程序方便地计算分子的各种热力学数据
- Shermo 程序的原文：Shermo: A general code for calculating molecular thermochemistry properties
EasyShermo 是 Kimariyb 开发的一款全自动批处理 Shermo 的 Python 脚本。目前已经在 Github 上开源，可以点此访问 EasyShermo 的 Github 主页。可以根据主页的 README 进行安装和使用，EasyShermo 可以轻松的瞬间调用 Shermo 批量处理几十、甚至上百个 Gaussian 或 Orca 的输出文件，并生成每一个文件对应的 Shermo 输出结果。
当前版本的 EasyShermo 仅支持读取 Gaussian 和 Orca 的单点输出文件，其中 Gaussian 只支持杂化、双杂化泛函以及 CCSD(T) 计算的单点输出文件，Orca 则全部支持。下一部分，笔者将使用 Gaussian、Orca、Shermo 以及 EasyShermo 结合不同计算级别，计算几个简单分子的热力学量，并对这些不同方法进行对比。
不同方法的对比
为了对常规方法和热力学组合方法进行比较，笔者分别使用热力学组合方法 CBS-QB3、G4(MP2) 以及 G4 计算了六个简单气体分子的生成焓。同时也使用了常规方法进行了计算，常规方法在 B3LYP/6-311G** 的级别下做的优化和振动分析（除 C(g)，C(g) 是单原子，因此只在相同级别下做了振动分析）。并在优化后的结构上，使用三种不同的级别分别计算了单点，再通过 Shermo 计算得到热力学量。所有热力学数据都是在 298.15 K、1 个标准压力下计算得到的，基态单质氧和基态气态碳的热力学数据都是在三重态下计算得到的。计算结果如下表所示：
值得一提的是，计算生成焓的策略就和生成焓的定义一致，由于生成焓是最稳定单质生成物质的热量变换，所以直接把最稳定单质和产物的热力学量算出来做差就行了。不过 C 比较特殊，C 的最稳定单质为石墨，Gaussian 里直接计算的是 C(g) 的热力学量，如果想根据定义得到合理的生成焓，就必须加上石墨的升华焓 715 kJ/mol。调用 Shermo 的命令如下所示。
1 | Shermo xxx.out -E xxx -ilowfreq 2 -sclZPE 0.9882 -sclheat 1 -sclS 1 -sclCV 1.0 -T 298.15 -P 1.0 |
在使用常规方法计算时，强烈推荐使用笔者开发的 EasyShermo 搭配 Shermo 使用，笔者要得到这些数据一共得用 12 个振动分析结果和 12 个单点结果，使用 EasyShermo 直接批量调用 Shermo 一秒就能得到热力学量的计算结果，并且输出在了 txt 文件中，方便检查，十分便捷。
参考文章
本文主要参考了以下 Sobereva 的博文：