GROMACS 水模型与虚拟位点

准备水溶液时，先选与研究问题及溶质、离子参数相容的水模型，再生成坐标与拓扑。本篇直接从3nm立方水盒开始，完整走一遍建模、最小化、平衡、生产和提取流程。

1. 三个原子，为什么有时是四个位点#

TIP3P、SPC/E采用三个相互作用位点；TIP4P-Ew和本例OPC增加一个无质量的电荷位点M。M由O、H的几何位置构造，不是额外的化学原子，也不增加一个普通有质量粒子的运动自由度。

这些数值来自各模型实际生成的体系。四点水在GRO里多出虚拟点，因此不能拿三点水的坐标数直接对应四点模型的拓扑。

2. 先准备拓扑，再加水#

下面是本例TIP3P水盒的完整拓扑：

#include "amber99sb-ildn.ff/forcefield.itp"
#include "amber99sb-ildn.ff/tip3p.itp"
[ system ]
Native water model and ion-count example
[ molecules ]
SOL               884

建模时先留空 [ molecules ]，gmx solvate -p 会补上实际的SOL数量。上面展示的是生成后的真实文件。

module load gromacs/2026.4-p100-threadmpi
gmx solvate -cs spc216.gro -box 3 3 3 \
  -o solv.gro -p topol.top

三点水使用 spc216.gro 作为初始填充坐标，随后按选定模型的拓扑约束与最小化处理。四点水用 tip4p.gro，并在拓扑里引用对应的 tip4pew.itp 或 opc.itp。坐标模板名称不是最终力场名称，真正的电荷、LJ参数、约束和虚拟点由拓扑决定。

资料包已经提供四套独立流程。在分配好的资源里可以分别执行：

bash run.sh water-tip3p
bash run.sh water-spce
bash run.sh water-tip4pew
bash run.sh water-opc

TIP3P、SPC/E、TIP4P-Ew水盒采用本包的AMBER99SB-ILDN目录对应水文件；OPC采用AMBER14SB目录。AMBER14SB这套输入库注明需要GROMACS2026或更高，本篇固定使用2026.4，不能直接将新库套进旧版本程序。

3. 先最小化，消除不合适接触#

gmx grompp -f em.mdp -c ions.gro -p topol.top -o em.tpr
gmx mdrun -deffnm em -ntmpi 1 -ntomp 2 -nb cpu -pme cpu

纯水没有加离子时，脚本将溶剂化坐标复制为 ions.gro，使后续流程统一。这里要求最小化实际达到最大力小于1000kJ/(mol·nm)，再开始动力学。

integrator = steep
nsteps = 3000
emtol = 1000
emstep = 0.01
nstenergy = 100
nstlog = 100
cutoff-scheme = Verlet
nstlist = 20
verlet-buffer-tolerance = -1
rlist = 1.1
coulombtype = PME
rcoulomb = 0.9
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 0.9
DispCorr = no
fourierspacing = 0.12
pme-order = 4
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1

4. NVT、NPT，再进入生产#

gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt -ntmpi 1 -ntomp 2 -nb gpu -pme gpu
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt \
  -p topol.top -o npt.tpr
gmx mdrun -deffnm npt -ntmpi 1 -ntomp 2 -nb gpu -pme gpu
gmx grompp -f md.mdp -c npt.gro -t npt.cpt \
  -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 2 -nb gpu -pme gpu

上面突出主要文件传递，下载脚本还显式固定键合、更新和绑核选项。NVT生成初速度；之后使用CPT传递完整精度的坐标、速度和盒子，不重新随机生成速度。各阶段单独从0计时，因此生产文件的20ps不是整个流程的总时长。

integrator = md
dt = 0.002
nsteps = 2500
continuation = no
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = no
gen-vel = yes
gen-temp = 300
gen-seed = 2612300
comm-mode = Linear
nstcomm = 100
nstenergy = 100
nstcalcenergy = 100
nstlog = 1000
nstxout = 0
nstvout = 0
nstfout = 0
nstxout-compressed = 0
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 20
verlet-buffer-tolerance = -1
rlist = 1.1
coulombtype = PME
rcoulomb = 0.9
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 0.9
DispCorr = no
fourierspacing = 0.12
pme-order = 4
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
ld-seed = 2613300

