在分子动力学模拟的后处理中，蛋白‑蛋白双链复合物体系的轨迹分析往往比单链蛋白或蛋白‑小分子体系更易出现隐藏错误。笔者在分析一组蛋白双链复合物的100 ns的轨迹时，遇到了自定义分组导致蛋白结构残基大片缺失、校正后轨迹中两条链周期性“飞动”、以及骨架RMSD异常剧烈跳跃扰动等这一系列的现象。经实验与排查，最终发现到两个独立的问题根源：GROMACS内部残基索引与原始PDB残基编号的混淆，以及周期性矫正脚本中步骤顺序和输入轨迹的错误。希望这个文章能够对分析双链体系的朋友有帮助。

1. 背景与初始问题

我们研究对象为包含A、B两条蛋白链的复合物，蛋白残基总数372，体系总原子数约6.6万。模拟已用GROMACS完成，拥有标准轨迹文件 md.xtc 和拓扑 md.tpr。分析目标为获得整体蛋白、A链和B链的骨架RMSD、RMSF、回旋半径，以及链间氢键、盐桥、DCCM和自由能形貌图。我们参照成熟的蛋白‑小分子复合物分析脚本，编写了双链体系自动分析脚本，并通过 gmx make_ndx 创建自定义索引以区分两条链。

脚本运行后，三个异常现象同时出现，并且相互矛盾：

症状	表现
结构缺失	用自定义 Chain_AB 组输出轨迹PDB，残基335‑520区域完全不可见；若改用系统默认 Protein 组输出，结构完整
轨迹飞动	校正后轨迹动画中，两条链在某时刻突然分离、飞向盒子两端，随即弹回原位，仿佛“镜像瞬移”
RMSD剧烈跳跃	A链骨架RMSD在约10 ns内从0.5 nm骤升至2.5 nm以上并持续波动，但肉眼观察轨迹时复合物结合稳定，未见解离

1.1 结构缺失的问题线索

我们首先尝试用自定义的索引组从周期性校正轨迹中提取代表性构象。执行命令：

gmx trjconv -s md.tpr -f md_center.xtc -n index_AB.ndx -o test.pdb -dump 50000

在交互界面中选择组 19（即 Chain_AB），程序输出：

Selected 19: 'Chain_AB'

然而在 PyMOL 中打开生成的 PDB 文件，发现蛋白质主链在中间位置出现大面积断裂，残基编号约 335 至 520 的区域完全不可见，仅余两个孤立的片段。当我们改用系统默认的 Protein 组（组 1）输出同一帧时，结构却完整无缺。进一步检查索引原子数，发现 Chain_AB 仅包含 3449 个原子，而 Protein 组有 5724 个原子，二者相差 2275 个原子，恰好对应那些“消失”的残基。（详细的错误建组流程在后面）

箭头所指区域就是我们当时出现残缺的氨基酸，并且仔细观察可以发现，体系中出现了一部分少量的水分子，但是我们选择的组别是Chain_AB，不应该包含水分子，这是疑点之一，并且蛋白质结构，红色加蓝色部分是我们完整的模型，而我们Chain_AB导出的构象仅仅含有蓝色部分结构和水分子。

分组创建的原始错误逻辑

为了理解这一问题的根源，有必要回溯我们最初创建索引分组的完整过程。由于原始PDB文件中明确标注了链标识符——A链为残基261‑558，B链为残基1‑74，我们很自然地认为GROMACS内部保留了这些原始编号。基于这一错误假设，我们使用以下命令创建了 index_AB.ndx：

gmx make_ndx -f md.tpr -o index_AB.ndx << EOF
ri 261-558
name 17 Chain_A
ri 1-74
name 18 Chain_B
17 | 18
name 19 Chain_AB
"Chain_A" & Backbone
name 20 Chain_A_bb
"Chain_B" & Backbone
name 21 Chain_B_bb
"Chain_AB" & Backbone
name 22 Chain_AB_bb
q
EOF

在交互执行过程中，gmx make_ndx 返回了以下关键信息：

