Moltemplate自动生成键角与二面角

建立有键分子时，不一定要手写每一个角和二面角。Moltemplate可以根据已有的连接关系，结合按类型定义的规则，生成相应相互作用。本篇用四个珠子的A–B–A–B链，演示从 Data Bond List 到3条键、2个角和1个二面角，再让LAMMPS实际读回能量和力。

这是自建教学势，A/B不是实际元素，质量为1/2，采用LAMMPS约化LJ单位。不要把坐标直接当成Å，或将结果解释成某种真实分子的性质。

1. 下载、运行与结果位置#

下载完整LT、LAMMPS输入、真实输出与绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd moltemplate-graph-types-charges

先查看已经计算的数据：

bash show-reference.sh
cat graph-energy-forces.csv
cat charge-fields.csv

重新计算时，在新解压目录提交：

sbatch job.slurm

CPU模板按56核整节点申请，本例短任务实际单进程执行；分区维护时等待可用资源。在已分配的资源中也可以用 bash run.sh。脚本依次完成：生成拓扑和参数 → 对最终键表赋一次电荷 → 含电荷、仅读DATA、改变平衡角三个 run 0 → 导出表格和图片。重新执行整套流程请换新目录，避免覆盖已有 work/。

使用现有Moltemplate2.22.5/Python3.11与LAMMPS30Sep2026 CPU版，后者以Intel oneAPI2024.2、Intel MPI和MKL构建。分析绘图通过 sisl/0.16.4 的Python完成，不用修改公共软件。

2. 四个珠子与三条连接#

inputs/molecule.lt 写坐标、类型和连接：

import "ff.lt"
Toy inherits ToyFF {
 create_var {$mol}
 write("Data Atoms") {
  $atom:a1 $mol @atom:A 0.0 0.0 0.0 0.0
  $atom:a2 $mol @atom:B 0.0 1.3 0.0 0.0
  $atom:a3 $mol @atom:A 0.0 1.3 1.3 0.0
  $atom:a4 $mol @atom:B 0.0 2.6 1.3 0.4
 }
 write("Data Bond List") {
  $bond:b1 $atom:a1 $atom:a2
  $bond:b2 $atom:a2 $atom:a3
  $bond:b3 $atom:a3 $atom:a4
 }
}

Data Bond List 每行是键变量、原子1、原子2，此时不直接写数字键类型。$atom:a1 是原子实例变量，@atom:A 是静态类型变量；它们在构建时得到相应编号。不要将分子内的变量名和最终DATA中的数字ID混为一谈。

顶层 inputs/system.lt：

import "molecule.lt"
m=new Toy
write_once("Data Boundary") {
 -10 20 xlo xhi
 -10 20 ylo yhi
 -10 20 zlo zhi
}

3. 按类型生成相互作用#

完整 inputs/ff.lt：

ToyFF {
 write_once("In Init") {
  units lj
  atom_style full
  pair_style zero 4.0
  bond_style harmonic
  angle_style harmonic
  dihedral_style opls
 }
 write_once("Data Masses") {
  @atom:A 1.0
  @atom:B 2.0
 }
 write_once("In Settings") {
  pair_coeff * *
  bond_coeff @bond:link 100.0 1.3
  angle_coeff @angle:turn 5.0 90.0
  dihedral_coeff @dihedral:torsion 0.2 0.1 0.05 0.0
 }
 write_once("Data Bonds By Type") {
  @bond:link @atom:* @atom:*
 }
 write_once("Data Angles By Type") {
  @angle:turn @atom:* @atom:* @atom:*
 }
 write_once("Data Dihedrals By Type") {
  @dihedral:torsion @atom:* @atom:* @atom:* @atom:*
 }
 write_once("Data Charge By Bond") {
  @atom:A @atom:B 0.2 -0.2
 }
}

这里 Data Bonds By Type 为连接匹配link类型；Data Angles By Type 匹配连续三个原子；Data Dihedrals By Type 匹配连续四个原子。星号用于本例通配匹配。真实力场通常需要更具体的类型组合和参数，不能直接把这组通配规则套到所有化学体系。

这条链实际生成的连接为：键1–2、2–3、3–4；角1–2–3和2–3–4；二面角1–2–3–4。并不是任意三原子或四原子组合都生成一项，需要满足连接图与匹配规则。

4. 单独生成并查看文件#

source ./env.sh
mkdir customer-build
cp inputs/*.lt customer-build/
cd customer-build
moltemplate.sh system.lt
grep -E 'atoms$|bonds$|angles$|dihedrals$' system.data
cat system.in.settings

本例输出4 atoms、3 bonds、2 angles、1 dihedrals。system.in.settings 写实际系数；只看原子数正确还不够，连接数、参数和电荷也要一起检查。按键赋电荷的具体处理见电荷与生成文件教程。

5. 赋电荷后让LAMMPS读回#

source ./env.sh
cd customer-build
bash ../assign-charges.sh
cp ../inputs/check0-with-charges.in .
lmp -var dumpfile customer-forces.dump -in check0-with-charges.in -log customer-check0.log

其中完整的读回输入是：

include system.in.init
read_data system.data
include system.in.settings
include system.in.charges
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes
thermo 1
thermo_style custom step atoms pe ke etotal ebond eangle edihed
thermo_modify norm no format float %.16g
run 0
write_dump all custom ${dumpfile} id mol type q x y z fx fy fz modify sort id format float %.16g

run 0 只计算当前位置的能量和力，不积分推进，也不优化结构。read_data 之后读参数，再读 system.in.charges，最后计算。少了电荷文件，运行中的电荷可能仍是初始DATA的0。

6. 改平衡角，观察非零角力#

原几何的两个角均为90°，最初平衡角也是90°，因此角能为0。配套的 check0-angle-shift.in 在相同坐标上加入：

angle_coeff * 5.0 105.0

在自己的副本中运行这一对照：

source ./env.sh
cd customer-build
cp ../inputs/check0-angle-shift.in .
lmp -var dumpfile customer-angle-forces.dump -in check0-angle-shift.in -log customer-angle.log

角势为 K(θ−θ0)^2，计算时角度差使用弧度；两个15°偏差合计得到约0.685389195的角能。几何没有改变，力发生变化，便于区分“模型位置”和“势函数平衡参数”。

7. 本例实际数据#

能量均为整体系约化单位。程序还逐分量对照了能量有限差分与实际LAMMPS力，最大差约7.1×10⁻¹⁰。这个短例子说明构建、读取和核对方法，没有做MD、几何优化或材料参数拟合。

官方规则说明：Moltemplate手册。