Zhang Xin's Tech Blog

Back
/ Update
13 min read
中文
分子动力学模拟 /
gromacs

分子动力学模拟：Gromacs教程1-水中的水解酶

本文是对Gromacs官方教程的中文翻译，并且包含了我的实操笔记
views | comments

本文是对Gromacs官方教程第一课Tutorial 1: Lysozyme in Water的中文翻译，是溶菌酶在水溶液中的模拟。包含了：

  • 蛋白质预处理和溶剂化
  • 能量最小化和平衡
  • 基本的MD模拟流程
  • 结果分析方法

第一阶段：打好基础#

1. 环境搭建#

我自己是用的Autodl的云服务器，你也可以按照如下的方法自己配置

操作系统：Linux 是首选（Ubuntu 是最友好的选择）。你可以在物理机安装、虚拟机安装，或者使用 Windows 的 WSL2。

  • 推荐：对于Windows用户，使用 WSL2 是最佳平衡方案。安装教程很容易在网上找到。

安装 GROMACS：

  • 对于 Ubuntu/WSL，最简单的方式是用 apt： 
    sudo apt update
sudo apt install gromacs
    bash
  • 如果你想用最新版，需要从官网下载源码编译（稍复杂，可后续再学）。

验证安装：

gmx -version
bash

如果显示出版本信息，恭喜你，第一步成功了！

2. 认识一个标准的模拟流程#

这是一个经典的 MD 模拟流程：

流程图
100%

正在加载流程图；也可展开下方源码查看。

查看流程图源码
graph TD
    A[蛋白质/小分子<br>初始结构] --> B[预处理<br>补全缺失原子/残基];
    B --> C[构建拓扑<br>定义分子连接与力场];
    D[力场文件<br>e.g., charmm36, amber99sb-ildn] --> C;
    E[溶剂盒子<br>e.g., 立方体， 十二面体] --> F[溶剂化<br>将分子放入盒子];
    C --> F;
    F --> G[添加离子<br>中和体系并设定离子浓度];
    G --> H[能量最小化<br>消除原子间的近距离接触];
    H --> I[NVT 平衡<br>稳定体系温度];
    I --> J[NPT 平衡<br>稳定体系密度和压力];
    J --> K[成品模拟<br>生产用于分析的数据];
    K --> L[轨迹分析<br>提取有意义的科学信息];

第二阶段：你的第一个模拟——水中的溶菌酶#

让我们按照上述流程，亲手完成一个最简单的模拟。这是 GROMACS 的 “Hello World”。

步骤 1: 准备结构文件#

  1. 从 PDB 数据库下载溶菌酶的pdb结构： 1AKI ↗
  2. 用文本编辑器（如 vim, nano）或可视化软件（如 PyMOL, VMD）打开它，删除所有的水分子和无关的离子。只保留蛋白质链 A。保存为 1aki_clean.pdb。

在这里我使用的是PyMOL得到干净的pdb文件

# 使用PyMOL打开1AKL.pdb文件
# 移除水分子
remove resn HOH
# 移除无关离子
remove not (polymer or resn UNK)
# 保存蛋白
save 1aki_clean.pdb
plaintext

步骤 2: 生成拓扑和力场#

拓扑文件描述了分子由哪些原子组成，它们如何连接，以及使用了什么力场。

# 使用 pdb2gmx 工具生成拓扑
gmx pdb2gmx -f 1aki_clean.pdb -o 1aki_processed.gro -p topol.top -water spce
bash
  • -f： 输入文件
  • -o： 处理后的结构文件（GROMACS 格式）
  • -p： 输出的拓扑文件
  • -water： 指定水模型（这里用 SPC/E）
  • 此时，终端会让你选择一个力场。对于初学者，输入 15 选择 CHARMM36 力场是一个不错的选择。

步骤 3: 定义模拟盒子并填充水#

# 创建一个立方体盒子，边界距离蛋白质至少 1.0 nm
gmx editconf -f 1aki_processed.gro -o 1aki_newbox.gro -c -d 1.0 -bt cubic

# 用水分子填充盒子
gmx solvate -cp 1aki_newbox.gro -cs spc216.gro -o 1aki_solvated.gro -p topol.top
bash
  • editconf -c： 将蛋白质置于盒子中心
  • solvate -cs： 指定水模型的结构文件（SPC 水）

步骤 4: 添加离子以中和电荷并设定盐浓度#

首先，需要生成一个 .tpr 输入文件，用于 grompp，然后使用 genion。

# 1. 准备一个离子参数文件（.mdp）。你可以从 GROMACS 官网教程或安装目录下找到。
# 下载 ‘ions.mdp’ 文件到当前目录。
wget http://www.mdtutorials.com/gmx/lysozyme/Files/ions.mdp

# 2. 生成 .tpr 文件
gmx grompp -f ions.mdp -c 1aki_solvated.gro -p topol.top -o ions.tpr

# 3. 添加离子
gmx genion -s ions.tpr -o 1aki_solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15
bash
  • -neutral： 添加钠离子或氯离子使体系整体电荷为中性。
  • -conc 0.15： 添加盐至生理浓度 0.15 M。
  • 此时，终端会问你“选择要替换的分子组”，输入 13（代表 “SOL”， 水分子）。程序会用离子替换掉一些水分子。

步骤 5: 能量最小化#

消除原子间的过度靠近，防止模拟开始时系统“爆炸”。

# 1. 准备能量最小化参数文件（minim.mdp）
wget http://www.mdtutorials.com/gmx/lysozyme/Files/minim.mdp
# 2. 生成 .tpr 文件
gmx grompp -f minim.mdp -c 1aki_solv_ions.gro -p topol.top -o em.tpr
# 3. 运行能量最小化
gmx mdrun -v -deffnm em
bash
  • -deffnm： 为所有输入和输出文件设置一个默认名称。