Reading structure file
Reading file md.tpr, VERSION 2024.2 (single precision)
Analysing residue names:
There are:   372    Protein residues
There are: 20091      Water residues
There are:    13        Ion residues
Analysing Protein...

  0 System              : 66010 atoms
  1 Protein             :  5724 atoms
  2 Protein-H           :  2925 atoms
  3 C-alpha             :   372 atoms
  4 Backbone            :  1116 atoms
  ...

> ri 261-558
Found 2285 atoms with resind.+1 in range 261-558

> ri 1-74
Found 1166 atoms with resind.+1 in range 1-74

> 17 | 18
Merged two groups with OR: 2285 1166 -> 3451

此时的输出已经出现了第一个明显的不协调：Chain_A（ri 261-558）包含了2285个原子，Chain_B（ri 1-74）包含了1166个原子，合并后的 Chain_AB 总原子数为3451。然而系统默认的 Protein 组显示有5724个原子——这意味着我们创建的 Chain_AB 组竟然缺少了2273个蛋白原子，缺失比例高达整个蛋白的40%。

接下来的布尔操作又引发了连锁问题：

> "Chain_A" & Backbone
Syntax error: "Backbone"

由于 make_ndx 不能直接使用用户自定义的带引号名称进行交集操作，我们只得改用组号：

> 17 & 4
Merged two groups with AND: 2285 1116 -> 222

这里 Chain_A_bb 仅含222个主链原子，而一个298残基的链，其主链原子数（N、Cα、C、O）应为298 × 4 = 1192个。222个原子的巨大偏差进一步暗示了分组选出的残基集合与预期严重不符。

从怀疑到验证

当我们将这个残缺不全的 Chain_AB 组用于后续的周期性矫正和轨迹输出时，所得到的结构自然是断裂的。但当时我们并不清楚问题是出在残基编号映射上，而是一度怀疑周期性边界条件处理不当导致了分子断裂。经过几轮实验排除了PBC的嫌疑后，我们决定从根本上检查索引文件本身。

在 make_ndx 交互环境中执行 l 命令，去输出md.tpr的全部残基，程序输出了以下信息：

> l

   1 ASN     2 PRO     3 MET     4 HIS     5 LYS     6 GLU     7 ILE     8 SER
   9 GLN    10 ARG    11 SER    12 THR    13 ALA    14 THR    15 MET    16 TYR
  17 ILE    18 ILE    19 GLY    20 GLY    21 TYR    22 TYR    23 TRP    24 HIS
  25 PRO    26 LEU    27 SER    28 GLU    29 VAL    30 HIS    31 ILE    32 TRP
  33 ASP    34 PRO    35 LEU    36 THR    37 ASN    38 VAL    39 TRP    40 ILE
  41 GLN    42 GLY    43 ALA    44 GLU    45 ILE    46 PRO    47 ASP    48 TYR
  49 THR    50 ARG    51 GLU    52 SER    53 TYR    54 GLY    55 VAL    56 THR
  57 CYS    58 LEU    59 GLY    60 PRO    61 ASN    62 ILE    63 TYR    64 VAL
  65 THR    66 GLY    67 GLY    68 TYR    69 ARG    70 THR    71 ASP    72 ASN
  73 ILE    74 GLU    75 ALA    76 LEU    77 ASP    78 THR    79 VAL    80 TRP
  81 ILE    82 TYR    83 ASN    84 SER    85 GLU    86 SER    87 ASP    88 GLU
  89 TRP    90 THR    91 GLU    92 GLY    93 LEU    94 PRO    95 MET    96 LEU
  97 ASN    98 ALA    99 ARG   100 TYR   101 TYR   102 HIS   103 CYS   104 ALA
  ...
 370 ILE   371 GLY   372 VAL
 373 - 20463 SOL     20464 - 20476 NA

从这份完整的残基列表中可以看到，蛋白质残基的 GROMACS 内部索引编号为 1 至 372，且是连续无缺失的。仔细观察残基名称的排列规律：A 链（原始编号 261‑558）在 pdb2gmx 时被首先写入拓扑，因此占据了索引 1‑298 的位置；而 B 链（原始编号 1‑74）紧随其后，位于索引 299‑372。

注意：A链在输出的列表前排布，而B链在列表后排布，可以手动对照，切忌勿要倒置，不然结果数据会出现逻辑性错误，但肉眼很难观察出。

