Packmol接GROMACS：水盒、拓扑与短MD

下载完整建模、拓扑、运行和绘图脚本。解压后进入water-md子目录。

1 坐标之外，还要准备拓扑#

Packmol负责摆放分子，不能替你生成任意分子的力场参数。本例完整搭建216个SPC/E水分子、648原子的周期盒，先装箱，再给GROMACS配套拓扑，最小化后运行20 ps NVT。所用Packmol 21.2.3为GFortran 13.3串行构建，GROMACS 2026.4使用8173m CPU thread-MPI模块；输入准备和画图用现有GPUMDkit Python环境。

水分子模板采用SPC/E几何：O–H为1.0 Å，H–O–H为109.47°。盒边长18.62578655 Å按216个水和初始1 g/cm³密度计算。目标密度是建模输入，未经过NPT求得。

2 生成初始盒与GRO坐标#

module purge
module load gpumdkit/1.5.8 packmol/21.2.3
python prepare.py
packmol < packmol.inp > packmol.log
python convert.py
cat PACKING.json

完整Packmol输入：

tolerance 2.0
filetype xyz
output packed.xyz
seed 27129
nloop 500
precision 0.000001
pbc 0 0 0 18.6257865512 18.6257865512 18.6257865512
structure water.xyz
 number 216
 inside box 0 0 0 18.6257865512 18.6257865512 18.6257865512
end structure

pbc检查周期盒两侧分子之间的距离。本次最短不同分子原子距离为2.198108 Å，大于2.0 Å下限。GRO写入时将Å换成nm，并按每个分子OW HW1 HW2排序、残基名写为SOL；这些名称和顺序必须与拓扑对应。

3 写配套拓扑#

topol.top完整内容：

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"

[ system ]
Packmol SPC/E water

[ molecules ]
SOL 216

GROMACS自带oplsaa.ff/spce.itp定义SOL的电荷和SETTLE刚性约束，O电荷−0.8476e、H各+0.4238e。SOL 216必须与坐标中实际水分子数量一致。这里是为水选取已知模型，不是自动为任意有机分子赋参。

4 最小化后再启动动力学#

em.mdp：

integrator=steep
nsteps=2000
emtol=1000
emstep=0.01
cutoff-scheme=Verlet
coulombtype=PME
rcoulomb=0.8
rvdw=0.8
pbc=xyz
constraints=none
nstenergy=1

module purge
module load gromacs/2026.4-8173m-threadmpi
gmx grompp -f em.mdp -c packed.gro -p topol.top -o em.tpr
gmx mdrun -deffnm em -ntmpi 1 -ntomp 2 -nb cpu

本次13步达到Fmax<1000 kJ mol⁻¹ nm⁻¹，最大力831.563 kJ mol⁻¹ nm⁻¹，势能−5829.6245 kJ/mol。这是初始去重叠标准，不是量子化学结构优化精度。

nvt20.mdp完整内容：

integrator=md
dt=0.001
nsteps=20000
continuation=no
cutoff-scheme=Verlet
coulombtype=PME
rcoulomb=0.8
rvdw=0.8
pbc=xyz
tcoupl=v-rescale
tc-grps=System
tau-t=0.1
ref-t=300
pcoupl=no
gen-vel=yes
gen-temp=300
gen-seed=27129
constraints=none
nstenergy=10
nstxout-compressed=100

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

水的刚性约束由拓扑中的SETTLE提供，constraints=none不表示这些水变成了柔性模型。时间步长0.001 ps，20000步共20 ps，v-rescale温控300 K；没有压力耦合。

5 提取实际曲线#

printf 'Potential\nTemperature\nPressure\n0\n' | gmx energy -f nvt20.edr -o thermo20.xvg

最初0.2 ps的试运行温度平均约451 K，说明随机装箱后的结构还在释放能量。延长后的10–20 ps区间温度为 298.51 K，时间序列标准差12.15 K；压力平均 309.66 bar、标准差962.25 bar。这里标准差描述瞬时波动，不是统计误差条。

温度接近300 K不等于压力已达1 bar。NVT固定体积，短小盒的压力波动也很大；需要目标压力与平衡密度时，应继续设置合适的NPT平衡，并检查尺寸和采样收敛。本例不用于报告扩散系数或水的热力学基准。

6 换成自己的体系时改什么#

同时修改分子模板、数量、盒子和拓扑。混合溶剂还要为每个分子选择兼容参数，并保持坐标顺序与[ molecules ]一致；只改Packmol数量、不改拓扑，不能直接运行。继续采样前检查能量、温度、压力、约束警告和结构是否合理。

附件run.sh还包含下一篇LAMMPS路线，会在CPU运行两套短工作流；只需要GROMACS时可按本篇命令逐步运行。参考：Packmol输入说明、GROMACS力场组织。