GROMACS 血红素辅因子与铁配位建模

含辅因子蛋白不能只生成标准氨基酸拓扑后直接运行。本例使用1A6N去氧肌红蛋白构建体，保留血红素并建立His93 NE2–Fe连接；用官方预定义CHARMM36参数和原始PHEM补丁表达经典配位模型。该模型不计算电子转移、氧化态变化、铁自旋或配体成键反应。

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的组合；硬件重复不是新增独立种子。

全新模拟需要官方 charmm36-feb2026_cgenff-5.0.ff.tgz。附件提供下载与63文件哈希检查脚本，没有重复打包整套力场。输入参考结果的分析不依赖这次下载。

module load ase/3.29.0 gromacs/2026.4-v100-threadmpi
python3 get_forcefield.py
ls forcefields/charmm36-feb2026_cgenff-5.0.ff
python3 prepare.py
cd cases/v100-myoglobin-s13579
sbatch job.slurm

登录环境无法访问外网时，可以在自己的电脑从CHARMM官方网站获取该固定版本，上传到实例包的 download-cache/，文件名与 CHARMM36-SOURCE.json一致，然后再次运行 python3 get_forcefield.py。本例原始压缩包SHA256是 c2b0c70e21cced150528ee25d205e5bd3367a5e55e85c2ff0403911980528151。

这里用CHARMM36m蛋白参数、预定义HEME和CHARMM修改的TIP3P。目录名含CGenFF5.0并不代表调用了CGenFF程序，也不代表未知配体已有参数。GROMACS2026.4使用GCC/G++13.3、CUDA12.8、AVX2编译；原子模型与粗粒化Martini不能混合。

2 His–Fe连接为何需要单独处理#

原始CHARMM PHEM规定一条His NE2–Fe键、六个跨组分角度和一条 His CD2–NE2–Fe–heme NA周期扭转，同时删除两个相对的卟啉N–Fe–N角，并明确不要在这个补丁之后自动生成角度和扭转。

附件 patch_heme.py严格针对本例His93/HEME拓扑执行这一连接定义：不改原子类型、电荷或已有键与1–4对，不自行增加零力常数；若原子身份或预期连接不匹配，就停止，让使用者先检查体系。原始CHARMM血红素拓扑与PHEM定义

任务内部先执行：

gmx pdb2gmx -f input.pdb -o processed.gro -p topol.top   -ff charmm36-feb2026_cgenff-5.0 -water tip3p -ignh -merge all
python3 ../../patch_heme.py
cat PHEM-PATCH.json
head -12 processed.gro

自动任务已经生成拓扑时不必再次运行这些命令。-merge all让蛋白与辅因子进入可包含共价连接的同一个分子；PHEM-PATCH.json记录对应连接项和参数来源。这里的脚本不是任意金属中心或配体的通用参数化工具，换配体、质子化或氧化态时必须重新选择匹配模型。

保留的周期扭转使用 k=0.2092 kJ/mol、相位0°、周期数4，即 V=k[1+cos(4φ)]。下图用同一真实构型，把血红素沿Fe–His轴刚性旋转0、15、30、45°，比较原生GROMACS的Proper Dih.能量变化与这一项的独立计算；不是结合自由能扫描。

参考构型和能量在 reference-controls/。在有GPU分配的Slurm任务中，可用原生命令重新提取：

gmx mdrun -s control.tpr -rerun rotations.trr -deffnm rerun   -ntmpi 1 -ntomp 1 -nb gpu -pme gpu -bonded cpu -update cpu
gmx energy -f rerun.edr -o proper-energy.xvg
# 选择Proper-Dih.，最后输入0

3 水、盐和动力学输入#

建盒、加水与四阶段命令顺序见二硫键完整流程。CHARMM离子名称在这里是 SOD和 CLA：

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 SOD -nname CLA -neutral -conc 0.15 -seed 13579
# 选择SOL

不同力场采用各自配套的原子类型、离子名称、水模型和非键处理。本例采用CHARMM的PME、1.2 nm截断、1.0–1.2 nm force-switch、不加色散尾修正；不要照搬上一节AMBER的1.0 nm设置。

最小化完整输入#

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.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = no
integrator = steep
emtol = 500
emstep = .001
nsteps = 5000
nstenergy = 100
nstlog = 100

emstep=0.001 nm是本例的初始最小化步长，用较小位移逐步松弛局部接触；它不是生产MD时间步。检查最小化的最大力、结构和约束输出后，再进入平衡。

生产完整输入#

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.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = no
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

NVT10 ps和NPT20 ps的完整输入分别在 mdp/myoglobin/nvt.mdp、npt.mdp，带蛋白重原子位置限制；50 ps生产不带位置限制。共价Fe连接和氢键约束仍由拓扑保留。仅仅看到Fe–His距离稳定，不能推断预测出了真实配位交换或结合强度。

4 提取铁配位距离与蛋白波动#

module load ase/3.29.0 gromacs/2026.4-v100-threadmpi
python3 ../../export_native.py
head -35 native-analysis/bond-distances.xvg
gmx distance -s md.tpr -f md.xtc   -select 'resnr 154 and name FE plus resnr 93 and name NE2'   -oall Fe-His93.xvg

本例血红素保留原始编号154；换结构时先查看自己的拓扑再改选择表达式。附件工具按真实原子ID提取Fe–NA、NB、NC、ND和His93 NE2的五条连接，避免默认组名遗漏HETATM辅因子。

gmx distance文本会舍入到0.001 nm；更细的比较可读取原始XTC或使用 gmx traj -ox导出坐标。MDAnalysis默认坐标单位为Å，转换成GROMACS的nm时除以10。RMSD和RMSF定义与二硫键分析一致，使用生产TPR参考。

四组数据来自两个速度种子在P100与V10016G上的重复流程。50 ps用于学习辅因子参数、拓扑连接和后处理，不是收敛的构象采样、结合能或金属化学研究。得到本例流程后，应根据实际研究的问题选择配位模型并延长采样。