适合对象：GROMACS 初学者、刚开始学习蛋白–配体分子动力学模拟的同学。
平台：AutoDL Linux 服务器
软件：GROMACS 2024.6 CUDA 版
任务：蛋白–配体复合物体系构建、GPU 分子动力学模拟与结果可视化分析

---

1. 实验目标

本次工作的目标是在 AutoDL 服务器上完成一次蛋白–配体复合物的分子动力学模拟。完整流程包括：

1. 准备蛋白受体和小分子配体文件；
2. 使用 CHARMM36 力场生成蛋白拓扑；
3. 使用 CGenFF 生成小分子配体参数；
4. 合并蛋白和配体坐标；
5. 构建水盒子并加入离子；
6. 进行能量最小化；
7. 进行 NVT 和 NPT 平衡；
8. 使用 GPU 进行正式 MD 模拟；
9. 对轨迹结果进行 RMSD、RMSF、回旋半径、SASA、氢键、PCA 和自由能形貌图分析。

本次最终使用的环境如下：

服务器平台：AutoDL
操作系统：Linux
模拟软件：GROMACS 2024.6 CUDA 版
GPU 后端：CUDA
蛋白力场：CHARMM36-feb2026_cgenff-5.0
水模型：CHARMM-modified TIP3P
配体参数：CGenFF
正式 MD 输出前缀：md_gpu
测试模拟时长：1 ns

---

2. 工作目录准备

在 AutoDL 中，建议将项目文件统一放在 /root/autodl-tmp 下，避免占用系统盘空间。

cd /root/autodl-tmp
mkdir -p gromacs_tutorial
cd gromacs_tutorial

本次使用的初始文件包括：

rep.pre.pdb     原始受体文件
rep.pdb         补全原子后的受体文件
lig.sdf         从 PubChem 下载的配体 SDF 文件

其中，rep.pdb 是后续用于 GROMACS 建模的受体文件，lig.sdf 是后续配体参数化的起点。

需要注意：PubChem 下载的配体结构通常是游离小分子构象，不一定处于蛋白结合口袋中。如果用于正式蛋白–配体复合物研究，建议先进行分子对接，然后使用对接后的配体构象。

---

3. 生成蛋白受体拓扑

使用 GROMACS 的 pdb2gmx 处理受体蛋白：

gmx pdb2gmx -f rep.pdb -o rep_gmx.gro -p topol.top -ignh

在力场选择界面中选择 CHARMM 力场：

CHARMM all-atom force field

实际使用的是：

Charmm36-feb2026_cgenff-5.0

水模型选择：

TIP3P，CHARMM-modified TIP3P water model

该步骤成功后会生成：

rep_gmx.gro     处理后的蛋白坐标文件
topol.top       系统拓扑主文件
posre.itp       蛋白位置限制文件

如果日志中出现类似 Long Bond 的提示，说明蛋白结构中可能存在异常长键、断链或不合理连接。为了流程复现，可以暂时使用 -maxwarn 1 继续；但如果用于正式科研，建议回到原始结构中修复该问题。

---

4. 配体结构准备与格式转换

原始教学数据中的 lig.mol2 可能来自 PyMOL 导出，常见问题包括：

1. 原子名大量重复，例如 C、C、C；
2. 残基名不规范；
3. 碳原子几乎全部被识别为 C.3；
4. 键级和芳香性识别不准确。

因此，本次改用 PubChem 下载的 SDF 文件：

lig.sdf

使用 Open Babel 转换为 CGenFF 可识别的 mol2 文件：

obabel lig.sdf -O lig_cgenff.mol2 -h --partialcharge gasteiger

该命令的作用是：

1. 将 SDF 转换为 MOL2；
2. 补全氢原子；
3. 生成初始 Gasteiger 部分电荷；
4. 尽量保留合理的键级、芳香性和原子类型。

转换后可以检查文件：

head -40 lig_cgenff.mol2

如果看到类似下面的原子类型，说明结构信息较为正常：

C.ar
O.2
O.3
C.3

---

5. 规范化配体 mol2 文件

为了让 CGenFF 更容易识别，需要将分子名、残基名和原子名统一为规范格式。

最终生成文件：

lig_cgenff_fixed.mol2

该文件中应满足：

分子名：LIG
残基名：LIG
原子名：O1、O2、C1、C2、H1 等唯一名称

检查文件：

head -40 lig_cgenff_fixed.mol2

应看到：

@<TRIPOS>MOLECULE
LIG

ATOM 部分应类似：

1 O1 ...
2 O2 ...
17 C1 ...
18 C2 ...

---

