GROMACS AWH输出与解析对照
原生提取、kT与摩擦单位、七列输出，以及16组谐势实算和解析自由能比较。
阅读本篇目录
下载完整输入、实际轨迹、原生AWH输出、CSV与绘图脚本。解压进入资料包目录，先处理已经算好的结果：
bash analyze.sh
图片写入 reference/figures，逐条件CSV、JSON和NPZ写入 reference/analysis。想重新计算，在新解压目录执行：
bash submit.sh
脚本提交22个条件的数组任务，全部完成后自动提取和画图。每个模拟任务申请1张P100、2个CPU核心；分析任务主要用CPU，附件模板仍申请GPU分区资源，按相应分区计费。不要在登录节点运行MD。新结果写入 work，已有 reference 保留。
使用 gromacs/2026.4-p100-threadmpi：GNU13.3、CUDA12.8、thread-MPI/OpenMP、混合精度；短程非键相互作用和PME走GPU，更新、约束等仍在CPU执行。分析使用ASE3.29.0环境中的Python、NumPy、SciPy和Matplotlib。这里调用GROMACS原生AWH，不需要加载PLUMED；MPI复制交换使用另一构建，见温度复制交换与轨迹拆分。
1 直接调用原生工具提取#
在已有结果的单个条件目录执行：
module load gromacs/2026.4-p100-threadmpi
cd reference/kp500-kc1500-convolved-s1
gmx awh -s md.tpr -f md.edr -o awh.xvg -more -fric friction.xvg
输出按时间分成 awh_t0.xvg、awh_t20.xvg，一直到 awh_t1000.xvg；对应摩擦数据为 friction_t时间.xvg。这不是把所有时间叠在一个文件里。
只提最后一帧，并把能量单位改成kT：
gmx awh -s md.tpr -f md.edr -b 1000 -e 1000 -o awh-kT.xvg -more -kt -fric friction-kT.xvg
NaCl生产到2000 ps，需把这里的1000改成2000。300 K时RT=2.494338785 kJ/mol；-kt 将PMF和Coord bias除以RT，也改变摩擦张量相应能量单位，不会把采样分布列一起除以RT。
想少输出一些时间点：
gmx awh -s md.tpr -f md.edr -b 200 -e 1000 -skip 4 -more -o awh-sparse.xvg
本例原始AWH帧间隔20 ps，这条命令实际得到200、280、360……1000 ps，每4个可用帧保留1个。-skip 4 不表示只处理4帧，也不改变原MD采样。
2 七列输出分别是什么#
| 列 | 原生名称 | 读法 |
|---|---|---|
| 1 | coordinate | 本例单位nm；是网格坐标，不是时间 |
| 2 | PMF | 原始反应坐标的自由能估计，默认kJ/mol；每帧减去了自身最小值 |
| 3 | Coord bias | 经过归零的卷积偏置表示；在convolved模式中与实际施加势符号相反，不能直接当成原始PMF |
| 4 | Coord distr | 原始坐标的相对采样权重 |
| 5 | Ref value distr | 辅助参考位置的相对采样权重 |
| 6 | Target ref value distr | 辅助参考位置的目标相对权重；本例均匀目标为1 |
| 7 | sqrt(friction metric) | 摩擦张量平方根按网格平均值归一化后的诊断量 |
第4–6列在原生一维均匀网格输出中平均值为1，总和等于网格点数；它们不是已经按坐标积分归一的概率密度。自行绘制概率密度时应明确重新归一化，不要把“1”理解成每个网格100%概率。
第3列在 umbrella 模式下仍表示卷积偏置，不能用它替代某一时刻实际辅助位置的谐势。不同时间的归零常数可变，比较两帧必须先使用相同参考位置。
摩擦文件默认的一维单位为kJ mol⁻¹ nm⁻² ps；配合 -kt 后为nm⁻² ps。它描述该受偏置采样的诊断特征，并不是水的宏观黏度或已经测得的无偏扩散系数。
3 PMF和辅助自由能不是同一个答案#
本例物理势只约束示踪粒子的z方向，其余相互作用解耦，因此原始坐标q的PMF有解析式：
Phi(q) = kp * q² / 2 + C
加入AWH辅助弹簧kc后，对原始坐标积分所得到的辅助自由能为：
F(lambda) = [kp * kc / (kp + kc)] * lambda² / 2 + C'
例如kp=500、kc=1500时，原始PMF曲率是500，辅助自由能曲率是375。原生 gmx awh 的PMF列估计前者，不能拿375作为它必须收敛到的目标。附件还用独立的高精度数值积分核对了这个卷积解析式。
图中黑线是原始PMF，灰虚线是辅助自由能；每条实算曲线仅在中心区拟合一个任意常数，不改变其形状或曲率。端点附近受有限范围和采样影响更明显。
4 实际16组结果#
对每个1 ns结果，在 |q|≤0.24 nm拟合二次曲率，并与物理解析PMF比较。下表RMSE只减去一个最佳常数，单位kJ/mol；它是这次模型对照的偏差，不是一般体系的统计置信区间。
| kp | kc | 模式 | 种子 | 实算拟合曲率 | 解析对照RMSE |
|---|---|---|---|---|---|
| 250 | 1500 | convolved | 1 | 247.01 | 0.213 |
| 250 | 1500 | convolved | 2 | 248.57 | 0.104 |
| 250 | 1500 | umbrella | 1 | 268.83 | 0.220 |
| 250 | 1500 | umbrella | 2 | 237.20 | 0.174 |
| 250 | 3000 | convolved | 1 | 238.74 | 0.234 |
| 250 | 3000 | convolved | 2 | 240.25 | 0.228 |
| 250 | 3000 | umbrella | 1 | 253.12 | 0.137 |
| 250 | 3000 | umbrella | 2 | 259.39 | 0.160 |
| 500 | 1500 | convolved | 1 | 482.16 | 0.304 |
| 500 | 1500 | convolved | 2 | 494.39 | 0.126 |
| 500 | 1500 | umbrella | 1 | 480.55 | 0.202 |
| 500 | 1500 | umbrella | 2 | 489.28 | 0.121 |
| 500 | 3000 | convolved | 1 | 509.13 | 0.146 |
| 500 | 3000 | convolved | 2 | 501.47 | 0.124 |
| 500 | 3000 | umbrella | 1 | 504.91 | 0.119 |
| 500 | 3000 | umbrella | 2 | 512.49 | 0.212 |
两个随机种子和不同模式仍有差别。不要因为某个单一条件更接近解析式，就声称它对所有体系更优；也不要用这一简单模型的误差替任意蛋白或离子对提供精度保证。
5 从坐标独立检查物理势能#
在同一条件目录提取原生坐标和能量：
gmx trajectory -s md.tpr -f md.xtc -select 'atomnr 1' -ox argon-position.xvg
printf 'Position-Rest.
Temperature
0
' | gmx energy -f md.edr -o energy.xvg
由XTC中z坐标计算q=z−1.409 nm，再由高频 md_pullx.xvg 计算 kp q²/2，与原生 Position Rest. 对照。图中残差主要处于混合精度和XVG打印精度量级。XTC坐标精度为0.001 nm，不能要求从它重建距离时得到任意多位小数。
批量CSV提取和上述图已包含在 bash analyze.sh 中。七列原始数值逐条件保存在CSV，不用从图片读数。对自己的体系，接着检查多次独立重复、时间演变、采样范围和隐藏的慢自由度，见NaCl实例。
输出选项见 gmx awh官方帮助。