MDAnalysis读取轨迹、周期边界与RMSD

下载本篇输入、完整脚本、数据与高清图。解压后在新目录运行，reference保存鸢算科技实际运行的参考结果。

1 准备拓扑和轨迹#

本例使用MDAnalysis 2.9.0及NumPy，在CPU上分析已有GROMACS TIP3P水轨迹：884个水分子、2652个原子，0–100 ps共1001帧，存帧间隔0.1 ps。下载包提供真实GRO/XTC文件、脚本与参考输出。

解压到新的工作目录后可一次运行全部分析：

bash run.sh

只运行MDAnalysis时：

mkdir mda
cp -a inputs/mda mda/input
cp mdanalysis_examples.py mda/
cd mda
/public/toolkit/miniconda3/envs/qvasp-2.26/bin/python mdanalysis_examples.py

这里直接调用已有Python环境，无需安装新依赖。脚本使用4线程以内资源；更大轨迹请在分配到的计算资源中处理。

2 读取文件与检查单位#

import MDAnalysis as mda
u = mda.Universe('input/water.gro', 'input/water.xtc')
print(len(u.atoms), len(u.residues), len(u.trajectory))
print(u.trajectory.dt, u.dimensions)

输出原子、残基、帧数应分别为2652、884、1001。MDAnalysis内部长度用Å、时间用ps；GROMACS原文件常用nm，读取时会转换。不能将读取后的30.13 Å盒长再当作30.13 nm。

本次GROMACS 2026的TPR版本138不能由MDAnalysis 2.9.0直接读取，因此示例使用GRO提供原子名称和残基映射。GRO没有完整力场拓扑、键和电荷；需要这些属性的分析必须补充对应拓扑或明确提供参数，不能把名称映射当成完整拓扑。

3 选择原子与周期边界距离#

oxygen = u.select_atoms('resname SOL and name OW')
hydrogen = u.select_atoms('resname SOL and name HW*')
print(len(oxygen), len(hydrogen))
from MDAnalysis.lib.distances import distance_array
d = distance_array(oxygen.positions[:10],
                   oxygen.positions[:10], box=u.dimensions)

得到884个O和1768个H。显式传入当前帧的box，距离才会使用周期最小镜像。脚本将一个氧整体移动一个周期盒后，距离矩阵最大差2.59×10⁻⁶ Å；这是坐标精度范围内的周期等价检查。

按残基对应OW与HW1，本次平均O–H距离为0.957247 Å。不要把盒两侧的一对成键原子直接做普通欧氏距离，否则可能得到接近盒长的“键”。

4 LAMMPS格式与GRO/DCD往返#

脚本将真实第一帧写为LAMMPS dump，明确设类型1为O、2为H，再由MDAnalysis读取，2652个坐标无差异。这个文件只是坐标数据，没有自动生成LAMMPS分子拓扑或势函数。

u.trajectory[0]
u.atoms.write('first.gro')
with mda.Writer('sample.dcd', n_atoms=len(u.atoms), dt=1.0) as w:
    for ts in u.trajectory[:101:10]:
        w.write(u.atoms)

DCD选取0–10 ps每1 ps一帧，共11帧，回读坐标最大差0 Å。GRO回读最大坐标差3.81×10⁻⁶ Å；不同输入精度时允许误差需按格式调整。DCD时间间隔要与采样间隔一致，格式转换不会自动保存全部力场信息。

5 单分子对齐、RMSD和RMSF#

先用周期最小镜像重建完整的第一个水分子，再移除整体平移、转动，对每1 ps的101帧计算RMSD与RMSF。完整实现见mdanalysis_examples.py的rms()：

from MDAnalysis.analysis.rms import RMSD, RMSF
from MDAnalysis.analysis.align import AlignTraj
r = RMSD(v, select='all').run()
AlignTraj(v, v, select='all', in_memory=True).run()
f = RMSF(v.atoms).run()
print(r.results.rmsd, f.results.rmsf)

其中v是脚本准备好的完整水分子轨迹。最大拟合RMSD为 0.007504 Å，O、H、H的RMSF约 0.003071、0.002646、0.002508 Å。

TIP3P本例是刚性水，小幅残差主要体现XTC坐标压缩精度，不代表真实键振动。分析蛋白质时应改用合适的主链选择，先处理分子跨边界，再对齐，最后解释柔性。

参考：MDAnalysis分析模块。