MDTraj入门：Trajectory对象创建，帧序号与时间信息，原子坐标访问

做分子动力学分析，第一步就是把轨迹文件读进内存。MDTraj 提供了非常简洁的接口，几行代码就能搞定。我个人习惯先加载拓扑再读轨迹，这样能避免很多莫名其妙的问题。

1. 创建 Trajectory 对象

MDTraj 支持两种方式创建 Trajectory 对象：

  • 从文件直接加载——最常用的姿势
  • 从内存数据构造——自己做预处理或拼接轨迹时用

先看第一种。你只要把轨迹文件路径和拓扑文件路径扔进去就行：

import mdtraj as md

# 从 xtc 加载，指定拓扑
traj = md.load('trajectory.xtc', top='topology.pdb')

print(traj)
# 输出示例：
# <mdtraj.Trajectory with 500 frames, 15246 atoms>

这里 top 参数可以是 PDB、GRO、HDF5 等格式。我遇到过一些初学者只传了 xtc 忘了给拓扑，结果 MDTraj 抛出异常说找不到原子信息。记住，拓扑信息必须提供，除非你的 xtc 本身自带拓扑（很少见）。

另一种情况是你要从内存创建对象。比如你从其他软件读到了坐标数据，想用 MDTraj 的函数做分析：

import numpy as np
import mdtraj as md

# 假设你有 100 帧，3 个原子
n_frames, n_atoms = 100, 3
xyz = np.random.randn(n_frames, n_atoms, 3) * 0.1   # 单位：纳米

# 构建一个简单的拓扑（3个碳原子）
topology = md.Topology()
chain = topology.add_chain()
residue = topology.add_residue('C1', chain)
for i in range(3):
    topology.add_atom('CA', element=md.element.carbon, residue=residue)

traj = md.Trajectory(xyz=xyz, topology=topology, time=np.arange(n_frames))
print(traj)

说实话，后一种场景我用得不多，但当你需要自定义测试数据、或者把多个轨迹片段合并时，这招特别管用。

2. 帧序号与时间信息

轨迹是一系列快照的集合。每个快照被称为一帧（frame），对应一个时间点。MDTraj 把这些信息都挂在了 Trajectory 对象上。

属性 类型 说明
n_frames int 帧总数
n_atoms int 原子总数
time ndarray, shape (n_frames,) 每帧对应的时间（皮秒）
timestep float 两帧之间的时间间隔（皮秒）

来看具体怎么用：

print('总帧数:', traj.n_frames)
print('时间点（前5个）:', traj.time[:5])
print('时间步长:', traj.timestep)

# 你也可以手动计算时间步长
dt = traj.time[1] - traj.time[0]
print('手动计算 dt:', dt)

我曾经帮一个同事排查过轨迹时间错位的问题。他加载的 xtc 里没有写入时间信息，MDTraj 默认 time 是从 0 开始的整数序列，但他以为时间是从 1000 ps 开始的。结果 RMSD 曲线怎么画都不对。后来我建议他先检查 traj.time[0]，才发现问题所在。所以你看，第一时间确认时间信息是否和你预期的匹配，这是血的教训。

我的小习惯： 加载后马上打印 traj 对象本身，看帧数、原子数、时间范围是否正确。这一步能帮你把 90% 的粗心错误扼杀在摇篮里。

3. 原子坐标访问

坐标的核心属性是 traj.xyz，一个三维数组。它的形状是 (n_frames, n_atoms, 3)。单位固定是纳米，这一点你一定要记住，不然和 VMD 里的单位混淆就麻烦了。

# 获取所有帧的坐标
all_coords = traj.xyz
print('坐标形状:', all_coords.shape)   # (500, 15246, 3)

# 访问第 10 帧的所有原子
frame_10 = traj.xyz[10]
print('第10帧形状:', frame_10.shape)   # (15246, 3)

# 访问第 10 帧的第 0 号原子坐标
atom0_frame10 = traj.xyz[10, 0]
print('原子0坐标:', atom0_frame10)     # [x, y, z] 纳米