6. 使用 CGenFF 生成配体参数

将 lig_cgenff_fixed.mol2 上传到 CGenFF 网站。

网页参数设置如下：

Use Version: 5.0
Job Name: LIG
Output Mol2: 开启
Include Warnings: 开启
Include Debug Info: 可开可不开
Copy Existing Params: 关闭
Net Charge: 0

运行 CGenFF 后下载结果文件，其中最重要的是：

lig_cgenff.str

上传回 AutoDL 后重命名为：

mv lig_cgenff.str lig.str

---

7. 使用 cgenff_charmm2gmx.py 转换配体拓扑

准备转换脚本：

unzip cgenff_charmm2gmx-main.zip
cp cgenff_charmm2gmx-main/cgenff_charmm2gmx.py .

执行转换：

FFDIR=/root/miniconda3/envs/gromacs/share/gromacs/top/charmm36-feb2026_cgenff-5.0.ff

python cgenff_charmm2gmx.py LIG lig_cgenff_fixed.mol2 lig.str "$FFDIR"

转换成功后生成：

lig.itp         配体拓扑文件
lig.prm         配体参数文件
lig.top         配体单独拓扑文件
lig_ini.pdb     配体坐标文件

将配体 PDB 坐标转换为 GROMACS 的 .gro 格式：

gmx editconf -f lig_ini.pdb -o lig.gro

---

8. 修改系统拓扑文件 topol.top

接下来需要把配体参数加入系统拓扑文件 topol.top。

首先确认配体分子名：

grep -A 5 "\[ moleculetype \]" lig.itp

结果应类似：

[ moleculetype ]
; Name     nrexcl
LIG        3

其中 3 是 nrexcl，表示排除近邻非键相互作用的级数，不是配体数量。

8.1 加入 lig.prm

在 topol.top 中找到：

#include "charmm36-feb2026_cgenff-5.0.ff/forcefield.itp"

在其下方加入：

#include "lig.prm"

最终应为：

#include "charmm36-feb2026_cgenff-5.0.ff/forcefield.itp"
#include "lig.prm"

lig.prm 必须放在任何 [ moleculetype ] 定义之前，因为其中包含配体的 atomtypes、bondtypes、angletypes 和 dihedraltypes 等参数。

8.2 加入 lig.itp

查找水模型位置：

grep -n "tip3p.itp" topol.top

在：

#include "charmm36-feb2026_cgenff-5.0.ff/tip3p.itp"

之前加入：

#include "lig.itp"

最终应为：

#include "lig.itp"
#include "charmm36-feb2026_cgenff-5.0.ff/tip3p.itp"

8.3 在 [ molecules ] 中加入配体

在 topol.top 文件底部找到 [ molecules ] 部分：

[ molecules ]
; Compound        #mols
Protein_chain_A     1

加入配体：

LIG                 1

最终应为：

[ molecules ]
; Compound        #mols
Protein_chain_A     1
LIG                 1

这里必须写 LIG 1，表示体系中有 1 个配体分子。

---

9. 合并蛋白和配体坐标

将蛋白坐标 rep_gmx.gro 和配体坐标 lig.gro 合并为复合物坐标文件 complex.gro。

创建脚本：

nano merge_gro.py

写入：

protein_gro = "rep_gmx.gro"
ligand_gro = "lig.gro"
output_gro = "complex.gro"

with open(protein_gro) as f:
    p_lines = f.readlines()

with open(ligand_gro) as f:
    l_lines = f.readlines()

p_atoms = p_lines[2:-1]
l_atoms = l_lines[2:-1]
box = p_lines[-1]

total_atoms = len(p_atoms) + len(l_atoms)

with open(output_gro, "w") as f:
    f.write("Protein-ligand complex\n")
    f.write(f"{total_atoms:5d}\n")
    for line in p_atoms:
        f.write(line)
    for line in l_atoms:
        f.write(line)
    f.write(box)

print(f"Generated {output_gro} with {total_atoms} atoms.")

运行：

python merge_gro.py

生成：

complex.gro

需要说明的是，这种方式只是简单拼接两个 .gro 文件。若进行严格蛋白–配体复合物模拟，配体坐标应来自分子对接结果，而不是 PubChem 游离构象。

---

10. 定义盒子、加水和加离子

10.1 定义盒子

gmx editconf -f complex.gro -o complex_newbox.gro -c -d 1.0 -bt cubic

10.2 加水

gmx solvate -cp complex_newbox.gro -cs spc216.gro -o complex_solv.gro -p topol.top

加水成功后，topol.top 中会自动加入水分子数量。例如：

SOL             33090

