VASP计算吉布斯自由能：从声子谱到化学反应方向判定

吉布斯自由能G = H – TS是判断化学反应方向和相稳定性的热力学核心量。在VASP框架下，G的计算需要将静态DFT能量与有限温度下的振动熵、电子熵等贡献组合，本质上是多步计算的系统集成。

吉布斯自由能的DFT分解

对于固体体系，吉布斯自由能可以写为：

G(T,P) = E_DFT + F_vib(T) + F_el(T) + PV

在常压条件下（P=1 atm ≈ 0.0001 GPa），PV项对固体几乎可以忽略（典型值<0.001 eV/atom），因此实际计算中：

G(T) ≈ E_DFT + F_vib(T) + F_el(T)

三个分量的计算方式和适用场景各不相同：

E_DFT：0K静态总能量，来自VASP结构优化后的单点计算。

F_vib(T)：振动自由能，来自声子谱计算。这是有限温度修正的主要贡献者。在简谐近似下，F_vib = k_BTΣ ln[2sinh(ℏω_i/2k_BT)]，需要完整的声子频率谱。实际操作中用Phonopy计算，输入VASP的力常数。

F_el(T)：电子自由能，来自电子能带结构的有限温度占据。对于绝缘体/半导体可忽略；对于金属，尤其在高温下不可省略。VASP在ISMEAR>0时自动计算电子熵，OUTCAR中”free energy TS(EENTROP)”一项就是-TS_el。要获得F_el(T)的温度依赖，需要在不同温度下做SCF计算（TEBEG参数），或后处理DOS积分。

完整计算流程

以相图计算为例，判断TiO₂在室温到2000K范围内金红石相与锐钛矿相的稳定性：

第一步：高精度结构优化

两相分别优化，参数必须完全一致：

PREC = Accurate
EDIFF = 1E-8
ENCUT = 600
ISIF = 3
IBRION = 2
NSW = 100
ISMEAR = 0
SIGMA = 0.01
ADDGRID = .TRUE.
LREAL = .FALSE.

EDIFF=1E-8至关重要——声子频率对力常数极度敏感，能量收敛不够直接产生虚频。ENCUT取600 eV（Ti的ENMAX约222 eV，O约400 eV，取最大值的1.3倍以上）确保Pulay应力可控。

第二步：声子谱计算

用Phonopy生成位移超胞，对每个超胞做VASP单点力计算。金红石TiO₂（6原子原胞）用2×2×2超胞（48原子），锐钛矿TiO₂（12原子原胞）用2×2×2超胞（96原子）。

位移超胞的INCAR：

IBRION = -1
NSW = 0
EDIFF = 1E-8
PREC = Accurate
ENCUT = 600
ISMEAR = 0
SIGMA = 0.01
LREAL = .FALSE.

第三步：F_vib(T)温度扫描

Phonopy的热力学属性功能直接输出F_vib(T)曲线，温度步长通常取10K：

phonopy-qha -e energy.dat -v volume.dat \
  --tmax=2000 --tstep=10 \
  -p -s > thermo_properties.yaml

输出文件中free_energy一列就是F_vib(T)。

第四步：电子自由能修正

TiO₂是半导体（带隙~3 eV），F_el在2000K以下可忽略。如果计算的是金属Ti，则需要用F_el(T) = -T*S_el，其中S_el通过DOS积分获得：

# VASP计算DOS
# INCAR: LORBIT=11, ICHARG=11, NEDOS=2001
# 后处理：用vaspkit或自写脚本积分Fermi-Dirac分布

第五步：组装G(T)并比较

G_rutile(T) = E_DFT(rutile) + F_vib(rutile, T) G_anatase(T) = E_DFT(anatase) + F_vib(anatase, T) ΔG(T) = G_rutile – G_anatase

当ΔG(T) < 0时金红石相稳定，ΔG(T) > 0时锐钛矿相稳定。

计算结果与实验对照

对于TiO₂体系，DFT-PBE计算给出：

• 0K下金红石比锐钛矿高约0.05 eV/f.u.（实验值0.04-0.06 eV）
• 约在900K处ΔG穿越零点，金红石开始变得稳定（实验相变温度~1000K）
• 误差来源：PBE系统性低估带隙（计算1.8 eV vs 实验3.0 eV），导致F_el偏大

加Hubbard U修正（U_eff=4.2 eV for Ti 3d）后，带隙修正到~2.7 eV，相变温度改善到950K，精度提升显著。

非谐效应与高温修正

简谐近似（HA）在高温下的失效是吉布斯自由能计算的主要误差来源。判断何时需要非谐修正的经验法则：

• T < 0.3*Θ_D（德拜温度）：HA足够精确，误差<0.5%
• 0.3*Θ_D < T < Θ_D：HA误差1-3%，QHA可以修正
• T > Θ_D：非谐效应显著，需要三阶力常数

对于声子谱计算中已经出现虚频的高温结构，QHA不再适用，此时应考虑分子动力学方法（AIMD）直接计算自由能，例如用 Thermodynamic Integration 或 Two-Phase Thermodynamic (2PT) 方法。

气相反应中的吉布斯自由能

涉及气体的反应（如催化、氧化反应），气体分子的G需要特别处理。以CO氧化反应CO + 1/2 O₂ → CO₂为例：

• 固体催化剂表面用VASP+Phonopy计算G_solid(T)
• 气体分子G_gas(T)从NIST-JANAF热化学表查取，或用VASP计算分子振动频率后按理想气体统计力学公式算
• ΔG_rxn(T) = G(CO₂, gas) – G(CO, gas) – 1/2*G(O₂, gas) + ΔG_surface

气体分子的平动和转动熵贡献巨大（室温下CO₂的S_trans≈155 J/mol·K），不能从DFT直接计算——这是初学者最容易犯的错误。