GROMACS BAR、TI与重叠分析
真实双种子数据，解释dhdl列、kT单位、BAR/TI、相邻重叠与连续块误差。
阅读本篇目录
下载完整输入、真实轨迹、dhdl、BAR结果与绘图脚本。解压进入资料包目录，先处理已有结果：
bash analyze.sh
reference/figures 保存数据图，reference/analysis 保存CSV、完整JSON与NPZ。想重新计算，在新解压目录只需：
bash submit.sh
脚本提交16个λ状态的数组任务，每个任务用1张P100、2个CPU核心顺序跑两组种子；全部结束后自动提交一个小约束对照与分析任务。不会让你手工串接一长串作业号。新结果在 work。这里是实际32个λ/种子条件加1个约束对照，不是33种独立性质。
本例使用 gromacs/2026.4-p100-threadmpi，GNU13.3、CUDA12.8、thread-MPI/OpenMP环境，生产阶段非扰动短程相互作用和PME使用GPU，扰动相关部分按原生程序的CPU路径执行，更新/约束使用CPU；不是所有计算都在GPU上。双精度或MPI/PLUMED构建与这里不同。数据处理使用ASE3.29.0环境中的Python、NumPy、SciPy与Matplotlib。不要在登录节点运行模拟；已有结果的提取与绘图无需重新跑MD。
1 先确定方向、时间和单位#
本资料包的λ从完整相互作用到完全去耦。每个条件生产1 ns，分析实际保留100–1000 ps的4501行dhdl。run_bar.py 先依据文件里每行的实际时间建立筛选副本，再调用原生BAR，原始 md.xvg 保留。
module purge
module load gromacs/2026.4-p100-threadmpi
python3 -B run_bar.py reference
脚本分别计算seed1、seed2和合并数据，结果在 reference/native-bar，选择记录在 SELECTION.json。你重跑的结果则将 reference 改为 work。
想手工处理一组，进入 reference/native-bar：
gmx bar -f selected/lambda*-seed1.xvg -temp 300 -prec 6 -o manual.xvg -oi manual-integral.xvg -oh manual-histogram.xvg
这些筛选文件已经包含正确时间区间，命令不再二次截时间。原生默认以5个时间块估计误差；输出文本最后一部分是kJ/mol。当前原生 -o XVG纵轴是kT，横轴是相邻状态序号的中点，不能直接当kJ/mol与真实λ。配套脚本根据图例确认单位，再乘RT并映射到实际λ数组。
2 BAR和TI的真实结果#
| 数据 | BAR去耦ΔG / kJ mol⁻¹ | TI梯形 / kJ mol⁻¹ | TI Simpson / kJ mol⁻¹ |
|---|---|---|---|
| seed1 | -9.407505 | -9.472582 | -9.413955 |
| seed2 | -9.505530 | -9.569208 | -9.418937 |
| pooled | -9.456421 | -9.520895 | -9.416446 |
BAR利用相邻两态的真实能差分布，TI积分每个λ的平均dH/dλ。配套程序分别实现BAR方程求根、非均匀λ梯形积分和Simpson积分，BAR逐段与原生结果对照。TI两种积分的差异保留，不能只挑最接近BAR的一个方法作为“正确答案”；应加密变化剧烈的λ区间，并检查采样。
本模型的合并BAR为−9.456421 kJ/mol，方向是去耦；反方向耦合为+9.456421 kJ/mol。这个数依赖当前固定密度、LJ截断与水模型，不是普适实验值。
3 重叠怎样看#
原生 -oh 写相邻态能差直方图。配套图将正向ΔU和反向数据取负后放在相同的 U(next)−U(current) 横轴上，不能把两种不同符号的工作直接叠在一起。
左图是本脚本明确定义的平均Fermi接受诊断：本例等样本数时，完全相同两态接近0.5；它不是所有软件通用的MBAR重叠矩阵，也不该套用别的重叠阈值。最小值约0.404，出现在λ=0.5→0.6。右图显示该对状态中央99.6%范围的实际直方图，原始尾部未从BAR计算中删除。
若相邻分布几乎不重叠，应优先增加中间λ状态；仅延长已有两端轨迹通常不能快速弥补很差的相空间覆盖。
4 两次种子、时间分块与误差#
两次独立重复相差约0.0980 kJ/mol。合并数据用50 ps连续块重采样100次，BAR标准差约0.1056 kJ/mol，TI约0.1130 kJ/mol；这些是当前短轨迹下的统计诊断，不是已经证明收敛的精确置信区间。
在每个状态内，同一时间块的全部外部能差列一起抽取，保持它们之间的相关性；不同状态、不同种子分别抽块。不能把每0.2 ps的一帧都当成独立样本，也不能把原生5块误差和自定义50 ps块误差误认为完全同一种估计。若慢过程超过块长，这两种误差都可能低估不确定性。
保留100 ps之后的200、400、650、900 ps时，合并BAR分别约−9.4973、−9.4241、−9.4747、−9.4564 kJ/mol。曲线看起来较平稳也不能单独证明盒大小、模型或λ离散误差已经收敛。
5 一次提取所有图和数据#
bash analyze.sh
head reference/analysis/free-energy-estimates.csv
head reference/analysis/state-sampling.csv
重跑结果用 bash analyze.sh work。完整逐段BAR、λ导数、重叠诊断与原始统计在 RESULTS.json；100次块结果在NPZ，时间窗口在 time-windows.json。数据不压成一张只有最终数字的截图，便于你按自己的体系修改时间和λ路径。
原生命令和统计解释参见 GROMACS gmx bar。本教程不计算电荷去耦、相变或配体结合；带电体系还需设计电荷路径与有限尺寸修正。