10.3 准备 ions.mdp

创建 ions.mdp：

integrator  = steep
emtol       = 1000.0
emstep      = 0.01
nsteps      = 50000

cutoff-scheme = Verlet
nstlist     = 10
rcoulomb    = 1.0
rvdw        = 1.0
pbc         = xyz

10.4 生成加离子输入文件

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

10.5 加离子

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

这里使用：

SOD：Na+
CLA：Cl-

因为当前体系采用 CHARMM36 力场。运行时选择 SOL 组，用离子替换部分水分子。

最终生成：

complex_ions.gro

---

11. 能量最小化

创建 minim.mdp：

integrator              = steep
emtol                   = 1000.0
emstep                  = 0.01
nsteps                  = 50000

cutoff-scheme           = Verlet
nstlist                 = 20
rlist                   = 1.2

coulombtype             = PME
rcoulomb                = 1.2

vdwtype                 = cutoff
vdw-modifier            = force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2

pbc                     = xyz

生成能量最小化输入文件：

gmx grompp -f minim.mdp -c complex_ions.gro -p topol.top -o em.tpr -maxwarn 1

这里使用 -maxwarn 1 是因为体系存在一个长键 warning。流程复现阶段可以继续，但正式模拟应修复结构。

运行能量最小化：

gmx mdrun -v -s em.tpr -deffnm em

完成后生成：

em.gro
em.edr
em.log
em.trr
em.tpr

检查势能：

gmx energy -f em.edr -o potential.xvg

选择：

Potential

---

12. NVT 平衡

创建 nvt.mdp：

define                  = -DPOSRES
integrator              = md
dt                      = 0.002
nsteps                  = 50000

nstxout-compressed      = 500
nstenergy               = 500
nstlog                  = 500

continuation            = no
constraint_algorithm    = lincs
constraints             = h-bonds
lincs_iter              = 1
lincs_order             = 4

cutoff-scheme           = Verlet
nstlist                 = 20
rlist                   = 1.2

vdwtype                 = cutoff
vdw-modifier            = force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2

coulombtype             = PME
rcoulomb                = 1.2
pme_order               = 4
fourierspacing          = 0.16

tcoupl                  = V-rescale
tc-grps                 = Protein Non-Protein
tau_t                   = 0.1 0.1
ref_t                   = 300 300

pcoupl                  = no
pbc                     = xyz
DispCorr                = no

gen_vel                 = yes
gen_temp                = 300
gen_seed                = -1

生成并运行 NVT：

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr -maxwarn 1
gmx mdrun -v -s nvt.tpr -deffnm nvt

本次 NVT 设置为：

50000 steps × 0.002 ps = 100 ps

NVT 的作用是在固定体积条件下让体系温度稳定到 300 K。

检查温度：

gmx energy -f nvt.edr -o temperature.xvg

选择：

Temperature

---

13. NPT 平衡

创建 npt.mdp：

define                  = -DPOSRES
integrator              = md
dt                      = 0.002
nsteps                  = 50000

nstxout-compressed      = 500
nstenergy               = 500
nstlog                  = 500

continuation            = yes
constraint_algorithm    = lincs
constraints             = h-bonds
lincs_iter              = 1
lincs_order             = 4

cutoff-scheme           = Verlet
nstlist                 = 20
rlist                   = 1.2

vdwtype                 = cutoff
vdw-modifier            = force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2

coulombtype             = PME
rcoulomb                = 1.2
pme_order               = 4
fourierspacing          = 0.16

tcoupl                  = V-rescale
tc-grps                 = Protein Non-Protein
tau_t                   = 0.1 0.1
ref_t                   = 300 300

pcoupl                  = C-rescale
pcoupltype              = isotropic
tau_p                   = 2.0
ref_p                   = 1.0
compressibility         = 4.5e-5

pbc                     = xyz
DispCorr                = no

gen_vel                 = no

生成并运行 NPT：

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr -maxwarn 1
gmx mdrun -v -s npt.tpr -deffnm npt

NPT 的作用是在恒温恒压条件下平衡体系密度和盒子大小。

完成后生成：

npt.gro
npt.cpt
npt.edr
npt.log
npt.xtc

---

14. 编译并使用 GPU 版 GROMACS

最初 conda 安装的 GROMACS 2026.0 是 OpenCL 版，不能在 AutoDL 上正常调用 NVIDIA GPU。因此后续编译了 GROMACS 2024.6 CUDA 版。

AutoDL 当前 CUDA 版本为：

nvcc --version

显示：

CUDA 11.8

GROMACS 2026.0 要求 CUDA 12.1 及以上，因此选择编译支持 CUDA 11.8 的 GROMACS 2024.6。

