LAMMPS分子拓扑与读写实操
用四点机械链操作分子模板、类型标签、添加和删除键，核对data、restart和原子编号重排后的能量与力。
阅读本篇目录
在LAMMPS中，原子坐标和分子拓扑是两部分：两个点靠得近，并不自动表示存在一条键。这里用一个明确给定的四点机械链，练习拓扑修改、模板导出和完整状态读写。它不是某个真实分子的已标定化学力场，不用于解释反应动力学。
1. 下载并运行#
下载输入、结果参考与绘图包，进入 lammps-topology-tables，在已分配的CPU作业内运行：
bash run.sh
也可 sbatch job.slurm 提交；CPU分区按56核整节点分配和计费，本包仅用一个线程。软件为现有 lammps/30Sep2026-8173m，依赖oneAPI2024.2和PLUMED2.9.4；本例不调用PLUMED，也不使用GPU。Python绘图由 sisl/0.16.4 环境提供。
2. 准备坐标与拓扑#
inputs/chain.data 包含4个原子、3条键、2个键角、1个二面角和1个improper。所有点质量12.01、电荷0，units real 下长度Å、能量kcal/mol。非键相互作用为zero，只演示这里明确指定的有键项。
cat inputs/chain.data
cat inputs/chain.mol
chain.data 是包含盒子、质量和原子属性的数据文件；chain.mol 是可复用的分子模板，包含坐标、类型、电荷和连接关系。二者格式不同，不能互相改扩展名代替。
3. 完整输入和操作顺序#
下面是资料包实际执行的 input.lmp。先记录原始状态，再添一条1–4键、删除所有指定类型的键，最后分别从restart、data恢复并重排ID。
展开完整LAMMPS输入
units real
atom_style full
boundary f f f
read_data chain.data extra/bond/per/atom 2 extra/special/per/atom 4
pair_style zero 8.0
pair_coeff * *
bond_style harmonic
bond_coeff 1 50.0 1.0
angle_style harmonic
angle_coeff 1 20.0 110.0
dihedral_style harmonic
dihedral_coeff 1 0.3 1 2
improper_style harmonic
improper_coeff 1 2.0 90.0
labelmap atom 1 C
labelmap bond 1 spring
molecule chain chain.mol
write_molecule chain chain-copy.mol
write_molecule chain chain-copy.json
thermo 1
thermo_style custom step atoms pe ebond eangle edihed eimp
thermo_modify norm no format float %.15g
run 0
print "STATE base $(pe:%.15g) $(ebond:%.15g) $(eangle:%.15g) $(edihed:%.15g) $(eimp:%.15g)"
write_dump all custom base.dump id type x y z fx fy fz modify sort id format float %.15g
write_data chain-rewrite.data
write_restart base.restart
bond_write 1 81 0.7 1.8 bond.table SPRING
angle_write 1 91 angle.table ANGLE
dihedral_write 1 120 dihedral.table TORSION
create_bonds single/bond 1 1 4
run 0
print "STATE added $(pe:%.15g) $(ebond:%.15g) $(eangle:%.15g) $(edihed:%.15g) $(eimp:%.15g)"
write_dump all custom added.dump id type x y z fx fy fz modify sort id format float %.15g
delete_bonds all bond 1 remove special
run 0
print "STATE deleted $(pe:%.15g) $(ebond:%.15g) $(eangle:%.15g) $(edihed:%.15g) $(eimp:%.15g)"
write_dump all custom deleted.dump id type x y z fx fy fz modify sort id format float %.15g
clear
read_restart base.restart
thermo_style custom step atoms pe ebond eangle edihed eimp
thermo_modify norm no format float %.15g
run 0
print "STATE restarted $(pe:%.15g) $(ebond:%.15g) $(eangle:%.15g) $(edihed:%.15g) $(eimp:%.15g)"
write_dump all custom restarted.dump id type x y z fx fy fz modify sort id format float %.15g
clear
units real
atom_style full
boundary f f f
pair_style zero 8.0
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style harmonic
read_data chain-rewrite.data
thermo_style custom step atoms pe ebond eangle edihed eimp
thermo_modify norm no format float %.15g
run 0
print "STATE reread $(pe:%.15g) $(ebond:%.15g) $(eangle:%.15g) $(edihed:%.15g) $(eimp:%.15g)"
write_dump all custom reread.dump id type x y z fx fy fz modify sort id format float %.15g
reset_atoms id sort yes
run 0
print "STATE reset_ids $(pe:%.15g) $(ebond:%.15g) $(eangle:%.15g) $(edihed:%.15g) $(eimp:%.15g)"
write_dump all custom reset_ids.dump id type x y z fx fy fz modify sort id format float %.15g
labelmap 把数值类型1映射成 C 和 spring 标签，它不改变质量和势参数。molecule 读入模板，write_molecule 可导出原生格式或JSON；JSON格式输出是此版本的功能，不能直接套用到不含该命令的旧版。
delete_bonds all bond 1 remove special 只针对这里指定的键类型。 本例删除键后，键角、二面角和improper仍然存在，因此势能不会变成零。实际要删除哪些拓扑，必须按建模目的分别处理，不能把此命令直接当作真实化学断键模型。
4. 查看六种状态与完整力#
grep '^STATE ' raw/native.out
head -n 14 raw/base.dump
head -n 20 raw/chain-copy.json
原始链势能为 2.89896473 kcal/mol，其中键0.93325507、键角0.49880996、二面角0.28918558、improper1.17771413。添加1–4键提高了键能；删除所有键后，其他三类项仍保留。这里只是给定模型的能量记账。
结构上的箭头来自实际原子力，显示长度经过统一缩放。右侧添加键的能量较大，其他项的精确值可从终端逐项读取。
write_restart/read_restart 保存并恢复本例所需的原生状态；write_data/read_data 可转移结构、拓扑和本例势参数。生产体系里的fix、外部模型、随机状态等是否需要重新定义，要按对应命令说明判断，不能仅凭本小例推断一切都会自动恢复。
reset_atoms id 可能改变原子编号。比较前后原子力时应按同一原子或坐标映射对应，不能只按文件行号硬减。
5. 重新提取与绘图#
source /etc/profile
module load sisl/0.16.4
python -B verify.py
本例从restart、data和重排编号后得到的完整分项能量与力一致。results/full-fields.npz 保存全数组，新图在 plots/；原生输出在 raw/。继续阅读有键势表与导数。命令用法见官方的molecule、write_molecule、delete_bonds。