📚 课程目录第3章轨迹预处理

3. 轨迹预处理：自定义topology文件的生成与检查（pdb2gmx）

好，咱们直接切入正题。上一章我们聊了怎么把原始结构文件里的残缺补全、加氢、理顺命名。但你有没有想过，为什么GROMACS非得让你跑一遍pdb2gmx，不能直接用你的结构？

说到底，GROMACS不认识什么PDB里的二面角、力常数。它只认自己的一套「暗号」——原子类型、电荷、成键参数。pdb2gmx这个工具，就是帮你的蛋白质翻译暗号的翻译官。我入行那年第一次跑流程，傻乎乎地跳过了这一步，结果mdrun直接报错「Fatal error: atomtype not found」。嗯，从那以后我就再也不敢省了。

一句话总结：pdb2gmx把一个静态结构坐标，翻译成了GROMACS能懂的拓扑力场文件。你后面所有的模拟、分析，根基就在这里。

3.1 核心命令：pdb2gmx参数扫盲

我建议你把这条命令存在笔记本里。每次新建项目，大概率都要跑一遍。

gmx pdb2gmx -f model.pdb -o complex.gro -p topol.top -ff amber99sb-ildn -water tip3p -ignh

参数逐个拆开说：

-f：输入结构文件，PDB或GRO格式都行。我个人习惯统一转成GRO后再喂进去，少踩格式兼容性的坑。
-o：输出的坐标文件，会补上氢原子和缺失的残基侧链坐标。
-p：输出的拓扑文件。
-ff：选力场。这个是关键中的关键，直接决定你后面模拟靠不靠谱。
-water：水模型。力场和模型要匹配，不能乱配。
-ignh：忽略PDB文件中原有的氢，让pdb2gmx自己重新加。这能避免力场和PDB中氢命名不一致的问题。

小技巧：如果你不确定PDB里有没有氢，直接加 -ignh 是安全牌。我曾经在一个结构里碰到PDB和力场对氢的命名规格差了一个字母，结果拓扑文件里一堆未定义的原子类型。

3.2 选择力场：用对的人，做对的事

你想想看，力场就是分子模拟里的物理规则。选错了，模拟跑得再长再稳也没用。我刚开始带团队时有个同事，拿CHARMM力场跑了GB1蛋白的折叠模拟，结果能量死活对不上——后来发现应该用CHARMM36m。

力场	适用场景	典型水模型
AMBER99SB-ILDN	球蛋白、酶，侧链参数优化得很好	TIP3P, TIP4P-Ew
CHARMM36m	膜蛋白、无序蛋白、核酸	CHARMM TIP3P/M
OPLS-AA/L	有机小分子、蛋白配体复合物	TIP3P, SPC
GROMOS 54A7	粗粒度和全原子互作，去折叠研究	SPC

警告：力场不要随便混搭。比如你用AMBER力场配了CHARMM的拓扑文件，后面的一个键角势项就能让你的模拟崩掉。遇到配体、非标准残基，老老实实用CGenFF或者ACPYPE单独生成。

3.3 自动化流程：你不需要每次都手打参数

实际项目里，你不会只跑一次pdb2gmx。跑突变扫描、做不同力场对比，一遍遍敲命令太傻了。写个脚本自动判断力场和水模型，才是工程师的做派。

#!/bin/bash
# auto_pdb2gmx.sh
# 自动判断力场和输出文件名

pdb_file=$1
system_name=$(basename "$pdb_file" .pdb)

# 根据残基数选择水模型和力场
residue_count=$(grep "ATOM" "$pdb_file" | awk '{print $5}' | sort -u | wc -l)

if [ "$residue_count" -gt 300 ]; then
  ff="amber99sb-ildn"
  water="tip3p"
else
  ff="charmm36m"
  water="charmm-tip3p"
fi

# 检查力场是否存在
if ! gmx pdb2gmx -h 2>&1 | grep -q "$ff"; then
  echo "[-] 力场 $ff 在系统上不可用"
  exit 1
fi

gmx pdb2gmx -f "$pdb_file" -o "${system_name}.gro" -p topol.top \
    -ff "$ff" -water "$water" -ignh

echo "[+] 完成: 输出 ${system_name}.gro, topol.top"
echo "[+] 使用力场: $ff, 水模型: $water"

