GROMACS — 高性能分子动力学模拟软件¶

一句话说明¶

GROMACS 是用于蛋白质、脂质膜、DNA 等生物大分子的分子动力学（MD）模拟软件，通过模拟原子在力场作用下随时间的运动来研究蛋白质折叠、构象变化、配体结合等生物过程。白话理解：用物理方程让蛋白质"动起来"，看它怎么运动。

安装与配置¶

# 方法1：conda 安装（推荐，含 GPU 版本）
conda install -c bioconda -c conda-forge gromacs  # CPU 版本
# GPU 版本（需要 CUDA）
conda install -c conda-forge gromacs cudatoolkit=11.8  # CUDA 11.8

# 确认版本（当前 v2024.x）
gmx --version  # 查看 GROMACS 版本信息

# 方法2：从源码编译（针对特定 GPU 优化，性能最好）
cmake .. -DGMX_BUILD_OWN_FFTW=ON \   # 自带 FFTW 库
         -DGMX_GPU=CUDA \             # 启用 CUDA GPU 加速
         -DGMX_MPI=off \              # 不用 MPI（单机模式）
         -DCMAKE_INSTALL_PREFIX=$HOME/gromacs2024  # 安装目录
make -j 16  # 16 线程编译
make install

# 测试安装
gmx grompp --help  # 测试 grompp 命令

核心用法¶

标准 MD 模拟流程（5步）¶

# === 步骤 1：准备拓扑文件 ===
# 从 PDB 文件生成 GROMACS 拓扑
gmx pdb2gmx \
  -f protein.pdb \        # 输入 PDB 文件
  -o processed.gro \      # 输出结构文件（GROMACS 格式）
  -p topol.top \          # 输出拓扑文件
  -water spce \           # 水模型：SPC/E（常用）
  -ff amber99sb-ildn      # 力场：AMBER99SB-ILDN（蛋白质常用）

# === 步骤 2：构建溶剂盒子 ===
# 在蛋白质外面加水盒子
gmx editconf \
  -f processed.gro \      # 输入结构
  -o box.gro \            # 输出加盒子后的结构
  -c \                    # 居中蛋白质
  -d 1.0 \                # 蛋白质到盒子边界距离 1.0 nm
  -bt cubic               # 立方体盒子

gmx solvate \
  -cp box.gro \           # 含盒子的结构
  -cs spc216.gro \        # 水分子库（GROMACS 自带）
  -o solvated.gro \       # 加水后的结构
  -p topol.top            # 拓扑文件（自动更新水分子数）

# 加离子（中和电荷）
gmx grompp -f ions.mdp -c solvated.gro -p topol.top -o ions.tpr  # 生成运行文件
gmx genion -s ions.tpr -o ionized.gro -p topol.top -pname NA -nname CL -neutral  # 加 Na+/Cl-

# === 步骤 3：能量最小化 ===
gmx grompp \
  -f em.mdp \             # 能量最小化参数文件
  -c ionized.gro \        # 输入结构
  -p topol.top \          # 拓扑文件
  -o em.tpr               # 输出运行文件

gmx mdrun \
  -v \                    # 显示进度
  -deffnm em              # 所有输出文件前缀为 em（em.edr/em.trr/em.log）

# === 步骤 4：平衡（NVT + NPT）===
# NVT 平衡（固定体积，平衡温度，100ps）
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt -ntmpi 1 -ntomp 8 -gpu_id 0  # 使用 GPU 0，8 个 CPU 线程

# NPT 平衡（固定压力，平衡密度，100ps）
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt -ntmpi 1 -ntomp 8 -gpu_id 0

# === 步骤 5：生产性 MD 模拟 ===
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 8 -gpu_id 0  # 正式模拟

参数详解¶

mdp 文件关键参数（md.mdp）¶

; 模拟时间控制
integrator  = md         ; 积分算法：md=蛙跳算法（最常用）
nsteps      = 5000000    ; 模拟步数（×dt = 总时间：5000000×0.002ps = 10ns）
dt          = 0.002      ; 时间步长：0.002ps（2飞秒，常用值）

