2.2 pdb2gmx 拓扑生成实操

2.2 pdb2gmx 拓扑生成实操

本节摘要：pdb2gmx 是 GROMACS 的拓扑生成器：读入结构，补氢、指派原子类型与电荷、写出全部成键与非键定义。本节以一条多肽走通它的完整交互流程，逐项解读端基、质子化与二硫键问答，再用原子数对账法验证输出；随后转向 pdb2gmx 认不了的分子——聚合物与脂质——给出 CGenFF 参数化与手写 itp 两条替代路线，并演示拓扑文件的骨架结构，为 2.3 的装配备好图纸。

2.1 把参数验对了，现在要把参数落到具体分子上。拓扑文件的发明理由很朴素：模拟器每步都要问"这堆原子谁连着谁、连的是哪种弹簧、电荷多少"，与其每次重算，不如在开工前登记成册。pdb2gmx 就是给蛋白、核酸这类"标准件分子"自动办册子的窗口；而非标准件（聚合物、脂质、绝大多数药物分子）得走别的窗口。两条路本节都走一遍。

标准流程：一条多肽的完整会话

先交代 pdb2gmx 实际做了四件事：加氢（PDB 结构通常只有重原子）；识别残基与端基并指派质子化态；查表指派原子类型与电荷；写出成键表（键、角、二面角、排除表）。命令一行，问答两三处：

$ gmx pdb2gmx -f peptide.pdb -o processed.gro -p topol.top -i posre.itp \
              -ff amber14sb -water tip3p -ignh
                  :-) GROMACS - gmx pdb2gmx, 2024.2 (-:
...
Opening force field file ./amber14sb.ff/forcefield.itp
Identified residue LYS2 as a starting terminus.
Identified residue LEU9 as an ending terminus.
Start terminus: NH3+          # 问你 N 端选哪种：带电 NH3+ 还是中性 NH2
End terminus: COO-            # C 端：带电 COO- 还是中性 COOH
...
Total charge in system: 1.000e+00   # 拓扑记账：净电荷 +1
​

三个问答的决策依据要讲清楚。端基：除非你的化学明确是封端修饰（如乙酰化 N 端），实验条件下端基几乎总是带电的——默认选带电态，别为了"好看"选中性。质子化：组氨酸最常被问（HIS 的 HISA/HISB/HISH 变体对应氮上氢的位置或带电态），pKa 附近的残基要结合体系 pH 用经验判断；拿不准时算完可以两种态各跑短模拟对比稳定性。-ignh：丢弃 PDB 自带的氢重新生成——实验结构里的氢位置本就不可靠，删掉让力场自己的规则来加，几乎总是更稳的选择。

输出对账是这一步的关键动作。pdb2gmx 会打印每种残基的原子数变化，你要做的是手工复算总账：

$ grep -c "^ATOM" peptide.pdb        # 原始重原子数，例如 128
$ grep -v "^#" processed.gro | head -n 2 | tail -n 1   # gro 第二行 = 总原子数，例如 261
# 差值 133 = 新加的氢 + 端基处理。数量级 sanity check：
# 每个残基平均带约 15-16 个氢，9 残基加两个端基氢，100 上下合理
​

对不上的情况只有两种来源：PDB 里有 pdb2gmx 不认识的残基名（报 "Residue X not found"），或结构里有链断裂（它会把链拆成多个分子块分别记账）。遇到前者去改 PDB 命名规范（如把 HSD 改成 HIS），遇到后者用 -chainsep 或手动补缺口。

二硫键的问答值得单独一提：两个 CYS 硫原子距离低于阈值时，pdb2gmx 会问"这算二硫键吗"。确认前先目测距离——真正的 S–S 键长约 0.203 nm；把两个只是靠得近的自由半胱氨酸误连成二硫键，是能把折叠方向整个带偏的灾难级手滑。

替代路线：pdb2gmx 不认识 PEO，怎么办

聚合物、脂质、药物小分子不在 pdb2gmx 的残基字典里，标准做法是外部参数化：CGenFF（CHARMM 生态）、GAFF（AMBER 生态）或 OpenFF，输入分子结构，输出带类型与电荷的参数，再翻译成 GROMACS 的 itp。流程以 CGenFF 为例：

# 在 CGenFF 服务或本地程序上提交 PEO 单体-链的 mol2 结构，得到 polyethylene_glycol.str
# 翻译为 GROMACS 格式（社区工具 cgenff_charmm2gmx 一类脚本）
$ python cgenff_charmm2gmx.py PEO polyethylene_glycol.mol2 \
            polyethylene_glycol.str charmm36.ff
# 产物：peo.itp（分子定义）+ peo.pdb（带类型标注的结构）
​

拿到 itp 后要过一遍"半自动质检"：CGenFF 会给每个片段打罚分（penalty），高于阈值（通常取 10）的片段表示外推成分大，参数可信度存疑——对聚合物而言整条链由重复单体拼成，只需把单体参数化一次，再用拓扑预处理器的多重复用语法把单元串起来。手写 itp 的骨架长这样：

; peo.itp —— 20 单体 PEO 链的手写骨架（联合原子示意）
[ moleculetype ]
PEO20  3                    ; 分子名，nrexcl=3 表示隔三个键内不算非键

[ atoms ]
;  nr  type  resnr residue  atom  cgnr  charge  mass
    1   CH3     1    PEO    CH3     1      0.0   15.035
    2   CH2     1    PEO    CH2     1     0.10   14.027
    3    OE     1    PEO     O      1    -0.20   15.999   ; 醚氧带负电
    ...                         ; 每单体三个格点，20 单体共 60 格点

[ bonds ]    ; 邻接格点两两成键，参数引用 2.1 核对过的表
[ pairs ]    ; 1-4 对：按 nrexcl 自动生成，或显式列出
[ angles ]
[ dihedrals ]

#ifdef POSRES                ; 平衡阶段用的位置限制开关，3.6 会用到
#include "posre_peo.itp"
#endif
​

手写路线的纪律只有一条：单体重复单元的参数绝不允许有两套。把 60 个格点逐个手写难免抄错一个电荷，而电荷总和错 0.05 e，在 2.3 的中和步骤就会现形。稳妥做法是写成"单体块 + 链接键"的生成脚本，链长变成一个参数。

拓扑的骨架与对账终点

无论哪条路线，最终都要汇进 topol.top。它的固定格式是"包含区 + 体系区"：

#include "charmm36.ff/forcefield.itp"     ; 力场总纲：函数形式与默认规则
#include "charmm36.ff/tip3p.itp"          ; 水模型
#include "peo.itp"                        ; 你的分子（2.1 验过参数，2.2 生成的定义）

#ifdef POSRES
#include "posre.itp"
#endif

[ system ]
PEO20 in water, 20 chains, 5 wt percent

[ molecules ]
; 分子名        数目
PEO20           20
SOL             115        ; 2.3 装水后会改写
​

最后的对账公式贯穿后续所有环节：总原子数 = Σ(分子原子数 × 分子数)，净电荷 = Σ(分子电荷 × 分子数)。2.3 加离子时改写的就是 [ molecules ] 区，改完各算一遍，账平了才许进第三章。

速记：pdb2gmx 只服务标准件分子，问答三处（端基、质子化、二硫键）每个都有物理含义；非标准件走外部参数化加 itp 引用；单体参数化一次、脚本复用，拒绝手抄；拓扑对账两本账——原子数与净电荷。图纸齐了，下一节搭台：盒子、水、离子。