# 更 Pythonic 的方式：直接用切片
first_atom_z = traj.xyz[:, 0, 2]       # 所有帧中第一个原子的 z 坐标
print('第一个原子的Z随时间变化:', first_atom_z[:5])

另外，MDTraj 还提供了索引访问的快捷方式。你可以像访问列表一样用方括号取单帧：

# 等价写法
frame_10 = traj[10]          # 返回一个新 Trajectory 对象（只有1帧）
frame_10_xyz = traj[10].xyz

# 多帧切片
sub_traj = traj[100:200:2]   # 从第100到199帧，步长2
print('子轨迹帧数:', sub_traj.n_frames)

注意 traj[10] 返回的是 Trajectory 对象，不是裸数组。如果你想直接拿 numpy 数组，还是用 traj.xyz[10]。个人建议在分析循环里一律用 traj.xyz 避免不必要的对象构造开销。

关键知识点：
  • xyz 形状 = (帧数, 原子数, 3)
  • 单位 = 纳米（不是 Å）
  • 索引帧时返回 Trajectory，索引 xyz 时返回 ndarray

4. 整个流程的关系图

下面这张图概括了 Trajectory 对象的内部结构与数据流向：

Trajectory 对象结构 .xtc / .dcd 拓扑文件 (.pdb/.gro) Trajectory 对象 n_frames, n_atoms, topology time (n_frames,) xyz (n_frames, n_atoms, 3) topology 核心操作 traj[frame] → Trajectory | traj.xyz[frame, atom, dim] → float

从图上可以看得清楚：文件加载 + 拓扑 → 构造 Trajectory → 通过 time 和 xyz 访问帧信息和坐标。如果你想用 traj[100] 取单帧，本质上还是从 xyz 里切片出来的。

5. 避坑指南

我曾经踩过的坑：
  • 单位混淆：MDTraj 内部坐标单位是纳米，但很多力场参数用的是 Å（0.1 nm）。如果你直接用 traj.xyz 计算距离然后跟文献里的 Å 数比，结果差 10 倍。建议统一一个单位，我习惯在程序头部就做转换。
  • 拓扑必须对应：如果拓扑文件中原子顺序和轨迹文件中的原子顺序不一致，加载不会报错，但后续分析全是错的。我会先跑一个 traj.topology 看看前几个原子的名字，跟预期对比一下。
  • 内存爆炸：一个 1000 帧、10 万原子的轨迹，traj.xyz 占内存约 1000×100000×3×8 ≈ 2.4 GB。如果你要处理超级大的轨迹，考虑用 md.iterload 分块读，这个后面课程会讲。

6. 一个完整的快速验证脚本

我每次拿到新轨迹都会跑下面这段代码，确认基本参数：

import mdtraj as md

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

print(f'帧数: {traj.n_frames}')
print(f'原子数: {traj.n_atoms}')
print(f'时间范围: {traj.time[0]:.1f} ~ {traj.time[-1]:.1f} ps')
print(f'坐标形状: {traj.xyz.shape}')

# 检查第一帧第一个原子的坐标
first_xyz = traj.xyz[0, 0]
print(f'第一个原子坐标: {first_xyz} nm')

# 检查时间步长一致性
dts = np.diff(traj.time)
if np.allclose(dts, dts[0]):
    print(f'时间步长均匀: {dts[0]:.4f} ps')
else:
    print('警告：时间步长不均匀！')

这条脚本能快速暴露八九成的问题，比如帧数不对、拓扑错位、时间不均匀等。我建议你把它加到自己的分析流水线开头。记住，花 30 秒确认基础数据正确，比花 3 小时分析出错结果划算得多。

好了，关于 MDTraj 的入门操作就聊到这里。咱们现在手里已经握住了轨迹的时间线和坐标数组，接下来就可以做 RMSD、回转半径、氢键分析这些真正有意思的事了。

← 上一章 📖 返回目录 下一章 →