📚 课程目录第2章轨迹文件格式与I/O工程

2. 轨迹文件格式与I/O工程

分子动力学模拟产生的轨迹文件，是后续所有分析工作的数据基础。一个典型的全原子体系在 1 ns 采样、每 1 ps 保存一帧时，就会产生上千帧坐标；若体系含数十万原子，单文件体积可达数十 GB。因此，理解轨迹格式的存储结构、掌握高效的读写与流式处理方式，是轨迹分析工程的第一道门槛。本章围绕 DCD、XTC、TRR、NetCDF、HDF5 五类主流格式，以及 MDAnalysis、MDTraj、GROMACS 三类工具的 I/O 能力展开。

一、格式的本质差异：存储模型决定访问方式

轨迹格式的差异，根源在于“如何组织坐标、时间、盒子与拓扑”这一存储模型。DCD 是 CHARMM 系格式，采用固定步长、定长记录的二进制布局，只存坐标（可选存盒子），不存拓扑、时间与速度，因此必须配合 PDB/PSF 等拓扑文件才能还原体系。XTC 是 GROMACS 的压缩轨迹格式，使用 3D 无损压缩，坐标精度约 0.001 nm，体积通常只有 DCD 的 1/3～1/5，但压缩/解压带来 CPU 开销。TRR 是 GROMACS 的全精度格式，保留坐标、速度、力与盒子，适合需要速度做动力学分析或续算的场景。NetCDF（AMBER 系）与 HDF5（如 NAMD、OpenMM 常用）都是自描述容器格式，可在同一文件内组织多组数据集，HDF5 还支持分块与并行 I/O。

这些差异直接决定了访问模式：定长记录格式（DCD）可以按帧号直接 seek 到任意帧，随机访问成本低；压缩格式（XTC）必须顺序解压，随机访问代价高；容器格式（HDF5）则可通过数据集路径与切片灵活读取。选择格式时，应在“体积、精度、随机访问、并行能力”之间权衡。

二、用 MDAnalysis 做统一读写

MDAnalysis 的核心抽象是“Universe = 拓扑 + 轨迹”。拓扑文件（PDB、PSF、GRO、PRMTOP）提供原子名、残基、质量、键连等信息，轨迹文件提供逐帧坐标。加载时把两者绑定，即可用统一接口访问任意格式。

import MDAnalysis as mda

# 拓扑 + 轨迹绑定，自动识别格式
u = mda.Universe("system.pdb", "traj.xtc")

print(u.trajectory.n_frames)      # 总帧数
print(u.trajectory.dt)            # 帧间隔（ps）

# 顺序遍历：每帧计算质心
for ts in u.trajectory:
    com = u.atoms.center_of_mass()
    # ts.frame, ts.time 可直接读取

# 随机访问：直接跳到第 500 帧
u.trajectory[500]
print(u.trajectory.time)

# 切片读取：每 10 帧取一帧
for ts in u.trajectory[::10]:
    pass

这里有一个工程要点：u.trajectory[i] 对 DCD 是真正的 seek，对 XTC 则可能触发从最近关键帧开始的顺序解压。因此对 XTC 做大量随机访问时，应尽量按帧号升序批量读取，或先用 gmx trjconv 转成 DCD/HDF5 再随机访问。

三、MDTraj 的轻量读取与流式处理

MDTraj 以“加载即得坐标数组”为设计取向，适合快速统计与机器学习特征提取。它支持 load 一次性载入，也支持 iterload 分块流式读取，后者是处理超大轨迹的关键。

import mdtraj as md

# 一次性载入（小轨迹）
traj = md.load("traj.dcd", top="system.pdb")
print(traj.n_frames, traj.n_atoms)

# 流式分块：每次 100 帧，内存占用可控
for chunk in md.iterload("traj.xtc", top="system.pdb", chunk=100):
    ca = chunk.atom_slice(chunk.top.select("name CA"))
    # 对 chunk 做逐块分析，如计算回旋半径
    rg = md.compute_rg(ca)

iterload 的价值在于：无论轨迹多大，内存中始终只保留一个 chunk。对于需要跨帧累积的统计量（如平均结构、直方图），应在循环外维护累加器，而不是把每帧结果都存进列表。

四、GROMACS 命令行 I/O 与格式转换

GROMACS 的 gmx trjconv 是格式转换与轨迹裁剪的主力工具。常见操作包括：把 XTC 转成 DCD 以便随机访问、抽取指定帧区间、去除周期性边界条件、把轨迹对齐到参考结构。

# XTC 转 DCD，便于随机访问
gmx trjconv -f traj.xtc -s topol.tpr -o traj.dcd

# 只取 1000–2000 ps，每 10 ps 一帧
gmx trjconv -f traj.xtc -s topol.tpr \
  -b 1000 -e 2000 -dt 10 -o slice.xtc

# 去除 PBC 并居中，输出到新轨迹
gmx trjconv -f traj.xtc -s topol.tpr \
  -pbc mol -center -o traj_nopbc.xtc

执行 trjconv 时，程序会交互式询问选择组（如 System、Protein、Backbone）。在脚本中可用 echo 管道传入选择，例如 echo "Protein System" | gmx trjconv ...，避免交互阻塞批处理。

五、流式处理的工程原则

面对超出内存的轨迹，核心原则是“按帧或按块处理，只保留必要状态”。具体做法包括：用 iterload 或 u.trajectory[::stride] 控制读取粒度；对每帧只提取分析所需的原子子集（如 CA 原子），而不是全原子；把逐帧结果写入磁盘（CSV、HDF5），而不是累积在内存列表；对可并行的逐帧计算，用多进程按帧区间分块，各进程独立打开轨迹文件。

格式选择建议：长期归档用 XTC 省空间；需要速度/力或续算用 TRR；需要随机访问与并行 I/O 用 HDF5；与 AMBER 生态交互用 NetCDF；与 CHARMM/NAMD 交互用 DCD。转换前先确认拓扑与轨迹的原子顺序一致，否则坐标会错位。

常见陷阱：DCD 不含时间信息，MDAnalysis 只能按帧号推断时间，若原始采样间隔不均匀会失真；XTC 压缩会引入约 0.001 nm 量化误差，对键长/键角级别的分析应改用 TRR 或 NetCDF；HDF5 文件若被多个进程同时写入，必须使用并行 HDF5 或分文件写入，否则会损坏。

掌握上述格式特性与工具接口后，就能根据分析目标选择合适格式，并用流式策略稳定处理任意规模的轨迹，为后续的 RMSD、RMSF、回旋半径、氢键、聚类等分析打下可靠的数据基础。