GROMACS 膜孔道建模与水占据数

本例学习怎样在膜上建立一个初始缺口，加入不重叠的水珠，再观察缺口周围的水和脂质怎样变化。采用 Martini3 粗粒化 DPPC 膜，不是蛋白质通道，也不把初始缺口当成稳定孔道。

下载完整输入、真实轨迹与绘图脚本。资料包约103 MiB，含8组生产轨迹及平衡阶段数据。解压进入资料包目录，先画已有数据：

bash view.sh

在 analysis/figures 得到4张高清图，各条件的时间序列在 analysis/条件名/pore-time-series.csv。画图只读取附件数据，不提交计算。

1 完整流程先后顺序#

1. Insane 建立 DPPC 双层膜，按网格移除孔口附近脂质。
2. gmx insert-molecules 尝试向孔口附近加入水珠，只保留实际放入且不重叠的珠子，更新拓扑中的分子数量与顺序。
3. 能量最小化，再做20 ps NVT、100 ps半各向异性NPT。
4. 采用平衡后盒子做1 ns固定盒子NVT，每5 ps保存一帧。
5. 用 gmx select 统计固定观察圆柱内的水珠，结合膜结构图、密度与尾链取向判断变化。

本例两种开孔参数1.6、2.4 nm，各用13579、24680两个建模/速度种子，分别在P100与V100 16G计算，共8组。两个平台重复同两个种子，不当成四个独立初始种子。

2 新建算例并提交#

python3 prepare.py
cd cases/p100-pore-r2p4-s13579
cat config.json
sbatch job.slurm
squeue -u "$USER"

V100 16G示例改为进入 cases/v100-pore-r2p4-s13579。每个任务使用1张对应GPU和1个CPU核心，模板为核时分区；套餐用户按自己的账户规则选择分区。已有结果目录不要覆盖，重新解压到新目录再准备。

GROMACS2026.4 thread-MPI，GNU13.3、CUDA12.8、FFTW3.3.10。短程非键计算使用GPU，键相互作用与坐标更新使用CPU。后处理环境含MDAnalysis2.9.0、NumPy和Matplotlib。

3 建模时哪些参数需要看懂#

build_membrane.py 调用附件内固定版本的 Insane，主要参数是：

-x 8 -y 8 -z 10 -sol W -hole 2.4 -l DPPC:1

单位为nm，初始盒子8×8×10 nm。-hole 在建模网格上选取要移除的脂质，不是最终孔径的测量值。本例网格孔中心在xy盒长的0.55、0.55处，不能直接假设它在0.5、0.5处。

小孔条件有198个DPPC，每侧99个；大孔条件182个，每侧91个。两种模型均保持电中性。本例没有离子、跨膜电压或膜蛋白。

Martini标准W珠代表四个水分子；图和表里的“W珠数”不直接写成原子级水分子数。使用Martini3基础参数与v2磷脂参数，不能混入Martini2拓扑。力场及Insane来源、许可证随附件提供。

Insane默认溶剂布置不会完整填充膜中间的缺口，所以脚本还执行：

gmx insert-molecules -f membrane-dry.gro -ci one-W.gro   -ip pore-water-positions.dat -nmol 40 -try 1 -dr 0 0 0   -rot none -seed 13579 -scale 1 -radius .235 -o membrane.gro

pore-water-positions.dat 给出40个候选点；实际加入多少，以 BUILD.json 和 topol.top 为准。小孔各加6个、大孔各加30个W珠。局部 vdwradii.dat 明确指定0.235 nm几何排斥半径，用于加水时检查重叠，不修改分子动力学力场参数。

cat BUILD.json
tail -n 12 topol.top

4 输入文件怎么写#

以下是附件中完整的各阶段MDP。脚本把速度种子替换后传给 grompp，处理后的实际输入也保留在每个结果目录中。

能量最小化：em.mdp#

cutoff-scheme = Verlet
nstlist = 20
coulombtype = Reaction-Field
rcoulomb = 1.1
epsilon-r = 15
epsilon-rf = 0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.1
rlist = 1.1
pbc = xyz
constraints = none
integrator = steep
emtol = 100
emstep = .01
nsteps = 5000
nstenergy = 100
nstlog = 100

20 ps NVT：nvt.mdp#

cutoff-scheme = Verlet
nstlist = 20
coulombtype = Reaction-Field
rcoulomb = 1.1
epsilon-r = 15
epsilon-rf = 0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.1
rlist = 1.1
pbc = xyz
constraints = none
integrator = md
nstenergy = 100
nstcalcenergy = 100
nstlog = 1000
nstxout-compressed = 100
compressed-x-precision = 1000
tcoupl = v-rescale
tc-grps = Bilayer Solvent
tau-t = 1 1
ref-t = 323 323
comm-mode = Linear
nstcomm = 100
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
dt = .005
nsteps = 4000
pcoupl = no
continuation = no
gen-vel = yes
gen-temp = 323
gen-seed = GENERATION_SEED

