GROMACS 距离约束与 NMR 输出分析
定义距离约束、将相关量写入EDR，提取瞬时距离、平均违反量与约束能，并对照真实坐标。
阅读本篇目录
下载本组输入、真实参考结果和绘图脚本，解压后进入资料包目录。inputs 是可以修改的输入，reference 是鸢算科技集群已经计算得到的参考结果。每次运行在 work 下建立新目录，路径写入 LAST-RUN.txt。计算请在已经分配的资源中执行；整组可用 sbatch run.sbatch 提交。
本例要得到什么#
使用柔性三原子水分子说明距离约束的操作：在两个 H 原子之间定义上限、运行短 MD、从 EDR 提取距离与违反量，并与 TRR 坐标独立比较。它是约束算法教学模型，不是实验 NOE 数据拟合。
bash run.sh nmr
如希望同时生成独立坐标核对和所有图，可用 bash run.sh all；脚本中的分析使用集群现有的 MDAnalysis 2.9.0 环境，不要求客户自行写 Python。
1. 约束写在拓扑的哪里#
完整 inputs/nmr/topol.top：
; Explicit flexible TIP3P-type teaching monomer; nm, kJ/mol, e units.
[ defaults ]
1 2 yes 0.5 0.8333333333
[ atomtypes ]
OW 8 15.9994 0 A 0.315061 0.636386
HW 1 1.008 0 A 0 0
[ moleculetype ]
SOL 2
[ atoms ]
1 OW 1 SOL OW 1 -0.834 15.9994
2 HW 1 SOL HW1 1 0.417 1.008
3 HW 1 SOL HW2 1 0.417 1.008
[ bonds ]
1 2 1 0.09572 502416
1 3 1 0.09572 502416
[ angles ]
2 1 3 1 104.52 628.02
[ distance_restraints ]
; ai aj funct index type low up1 up2 fac
2 3 1 0 1 0.0 0.149 0.18 1.0
[ system ]
Flexible water monomer
[ molecules ]
SOL 1
两个原子编号为 2、3，单位 nm。up1=0.149 表示上限开始，up2=0.18 决定上限之外势函数的分段；本例采样距离低于 up2。这属于有势能惩罚的 distance restraint，不是将距离固定不变的 LINCS/SETTLE constraint。
2. 准备并运行 MD#
完整 restraints.mdp：
cutoff-scheme = Verlet
nstlist = 1
rlist = 1.0
coulombtype = Cut-off
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
pbc = xyz
constraints = none
verlet-buffer-tolerance = -1
nstenergy = 1
nstcalcenergy = 1
integrator = md
dt = 0.0002
nsteps = 200
nstxout = 1
nstvout = 1
comm-mode = Linear
nstcomm = 1
gen-vel = yes
gen-temp = 300
gen-seed = 20261028
disre = simple
disre-weighting = conservative
disre-fc = 1000
disre-tau = 0
nstdisreout = 1
nstdisreout=1 将距离约束相关量逐步写入 EDR；普通能量项存在不代表所有 NMR 分析需要的量都已保存。disre-tau=0 关闭时间平均，本例更容易从瞬时坐标核对。
module purge
module load gromacs/2026.4-p100-threadmpi
mkdir nmr-manual
cp inputs/nmr/* nmr-manual/
cd nmr-manual
gmx grompp -f restraints.mdp -c start.gro -p topol.top -o restraints.tpr
gmx mdrun -deffnm restraints -ntmpi 1 -ntomp 2 -nb cpu -bonded cpu -update cpu
本例使用 CPU 执行路径，200 步×0.2 fs=40 fs，输出包含初始时刻在内的 201 帧。
3. 从 EDR 提取约束数据#
printf '1\n0\n' | gmx nmr -f restraints.edr -s restraints.tpr -viol nmr-violations.xvg -pairs nmr-pairs.xvg
printf 'Dis.-Rest.\n0\n' | gmx energy -f restraints.edr -o restraint-energy.xvg
本例 NMR 提示中选择 1 分析该约束，再用 0 结束选择。实际菜单以程序显示为准；其他拓扑可能有多个约束。nmr-pairs.xvg 保存按帧的成对距离相关数据，nmr-violations.xvg 与日志给出平均相关结果。
4. 对照真实数值#
| 量 | 本例结果 |
|---|---|
| H—H 瞬时距离范围 | 0.1430646–0.1578624 nm |
<r⁻³>⁻¹/³ 有效距离 |
0.1501794 nm |
| 有效距离相对上限的平均违反量 | 0.00117936 nm |
| 平均约束能 | 0.0101703 kJ/mol |
本例只有一个距离对，平均分析使用 <r⁻³>⁻¹/³，不要把它当作普通距离算术平均。瞬时距离超过 0.149 nm 并不意味着整体轨迹报错；约束势会产生恢复作用。
当前距离区间的独立能量公式为：
violation = max(rHH - 0.149, 0)
V = 0.5 × 1000 × violation² # kJ/mol，距离用 nm
资料包从 TRR 独立逐帧读坐标，能量与 GROMACS 原生能量最大差约 5.4e-7 kJ/mol。原生距离文件打印位数有限，逐帧距离最大差小于 5e-5 nm。这些是短流程的数值对照，不是约束参数已适合任意真实体系。
采用实验距离约束时还需核对原子映射、等价原子、约束分组、时间平均和力常数；不应把本例水分子的原子编号及阈值直接套给蛋白。
参数见 GROMACS nmr 官方说明。
从已保存轨迹逐帧检查#
已有TPR/TRR时，可按轨迹距离约束分析使用独立 gmx-disre-compat。该篇包括单距离对、周期、宽阈值与含等价对的分组控制，并说明它与EDR平均分析的区别。