pymatgen：读取真实NEB路径和能垒

本篇用pymatgen 2026.9.24读取本集群已完成的Al(100)表面直接跳跃路径。模型为2×2三层Al表面12原子加1个吸附Al，底部8原子固定；两个独立优化端点和1个CI中间图像，共3个图像。PBE PAW、ENCUT=450 eV、EDIFF=1E-6、2×2×1 Gamma网格，优化力0.02 eV/Å，实际VASP 6.6.1 GPU计算。完整计算步骤见CI-NEB表面迁移。

这里仅用已有输出做新的读取与曲线分析，不重新提交DFT。它是一条直接跳跃候选路径，不证明它比交换机制等所有路径的能垒都低，也不提供收敛扩散系数。

下载输入、真实参考表与操作脚本。

unzip pymatgen-NEB-wave666-*.zip
cd pymatgen-NEB
bash results.sh
bash plot.sh

reference-data/保存实际计算的数值，figures/reference.png为重绘图；读取和绘图不会启动VASP。若要运行本篇全部后处理，普通账号提交sbatch job.slurm。脚本申请cpu-8173m独占整节点1任务×56核，小例内部单线程，计费仍按整节点；不要在管理节点启动计算。

1. 自己的输出应有什么#

own-neb/
├── 00/  POSCAR  OUTCAR
├── 01/  CONTCAR OUTCAR
└── 02/  POSCAR  OUTCAR

端点读取POSCAR，中间图像读取最终CONTCAR；所有图像需要各自OUTCAR。多个中间图像时按00、01、02……排序，不能只检查一个图像的收敛或把不同路径混成同一序列。

source env.sh
python read_own_neb.py /自己的完整NEB目录

请将中文示意路径换成实际目录。脚本调用原生接口并把提取的数值与结构写到own-NEB-numeric-input.json，原始OUTCAR保持不变：

from pymatgen.analysis.transition_state import NEBAnalysis
a = NEBAnalysis.from_dir('own-neb')
print(a.r)
print(a.energies)
print(a.relative_energies)
print(a.forces)

若端点能量放在独立优化目录，使用relaxation_dirs=('start-relax','end-relax')指定对应端点。结构序号、原子数与元素顺序仍须和路径一致。

2. 附件里的真实例子怎么运行#

公开包提供从本集群原始输出提取的数值与结构，不分发全文OUTCAR。calculate.py用这些明确标注的真实数值重建原生NEBAnalysis对象：

sbatch job.slurm
head results/p100-curve.csv

核心操作为：

from monty.serialization import loadfn
from pymatgen.analysis.transition_state import NEBAnalysis
d = loadfn('reference-data/p100-NEB-numeric-input.json')
a = NEBAnalysis(d['r_A'], d['energies_eV'],
                d['tangent_forces_eV_A'], d['structures'])
print('discrete barrier / eV:', max(a.relative_energies))

这一步没有伪造VASP输出；输入是原始OUTCAR中实际的energy(sigma->0)、切向力和最终结构，并与原生from_dir读取逐项对应。

3. Eσ→0、路径长度与能垒#

P100与V100 16G原始路径独立运行，本次读取结果为：

这里比较的是同一个小模型的实际路径数值，不是两种GPU的性能结论。原生接口读取energy(sigma->0)，并减去起点能量，不能把混用TOTEN和Eσ→0造成的差当成能垒变化。

当前实现对相邻图像计算各原子最小周期距离，再取√Σᵢ|Δrᵢ|²，沿图像累加。它没有再除以原子数，所以不要把这条配置空间坐标叫成“单个迁移原子的位移”或每原子RMS距离。

4. 接下来如何理解曲线#

插值只连接现有图像，没有新增电子计算点。NEB样条的实际约定继续说明默认端点斜率、内部切向力和极值采样；这些差别在图像很少时尤其重要。