2. 频率分析 (3)：热力学能矫正#
创建时间：2021-06-18
在这份文档中，我们将简单地讨论从 Gaussian 生成的 formated checkpoint 文件 (fchk 或 fch 后缀名)，产生热力学能量矫正。
我们计算的分子与 频率分析 (1) 相同，为没有优化到能量最低结构的 C2O4H+ 离子。
该文档的重要的参考资料是 Gaussian 的白皮书 Thermochemistry in Gaussian [1]。我们所使用的公式也 (很无耻地) 会与该文档几乎完全相同；并且作为程序实现笔记，也不会对具体的公式推导作讨论。
警告
不处于能量最低结构的分子一般来说不适合用作频率分析。此时不仅分子光谱、同时热力学矫正从理论上也是不允许的。
这份文档尽管使用了有虚频的分子，但若要进行有价值的热力学矫正，仍然需要先对分子的结构进行优化。
分子对应的输入卡 C2O4H.gjf、输出文件 C2O4H.out 与 fchk 文件 C2O4H.fchk 在链接中可供下载。这份文档的目标将是重复输出文件中的热力学矫正量。其一部分输出是：
with open("C2O4H.out", "r") as f:
while "Zero-point correction" not in f.readline(): continue
for _ in range(23): print(f.readline()[:-1])
Thermal correction to Energy= 0.030349
Thermal correction to Enthalpy= 0.031293
Thermal correction to Gibbs Free Energy= -0.002096
Sum of electronic and zero-point Energies= -370.093195
Sum of electronic and thermal Energies= -370.088816
Sum of electronic and thermal Enthalpies= -370.087872
Sum of electronic and thermal Free Energies= -370.121261
E (Thermal) CV S
KCal/Mol Cal/Mol-Kelvin Cal/Mol-Kelvin
Total 19.044 14.495 70.273
Electronic 0.000 0.000 0.000
Translational 0.889 2.981 39.370
Rotational 0.889 2.981 26.171
Vibrational 17.267 8.534 4.732
Vibration 1 0.683 1.701 1.500
Vibration 2 0.731 1.565 1.145
Vibration 3 0.897 1.158 0.562
Vibration 4 0.941 1.067 0.479
Q Log10(Q) Ln(Q)
Total Bot 0.920866D+01 0.964196 2.220144
Total V=0 0.811762D+13 12.909429 29.725059
Vib (Bot) 0.238627D-11 -11.622280 -26.761289
我们的文档的内容原则上可以重现所有的上述数据。
备注
热力学能量矫正本质上是统计热力学的初步应用。我们这里遇到的分子是最为常见的没有对称性、单重态的分子。但一些特殊的情况，譬如具有对称性、多重态等情况，需要对下述代码作改动。这些改动我们不在文档中作补充，因此读者还是需要参考 Gaussian 白皮书 [1]。
2.1. 准备工作#
频率分析部分是由 freqanal.py 完成的；它还需要调用读取 Gaussian formchk 文件的小程序 formchk_interface.py。这些程序可以下载。频率分析的具体做法已经在 频率分析 (1) 有所陈述；对于比较一般的非线性分子，它应当可以输出与 Gaussian 近乎一致的频率结果。
from freqanal import FreqAnal
import numpy as np
np.set_printoptions(5, linewidth=150, suppress=False)
为了能进行单位换算，我们还需要定义一些常数。这里对它们在这份文档中的变量名与符号作说明，并给出大致的数量。
E_h\(E_\mathrm{H}\)：Hartree 能量 \(4.360 \times 10^{-18} \ \mathrm{J}\)；a_0\(a_0\)：Bohr 半径 \(5.292 \times 10^{-11} \ \mathrm{m}\)N_A\(N_\mathrm{A}\)：Avogadro 常数 \(6.022 \times 10^{23} \ \mathrm{mol}^{-1}\)c_0\(c\)：真空光速 \(2.998 \times 10^{8} \ \mathrm{m} \ \mathrm{s}^{-1}\)k_B\(k_\mathrm{B}\)：Boltzmann 常数 \(1.381 \times 10^{-23} \ \mathrm{J} \ \mathrm{K}^{-1}\)R\(R\)：Mole 气体常数 \(8.314 \ \mathrm{J} \ \mathrm{mol}^{-1} \ \mathrm{K}^{-1}\)h\(h\)：Planck 常数 \(6.626 \times 10^{-34} \ \times{J} \ \mathrm{s}\)P_0\(P_0\)：标准大气压 \(101325. \mathrm{kg} \ \mathrm{m}^{-1} \mathrm{s}^{-2}\)amu\(m_\mathrm{u}\)：原子质量单位 \(1.661 \times 10^{-27} \ \mathrm{kg}\)
# https://docs.scipy.org/doc/scipy/reference/constants.html
from scipy import constants
from scipy.constants import physical_constants
E_h = physical_constants["Hartree energy"][0]
a_0 = physical_constants["Bohr radius"][0]
N_A = physical_constants["Avogadro constant"][0]
c_0 = physical_constants["speed of light in vacuum"][0]
k_B = physical_constants["Boltzmann constant"][0]
R = physical_constants["molar gas constant"][0]
h = physical_constants["Planck constant"][0]
P_0 = physical_constants["standard atmosphere"][0]
amu = physical_constants["atomic mass constant"][0]
Calorie 与 Joule 的换算比例 4.184
pi\(\pi\)：圆周率
calorie = 4.184
pi = np.pi
我们在整个文档中，使用的温度是 \(T = 298.15 \ \mathrm{K}\)
T = 298.15
对 C2O4H+ 离子的频率分析对象储存在变量 fa 中。
fa = FreqAnal("C2O4H.fchk")
整个热力学分析过程需要分为平动 (translation)、电子态 (electronic)、转动 (rotation)、振动 (vibrational) 四部分。需要计算的基础热力学矫正量是熵 \(S\) (entropy)、内能 \(E\) (thermo energy)、热容 \(C\) (heat capacity)。这里的热容是指恒容热容。在计算熵时，我们还需要给出配分函数 \(q\) (partition function)。
2.2. 平动 (Translation)#
配分函数
\[ q_\mathrm{t} = \left( \frac{2 \pi m k_\mathrm{B} T}{h^2} \right)^{3/2} \frac{k_\mathrm{B} T}{P_0} \]m\(m\) 指分子质量，Gaussian 的输出是原子质量单位。它需要先转为 SI 单位制再进行计算。
m = fa.mol_weight.sum() * amu
q_t = (2 * pi * m * k_B * T / h**2)**(3/2) * (k_B * T / P_0)
q_t
32994942.55727539
熵 (\(\mathrm{cal} \ \mathrm{mol}^{-1} \ \mathrm{K}^{-1}\))
\[ S_\mathrm{t} = R (\ln q_\mathrm{t} + 1 + 3/2) \]
S_t = R * (np.log(q_t) + 1 + 3/2) / calorie
S_t
39.37022220447884
内能 (\(\mathrm{kcal} \ \mathrm{mol}^{-1}\))
\[ E_\mathrm{t} = \frac{3}{2} R T \]
E_t = 3/2 * R * T / 1000 / calorie
E_t
0.8887274245542662
热容 (\(\mathrm{cal} \ \mathrm{mol}^{-1} \ \mathrm{K}^{-1}\))
\[ C_\mathrm{t} = \frac{3}{2} R \]
C_t = 3/2 * R / calorie
C_t
2.9808063879063096
2.3. 电子态 (Electronic)#
警告
电子态的配分函数 \(q_\mathrm{e}\) 受分子多重度的影响；这也会同时影响到熵矫正 \(S_\mathrm{e}\)。因此对于多重度不为 1 的分子，下述代码将需要作修改。具体地来说，需要将 ω_0 \(\omega_0\) 改为多重度的数值。
配分函数
\[ q_\mathrm{e} = \omega_0 \]ω_0\(\omega_0\) 指分子多重度。
ω_0 = 1
q_e = ω_0
q_e
1
熵 (\(\mathrm{cal} \ \mathrm{mol}^{-1} \ \mathrm{K}^{-1}\))
\[ S_\mathrm{e} = R \ln q_\mathrm{e} \]
S_e = R * np.log(q_e) / calorie
S_e
0.0
内能 \(E_\mathrm{e}\) 与热容 \(C_\mathrm{e}\) 取零值。
E_e = C_e = 0
2.4. 转动 (Rotation)#
警告
分子转动的配分函数、熵、内能与热容计算都受分子本身构型而影响。
我们这里计算的是不具有对称性的分子；
对于线性分子，所有物理量的计算都将发生变化，这里的代码将不能使用；
对于具有对称性的分子需要修改代码；具体地来说，需要将
σ_r\(\sigma_\mathrm{r}\) 设置为分子的对称数；对称数应当指分子的一些原子置换，但通过一系列旋转或反映操作可以重现原来分子的总数量，对于水为 2、氨气为 3、甲烷或苯为 6。
q_r = 1.16957e5
S_r = R * (np.log(q_r) + 3/2) / calorie
S_r
26.17060894487201
转动惯量
Ixyz\(I_x, I_y, I_z\) (\(\mathrm{kg} \ \mathrm{m}^2\))。需要注意，这里并非真的是绕 \(x, y, z\) 轴转动的惯量，即\[ I_x \neq \sum_{A} m_\mathrm{A} r_{Ax}^2 \]它需要通过对 \(3 \times 3\) 的转动惯量矩阵 \(I_{xx}, I_{xy}, \cdots, I_{zz}\) 作对角化得到。对角化同时会得到三个转动主轴；由于这三个转动主轴相互垂直，确实地可以构建坐标系，因此才称为 \(I_x, I_y, I_z\)，但这里的 \(x, y, z\) 相对于输入卡的坐标系一般有旋转。
Ixyz = fa.rot_eig * amu * a_0**2
Ixyz
array([1.33218e-45, 2.34564e-45, 3.43473e-45])
转动特性温度
Θxyz_r\(\Theta_{\mathrm{r}, x}, \Theta_{\mathrm{r}, y}, \Theta_{\mathrm{r}, z}\) (\(\mathrm{K}\))\[ \Theta_{r, x} = \frac{h^2}{8 \pi^2 I_x k_\mathrm{B}} \]
Θxyz_r = h**2 / (8 * pi**2 * Ixyz * k_B)
Θxyz_r
array([0.30233, 0.1717 , 0.11726])
配分函数
\[ q_\mathrm{r} = \frac{\pi^{1/2}}{\sigma_r} \left( \frac{T^3}{\Theta_{\mathrm{r}, x} \Theta_{\mathrm{r}, y} \Theta_{\mathrm{r}, z}} \right)^{1/2} \]其中 \(\sigma_\mathrm{r}\) 为分子的对称数。
σ_r = 1
q_r = (pi**(1/2) / σ_r) * (T**3 / Θxyz_r.prod())**(1/2)
q_r
116957.22677656508
熵 (\(\mathrm{cal} \ \mathrm{mol}^{-1} \ \mathrm{K}^{-1}\))
\[ S_\mathrm{r} = R (\ln q_\mathrm{r} + 3/2) \]
S_r = R * (np.log(q_r) + 3/2) / calorie
S_r
26.170612798005376
内能 (\(\mathrm{kcal} \ \mathrm{mol}^{-1}\))
\[ E_\mathrm{r} = \frac{3}{2} R T \]
E_r = 3/2 * R * 298.15 / 1000 / calorie
E_r
0.8887274245542662
热容 (\(\mathrm{cal} \ \mathrm{mol}^{-1} \ \mathrm{K}^{-1}\))
\[ C_\mathrm{r} = \frac{3}{2} R \]
C_r = 3/2 * R / calorie
C_r
2.9808063879063096
2.5. 振动 (Vibrational)#
备注
热力学问题的一般都要求假设分子处于稳定构型 (因此才能有热力学稳定的状态，各种热力学能量作为温度的状态函数才能成立)。因此，具有振动虚频的分子一般认为是不能进行振动分析的。
我们这里的计算单纯地是 Gaussian 的结果作比对。Gaussian 在计算热力学能时，对虚频的处理是直接忽视。
振动特征温度 (\(\mathrm{K}\))
\[ \Theta_{\mathrm{v}, K} = h \nu_K = \varpi_K h c \]其中，\(K\) 是指震动模式，\(\varpi_K\) 是以长度为量纲的振动频率。
ΘK_v = fa.freq[fa.freq > 0] / 1e-2 * h * c_0 / k_B
配分函数
\[ q_\mathrm{v} = \prod_K \frac{e^{- \Theta_{\mathrm{v}, K} / 2 T}}{1 - e^{- \Theta_{\mathrm{v}, K} / T}} \]这里采用的是与 Gaussian 的最终热力学矫正一致的输出 (即
BOT结果)。
qK_v = np.exp(-ΘK_v / (2 * T)) / (1 - np.exp(-ΘK_v / T))
q_v = qK_v.prod()
q_v
2.386188575409033e-12
熵 (\(\mathrm{cal} \ \mathrm{mol}^{-1} \ \mathrm{K}^{-1}\))
\[ S_\mathrm{v} = R \sum_K \left( \frac{\Theta_{\mathrm{v}, K} / T}{e^{\Theta_{\mathrm{v}, K} / T} - 1} - \ln (1 - e^{- \Theta_{\mathrm{v}, K} / T}) \right) \]
SK_v = R * (ΘK_v / T / (np.exp(ΘK_v / T) - 1) - np.log(1 - np.exp(-ΘK_v / T))) / calorie
S_v = SK_v.sum()
S_v
4.732478685230292
内能 (\(\mathrm{kcal} \ \mathrm{mol}^{-1}\))
\[ E_\mathrm{v} = R \sum_K \Theta_{\mathrm{v}, K} \left( \frac{1}{2} + \frac{1}{e^{\Theta_{\mathrm{v}, K} / T} - 1} \right) \]
EK_v = R * (ΘK_v * (1/2 + 1 / (np.exp(ΘK_v / T) - 1))) / 1000 / calorie
E_v = EK_v.sum()
E_v
17.266670082670938
热容 (\(\mathrm{cal} \ \mathrm{mol}^{-1} \ \mathrm{K}^{-1}\))
\[ C_\mathrm{v} = R \sum_K e^{- \Theta_{\mathrm{v}, K} / T} \left( \frac{\Theta_{\mathrm{v}, K} / T}{e^{- \Theta_{\mathrm{v}, K} / T} - 1} \right)^2 \]这里 Gaussian 白皮书的记号有略微的错误。
CK_v = R * (np.exp(-ΘK_v / T) * (ΘK_v / T / (np.exp(-ΘK_v / T) - 1))**2) / calorie
C_v = CK_v.sum()
C_v
8.533482711830887
2.6. 总热力学矫正量#
最终的热力学矫正量直接通过将四部分 (平动、电子态、转动、振动) 相加即可。对于配分函数，其结果通过相乘得到。
E_thermal = E_t + E_e + E_r + E_v
E_thermal
19.04412493177947
C_thermal = C_t + C_e + C_r + C_v
C_thermal
14.495095487643507
S_thermal = S_t + S_e + S_r + S_v
S_thermal
70.2733136877145
q_total = q_t * q_e * q_r * q_v
q_total
9.208294504188077
2.7. 最终热力学矫正#
化学中关心的通常是零点能 (zero point energy, ZPE)、焓 (Enthalpy)、Gibbs 自由能 (Gibbs free energy)。这些量不难通过上面的结果给出。
零点能 (zero point energy, ZPE) (\(E_\mathrm{H}\), Hartree)
\[ E_\mathrm{ZPE} = \sum_{K} \frac{1}{2} \Theta_{\mathrm{v}, K} k_\mathrm{B} \]
corr_zpe = ΘK_v.sum() / 2 * k_B / E_h
corr_zpe
0.025969755181835332
内能矫正 (\(E_\mathrm{H}\), Hartree) 实质上就是刚才结果的单位变换而已。
corr_thermal = E_thermal / (E_h * N_A / 1000 / calorie)
corr_thermal
0.03034874486989149
焓 (Enthalpy) 矫正 (\(E_\mathrm{H}\), Hartree)
\[ H_\mathrm{corr} = E_\mathrm{corr} + k_\mathrm{B} T \]
corr_enthalpy = corr_thermal + k_B * T / E_h
corr_enthalpy
0.03129292973753578
Gibbs 自由能 (Gibbs free energy) 矫正 (\(E_\mathrm{H}\), Hartree)
\[ G_\mathrm{corr} = H_\mathrm{corr} - T S_\mathrm{tot} \]
corr_gibbs = corr_enthalpy - T * S_thermal / (E_h * N_A / calorie)
corr_gibbs
-0.002096189219235073