VASP计算自由能：从振动频率到声子谱的完整方法

在第一性原理计算中，自由能是连接电子结构层次与热力学宏观性质的核心桥梁。VASP本身不直接输出”自由能”这一个量，而是通过振动频率分析或声子谱计算，将有限温度下的振动熵贡献纳入能量修正，最终得到亥姆霍兹自由能F(T)。

自由能的物理内涵与计算路径

晶体的亥姆霍兹自由能可以分解为三部分：

F(T) = E_DFT + F_vib(T) + F_el(T)

其中E_DFT是0K下的静态DFT总能量，F_vib是晶格振动（声子）贡献的振动自由能，F_el是电子热激发贡献（金属体系不可忽略）。对于绝缘体和半导体，F_el通常可以忽略，重点在于F_vib的计算。

VASP计算F_vib有两条路径，选择哪条取决于体系大小和精度需求：

路径一：有限位移法（ Phonopy + VASP）

这是目前最主流的方案。Phonopy作为外部工具驱动VASP完成超胞力计算，流程清晰：

1. 用VASP对原始cell做高精度结构优化（EDIFF=1E-8, EDIFFG=-1E-4）
2. 用Phonopy生成位移超胞（通常位移量±0.01 Å）
3. 对每个位移超胞做VASP静态计算（EDIFF=1E-8, PREC=Accurate）
4. 用Phonopy提取力常数，计算声子色散和态密度
5. 在准谐近似（QHA）下，对不同体积重复上述步骤，得到F_vib(T,V)

关键参数设置上，超胞大小直接影响声子谱的精度。对于简单金属（如Al、Cu），2×2×2超胞通常足够；对于低对称性晶体或含轻元素体系（如氢化物），建议3×3×3超胞并加密k网格。

路径二：密度泛函微扰理论（DFPT）

VASP内置了DFPT功能（LEPSILON=.TRUE.或IBRION=6/7/8），可以直接在原始cell中计算力常数矩阵，无需超胞。优势在于计算量小，尤其适合大体系。但DFPT对k网格密度更敏感，且部分泛函（如HSE06）不支持线性响应。

INCAR参数配置要点

以Phonopy+VASP的有限位移法为例，位移超胞的INCAR核心参数：

PREC = Accurate
EDIFF = 1E-8
ENCUT = 520
IBRION = -1
NFREE = 0
NSW = 0
ISMEAR = 0
SIGMA = 0.01
ADDGRID = .TRUE.
LREAL = .FALSE.

这里IBRION=-1且NSW=0表示不做结构优化，只计算单点力。LREAL=.FALSE.强制在倒空间投影，保证力计算的精度。EDIFF=1E-8比常规计算严格两个量级，因为声子频率对力常数非常敏感，能量收敛不够会直接导致虚频。

如果用DFPT路径，INCAR需要改为：

IBRION = 8
LEPSILON = .TRUE.
PREC = Accurate
EDIFF = 1E-8
NSW = 1

IBRION=8对应有限电场计算，会自动计算Born有效电荷和介电常数，适合后续做热力学修正。

从声子谱到自由能的转换

拿到声子频率后，振动自由能在简谐近似下的表达式为：

F_vib(T) = k_B * T * Σ_i ln[2 * sinh(ℏω_i / 2k_B * T)]

Phonopy可以直接输出这个量的温度依赖曲线。需要特别注意的几个实战问题：

虚频处理：如果声子谱出现虚频（Γ点附近的ω²<0），说明结构没有在真正的能量极小点上。常见原因是结构优化不够充分，或者ENCUT偏低导致Pulay应力。解决办法是重新优化，增大ENCUT到至少1.3倍ENMAX，并用ADDGRID=.TRUE.。

QHA适用范围：准谐近似假设振动频率只是体积的函数，温度效应通过体积变化间接体现。在高温下（德拜温度的2倍以上），非谐效应显著增强，QHA会系统性偏高自由能。此时需要考虑声子-声子相互作用的三阶力常数，用Phonopy-3ph或ALAMODE做非谐修正。

实际案例：MgO体系的热膨胀

以MgO为例说明完整流程。MgO是离子晶体，常温德拜温度约743K，QHA在1000K以下适用良好：

• 结构优化后用2×2×2超胞（64原子）做有限位移计算
• Phonopy输出声子色散，Γ点光学模频率约383 cm⁻¹，与实验值382 cm⁻¹吻合
• 在5个体积点（±3%应变）重复计算，得到F(T,V)曲面
• 对F(T,V)做Birch-Murnaghan状态方程拟合，提取热膨胀系数
• 300K下得到α=3.1×10⁻⁵ K⁻¹，实验值3.1×10⁻⁵ K⁻¹，误差<1%

这个案例说明，只要参数设置得当，VASP+Phonopy计算自由能的精度可以完全达到实验测量水平。

计算成本与效率权衡

对于不同体系规模的计算量估算：32原子的超胞，单个位移计算（VASP静态）在16核节点上约2-4小时；2×2×2超胞通常需要12-24个位移构型，总计算时间约24-96核时。如果体系更大（如含64+原子的超胞），建议用DFPT路径减少计算量。

自由能计算的核心在于力常数矩阵的精度，而不是泛函的选择。PBE+PAW在大多数体系已经足够，不需要换用更贵的HSE06。把算力集中在ENCUT、EDIFF和k网格上，比换泛函的性价比高得多。