1.原理
对于极小体系真空下的热力学计算，要考虑的体系热力学能量可以分为电子能量和热力学校正两部分。前者对应势能面极小值点的能量。
电子能量在用Hartree-Fock和post Hartree-Fock方法计算的时候，可以分为自洽场能量即HF能量和电子相关能校正两部分。
电子能量一个维度上受电子相关作用考虑程度的影响，一般来说越高级的电子相关post HF方法此维度上会越精确，但是计算量也会越大。另一个维度上受所用的基函数组的完备性的影响，越完备误差越小，但是也越昂贵。
对于Dunning的相关一致性基组cc-pVxZ，随着zeta数增加基组越来越大，自洽场能量和电子相关校正能量都有很好的收敛性。有人事先拟合了一些公式和参数，按照这些方法可以在较小基组下获得较大基组的近似结果，这称为基组外推。
CCSD(T)方法对动态电子相关考虑已经很好，但是计算量太大。MP2方法考虑一般，但是计算量较小。利用电子相关能近似可加性，可以用MP2下小基组和大基组的电子相关能能量差来估计CCSD(T)下小基组和大基组的电子相关能能量差，这样可以用MP2下大基组的结果和CCSD(T)下小基组的结果(后者计算过程中会顺便得到MP2的结果)来估计CCSD(T)下大基组的电子相关能。
使用2-zeta，3-zeta和4-zeta进行两者的外推是比较合理的。
本文的脚本没有弥散函数，如果对于阴离子体系，可以使用月份基组加上一定程度的弥散函数。
热力学校正部分能量在0K下为零点能，即最低振动能级和势能面极小值点的能量差。此温度下焓和吉布斯自由能的校正都是0。在室温下，总的吉布斯自由能在零点能的基础上额外还要加上吉布斯自由能随着升温的变化。这可以由Gibbs-Helmholtz方程用升温导致的焓的校正和熵的校正计算出来。
由于非谐振模型计算极为耗时，常用的是谐振模型，这会引入一定系统偏差。为此有人提出了频率校正因子，分别对基频、零点能、升温的焓校正以及熵的校正做了拟合，将频率分别乘以对应的校正因子之后再用配分函数等算出的这类热力学校正能量更接近实验值。
所有计算前必须进行几何优化，几何优化和频率计算必须在同一水平。本程序采用B3LYP-D3(BJ)/6-31G*进行优化和频率计算。这个精度对于这两个目的已经足够，而且其频率校正因子十分常见。电子能量部分分别计算CCSD(T)/cc-pVTZ、CCSD(T)/cc-pVDZ以及MP2/cc-pVQZ的能量。注意MP2的能量计算会顺便给出该基组的HF能量，CCSD(T)的能量会顺便给出该基组的MP2能量和HF能量。
本方法最大的不足之一在于对大量涉及氢振动的低频模式，缺乏隧道效应的校正，会有精度的损失。但是由于热力学校正部分能量影响本来占的比例就较小，本文目的也只是学习为主，所以问题不大。
其实Gaussian有自带的更好的高精度热力学组合方法使用。本文只是实现了一个低配的热力学组合方法。
注意化学研究的是变化，只有相对数值有意义。直接用此程序算一个分子的热力学能量什么都不能说明。但是对于两个元素相同的反应，比如丙二烯到丙炔的异构化，用该方法分别计算两个分子的能量，然后相减，用得到的吉布斯自由能差通过van't Hoff等温式换算成标准平衡常数，是有意义的。

2.程序说明。
本程序提供了一个Gaussian的输入文件样板，对Gaussian 16 B.01的Windows版和Linux版测试过。
本文的程序可在下面的链接下载。如果链接挂了请评论通知。
带Windows字样的压缩包是Windows的64位版本，带Linux字样的压缩包是Linux的64位版本。
带source字样的是备份的样板文件、配置文件和C语言的源代码。那个不用管。
解压后用文本编辑器打开template.gjf，将最前面的坐标和电荷、自旋多重度改成你要计算的分子，别的什么都不要动，保存。然后自己调用Gaussian计算这个文件。Linux下输出文件默认是template.log，手动重命名成template.out。Windows下默认输出文件就是template.out。确保输出文件也在程序的文件夹中。
找到template.chk文件。Windows版可能在安装目录的Scratch目录。命令行调用就会生成在本地。也可以在template.gjf计算前编辑里面所有的%chk手动指定路径，但是文件名一定最后要改成template.chk放在程序的文件夹中。
可以手动调用formchk工具把template.chk转化成template.fch，Linux下默认转化出来叫template.fchk，手动重命名成template.fch，把这个文件也放入程序的文件夹中。如果已经有了fch文件，就可以删掉chk文件节省空间。你也可以仅保留chk文件，程序会自动转化。
确定环境变量设置正确，Gaussian的安装目录必须在PATH中，比如C:\G16W或者/home/gaussian/g16这样。PATH变量里面Windows下如何编辑请自行用搜索引擎查看。
GAUSS_EXEDIR最好也设置成Gaussian的安装目录。
编辑配置文件settings.ini，里面的各个参数是计算过程中要用到的。可以手动编辑你的温度，但是明显偏离室温误差会很大。其他参数一般不用动，除非文献有更新。如果找不到settings.ini，程序会用内置的默认参数，和settings.ini一开始的一样。
运行HighPreciseThermalEnergy的可执行文件，Linux下就叫这个，Windows下叫HighPreciseThermalEnergy.exe。WIndows下双击即可。等几秒钟，程序就会输出一些必要的文件和所有包含所有结果的Results.txt文件。这个文件包括提取和计算出的所有有用信息以及公式。最后一行的数据就是最终总的热力学能量。所有能量单位都是Hartree也就算a.u.，换算成J/mol要乘上换算因子((m_e)*e^4)/(4*(epsilon_0)^2*h^2)*(N_A)，其中m_e是电子质量9.10956E-31 kg，e是电子电荷-1.602176634E-19 C，epsilon_0是真空介电常数8.854187817E-12 F/m，h是普朗克常数6.62607015E-34 J·s，N_A是阿伏伽德罗常数6.02214076E23 /mol。换算成kJ./mol要再除以1000 J/kJ，换算成kcal/mol要再除以指定的水的比热容4.184 J/(kg·K)。

http://tieba.baidu.com/p/6471669758?pid=129707226529&cid=0#0