LAMMPS Si共价势与能量力实操

下载完整输入、势文件、实际轨迹和绘图脚本，解压进入目录，先查看已有结果：

bash run.sh

输出 analysis/ 和 figures/。重新计算时先核对脚本里的分区适合自己的账户，再提交：

sbatch submit.slurm

这套脚本申请1个P100和1个CPU，依次运行普通计算路径、适用的GPU/Kokkos路径，再分析自己的结果。已有 work/analysis/figures 时请另解压一份，不要覆盖。直接执行 bash run.sh compute 或文中的 lmp 命令时，必须已经处于分配好的计算资源中。

实际软件为 LAMMPS 30Sep2026，GCC 13.3、OpenMPI 5.0.10、CUDA 12.8；GPU包采用Mixed Precision，Kokkos 5.2.1采用双精度。程序链接现有PLUMED 2.9.4，但本文不使用偏置。数据后处理使用ASE 3.29环境中的NumPy与Matplotlib。普通路径的数值对照在同一P100计算资源的CPU上完成，不是Intel Xeon 8173M性能测试。

1. 先明确这三个势各自做什么#

本篇用64原子金刚石结构Si演示SW、Tersoff和COMB。前两种通过局域环境/键角描述共价键；COMB还需要可变电荷与对应的电荷平衡。三者是三套独立参数化，不能把不同势的绝对能量直接拼到同一物性公式里，也不能用单个晶格点的能量排名判断哪一套普遍更准确。

输入从5.43 Å常规晶胞开始，沿三个方向各重复2次，得到64原子。固定晶胞尺度依次乘0.98、0.99、1、1.01、1.02。这里是固定晶胞模型对照，没有把原子和晶胞全部优化后的曲线冒充严格状态方程。

2. SW完整输入：先评价能量和力#

inputs/sw-force-0-cpu/in.lammps：

newton on
units metal
atom_style atomic
boundary p p p
lattice diamond 5.43
region box block 0 2 0 2 0 2
create_box 1 box
create_atoms 1 box
mass 1 28.0855
pair_style sw
pair_coeff * * potentials/Si.sw Si
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
group probe id 1
displace_atoms probe move 0.02 0 0 units box
compute pea all pe/atom
compute pesum all reduce sum c_pea
timestep 0.001
thermo 100
thermo_style custom step time temp pe ke etotal press vol c_pesum
thermo_modify norm no format float %.16g
run 0
write_dump all custom final.dump id type x y z vx vy vz fx fy fz c_pea modify sort id format float %.16g
print "$(pe:%.16g) $(vol:%.16g) $(c_pesum:%.16g)" file energy.dat screen no

units metal 表示长度Å、能量eV、时间ps。run 0只评价当前坐标的能量与力，不推进动力学；没有时间积分fix的提示在这种单点操作里不代表漏算。示例把1号原子沿x移动0.02 Å，避免在高对称静止点只得到一堆零力。

module purge
module load lammps/30Sep2026-p100-gpu
cd inputs/sw-force-0-cpu

输入中的 potentials/ 由整套运行脚本放到相应工作目录；最简单的原生复算仍是上面的 sbatch submit.slurm。要单独运行这个输入，先把包中的势目录拷到当前目录，再调用：

cp -r ../../potentials .
lmp -in in.lammps
cat energy.dat

pair_coeff 最后的Si是势文件中的元素名，对应LAMMPS的type 1；它不是材料文件名或随机标签。文件前的 * * 表示为所有原子类型建立映射。

3. Tersoff及GPU输入#

下面是同一64原子结构的实际GPU包输入：