最终安装目录为：

/root/autodl-tmp/gromacs_cuda_2024

每次使用 GPU 版 GROMACS 前，需要加载环境：

source /root/autodl-tmp/gromacs_cuda_2024/bin/GMXRC
hash -r
which gmx
gmx --version

应确认：

Executable: /root/autodl-tmp/gromacs_cuda_2024/bin/gmx
GROMACS version: 2024.6
GPU support: CUDA

由于 GROMACS 2024 和 2026 的 .tpr 文件不建议混用，所以正式 MD 前重新生成 2024 版输入文件：

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_2024.tpr -maxwarn 1

如果 GROMACS 2024 找不到 charmm36-feb2026_cgenff-5.0.ff，需要将该力场目录复制到当前项目目录：

cp -r /root/miniconda3/envs/gromacs/share/gromacs/top/charmm36-feb2026_cgenff-5.0.ff .

---

15. 正式 MD 模拟

正式 MD 使用 md.mdp，本次先进行 1 ns 测试模拟：

dt = 0.002 ps
nsteps = 500000
总时长 = 500000 × 0.002 ps = 1000 ps = 1 ns

运行 GPU 版正式 MD：

gmx mdrun -v -s md_2024.tpr -deffnm md_gpu -nb gpu

日志显示：

1 GPU selected for this run.
PP:0,PME:0
PP tasks will do short-ranged interactions on the GPU
PP task will update and constrain coordinates on the GPU
PME tasks will do all aspects on the GPU

说明正式 MD 已经成功使用 GPU 加速。CPU 仍参与调度和部分辅助计算，但主要计算任务已经分配到 GPU 上。

运行完成后生成：

md_gpu.xtc     轨迹文件
md_gpu.gro     最后一帧结构
md_gpu.edr     能量文件
md_gpu.log     日志文件
md_gpu.cpt     断点续跑文件

检查是否完成：

tail -40 md_gpu.log
grep -i "writing final coordinates\|performance\|finished" md_gpu.log

如果日志中出现：

Writing final coordinates
Performance

说明正式 MD 已完成。

---

16. 轨迹预处理

创建分析目录：

mkdir -p analysis

去除周期性边界跳跃：

gmx trjconv -s md_2024.tpr -f md_gpu.xtc -o analysis/md_nojump.xtc -pbc nojump

选择：

System

蛋白居中：

gmx trjconv -s md_2024.tpr -f analysis/md_nojump.xtc -o analysis/md_center.xtc -pbc mol -center -ur compact

选择：

center group: Protein
output group: System

后续所有分析主要使用：

analysis/md_center.xtc

---

17. 创建 index 文件

为了后续选择配体分组，创建 index 文件：

gmx make_ndx -f md_gpu.gro -o analysis/index.ndx

交互输入：

r LIG
q

这样可以确保后续分析中有 LIG 分组。

---

18. RMSD 分析

蛋白 RMSD：

gmx rms -s md_2024.tpr -f analysis/md_center.xtc -o analysis/rmsd_protein.xvg -tu ns

选择：

Least-squares fit group: Backbone
RMSD calculation group: Backbone

配体 RMSD：

gmx rms -s md_2024.tpr -f analysis/md_center.xtc -n analysis/index.ndx -o analysis/rmsd_ligand.xvg -tu ns

选择：

Least-squares fit group: Protein
RMSD calculation group: LIG

蛋白 RMSD 用于判断蛋白整体结构是否稳定；配体 RMSD 用于判断配体是否保持相对稳定的结合构象。

---

19. RMSF 分析

gmx rmsf -s md_2024.tpr -f analysis/md_center.xtc -o analysis/rmsf_residue.xvg -res

选择：

Protein

RMSF 用于观察各个残基的柔性。峰值较高的区域通常对应末端、loop 区域或运动较大的结构片段。

---

20. 回旋半径分析

gmx gyrate -s md_2024.tpr -f analysis/md_center.xtc -o analysis/gyrate.xvg

选择：

Protein

回旋半径用于评估蛋白整体紧凑程度。如果曲线较稳定，说明模拟过程中蛋白没有明显展开或塌缩。

---

21. SASA 分析

gmx sasa -s md_2024.tpr -f analysis/md_center.xtc -o analysis/sasa.xvg

选择：

Protein

SASA 表示溶剂可及表面积，用于评估蛋白暴露面积随时间的变化。

---

22. 蛋白–配体氢键分析

gmx hbond -s md_2024.tpr -f analysis/md_center.xtc -n analysis/index.ndx -num analysis/hbnum_protein_ligand.xvg