换言之，GROMACS 在 pdb2gmx 步骤中已经将整个蛋白的残基按照在 PDB 文件中出现的物理顺序重新编号为 1‑372 的连续序列，而原始 PDB 中的残基编号虽然被记录在拓扑文件中，却不能被 make_ndx 的 ri 命令所使用。

反观我们当时的命令 ri 261-558，程序忠实地返回了 “Found 2285 atoms with resind.+1 in range 261‑558”——它确实找到了索引 261 至 558 的约 298 个残基，但这些残基已经远超出了 A 链的范围，大部分落入了 B 链甚至水分子、离子区域，选出的原子自然不属于目标 A 链，这就是为什么我们上面出现了一部分的水分子。同理，ri 1-74 虽然选出了 A 链的前 74 个残基，但依然不是完整的 B 链。两组选出的原子集合与真正的 A、B 链毫无对应关系，合并后便造成了视觉上的结构断裂和原子数严重缺失。
v的提取验证**

选取 RMSD 平稳的时间段（90‑100 ns），从校正轨迹中提取动画。播放发现飞动完全消失，蛋白保持稳定构象。这进一步证实飞动与 RMSD 跳跃发生在同一时间窗口，且二者均与矫正处理密切相关。

通过这三阶段实验，我们排除了原始模拟数据本身异常的假设，将问题根源锁定在两个独立的环节：索引文件创建时的残基编号错误和周期性矫正脚本中步骤顺序与输入轨迹的逻辑错误。后续章节将逐一剖析这两个技术陷阱及其修复方法。

1.2 轨迹飞动与RMSD剧烈跳跃：独立的症状，相同的指向

在修正索引分组之后，我们曾预期所有异常会一并消失。但实际情况远非如此——RMSD曲线依然剧烈跳跃，而更诡异的是，轨迹动画中出现了两条链的“瞬移”现象。这两个问题几乎同时暴露，但它们的表现方式和时间特征并不相同也不具有耦合性。

我们首先从校正后的轨迹中提取了70-80 ns的多帧PDB生成动画：

gmx trjconv -s md.tpr -f md_center.xtc -n index_AB.ndx -o anim.pdb -b 70000 -e 80000 -dt 100 -pbc mol

选择组19（Chain_AB）后，程序输出确认：

Selected 19: 'Chain_AB'

在PyMOL中播放该动画时，令人震惊的画面出现了。在约75 ns附近，A、B两条链突然朝相反方向飞离，间距瞬间拉开至接近盒子尺寸，随后在短短几帧之内又弹回原位，恢复紧密接触。这一“飞动—弹回”的瞬态过程在我们截取的70-80 ns区间内反复出现了数次，像是蛋白在盒子内“瞬移”。然而，在其余部分时间内，复合物构象表现得相当稳定，链间结合紧密，未见任何解离或大尺度构象变化。

校正后轨迹动画截图对比——分别展示飞动瞬间ABC三帧的构象，明显在A帧时正常，下一帧B帧时突然瞬移，C帧又恢复正常。

与此同时，我们使用以下两条命令分别从原始轨迹和校正后轨迹计算了A链骨架RMSD：

gmx rms -s md.tpr -f md.xtc -n index_AB.ndx -o rmsd_raw.xvg -tu ns
gmx rms -s md.tpr -f md_center.xtc -n index_AB.ndx -o rmsd_center.xvg -tu ns

原始轨迹与校正轨迹的A链骨架RMSD曲线对比。横轴为时间（0-100 ns），纵轴为RMSD（nm）。两条曲线几乎完全重合，均显示0-100 ns出现RMSD的扰动，数值骤升至2.5 nm以上的持续平台段。]

两条RMSD曲线几乎完全重合，展现出相同的异常模式：部分时间内，RMSD平稳维持在0.4-0.6 nm，表明蛋白质主链构象在这些时间段内稳定；然而到了例如80-90 ns附近，RMSD突然飙升至2.5 nm以上，之后在1.5-3.0 nm之间持续剧烈波动，又再次回到平稳基线。在正常的蛋白MD模拟中，2.5 nm量级的骨架RMSD变化通常意味着蛋白发生了完全解折叠，或者两条链彻底解离。但我们在PyMOL中逐帧检查这段轨迹时，复合物明明始终保持紧密结合，未出现任何显著的构象重排或链间分离——飞动只是瞬态的“弹跳”，并非真实的解离。

