VASP 状态方程与高压相比较

本篇解决两个实际问题：“材料压缩起来有多难？”以及“给定压力下，两种结构哪个更有利？”前者从体积—能量曲线拟合，后者比较同一压力下的焓。

本组教程：1. 结构优化 · 2. 能量差 · 3. 相稳定性 · 4. 状态方程与压力 · 5. 表面与界面 · 6. 结构搜索。

1. 状态方程怎样算#

以8原子金刚石Si晶胞为例：先准备接近平衡的结构，再取参考体积的0.88、0.94、1.00、1.06、1.12倍。体积乘f时，三条晶胞矢量应乘 f的立方根，不是每条都乘f。

1. 复制出五个目录，各自按体积比例缩放晶胞和原子。
2. 保持每个目录的目标体积不变。一般材料可以固定体积优化内部位置及形状；本例理想金刚石Si的内部位置由对称性确定，采用固定形状静态计算。
3. 所有目录统一PBE、ENCUT=550 eV、EDIFF=1E-6和Gamma中心7×7×7网格，设置 NSW=0、IBRION=-1。
4. 提取每个目录的最终体积V和同一种能量E，统一到每原子。
5. 画出E(V)，确认最低点在取样区间内，再拟合Birch–Murnaghan状态方程。

下载本组完整算例与数据（ZIP） · 下载实算结果表（CSV）。

对应目录为 eos-diamond-88 至 eos-diamond-112。准备好POTCAR后，在每个目录中执行 sbatch run.sh。五个任务彼此独立。

2. 本次Si的体积点与拟合#

拟合得到：

五点拟合适合入门；B₀′通常比V₀更敏感。如果需要可信的压力导数或材料间的小差别，应增加体积点并检查拟合区间、k点和ENCUT。不要仅凭拟合误差很小就认为全部物理误差都小。

左图同时给出β-Sn型Si的固定形状体积扫描。它是指定结构分支的曲线，不能仅凭两个零压最低能量判断高压稳定性；高压比较还要考虑PV及晶胞随压力的调整。

3. 给定压力下，应该比较焓#

零温静态近似下：

$$H(P)=E+PV.$$

E用eV、V用ų、P用GPa时，PV的eV值为 P×V/160.21766208。比较不同原子数的晶胞时，将E、V和H全部除以原子数。

本次分别在0、10、20 GPa优化金刚石型和β-Sn型Si。对应优化目录为 pressure-diamond-0/10/20 和 pressure-beta-0/10/20。仍使用EDIFF=1E-6、EDIFFG=-0.02、ENCUT=550 eV、7×7×7网格。

以10 GPa为例：

ISIF = 3
IBRION = 2
PSTRESS = 100

PSTRESS的单位是kbar，10 GPa应填100，而不是10。 优化后用CONTCAR建立相应的 -static 目录，固定几何、去掉PSTRESS，重新计算不含外加PV项的电子能，再手动加PV。这样可以避免把压力优化输出里已包含的PV项再加一次。

操作顺序是：压力优化 → 检查结构与残余应力 → 复制最终结构 → 零PSTRESS静态计算 → 统一每原子焓 → 同压比较。

4. 两种Si结构的实算比较#

ΔH定义为H(β-Sn型)−H(金刚石型)：正数表示本次比较中金刚石型更低，负数表示β-Sn型更低。本次0 GPa时金刚石型更低，10与20 GPa时β-Sn型更低，因此这两条结构分支的焓交叉位于0–10 GPa之间。三个点只能限定区间；要精确定位，需要在区间内增加压力点。

这不是完整高压相图。 本例只比较两种起始结构，未包含所有高压相，也没有声子稳定性、零点能和温度贡献。结构优化过程中还应确认最终结构是否保留预期相；如果它变成另一种结构，应按最终结构重新命名和分析，不能只凭目录名判断物相。

5. 常见操作错误#

官方说明：体积优化与EOS、PSTRESS。