LAMMPS键、键角与二面角势表

当你想确认势参数的形状、平衡位置和量纲，可以先导出势表。本篇继续使用四点机械链，介绍LAMMPS30Sep2026的 bond_write、angle_write 和 dihedral_write。

1. 跑通共同资料包#

下载三种势表、输入和绘图资料包，进入 lammps-topology-tables，在CPU作业中运行：

bash run.sh

如果已经完成前篇，不必在原目录再次运行原生计算，直接读取 raw/。要做一次全新运行则解压到新目录。

2. 导出命令和各自的坐标范围#

本例完整输入里实际执行：

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

键长表有81点，0.7–1.8 Å；键角表91点，0–180°；二面角表120点，0–357°，360°与0°等价，不另存一份。这些是势函数扫描点，不是运行了81/91/120个独立分子动力学任务。

head -n 12 raw/bond.table
head -n 12 raw/angle.table
head -n 12 raw/dihedral.table

四列分别为行号、坐标、能量和该坐标下的负导数。force 列不是某一个原子的笛卡尔力。当前本例键表按kcal/mol/Å核对；两个角度表按对degree求导核对，转换成对radian求导要乘180/π，不能混用。

3. 用本例解析式核对#

本例三种势能分别为：

bond:     E = 50 (r − 1)^2
angle:    E = 20 (θrad − 110π/180)^2
dihedral: E = 0.3 [1 + cos(2φrad)]

harmonic的系数已经按该style约定包含所需因子，不要再自行补一个1/2。键的负导数为−100(r−1)；角度导数先对弧度求，再按角度坐标换算。

head -n 8 results/potential-tables.csv
tail -n 3 raw/angle.table

CSV同时保留 native_energy_kcalmol、analytic_energy_kcalmol、native_derivative、analytic_derivative、导数单位和差值。解析参考是单独计算的列，没有替换原生表里的数字。

左列比较原生能量与解析式，右列比较负导数。曲线来自实际导出数据，不是从网上势函数示意图描画。

4. 注意共线端点和表格行号#

angle_write 通过构造三原子几何并做数值微分得到导数。0°和180°是共线退化几何，本版本该例端点导数与解析值最大相差约 4.43×10⁻⁴ kcal/mol/degree，不适合当作高精度解析导数。

当前输出的180°最后一行，其行号还可能重复为1。读取时按角度列识别位置，不按这个行号拼接；资料包保留原始表，CSV另列明确坐标。需要精确端点导数时使用已确认的解析式，或避开退化端点重新选择分析范围。

这些导出器也不支持所有势类型。带跨项、依赖其他几何量、不能写出系数或某些hybrid/table形式，需要按官方限制另选方法。将表重新用于动力学时，还要分别检查插值误差、表格style的格式及单位，不能只看曲线相似。

source /etc/profile
module load sisl/0.16.4
python -B verify.py

新图在 plots/，原生势表在 raw/，完整数值数组在 results/。官方说明：bond_write、angle_write、dihedral_write。这里的四点机械模型只用于输入和结果解释，不是特定真实分子的验证参数集。