pymatgen：分子RRHO与准RRHO自由能
从实际水和硫化氢的结构、电子能与频率出发，计算气相自由能并比较温度、压力和转动对称数。
阅读本篇目录
分子优化得到的电子能，不直接等于某个温度和压力下的气相Gibbs自由能。还需要平动、转动、振动与电子简并度的热力学校正。本例用pymatgen处理已经实际算出的水和硫化氢数据，不重新计算电子结构。
unzip pymatgen-qrrho-wave710-*.zip
cd pymatgen-qrrho
bash results.sh
bash plot.sh
1. 输入来自哪里#
两个孤立分子来自CREST频率实操的实际计算：CREST3.0.2内置tblite/GFN2，先优化，再计算Hessian。这里沿用原有结果，不把GFN2数据说成高层级从头算结果。
| 分子 | 电子能/Hartree | 三个正振动频率/cm⁻¹ |
|---|---|---|
| H₂O | −5.070544450296 | 1539.58，3642.91，3651.05 |
| H₂S | −4.256278831653 | 1157.64，2475.73，2480.68 |
两者均为非线性三原子分子，内部振动数为3N−6=3。原Hessian中的六个刚体平动/转动零模不作为振动输入；不能把真实的虚频取绝对值后当成稳定分子。
cat inputs/water.xyz
cat inputs/hydrogen-sulfide.xyz
cat inputs/molecules.json
2. 加载环境和随包入口#
后处理环境为pymatgen2026.9.24、Python3.11、NumPy/SciPy和Matplotlib。本例compat.py在当前Python进程中统一准RRHO的频率单位：有效转动惯量用物理频率，阻尼阈值与输入波数保持一致；普通RRHO部分保持原接口的计算方式。公共安装不修改。
module load pymatgen/2026.9.24
bash run.sh
脚本检查对应版本源码后使用本地入口。换pymatgen版本应先核对接口和单位，不能删掉版本检查直接运行。
3. 单个分子的完整调用#
以下代码在资料包根目录运行：
import json
from pymatgen.core import Molecule
from compat import QuasiRRHO
r = json.load(open('inputs/molecules.json'))['water']
molecule = Molecule(r['symbols'], r['positions_A'])
q = QuasiRRHO(
mol=molecule,
frequencies=r['vibrational_frequencies_cm_1'],
energy=r['electronic_Eh'],
mult=1, sigma_r=2,
temp=298.15, press=101325.0, v0=100.0)
print('G / Hartree:', q.free_energy_quasiRRHO)
print('H correction / Hartree:', q.h_corrected)
print('S / J mol-1 K-1:', q.entropy_quasiRRHO * 4.184)
energy单位Hartree，频率cm⁻¹，坐标Å，温度K，压力Pa。mult=1表示单重态，sigma_r=2是这两个C₂v分子的转动对称数。v0=100是准RRHO阻尼波数，不是温度。
这个安装版的熵属性数值按**cal/(mol·K)**给出，导出J/(mol·K)时乘4.184；不要把属性说明里的原子单位直接套在这个数值上。能量属性仍按Hartree使用。
4. 实际结果#
在1atm、对称数2时：
| 分子 | 温度/K | S/J·mol⁻¹·K⁻¹ | G−E电子/kJ·mol⁻¹ |
|---|---|---|---|
| H₂O | 298.15 | 188.504623 | 6.560416 |
| H₂O | 400 | 198.451596 | −13.168256 |
| H₂S | 298.15 | 205.376018 | −14.694406 |
| H₂S | 400 | 215.641269 | −36.156566 |
这些是各分子相对自身电子能的校正，不用不同组成分子的总G高低直接比较稳定性。反应自由能需要按配平方程组合相同方法、相同标准态的各物种。
cat analysis/Molecular-quasiRRHO-thermochemistry.csv
本例内部频率均较高，所以准RRHO与普通谐振子结果很接近。对有真实低频扭转和软模式的柔性分子，应更认真检查低频处理、构象采样和阻尼参数；这两个小分子不代表柔性体系的全部误差。
5. 压力和对称数为什么要写清楚#
完整脚本实际比较了298.15/400K、101325/100000Pa、转动对称数1/2，共16组后处理条件。
对于同一理想气体分子，压力从P₁改为P₂，自由能变化应满足RT ln(P₂/P₁)。1atm和1bar分别是101325Pa与100000Pa，不能写了1bar却仍沿用1atm数值。
对称数从1改为2，转动熵减少R ln2，自由能增加RT ln2。对称数不是原子数，也不是电子自旋多重度；复杂分子需要根据实际构型判断。
气相1atm/1bar结果也不能直接当成溶液1mol/L标准态自由能。若研究溶液反应，应一致处理溶剂化、标准态和构象等贡献。
6. 换自己的分子#
准备同一结构、同一方法的优化几何、电子能和频率；确认目标结构的虚频情况，再设置电荷/自旋多重度、转动对称数及标准态。不要把一个构象的电子能与另一个构象的频率拼起来。
ls inputs
bash run.sh
run.sh使用新的解压目录，结果写到analysis。只读参考结果和重画已有图用results.sh、plot.sh即可。
原始接口见pymatgen准RRHO说明，频率与有效转动惯量约定可参照ASE的RRHO实现及其引用的Grimme方法。