ASE分子热化学与自由能修正

下载全部输入、可运行脚本、实际数据和四张高清图。在已分配CPU资源内，使用新解压目录执行 bash run.sh；ASE3.29.0、单进程单线程，生成 work/。本包不调用GPU或外部DFT程序。

1 先确定你要哪一种自由能#

气相分子的自由能含平移、转动和振动；吸附态常用的谐振近似不直接加入气相平移/转动。两者物理模型不同，不能仅因为都由频率算出，就混用修正值。

本例复用集群已完成的GFN2-xTB6.7.1水优化和Hessian，原始结构、JSON与频率表在包内 inputs/，不重新运行xTB：

bash run.sh
cat inputs/vibspectrum
cat work/phonon-output/water-thermo.csv

三个真实振动频率为 1539.649482、3642.449480、3650.577078cm⁻¹。另六个接近零的模对应整体平移/转动的数值残差，不作为分子内部振动加入。

2 用原始数据建立理想气体模型#

import json
import numpy as np
from ase.io import read
from ase.units import invcm, Hartree
from ase.thermochemistry import IdealGasThermo

raw = json.load(open('inputs/xtbout.json'))
frequency = np.array(raw['vibrational frequencies / rcm'])
vib_energy = frequency[frequency > 10] * invcm
water = read('inputs/xtbopt.xyz')
thermo = IdealGasThermo(vib_energies=vib_energy,
                       geometry='nonlinear',
                       potentialenergy=raw['total energy'] * Hartree,
                       atoms=water, symmetrynumber=2, spin=0)
G = thermo.get_gibbs_energy(298.15, 100000, verbose=False)
print(G)

invcm将cm⁻¹转换为eV；原始总能以Hartree计，乘 Hartree后才与振动能量相加。geometry='nonlinear'适用于非线性水分子；对称数2，自旋量子数0。其他分子不能照搬水的对称数、线性/非线性选项和自旋。

压强单位是 Pa。100000Pa为1bar，101325Pa为1atm；把1输入为压强会得到1Pa条件，并非1bar。

3 比较纯振动谐振模型#

from ase.thermochemistry import HarmonicThermo
harmonic = HarmonicThermo(vib_energies=vib_energy,
                         potentialenergy=raw['total energy'] * Hartree)
F = harmonic.get_helmholtz_energy(298.15, verbose=False)
print(F)

HarmonicThermo这里只使用三个振动能量；IdealGasThermo还根据分子结构、温度、压强和对称性加入气体平移/转动项。

左图比较同一套水振动数据的两种模型；右图是独立Cu EMT晶体示例，参考电子能和归一单位不同，不直接比较两幅图的绝对值。

4 自己查看温度和压强影响#

附件输出100、200、298.15、400、600K的H、S、G与修正值：

head work/phonon-output/water-thermo.csv

其中H和G为eV/分子，S为eV/(分子·K)。对理想气体，同温下从1bar升到2bar应有 ΔG=kBT ln2；附件用这个关系检查压强设置与单位。比较反应自由能时，对反应物和产物使用一致方法、标准态、温度及化学计量系数。

5 如何用于自己的体系#

先优化对应电荷与自旋态，确认频率没有真实虚频，再确定哪些自由度进入模型。低频模的谐振近似可能带来很大熵贡献；约束、表面吸附、溶剂及标准态变化需要相应方法，不能随意删频率或把全部模取绝对值。

本例仅排除六个数值接近零的整体模；对真实软模、虚频或受约束结构应逐个判断。GFN2-xTB和PBE等方法的电子能与频率也不应随意混搭。参考：ASE热化学说明。