LAMMPS键与角势能表实操
原生bond_write、angle_write、dihedral_write输出，实际曲线与解析能量斜率检查。
阅读本篇目录
设置分子力场后，先把势能随键长、键角和二面角的变化画出来，有助于核对参数和单位。本篇用一个四站点教学结构，实际输出三张原生势能表并检查能量与斜率。
下载六组完整输入、实际原生输出、势能表、JSON与绘图脚本，解压进入run.sh目录。
1 提交与目录#
sbatch submit.slurm
模板为核时P100分区1GPU1CPU，顺序完成六组很小的原生工作流，再自动提取与绘图。这里调用lammps/30Sep2026-p100-gpu，但未启用加速后缀，演示的是其CPU执行路径；这不是GPU性能测试。构建为GNU C++13.3.0/C++17、Open MPI5.0.10，分析使用ASE3.29.0中的NumPy与Matplotlib环境。
套餐用户用yskit生成所属分区脚本，在分配资源中执行bash run.sh；不要在登录节点运行计算。inputs是干净输入，reference是本文实际输出，work是新结果；已有work重新解压。
只想手动练习一个输入时，在已分配资源中复制对应目录，再进入它：
cp -r inputs/bonded-tables my-tables
cd my-tables
module load lammps/30Sep2026-p100-gpu
export OMP_NUM_THREADS=1
lmp -in in.lammps
单独目录不经过整套分析脚本；完整绘图与检查使用前面的提交方式。
2 完整输入#
units real
atom_style full
boundary f f f
region box block -10 20 -10 20 -10 20
create_box 1 box bond/types 1 angle/types 1 dihedral/types 1 improper/types 1 extra/bond/per/atom 4 extra/angle/per/atom 4 extra/dihedral/per/atom 2 extra/improper/per/atom 2 extra/special/per/atom 8
mass 1 12.0
pair_style zero 8.0
pair_coeff * *
bond_style harmonic
bond_coeff 1 100.0 1.5
angle_style harmonic
angle_coeff 1 50.0 109.47
dihedral_style opls
dihedral_coeff 1 1.0 0.5 0.2 0.0
improper_style harmonic
improper_coeff 1 10.0 0.0
create_atoms 1 single 0 0 0
create_atoms 1 single 1.5 0 0
create_atoms 1 single 2.0 1.4 0
create_atoms 1 single 2.7 1.7 1.2
set atom * mol 1
create_bonds single/bond 1 1 2
create_bonds single/bond 1 2 3
create_bonds single/bond 1 3 4
create_bonds single/angle 1 1 2 3
create_bonds single/angle 1 2 3 4
create_bonds single/dihedral 1 1 2 3 4
create_bonds single/improper 1 1 2 3 4
thermo_style custom step atoms pe ebond eangle edihed eimp
thermo_modify norm no
run 0
bond_write 1 401 0.5 2.5 bond.table HarmonicBond
angle_write 1 361 angle.table HarmonicAngle
dihedral_write 1 720 dihedral.table OPLSTorsion
write_data structure.data
write_dump all custom atoms.dump id type x y z fx fy fz modify sort id
units real使用Å和kcal/mol。harmonic键为E=K(r−r0)²，K=100、r0=1.5Å；harmonic键角为E=K(θ−θ0)²，内部角度使用弧度，输入平衡角为109.47度。不要再给这两个K额外乘一个1/2。
本例还显式创建了3条键、2个键角、1个二面角和1个improper。create_bonds single/dihedral不会顺带创建前三条键，必须分别定义；所有力场系数在写表前设置齐全。
3 三条写表命令#
bond_write 1 401 0.5 2.5 bond.table HarmonicBond
angle_write 1 361 angle.table HarmonicAngle
dihedral_write 1 720 dihedral.table OPLSTorsion
第一条对1号键类型在0.5–2.5Å取401点；第二条取0–180度共361点；第三条对0–360度取720点，360度与0度重复，所以最后一个实际角度为359.5度。
head -n 12 work/cases/bonded-tables/bond.table
head -n 12 work/cases/bonded-tables/angle.table
head -n 12 work/cases/bonded-tables/dihedral.table
表格每个数据行依次为序号、坐标、能量、负能量斜率。键表最后一列按Å求导；本次angle_write和dihedral_write原生角度表最后一列，实际与按“度”求导的负能量斜率一致。分析脚本同时代入解析表达式核对，不能把这一列直接当作某个原子的笛卡尔力，也不要与按弧度表示的力矩混用。
angle_write在0和180度附近构造的三原子几何有端点数值敏感性，脚本保留完整原生表，但解析斜率检查选取1–179度。二面角采用OPLS余弦表达式，本例K1/K2/K3/K4=1/0.5/0.2/0，单位kcal/mol。
4 实际曲线#
三张表共1482个原生采样点。键长最低点在1.5Å；键角最低点接近109.47度；这组OPLS系数的最低点在180度。曲线只说明输入势函数的形状，不表示相应真实分子的能垒或热力学性质。
5 换用自己的参数#
先按所用力场核对units和K的定义，再修改style与coeff；确认写表命令支持该样式。angle_write和dihedral_write不是任意混合/多几何依赖势都能直接使用，例如hybrid、已有table及某些交叉项有专门限制。新表应先核对解析形式或数值导数，再用于后续工作。
提取与重画#
module load ase/3.29.0
python analyze.py work
python plot_results.py work
python check_results.py work
ls work/figures
cat work/analysis/coulomb.csv
分析检查原生日志是否完成、势能表与解析表达式、拓扑数量和实际电荷对的能量/力。查看本文参考数据时把work改成reference。四站点模型是人为设置的教学结构，参数用于解释命令，没有拟合为某个真实分子的完整力场。
参数与适用条件：bond_write、angle_write、dihedral_write。