LAMMPS自旋优化与力矩核对

普通NEB优化原子位置，自旋NEB优化磁矩方向。本例原子始终固定在原点，只让一个单位磁矩在外场和单轴各向异性能量下翻转，实际比较 spin、spin/cg、spin/lbfgs 三种优化器与7个路径图像。

这是一个有解析能量基准的经典单自旋教学模型，不是某种铁、镍或其他材料的DFT翻转能垒。它便于先弄清文件格式、固定模长、路径、力矩和单位，再换成合适的真实磁性模型。

下载完整输入、21个实际图像、角度差分数据与图片，在个人目录执行：

mkdir spin-neb
cd spin-neb
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-spin-neb-wave82-2c81dd6c57bd.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

脚本申请7个CPU核、一个P100资源，每个图像1个MPI进程，顺序完成三种优化器和单点检查。实际GNEB使用CPU/MPI，没有GPU力计算或速度对比。模块为 lammps/30Sep2026-p100-gpu（GCC13.3、OpenMPI5.0.10、CUDA12.8）；分析用 ase/3.29.0 的NumPy/Matplotlib。

完成后：

module purge
module load ase/3.29.0
python analyze.py work
python export_spin_path.py work
python plot_results.py work
head -10 work/analysis/spin-neb.csv

已有实算在 reference/，可将后处理参数 work 换为 reference 先看结果。新运行换一个解压目录，不覆盖自己的work。

1 三种优化器怎样设置#

spin 使用阻尼自旋动力学；spin/cg 和 spin/lbfgs 提供其他自旋优化方案。GNEB中的图像约束不是普通保守力，后两者在这个任务中关闭线搜索：

min_style spin/cg
min_modify discrete_factor 20
min_modify line spin_none

把第一行改为 spin/lbfgs 就是另一例。普通 spin 输入如下，完整文件可直接从资料包使用：

variable image uloop 7
units metal
atom_style spin
atom_modify map array
region box block -5 5 -5 5 -5 5 units box
create_box 1 box
create_atoms 1 single 0 0 0 units box
mass 1 1
set atom 1 spin 1 0.05789010000000001 0 0.9983229619326554
pair_style zero 3.0
pair_coeff * *
fix external all precession/spin zeeman 10 1 0 0 anisotropy .01 0 0 1
fix_modify external energy yes
fix band all neb/spin 1.0
compute magnetic all property/atom sp spx spy spz fmx fmy fmz
thermo_modify norm no format float %.16g
min_style spin
min_modify discrete_factor 20
neb/spin 0.0 1.0e-7 4000 4000 20 each spin-${image}.coords
write_dump all custom neb-${image}.dump id type x y z c_magnetic[*] 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 c_magnetic[*] modify sort id format float %.16g

discrete_factor 控制自适应步长的离散程度，不是实际MD输出时间。所有方案都只优化自旋方向，不移动晶格原子。

2 力矩不等于磁进动向量本身#

compute property/atom 导出的fmx/fmy/fmz在这里是磁进动向量ω。实际驱动方向改变的量包含 s×ω。即使端点ω不为0，只要它与s平行，力矩仍为0。

当前30Sep2026、metal单位下，最小化停止条件使用 ħ·||s×ω|| 的角能量力矩，单位是eV尺度。程序实际ħ为 4.135667403e-3/(2π) eV·ps。本例阈值是 1e-7 eV。

全局日志有不同的原生量：MaxReplicaTorque 是未乘ħ的频率力矩平方和，MaxAtomTorque 是频率力矩最大分量。它们的范数与单位不同，不能把原始列值直接与eV阈值比较。配套脚本从向量重新计算，并将平方量开方后乘ħ核对。

tail -4 work/cases/spin/log.lammps
head -12 work/cases/spin/neb-2.dump
head -12 work/cases/spin/physical-2.dump

neb-编号.dump 是NEB约束尚在时的磁向量；随后 unfix band、run 0 重算普通磁向量，写到 physical-编号.dump。中间图像的普通物理力矩可沿路径保持非零，不能直接用它代替NEB修正力矩判断收敛。

3 用真实旋转能量差分核对#

单点控制取实际优化路径的第2图像，分别绕x/y/z轴旋转±1e-5 rad，用原生LAMMPS重新计算六个能量。负旋转能量导数与 ħ(s×ω) 的对应分量一致。

cat work/cases/control/center.dat
cat work/cases/control/rotation-energy.dat
module load ase/3.29.0
python - <<'PY'
from pathlib import Path
rows=[s.split() for s in Path('work/cases/control/rotation-energy.dat').read_text().splitlines()]
for axis in ['x','y','z']:
    pair=[r for r in rows if r[1]==axis]
    torque=-(float(pair[1][3])-float(pair[0][3]))/(2e-5)
    print(axis,'angular-energy torque / eV:',torque)
PY

左图转换到与停止阈值相同的eV尺度；右图是实际单点旋转能量差分与磁向量的独立对照，没有用原子平移力来替代自旋力矩。

4 导出方向向量与查看数据#

python export_spin_path.py work
head -9 work/analysis/spin-spin-path.xyz
head -10 work/analysis/spin-neb.csv
python check_results.py work

扩展XYZ有固定位置和spin方向属性，物种X是这个教学磁矩的占位标签。可用配套绘图脚本直接画方向与能量；这些帧是路径图像，不是MD时间，也不是原子迁移轨迹。

更改磁矩模长、场强、各向异性或单位后，同步更新解析基准和检查脚本。不要把本例单自旋模型的能垒或误差当作其他材料的结果。返回完整流程。

优化器约定见官方min_style spin说明，GNEB输入见官方neb/spin说明。