VASP单坐标非谐自由能示例
从实际势能曲线建立一维量子模型，比较谐近似与非谐振动自由能。
阅读本篇目录
有了势能曲线，还可以计算指定模型的振动能级和自由能。本例接着硅冻结位移的9组VASP结果，建立一个两原子周期原胞的单坐标量子模型，帮助理解非谐自由能的计算步骤。
1. 明确这是哪个模型的自由能#
只保留两个Si沿x反向移动的坐标u，其余自由度固定。势能采用前篇的偶次多项式，有效质量m_eff=2mSi=56.170 u。
这是一维约化模型，不是整块硅的振动自由能。真实晶体需要考虑全部声子波矢、分支及其耦合；不能把这一原胞坐标的自由能直接乘以晶胞数，宣称得到体材料的非谐自由能。
2. 先完成VASP，再解一维方程#
按前篇完成9个单点计算，在附件根目录运行extract_anharmonic_mode.py。脚本用有限差分求解：
[−ℏ²/(2m_eff) d²/du² + V(u)] ψₙ = Eₙ ψₙ。
计算范围u∈[−0.16,+0.16] Å，两端波函数取零，保留最低48个能级。分别使用801和1201个网格点作数值对照。这个边界条件属于约化模型的数值设置，因此需要检查波函数在边界附近是否足够小。
3. 从能级计算自由能#
T>0时F=−kBT ln Σₙ exp(−Eₙ/kBT)；T=0取最低能级，包含这个坐标的零点能。谐参照使用相同a₂与有效质量：F_harm=ℏω/2+kBT ln[1−exp(−ℏω/kBT)]。
以下能量都相对于静态u=0势能，不包含其他自由度：
| T / K | 非谐模型F / meV | 谐模型F / meV | 差值 / meV | RMS位移u / Å |
|---|---|---|---|---|
| 0 | 31.1300 | 31.2508 | -0.1208 | 0.02449 |
| 300 | 28.6648 | 28.8374 | -0.1726 | 0.02686 |
| 600 | 12.5089 | 12.9166 | -0.4077 | 0.03360 |
801与1201点给出的最大自由能差为0.000497 meV；1201点计算中|u|>0.14 Å区域的最大总概率为2.442e-05。这些只检查离散网格和所取边界附近的波函数分布，不替代VASP采样密度、拟合阶数与其他振动自由度的检验。
4. 怎样扩展到实际材料问题#
本例输出能级、温度—自由能表、RMS位移和对应图片。它展示非谐势能怎样改变指定坐标的量子统计，而非给出硅的完整非谐热容、相稳定性或相变温度。
若要比较两个相的稳定性，需要每个相一致的静态能量、全部相关自由度、体积与压力项，并选择合适的非谐自由能方法。不能从这条单坐标曲线推算熔点，也不能把有限盒子的高激发态解释成真实晶体的全部声子态。
数值方法参考：SciPy对称三对角矩阵求解。