这里有一个关键的观察：飞动现象出现在约大部分偶然时间帧的区间，而RMSD的剧烈跳跃则最一定的长时间周期内才开始，二者在时间单位上并不严格吻合。飞动发生的时段，RMSD尚处于平稳期；而RMSD开始跳跃后，飞动反而减弱了。这意味着飞动和RMSD跳跃并不是同一个物理效应的两种表现，而是两个独立的异常信号。但它们有一个共同的指向——都与周期性矫正处理密切相关：飞动是矫正后轨迹才有的，而RMSD跳跃则在原始轨迹中就已存在，说明原始轨迹的坐标记录本身在后期存在某种问题，而矫正步骤非但没有修复它，反而在另一个时间窗口引入了新的畸变。

三阶段对照实验

修正索引之后，Chain_AB 组的原子数已与 Protein 组完全一致，从校正轨迹中输出的蛋白结构也不再残缺。然而，RMSD曲线依然剧烈跳跃，轨迹动画中仍然出现两链“瞬移”的诡异现象。这两个问题几乎同时暴露，但它们的表现方式和时间特征并不相同，且都指向了周期性矫正这一共同环节。为了厘清根源，我们分三个阶段进行了系统的对照实验。

第一阶段：原始轨迹中不同分组的一致性验证

在怀疑矫正流程之前，我们必须首先确认目前使用的索引本身是否会在计算中引入错误。最直接的验证方式，是在未经任何矫正处理的原始轨迹 md.xtc 上，分别用系统默认的 Protein 组和我们自定义的 Chain_AB 组提取同一时刻的结构并计算RMSD。

首先，从原始轨迹中提取单帧结构（以t = 50 ns为例）进行比较：

# 使用Protein组输出
gmx trjconv -s md.tpr -f md.xtc -o raw_protein.pdb -dump 50000 -pbc mol <<EOF
1
EOF

# 使用Chain_AB组输出
gmx trjconv -s md.tpr -f md.xtc -n index_AB.ndx -o raw_chainAB.pdb -dump 50000 -pbc mol <<EOF
19
EOF

在PyMOL中将两个PDB文件叠合，结果完全重叠，原子坐标没有任何偏差。这初步说明，Chain_AB 组选出的原子在原始轨迹中与 Protein 组对应的原子完全一致。

接下来，我们用同样的两组索引分别计算A链骨架RMSD，以验证数值计算结果是否也一致：

# 使用Protein组（组1）做拟合，Chain_A_bb（组21）做计算
gmx rms -s md.tpr -f md.xtc -n index_AB.ndx -o rmsd_raw_viaProtein.xvg -tu ns <<EOF
1
21
EOF

# 使用Chain_AB组（组19）做拟合，Chain_A_bb（组21）做计算
gmx rms -s md.tpr -f md.xtc -n index_AB.ndx -o rmsd_raw_viaChainAB.xvg -tu ns <<EOF
19
21
EOF

两条RMSD曲线几乎完全重合，表明无论用 Protein 组还是 Chain_AB 组作为拟合参考，RMSD计算结果都一致。这确认了索引本身的计算功能没有问题。

然而，正是在这条RMSD曲线上，我们看到了一个棘手的事实：即使在原始轨迹上，A链骨架RMSD从约80 ns开始仍然出现了剧烈跳跃，从0.5 nm飙升至2.5 nm以上。这说明，飞动现象虽未在原始轨迹md.xtc的可视化中出现，但RMSD异常在原始数据中就已存在。

原始轨迹中分别使用Protein组和Chain_AB组计算的A链RMSD，两条曲线完全重合，均显示剧烈跳跃。

第二阶段：错误矫正轨迹中的分组对比与飞动定位

第一阶段的结果证实了我们的index_AB.ndx索引的可靠性，也表明原始轨迹中的RMSD跳跃尚未解决。接下来，我们将目光转向那条由旧版错误脚本生成的校正轨迹 md_center.xtc。在这条轨迹上，我们重复了与第一阶段类似的分组对比，但更重要的是直接观察轨迹动画本身。

首先，同样提取t = 50 ns的单帧结构进行比较：

