VASP非谐振动势能曲线
用硅的冻结光学位移，比较二次、四次与六次项，并用原子力检查拟合。
阅读本篇目录
谐近似把振动势能写成二次函数。位移增大后，高阶项可能影响振动频率和热力学。本例选择硅两原子原胞中的一个光学位移坐标，直接计算势能曲线，说明如何从VASP数据识别非谐性。
1. 先定义位移，不做结构优化#
使用既有优化硅原胞，固定晶胞。两个Si分别沿笛卡尔x移动+u、−u，其他坐标固定。u取−0.16至+0.16 Å，间隔0.04 Å，共9组。
这些是指定坐标上的单点采样，若重新优化原子，就会丢失要研究的势能曲线。PBE PAW64、ENCUT=550 eV、EDIFF=1E-6、9³Gamma网格、ISMEAR=0、SIGMA=0.05 eV、ISYM=0。使用vasp/6.6.1-p100，NVHPC 25.3/CUDA 12.8构建。
附件calculations目录包含所有输入。按POSCAR准备合法Si POTCAR，分别sbatch run.sh，得到能量和力。
2. 提取能量与对应广义力#
完成后在附件根目录执行：
module load ase/3.29.0
python extract_anharmonic_mode.py
能量标签采用与原子力一致的VASP自由能F，并减去u=0参照。因为两个原子反向移动，广义力为Fᵤ=F₁ₓ−F₂ₓ=−dE/du，不能只取一个原子的力。
| u / Å | 相对能量 / meV原胞⁻¹ | 广义力 / eV Å⁻¹ |
|---|---|---|
| -0.160 | 618.0132 | +7.23561 |
| -0.120 | 357.7725 | +5.72412 |
| -0.080 | 163.5967 | +3.98284 |
| -0.040 | 41.6158 | +2.07426 |
| +0.000 | 0.0000 | +0.00000 |
| +0.040 | 41.6158 | -2.07426 |
| +0.080 | 163.5967 | -3.98284 |
| +0.120 | 357.7725 | -5.72413 |
| +0.160 | 618.0131 | -7.23561 |
3. 用偶次多项式拟合#
先比较+u与−u。该方向按对称性应有相同势能；本例最大奇对称能量部分为0.000005 meV。对正负两侧能量取平均后拟合：
E(u)−E(0)=a₂u²+a₄u⁴+a₆u⁶。
| 系数 | 本例结果 |
|---|---|
| a₂ / eV Å⁻² | 26.2460128 |
| a₄ / eV Å⁻⁴ | -116.4832676 |
| a₆ / eV Å⁻⁶ | 1338.2543505 |
| 能量拟合RMSE / meV原胞⁻¹ | 0.045214 |
| 未参与拟合的广义力RMSE / eV Å⁻¹ | 0.047818 |
这里系数定义不含1/2!、1/4!、1/6!。若需要对应的势能导数，应乘相应阶乘。只用能量拟合后，再用实际力检验导数，能发现只把几个能量点拟合得很好、力却不匹配的情况。
4. 最后得到什么#
得到这一个坐标的非谐势能曲线、拟合系数、能量残差和力检验。a₂对应的单坐标谐振动频率约15.1128 THz；本例两个原子反向移动，有效质量是2mSi。
这不是完整四阶力常数张量，也没有计算四声子散射率或四声子热导。完整晶体非谐输运需要多个原子、方向、波矢之间的耦合与更充分采样。这里的多项式只在已计算位移范围内使用，不向任意大位移外推。进一步可用它演示单坐标非谐自由能。
数值方法参考：SciPy对称三对角矩阵求解。