4. 拓扑与原子选择工程
在分子动力学轨迹分析中，轨迹文件（如 XTC、DCD、TRR）通常只保存坐标、速度或力，而不保存原子名称、残基编号、链标识和成键信息。因此，任何一次分析都必须先把拓扑文件与轨迹配对，才能把“第 1732 号原子”翻译成“A 链第 87 号残基的 CB 原子”。拓扑与原子选择工程的核心任务，就是建立并维护这套映射关系，使选择语句、索引和实际原子三者始终一致。
拓扑文件解析：从文件到内存对象
常见拓扑来源有三类：PDB 文件（含坐标与部分连接信息）、结构文件加力场拓扑（如 GROMACS 的 .tpr、AMBER 的 prmtop）、以及 MDAnalysis 等库直接读取的 PSF、GRO、PDB。以 MDAnalysis 为例，加载时执行：
import MDAnalysis as mda
u = mda.Universe("system.prmtop", "traj.dcd")
print(u.atoms.n_atoms, u.residues.n_residues, u.segments.n_segments)
解析后，拓扑被组织为层级对象：Universe → segments（链）→ residues（残基）→ atoms（原子）。每个原子携带 index（全局 0 基索引）、id、name、type、resname、resid、segid、mass、charge 等属性。需要特别注意的是，index 是解析顺序决定的连续整数，而 id 来自文件中的原子序号，两者在多链或配体插入后往往不一致。分析脚本中若混用二者，会出现“选错原子但程序不报错”的隐蔽错误。
原子组选择语言：把化学意图写成表达式
选择语言的作用是把“我想分析蛋白主链”“我想看配体周围 5 Å 内的水”这类意图，转成可复现的原子集合。MDAnalysis 的选择语法以 select_atoms() 为入口，常用关键字包括 name、resname、resid、segid、protein、backbone、water、around、byres、same residue as 等。例如：
bb = u.select_atoms("backbone")
lig = u.select_atoms("resname LIG")
pock = u.select_atoms("protein and around 5.0 resname LIG", lig=lig)
选择语言的关键工程习惯是：先验证数量，再使用结果。执行 print(bb.n_atoms) 应等于“残基数 × 4（N、CA、C、O）”，若明显偏大，说明 backbone 匹配到了非标准残基或配体。另一个高频陷阱是 around 默认按原子距离筛选，若希望保留整个残基，应加 byres，即 "byres (protein and around 5.0 resname LIG)"，否则会出现残基被“切一半”的原子组。
残基与原子索引映射：把选择结果落到数组
选择语句返回的是 AtomGroup，其底层是 indices 数组。真正做数值计算时，应显式取出索引并核对：
idx = bb.indices
print(idx[:5], idx.shape)
# 用索引直接取坐标，避免重复选择
coords = u.trajectory[0].positions[idx]
残基级映射则通过 residues 与 atoms 的对应关系完成。例如要统计每个残基的 Cα 坐标，可先选 name CA，再按 resids 分组。若拓扑来自 PDB，残基编号可能因缺失残基而跳号，此时 resid 不连续，但 resindex 始终连续。做逐残基循环时，用 resindex 更安全；做与实验编号对照时，用 resid 更直观。两者不可混用。
多链与配体处理
多链体系最常见的问题是残基编号在链间重复。假设 A 链和 B 链各有 resid 1–120，执行 u.select_atoms("resid 50") 会同时选中两条链的第 50 号残基，原子数翻倍。正确做法是加上链标识："segid A and resid 50"。若拓扑来自 PDB 且没有 segid，可用 chainID 或先按 segments 拆分再选择。
配体处理有三个要点。第一，配体通常作为独立 segment 或独立 residue 存在，选择时应同时限定 segid 与 resname，避免与标准残基重名。第二，配体的原子名可能包含数字或特殊字符，选择时建议用 resname LIG 整体选中，再按 name 细分。第三，若配体在轨迹中发生键断裂或质子化状态变化，拓扑中的成键信息可能与实际不符，此时应避免依赖 bonded to 选择，改用距离或残基整体选择。
selections.py 或配置字典中，分析脚本只引用名称。这样当拓扑更新、链编号变化时，只需改一处，避免散落在各脚本中的硬编码索引失效。一个可复用的检查清单
| 检查项 | 命令/方法 | 预期 |
|---|---|---|
| 原子总数 | u.atoms.n_atoms | 与拓扑文件一致 |
| 链数 | u.segments.n_segments | 蛋白链数 + 配体 + 水/离子 |
| 主链原子数 | u.select_atoms("backbone").n_atoms | ≈ 残基数 × 4 |
| 配体原子数 | u.select_atoms("resname LIG").n_atoms | 与配体化学式一致 |
| 残基编号连续性 | u.residues.resids | 允许跳号，但需记录 |
| index 与 id 差异 | (u.atoms.indices != u.atoms.ids).sum() | 多链/配体后通常 > 0 |
完成上述检查后，再进入轨迹读取与坐标计算，才能保证后续每一帧取到的原子与拓扑定义严格对应。拓扑与选择工程做扎实，后面的 RMSD、氢键、接触分析才有可靠基础。