VASP 准谐自由能与热膨胀入门
在五个体积分别算硅的声子，理解自由能、平衡体积和热膨胀的关系。
阅读本篇目录
单个体积的声子可以给出定容热力学；想知道升温后晶胞如何变化，需要比较多个体积。准谐近似把温度效应通过声子自由能带入体积选择。
1. 整个工作流#
- 取优化体积V₀，构造0.96、0.98、1.00、1.02、1.04倍的体积。各晶格矢量乘以体积比例的立方根，不能直接乘体积比例。
- 在每个固定体积下优化必要的内部原子。本例高对称硅的内部位置由对称性固定。
- 每个体积做静态能量，另做一套有限位移声子计算。
- Phonopy给出各体积的振动自由能Fvib(V,T)。统一单位后加上静态E(V)。
- 固定一个温度，拟合E+Fvib随体积的曲线，寻找极小值；换下一个温度重复。
- 得到V(T)，进而计算体膨胀系数β=(1/V)dV/dT。立方材料的线膨胀系数约为β/3。
本例静态原胞采用9×9×9 k网格，声子采用2×2×2超胞、4×4×4 k网格和16³ q网格。五组数据必须对应同一种泛函、赝势与数值设置。
2. 单位与拟合#
本例Fvib以每摩尔2Si原胞的kJ/mol输出，除以96.485332转换为每原胞eV；静态能量也必须是同一个2原子原胞。不能把16原子超胞的总能量直接加到2原子的自由能上。
这里在窄体积区间用五点二次拟合展示流程。正式计算可用Phonopy的QHA工具拟合状态方程，并扩大体积采样检查拟合是否稳定。有外压时加入pV项，注意GPa·Å³到eV的转换。
3. 实际结果#
| T / K | 平衡体积 / ų（2 Si） | 曲率得到的B / GPa |
|---|---|---|
| 0 | 41.12634 | 89.926 |
| 100 | 41.11640 | 89.528 |
| 300 | 41.12380 | 87.766 |
| 500 | 41.17962 | 86.154 |
| 800 | 41.29121 | 83.892 |
所有列出的自由能极小值都在采样范围内。0 K结果含零点振动，所以不要求与纯静态E(V)的最小值一致。本例低温有轻微体积回缩，随后随温度升高膨胀；不能仅凭这组小超胞数据宣称高精度复现实验热膨胀。
4. 哪些情况要换方法#
如果某体积有明显不稳定声子，不能把其热力学直接塞入普通QHA。结构相变、软模、强非谐、高温扩散和熔化也可能超出准谐近似。Cp与Cv不同；需要Cp时用完整QHA的温度导数与热膨胀关系，不把本页上一教程的Cv改个名称就交付。
方法参考：Phonopy准谐近似。
计算环境与资料下载#
本页数据来自鸢算科技集群实际计算。VASP 6.6.1 GPU版采用 NVHPC 25.3、CUDA 12.8、HPC-X 2.22.1 与 MKL 2023.2；PBE、PAW PBE 版本64、ENCUT=550 eV、EDIFF=1E-6。需要优化的自由原子采用0.02 eV/Å；静态位移点、外加应变点和动力学轨迹不按结构优化阈值判断。
下载输入、分析脚本和实算数据（ZIP） · 下载结果表（CSV）
解压后先读 README，选择本文对应目录，按POSCAR元素顺序准备本人有权使用的POTCAR，再运行 sbatch run.sh。POTCAR、WAVECAR及电荷密度大文件不随资料包分发。各目录的原始输入和实际结果表用于复现；更换体系后应重新提取自己的结果。