LAMMPS NEB能垒计算实操

本例先用一个能直接写出能量和力的弯曲双阱模型，演示 LAMMPS 的普通 NEB 与爬山图像流程。移动粒子从 (-1,0,0) 到 (1,0,0)，真实计算后路径向上弯曲，穿过 (0,0.5,0) 的鞍点。另一个粒子固定在 (0,3,0)，只作几何参照，不与移动粒子相互作用。

这是约化单位下的外场教学模型，不是某种材料的迁移能垒。它的优势是每个输入参数、路径、能量和力都能独立核对，再按相同流程换成自己的势函数和结构。

下载5/9图像完整输入、原始日志与分析脚本，在个人目录执行：

mkdir neb-demo
cd neb-demo
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-neb-wave78-6256297f06a4.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

脚本申请9个CPU核及一个P100资源，先跑5个、再跑9个分区，每个图像用1个MPI进程。本例原生NEB走CPU/MPI路径，没有GPU力计算，也不作CPU/GPU性能对比。使用现有 lammps/30Sep2026-p100-gpu 模块加载 LAMMPS 30Sep2026 及配套MPI环境；后处理用 ase/3.29.0 的NumPy/Matplotlib。

计算结束后：

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/neb-images.csv

reference/ 保存已有实算，可把后处理命令的 work 换为 reference 先看配套结果。新计算使用新解压目录，不覆盖自己的既有结果。

1 明确初态、末态与能量模型#

初态在输入中建立。末态文件 final.coords 第一行写本文件包含的原子数，随后每行是 原子ID x y z；不是POSCAR或完整data文件。本例只移动ID为2的粒子：

1
2 1.0 0.0 0.0

final 模式根据初末态线性插值中间图像。所有分区必须保持相同粒子ID、类型、盒子和势函数；真实材料先分别充分优化两个端点，再建立路径，不要让端点未优化的残余力影响能垒。

模型能量为 U=(x²−1)²+5(y−0.5+0.5x²)²+5z²。两个最低点能量均为0，鞍点能量为1。输入用 fix addforce ... energy 同时给出外场力和对应能量，并用 fix_modify field energy yes 把这部分能量计入势能；缺少这一步会让优化所见能量与力不一致。

2 读懂完整输入#

variable image uloop 9
units lj
atom_style atomic
atom_modify map array
boundary f f f
region box block -4 4 -4 4 -4 4 units box
create_box 1 box
mass 1 1
create_atoms 1 single 0 3 0 units box
create_atoms 1 single -1 0 0 units box
group anchor id 1
group mobile id 2
pair_style zero 2.0
pair_coeff * *
fix fixed anchor setforce 0 0 0
variable gap atom y-0.5+0.5*x*x
variable fieldenergy atom (x*x-1)*(x*x-1)+5*v_gap*v_gap+5*z*z
variable fieldx atom -4*x*(x*x-1)-10*x*v_gap
variable fieldy atom -10*v_gap
variable fieldz atom -10*z
fix field mobile addforce v_fieldx v_fieldy v_fieldz energy v_fieldenergy
fix_modify field energy yes
fix band mobile neb 1.0
thermo_modify norm no
timestep 0.02
min_style quickmin
neb 0.0 1.0e-6 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 "FINAL ${image} $(pe:%.16g) $(x[2]:%.16g) $(y[2]:%.16g) $(z[2]:%.16g) $(fx[2]:%.16g) $(fy[2]:%.16g) $(fz[2]:%.16g)" file final-${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

• group mobile id 2 选择参与NEB弹簧连接的原子；固定参照粒子不属于这个组。
• fix band mobile neb 1.0 设置图像间的弹簧强度。弹簧是数值路径约束，不属于材料真实势能。
• min_style quickmin 与 timestep 0.02 设置阻尼动力学最小化；时间步影响收敛，不是本例真实动力学采样时间。
• neb 0.0 1.0e-6 4000 4000 20 final final.coords 把能量停止条件关闭，以力阈值判断；普通阶段和爬山阶段各最多4000步，每20步输出路径。达到上限不等于收敛，要看力和日志。
• variable image uloop 9 为9个分区分别赋予图像编号，输出文件因此不会相互覆盖。

3 MPI分区数怎样对应图像数#

提交整包最简单的是 sbatch submit.slurm。下面展示计算资源已经分配后的9图像原生命令，便于理解作业脚本：

module purge
module load lammps/30Sep2026-p100-gpu
mkdir my-neb
cp inputs/images-9/* my-neb/
cd my-neb
mpirun --bind-to none -np 9 lmp -partition 9x1 -in in.lammps

9x1 是9个分区、每个1个MPI进程。修改成5图像时，要同时改分区数、MPI总进程数、uloop 数和作业资源；不要只把输入文件的9改成5。资料包已经分别配好两个输入，由 run_one.py 根据实际分配的CPU核顺序执行。

4 从输出得到能垒与路径#

全局 log.lammps 的 EBF/EBR 是正向/反向能垒，RD1... 是归一化路径位置，PE1... 是各图像能量；分区日志为 log.lammps.0 等。每个分区有 final-1.dat 等最终数据。

tail -4 work/cases/images-9/log.lammps
cat work/cases/images-9/final-5.dat
head -10 work/analysis/neb-images.csv

图中背景是本例已知解析势能，折线和节点来自实际LAMMPS输出。直线初始路径被优化为弯曲路径；9个图像更细致地描述几何形状。

初始直线路径的最高能量为2.25，优化后两种图像数都得到1.00。这个对称模型的中间图像正好经过鞍点，因而能垒一致；两条折线路径长度仍不同，不能凭这两个能垒相等就宣称路径或真实材料计算已经收敛。

5 换成自己的结构#

把 create_atoms 与外场部分替换为实际结构、势函数和参数；按相同原子ID准备末态坐标，保留 fix neb、合适的最小化器和分区方式。先检查端点及中间图像是否出现原子重叠、错误周期跨越、断键或不合理近邻，再增加图像数、收紧力阈值，并尝试不同初始路径。

NEB这里得到的是势能面路径与能垒，不自动等于有限温度的激活自由能。下一篇：收敛力与独立核对说明为什么中间图像的普通物理力可以不为零。

官方操作说明：neb、fix neb。

继续：Ni空位与副本交换#

Ni107空位NEB流程及端点、力和P100单点检查把建模与真实材料势衔接；NVT/NPT副本交换、gREM广义系综和按状态整理与续算提供四副本完整输入及工具。原生通用命令104/114有成功实例；短模型与低交换率的限制见各篇。