integrator = md
dt = 0.002
nsteps = 5000
continuation = yes
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = C-rescale
pcoupltype = isotropic
tau-p = 5
ref-p = 1
compressibility = 4.5e-5
gen-vel = no
comm-mode = Linear
nstcomm = 100
nstenergy = 100
nstcalcenergy = 100
nstlog = 1000
nstxout = 0
nstvout = 0
nstfout = 0
nstxout-compressed = 0
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 20
verlet-buffer-tolerance = -1
rlist = 1.1
coulombtype = PME
rcoulomb = 0.9
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 0.9
DispCorr = no
fourierspacing = 0.12
pme-order = 4
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
ld-seed = 2613300

integrator = md
dt = 0.002
nsteps = 10000
continuation = yes
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = no
gen-vel = no
comm-mode = Linear
nstcomm = 100
nstenergy = 100
nstcalcenergy = 100
nstlog = 1000
nstxout = 100
nstvout = 100
nstfout = 0
nstxout-compressed = 100
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 20
verlet-buffer-tolerance = -1
rlist = 1.1
coulombtype = PME
rcoulomb = 0.9
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 0.9
DispCorr = no
fourierspacing = 0.12
pme-order = 4
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
ld-seed = 2613300

5. 最后拿到什么#

md.tpr是生产输入，md.cpt用于续算，md.xtc保存压缩坐标，md.trr保存坐标和速度，md.edr保存能量等原生量。

gmx check -f reference/water-opc/md.trr
printf 'Temperature\nPressure\nVolume\nDensity\n0\n' | \
  gmx energy -f reference/water-opc/npt.edr -o npt-state.xvg
bash analyze.sh

提取脚本还按拓扑中的实际质量和盒体积计算密度，单独检查SETTLE几何以及M位点构造。四点水仍只有三个有质量原子；刚性纯水的自由度在去除整体平动后为 6×水分子数−3，不能按GRO中的四个位点当成四个自由运动原子。

官方参考：力场与水模型。

下载、提交与结果位置#

下载全部输入、真实轨迹、拓扑、能量表和五张图的绘图脚本。解压进入资料包目录，先查看现成结果：

bash analyze.sh

图片在 reference/figures，总表在 reference/analysis/conditions.csv，各个实际条件的完整文件在 reference/条件名。

从头运行一个水盒示例：

sbatch run.slurm
squeue -u "$USER"

run.slurm 默认选择 water-tip3p。需要其他条件时，修改其中的 CASE_NAME 一行；所有条件名见 CASES.json。新结果写入 work/条件名。完整18条件对照可在新解压目录执行 bash submit.sh；只生成脚本先看、不提交，用 bash submit.sh --dry-run。已有工作目录不会被覆盖。

模板默认核时分区 gpu-8p100-hs，一张 P100、两个 CPU 核；套餐用户按自己可用的分区修改 #SBATCH --partition，例如 P100 套餐分区 gpu-8p100-tc。正文中的单条 bash run.sh ... 仅在已分配的任务资源里使用，不在登录/管理节点运行 MD。

统一环境：gromacs/2026.4-p100-threadmpi，GROMACS2026.4、GNU13.3、CUDA12.8、FFTW3.3.10、thread-MPI/OpenMP。建模、预处理和最小化使用 CPU；动力学阶段短程非键和 PME 使用 GPU，更新/约束使用 CPU。提取使用原生 GROMACS，数值与结构处理使用 NumPy、MDAnalysis2.9.0，图片由 Matplotlib 生成。

本组每个条件都从建模开始，完成最小化、5ps NVT、10ps NPT和20ps生产NVT。它们是短时操作实例；固定0.9nm截断、DispCorr=no和短平衡不用于给水模型、力场精度或性能排名。正式研究按所选参数集文献设置相容的截断/长程处理，并充分平衡与采样。