GROMACS 原生AWH实操
从解析谐势模型到真实离子对，准备输入、平衡、提交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 先理解这套例子要做什么#
AWH（Accelerated Weight Histogram）沿选定坐标逐步调整偏置，让原本难以访问的区间获得更多采样，再估计原始体系的PMF。它需要你先定义坐标、范围、耦合弹簧与采样参数，并不是打开一个开关就能得到可靠自由能。
这套工作流先用一个有解析答案的模型认识输入和输出，再做真实溶剂中的NaCl离子对：
| 模型 | 坐标与参考答案 | 实际运行 |
|---|---|---|
| 577个SPC/E水分子中的一个解耦示踪粒子 | q=z−1.409 nm；物理势与原始PMF均为kp q²/2 | kp=250/500、kc=1500/3000、两种偏置模式、两个种子，共16组，各20 ps平衡+1 ns |
| 576个SPC/E水分子中的Na⁺和Cl⁻ | 三维最小镜像距离r；没有简单解析答案 | 三档kc/范围，每档两个种子，共6组，各20 ps平衡+2 ns |
这22个条件服务于同一个AWH教学工作流。它们不是22种独立材料性质，也不是收敛自由能数据库。
结构图来自实际输入坐标，水氧为投影。离子连线表示所选距离坐标，不表示添加化学键。
2 准备拓扑，分清两个弹簧#
先看 inputs/kp500-kc1500-convolved-s1。这里的示踪粒子使用Ar质量，但电荷和LJ相互作用都设为0，目的是隔离出一个能手算的自由能问题，不能把它的结果当成真实氩的水合性质。
完整拓扑：
#include "oplsaa.ff/forcefield.itp"
[ atomtypes ]
ARLJ 18 39.948 0.0 A 0.340 0.0
[ moleculetype ]
AR 1
[ atoms ]
1 ARLJ 1 AR AR 1 0.0 39.948
[ position_restraints ]
1 1 0 0 500
#include "oplsaa.ff/spce.itp"
[ system ]
One neutral Lennard-Jones Ar in SPC/E water
[ molecules ]
AR 1
SOL 577
[ position_restraints ] 中z方向的500是物理势的力常数kp，x、y为0。MDP中的 awh1-dim1-force-constant = 1500 是另一根AWH辅助弹簧的kc，两者用途不同，不能只看数值把它们当成同一参数。单位均为kJ mol⁻¹ nm⁻²。
索引里 AR 包含第1个原子，另有包含全部1732个原子的 System 组。完整索引已在附件中提供；自定义索引时要同时保留温控和其他配置引用的组。
3 生产阶段的完整MDP#
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
comm-mode = none
nstcomm = 0
pull = yes
pull-ngroups = 1
pull-ncoords = 1
pull-group1-name = AR
pull-coord1-type = external-potential
pull-coord1-potential-provider = awh
pull-coord1-geometry = direction
pull-coord1-groups = 0 1
pull-coord1-origin = 0 0 1.409
pull-coord1-vec = 0 0 1
pull-coord1-dim = N N Y
pull-nstxout = 100
pull-nstfout = 100
pull-print-ref-value = no
awh = yes
awh-potential = convolved
awh-seed = 3012601
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 = 1500
awh1-dim1-start = -0.3
awh1-dim1-end = 0.3
awh1-dim1-diffusion = 0.05
nsteps = 500000
ld-seed = 2922601
continuation = yes
gen-vel = no
nstxout-compressed = 500
关键项按下面理解：
| 配置 | 在这个算例中的作用 |
|---|---|
external-potential、potential-provider = awh |
将该pull坐标交给原生AWH提供偏置 |
geometry = direction、groups = 0 1 |
相对固定原点取方向投影，本例只取z方向 |
origin = 0 0 1.409、vec = 0 0 1 |
q=z−1.409 nm，不是两粒子间的三维距离 |
start = -0.3、end = 0.3 |
所采样的q区间；粒子涨落可到区间外，不是硬墙 |
awh-potential = convolved |
使用卷积偏置；其他对照使用 umbrella 模式 |
nstsample = 10、nsamples-update = 10 |
每10步取样，每10个样本更新；本例2 fs步长，对应0.02 ps和0.2 ps |
awh-nstout = 10000 |
每20 ps把AWH数据写入能量文件 |
diffusion = 0.05 |
用于算法的粗略坐标扩散估计，单位nm²/ps；不是本次测得的材料扩散系数 |
error-init = 10、growth = exp-linear |
初期适应参数；输出的目标误差不能代替独立统计误差 |
comm-mode = none |
解析对照中有外加绝对位置势，不再扣除整体质心运动 |
umbrella 模式通过辅助参考位置的采样施加谐偏置；convolved 模式把这些辅助位置的贡献卷积起来。两者面向相同的平衡自由能目标，实际短时轨迹与收敛速度可以不同。
4 先平衡，再启动AWH#
平衡用 eq.mdp，不启用AWH；生产接平衡后的坐标和速度。已分配计算资源中可以只跑第0个小对照：
bash run.sh 0
不要在同一目录重复执行去覆盖已有 work。要重算请重新解压到新目录。工作流内部实际使用以下顺序：
gmx grompp -f eq.mdp -c start.gro -r start.gro -p topol.top -n index.ndx -o eq.tpr
gmx mdrun -deffnm eq -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu
gmx grompp -f md.mdp -c eq.gro -r start.gro -t eq.cpt -p topol.top -n index.ndx -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu
这些命令在单个条件的工作目录中执行。-r start.gro 给物理位置约束提供参考，-t eq.cpt 传递平衡状态；生产设置 continuation=yes、gen-vel=no，不再重新生成速度。
5 得到哪些文件，接着看什么#
| 输出 | 用途 |
|---|---|
md.edr、md.tpr |
原生AWH后处理的主要输入；仅有XTC不足以重新提取AWH历史 |
md_pullx.xvg、md_pullf.xvg |
实际坐标和偏置力随时间变化 |
awh_t1000.xvg |
1 ns解析对照的最终AWH输出；NaCl对应 awh_t2000.xvg |
friction_t1000.xvg |
原生摩擦张量诊断 |
md.log |
时间步、初期阶段是否结束及运行过程 |
analysis/*.csv、figures |
已转换的七列结果和高清图 |
先读AWH输出与解析对照，理解每列含义，再读离子对AWH与采样检查。这里的随机动力学和自适应偏置用于平衡分布采样，不能用轨迹直接报告真实跃迁速率或动力学扩散系数。
原生算法和参数说明见 GROMACS AWH参考手册。