LAMMPS空位NEB结果核对

本篇把NEB流程放到一个实际金属模型中：3×3×3的FCC镍共有108个晶格位置，删除一个原子形成空位，余下107个原子。相邻镍原子跳入空位，空位便移到相邻位置。我们实际完成了两个端点优化、5/9图像NEB、能量与力检查，以及一个P100单点对照。

这里得到的是特定EAM势、固定10.56 Å立方盒中的教学结果，不是经过尺寸、势函数和温度收敛的体相实验预测。它适合作为“自己的结构该怎样接入NEB”的操作参考。

下载结构生成器、完整输入、势文件、实际数据与图片，在自己的工作目录执行：

mkdir ni-vacancy
cd ni-vacancy
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-ni-vacancy-wave79-f965ac3b49ce.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

作业脚本自动生成结构并顺序执行全部步骤，申请9个CPU核、一个P100。NEB用CPU/MPI计算；P100用于最后同一结构的 eam/gpu 单点核对，本例没有测试GPU版NEB加速比。LAMMPS模块为 lammps/30Sep2026-p100-gpu：30Sep2026、GCC13.3.0、OpenMPI5.0.10、CUDA12.8、mixed GPU精度；结构与绘图使用 ase/3.29.0。

先读作业输出，再提取结果：

module purge
module load ase/3.29.0
python analyze.py work
python export_path.py work
python plot_results.py work
head -10 work/analysis/ni-vacancy-neb.csv

已有实算位于 reference/，可将后处理命令中的 work 换成 reference 直接查看。重新计算换一个新目录；整包提交会自行创建work，不要先手动生成后又在同一目录重复提交。

1 核对空位与端点#

输入有107个原子、108个理想晶格位置。分析程序用周期性最小镜像距离找出缺原子的位置：初态在 (3.52,3.52,3.52)，末态在 (3.52,5.28,5.28)；其余107个晶格位置都有附近原子。整个路径保留原子ID和类型，也检查最近粒子距离，避免画出一条看似平滑却映射错误的路径。

完整9图像NEB输入如下，便于对照自己的流程；端点优化的两份输入见Ni空位完整流程教程：

variable image uloop 9
units metal
atom_style atomic
atom_modify map array
boundary p p p
read_data initial.data
mass 1 58.6934
pair_style eam
pair_coeff * * ../DATA/Ni_u3.eam
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
thermo_modify norm no
fix band all neb 1.0
min_style quickmin
timestep .002
neb 0.0 1.0e-3 4000 4000 20 final final.coords verbosity default
write_dump all custom neb-${image}.dump id type x y z fx fy fz modify sort id format float %.16g
unfix band
run 0
print "ENERGY ${image} $(pe:%.16g)" file energy-${image}.dat screen no
write_dump all custom physical-${image}.dump id type x y z fx fy fz modify sort id format float %.16g

2 看NEB修正力，不把普通力混在一起#

每个分区的普通阶段与爬山阶段都达到力阈值。neb-编号.dump 保存NEB修正力，随后移除 fix band、执行 run 0，physical-编号.dump 保存同一位置的普通EAM物理力。

中间路径图像沿路径有势能斜坡，普通力可能不为零；比较输入阈值时使用NEB修正后的力。端点应先优化，爬山图像则需要靠近候选过渡态。

grep 'Stopping criterion' work/cases/images-9/log.lammps.4
head -12 work/cases/images-9/neb-5.dump
head -12 work/cases/images-9/physical-5.dump

右图把单个跃迁原子的进度与全体系路径坐标对照，邻居也在调整，所以两个坐标定义不必相同。

3 独立位移能量差分核对力#

对9图像中间结构，程序用原生 read_dump 重放坐标，在三个笛卡尔方向分别把跃迁原子位移 ±0.0001 Å，独立计算六个能量，再用中心差分得到力。普通EAM力与能量差分最大差约 3.34e-9 eV/Å。

查看实际数据并直接提取：

cat work/cases/control/center.dat
cat work/cases/control/energy-fd.dat
module load ase/3.29.0
python - <<'PY'
from pathlib import Path
import numpy as np
p=Path('work/cases/control')
center=np.loadtxt(p/'center.dat',usecols=(1,2,3,4))
rows=[s.split() for s in (p/'energy-fd.dat').read_text().splitlines()]
for i,axis in enumerate(['x','y','z']):
    pair=[r for r in rows if r[1]==axis]
    force=-(float(pair[1][3])-float(pair[0][3]))/0.0002
    print(axis,'energy-FD force=',force,'native=',center[i+1])
PY

程序还沿理想跃迁方向 (0,1,1)/√2 做了两个位移，得到负曲率约 −2.466 eV/Ų。这只检查一个方向，不是完整Hessian，也没有据此证明全部振动模式中只有一个虚频。需要严格过渡态确认时，还应补相应振动或Hessian分析。

4 同一结构的P100单点精度#

使用同一中间结构，分别调用CPU eam 与P100混合精度 eam/gpu，不再移动原子。实际P100输入如下：

package gpu 1 neigh no
newton off
units metal
atom_style atomic
atom_modify map array
boundary p p p
read_data initial.data
read_dump saddle.dump 110 x y z box no replace yes
mass 1 58.6934
pair_style eam/gpu
pair_coeff * * ../DATA/Ni_u3.eam
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
thermo_modify norm no format float %.16g
run 0
print "GPU $(pe:%.16g)" file gpu-energy.dat screen no
write_dump all custom gpu.dump id type x y z fx fy fz modify sort id format float %.16g

这些差异属于这个势、结构和混合精度构建的实际观察；不要用一个结构的精度对照推断所有体系，更不能把它当成GPU版NEB的性能测试。本例NEB一直使用CPU/MPI路径。

5 保存哪些结果#

保留初末态data、final.coords、全部图像与分区日志、势文件版本及哈希、输入脚本和分析CSV。XYZ仅方便看结构，不能代替原生数据、周期映像或重启文件。图件来自实际输出，可点击原图放大。

python export_path.py work
head -8 work/analysis/images-9-Ni-path.xyz
python check_results.py work

返回完整Ni空位计算流程。收敛与力的定义见官方fix neb，EAM格式见官方EAM说明。