GROMACS PMF采样与偏置检查
实际验证偏置力和周期距离,比较种子、采样时间与连续块统计,识别尚未稳定的接触区。
阅读本篇目录
下载完整输入、实际拉伸与窗口轨迹、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 先分清施加的偏置与要计算的PMF#
每个窗口施加已知谐偏置:
U_bias(r) = k(r-r0)²/2
F_bias(r) = -k(r-r0)
偏置使轨迹在r0附近采样;WHAM要去除它,得到原来体系的分布。因此,画出一条输入抛物线并不等于得到了目标PMF。

图使用中心0.73 nm窗口的真实坐标、拉力和原生 COM Pull En.。距离、力与能量的XVG打印精度有限,逐点重建应考虑它和混合精度算术,而不是要求每位完全相等。
2 用终端分别提取距离和偏置能#
在已完成窗口目录里:
gmx trajectory -s md.tpr -f md.xtc -select 'atomnr 1 or atomnr 2' -ox ion-coordinates.xvg
gmx energy -f md.edr -o bias-energy.xvg
能量菜单选择 COM-Pull-En.、Temperature、Pressure,最后输入0。程序按完整名称选择偏置能,不用固定菜单数字;读取时核对XVG图例。
周期盒内的位置不能直接相减当作距离。先按盒长对离子差向量取最小镜像,再求长度;配套脚本将实际XTC离子坐标与原生pullx对照。分析XTC还要考虑默认0.001 nm的坐标精度。
3 两次种子与时间窗口都要看#
本例每个窗口做两次独立初速度/随机序列,另外比较保留100 ps之后的200、400和900 ps。它们使用相同的bin与PMF零点。

以r=0.2825 nm的接触区为例,所保留的200、400、900 ps分别得到:
| 每窗/种子保留时间 / ps | w(0.2825 nm) / kJ mol⁻¹ |
|---|---|
| 200 | -8.5565 |
| 400 | -8.1027 |
| 900 | -7.3278 |
接触区在这个短测试中仍有约1 kJ/mol的时间变化;一次曲线看起来平滑,不表示采样已经足够。面向自己的体系,应延长关键窗口,增加独立重复,并检查离子水合或分子构象等慢变量。
4 时间块误差与整条直方图误差#
资料包提供两种不同的诊断,各重采样50次:
- 原生WHAM对整条窗口直方图赋随机权重的Bayesian bootstrap。
- 配套程序在每个种子、每个窗口内,按50 ps连续时间块重采样,然后重新做WHAM。
每条轨迹保留900 ps,共18块;没有把每0.2 ps的一帧都视为独立样本。两种方法的误差不必相同,也不应相加成一个“更准确”的误差条。若慢过程超过块长或两次种子都没有访问另一个构象,它们仍可能低估真实不确定性。
图中的阴影是相同PMF零点下的块标准差,不包括力场模型、盒大小、反应坐标或未采样构象的系统误差。
5 为什么单次拉伸的功不当作自由能#
750 ps非平衡拉伸的实际协议功积分约8.5163 kJ/mol。它包括具体速率和路径下的非平衡过程;这里用这条轨迹取窗口起点,再通过平衡分布重建PMF。不能将这个功值直接改名为离子对分离自由能。
如果关心真实解离动力学,还需使用适当的动力学方案;本例随机动力学热浴轨迹用于分布采样,不直接提供实验离子对寿命。
6 数据保留到哪一层#
bash analyze.sh
python3 -m json.tool reference/analysis/RESULTS.json | less
head reference/analysis/PMF-profile.csv
重新计算的结果用 bash analyze.sh work。JSON保存每个窗口的温度、平均距离、重叠与时间窗口,NPZ保留完整WHAM、径向修正和每次块重采样曲线。原始pullx/pullf、XTC、TPR和能量随包提供,可以更改bin、时间区间和块长继续分析。
统计方法和限制参见 GROMACS WHAM官方说明。这批输入提供普通用户可复用的多步流程;正式缔合自由能、精确PMF或速率仍需按目标体系完善模型与收敛检查。