# 使用Protein组输出
gmx trjconv -s md.tpr -f md_center.xtc -o center_protein.pdb -dump 50000 -pbc mol <<EOF
1
EOF

# 使用Chain_AB组输出
gmx trjconv -s md.tpr -f md_center.xtc -n index_AB.ndx -o center_chainAB.pdb -dump 50000 -pbc mol <<EOF
19
EOF

叠合结果显示，两者在重合原子区域的坐标仍然一致。这说明，即使经过了错误的矫正处理，Chain_AB 组与 Protein 组之间的相对关系并未被破坏，飞动不是由于索引选取了错误的原子或原子顺序混乱造成的，这一点很重要，我们排除了索引造成的差错。

接下来，我们直接从这条校正轨迹中提取70-80 ns的多帧动画，以检验飞动现象：

gmx trjconv -s md.tpr -f md_center.xtc -n index_AB.ndx -o anim_center.pdb -b 70000 -e 80000 -dt 100 -pbc mol <<EOF
19
EOF

在PyMOL中播放该动画时，预期的结果出现了。在约75 ns附近帧左右，A、B两条链突然朝相反方向飞离，间距瞬间拉开至接近盒子尺寸（约10 nm），随后在短短几帧之内又弹回原位，恢复紧密接触。这一过程在70-80 ns区间内反复出现，但在其余时间段，复合物构象表现得相当稳定。

错误校正轨迹的动画截图——图为飞动瞬间（约75 ns），两链分离至盒子两端；。

与此同时，我们计算了这条校正轨迹的A链RMSD，发现其与原始轨迹的RMSD曲线几乎重合，RMSD的跳跃扰动特征完全一致。也就是说，错误的矫正流程既没有修复原始轨迹中已有的RMSD异常，还在特定时间窗口引入了轨迹可视化肉眼可见的飞动假象。飞动和RMSD跳跃虽然在时间上并不严格同步，但它们的共同点是：都与周期性矫正的处理方式密切相关。

错误校正轨迹与原始轨迹的A链RMSD，两条曲线几乎重合，均发生剧烈跳跃扰动。]

第三阶段：更改矫正流程和验证

前两个阶段的实验形成了一个清晰的逻辑闭环：原始轨迹中，不同索引分组的结果一致，说明索引正确，但RMSD本身存在跳跃；错误矫正后的轨迹中，分组对比仍然一致，说明飞动并非索引错误，而是矫正过程引入了新的全局坐标畸形变。那么，如果恢复标准的矫正顺序，是不是能否同时解决这两个问题？

我们回到此前我们成熟的蛋白-小分子复合物分析脚本中采用周期性矫正的“标准三步法”（这一部分内容可以前去这篇我们之前的博文查看CSDN文章链接、博客链接）：

1. 从原始轨迹 md.xtc 执行 -pbc nojump，输出 md_nojump.xtc；
2. 从原始轨迹 md.xtc 执行 -pbc cluster -center -ur compact，输出 md_DH.xtc；
3. 从 md_nojump.xtc 执行 -fit rot+trans -center，输出 md_center.xtc。

所有居中与拟合的参考组均使用系统默认的 Protein 组（组1），输出均为整个 System（组0）。

# 步骤1：去周期跳跃，输出全系统
gmx trjconv -s md.tpr -f md.xtc -o md_nojump.xtc -pbc nojump <<EOF
0
EOF
# 步骤2：从原始轨迹做团簇化、居中、紧凑化，输出全系统，这里我们并没有用这一个，在这里写这一个的原因在后面
gmx trjconv -s md.tpr -f md.xtc -o md_DH.xtc -pbc cluster -center -ur compact <<EOF
1
1
0
EOF

# 步骤2：从 nojump 轨迹做旋转平移拟合 + 居中，输出全系统
gmx trjconv -s md.tpr -f md_nojump.xtc -o md_center_new.xtc -fit rot+trans -center <<EOF
1
1
0
EOF

用这条新生成的 md_center_new.xtc 重新计算A链骨架RMSD：

gmx rms -s md.tpr -f md_center_new.xtc -n index_AB.ndx -o rmsd_new.xvg -tu ns <<EOF
21
21
EOF

