01 — 蛋白质-配体分子动力学模拟完整流程

本文档详细介绍从 PDB 文件到完成 MD 模拟的全流程，包括：拓扑构建、小分子处理、溶剂化、能量最小化、平衡模拟和生产模拟。

1. 准备工作

所需文件

文件	说明	位置
`protein.pdb`	蛋白质结构文件（来自 PDB 数据库或同源建模）	工作目录
`unk.mol2`	小分子配体结构文件（来自药物设计软件）	工作目录
`em.mdp`	能量最小化参数文件	`docs/mdp配置文件/em.mdp`
`ions.mdp`	离子添加参数文件	`docs/mdp配置文件/ions.mdp`
`nvt.mdp`	NVT 平衡参数文件	`docs/mdp配置文件/nvt.mdp`
`npt.mdp`	NPT 平衡参数文件	`docs/mdp配置文件/npt.mdp`
`md.mdp`	生产模拟参数文件	`docs/mdp配置文件/md.mdp`

目录结构建议

project/
├── protein.pdb          # 蛋白质结构
├── unk.mol2             # 配体结构
├── em.mdp               # MDP 参数文件
├── ions.mdp
├── nvt.mdp
├── npt.mdp
├── md.mdp
├── charmm36-jul2022.ff/ # 力场文件夹
├── cgenff_charmm2gmx_py3_nx2.py  # 配体拓扑转换脚本
└── sort_mol2_bonds.pl            # MOL2 文件修复脚本

力场文件夹和 Python/Perl 脚本位于本指南的 important_file/ 目录下，请复制到项目工作目录中使用。

2. 构建蛋白质拓扑

# 激活 GROMACS 环境
conda activate gmx

# 由 PDB 生成 GROMACS 拓扑文件
gmx pdb2gmx -f protein.pdb -o protein_processed.gro -ignh -ter

参数说明

参数	含义
`-f protein.pdb`	输入的蛋白质结构文件
`-o protein_processed.gro`	输出的 GROMACS 结构文件
`-ignh`	忽略 PDB 中的氢原子，由 GROMACS 重新添加
`-ter`	交互式指定末端残基的质子化状态

交互步骤

执行后会提示选择力场和末端状态：

# 步骤 1：选择力场
# 输入 1 选择 CHARMM36 all-atom force field（推荐）
1

# 步骤 2-3：指定 N 末端状态
#  1: NH3+ (带电, 默认)
#  0: NH2 (中性)
1

# 步骤 4-5：指定 C 末端状态
#  1: COO- (带电, 默认)
#  0: COOH (中性)
0

生成的文件

文件	说明
`protein_processed.gro`	处理后的蛋白质结构
`topol.top`	体系拓扑文件（后续会持续修改）
`posre.itp`	蛋白质位置限制文件

3. 处理小分子配体

3.1 修复 MOL2 文件

CGenFF 对 MOL2 文件格式有严格要求，需要先修复：

perl sort_mol2_bonds.pl unk.mol2 unk_fix.mol2

3.2 在线生成配体拓扑

打开网站 https://cgenff.com/
上传 unk_fix.mol2 文件
下载生成的 unk.str 文件（包含该小分子的力场参数）

3.3 转换为 GROMACS 兼容格式

python cgenff_charmm2gmx_py3_nx2.py UNK unk_fix.mol2 unk_fix.str charmm36-jul2022.ff

注意：此脚本需要较高配置（建议 3080 及以上显卡），因为 CGenFF 参数计算需要大量计算资源。

生成的文件

文件	说明
`unk.itp`	配体拓扑文件（原子类型、成键信息）
`unk.prm`	配体参数文件（键参数、角度参数等）
`unk_ini.pdb`	配体初始结构

# 将配体 PDB 转换为 GRO 格式
gmx editconf -f unk_ini.pdb -o unk.gro

4. 组装复合物体系

4.1 手动合并蛋白质和配体

需要手动创建 complex.gro 文件：

打开 protein_processed.gro：保留所有行
打开 unk.gro：只保留原子坐标行（从第 3 行开始，到倒数第 2 行结束）
合并：在 protein_processed.gro 的倒数第二行（最后一行是盒子信息上方）按回车插入新行，粘贴配体的原子坐标行
修改原子总数：将第 2 行的原子数改为蛋白质原子数 + 配体原子数

示例：蛋白质有 5000 个原子，配体有 50 个原子，则第 2 行改为 5050

4.2 修改拓扑文件

在 topol.top 文件中添加配体相关信息：

# 编辑 topol.top，在合适位置插入以下内容

插入位置 1 — 在 [ moleculetype ] 上方（在 #include "forcefield.itp" 之后）：

; Include ligand parameters
#include "unk.prm"