生产建议：我习惯每次跑之前先加个 -apnd 参数，把pdb2gmx输出的内部诊断信息打印出来，快速扫一遍有没有警告。特别是「Atom X in residue Y not found in force field」这种，必须排查清楚再往下走。

3.4 常见陷阱：我踩过的那些坑

这么多年下来，pdb2gmx报的错我基本都见过。列几个高频的，你们提前避一避：

残基命名不对口：PDB里写HIE、HID、HIP，力场里只认HSD、HSE、HSP。跑之前先查查力场文档。
非标准残基：磷酸化、乙酰化、糖基化。pdb2gmx不认，需要先用CGenFF生成小分子参数，再用gmx的配体处理工具导入。
二硫键：PDB里两个CYS离得近不代表形成了二硫键。必须手动指定：在 -ss 后面写 "1 25"、"32 67" 这样成对的残基号。不指定的话，拓扑文件里没有桥键，模拟时两个巯基强行靠近，能量会炸。
末端电荷：默认配二末端NTER/CTER正常，但如果你做的是短肽片段，可能需要改电荷态——比如手动加上ACE/NME封端。

我曾经犯过一个低级错误：用pdb2gmx处理一个含有三个二硫键的小蛋白，忘了加 -ss 参数。模拟跑了20 ns，距离RMSD始终降不下来。最后查了轨迹才发现四个CYS的磺原子一直在来回晃动打转，压根没配对。加回二硫键约束再跑，10 ns就收敛了。

3.5 生成的文件，到底有啥用？

跑完之后你得到三个核心文件：

conf.gro：重新加氢、补残基后的坐标。直接用作后续模拟的初始结构。
topol.top：主拓扑文件。包含原子类型、电荷、质量、成键参数、二面角、非键参数、位置约束。
posre.itp：位置约束参数文件。对蛋白骨架（Cα、C、N）施加2000 kJ/(mol·nm²)的约束力，防止在能量最小化时蛋白飞到盒子外。你也可以在topol.top里include这个文件。

; topol.top 核心部分节选
[ atomtypes ]
; 原子类型定义
OPLS_AA_N  14.0067   -0.5000   A   3.50e-01  3.11e-01

[ moleculetype ]
; 分子类型名  接序数
Protein_chain_A        3

[ atoms ]
;   nr    type  resnr residue  atom   cgnr     charge    mass
     1  OPLS_AA_N       1    LYS      N      1    -0.5000  14.0067
     2  OPLS_AA_H       1    LYS     H1      2     0.3000   1.0080
; ... 后面还有成键、角度、二面角、位置约束定义

[ system ]
; 系统名
My simulation system

[ molecules ]
; 分子类型及数量
Protein_chain_A      1
SOL              10236
NA+                  18
CL-                  14

3.6 质量检查：跑完pdb2gmx后该看什么

别急着跳到下一步。花几分钟做三层检查：

检查拓扑文件的原子总数：和PDB实际的残基数核对。比如说你100个氨基酸残基，水模型不算，原子数大致在100*15到100*20之间。偏差太大说明有问题。
检查位置约束：看看posre.itp里约束了哪些原子。我遇到过一次因为残基顺序变了，位置约束约束错了残基，结果模拟中端到端距离完全下不来。
检查日志末尾：选-force_field 参数后，pdb2gmx会输出一段总结日志。里面会明确告诉你一共处理了多少残基、多少原子、用了什么力场和水模型。

# 快速检查示例
# 从topol.top里提取总原子数
grep -c "\[ atoms \]" topol.top   # 有几个分子类型
grep -c "^\s*[0-9]" topol.top    # 找原子定义行

# 查看结构是否合理
gmx check -f conf.gro -s topol.top | grep "Number of atoms"

最后提醒一句：别把pdb2gmx的输出当成真理。它只是根据力场把你给的坐标做了最合理的假设。真实的蛋白可能有修饰、有配体、有金属离子——这些都得你手动补。pdb2gmx不是魔法，它是你的得力帮手。

好，到这一步，你的拓扑文件已经生成、检查完毕，结构也从PDB的残缺状态变成了完整的、加了氢的、带全参数的系统。接下来就可以往溶液或者膜里扔进去了。