4、读取轨迹文件：mdtraj.load() 的十八般武艺

读取轨迹，是分子模拟分析的第一步。说白了，这一步没走对，后面全是白费功夫。我刚开始用 MDTraj 时，就被各种文件格式和 topology 参数折磨过。

今天我们来聊聊 mdtraj.load() 这个核心函数。它到底能读什么？怎么读？topology 参数又是干嘛的？

4.1 支持的格式：比你想象的要多

MDTraj 支持的轨迹格式相当丰富。我粗略数过，大概有20多种。这里列出最常用的：

格式	扩展名	典型来源	备注
PDB	.pdb	Amber, GROMACS, CHARMM	可含多帧，但我不建议用它存大量轨迹
XTc	.xtc	GROMACS	压缩率高，我最常用的格式之一
TRR	.trr	GROMACS	全精度，包含速度/力
DCD	.dcd	CHARMM, NAMD	老牌格式，生物领域常见
NetCDF	.nc	Amber	Amber 用户的标配
HDF5	.h5	MDTraj 原生格式	我个人很推荐，读写都很快
LAMMPS dump	.dump	LAMMPS	需要指定原子类型映射

你知道吗？mdtraj.load() 其实不需要你手动指定格式。它会根据文件扩展名自动判断。但有些场景下，扩展名不标准或文件没有后缀名，这时你就得用 format= 参数了。

4.2 topology 参数：为什么必须要有它？

轨迹文件存的是原子坐标，但它不存原子名称、残基信息、成键关系。这些信息存放在拓扑文件里。

拿 XTC 举例。一个标准的 XTC 文件只存坐标和时间。如果你直接用 load('traj.xtc')，MDTraj 会抛出异常：

ValueError: No topology provided and trajectory format does not contain one.
Please provide the 'top' argument.

嗯，这里的教训很直接。你得给它一个拓扑文件：

import mdtraj as md

# 正确方式：提供拓扑文件
traj = md.load('traj.xtc', top='system.pdb')

# 或者用GroTop文件
traj = md.load('traj.xtc', top='system.gro')

这里有个细节：拓扑文件可以是 PDB、GRO、PRMTOP、PSF 等格式。MDTraj 会从这些文件中提取出完整的拓扑信息。

我曾经踩过一个坑：用 Amber 的 prmtop 作为拓扑，读 XTC 轨迹。结果 MDTraj 直接挂掉了。为什么？因为 prmtop 需要配合 NetCDF 轨迹的单元盒信息，而 XTC 的盒子表示方式不同。后来我改用 PDB 作为拓扑文件，才解决问题。

4.3 原子选择：有了拓扑，你就是上帝

一旦有了拓扑，你就获得了原子选择的能力。这是 MDTraj 最让我喜欢的地方。

看下面这段代码：

import mdtraj as md

traj = md.load('traj.xtc', top='system.pdb')

# 选择所有蛋白质原子
protein = traj.topology.select('protein')

# 选择所有C-alpha原子
ca_indices = traj.topology.select('name CA')

# 选择残基1到10的骨架
backbone = traj.topology.select('resid 1 to 10 and backbone')

# 选择水分子中的氧原子
water_oxygen = traj.topology.select('resname HOH and name O')

# 获取子轨迹
protein_traj = traj.atom_slice(protein)

选择语法和 VMD 类似，但更简洁。常用的关键词有：

选择器	作用	示例
protein	选择所有蛋白质原子	top.select('protein')
name CA	按原子名称选择	top.select('name CA or name CB')
resid 1 to 20	按残基序号范围选择	top.select('resid 1 to 20')
resname ALA	按残基类型选择	top.select('resname ALA')
backbone	骨架原子（N, CA, C）	top.select('backbone')
sidechain	侧链原子	top.select('sidechain')
water	水分子	top.select('water')
not water	除水外的所有原子	top.select('not water')
element C	按元素选择	top.select('element C')

选择器支持逻辑运算：and、or、not。也支持括号分组。

举个例子，我想选择亮氨酸和异亮氨酸的侧链原子：

target = traj.topology.select('(resname LEU or resname ILE) and sidechain')

4.4 实战：从 PDB 开始，一步步来

我建议新手先从 PDB 文件开始。因为 PDB 自带拓扑，不需要额外传 top 参数。

import mdtraj as md

# 单帧PDB
traj = md.load('protein.pdb')
print(traj.n_frames)   # 大概率是1
print(traj.n_atoms)    # 原子数
print(traj.n_residues) # 残基数

# 多帧PDB（有些MD程序会输出多帧PDB）
traj_multi = md.load('trajectory.pdb')
print(traj_multi.n_frames)  # 帧数

# 混合：拓扑来自PDB，轨迹来自XTC
traj = md.load('traj.xtc', top='ref.pdb')

我习惯把参考结构存成 PDB，然后用它作为拓扑文件。这样即使原始轨迹文件被删了，我还能靠 PDB 重建拓扑关系。

4.5 避坑指南

这里列几个我亲眼见过的坑：

• 拓扑不匹配：拓扑文件的原子数必须和轨迹文件一致。不一致的话，MDTraj 会直接报错，不会含糊。
• 单元盒信息丢失：有些轨迹格式不存盒信息。比如不指定盒子的 PDB。这时需要手动设置，用 traj.unitcell_vectors = ... 赋值。
• 原子顺序混乱：不同程序输出轨迹的原子顺序可能不同。我遇到过从 GROMACS 转成 PDB 后，水分子顺序变了的案例。一定要用 topology.select() 确认一下。
• 大文件读取太慢：几百 MB 的 XTC 还好。但如果是几千帧的 PDB，读取速度会让人崩溃。我会先用 mdtraj.load_frame() 读单帧测试，确认没问题再批量读取。

4.6 一手经验：如何判断读取成功了？

我每次加载轨迹后，都会快速检查几个关键属性：

traj = md.load('production.xtc', top='start.pdb')

# 快速检查
print('帧数:', traj.n_frames)
print('原子数:', traj.n_atoms)
print('残基数:', traj.n_residues)
print('盒子向量:', traj.unitcell_vectors.shape)

# 查看第一帧的坐标范围
import numpy as np
print('坐标范围:', np.min(traj.xyz[0]), '~', np.max(traj.xyz[0]))

如果坐标范围异常大（比如超过 10 nm），通常说明拓扑文件和轨迹文件不匹配，或者单位不统一。MDTraj 默认用纳米为单位。如果你的轨迹是埃的，需要做转换。

说到单位，我想到一个案例。有一次我读了一个 LAMMPS 的 dump 文件，坐标直接飞到了几千。后来发现 LAMMPS 输出的单位是埃，而 MDTraj 默认是纳米。乘以 0.1 就恢复正常了。

# 如果坐标单位是埃，手动转成纳米
traj.xyz = traj.xyz * 0.1

4.7 本章要点总结

读取轨迹这件事，核心就是三点：

1. 选对方法：mdtraj.load() 是万能入口。遇到特殊格式，可以用 format= 提示。
2. 配好拓扑：没有拓扑，MDTraj 就无法知道原子的归属。PDB 是最保险的拓扑文件来源。
3. 用好选择器：topology.select() 的语法很直观。花几分钟记住常用关键词，后续分析效率翻倍。

好了，这章就到这里。下次我们聊轨迹的帧操作——怎么切帧、怎么对齐、怎么把轨迹拆成片段。这些都是实打实的实用技能。