插入位置 2 — 在 ; Include water topology 上方：

; Include ligand topology
#include "unk.itp"

修改文件末尾 — 在 [ molecules ] 部分添加：

UNK               1

验证：确保 SOL 和 UNK 不在同一行，如果它们在同一行，在 SOL 之后按回车换行。

4.3 生成复合物 GRO 文件

gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0

参数	含义
`-bt dodecahedron`	十二面体盒子（比立方体节省约 30% 溶剂分子）
`-d 1.0`	溶质与盒子边界的距离（nm）

5. 溶剂化与中和

5.1 添加溶剂（水分子）

gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

参数	含义
`-cs spc216.gro`	使用 SPC 水模型（与 CHARMM36 兼容）
`-p topol.top`	自动更新拓扑文件中的 SOL 数量

5.2 添加离子中和电荷

首先生成离子添加的 TPR 文件：

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr

然后添加 Na⁺/Cl⁻ 中和体系电荷并设置生理浓度：

gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname SOD -nname CLA -neutral

交互步骤：

# 选择溶剂（水分子）组来替换为离子
# 输入 SOL 对应的组号（通常是 15）
15

参数	含义
`-pname SOD`	正离子名称（Na⁺）
`-nname CLA`	负离子名称（Cl⁻）
`-neutral`	自动中和体系总电荷

6. 能量最小化

对溶剂化体系进行能量最小化，消除不合理的原子接触：

gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr

# 如果遇到警告，可尝试：
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr -maxwarn 1

# 执行能量最小化
gmx mdrun -v -deffnm em

# 如果遇到 CPU 调用错误（如 3080 及以上显卡），使用：
gmx mdrun -v -deffnm em -ntmpi 1 -ntomp 0

参数	含义
`-v`	显示详细输出
`-deffnm em`	所有输出文件前缀为 em
`-ntmpi 1`	使用 1 个 MPI 进程
`-ntomp 0`	自动设置 OpenMP 线程数

判断是否收敛

能量最小化完成后，检查最后一步的势能： - 查看 em.log 中的 Potential Energy，应在负值范围内（通常 -10⁵ ~ -10⁶ kJ/mol） - 或查看 em.gro 是否存在且不为空

7. 位置限制

对配体施加位置限制，防止在平衡阶段发生明显位移。

7.1 创建配体索引

gmx make_ndx -f unk.gro -o index_unk.ndx

交互：

0 & ! a H*    # 选择配体中非氢原子
q             # 退出

7.2 生成位置限制文件

gmx genrestr -f unk.gro -n index_unk.ndx -o posre_unk.itp -fc 1000 1000 1000

交互：选择配体组（通常为 3）

7.3 在拓扑中添加位置限制

编辑 topol.top，在 ; Include water topology 上方插入：

; Ligand position restraints
#ifdef POSRES
#include "posre_unk.itp"
#endif

7.4 创建完整索引文件

gmx make_ndx -f em.gro -o index.ndx

交互：

1 | 13      # 合并蛋白质(1)和小分子(13)为 Protein_UNK
14 | 15     # 合并 SOL(14) 和 SOD(15) 为 SOD_Water
q

注意：组号可能因体系不同而变化，请根据实际编号调整。

8. NVT 平衡

在恒定粒子数、体积、温度下进行 100 ps 的平衡模拟：

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

重要：修改 NVT 参数文件

编辑 nvt.mdp 中的 tc-grps 行，确保温度耦合组与索引文件中的组名匹配。根据上一步创建的索引：

tc-grps                 = Protein_UNK SOD_Water

必须将默认的组名改为索引文件中实际存在的组名，否则会报错。

# GPU 加速执行 NVT 平衡
gmx mdrun -v -deffnm nvt

9. NPT 平衡

在恒定粒子数、压力、温度下进行 100 ps 的平衡模拟：

gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr

# GPU 加速执行
gmx mdrun -v -deffnm npt

参数	含义
`-t nvt.cpt`	从 NVT 的检查点文件继续
`-r nvt.gro`	参考结构（用于位置限制）

10. 生产模拟 (MD)

10.1 运行生产模拟

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md100.tpr

# 前台运行
gmx mdrun -v -deffnm md100

md.mdp 中的 nsteps = 50000000，配合 dt = 0.002 ps，总模拟时长为 100 ns。

10.2 防止断连（推荐）

使用 nohup 在后台运行，即使关闭终端也不会中断：

nohup gmx mdrun -v -deffnm md100 > md100.log 2>&1 &

10.3 监控进度

# 实时查看模拟日志
tail -f md100.log

# 检查进程是否在运行
ps aux | grep mdrun

11. 附录：MDP 配置文件说明

