GROMACS WHAM与径向PMF实操
对应TPR与窗口轨迹，查看重叠、去除谐偏置，并解释三维距离的径向几何修正。
阅读本篇目录
下载完整输入、实际拉伸与窗口轨迹、WHAM和绘图脚本。解压进入资料包目录，先处理已有结果：
bash analyze.sh
图片在 reference/figures，CSV和完整统计在 reference/analysis。想重新计算整套流程，在新解压目录执行：
bash submit.sh
脚本依次提交拉伸取结构、16个窗口数组任务、WHAM与绘图；每个窗口任务顺序跑两组种子。每个模拟任务申请1张P100、2个CPU核心，不要在登录节点运行模拟。最后的分析任务主要用CPU，附件模板仍按相应GPU分区分配资源与计费；已有结果可以直接提取，无需重新跑MD。
使用 gromacs/2026.4-p100-threadmpi，GNU13.3、CUDA12.8、thread-MPI/OpenMP、混合精度构建；数据分析使用ASE3.29.0环境中的Python、NumPy/SciPy和Matplotlib。生产阶段短程相互作用与PME使用GPU，坐标更新、约束和部分其他计算仍在CPU执行。这里的750 ps拉伸和32个窗口/种子条件组成一个多步工作流，不是33种独立性质。
1 WHAM输入为什么要一一对应#
每条窗口轨迹来自各自的偏置中心和力常数，WHAM需要同时读取TPR与pullx。顺序错了，即使程序正常结束，所去除的偏置也可能不对。
本资料包分析前，按每行实际时间保留100–1000 ps的4501条距离记录；原始 md_pullx.xvg 完整保留。操作只需：
module purge
module load gromacs/2026.4-p100-threadmpi
python3 -B run_wham.py reference
重跑数据将 reference 换成 work。脚本建立三组TPR/pullx列表：seed1、seed2、pooled，输出在 native-wham；列表行顺序来自同一份窗口清单。
2 原生WHAM命令#
在 reference/native-wham 中，单组命令可以写成：
gmx wham -it seed1-tpr.dat -ix seed1-pullx.dat \
-b 0 -min 0.2 -max 1.2 -bins 200 -temp 300 \
-tol 1e-8 -zprof0 0.9 -unit kJ \
-o manual.xvg -hist manual-hist.xvg
先在资料包根目录运行 bash analyze.sh，列表便会按当前解压位置重建。上述命令与资料包使用相同的0.2–1.2 nm范围、200个bin和接近0.9 nm的零点，便于不同种子和采样时间比较。筛选文件已经保留100–1000 ps，所以命令不再第二次截掉100 ps。
WHAM从TPR读取谐偏置，按各窗口的采样数重建未加偏置的分布。它不能自动补足没有访问的构象，也不能用程序收敛替代分子动力学采样收敛。
3 先看直方图，再读曲线#
左图是两次种子合并后的真实有偏分布。中心间距0.05 nm，输入k=1500 kJ mol⁻¹ nm⁻²，使相邻窗口仍有重叠。右图的诊断定义为两个归一化直方图逐格取较小值再求和；完全相同分布给1，本例最小约0.382。它不是通用MBAR重叠矩阵，不能与别的诊断阈值直接混用。
若相邻窗口之间几乎没有重叠，应增加中间窗口、调整合适的力常数，并检查起点和慢构象。曲线端点尤其容易受到少量样本支配。
4 距离坐标自由能与两粒子径向平均力势#
对本例三维距离坐标，直接重建的是：
F(r) = -RT ln P(r) + constant
球壳的可及体积随r²增大。若要去掉这部分几何熵，得到两粒子的径向PMF，可写为：
w(r) = -RT ln[P(r)/(4πr²)] + constant
= F(r) + 2RT ln r + constant
这里是加 2RT ln r。对Cartesian的z坐标等其他反应坐标，几何因子不同，不能机械套用这条公式。
图统一在r=0.8975 nm的bin中心设为0；PMF的加性常数是任意的。下表列实际合并数据的径向PMF和50 ps块重采样标准差：
| r / nm | w(r) / kJ mol⁻¹ | 块标准差 / kJ mol⁻¹ |
|---|---|---|
| 0.2825 | -7.3278 | 0.5438 |
| 0.3825 | 4.1186 | 0.3695 |
| 0.4825 | -2.4390 | 0.3137 |
| 0.6325 | -0.2776 | 0.2710 |
| 0.8325 | -0.0153 | 0.1534 |
| 0.8975 | 0.0000 | 0.0000 |
| 1.0325 | 0.1385 | 0.1675 |
这些是当前短轨迹和模型的结果；接触区仍随采样时间变化，不能把四位小数理解为真实物理精度。它也没有给出标准态缔合自由能：若要从径向PMF求缔合常数，还需定义结合区域、做带球壳体积权重的积分并处理标准态。
5 空bin与数据精度怎样处理#
本例200个bin中，合并数据有187个bin实际出现样本。原生输出在零直方图区间也可能显示普通数字，不能把这些数字解释成物理平台或势垒。配套CSV将未采样区标为无数据，图仅展示合并样本数不少于20的bin；完整原始直方图和原生曲线保留。
脚本按当前混合精度构建的分箱约定核对原生结果，并另存普通双精度分箱的对照。刚好位于bin边界的有限打印精度数据可落入相邻格；正式比较时要统一边界、bin宽和输出精度。
6 导出曲线与全部统计#
bash analyze.sh
head reference/analysis/PMF-profile.csv
head reference/analysis/umbrella-states.csv
CSV同时保留径向坐标F、修正后的w、样本数和时间块误差，不能只保留一条平滑的最终曲线。原生用法见 gmx wham官方说明；方法原文见 Hub等，2010。