选择：

group 1: Protein
group 2: LIG

如果新版 gmx hbond 命令不兼容，可使用：

gmx hbond-legacy -s md_2024.tpr -f analysis/md_center.xtc -n analysis/index.ndx -num analysis/hbnum_protein_ligand.xvg

氢键数量可以反映蛋白和配体之间相互作用是否持续稳定。

---

23. PCA 和自由能形貌图

先计算协方差矩阵和特征向量：

gmx covar -s md_2024.tpr -f analysis/md_center.xtc -o analysis/eigenval.xvg -v analysis/eigenvec.trr -av analysis/average.pdb

选择：

fit group: C-alpha 或 Backbone
analysis group: C-alpha 或 Backbone

然后生成 PC1 和 PC2 的二维投影：

gmx anaeig -s md_2024.tpr -f analysis/md_center.xtc -v analysis/eigenvec.trr -first 1 -last 2 -2d analysis/pc12_2d.xvg

使用 Python 脚本将 pc12_2d.xvg 转换为自由能形貌图：

analysis/free_energy_landscape.png

自由能形貌图基于：

F = -kBT ln P

其中低自由能区域代表模拟中更常出现、更稳定的构象状态。

---

24. 结果图整理

最终分析结果保存在：

analysis/

主要结果包括：

rmsd_protein.png
rmsd_ligand.png
rmsf_residue.png
gyrate.png
sasa.png
hbnum_protein_ligand.png
eigenval.png
free_energy_landscape.png

打包下载：

zip -r md_analysis_results.zip analysis

---

25. 本次流程的关键经验

本次流程的核心是：先用 CHARMM36 处理蛋白，再用 CGenFF 处理配体，最后在 GROMACS 中合并体系并完成 MD 模拟。对于初学者来说，最容易出错的地方主要有三个。

第一，配体 mol2 文件必须可靠。PyMOL 导出的 mol2 往往存在原子名重复、键级不准确和芳香性错误等问题，不适合直接用于 CGenFF。更稳妥的方法是从 PubChem 下载 SDF，再用 Open Babel 转为 mol2。

第二，蛋白和配体的力场要尽量一致。本次蛋白使用 CHARMM36，因此配体也使用 CGenFF，而不是混用 GAFF 或其他力场。

第三，GROMACS 版本和 CUDA 版本要匹配。AutoDL 当前 CUDA 是 11.8，因此 GROMACS 2026.0 CUDA 版无法编译，最终选择 GROMACS 2024.6 CUDA 版，并成功调用 GPU。

---

26. 本次模拟的局限性

本次流程已经完整跑通，但如果用于正式科研，还需要进一步改进。

首先，目前正式 MD 只跑了 1 ns，更适合作为流程验证，不适合得出强科学结论。正式分析通常需要更长时间，例如 10 ns、50 ns 或更长，并建议设置重复模拟。

其次，当前配体坐标来自 PubChem 游离构象，不一定在蛋白结合口袋中。正式蛋白–配体 MD 应先进行分子对接，将配体放入结合口袋，再进行动力学模拟。

第三，蛋白结构存在 long bond warning。复现流程时可以使用 -maxwarn 1 暂时继续，但正式研究应检查对应残基和原子，修复断链、缺失残基或不合理连接。

---

27. 总结

本次在 AutoDL 平台上完成了基于 GROMACS 的蛋白–配体分子动力学模拟。首先对蛋白受体进行预处理，并采用 CHARMM36/CGenFF 力场和 CHARMM-modified TIP3P 水模型生成蛋白拓扑。随后从 PubChem 获取配体结构，经 Open Babel 转换和格式规范化后，使用 CGenFF 生成配体参数，并通过 cgenff_charmm2gmx.py 转换为 GROMACS 可用的 lig.itp 和 lig.prm 文件。

之后，将蛋白和配体坐标合并，构建复合物体系，依次完成加盒子、加水、加离子、能量最小化、NVT 平衡和 NPT 平衡。最后，针对 AutoDL 的 CUDA 11.8 环境，编译并使用 GROMACS 2024.6 CUDA 版完成 1 ns GPU 加速正式 MD 模拟。

模拟结束后，对轨迹进行了周期性边界校正和蛋白居中处理，并进一步完成 RMSD、RMSF、回旋半径、SASA、蛋白–配体氢键、PCA 和自由能形貌图等后处理分析。该流程完整覆盖了蛋白–配体分子动力学模拟从体系构建到结果可视化的主要步骤，为后续开展更长时间尺度、更严格的酶–底物复合物动力学分析奠定了基础。