VASP 刚性解理：从结构到能量曲线
用8原子Si的七个实际DFT单点，演示开裂间距扫描、Slurm数组、能量提取和面积换算。
阅读本篇目录
下载七点结构、输入、实际能量和作图脚本，保存为 example.zip：
python3 -m zipfile -e example.zip .
cd vasp-rigid-cleavage-examples
bash results.sh
bash plot.sh
本例使用VASP 6.6.1，NVHPC 25.3 / CUDA 12.8及配套MPI构建，在P100实际计算七个间距点，并用V100 16G计算同一零间距参照。每个任务一张卡、一个MPI进程，Slurm数组负责调度。参考数据在 reference-data，重算目录为 work。
1. 这个例子计算什么#
取5.43 Å的金刚石Si常规晶胞，共8个原子。保持原子之间的笛卡尔位置关系不变，沿z增加晶胞长度，在周期像之间打开(001)方向的间隙。七个间距为0、0.5、1、2、4、8、12 Å。
每次还沿z整体平移间距的一半，让薄层位于晶胞中部；这次整体平移不改变内部几何。没有对表面原子做弛豫，也没有让晶胞重新优化。
这是学习刚性解理路径的极小模型。它不包括裂尖、裂纹扩展、塑性、表面重构或温度采样，不能输出断裂韧性KIC。薄层只有四个原子层，k网格也只用于本次流程教学。
2. 查看输入结构与计算参数#
零间距POSCAR为：
Own Si8 rigid cleavage gap 0 A
1.0
5.43 0 0
0 5.43 0
0 0 5.43
Si
8
Cartesian
0.000000000000 0.000000000000 0.000000000000
0.000000000000 2.715000000000 2.715000000000
2.715000000000 0.000000000000 2.715000000000
2.715000000000 2.715000000000 0.000000000000
1.357500000000 1.357500000000 1.357500000000
1.357500000000 4.072500000000 4.072500000000
4.072500000000 1.357500000000 4.072500000000
4.072500000000 4.072500000000 1.357500000000
开裂2 Å时，c改为7.43 Å，所有原子z坐标同时增加1 Å，其余坐标不动。可以直接比较输入：
cat inputs/gap-00/POSCAR
cat inputs/gap-03/POSCAR
diff -u inputs/gap-00/POSCAR inputs/gap-03/POSCAR
本例使用Cartesian坐标，避免增加c时把内部层间距一起拉伸。若只改Direct坐标结构的c而保留相同分数坐标，得到的是整体拉伸，物理路径已经改变。
实际INCAR：
SYSTEM = gap-00
ENCUT = 450
PREC = Accurate
EDIFF = 1E-6
IBRION = -1
NSW = 0
ISIF = 2
ISMEAR = 0
SIGMA = 0.05
ISPIN = 1
ISYM = 0
ALGO = Normal
NELM = 180
LREAL = .FALSE.
LASPH = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
EDIFF=1E-6控制电子收敛。NSW=0、IBRION=-1做固定结构单点，ISIF=2输出应力；这里没有离子优化步骤，因此没有使用0.02 eV/Å作为终止条件。若希望得到约束下弛豫后的解理曲线，应另建任务，明确哪些原子能动，再用 EDIFFG=-0.02，不能直接给本结果改名。
实际KPOINTS：
Own fixed in-plane mesh;teaching only
0
Gamma
5 5 1
0 0 0
3. 准备POTCAR并提交七个单点#
在解压目录根部放好自己有权使用、与输入匹配的PBE Si POTCAR。本包不包含POTCAR和波函数。然后执行：
python3 prepare.py
sbatch job.slurm
squeue -u "$USER"
prepare.py把七套输入复制到 work/gap-00至 work/gap-06。为避免覆盖结果，已有 work 时脚本会停止。
提交脚本的关键资源设置如下：
#SBATCH --partition=gpu-8p100-hs
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --gres=gpu:p100-16g:1
#SBATCH --array=0-6%3
每个数组任务运行一个间距点，最多三个同时运行，实际开始时间由队列资源决定。每个任务内部调用 mpirun -np 1 vasp_std，不是把多个数组任务挤到同一张卡。
4. 确认完成，再提取数据#
等七个任务都结束后执行：
python3 collect.py
cat results/Si-rigid-cleavage.csv
bash plot.sh --data results/Si-rigid-cleavage.csv
提取脚本要求每个OUTCAR都包含正常结束和达到EDIFF的记录。这里采用最终的 energy(sigma->0)，全部点采用同一能量定义，不能把某些点的TOTEN与其他点的E₀混用。
| 开裂间距/Å | E₀/eV（8 Si） | 相对能量/eV | 刚性解理功/J·m⁻² |
|---|---|---|---|
| 0 | -43.40461216 | 0.00000000 | 0.000000 |
| 0.5 | -42.29410634 | 1.11050582 | 0.603436 |
| 1 | -40.03299782 | 3.37161434 | 1.832098 |
| 2 | -36.99289872 | 6.41171344 | 3.484054 |
| 4 | -36.07578313 | 7.32882903 | 3.982404 |
| 8 | -36.13667354 | 7.26793862 | 3.949317 |
| 12 | -36.13709293 | 7.26751923 | 3.949089 |
本次八个任务合计约120 GPU秒，单个实际执行约10–24秒；这是本例规模和当时环境下的耗时，不是其他体系或排队时间的承诺。P100与V10016零间距参照在OUTCAR打印的能量精度内一致。
5. 面积、两张表面和单位#
本例截面积 A=|a×b|=29.4849 Ų。刚性解理功为 W(δ)=[E(δ)−E(0)]/A，换算使用 1 eV/Ų=16.02176634 J/m²。
12 Å点得到 3.949089 J/m²。这是本次未弛豫路径的能量增量除以面积。裂开会形成两张表面；若要在合适条件下讨论单面表面能，还要除以2，并处理表面弛豫、厚度、k点和参考态。不能把本例的W直接标成单面表面能γ。
4到8 Å的曲线有小幅回落，8到12 Å变化较小。这里保留真实点，没有强行拟合成单调曲线。小晶胞、固定原子、有限真空与k采样都会影响数值；要研究材料定量结果，应增加厚度与数值收敛比较。
6. 与拉伸和滑移有什么区别#
整体拉伸会改变层内、层间间距；刚性解理只打开选定界面；滑移主要改变界面两侧的平行位移。三种路径的受限自由度不同，不能只换图标题。继续看二维应变加载、刚性滑移与力学和有限温度工作流。晶胞与原子自由度见 VASP ISIF官方说明。