开源量子化学软件PSI4[9]:用振动频率分析给总能加上零点能校正

往期回顾：

开源量子化学软件PSI4[1]：综述

开源量子化学软件PSI4[2]：Linux终端下使用 Anaconda 安装 PSI4

开源量子化学软件PSI4[3]：用一个例子讲解PSI4的输入文件

开源量子化学软件PSI4[4]：利用DFT计算孤立氢原子（H）的总能

开源量子化学软件PSI4[5]：利用DFT计算孤立质子（H+）的总能

开源量子化学软件PSI4[6]：利用DFT计算孤立氢负离子（H-）的总能

开源量子化学软件PSI4[7]：利用DFT优化氢气分子的结构

开源量子化学软件PSI4[8]：计算水分子分解为氧气和氢气的形成能

写在前面：第[8]篇那个数，为什么和实验对不上

上一篇《开源量子化学软件PSI4[8]：计算水分子分解为氧气和氢气的形成能》里，我们把水、氧气、氢气三个体系的总能代进同一个式子，得到：

-76.4206350591150851 - 0.5*(-150.3340408515116167) - (-1.1736020784246788) = -0.080012554935 (Hartree) = -2.177252535609 (eV)

这个差值是 E(H₂O) − E(H₂) − ½E(O₂)，即水由单质生成时的能量变化；反过来把一个水分子拆开，需要吸收同样大小、符号相反的能量：+0.080013 Hartree = +2.177253 eV（换算成摩尔量是 210.07 kJ/mol）。

但翻开热化学数据表，298.15 K 下气态水分解的标准焓变 ΔH° 是 +241.826 kJ/mol，约 2.506 eV；标准自由能变 ΔG° 是 +228.57 kJ/mol，约 2.369 eV（NIST Chemistry WebBook，CODATA 评审值）。与实验焓变之间差了 31.75 kJ/mol，也就是 0.329 eV。

这不一定是“算错了”，更根本的原因是两者本来就不是同一个物理量。前面几篇计算里从输出文件读出的 Total Energy，是把原子核固定在平衡位置、温度 0 K、完全没有核运动时的纯电子能。要跟量热实验站在同一把尺子上，至少还要补三件事：零点振动能、有限温度的热焓校正、以及熵的贡献。这一讲先讲清楚怎么用 PSI4 的 frequency() 补上第一件事，并把 Hartree、eV、kJ/mol 三种单位的换算一次讲清楚。

一、从电子能到反应焓：四级台阶

图 1　从纯电子能走到反应自由能的四级校正（示意）

第一级 ΔE_elec 就是前面几篇计算里的总能。第二级是零点振动能：即使在 0 K，原子核也不会静止，谐振近似下每个振动模式贡献 ½hν_i，所以 ZPE = ½Σhν_i，这里的 ν_i 正是频率计算给出的振动波数。第三级把温度升到 298.15 K，按刚性转子—谐振子近似补上平动、转动、振动的热校正，得到 ΔH。第四级再减去 TΔS，得到能和平衡常数、电极电势直接对照的 ΔG。

对 H₂O → H₂ + ½O₂ 这类“分子数变多”的反应，平动熵和转动熵的贡献尤其不能忽略：3 个原子组成的非线性水分子有 3 个振动模式，而双原子的 H₂ 和 O₂ 各只有 1 个。反应后平动自由度增加，熵项会把 ΔG 明显往下拉。

二、输入文件 H2O_freq.in

分子坐标直接沿用第[8]篇 PSI4 优化后输出的几何，保证从同一个极小点出发。把下面这段存成 H2O_freq.in：

# 设置计算的分子（几何取自第[8]篇 PSI4 优化后的 H2O） molecule h2o { 0 1 O   -0.000000000000   0.000000000000  -0.067689544787 H    0.000000000000  -0.756603946757   0.537141098344 H   -0.000000000000   0.756603946757   0.537141098344 units angstrom }  # 设置计算方法（与前八篇保持一致） set { basis = cc-pVDZ scf_type = pk }  # 先优化到稳定结构，再做振动频率分析 optimize('b3lyp-d3') frequency('b3lyp-d3')

第 1 行：# 后面是注释，PSI4 不执行。

第 2～8 行：molecule h2o { ... } 定义分子；0 1 表示电荷 0、自旋多重度 1（闭壳层水分子）；三行是 O、H、H 的笛卡尔坐标；units angstrom 指定单位为埃。这几行的写法在第[3]篇已经逐行讲过。

第 11～14 行：set 块沿用系列一贯的设置，basis = cc-pVDZ、scf_type = pk，方便和第[4]～[8]篇的结果直接比较。

第 17 行：optimize('b3lyp-d3') 先做几何优化。这一步不能省：频率分析默认在给定几何处求 Hessian，如果结构没落到极小点，就会算出虚频。

第 18 行：frequency('b3lyp-d3') 是本讲的新东西。它会在同一方法、同一基组下计算 Hessian 矩阵，对角化后给出振动频率，同时把 ZPE 和热力学校正量一并打印出来。注意它比单点能贵得多——能量对坐标的二阶导数要算很多次，耗时通常是单点能的数十倍，建议先用小基组试通再上正式计算。

三、运行，以及输出文件里该看哪几段

psi4 H2O_freq.in H2O_freq.out

输出文件很长，按下面这个顺序找关键段就够了：

① ==> Geometry <==：确认程序读进去的确实是优化后的结构，坐标和输入文件一致。

② Harmonic Vibrational Analysis：频率表。水分子有 3 个振动模式，三个频率都应该是正数；一旦出现负频（虚频），说明这不是极小点，先回头检查优化是否收敛、初始几何是否对称得过头。

