CP2K 分子零点能与气相热修正

电子能量不等于有限温度自由能。得到分子驻点和谐振频率后，可以加入零点能、振动、平动与转动的热贡献。本篇沿用真实H₂O、D₂O和CO₂振动，给出298.15 K、101325 Pa下的理想气体RRHO流程。

下载完整输入、原始输出、频率表和绘图脚本，解压后进入 run.sh 所在目录。本文在鸢算科技P100上用CP2K2026.2双精度 cp2k.psmp 实算；构建采用GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3。分析环境为ASE3.29.0。

1 提交与查看结果#

sbatch submit.slurm water

模板申请四张P100、四个MPI进程，每进程一个CPU核和一个OpenMP线程。NPROC_REP 1让不同进程计算不同位移，脚本按MPI全局编号分配GPU。在已分配的计算资源里也可以执行 bash run.sh water。套餐用户用 yskit 生成对应分区的提交脚本，再保留这条运行命令；不要在登录节点直接运行。

脚本依次执行“优化结构 → 提取最终XYZ → 有限位移力、偶极与极化率 → 振动模式 → 数据提取 → 画图”。inputs为可编辑输入，reference为本文结果，work为本次新结果。每个条件请新解压一份，脚本不会覆盖已有work。

2 明确热力学模型#

H₂O与D₂O是非线性分子，有三个转动自由度；CO₂是线性分子，有两个。三者转动对称数都取2，自旋取0；重水原子质量与振动计算一致。振动频率必须来自相应稳定结构，不能把过渡态的虚频直接塞进稳定分子公式。

CP2K输入开启 THERMOCHEMISTRY TRUE 时会写一份原生热化学表；本篇使用ASE3.29.0显式指定几何与对称数的RRHO结果。两份表采用的转动/对称处理不同，附件保留原生表便于追溯，正文表格统一使用下列ASE计算，不能从两份表中挑选各一列拼起来。

3 提取频率，算H、S、G#

公开的 analyze.py 已完成该流程。以下代码可以在算完water后直接执行，用实际日志和最终结构得到同样的气相修正：

module load ase/3.29.0
python - <<'PY'
from pathlib import Path
import re
from ase.io import read
from ase.units import Hartree, invcm, kJ, mol
from ase.thermochemistry import IdealGasThermo
d = Path('work/cases/water-ir-raman')
text = (d / 'output.log').read_text()
freq = [float(x) for line in re.findall(
    r'VIB\|Frequency \(cm\^-1\)\s*([^\n]+)', text) for x in line.split()]
opt_text = (d / 'opt' / 'output.log').read_text()
energy = float(re.findall(
    r'ENERGY\| Total FORCE_EVAL \( QS \) energy \[hartree\]\s+([-+0-9.Ee]+)', opt_text)[-1]) * Hartree
atoms = read(d / 'optimized.xyz')
atoms.set_masses([15.9994 if s == 'O' else 1.00794
                  for s in atoms.get_chemical_symbols()])
gas = IdealGasThermo(vib_energies=[v * invcm for v in freq],
    geometry='nonlinear', potentialenergy=energy, atoms=atoms,
    symmetrynumber=2, spin=0, vib_selection='exact')
T, pressure = 298.15, 101325.0  # K、Pa，显式传入
H = gas.get_enthalpy(T, verbose=False)
S = gas.get_entropy(T, pressure, verbose=False)
G = gas.get_gibbs_energy(T, pressure, verbose=False)
factor = 1 / (kJ / mol)
print('ZPE kJ/mol:', sum(freq) * invcm / 2 * factor)
print('H-E kJ/mol:', (H-energy) * factor)
print('S J/mol/K:', S * factor * 1000)
print('G-E kJ/mol:', (G-energy) * factor)
PY

CP2K能量先由Hartree转为eV，频率由cm⁻¹转为eV；ASE返回的每分子eV再转成kJ/mol。此处ASE3.29.0的熵接口显式传压力数值 101325.0 Pa，不要再乘一个原子单位压力换算常数。CO₂应使用 geometry='linear' 和四个内部频率，D₂O应采用对应频率与2.01410177812 u的氢质量。

4 本例结果#

H−E电子已包含零点能和温度贡献，G−E电子又包含−TS，不能在G修正之外再加一次ZPE。图表中的量是修正项；不同分子的总能量需要按反应方程式配平之后才能求反应热/自由能。

5 从修正到反应自由能#

对每个独立气相物种分别算电子能和一致定义的修正，再按化学计量数求差：ΔG=Σνᵢ[E电子,ᵢ+(G−E电子)ᵢ]。温度、压力、质量、对称数和计算方法需要一致；改变气相标准态会改变熵与自由能。

这里没有计算具体反应的ΔG，也没有把这些气相平动/转动项加到吸附态。柔性分子低频、受阻转动、溶液标准态、表面吸附和非谐效应需要另外的物理处理。低频模式出现时，先判断模型是否适合RRHO。

提取与重画#

module load ase/3.29.0
python analyze.py work
python plot_results.py work
ls work/analysis
ls work/figures

CSV给出模式编号、频率、IR强度和Raman活性；JSON保留质量、最终结构和热修正。原生 water-VIBRATIONS-1.mol 是Molden振动文件，可以用支持该格式的分子查看器检查模式位移。项目名前缀water是本包统一输出名，不代表CO₂输入成了水。

只重画本文结果，把上述命令中的work换成reference。脚本用实际模式的线位置与活性叠加高斯曲线，展示宽度σ=10 cm⁻¹；这是绘图展宽，不是计算得到的声子寿命。

模型说明见ASE热化学文档与CP2K振动/热化学定义。