结果令人振奋：RMSD曲线在整个0-100 ns范围内都保持平滑，波动范围稳定在0.3-0.45 nm，异常跳跃扰动完全消失。同时，从新轨迹中提取70-80 ns的动画也不再出现任何飞动，蛋白全程相对正常。

[正确矫正后的A链RMSD曲线，显示全程平滑，无异常跳跃。

将旧版错误脚本的矫正命令与标准三步法进行对比，可以清楚地看到两处致命差异。其一，旧脚本第一步使用了 -pbc mol -center，而标准流程是 -pbc nojump。-pbc mol 只能保证分子完整性，无法消除原子跨盒子时的坐标跳变，这一步的失当直接导致原始轨迹中的RMSD跳跃无法被修复。其二，旧脚本将团簇化步骤的输入从原始轨迹错误地改为了已处理过的 md_nojump.xtc，而 cluster 算法需要原始的分子拓扑信息才能正确团簇化。正是这一逻辑混乱，在局部区域产生了坐标失真，表现为轨迹中的“瞬移”飞动。当矫正顺序恢复标准三步法后，两个问题被一并根除。

2. 根源一：残基编号与内部索引的错位

错误表现：原始 PDB 中 A 链编号为 261‑558，B 链编号为 1‑74。最初使用 ri 261‑558 和 ri 1‑74 创建链组，得到的 Chain_AB 仅含 3449 个原子，而系统默认的 Protein 组含 5724 个原子，相差超过 2000 个原子。用该组输出校正轨迹 PDB 时，蛋白中间区域出现大片断裂，残基 335‑520 区域完全不可见，同时结构中出现了不应存在的水分子。而改用 Protein 组输出则结构完整。这些现象直接表明索引选出的原子集合与真实蛋白链严重不符。

原因：GROMACS 在 pdb2gmx 步骤中会按照残基在 PDB 文件中出现的物理顺序，将所有蛋白残基重新编号为从 1 开始的连续整数，这个内部序号称为残基索引。make_ndx 的 ri 命令使用的正是残基索引，而非原始 PDB 中的残基编号。在本体系中，A 链（原始编号 261‑558）被首先写入拓扑，占用了残基索引 1‑298；B 链（原始编号 1‑74）紧随其后，位于索引 299‑372。因此，执行 ri 261‑558 时，程序实际选取的是索引 261 至 558 的残基，而这些索引位置已经远远超出 A 链的范围，大部分落入了 B 链以及随后的水分子和离子区域。ri 1‑74 则选出了 A 链的前 74 个残基，而非完整的 B 链。两组选出的原子集合与真正的 A、B 链完全错位，导致合并后的 Chain_AB 原子数异常、结构残缺，且混入了非蛋白原子。

检查：在 make_ndx 交互环境中执行 l 命令列出全部残基，可获得残基索引号、残基名称和原始编号的映射表。根据列表中残基名称的排列规律和原始编号的跳变，能准确判断出每条链对应的真实索引范围。正确的创建命令为：

gmx make_ndx -f md.tpr -o index_AB.ndx << EOF
ri 1-298
name 17 Chain_A
ri 299-372
name 18 Chain_B
17 | 18
name 19 Chain_AB
17 & 4
name 20 Chain_A_bb
18 & 4
name 21 Chain_B_bb
19 & 4
name 22 Chain_AB_bb
q
EOF

修复后，Chain_AB 原子数与 Protein 组完全一致，输出的轨迹结构恢复完整。

注意：在任何体系上创建自定义索引前，必须先用 l 命令确认真实的残基索引映射关系，绝不可假设原始 PDB 编号在 GROMACS 内部保持不变。同时需要注意链在输出列表中的先后顺序，避免 A、B 链倒置导致的逻辑性数据错误。

3. 根源二：周期性矫正步骤的顺序错误

错误表现：即使修复了索引，RMSD 曲线仍然剧烈跳跃，且校正后轨迹中依然存在蛋白链飞动的现象。然而，直接从原始轨迹 md.xtc 提取的动画完全稳定，飞动仅出现在经周期性校正的轨迹中；同时，原始轨迹与校正轨迹的 RMSD 曲线几乎完全重合，均出现剧烈跳变。这表明原始轨迹的坐标记录本身在后期即已存在问题，而矫正步骤非但未能修复，反而在特定时间窗口引入了新的全局坐标畸变。

矫正逻辑对比：此前成熟的蛋白‑小分子复合物分析脚本采用标准三步法，而双链体系脚本被意外错误改写。两套命令序逻辑顺序的关键差异如下：

步骤	正确脚本（标准三步法）	错误脚本（双链初版）
第一步	-pbc nojump，从 md.xtc 读，输出 md_nojump.xtc	-pbc mol -center，从 md.xtc 读，输出 md_nojump.xtc
第二步	-pbc cluster -center -ur compact，从 md.xtc 读，输出 md_DH.xtc	-pbc mol -center，从 md_nojump.xtc 读，输出 md_DH.xtc
第三步	-fit rot+trans -center，从 md_nojump.xtc 读，输出 md_center.xtc	-fit rot+trans -center，从 md_nojump.xtc 读，输出 md_center.xtc

错误的脚本包含两个致命缺陷。其一，用 -pbc mol 取代了 -pbc nojump。-pbc mol 只能保证单分子完整，并不能消除原子跨越周期边界时产生的坐标跳变，导致拟合时坐标仍不连续，原始轨迹中本应被修复的 RMSD 跳跃因此得以存留。其二，团簇化步骤的输入从原始轨迹错误地改为已处理过的 md_nojump.xtc。cluster 算法需要原始的分子拓扑信息才能正确实现团簇划分和居中，用已去跳变的轨迹作为输入会破坏这一前提，致使部分原子的居中参考发生偏离，在局部区域产生坐标失真，这正是轨迹动画中出现链瞬移飞动的来源。

# 第一步：去周期跳跃，输出全系统
gmx trjconv -s md.tpr -f md.xtc -o md_nojump.xtc -pbc nojump <<EOF
0
EOF

# 第二步：从原始轨迹做团簇化、居中、紧凑化，输出全系统
gmx trjconv -s md.tpr -f md.xtc -o md_DH.xtc -pbc cluster -center -ur compact <<EOF
1
1
0
EOF

# 第三步：从 nojump 轨迹做旋转平移拟合 + 居中，输出全系统
gmx trjconv -s md.tpr -f md_nojump.xtc -o md_center.xtc -fit rot+trans -center <<EOF
1
1
0
EOF

修正后，RMSD 曲线在整个 100 ns 范围内恢复平滑，波动范围稳定在 0.3‑0.45 nm，飞动现象完全消失。这一流程的恢复也再次印证了周期性矫正顺序在 GROMACS 轨迹后处理中的核心地位——输入轨迹的选择与操作符的匹配必须严格遵循物理逻辑，任何偏离都可能引入难以察觉但后果严重的数据偏差。

4. 通用稳健分析流程

最后基于上述排错经验，笔者整出一套适用于任意多链蛋白复合物MD分析的高鲁棒性流程，核心规则如下表：

阶段	关键操作	铁律
索引创建	make_ndx 中用 l 列出残基，确定链的残基索引范围	绝不用原始PDB编号直接 ri，必须验证索引映射
周期性矫正	严格三步骤：① -pbc nojump (读 md.xtc) → ② -pbc cluster -center -ur compact (读 md.xtc) → ③ -fit rot+trans -center (读 md_nojump.xtc)	输出全系统（组0），居中/拟合参考用 Protein（组1）
分析计算	所有 rms、rmsf、gyrate、hbond 等命令通过 -n index.ndx 选择子组	不要用索引直接输出分析子集轨迹，保持轨迹完整
子集轨迹生成	仅当外部程序要求匹配原子数时（如DCCM），先用 convert-tpr 生成匹配拓扑，再从全系统轨迹提取	确保拓扑与轨迹原子数一致

遵循该流程，我们已在当前两个不同体系上获得正常数据结果、完整结构和可重复的分析结果。

5. 最后

在MD后数据处理中，一个残基编号的误用、一条命令顺序的颠倒，都可能导致面目全非的分析结果。在计算化学公社论坛笔者也发现许多类似关于多链蛋白数据分析的问题，这里也希望这一经验能帮助大家在分析多链蛋白体系时少走弯路，确保计算结果真实反映模拟物理过程。

个人博客地址：