Moltemplate的MOL2导入与PDB整理

Moltemplate可以把MOL2中的坐标、电荷和连接读成LT，也提供PDB排序与坐标提取工具。本篇使用自建四珠子MOL2链和五条PDB格式记录，演示具体命令与需要检查的字段。它们是格式和建模示例，不是实验蛋白或已经参数化的化学分子。

1. 下载与运行#

下载控制点、LT模板、MOL2/PDB、真实结果及绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd moltemplate-curve-import-tools

先查看已有结果：

bash show-reference.sh
cat curve-import-checks.csv
head -n 5 curve-geometry.csv

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

sbatch job.slurm

CPU模板申请56核整节点，实际短示例单进程执行；分区维护时等待可用资源。在已分配资源中可用 bash run.sh。它依次生成曲线和LT、加入指定交联键、导入两种MOL2电荷、用LAMMPS读回四个模型，再提取PDB坐标并画图。重复整套流程请换新目录，避免覆盖已有 work/。

本例使用Moltemplate2.22.5/Python3.11，以及Intel oneAPI2024.2、Intel MPI/MKL构建的LAMMPS30Sep2026 CPU版；分析绘图通过 sisl/0.16.4 的Python完成。无需修改公共环境。

2. 先看MOL2里面有什么#

完整 inputs/toy.mol2：

@<TRIPOS>MOLECULE
OwnFourSiteTeachingChain
4 3 1 0 0
SMALL
USER_CHARGES

@<TRIPOS>ATOM
1 B1 0.000000 0 0 bead 1 TOY 0.100000
2 B2 1.400000 0 0 bead 1 TOY -0.100000
3 B3 2.800000 0 0 bead 1 TOY -0.100000
4 B4 4.200000 0 0 bead 1 TOY 0.100000
@<TRIPOS>BOND
1 1 2 1
2 2 3 1
3 3 4 1
@<TRIPOS>SUBSTRUCTURE
1 TOY 1

四个珠子名B1到B4，原子类型为本例力场定义的 bead。坐标沿x方向排列，电荷为+0.1、−0.1、−0.1、+0.1，三条键连接相邻珠子。实际化学体系应使用对应力场认可的原子类型和电荷，MOL2格式转换本身不会替你完成可靠参数化。

3. 转成LT，再生成LAMMPS输入#

source ./env.sh
mkdir customer-mol2
cp inputs/chainff.lt inputs/toy.mol2 inputs/override-charges.txt customer-mol2/
cp inputs/mol2-system.lt customer-mol2/system.lt
cd customer-mol2
mol22lt.py --in toy.mol2 --out molecule.lt --name ImportedChain --ff ChainFF --ff-file chainff.lt
moltemplate.sh system.lt
grep -E 'atoms$|bonds$' system.data

--name 指定LT对象名；--ff 与 --ff-file 指向本例实际存在的力场。MOL2读取后的连接进入 Data Bond List，这里通过 Data Bonds By Type 规则确定谐键类型。

安装版本的 mol22lt.py 明确属于实验性工具。它读取ATOM/BOND的相应字段，但不保留MOL2第四列所写的化学键级含义，也不自动证明多链、特殊子结构或所有第三方MOL2写法兼容。导入后要核对全部原子、拓扑、类型和参数；不要把“生成了LT”理解为力场已完整可用。

4. 需要替换电荷时#

本例外部电荷文件一行一个值：

0.2
-0.2
-0.2
0.2

在独立目录里带上 --charges：

source ./env.sh
mkdir customer-mol2-override
cp inputs/chainff.lt inputs/toy.mol2 inputs/override-charges.txt customer-mol2-override/
cp inputs/mol2-system.lt customer-mol2-override/system.lt
cd customer-mol2-override
mol22lt.py --in toy.mol2 --out molecule.lt --name ImportedChain --ff ChainFF --ff-file chainff.lt --charges override-charges.txt
moltemplate.sh system.lt
sed -n '/Atoms/,+6p' system.data

本例电荷变为+0.2、−0.2、−0.2、+0.2，总电荷仍为0。这里 pair_style zero 不计算库仑相互作用，因此两组电荷下的势能相同；这是电荷字段导入示例，不能当成静电能计算。真实带电体系还要选择并配置相应势函数与长程处理。

5. PDB先排序，再提取坐标#

PDB不是按空白任意分列的坐标表。链ID、残基号、插入码、原子序号和坐标都有固定列。本例故意将五条ATOM记录打乱，再按当前工具的规则整理：

source ./env.sh
pdbsort.py < inputs/unsorted.pdb > customer-sorted.pdb
pdb2crds.awk < customer-sorted.pdb > customer-sorted.raw
cat customer-sorted.raw

实际排序后的原始原子序号为4、5、3、2、1，对应先链ID、再残基号、插入码和原子序号；五行坐标均按PDB固定列提取。pdb2crds.awk 只输出坐标，不生成类型、键、电荷或力场。

当前 pdbsort.py 会把非ATOM/HETATM行输出在原子记录之前，例如本例END也会出现在前面。因此这里的输出仅用于后续坐标提取，不能直接当作保留完整记录顺序的PDB归档文件。多模型、替代构象等情形还需按自己的目的另行选择，原始PDB应始终保留。

6. 实际结果#

两组结构都有4珠子、3键，实际LAMMPS run 0 的总势能为3，与三条K=100、r0=1.3、实际键长1.4的谐键独立求和一致。全部坐标、质量、电荷和连接也逐项核对。这里只演示转换和读回，没有进行结构优化或MD。

需要把已有DATA再还原成LT，可继续看模板还原与类型清理；只处理轨迹快照，可看轨迹与DATA转换。官方说明：mol22lt。