5.1 轨迹预处理与周期性边界

5.1 轨迹预处理与周期性边界

本节摘要：预处理是分析前的修图工序：周期性边界会把跨盒分子"切"成两半、让质心在盒子间跳动，直接拿原始轨迹统计必然得到伪影。本节用 trjconv 的 -pbc mol/whole/nojump 与 -center 组合，按分析目标给出三种标准配方，并用一帧示意图解释每种模式修的是什么病。记住本章的开工口令：先修几何，再算物理。

4.2 把轨迹入了仓，本节开工处理。先讲病根。2.3 说过，模拟里的分子可以在右壁出去、左壁进来，坐标本身永远是"盒子内的那一份"。于是轨迹文件里常看到两类怪象：分子断裂——一条 PEO 链的端基出现在盒子对侧，可视化时链被切成两截；整体跳动——分子质心恰好走到边界附近时，整个分子从盒子一头"瞬移"到另一头（坐标被折回 0 到 L 区间）。这两类怪象在模拟物理里完全无害（最小镜像保证相互作用正确），在统计分析里却是毒药——任何基于"相邻帧位移"或"分子完整性"的计算都会被污染。

三种 -pbc 模式各修什么

三种模式处理的是三种不同的病，常要组合使用。标准配方一（可视化与通用分析）：先 whole 拼整分子，再 nojump 消跳动，再 center 居中：

$ printf "PEO\n" | gmx trjconv -s md.tpr -f md.xtc -o fix1.xtc -pbc whole
$ printf "PEO\n" | gmx trjconv -s md.tpr -f fix1.xtc -o fix2.xtc -pbc nojump
$ printf "PEO\nSystem\n" | gmx trjconv -s md.tpr -f fix2.xtc -o fix3.xtc \
        -pbc mol -center
​

（交互式选组时第一个提示选 PEO 作为居中参考组、第二个选 System 作为输出组；脚本里用 printf 管道喂选项，保持可复现。）配方二（扩散系数）：whole 加 nojump 即可，居中与否不影响均方位移。配方三（脂膜双层）：法向居中到双层中间平面——参考组选磷原子，这样上下叶始终在各自的半盒里，7.2 的厚度与序参数分析全靠这一步打底。

图 5-1 回卷与居中：一帧轨迹的三种形态

两个深坑：nojump 的起点与各向异性盒子

nojump 的实现是"让相邻帧位移取最小镜像"，因此它有一个隐蔽前提：处理必须从轨迹的第一帧开始。若你只截取后 50 ns 做 nojump，前 50 ns 累积的跨盒位移会在截取起点处被错误地"弹回"，原子被强制拉回第一帧的位置附近——扩散系数瞬间变成零附近的无意义小数。规矩：nojump 永远对完整轨迹做，做完了再按时间段截取。

第二个坑是各向异性盒子（脂膜模拟的常态，法向高度大于面内边长）：trjconv 的某些组合在非立方盒子上会把分子放进"错误的周期性盒子表示"。保险做法是处理脂膜时始终显式 -center 且参考组选双层中部（磷原子组），并抽查几帧可视化确认上下叶没有互换——上下叶互换这个错误不会报错，只会让你的序参数符号相反（7.2 详述）。

预处理的验收动作

修完不等于修好，验收两步。第一步抽三帧可视化（首、中、尾），确认没有断裂分子、没有贴壁怪异构象。第二步跑一个最小统计自检——对修复前后各算一次聚合物回转半径或双层厚度，若两值系统性不同（差超过涨落），说明修复改变了物理，配方用错了。这两步花不掉五分钟，却挡得住"分析完美、输入有毒"的最冤枉错误。

演算：量化"这条轨迹病得多重"

预处理不必盲改，先量化病情。下面的脚本统计原始轨迹里"断裂分子"的规模：每帧统计每个分子的原子中落在"距质心超过半盒宽"的比例，比例高即该分子正被边界切开。轨迹读取用 MDAnalysis 一类库，核心判据逻辑独立展示：

import numpy as np

def broken_fraction(frame_pos, bonds, box):
    """frame_pos: N×3；bonds: M×2 键表；box: 盒边长。
    返回"至少有一对键接原子相距超过半盒宽"的分子占比（粗判断裂）。"""
    half = box / 2
    d = frame_pos[bonds[:, 0]] - frame_pos[bonds[:, 1]]
    d -= box * np.round(d / box)          # 最小镜像还原真实位移
    stretched = np.linalg.norm(d, axis=1) > half
    broken_atoms = np.unique(bonds[stretched])          # 被切开的原子
    # 粗略折算成分子数：每条被切开的链至少贡献 2 个端点原子
    return len(broken_atoms) / frame_pos.shape[0]

# 演示：一条 60 格点的链放在 4 nm 盒子里，人为跨边界
chain = np.cumsum(np.full((60, 3), 0.13), axis=0)       # 直线链，长 7.7 nm
chain %= 4.0                                            # 折进盒内，模拟 PBC 写盘效果
bonds = np.array([[i, i + 1] for i in range(59)])
frac = broken_fraction(chain, bonds, 4.0)
print(f"被切开原子的占比：{frac:.2f}（0 说明该分子完好）")
​

对真实轨迹逐帧跑这个量，你会得到一条"病情曲线"：断裂帧占比几成、集中在哪些时段。它回答了两个实际问题——配方选多重的（断裂常见就必须 whole + nojump）、修复后要不要复核（修复后该占比必须归零，否则配方没起效）。

两个高频疑问

轨迹已经是 .xtc 压缩精度了，修复会不会放大误差？ 不会。whole 与 nojump 是整分子的刚体重排与最小镜像选择，不引入新的量化误差；0.001 nm 的坐标量化在键长尺度（0.1 nm 级）上无害。唯一要避免的是"反复转码"——修复链路上每多一次 trjconv 都是一次重写，配方一次写好、一条管道跑完，别修十轮。

分析时该用哪一版轨迹发给合作者？ 发修复版加它配套的 tpr，并附一句处理命令记录（哪几条 trjconv、什么顺序、参考组是谁）。"原始未修"版自己留档即可——合作者拿原始轨迹自己乱修，是合作破裂的经典导火索。

trjconv 常用旗标组合速查

把本节散落的命令收拢成一张卡，用时直接抄：

末行值得强调：先修复、再挑组——很多人为省事直接对含水轨迹做 whole 后切溶质，再把半成品当修复版用；切组会丢掉做 whole 所需的邻居信息，顺序反了就修不回去了。

速记：原始轨迹的断裂与跳动是几何病不是物理病；whole 修断裂、nojump 修跳动、mol 配 center 管完整性；扩散用前两样、双层补法向居中；nojump 必须从第一帧连续处理；修完抽帧目检加最小统计复核。几何修好了，下一节开始问第一个正经问题：体系稳定吗。