③ 文件末尾的热力学校正表，大致长这样（~ 处是你自己运行时才有的数值）：

  ==> Harmonic Vibrational Analysis <==    Irrep   Harmonic Frequency (cm-1)   IR Intensity   -----   -------------------------   ------------    A1            ~                        ~    A1            ~                        ~    B2            ~                        ~    Zero-Point Energy (Hartree)             :  ~   Thermal correction to Energy (Hartree) :  ~   Thermal correction to Enthalpy         :  ~   Thermal correction to Gibbs Free Energy:  ~   Total Energy (Hartree)                 :  ~   ZPE corrected energy (Hartree)         :  ~   Enthalpy (Hartree)                     :  ~   Gibbs Free Energy (Hartree)            :  ~

几个量的含义：Zero-Point Energy 是 ½Σhν_i；Thermal correction to Enthalpy 把 0 K 的电子能推到 298.15 K 的焓；Thermal correction to Gibbs Free Energy 再扣掉熵项；最后三行分别是含 ZPE 校正的能量、总焓与总自由能。

反应能要按物种分别算，再做差：ΔE₀ = ΔE_elec + [ ZPE(H₂) + ½ZPE(O₂) − ZPE(H₂O) ]。所以 H₂O、H₂、O₂ 三个体系各要跑一次 frequency()，把各自的 ZPE 代进方括号，才是 0 K 下含零点能的分解能。

这里刻意不替你贴出具体数值：频率与 ZPE 会随 PSI4 版本、基组、积分网格和收敛限出现末位差异，凭空写一个数反而会误导。请以你自己 H2O_freq.out 的真实输出为准，把三个物种各跑一遍再代入上式——这与本系列前八篇“贴输出文件、逐行读”的做法是一致的。

四、单位换算：把第[8]篇的数字放到同一把尺子上

前八篇里 Hartree、eV、kJ/mol、Å 混着出现，最容易被忽略的就是单位。记住 CODATA 2018 的四个换算常数：

最后一行是读频率表时最常用的：把 cm⁻¹ 的振动波数换成 eV。把第[8]篇实测的 0.080013 Hartree 代进去，就是下面这段可以直接跑的换算：

HA_TO_EV = 27.211386245988   # CODATA 2018 HA_TO_KJ = 2625.4996394799   # Hartree -> kJ/mol dE_ha = 0.080012554935       # 第[8]篇实测形成能的绝对值  print('eV     =', dE_ha * HA_TO_EV) print('kJ/mol =', dE_ha * HA_TO_KJ)

实跑结果：2.177253 eV、210.07 kJ/mol。前一个数与第[8]篇正文的 -2.177252535609 eV 取绝对值后一致（末位差异来自舍入），可以当作一次交叉校验。

图 2　水分解能量：第[8]篇两种实测电子能，与 NIST 气相热化学值在同一标尺下的对比（本次实验实测换算）

四根柱子讲的不是“谁对谁错”，而是三种不同的东西：左边两根是纯电子能（第[8]篇 PSI4 与 VASP 各自的实测输出），右边两根是含核运动与温度效应的实验焓变与自由能变。PSI4 的 2.177 eV 与实验 ΔH° 的 2.506 eV 相差 31.75 kJ/mol，缺口里既有零点能与热校正，也有 cc-pVDZ 基组、B3LYP-D3 泛函带来的方法误差；VASP 的 -2.53122095 eV（分解反应取绝对值即 2.531 eV）看似更贴近实验，但那是另一套赝势与平面波基组下的结果，不能据此断言某一个“更准”。把这些贡献一项项拆开，正是本篇和下一篇要做的事。

五、常见报错与排查

出现虚频：先看 optimize 是否真正收敛；仍不收敛时把坐标轻微扰动后再优化，或降低初始对称性。频率分析必须在极小点上做，否则后面的 ZPE 没有意义。

SCF 不收敛：在 set 块里加 maxiter 200、damping_percentage 20，或改用 guess sad 初始猜测；开壳层体系（如 O₂ 三重态）别忘了 reference = uhf 与正确的自旋多重度。

太慢或内存不足：加 memory 4 GB；用 scf_type df（密度拟合）能明显提速，但它与系列一直沿用的 pk 在数值末位会有细微差别，正式对比时请保持设置一致。

单位看错：频率表头写的是 cm⁻¹，能量表写的是 Hartree，热校正表可能同时出现 kcal/mol 和 kJ/mol。抄数之前先读表头，这是最容易出错、也最不容易被发现的一类问题。

六、本篇小结与下一篇预告

① 前面几篇计算输出里的 Total Energy 是 0 K、核固定条件下的纯电子能；② frequency() 一次给出 Hessian、振动频率、ZPE 与热力学校正，虚频是结构没落到极小点的信号；③ 反应能必须对每个物种各算一次 ZPE 后再相减；④ 记住 1 Hartree = 27.211386245988 eV = 2625.4996394799 kJ/mol，就能把计算值和热化学数据放到同一把尺子上。

下一篇《开源量子化学软件PSI4[10]》打算讲基组收敛性检验：把第[8]篇的三个体系从 cc-pVDZ 依次换到 cc-pVTZ、cc-pVQZ，看那 0.080013 Hartree 的形成能怎么一步步移动，把“方法和基组带来的误差”从上面那个缺口里单独拆出来。

你跑 PSI4 时最常卡在哪一步——安装、SCF 不收敛，还是频率算出的虚频？欢迎在留言区把报错贴出来，我们会挑典型问题放到后续篇里专门讲。