做催化的课题组对自由能图（Free Energy Diagram）再熟悉不过：横轴是反应坐标，纵轴是吉布斯自由能，每一个台阶代表一个中间态，每一个峰代表一个过渡态。图的最后一行给出决速步的自由能垒——数值超过0.75 eV，这条反应路径在室温下基本走不通。
问题在于：VASP算出来的直接输出是电子能量E_elec，不是自由能G。从E_elec到G，中间隔着振动零点能、有限温度的热焓修正、以及构型熵。每一步修正的物理意义不同，忽略任何一项都可能让自由能的排序出错。
在DFT框架下，一个吸附态或中间态的吉布斯自由能写作：
G = E_elec + ZPE + ∫C_p dT − TS
四项的物理来源分别是：
E_elec（电子能量）。 VASP SCF收敛后的总能量。这是T=0 K下不考虑核运动的静态能量。对于反应自由能来说，不同中间态的E_elec差值通常贡献了G差值的主体（比如80-90%），但剩下的10-20%来自振动修正——在室温附近，这个量级（0.1-0.3 eV）足以颠覆反应的定性结论。
ZPE（零点能）。 即使在0 K，原子也在平衡位置附近做量子力学的零点振动。这个振动能量不能从经典统计力学得出——它是量子谐振子的基态能量：ZPE = Σ(½hν_i)，对所有振动模式i求和。对于吸附在金属表面的小分子（CO、H₂O、NH₃），ZPE通常在0.1-0.5 eV量级。注意：ZPE对于反应物和生成物不是对称的——反应前后化学键数目和类型改变了，振动频率分布就变了，零点能差值（ΔZPE）必须单独计算。
∫C_p dT（焓温修正）。 从0 K升温到有限温度T，体系的焓增加了平移、转动和振动的热能贡献。对于周期性DFT计算，表面吸附物种的平移和转动自由度被冻结（或转化为受阻平动/转动的低频振动模），因此焓修正主要来自振动部分：H_vib(T) = Σ[hν_i / (e^(hν_i/kT) – 1)]。
−TS（熵贡献）。 吉布斯自由能和亥姆霍兹自由能的差值就这一项。分子的气相熵可以很大——室温下气相H₂O的平动熵约140 J/(mol·K)，对应TΔS约0.52 eV。但如果H₂O吸附在表面上，平动自由度变成了表面上的低频振动，熵贡献大幅降低（通常<0.1 eV）。气相到吸附相的熵损失（ΔS_ads）是很多反应自由能计算中最容易被低估的修正项。
计算ZPE和振动热力学修正的前提是拿到振动频率。VASP提供了两种频率计算方法：
IBRION=5（数值差分法）。 VASP对每个被标记的原子施加±位移（POTIM控制位移量，默认0.015 Å），用有限差分算Hessian矩阵。每个位移方向跑一次SCF——所以N个原子需要6N+1次SCF计算（±x、±y、±z各一次+参考点）。这是最可靠的方法，但成本高昂：一个100原子的slab模型，用数值差分做频率计算需要约600次SCF计算，在64核上可能要跑一周。
策略性地缩减计算量：不要对所有原子做频率计算。只对被吸附分子和直接键合的表面原子做频率分析，其余基底原子固定（Selective Dynamics + .FALSE.）。这个近似的物理基础是：基底原子的振动与表面吸附物种的振动频率差距很大（基底<300 cm⁻¹，吸附分子的伸缩/弯曲模在500-4000 cm⁻¹），解耦处理对ZPE的误差通常<0.05 eV。
IBRION=6（密度泛函微扰论，DFPT）。 DFPT用线性响应理论直接计算Hessian矩阵，不需要构造超胞和位移原子。它的计算量和SCF收敛一次的量级相同，比数值差分快一个数量级以上。但DFPT对两个东西敏感：一是体系必须是绝缘体（金属体系的DFPT会因费米面附近的占据数不连续而失败），二是k点密度必须足够（否则低频声学模的频率会虚）。
选择路线：绝缘体+大体系→DFPT优先；金属表面+小分子吸附→数值差分优先。如果DFPT给出的最低频率是虚数（负的平方值），说明初始结构不够稳定——需要重新优化到更严格的力收敛标准（EDIFFG ≤ -0.01 eV/Å）。
拿到频率后，VASP的OUTCAR文件末尾会自动输出ZPE和有限温度的热力学修正（在T=0到你设定的温度范围内）。你也可以用VASPKIT或phonopy做后处理，提取每个振动模的贡献并绘制声子态密度。
以HER（析氢反应）在Pt(111)上的自由能计算为例，三步走：
第一步，计算每个中间态的电子能量。 H⁺ + e⁻ → H*（吸附氢中间体）。气相H₂的能量可以直接用VASP算（一个H₂分子放在15×15×15 Å的格子里），吸附态H在Pt(111) slab上计算。电子能量差ΔE = E(H) − E(Pt slab) − ½E(H₂)。
第二步，对每个态做频率计算。 气相H₂的振动频率约4400 cm⁻¹（实验值），ZPE约0.28 eV。吸附态H在Pt表面有一个垂直于表面的Pt-H伸缩振动（约2000 cm⁻¹），ZPE约0.12 eV。零点能差值ΔZPE = ZPE(H) − ½ZPE(H₂) ≈ −0.02 eV——很小，但不能假设它总是可忽略。
第三步，组装自由能。 ΔG(H*) = ΔE + ΔZPE + Δ∫C_p dT − TΔS。气相H₂在300 K的标准熵是130.7 J/(mol·K)，对应TΔS约0.41 eV。H*是吸附态，熵小得多，通常用½TΔS(H₂)估算——它并不精确，但对于催化剂活性趋势的筛选（火山图的横轴）已经够用。
理想HER催化剂的ΔG(H*)应该接近0 eV——吸附不太强（否则H堵住活性位不放H₂），也不太弱（否则H⁺不愿上表面）。Pt(111)的ΔG(H)实验值约−0.09 eV，DFT计算结果通常在−0.05到−0.15 eV之间，取决于泛函选择和slab厚度。这个精度——约0.1 eV的误差——在室温下对应的反应速率差异约为一个数量级（exp(−ΔG/kT)），足以区分活性和非活性催化剂。
CI-NEB（Climbing Image Nudged Elastic Band）给出的是势能面上的最小能量路径（MEP）——它的能垒是势能垒Eₐ，不是自由能垒ΔG‡。
从Eₐ到ΔG‡的修正逻辑和中间态一样：在过渡态（TS）和初始态（IS）分别做频率计算，各算各的ZPE和热力学修正。一个确保质量的操作：过渡态的频率计算必须给出有且仅有一个虚频——如果虚频率不止一个，NEB找到的不是真正的鞍点；如果没有虚频，那这个态不是过渡态而是中间态。
反应速率常数k的计算最终由过渡态理论（TST）给出：k = (k_BT/h) × exp(−ΔG‡/RT)。这个公式里的ΔG‡不是势能垒Eₐ，而是包含了振动修正的自由能垒——对于室温下的表面反应，ΔZPE的修正量级在0.05-0.2 eV，对应的速率修正倍数是3-50倍。
很多DFT反应机理的文章在自由能数据上经不住推敲——要么是频算的原子范围选得太小（ZPE被低估），要么是气相和吸附相的熵差靠猜的（TΔS被高估或低估）。看一眼Supplementary Information里的振动频率数据和熵估算方法，基本就能判断结果的可信度。
科研学术网（https://www.keyanxueshu.com）参数库、过渡态搜索策略和自由能计算的标准操作流程。
DFT计算吉布斯自由能从来不是算一个数填到文章里就完了。它是一张图——每一个自由能台阶上都写着”这个吸附构型稳还是不稳”，每一个自由能峰上都标着”这条反应路径走得通还是走不通”。做对了，它就是比实验更快的高通量筛选工具。
