GROMACS 二硫键蛋白的完整模拟流程

本篇用1CRN的46残基crambin演示蛋白质的完整计算顺序。结构整理见上一节。二硫键分别是Cys3–Cys40、Cys4–Cys32、Cys16–Cys26；这里模拟的是含三条共价二硫键的氧化形式，不是可自动发生氧化还原反应的模型。

1 获取实例并准备新作业#

下载蛋白质、二硫键与血红素完整实例包，包含8组真实短模拟的输入、生产轨迹、TPR、EDR和绘图脚本。

unzip gromacs-protein-disulfide-heme-wave232*.zip -d protein-example
cd protein-example
sha256sum -c MANIFEST.sha256
bash view.sh
ls analysis/figures

view.sh只处理参考数据，不提交计算，也不需要下载力场。它使用集群现有的MDAnalysis2.9.0、NumPy、Matplotlib环境；可用 ANALYSIS_PYTHON 指定自己的Python。环境路径中带有qvasp名称，但这里不调用qvasp或VASP。8个实例是两种体系、两个种子、两类GPU的组合；硬件重复不是新增独立种子。

运行新模拟前加载环境，再生成案例目录：

module load ase/3.29.0 gromacs/2026.4-p100-threadmpi
python3 prepare.py
cd cases/p100-crambin-s13579
cat config.json
cat job.slurm
sbatch job.slurm
squeue -u "$USER"

prepare.py同时准备蛋白和血红素两套示例，会按血红素教程获取官方CHARMM36外部依赖；查看参考结果不需要这一步。每个Slurm任务一张GPU、一个CPU核，任务内部依次完成四个阶段。模板用GPU核时分区；套餐账户通过 yskit 选适合自己的分区，保留对应GROMACS程序入口。

使用GROMACS2026.4混合精度GPU thread-MPI构建，GCC/G++13.3、CUDA12.8、AVX2。AMBER99SB-ILDN配TIP3P，静电PME与1.0 nm截断；不要把CHARMM、GROMOS或其他模型的水和非键设置随意混合。

2 任务真正执行的建模命令#

以下命令在新案例目录中依次执行；自动任务已经执行过时无需再手动覆盖：

gmx pdb2gmx -f input.pdb -o processed.gro -p topol.top   -ff amber99sb-ildn -water tip3p -ignh -merge all
gmx editconf -f processed.gro -o box.gro -c -d 1.0 -bt cubic
gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top
gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o ions.gro -p topol.top   -pname NA -nname CL -neutral -conc 0.15 -seed 13579
# genion交互选择SOL

spc216.gro在这里提供初始水分子的空间坐标，真正的水相互作用由 topol.top 包含的TIP3P参数定义。genion替换水分子，添加中和离子及目标0.15 mol/L的盐；最终原子数以生成结构为准。

pdb2gmx应识别三对SG–SG连接，并将六个Cys按氧化形式处理；本例拓扑注释标为CYX，SG上没有巯基氢。不要一边保留SG–H，一边又手工加SG–SG键。

3 四个阶段如何衔接#

NVT和NPT使用 define=-DPOSRES，参考坐标固定为最小化的 em.gro；NPT用 refcoord-scaling=com配合盒缩放。生产输入去掉POSRES；氢键约束和水SETTLE仍保留，这与蛋白位置限制是不同概念。

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

完整模拟通过Slurm运行。上一阶段的 .cpt传递速度和耦合状态，-c传递坐标；不要在连续阶段重新生成速度。拓扑或缺失参数问题应先修正，不能用 -maxwarn 直接跳过。

EM完整输入#

cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
epsilon-r = 1
pme-order = 4
fourierspacing = .12
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
rcoulomb = 1.0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.0
DispCorr = EnerPres
integrator = steep
emtol = 500
emstep = .01
nsteps = 5000
nstenergy = 100
nstlog = 100

NVT完整输入#

cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
epsilon-r = 1
pme-order = 4
fourierspacing = .12
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
rcoulomb = 1.0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.0
DispCorr = EnerPres
integrator = md
nstenergy = 50
nstcalcenergy = 50
nstlog = 1000
nstxout-compressed = 500
compressed-x-precision = 1000
tcoupl = v-rescale
tc-grps = System
tau-t = .5
ref-t = 300
nsttcouple = 10
comm-mode = Linear
nstcomm = 100
dt = .001
nsteps = 10000
pcoupl = no
continuation = no
gen-vel = yes
gen-temp = 300
gen-seed = 13579
define = -DPOSRES

NPT平衡完整输入#

cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
epsilon-r = 1
pme-order = 4
fourierspacing = .12
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
rcoulomb = 1.0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.0
DispCorr = EnerPres
integrator = md
nstenergy = 50
nstcalcenergy = 50
nstlog = 1000
nstxout-compressed = 500
compressed-x-precision = 1000
tcoupl = v-rescale
tc-grps = System
tau-t = .5
ref-t = 300
nsttcouple = 10
comm-mode = Linear
nstcomm = 100
pcoupl = C-rescale
pcoupltype = isotropic
tau-p = 5
nstpcouple = 10
ref-p = 1
compressibility = 4.5e-5
dt = .002
nsteps = 10000
continuation = yes
gen-vel = no
define = -DPOSRES
refcoord-scaling = com

生产完整输入#

cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
epsilon-r = 1
pme-order = 4
fourierspacing = .12
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
rcoulomb = 1.0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.0
DispCorr = EnerPres
integrator = md
nstenergy = 50
nstcalcenergy = 50
nstlog = 1000
nstxout-compressed = 500
compressed-x-precision = 1000
tcoupl = v-rescale
tc-grps = System
tau-t = .5
ref-t = 300
nsttcouple = 10
comm-mode = Linear
nstcomm = 100
pcoupl = C-rescale
pcoupltype = isotropic
tau-p = 5
nstpcouple = 10
ref-p = 1
compressibility = 4.5e-5
dt = .002
nsteps = 25000
continuation = yes
gen-vel = no

4 提取RMSD、RMSF与二硫键长度#

任务完成后，在案例目录调用附件工具：

module load ase/3.29.0 gromacs/2026.4-p100-threadmpi
python3 ../../export_native.py
ls native-analysis
cat analysis-index.ndx

也可逐条使用原生命令，交互选择索引组 CA：

gmx rms -s md.tpr -f md.xtc -n analysis-index.ndx -o CA-rmsd.xvg -nomw
# 拟合组、计算组都选CA
gmx rmsf -s md.tpr -f md.xtc -n analysis-index.ndx -o CA-rmsf.xvg -res
# 选CA
gmx distance -s md.tpr -f md.xtc -select 'resnr 3 40 and name SG' -oall SS3-40.xvg
gmx energy -f md.edr -o temperature-density.xvg
# 选Temperature、Density，最后输入0

RMSD先去掉整体平移和旋转，再比较Cα与生产TPR中的参考坐标；本例参考是开始生产时的NPT末态，不是原始晶体结构。RMSF是同一拟合约定下各残基Cα围绕时间平均位置的均方根波动。周期边界下先保持分子完整；不要直接对被盒边切开的蛋白拟合。GROMACS rms说明、rmsf说明

5 八个实例实际得到什么#

生产每1 ps存一帧，0–50 ps共51帧；每0.1 ps记录一次能量，共501条。下面同时列出crambin和血红素模型，便于查对应目录；两种力场的能量不能直接用于材料或蛋白稳定性排名。

结构和数据图片来自附件中的实际轨迹。50 ps只用于展示输入传递与分析过程，不能据此判断蛋白长期稳定性、收敛构象分布或结合强度。短时曲线正常后，还需结合研究目标延长平衡和生产、安排独立重复。画新结果可在包根目录执行：

/public/toolkit/miniconda3/envs/qvasp-2.26/bin/python analyze.py   --cases cases --case p100-crambin-s13579 --output analysis-new
ls analysis-new/figures