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 先看配套结果。新计算使用新解压目录，不覆盖自己的既有结果。

先看输入与最终文件#

完整9图像输入如下，5图像版本只改变图像数及配套启动资源：

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

1 分清两种力#

在NEB优化中，普通势能梯度沿路径方向的分量会被处理，再加入图像间弹簧力。最小化器收敛的是这个NEB修正后的力；中间图像本来就在斜坡上，其普通物理力可以沿路径方向保持非零。端点应是稳定态，爬山图像应靠近鞍点。

本例先输出 neb-编号.dump 保存修正力，然后 unfix band、run 0 重新计算，再输出 physical-编号.dump 和 final-编号.dat 保存普通力。位置没有在这个零步计算中变化，因此可以直接对照。

head -12 work/cases/images-9/neb-3.dump
head -12 work/cases/images-9/physical-3.dump
cat work/cases/images-9/final-3.dat

9图像例中最大普通力约1.3681，而最终最大NEB修正力约 7.12e-7，达到输入 1e-6 阈值。不要把普通力1.3681直接当作NEB没有收敛，也不要只检查中间鞍点恰好为零就忽略其他图像。

左图取原生全局日志的 MaxReplicaForce，右图取移除NEB前后实际dump中的力；比较的是两种定义不同的量。

2 用能量差分独立核对外场力#

记 g=y−0.5+0.5x²，本例普通力应为 Fx=−4x(x²−1)−10xg、Fy=−10g、Fz=−10z。分析工具同时用每个图像位置前后 1e-5 的能量中心差分核对，而不是只把输入中的力公式重算一遍。

下面这段可以直接在整包目录执行，使用真实输出CSV：

module load ase/3.29.0
python - <<'PY'
import csv
import numpy as np
def energy(p):
    x,y,z=p
    return (x*x-1)**2+5*(y-.5+.5*x*x)**2+5*z*z
errors=[]
for row in csv.DictReader(open('work/analysis/neb-images.csv')):
    p=np.array([float(row[k]) for k in ['x','y','z']])
    actual=np.array([float(row[k]) for k in ['physical_fx','physical_fy','physical_fz']])
    h=1e-5
    fd=np.array([-(energy(p+np.eye(3)[i]*h)-energy(p-np.eye(3)[i]*h))/(2*h) for i in range(3)])
    errors.append(max(abs(fd-actual)))
print('Maximum finite-difference force error:',max(errors))
PY

实际最大误差约 9.0e-10。同时，原生势能与解析模型逐图像一致。真实材料势函数也应在可行时做独立能量/力检查，尤其使用自定义外场或插件时。

3 端点、鞍点和图像数#

本例端点坐标保持 (-1,0,0) 与 (1,0,0)，能量均为0；5/9图像中心均达到 (0,0.5,0)，能量为1。该点解析Hessian特征值为 −4,10,10，有且只有一个负曲率方向，与一阶鞍点一致。

两组都通过力阈值、端点、对称性、等间距、能垒和能量梯度核对。5图像折线路径长度2.30289，9图像2.32882，提醒我们图像数仍影响路径几何。这个例子不提供真实材料的图像数收敛结论。

4 用工具查看结构与数据#

export_path.py 导出 images-5-path.xyz 和 images-9-path.xyz，帧顺序是路径图像顺序，不是动力学时间。可用OVITO打开查看，但本篇配图由配套Matplotlib脚本直接读取实际数据生成，不依赖手工截屏。

python export_path.py work
head -12 work/analysis/images-9-path.xyz
python plot_results.py work
python check_results.py work

将自己的计算整理成相同输出格式后，可以继续复用CSV、XYZ和绘图思路；解析双阱能量/力公式只适用于本例，不能原封不动用来检查其他势函数。

返回NEB搭建流程。外场能量要求参照官方fix addforce说明，路径修正力参照官方fix neb说明。

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

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