所有 MDP 文件位于 docs/mdp配置文件/ 目录下。

`em.mdp` — 能量最小化

integrator  = steep       # 最陡下降法
emtol       = 1000.0      # 最大力收敛阈值 (kJ/mol/nm)
emstep      = 0.01        # 初始步长 (nm)
nsteps      = 50000       # 最大步数

`ions.mdp` — 离子添加（仅用于 grompp）

; 极简配置，仅用来生成 TPR 以便 genion 识别体系
coulombtype = cutoff      # 快速 cutoff 即可，此处不重要

`nvt.mdp` — NVT 平衡

define      = -DPOSRES    # 启用位置限制
integrator  = md          # Leap-frog 积分器
nsteps      = 50000       # 100 ps (0.002 ps/step × 50000)
dt          = 0.002       # 积分步长 2 fs
tcoupl      = V-rescale   # 速度重缩放恒温器
tc-grps     = Protein_UNK SOD_Water
ref_t       = 300  300    # 目标温度 (K)

`npt.mdp` — NPT 平衡

continuation = yes        # 从 NVT 继续
pcoupl      = Parrinello-Rahman  # 各向同性压力耦合
tau_p       = 2.0         # 压力耦合时间常数 (ps)
ref_p       = 1.0         # 目标压力 (bar)
compressibility = 4.5e-5  # 水的等温压缩率 (bar⁻¹)

`md.mdp` — 生产模拟

nsteps      = 50000000    # 100 ns (0.002 ps × 50,000,000)
nstenergy   = 5000        # 每 10 ps 输出一次能量
nstxout-compressed = 5000 # 每 10 ps 输出一次压缩轨迹
continuation = yes        # 从 NPT 继续

调整模拟时长：修改 md.mdp 中的 nsteps。例如 50 ns → nsteps = 25000000，200 ns → nsteps = 100000000。

完整执行流程速查

# 1. 蛋白质拓扑
gmx pdb2gmx -f protein.pdb -o protein_processed.gro -ignh -ter

# 2. 配体处理
perl sort_mol2_bonds.pl unk.mol2 unk_fix.mol2
# [上传 unk_fix.mol2 到 https://cgenff.com/，下载 unk.str]
python cgenff_charmm2gmx_py3_nx2.py UNK unk_fix.mol2 unk_fix.str charmm36-jul2022.ff
gmx editconf -f unk_ini.pdb -o unk.gro

# 3. 手动合并 complex.gro + 修改 topol.top（见第 4 节）

gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname SOD -nname CLA -neutral

# 4. 能量最小化
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

# 5. 配体位置限制
gmx make_ndx -f unk.gro -o index_unk.ndx
gmx genrestr -f unk.gro -n index_unk.ndx -o posre_unk.itp -fc 1000 1000 1000
# 编辑 topol.top 添加 posre_unk.itp 引用
gmx make_ndx -f em.gro -o index.ndx

# 6. NVT 平衡
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -v -deffnm nvt

# 7. NPT 平衡
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr
gmx mdrun -v -deffnm npt

# 8. 生产模拟
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md100.tpr
nohup gmx mdrun -v -deffnm md100 > md100.log 2>&1 &

01 — 蛋白质-配体分子动力学模拟完整流程

目录

1. 准备工作

所需文件

目录结构建议

2. 构建蛋白质拓扑

参数说明

交互步骤

生成的文件

3. 处理小分子配体

3.1 修复 MOL2 文件

3.2 在线生成配体拓扑

3.3 转换为 GROMACS 兼容格式

生成的文件

4. 组装复合物体系

4.1 手动合并蛋白质和配体

4.2 修改拓扑文件

4.3 生成复合物 GRO 文件

5. 溶剂化与中和

5.1 添加溶剂（水分子）

5.2 添加离子中和电荷

6. 能量最小化

判断是否收敛

7. 位置限制

7.1 创建配体索引

7.2 生成位置限制文件

7.3 在拓扑中添加位置限制

7.4 创建完整索引文件

8. NVT 平衡

重要：修改 NVT 参数文件

9. NPT 平衡

10. 生产模拟 (MD)

10.1 运行生产模拟

10.2 防止断连（推荐）

10.3 监控进度

11. 附录：MDP 配置文件说明

em.mdp — 能量最小化

ions.mdp — 离子添加（仅用于 grompp）

nvt.mdp — NVT 平衡

npt.mdp — NPT 平衡

md.mdp — 生产模拟

完整执行流程速查

💬 留言

`em.mdp` — 能量最小化

`ions.mdp` — 离子添加（仅用于 grompp）

`nvt.mdp` — NVT 平衡

`npt.mdp` — NPT 平衡

`md.mdp` — 生产模拟