VASP计算焓值：从总能量到生成焓的完整推导

焓（Enthalpy, H）在第一性原理计算中并不是VASP直接输出的物理量，而是通过DFT总能量加上压力-体积功得到的导出量。理解这一点，是正确使用VASP做热力学分析的前提。

焓的DFT定义与计算原理

在0K、零压条件下，晶体的焓近似等于DFT总能量：

H(0K, 0Pa) ≈ E_DFT

有限温度下，焓的定义为H = U + PV，其中U是内能，P是压力，V是体积。在VASP的框架中，U由静态DFT能量E_DFT加上振动内能U_vib(T)构成，P来自应力张量的迹。因此：

H(T) = E_DFT + U_vib(T) + P*V

对于常压下的固体（P≈0），PV项通常可以忽略（典型金属的PV约0.001-0.01 eV/atom），焓的计算简化为H(T) ≈ E_DFT + U_vib(T)。振动内能U_vib通过声子谱计算获得，与自由能共享同一套力常数数据。

VASP输出中与焓相关的量

VASP的OUTCAR文件中有几个关键输出需要正确理解：

sigma→0后的能量（E0）：OUTCAR中”Free energy of the ion-electron system (eV)”这一行给出的是E_DFT – T*S_el，其中S_el是电子熵。对于绝缘体，S_el≈0，E0就是纯DFT总能量。对于金属，ISMEAR>0时电子熵不可忽略，应该取sigma→0外推后的E0值。

应力张量（stress tensor）：OUTCAR中”in kB”标记的3×3矩阵就是应力张量。对角元的平均值乘以-1就是压力P，用于计算PV项。在结构优化中（ISIF=3），VASP会自动调整cell使应力张量趋近于目标压力（默认0 kB）。

总能量（TOTEN）：OSZICAR中每步的E0值。最终收敛后的E0就是H(0K, 0Pa)的近似值。

生成焓的计算方法

实际研究中最常用的是生成焓（Formation Enthalpy），即化合物从单质形成时的焓变。以二元化合物A_xB_y为例：

ΔH_f = H(A_xB_y) – xH(A) – yH(B)

其中H(A_xB_y)是化合物的DFT总能量，H(A)和H(B)是相应单质在最稳定结构下的DFT总能量。所有能量都是在各自平衡体积下计算的，且必须使用相同的ENCUT、泛函和精度参数。

实操中的坑：

1. 单质参考态的选择：必须用最稳定的同素异形体。比如氧的参考态是O₂分子（不是O原子），需要在一个足够大的盒子（至少10Å×10Å×10Å）中做分子计算，并加偶极修正（IDIPOL=3）。铁在室温下是BCC（α-Fe），不是FCC。
2. 原子数归一：比较不同化学计量比的化合物时，需要将能量归一化到每原子的焓值。比如Al₂O₃的生成焓要除以5（2个Al+3个O），才能和MgO（除以2）比较。
3. 磁矩处理：含过渡金属的体系，必须正确设置初始磁矩（MAGMOM），否则能量可能偏高0.1-0.5 eV/atom，直接导致生成焓符号翻转。

INCAR参数配置

计算焓值不需要特殊INCAR设置，核心在于高精度结构优化：

PREC = Accurate
EDIFF = 1E-6
ENCUT = 520
ISIF = 3
IBRION = 2
NSW = 60
ISMEAR = 0
SIGMA = 0.05
ADDGRID = .TRUE.
LREAL = .FALSE.

ISIF=3同时优化原子位置和cell形状/体积。对于金属体系，ISMEAR=1（Methfessel-Paxton）比ISMEAR=0更稳定，但最终静态计算时换回ISMEAR=0做精确能量。两步法——先ISIF=3优化，再复制CONTCAR为POSCAR做ISIF=2静态计算——可以消除残余应力对能量的影响。

反应焓的实战计算

以LiCoO₂正极材料的脱锂反应为例：

LiCoO₂ → CoO₂ + Li

反应焓ΔH_rxn = H(CoO₂) + H(Li) – H(LiCoO₂)

• LiCoO₂（R-3m空间群）：优化后E_DFT = -39.85 eV/f.u.
• CoO₂（同一结构去掉Li）：E_DFT = -32.71 eV/f.u.
• Li金属（BCC）：E_DFT = -1.90 eV/atom
• ΔH_rxn = -32.71 + (-1.90) – (-39.85) = 5.24 eV

这个正值表示脱锂反应吸热，与实验测量的4.9-5.3 eV范围一致。差异来源主要是DFT对Co³⁺/Co⁴⁺局域态的描述偏差，如果加Hubbard U（U_eff=3.3 eV for Co d轨道），精度可以进一步提升到4.95 eV。

常见错误与诊断

能量未收敛导致的焓值漂移：EDIFF=1E-4在结构优化中够用，但做能量比较时必须收紧到1E-6。两个结构的能量差可能只有0.01-0.05 eV，如果EDIFF太松，能量噪声会淹没真实信号。

k网格不一致：不同结构的k网格密度必须一致（至少同一密度量级），否则截断误差不一致。一个简单规则：k1×a1 ≈ k2×a2 ≈ k3×a3 ≈ 40-50（中等精度）或60-80（高精度）。

未考虑零点能：严格的热力学分析应该加零点振动能（ZPE）修正。ZPE = 1/2 * Σ ℏω_i，对于含氢体系可以达到0.1-0.3 eV，不能忽略。