GROMACS 约束与标准态修正
实际单原子谐约束轨迹，验证位移分布并讲清体积、单位与符号。
阅读本篇目录
下载完整输入、真实轨迹、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 去耦自由能为什么不是结合自由能#
水中小溶质去耦只比较溶质与溶剂相互作用开关。蛋白配体结合还涉及结合态和溶液态之间的转移、标准浓度、位置与取向约束，以及约束施加和释放。不能把 小溶质BAR结果直接改名为绝对结合自由能。
这里额外实际计算一个完全去耦的单原子，在盒内加各方向1000 kJ mol⁻¹ nm⁻²的谐位置约束。它没有分子取向自由度，适合单独解释平移体积修正；不是完整的蛋白配体双去耦算例。
2 约束输入和实际操作#
在AR的 [ atoms ] 之后、水分子定义之前加入：
[ position_restraints ]
1 1 1000 1000 1000
这是AR分子内部编号1号原子，三个数分别是x/y/z力常数。该对照同时将ARLJ的ε设为0、保持电荷0，完全关闭粒子与水的非键相互作用；因此不需要再启用λ变换。约束中心来自 start.gro 的实际参考坐标，不能随意假设总是盒中心。
完整工作流在所有λ任务结束后自动运行本例。在已分配资源中对应的原生命令为：
gmx grompp -f eq.mdp -c start.gro -r start.gro -p topol.top -o eq.tpr
gmx mdrun -deffnm eq -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu -pin off
gmx grompp -f md.mdp -c eq.gro -t eq.cpt -r start.gro -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu -pin off
gmx trajectory -s md.tpr -f md.xtc -select 'atomnr 1' -ox argon-position.xvg
-r start.gro 给位置约束参考坐标，生产阶段仍使用同一个中心。省掉或误改参考结构会改变所算问题。轨迹每0.1 ps输出，共10001帧，分析舍弃前100 ps。
3 用真实轨迹验证谐约束#
各方向理论位移方差为 RT/k，三维平均谐势能为 3RT/2。300 K和上述k下，理论方差0.00249434 nm²；真实轨迹三个方向分别为：
| 方向 | 实际方差 / nm² |
|---|---|
| x | 0.00230736 |
| y | 0.00251693 |
| z | 0.00267442 |
平均谐势能由实际坐标重建为3.749433 kJ/mol，理论为3.741508 kJ/mol。单帧坐标与原生约束能重建的差别受XTC默认0.001 nm精度影响，不能要求逐帧严格相等。
4 从约束体积换算标准态#
对于这里这个三维单原子谐约束，边界远离主要高斯分布区域，有效平移体积为：
Vrest = (2πRT/k)^(3/2)
ΔG_apply = -RT ln(Vrest / Vreference)
所有体积必须使用相同单位；这里是每个粒子的nm³，能量是每摩尔kJ。1 mol/L标准态每粒子体积为 10^24 / NA = 1.660539 nm³。
| 参考体积 | 体积 / nm³ | 施加该约束的ΔG / kJ mol⁻¹ |
|---|---|---|
| 模拟盒 | 17.576000 | 22.699271 |
| 1 mol/L标准态 | 1.660539 | 16.814148 |
| 本约束有效体积 | 0.001962018 | — |
施加约束使粒子的可及体积缩小，自由能升高；释放到1 M参考体积的修正符号相反，为−16.814148 kJ/mol。这个符号跟“从哪里到哪里”的方向一起使用，不能只记一个正负号套进不同的热力学循环。
5 哪些部分仍需为真实结合问题单独计算#
真实分子有转动和构象自由度，常用的距离/角度/二面角约束与这里的单原子三维位置约束不同。完整结合流程需匹配结合态和溶液态的相互作用路径，处理真实体系中的约束附加/释放、平移与取向标准态修正，必要时还要处理电荷有限尺寸与多结合姿态。这里的22.699或16.814都不能当作任意配体的通用修正值。
bash analyze.sh
python3 -m json.tool reference/analysis/TRAP-RESULTS.json | less
脚本用实际参考坐标读取约束中心，并从真实轨迹重建位移分布与约束能，再计算体积公式；短轨迹的50 ps块方差同图给出。正式应用仍需验证自己的约束定义、采样和完整热力学循环。
位置约束和自由能参数见 GROMACS官方MDP选项。