; 输出控制
nstxout-compressed = 5000   ; 每5000步输出一帧压缩轨迹（xtc格式）
nstenergy = 1000            ; 每1000步输出一次能量

; 温度控制
tcoupl      = V-rescale    ; 恒温算法：V-rescale（推荐）
tc-grps     = Protein Non-Protein  ; 分组控温
tau_t       = 0.1 0.1      ; 耦合时间常数（ps）
ref_t       = 310 310      ; 目标温度：310K（37°C，体温）

; 压力控制
pcoupl      = Parrinello-Rahman   ; 恒压算法（NPT 模拟时用）
pcoupltype  = isotropic           ; 各向同性压力
tau_p       = 2.0                 ; 压力耦合时间常数（ps）
ref_p       = 1.0                 ; 目标压力：1bar

; 非键相互作用
cutoff-scheme = Verlet     ; 截断方案（推荐）
coulombtype   = PME        ; 静电：粒子网格 Ewald（长程精确）
rcoulomb      = 1.2        ; 库仑截断距离：1.2nm
rvdw          = 1.2        ; 范德华截断距离：1.2nm

实战案例¶

案例：分析 MD 轨迹（RMSD + RMSF）¶

# RMSD：结构与初始状态的均方根偏差（评估稳定性）
gmx rms \
  -s md.tpr \          # 参考结构（初始拓扑）
  -f md.xtc \          # 轨迹文件
  -o rmsd.xvg \        # 输出 xvg 格式数据
  -tu ns \             # 时间单位：纳秒
  -select "Backbone"   # 只分析骨架原子（C-Cα-N）
# 选择 Group: 1 (Backbone) 回车两次

# RMSF：各残基的均方根波动（评估灵活性）
gmx rmsf \
  -s md.tpr \
  -f md.xtc \
  -o rmsf.xvg \
  -res \               # 按残基计算
  -select "Backbone"   # 分析骨架

# 用 Python 绘制 RMSD/RMSF 曲线
python3 << 'EOF'
import numpy as np
import matplotlib.pyplot as plt

# 读取 xvg 格式（跳过注释行）
data = np.loadtxt("rmsd.xvg",       # 读取 RMSD 数据
                  comments=["#", "@"])  # 跳过注释行
time = data[:, 0]                    # 第一列：时间（ns）
rmsd = data[:, 1] * 10               # 第二列：RMSD（nm→Å，×10）

plt.figure(figsize=(10, 4))
plt.plot(time, rmsd, color="steelblue", linewidth=1)
plt.xlabel("Time (ns)")
plt.ylabel("RMSD (Å)")
plt.title("Backbone RMSD")
plt.tight_layout()
plt.savefig("rmsd_plot.pdf")
print(f"平均 RMSD: {rmsd.mean():.2f} Å")  # 打印平均值
EOF

常见报错与解决¶

报错	原因	解决方法
`Fatal error: Atom not found`	PDB 有缺失原子/非标准残基	用 `pdb4amber` 或手动修复 PDB
`WARNING: clashes`	加水后原子重叠	能量最小化步数不够；减小 `-d` 参数
`LINCS warning`	时间步长太大	把 `dt` 从 0.002 降到 0.001 ps
GPU 内存不足	体系太大	增大 PME 网格分工；用多 GPU
模拟温度异常升高	初始结构有原子重叠	延长能量最小化步数

速查表¶

# 安装
conda install -c conda-forge gromacs

# 标准流程命令
gmx pdb2gmx   # PDB → 拓扑
gmx editconf  # 构建盒子
gmx solvate   # 加水分子
gmx genion    # 加离子
gmx grompp    # 生成运行文件（.tpr）
gmx mdrun     # 执行模拟

# 轨迹分析工具
gmx rms       # RMSD
gmx rmsf      # RMSF（残基灵活性）
gmx gyrate    # 旋转半径（蛋白质紧密程度）
gmx hbond     # 氢键分析
gmx energy    # 提取能量项
gmx trjconv   # 轨迹格式转换/中心化/PBC修正

# GPU 加速运行（单机单GPU）
gmx mdrun -ntmpi 1 -ntomp 8 -gpu_id 0

# 当前版本：GROMACS 2024.x