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 生成含空位结构#

generate_structure.py 用ASE生成FCC镍，再删除 (3.52,3.52,3.52) 位置的原子。初始相邻原子在 (3.52,5.28,5.28) 附近，它在删除原子后对应的LAMMPS ID是53。ID来自程序实际映射，不是随意选择的“第53行”。

from pathlib import Path
import os,json,hashlib
import numpy as np
from ase.build import bulk
from ase.io import write
B=Path(__file__).resolve().parent;W=B/'work';assert not W.exists(),'work已存在，请重新解压';W.mkdir();atoms=bulk('Ni','fcc',a=3.52,cubic=True).repeat((3,3,3));vac=np.array([3.52,3.52,3.52]);hop=np.array([3.52,5.28,5.28]);vi=int(np.argmin(np.linalg.norm(atoms.positions-vac,axis=1)));hi=int(np.argmin(np.linalg.norm(atoms.positions-hop,axis=1)));assert np.linalg.norm(atoms.positions[vi]-vac)<1e-12 and np.linalg.norm(atoms.positions[hi]-hop)<1e-12;del atoms[vi];hid=hi+1-(hi>vi);write(W/'initial.data',atoms,format='lammps-data',atom_style='atomic',units='metal',masses=True);model=json.loads((B/'MODEL.json').read_text());model.update(atoms=len(atoms),hop_atom_ID=hid,cell_A=atoms.cell.array.tolist(),vacancy_initial_A=vac.tolist(),vacancy_final_ideal_A=hop.tolist());(W/'MODEL.json').write_text(json.dumps(model,indent=2)+'\n');assert hashlib.sha256((W/'initial.data').read_bytes()).hexdigest()==hashlib.sha256((B/'seeds/initial.data').read_bytes()).hexdigest();print('GeneratedFCCNi108minusonevacancy=107atoms;hopID',hid)

整包任务已调用这个程序，无需另外重复运行。生成 work/initial.data 和 work/MODEL.json；data完整107原子列表由程序写出。可在提交前查看固定参考结构：

head -24 seeds/initial.data
cat MODEL.json
head -8 DATA/Ni_u3.eam

势函数是 Ni_u3.eam 的Ni funcfl，保留原作者与文献头，来源为Foiles等人的EAM参数（Phys. Rev. B 33, 7983, 1986），并保留SHA-256。eam 的funcfl会读取势文件自带质量58.71；写在 pair_coeff 前面的 mass 会被它替代。若确需更改质量，应在势参数读取之后设置，并考虑动力学和优化算法的影响。

2 分别优化两个端点#

初态完整输入如下；固定晶胞，仅优化原子坐标：

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 100
thermo_style custom step atoms pe fmax fnorm press
thermo_modify norm no format float %.16g
min_style fire
timestep .001
minimize 0.0 1.0e-6 4000 20000
write_data endpoint.data nocoeff
print "ENDPOINT $(pe:%.16g) $(fmax:%.16g) $(fnorm:%.16g)" file endpoint.dat screen no

末态从初态已优化的 endpoint.data 出发，把同一个ID53移到原空位附近，再独立优化。生成的末态输入如下：

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 100
thermo_style custom step atoms pe fmax fnorm press
thermo_modify norm no format float %.16g
min_style fire
timestep .001
set atom 53 x 3.52 y 3.52 z 3.52
minimize 0.0 1.0e-6 4000 20000
write_data endpoint.data nocoeff
print "ENDPOINT $(pe:%.16g) $(fmax:%.16g) $(fnorm:%.16g)" file endpoint.dat screen no

这里 minimize 0.0 1.0e-6 4000 20000 使用力条件，端点实际总力范数约 8e-7 eV/Å，两端能量分别约 −474.5177831570 eV，差约 2e-11 eV。等价端点能量相近符合这个对称小模型；真实不等价位点的两端能量不必相同。

流程会按原子ID排序读取末态坐标，生成107行 final.coords，保留相同盒子、原子数、类型和ID。

cat work/cases/initial/endpoint.dat
cat work/cases/final/endpoint.dat
head -6 work/cases/images-9/final.coords

3 设置普通NEB与爬山图像#

9图像完整输入：

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

fix band all neb 1.0 将所有107个原子纳入路径，允许邻近原子一起调整；并非只移动ID53而固定全部邻居。neb ... 4000 4000 分别限制普通与爬山阶段，1e-3 eV/Å 是这次教学NEB的力阈值。每个图像1个MPI进程。

提交脚本会顺序完成5与9图像，不需要手动敲很长的MPI命令。若要查看启动关系，在已经分配的计算资源内，9图像原生调用是：

mpirun --bind-to none -np 9 lmp -partition 9x1 -in in.lammps

图像数、分区数和MPI资源必须匹配。换成自己的体系时，先检查初末态映射和初始插值是否穿过不合理近邻，不要直接靠增加迭代次数处理原子重叠。

4 得到路径、能垒与结构#

图中展示实际x≈3.52 Å原子层的yz投影，橙色始终是同一个跃迁原子，虚线圆为两个理想晶格位置。空位没有原子标识，而是从原来缺原子的位置移到跃迁原子留下的位置。

势能垒约1.06085 eV。两组能垒差很小，但折线路径长度仍有差异；材料定量研究还需检查更大超胞、更多图像、不同初始路径、势函数适用范围及更严格的过渡态分析。

work/analysis/ni-vacancy-neb.csv 有每个图像的相对能量、集体路径位置、跃迁原子坐标与修正力。export_path.py 导出每帧107原子的XYZ，可用OVITO等工具打开；帧序号代表路径位置，不是时间轨迹。

5 怎样迁移到自己的材料#

本包默认生成器固定为107原子Ni，并核对参考结构哈希。建立其他体系时，准备自己的完整原生data、势文件和已优化端点；同步修改元素映射、原子ID、盒子、边界和分析假设。不要只更换元素名字，也不要把本例1.06085 eV当作其他Ni尺寸或其他金属的默认数值。

下一篇：端点、力和P100单点核对给出独立能量差分、收敛曲线与真实精度数据。解析模型NEB入门可用于先理解分区和爬山流程。

官方输入说明：neb、EAM势文件。