GROMACS 温度复制交换与轨迹拆分

下载输入、实际结果和绘图资料包。解压后进入资料包目录，inputs 是可修改输入，reference 是鸢算科技集群的参考结果。脚本每次在 work 下新建目录，结果位置写入 LAST-RUN.txt。计算请放在已分配的资源内执行。

先理解本篇要得到什么#

温度复制交换（temperature REMD）同时运行多份相同体系，在指定时刻尝试交换构型所在的温度。它用于帮助跨越能垒；是否改善采样，需要更长轨迹、温度梯度和目标体系共同判断。

本例是 216 个 SPC/E 水分子、648 个原子、300 K 与 310 K、各 4 ps。完整流程是准备两套相同拓扑和坐标 → 分别生成 TPR → 两个 MPI 进程协同运行 → 提取能量并核对交换 → 拆出随构型行走的轨迹。这个短算例便于掌握操作，不用于给出收敛自由能。

1. 准备两套目录#

直接运行整组：

sbatch run.sbatch

只做本篇时，在分配了 2 个 MPI 进程、每进程 1 个线程 的资源内运行：

bash run.sh replica
cat LAST-RUN.txt

手动操作可以从包内输入开始：

module purge
module load gromacs/2026.4-p100-mpi-plumed
mkdir my-remd
cp -r inputs/replica-300 inputs/replica-310 my-remd/
cd my-remd

两套目录都有 start.gro、topol.top、md.mdp。SPC/E 水的拓扑与初始结构已经配好，不需要另行生成结构。topol.top 中使用相同分子类型与数量；不要让两份体系的原子顺序、约束或自由度不同。

300 K 的完整 md.mdp：

integrator = md
dt = 0.001
nsteps = 4000
continuation = no
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 0.8
rvdw = 0.8
pbc = xyz
tcoupl = v-rescale
tc-grps = System
tau-t = 0.1
ref-t = 300
pcoupl = no
gen-vel = yes
gen-temp = 300
gen-seed = 20261328
ld-seed = 20261329
constraints = none
nstenergy = 10
nstxout-compressed = 50
fourierspacing = 0.12

310 K 目录主要改变 ref-t、gen-temp，并使用各自的固定随机种子。MD 步长 0.001 ps 即 1 fs，4000 步为 4 ps；能量每 10 步输出，压缩轨迹每 50 步输出。温度由整个 System 耦合，不把单个水原子分成独立温控组。

2. 生成 TPR，再一起运行#

for temperature in 300 310; do
  cd replica-$temperature
  gmx_mpi grompp -f md.mdp -c start.gro -p topol.top -o md.tpr
  cd ..
done

export OMP_NUM_THREADS=1
mpirun -np 2 --bind-to none --mca pml ob1 --mca btl self,sm,tcp \
  gmx_mpi mdrun -multidir replica-300 replica-310 -deffnm md \
  -replex 100 -reseed -270409831 -ntomp 1 \
  -nb cpu -bonded cpu -update cpu -pme cpu -pin off

-multidir 把目录作为一个耦合任务运行，不能分别启动两个普通 mdrun 来代替。-replex 100 每 100 步进入一次交换周期；-reseed 固定交换随机数种子。本例使用 GROMACS 2026.4 与 OpenMPI 5.0.10，采用 CPU 计算路径；提交模板分配的 GPU 不参与该例计算。

不要把本例两套温度直接照搬为任意蛋白的正式温度梯度。真实体系要检查相邻能量分布重叠、交换接受率和 walker 往返情况。

3. 提取能量和检查交换#

for temperature in 300 310; do
  cd replica-$temperature
  printf 'Potential\nTemperature\nKinetic-En.\n0\n' | \
    gmx_mpi energy -f md.edr -o thermodynamics.xvg
  cd ..
done

grep -E 'Replica exchange at step|Repl ex|Repl pr' replica-300/md.log

每个温度目录得到 401 个能量时间点、81 帧坐标、末时刻 4 ps。各列含义应读取 XVG 的 legend，不能假定工具按输入选择顺序排列输出列。

以两个温度的同一时刻势能计算：

import math
R = 0.008314462618  # kJ/(mol K)
delta = (1/300 - 1/310) * (U310 - U300) / R
probability = min(1.0, math.exp(-delta))

U310、U300 单位为 kJ/mol。集群参考第 100 步的无量纲能量项约 −0.6842，概率为 1；第 300 步约 0.4883，概率约 0.61。概率大于零不代表一定接受。

本次有 39 个交换周期、20 次实际相邻对尝试、6 次接受，接受率 30%。两份复制的奇偶配对轮次不同，有些周期没有实际成对尝试，不能用 6/39 当接受率。不同硬件或浮点路径可能改变后续交换轨迹，重算应检查完整输出与概率公式，而非只匹配一个接受次数。

4. 区分温度轨迹与 walker 轨迹#

replica-300/md.xtc 是固定温度槽位的轨迹，交换后其中的构型来源会变化。要跟踪同一 walker 在两温度之间的行走，先生成映射，再用原生工具拆轨迹：

mkdir demux
cd demux
demux.pl ../replica-300/md.log
gmx_mpi trjcat -f ../replica-300/md.xtc ../replica-310/md.xtc \
  -demux replica_index.xvg -o walker.xtc -nosort

输出 0_walker.xtc 与 1_walker.xtc，各 81 帧。replica_temp.xvg 展示温度路径，replica_index.xvg 用于拆分。这里直接调用 demux.pl 日志，不加额外时间间隔参数。

在这套输出中，恰好位于交换时刻的坐标帧仍采用交换前来源；随后帧才采用新的对应关系。资料包已经用真实坐标逐帧核对原生拆分结果。

5. 最后拿到哪些结果#

• 两个温度槽位的 md.edr、md.log、md.xtc 和 thermodynamics.xvg。
• 两个 walker 的连续来源轨迹、温度路径和映射表。
• 逐次交换能量项、概率以及独立坐标核对结果。

重新生成包内图：

bash run.sh plots

若要进一步做偏置坐标采样，见 AWH 距离采样与结果提取。AWH 与本篇 REMD 是不同方法，不需要为了做本例启用 PLUMED。