GROMACS 离子对AWH与采样检查

下载完整输入、实际轨迹、原生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 用同一离子对比较两种采样方法#

本例使用Na⁺和Cl⁻各一个、576个SPC/E水分子，2.6 nm立方周期盒。离子来自本地OPLS-AA参数，PME静电、0.9 nm势平移LJ截断、不加色散尾修正；与伞形采样实例保持相同拓扑，便于理解方法之间的差别。

每个AWH条件先在约0.68 nm参考距离做20 ps固定伞形平衡，再切换原生AWH运行2 ns，300 K、2 fs、固定盒体积、随机动力学。各条件从同一份实际前序MD结构开始，两个种子分别改变随机序列；它们不是从两个完全独立的溶剂构象集合开始。

这里给出六个短算例的真实结果，主要用来学习配置、提取和检查采样。不要把本页数值当成已收敛的实验离子缔合自由能。

2 反应坐标改成三维最小镜像距离#

索引前两个单原子组为 NA（1）和 CL（2），另有包含全部1730个原子的System。拓扑为：

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"
#include "oplsaa.ff/ions.itp"
[ system ]
One NaCl pair in SPC/E water
[ molecules ]
NA 1
CL 1
SOL 576

重点是 geometry = distance、groups = 1 2、dim = Y Y Y。它采样两个单原子组的三维最小镜像距离r，不是z轴投影，也不是固定原点到一个粒子的坐标。

细网格条件的完整生产配置：

integrator = sd
dt = 0.002
tcoupl = no
tc-grps = System
tau-t = 1.0
ref-t = 300
pcoupl = no
nstcalcenergy = 100
nstenergy = 1000
nstlog = 10000
nstxout = 0
nstvout = 0
nstfout = 0
cutoff-scheme = Verlet
nstlist = 20
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
fourierspacing = 0.12
pme-order = 4
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 0.9
DispCorr = no
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-iter = 1
lincs-order = 4
continuation = yes
gen-vel = no
nsteps = 1000000
ld-seed = 2921851
nstxout-compressed = 500
pull = yes
pull-ngroups = 2
pull-ncoords = 1
pull-group1-name = NA
pull-group2-name = CL
pull-coord1-type = external-potential
pull-coord1-potential-provider = awh
pull-coord1-geometry = distance
pull-coord1-groups = 1 2
pull-coord1-dim = Y Y Y
pull-nstxout = 100
pull-nstfout = 100
pull-print-ref-value = no
awh = yes
awh-potential = convolved
awh-seed = 3011851
awh-nstout = 10000
awh-nstsample = 10
awh-nsamples-update = 10
awh1-error-init = 10
awh1-growth = exp-linear
awh1-ndim = 1
awh1-dim1-coord-index = 1
awh1-dim1-force-constant = 12000
awh1-dim1-start = 0.20
awh1-dim1-end = 1.10
awh1-dim1-diffusion = 0.003

awh1-dim1-diffusion = 0.003 是粗略的算法估计。改变离子、溶剂或反应坐标后应重新考虑，不要直接视为该离子的扩散测量结果。AWH范围应同时考虑体系有意义的坐标区间、盒子半长和坐标涨落；本例全部实际距离均小于1.3 nm最小镜像边界。

3 为什么做三档网格/范围对照#

离子接触势阱很窄，粗网格可能不能清楚表示其形状。提高kc会缩小辅助弹簧的热涨落尺度，原生程序相应生成更细网格；但它也会改变采样特征，不保证更快收敛。这里从1500到6000同时扩大了范围，因此比较的是“分辨率与范围敏感性”，不能把所有差异只归因于网格。

4 提交、提取和画图#

整套重新运行：

bash submit.sh

只重新运行第16个条件（初始NaCl第一组种子），在已分配的资源中执行：

bash run.sh 16

提取某个NaCl条件的最终结果：

cd reference/nacl-kc12000-s1
gmx awh -s md.tpr -f md.edr -b 2000 -e 2000 -o final.xvg -more
gmx trajectory -s md.tpr -f md.xtc   -select 'atomnr 1 or atomnr 2' -ox ion-coordinates.xvg

gmx trajectory 的输出能用于独立检查最小镜像距离。不要直接用跨周期边界的两组显示坐标相减后求长度，必须先按盒长折回最小镜像。附件分析脚本已经完成这一检查。

5 径向自由能和离子对PMF分开画#

原生AWH给的是所选标量r的自由能 F(r)=−RT ln P(r)+C。三维距离的几何体积元包含r²，若想去掉这部分几何贡献，画离子对势均力，应使用：

w(r) = F(r) + 2 RT ln(r / r_ref) + C

本例把r_ref=0.9 nm处归零。这里是加号；固定z投影或前面的笛卡尔谐势对照不能套用这一径向修正。改变自由能零点只是加常数，不改变势垒和势阱的形状。

上排为径向坐标F，下排为修正后的pair PMF。灰线是此前独立伞形/WHAM短样本，只有有实际支持的bin进入比较；它自身也有约1 kJ/mol的接触区时间变化，不是精确真值。所有曲线在0.9 nm使用同一插值归零，原始各自零点仍保留在CSV。

6 这次短算例的数据告诉我们什么#

0.28 nm处数值由各自实际网格线性插值得到，粗网格可能错过窄势阱的最低点。不同条件、不同种子仍存在数kJ/mol的差异。例如kc=12000时两组种子在2 ns得到约−5.135与−8.473 kJ/mol，不能取一条看起来平滑的线就宣称收敛。

初期阶段结束只说明AWH内部增长阶段转换，不代表所有慢自由度已经平衡，也不代表最终自由能达到某个精度。更细网格能提高曲线分辨率，不能替代足够长的采样。

例如第一组细网格在500、1000、1500、2000 ps的w(0.28 nm)约为−8.011、−4.562、−6.308、−5.135 kJ/mol。这里直接展示变化，不把算法报告的目标误差当作这些变化的独立误差条。

7 换成自己的体系，按什么顺序检查#

1. 先确认原子组、反应坐标方向和周期边界。对多个原子的组，要额外检查质心和周期跨越。
2. 用短任务确认选定范围能访问、没有结构或约束异常，再考虑加长。
3. 同时看原始坐标分布、参考位置分布、时间演变和不同随机种子；只看最终PMF外观不够。
4. 对较窄势阱检查kc和实际网格分辨率。改变范围或kc后，重新比较采样与曲线，不保证参数越大越好。
5. 用多个更长独立样本检查结果稳定性，再讨论物理结论。不要沿用本页2 ns作为所有体系的固定时长。

本例没有积分出标准态缔合常数、没有额外1 M修正，也没有计算真实动力学速率。若要做结合自由能，还需要明确结合区域、体积和标准态定义，以及其他约束的处理。更多输出细节见AWH输出与解析对照，窗口方法见WHAM与径向PMF。