GROMACS 离子对伞形采样实操

下载完整输入、实际拉伸与窗口轨迹、WHAM和绘图脚本。解压进入资料包目录，先处理已有结果：

bash analyze.sh

图片在 reference/figures，CSV和完整统计在 reference/analysis。想重新计算整套流程，在新解压目录执行：

bash submit.sh

脚本依次提交拉伸取结构、16个窗口数组任务、WHAM与绘图；每个窗口任务顺序跑两组种子。每个模拟任务申请1张P100、2个CPU核心，不要在登录节点运行模拟。最后的分析任务主要用CPU，附件模板仍按相应GPU分区分配资源与计费；已有结果可以直接提取，无需重新跑MD。

使用 gromacs/2026.4-p100-threadmpi，GNU13.3、CUDA12.8、thread-MPI/OpenMP、混合精度构建；数据分析使用ASE3.29.0环境中的Python、NumPy/SciPy和Matplotlib。生产阶段短程相互作用与PME使用GPU，坐标更新、约束和部分其他计算仍在CPU执行。这里的750 ps拉伸和32个窗口/种子条件组成一个多步工作流，不是33种独立性质。

1 本例的体系和反应坐标#

在2.6 nm立方周期盒中放Na⁺、Cl⁻各一个和576个SPC/E水分子。离子参数来自本地OPLS-AA的 ions.itp，采用0.9 nm势平移LJ截断、PME静电、不加色散尾修正。这是一组用于学习流程的模型组合，不作为实验离子缔合数据。

坐标r是两个单原子组之间的三维最小镜像距离，不是只取z方向。盒子半长1.3 nm，中心范围0.28–1.03 nm；设计时同时为热涨落保留距离边界余量。各窗口保持相同水数、拓扑和盒体积，不能为了避免某个起点重叠而随意删掉不同数量的水。

图中离子连线表示距离坐标，不是力场中新增的一条化学键；为展示周期结构，对整盒作了平移并按周期等价位置显示。

2 先建立水盒和索引#

拓扑主体：

#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

初始模板只列NA和CL，gmx solvate -p 自动补水数；上面展示补齐后的版本，不要重复添加SOL。索引前两组如下：

[ NA ]
1
[ CL ]
2

实际文件还包含全部1730个原子的 System 组，用于温控；脚本按真实溶剂化后的原子数生成，每行15个编号。提供自定义索引时，要把所用的温控/输出组也包含进去。

3 先拉伸，再从真实轨迹取窗口起点#

拉伸坐标核心设置如下，其他完整配置在 steer-input/steer.mdp：

pull = yes
pull-ngroups = 2
pull-ncoords = 1
pull-group1-name = NA
pull-group2-name = CL
pull-coord1-type = umbrella
pull-coord1-geometry = distance
pull-coord1-groups = 1 2
pull-coord1-dim = Y Y Y
pull-coord1-start = no
pull-coord1-init = 0.28
pull-coord1-rate = 0.001
pull-coord1-k = 1500
pull-nstxout = 100
pull-nstfout = 100
pull-print-ref-value = yes

rate 单位nm/ps；参考距离从0.28 nm在750 ps内移动到1.03 nm，实际距离会在参考值附近涨落。单条非平衡拉伸轨迹用于取结构，不能直接把沿途功积分当作平衡PMF。

在已经分配的计算资源中，准备阶段可单独执行：

bash prepare.sh

它先能量最小化、50 ps固定参考平衡，再做实际750 ps拉伸，用原生 gmx trjconv 每50 ps取一个起点。例如工作目录中的实际命令：

gmx trjconv -s steer.tpr -f steer.xtc -dump 350 -o window07.gro

选择 System 输出全体系。提取后的16帧保持相同原子顺序，程序据此建立 runtime-inputs；原来的 inputs 作为参考模板保留。不要把拉伸过程中显示的参考距离当作每帧实际离子间距。

4 固定窗口里的完整生产配置#

第0窗口第一组种子的实际生产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
pull = yes
pull-ngroups = 2
pull-ncoords = 1
pull-group1-name = NA
pull-group2-name = CL
pull-coord1-type = umbrella
pull-coord1-geometry = distance
pull-coord1-groups = 1 2
pull-coord1-dim = Y Y Y
pull-coord1-start = no
pull-coord1-init = 0.28
pull-coord1-rate = 0
pull-coord1-k = 1500
pull-nstxout = 100
pull-nstfout = 100
pull-print-ref-value = yes
nsteps = 500000
ld-seed = 2921800
continuation = yes
gen-vel = no
nstxout-compressed = 2500

每个窗口把 pull-coord1-init 改成中心值，并将 rate 固定为0；16个中心每0.05 nm递增。两个种子分别改变初速度和随机动力学序列。每窗20 ps平衡后接续1 ns生产，300 K、2 fs、固定盒体积。

integrator = sd 用于这里的平衡分布采样，轨迹不能直接解释为真实动力学输运过程。本例平均温度299.09–301.16 K，各条件平均压力约−111到12 bar；没有做1 bar恒压控制。

5 实际窗口命令与输出#

数组任务内部依次执行：

gmx grompp -f eq.mdp -c 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 -pin off
gmx grompp -f md.mdp -c eq.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 -pin off

md_pullx.xvg 给时间、实际距离和参考距离，md_pullf.xvg 给偏置力；md.tpr保存中心和力常数等运行参数，WHAM要使用与每条pullx对应的TPR。md.xtc保存真实结构轨迹。

下一步看 WHAM、窗口重叠与径向PMF；短采样的重复和时间误差见 PMF采样与偏置检查。拉伸参数定义见 GROMACS官方MDP选项。