LAMMPS 镍位错建模与松弛
用Atomsk构建4320原子刃型和螺型位错，经GPU能量最小化，再观察实际核心变化与残余力。
阅读本篇目录
本篇把 Atomsk 建模 → LAMMPS 松弛 → 结构与数据输出 串起来。下一篇继续用 OVITO 提取位错线和 Burgers 矢量。全部结构与结果来自鸢算科技集群实际运行。
本例是有限尺寸、固定外边界的经典 EAM 镍模型，不涉及第一性原理磁性计算。它帮助准备位错输入和检查松弛，不能把两份有限胞能量之差直接当作位错核心能。
1. 下载与建立结构#
python3 -m zipfile -e lammps-Ni-dislocation-wave542-05f1cf979c2d.zip Ni-dislocation
cd Ni-dislocation
bash build.sh
这一步在登录节点用一个 CPU 线程运行 Atomsk，生成四个 4320 原子结构：刃型方向的完美晶体与位错胞，螺型方向的完美晶体与位错胞。此时还没有提交 LAMMPS 松弛。
head -n 18 work/models/edge-dislocation.lmp
head -n 18 work/models/edge-dislocation.data
.lmp 保留原生 Atomsk 输出。.data 只把非周期 x、y 导入盒边界扩到足以容纳位移后的原子；原子坐标、ID、类型及 z 周期均不改变。对自由表面模型，这是容纳位移后的边缘，不是再次施加应变。
2. 晶向和 Burgers 矢量怎样选#
FCC Ni 晶格常数取 a=3.52 Å，完整 1/2⟨110⟩ 位错的 Burgers 矢量大小为 a/√2=2.48901587 Å。两个晶向设置为：
| 模型 | x 轴 | y 轴 | z 轴/位错线 | 复制数 |
|---|---|---|---|---|
| 刃型 | [110] | [−111] | [1−12] | 18×10×4 |
| 螺型 | [1−12] | [1−1−1] | [110] | 12×10×6 |
刃型的 Burgers 矢量沿 x，垂直于 z 位错线；螺型沿 z，与位错线平行。各自的完美晶体使用相同晶向、复制数和原子数。
建模脚本实际调用的刃型命令包括：
atomsk --create fcc 3.52 Ni orient [110] [-111] [1-12] -duplicate 18 10 4 edge-perfect.lmp
atomsk edge-perfect.lmp -dislocation 0.51*box 0.51*box edge Z Y 2.489015869776647 0.31 edge-dislocation.cfg
atomsk edge-dislocation.cfg edge-dislocation.lmp
位错位置取 x、y 盒长的 0.51，避免直接落在原子上。0.31 是构造初始各向同性弹性位移场使用的泊松比，不是这次由 Ni 势计算出的弹性常数。螺型命令使用 screw Z Y 和相同的 |b|，不传泊松比。edge 保持原子数，不等于插入半原子面的 edge_add。
图中截取每个周期胞 z 范围的前四分之一并投影到 xy，优先显示 HCP/其他局部环境以便观察核心；原始三维坐标全部保留。颜色是 OVITO DXA 的局部结构分类，表面或核心的 Other 标签不自动表示一种新体相。
3. 提交 GPU 松弛#
bash run.sh
完整对照流程需要 P100 和 V100 16G 分区权限。脚本经普通账号 Slurm 提交八次松弛：四个结构分别使用两个 Kokkos GPU 版本。每个任务一张 GPU、两个申请的 CPU 核心、一个 MPI 进程和一个 OpenMP 线程。随后还有一个同几何单点与能量有限差分任务；都经 Slurm 执行。
work/SUBMISSIONS.json 与 work/CROSS-SUBMISSION.json 保存作业号。已有提交记录时不会重复提交；终端中断后先核对作业状态。想重新开始，请使用新解压目录。
cat work/SUBMISSIONS.json
head -n 5 tables/Ni-minimization.csv
4. 完整输入与固定边界#
下面是刃型位错的输入形式。fixed-group.in 由脚本根据建模阶段确定的原子 ID 写出，固定原子是实际 x、y 外边缘 5 Å 范围，z 方向仍周期。
units metal
atom_style atomic
boundary s s p
read_data ../models/edge-dislocation.data
pair_style eam
pair_coeff * * ../Ni_u3.eam
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
include fixed-group.in
group mobile subtract all fixed
fix raw all store/force
fix hold fixed setforce 0 0 0
compute virial all stress/atom NULL pair
variable fm atom sqrt(fx*fx+fy*fy+fz*fz)
compute fm mobile reduce max v_fm
thermo 100
thermo_style custom step atoms pe c_fm
thermo_modify norm no format float %.16g
dump d all custom 1 initial.dump id type x y z fx fy fz f_raw[1] f_raw[2] f_raw[3] c_virial[1] c_virial[2] c_virial[3] c_virial[4] c_virial[5] c_virial[6]
dump_modify d sort id format float %.16g
run 0
undump d
print "$(pe:%.16g) $(c_fm:%.16g)" file initial-properties.dat screen no
min_style cg
min_modify dmax 0.05
minimize 0 0.002 5000 50000
dump d all custom 1 final.dump id type x y z fx fy fz f_raw[1] f_raw[2] f_raw[3] c_virial[1] c_virial[2] c_virial[3] c_virial[4] c_virial[5] c_virial[6]
dump_modify d sort id format float %.16g
run 0
undump d
print "$(pe:%.16g) $(c_fm:%.16g)" file final-properties.dat screen no
write_data final.data nocoeff
minimize 0 0.002 5000 50000 不使用能量变化停止条件，目标是力二范数达到 0.002 eV/Å，最多 5000 次迭代和 50000 次力评估。力二范数与单原子最大力不同；本包同时导出二者。这里只进行静态最小化，没有温度、恒温器或有限温度 MD。
5. 实际结果怎么看#
P100 刃型位错胞能量从 −18830.463253 eV 降到 −18840.691169 eV，降低约 10.227916 eV；螺型胞降低约 4.610919 eV。末态 mobile 原子的最大力约为 1.4×10⁻⁴、1.5×10⁻⁴ eV/Å。两种晶向的表面与固定边界不同，不能通过比较它们的总能量直接判断刃型和螺型哪种核心更稳定。
结构位于 work/*/final.data、final.dump，初末能量在各目录的 initial-properties.dat 与 final-properties.dat。CSV 保留完整数值；图中完美对照很小的能量变化没有被当作材料常数。
环境：Atomsk 0.13.1（GFortran 13.3）；LAMMPS 30Sep2026 Kokkos（GCC 13.3、CUDA 12.8、OpenMPI 5.0.10、Kokkos 5.2.1）；OVITO 3.16.1。Ni_u3.eam 为 LAMMPS 附带的 Foiles 等 Ni 势，出处见文件首行，Phys. Rev. B 33, 7983 (1986)。建模参考Atomsk 官方位错教程，最小化参数参考LAMMPS 文档。