package gpu 1 neigh no
newton off
units metal
atom_style atomic
boundary p p p
lattice diamond 5.43
region box block 0 2 0 2 0 2
create_box 1 box
create_atoms 1 box
mass 1 28.0855
pair_style tersoff/gpu
pair_coeff * * potentials/Si.tersoff Si
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
group probe id 1
displace_atoms probe move 0.02 0 0 units box
compute pea all pe/atom
compute pesum all reduce sum c_pea
timestep 0.001
thermo 100
thermo_style custom step time temp pe ke etotal press vol c_pesum
thermo_modify norm no format float %.16g
run 0
write_dump all custom final.dump id type x y z vx vy vz fx fy fz c_pea modify sort id format float %.16g
print "$(pe:%.16g) $(vol:%.16g) $(c_pesum:%.16g)" file energy.dat screen no

这里 package gpu 1 neigh no、newton off 和 tersoff/gpu 配套使用。GPU包中支持某个势，不能推出所有势都有 /gpu 变体。

Kokkos则使用其单独的模块和输入，包含 atom_style atomic/kk、相应 /kk 势以及 run_style verlet/kk。例如在已分配GPU的资源中：

module purge
module load lammps/30Sep2026-p100-kokkos
cd inputs/tersoff-force-0-kk
cp -r ../../potentials .
lmp -k on g 1 -in in.lammps

输入已经显式选择后端，避免只加载模块却仍执行普通势。更多包与势的查找方法见实际style查询。

4. COMB：电荷属性与电荷平衡#

本例COMB使用 ffield.comb 中的Si参数；不要把COMB3的参数文件交给COMB读取。

newton on
units metal
atom_style charge
boundary p p p
lattice diamond 5.43
region box block 0 2 0 2 0 2
create_box 1 box
create_atoms 1 box
mass 1 28.0855
set group all charge 0
pair_style comb
pair_coeff * * potentials/ffield.comb Si
fix charges all qeq/comb 1 1e-8
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
group probe id 1
displace_atoms probe move 0.02 0 0 units box
compute pea all pe/atom
compute pesum all reduce sum c_pea
compute qatom all property/atom q
compute qsum all reduce sum c_qatom
timestep 0.001
thermo 100
thermo_style custom step time temp pe ke etotal press vol c_pesum c_qsum
thermo_modify norm no format float %.16g
run 0
write_dump all custom final.dump id type x y z q vx vy vz fx fy fz c_pea modify sort id format float %.16g
print "$(pe:%.16g) $(vol:%.16g) $(c_pesum:%.16g) $(c_qsum:%.16g)" file energy.dat screen no

关键区别是 atom_style charge、初始电荷和 fix qeq/comb。要统计总电荷，先用 compute property/atom q 取得每原子电荷，再将该向量求和。compute pe/atom 的和可与全胞势能对照。

5. 实际得到什么#

尺度为1、未移动原子的原生能量为：

energy.dat保存的是全胞总势能、体积和原子势能总和；表中每原子值是在后处理中除以64。final.dump保存原子坐标和力；完整包还有三个位置偏移的单点数据，用能量中心差分比较x方向的力。

SW另外按势文件中的两体项和三体项直接重建能量，并从重建的能量取差分；本例1号原子的非零力与原生结果一致。普通路径与Kokkos的有限差分误差约10⁻⁹ eV/Å量级；GPU混合精度以同坐标力对照为主，不用很小位移下的两个混合精度能量相减硬凑高精度导数。

6. 再看短时NVE#

下载包的 *-md-0-* 输入在初始100 K下演示1000步NVE。SW/Tersoff/COMB的时间步是0.001 ps。NVE并不把温度锁定在100 K；初始动能会与势能交换，非平衡小体系的瞬时温度也会变化。

本例COMB在1 ps内的最大总能变化为约0.00280735 eV/atom，明显大于这里的SW/Tersoff。它只用来演示如何读能量漂移，不能据此直接开展定量长MD；正式任务应继续收敛时间步、电荷求解条件，先优化结构并完成相应平衡。本文不据这些短轨迹报告热导或反应速率。

接着可看碳体系AIREBO和ReaxFF/COMB3可变电荷。官方输入说明：SW、Tersoff、COMB。