检查结果：

# 查看能量最小化是否收敛
gmx energy -f em.edr -o potential.xvg
# 当提示时，输入 `10 0` 选择 “Potential” 和 “零”。
# 用 xmgrace 或 gnuplot 绘图查看能量是否下降并趋于平稳。
bash

步骤 6 & 7: 平衡——NVT 和 NPT#

让系统在设定的温度和压力下达到平衡。

# NVT 平衡 (恒定粒子数、体积、温度)
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt

# 检查温度平衡
gmx energy -f nvt.edr -o temperature.xvg

# NPT 平衡 (恒定粒子数、压力、温度)
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt

# 检查压力平衡
gmx energy -f npt.edr -o pressure.xvg
# 检查密度平衡
gmx energy -f npt.edr -o density.xvg
bash

步骤 8: 成品 MD 模拟#

现在系统已经平衡好，可以开始长时间模拟来收集数据了。

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
gmx mdrun -deffnm md_0_1
bash

这个步骤可能会运行几分钟到几周，取决于你的体系大小和计算机性能。

第三阶段：分析结果#

模拟跑完了，我们得到了轨迹文件（.trr 或 .xtc）和能量文件（.edr）。现在来分析它。

  1. 回看轨迹:

    # 注意，这个需要图形界面支持
gmx view -f md_0_1.xtc -s md_0_1.tpr
    bash

  2. 计算 RMSD（蛋白质骨架稳定性）:

    gmx rms -s md_0_1.tpr -f md_0_1.xtc -o rmsd.xvg -tu ns
# 选择 “Backbone” 作为计算和拟合的组 (两次都选 4)
    bash

  3. 计算 RMSF（残基柔性）:

    gmx rmsf -s md_0_1.tpr -f md_0_1.xtc -o rmsf.xvg -res
# 选择 “Backbone” (4)
    bash

  4. 计算回旋半径（衡量蛋白质紧实度）:

    gmx gyrate -s md_0_1.tpr -f md_0_1.xtc -o gyrate.xvg
# 选择 “Protein” (1)
    bash

课后说明#

.mdp 参数文件： 我们简要说明mdp文件的含义

案例1：ions.mdp#

; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme	= Verlet    ; Buffered neighbor searching 
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = cutoff    ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions
txt

mdp文件以;符号为注释，

这个ions.mdp文件是用于生成离子添加过程的输入参数文件。让我详细解释每个参数的含义：

基本最小化参数#

integrator  = steep         ; 使用最陡下降法进行能量最小化
emtol       = 1000.0        ; 当最大力 < 1000.0 kJ/mol/nm时停止最小化
emstep      = 0.01          ; 最小化步长 (nm)
nsteps      = 50000         ; 最大最小化步数
bash

说明：这里使用了相对宽松的收敛标准（1000 kJ/mol/nm），因为目的只是生成一个合理的结构来添加离子，而不是完全的能量最小化。

邻居列表和相互作用参数#

nstlist         = 1         ; 每1步更新邻居列表
cutoff-scheme   = Verlet    ; 使用Verlet缓冲邻居搜索
ns_type         = grid      ; 使用网格方法搜索邻居原子
bash

非键相互作用参数#

coulombtype = cutoff        ; 使用截断方法处理静电相互作用
rcoulomb    = 1.0           ; 静电相互作用的截断半径 (1.0 nm)
rvdw        = 1.0           ; 范德华相互作用的截断半径 (1.0 nm)
bash

周期性边界条件#

pbc = xyz                   ; 在x、y、z三个方向都使用周期性边界条件
bash

这个配置文件的特殊用途#

这个ions.mdp文件不是用于真正的能量最小化，而是为了：

  1. 生成.tpr文件：grompp需要这个文件来创建GROMACS的输入文件
  2. 定义力场参数：为genion提供必要的力场和模拟参数
  3. 设置基本环境：定义截断半径、周期性边界条件等

案例2：minim.mdp#

这个minim.mdp文件是用于能量最小化的输入参数文件，用于执行体系的能量最小化，消除原子间的冲突和不利接触。

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions
txt

最小化算法参数#

integrator = steep     # 使用最陡下降法 - 简单稳定，适合初始最小化
emtol = 1000.0        # 收敛标准：最大力 < 1000 kJ/mol/nm
emstep = 0.01         # 步长0.01 nm - 保守值，保证稳定性
nsteps = 50000        # 最大步数 - 提供充足的优化空间
bash

邻居搜索设置#

nstlist = 1           # 每步都更新邻居列表 - 确保准确性
cutoff-scheme = Verlet # 现代Verlet缓冲搜索算法
ns_type = grid        # 网格搜索方法 - 计算效率高
bash

相互作用处理#

coulombtype = PME     # 使用粒子网格Ewald方法 - 精确处理长程静电
rcoulomb = 1.0        # 静电截断半径1.0 nm
rvdw = 1.0           # 范德华截断半径1.0 nm
pbc = xyz            # 三维周期性边界条件
bash

原文内容#

http://www.mdtutorials.com/gmx/ ↗

分子动力学模拟：Gromacs教程1-水中的水解酶
https://blog.zxmath.top/blog/gromacs1
Author Zhang Xin
Published at 2025年11月30日
Copyright CC BY-NC-SA 4.0
Buy me a cup of coffee ☕.

分子动力学模拟：Gromacs教程3-蛋白配体复合物

Linux运维：服务器使用简明教程
Comment seems to stuck. Try to refresh?✨