100 ps NPT：npt.mdp#

cutoff-scheme = Verlet
nstlist = 20
coulombtype = Reaction-Field
rcoulomb = 1.1
epsilon-r = 15
epsilon-rf = 0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.1
rlist = 1.1
pbc = xyz
constraints = none
integrator = md
nstenergy = 100
nstcalcenergy = 100
nstlog = 1000
nstxout-compressed = 100
compressed-x-precision = 1000
tcoupl = v-rescale
tc-grps = Bilayer Solvent
tau-t = 1 1
ref-t = 323 323
comm-mode = Linear
nstcomm = 100
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
pcoupl = C-rescale
pcoupltype = semiisotropic
tau-p = 12
ref-p = 1 1
compressibility = 3e-4 3e-4
dt = .01
nsteps = 10000
continuation = yes
gen-vel = no

1 ns固定盒子NVT：md.mdp#

cutoff-scheme = Verlet
nstlist = 20
coulombtype = Reaction-Field
rcoulomb = 1.1
epsilon-r = 15
epsilon-rf = 0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.1
rlist = 1.1
pbc = xyz
constraints = none
integrator = md
nstenergy = 250
nstcalcenergy = 250
nstlog = 1000
nstxout-compressed = 250
compressed-x-precision = 1000
tcoupl = v-rescale
tc-grps = Bilayer Solvent
tau-t = 1 1
ref-t = 323 323
comm-mode = Linear
nstcomm = 100
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
dt = .02
nsteps = 50000
continuation = yes
gen-vel = no
pcoupl = no

NPT的 pcoupltype = semiisotropic 让xy方向一起缩放，z方向独立缩放。生产阶段改为 pcoupl = no，让观察圆柱、二维网格和三维体素都对应同一个固定盒子。1 ns是本例的短操作演示，不表示孔道寿命或膜性质已经收敛。

在已分配的GPU资源内，脚本按以下顺序执行；登录节点使用 sbatch job.slurm 提交即可：

gmx grompp -f em.mdp -c membrane.gro -p topol.top -n index.ndx -o em.tpr
gmx mdrun -deffnm em -ntmpi 1 -ntomp 1 -nb gpu -bonded cpu -update cpu
gmx grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -deffnm nvt -ntmpi 1 -ntomp 1 -nb gpu -bonded cpu -update cpu
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr
gmx mdrun -deffnm npt -ntmpi 1 -ntomp 1 -nb gpu -bonded cpu -update cpu
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 1 -nb gpu -bonded cpu -update cpu

5 计算完以后提取什么#

在新计算的条件目录执行轻量后处理：

module load gromacs/2026.4-p100-threadmpi
python3 ../../export_native.py
cat cylinder-definition.json

V100选择对应模块。程序从实际TPR读取完整精度盒长，观察圆柱半径0.5 nm、总高1 nm，位于建模孔轴与生产盒子中平面交点。圆柱固定不移动，统计的是这个空间区域的W珠数。它不等于整个孔中总水量，也不能凭此计算通透率。

回到资料包根目录，分析刚完成的条件：

cd ../..
/public/toolkit/miniconda3/envs/qvasp-2.26/bin/python analyze.py   --cases cases --case p100-pore-r2p4-s13579 --output analysis-new
cat analysis-new/p100-pore-r2p4-s13579/metrics.json

若使用自己的Python环境，需要MDAnalysis、NumPy、Matplotlib。不要将不同条件的轨迹或能量文件混在同一目录。

6 八组实际结果怎么看#

灰色为脂质珠，蓝色为水珠，橙圈标出固定观察圆柱的xy投影。上排仅显示膜中平面薄层中的水珠，下排显示固定宽度侧面切片。绘图使用真实保存坐标。

本例大孔生产起点圆柱内还有2–6个W珠，末帧均为0；最后500 ps多数条件持续为0，一组均值0.020。小孔在生产起点已干燥，变化可能发生在前面的平衡阶段，所以附件保留NVT/NPT轨迹供继续观察。这里展示的是无额外稳定约束的短时缺口弛豫，不能称为稳定开放通道，也不能把某一帧有水当成跨膜输运证据。

要判断真实材料中的稳定孔、膜蛋白通道或通透性，需要合适的全原子/粗粒化模型、充分平衡、独立重复、自由能与传输事件分析。改变体系后重新确定几何观察区域，不能直接照抄本例坐标。

继续查看二维密度与三维空间分布；普通完整膜的建模、面积和厚度见GROMACS粗粒化膜建模。