GROMACS 脂质膜建模与半各向异性平衡
用固定版本Insane和Martini3脂质参数构建DPPC、DOPC及混合双层膜，展示最小化、平衡和1ns生产的完整命令与MDP。
阅读本篇目录
本篇从建膜开始，完成最小化、NVT平衡、半各向异性NPT和生产，再进入膜面积、厚度与取向分析。使用Martini3粗粒化模型，适合先掌握建模与分析流程；不能把这些短轨迹直接当成平衡膜材料数据。
1 下载与查看#
下载DPPC、DOPC和混合膜输入、完整生产轨迹与分析脚本，约200 MiB，包含12个真实实例：
unzip gromacs-Martini3-bilayers-wave231*.zip -d membrane-example
cd membrane-example
sha256sum -c MANIFEST.sha256
bash view.sh
ls analysis/figures
view.sh只处理参考数据，不提交模拟。它使用集群现有的MDAnalysis2.9.0、NumPy和Matplotlib环境；该Python环境路径名称中带有qvasp，但这里不调用qvasp，也不运行VASP。其他机器可通过 ANALYSIS_PYTHON 指向具有相同依赖的Python。
2 本例到底有多少粒子#
每个膜有200个脂质，上下叶层各100个。分别计算纯DPPC、纯DOPC及1:1混合膜，两组建模/初速度种子各在P100和V100 16G运行。初始盒长8×8×10 nm，膜法向为z；生产目标323 K、1 bar。
| 建模种子 | 脂质数 | 粗粒化水珠数 | 总相互作用位点 |
|---|---|---|---|
| 13579 | 200 | 3032 | 5432 |
| 24680 | 200 | 3040 | 5440 |
随机建模会影响溶剂格点的填充数量，应该读取各自拓扑，不能假设所有种子的粒子数都相同。一个标准Martini粗粒化水珠代表4个水分子；表里的3032不是全原子水分子数量。
使用GROMACS2026.4混合精度GPU thread-MPI版本，以GCC/G++13.3、CUDA12.8编译，CPU指令为AVX2；固定Insane源码提交 e23a501036c6，Martini3相互作用参数与2025脂质v2参数。力场文件、来源、许可证和哈希均在附件中。不要把Martini2脂质拓扑与这里的Martini3相互作用参数混合。Insane与Martini3建膜说明
3 生成目录并提交#
python3 prepare.py
cd cases/p100-DPPC-s13579
cat config.json
cat job.slurm
sbatch job.slurm
squeue -u "$USER"
一次作业依次完成建膜和四个计算阶段，一张GPU、一个CPU核。模板对应GPU核时分区；套餐用户用 yskit 获取适合自己账户的分区，保留程序入口。
run_native.py会调用附件里的固定版本Insane，生成 membrane.gro 与原始分子数量，再按匹配的Martini3文件整理 topol.top，建立Bilayer和Solvent索引组。无需自行寻找不同年代的脂质参数。查看已完成实例的建模输入和结构：
cat BUILD.json
head -8 membrane.gro
head -15 topol.top
cat index.ndx | head -8
INSANE! Membrane UpperLeaflet>DPPC=1.0 LowerLeaflet>DPPC=1.0
5432
1DPPC NC3 1 2.971 1.991 6.950
1DPPC PO4 2 2.981 1.995 6.650
1DPPC GL1 3 3.013 1.978 6.350
1DPPC GL2 4 2.880 2.059 6.350
1DPPC C1A 5 2.983 2.009 6.050
1DPPC C2A 6 2.984 2.003 5.750
BUILD.json保存真正调用Insane的参数与种子。当前配置的核心等价于 insane -l DPPC:1 -x 8 -y 8 -z 10 -sol W；实际入口通过附件内的Insane源码调用，避免依赖另一个公共安装版本。混合膜给出DPPC和DOPC两个 -l 参数，各占一份。
想在运行前查看结构，可直接打开 reference-results 中对应实例的完整建模文件。新案例目录的结构由任务生成，避免先手动生成一遍、随后让任务重复覆盖。
4 四个阶段分别干什么#
| 阶段 | 目标 | 时间步 | 长度 |
|---|---|---|---|
| EM | 缓解初始接触 | 最陡下降 | 最多5000步，力阈值100 kJ/(mol·nm) |
| NVT | 初始化温度 | 5 fs | 20 ps |
| NPT平衡 | 开始调整膜平面面积与盒高 | 10 fs | 100 ps |
| NPT生产 | 保存轨迹、面积与厚度数据 | 20 fs | 1 ns |
粗粒化模型可以使用与全原子体系不同的步长，不能将20 fs照搬给普通全原子蛋白或脂质。初始构型若存在很大力，应先改善最小化或使用更短的平衡步长。
最小化输入#
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
生产输入#
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 = .02
nsteps = 50000
continuation = yes
gen-vel = no
semiisotropic让x和y一起缩放，z单独缩放，适合平面膜与水层分别响应压强。这里用C-rescale控压，tau-p=12 ps；两个压强目标分别对应xy平面和z方向。Bilayer、Solvent分别控温到323 K，避免只看总温度而忽略组间温差。
静电使用Reaction-Field、相对介电常数15，非键截断1.1 nm。脂质头基含正负电珠，但每个脂质及整个模型净电荷为零。这些是所用Martini3参数体系的一部分，不是通用全原子MD推荐值。
5 原生GROMACS命令顺序#
任务实际按以下顺序传递坐标和检查点；不要给 grompp 随意添加 -maxwarn 略过拓扑问题。
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
以上完整模拟应通过Slurm运行；轻量提取和绘图可在登录节点完成。生成的 *-processed.mdp、md-processed.top用于确认程序真正读到的参数与原子顺序，不只检查手写输入。
6 混合膜与下一步#
图中来自1 ns末帧，不是人为绘制的相分离示意。1 ns仍可能保留建模格点松弛和面积漂移，不能因图上出现某个局部聚集就判断发生相分离。继续膜分析教程，提取时间曲线后再决定如何延长平衡和采样。
附件脚本用于复现这组固定参数实例。修改体系时，要同步调整Insane建模参数、三份动力学MDP及分析定义，温度以实际MDP为准；不要仅修改JSON里的说明字段就认为所有物理输入都随之改变。