GROMACS 分子动力学

GPU 加速的分子动力学模拟——从蛋白质折叠到药物发现

GROMACS（GROningen 化学模拟软件）是世界上使用最广泛的分子动力学模拟软件包。最初在格罗宁根大学开发，如今由全球社区维护，是计算化学和结构生物学实验室的主力工具。

借助 GPU 加速，GROMACS 可以以在仅使用 CPU 的硬件上需要数周时间才能完成的速度模拟百万级原子系统。Clore.ai 的经济实惠 GPU 租赁使大规模 MD 模拟对个人研究人员和小型实验室变得可及。

---

可以模拟什么？

• 蛋白折叠与动力学 —— 观察纳秒到微秒尺度的构象变化

• 药物-蛋白结合 —— 计算结合自由能以用于药物发现

• 膜模拟 —— 脂质双层、膜蛋白、离子运输

• 蛋白-蛋白相互作用 —— 研究复合物形成与界面动力学

• 材料科学 —— 聚合物、纳米粒子、水模型

• 自由能计算 —— 变换态（alchemical 转换）、PME

---

先决条件

• 具有 GPU 租用的 Clore.ai 帐户

• 基本的 Linux 命令行知识

• 分子体系文件（拓扑 + 坐标），或使用示例体系

• 可选：在本地安装 GROMACS 以便可视化（VMD、Pymol）

---

为何使用 GPU 加速的 GROMACS？

启用 GPU 卸载的 GROMACS 可带来显著提速：

系统规模	仅 CPU（ns/天）	单块 A100（ns/天）	加速比
25K 原子	~50	~800	约 16 倍
100K 原子	~15	~400	约 27 倍
500K 原子	~3	~150	约 50 倍
1M 原子	~1	~80	约 80 倍

GPU 加速对大型系统（>100K 原子）最为有益。 对于小型测试系统，由于数据传输开销，CPU 性能可能具有可比性。

---

步骤 1 — 在 Clore.ai 上租用 GPU

1. 前往 clore.ai → 市场

2. 按 GPU 筛选： A100、RTX 4090 或 RTX 3090 推荐

3. 对于大型系统（>500K 原子）：选择 A100 40GB 或 80GB

4. 对于标准模拟：RTX 4090 或 RTX 3090 性价比极佳

推荐配置：

• GPU：A100 40GB 或 RTX 4090

• CPU：16+ 核（GROMACS 在非键合相互作用上使用多核）

• 内存：32GB+

• 磁盘：50GB+（轨迹文件可能很大）

---

步骤 2 — 部署 GROMACS 容器

使用 NVIDIA 的官方 HPC GROMACS 镜像——针对支持 CUDA 的 NVIDIA GPU 做了优化：

Docker 镜像：

nvcr.io/hpc/gromacs:2023.2

暴露端口：

22

环境变量：

NVIDIA_VISIBLE_DEVICES=all
NVIDIA_DRIVER_CAPABILITIES=compute,utility
GMX_GPU_DD_COMMS=true
GMX_GPU_PME_PP_COMMS=true
GMX_FORCE_UPDATE_DEFAULT_GPU=true

GROMACS 的 NVIDIA 环境变量：

• GMX_GPU_DD_COMMS=true —— 启用基于 GPU 的域分解通信

• GMX_GPU_PME_PP_COMMS=true —— 启用基于 GPU 的 PME-PP 通信

• GMX_FORCE_UPDATE_DEFAULT_GPU=true —— 强制 GPU 坐标更新（显著加速）

---

第 3 步 — 连接并验证

ssh root@<server-ip> -p <ssh-port>

# 检查 GROMACS 版本
gmx --version

# 检查 GPU 是否可用
nvidia-smi

# 验证 GROMACS 能否识别 GPU
gmx mdrun -h 2>&1 | grep -i gpu

预期输出 gmx --version 应显示：

GROMACS 版本：2023.2
CUDA 版本：11.x 或 12.x
GPU 支持：CUDA

---

步骤 4 — 准备你的体系

使用示例体系（溶液中的溶菌酶）

这是经典的 GROMACS 教程体系——非常适合测试你的环境：

# 创建工作目录
mkdir -p /workspace/lysozyme && cd /workspace/lysozyme

# 下载溶菌酶 PDB 结构
wget https://files.rcsb.org/download/1AKI.pdb -O 1AKI.pdb

# 从晶体结构中去除水分子
grep -v HOH 1AKI.pdb > 1AKI_clean.pdb

# 使用 AMBER99SB 力场生成拓扑
gmx pdb2gmx \
    -f 1AKI_clean.pdb \
    -o processed.gro \
    -water spce \
    -ff amber99sb-ildn

在提示选择力场时，选择 amber99sb-ildn （通常是选项 6）。

---

步骤 5 — 构建模拟盒

# 定义模拟盒（十二面体，距蛋白 1.0 nm）
gmx editconf \
    -f processed.gro \
    -o boxed.gro \
    -c \
    -d 1.0 \
    -bt dodecahedron

# 用水溶剂化盒子
gmx solvate \
    -cp boxed.gro \
    -cs spc216.gro \
    -o solvated.gro \
    -p topol.top

# 添加离子以中和电荷
# 首先为添加离子创建 tpr 文件
gmx grompp \
    -f /usr/local/gromacs/share/gromacs/top/em.mdp \
    -c solvated.gro \
    -p topol.top \
    -o ions.tpr

# 添加 Na+ 和 Cl- 离子（0.15M NaCl）
gmx genion \
    -s ions.tpr \
    -o ionized.gro \
    -p topol.top \
    -pname NA \
    -nname CL \
    -neutral \
    -conc 0.15
# 在提示时选择组 13（SOL）

---

步骤 6 — 能量最小化

# 创建能量最小化 MDP 文件
cat > em.mdp << 'EOF'
; 能量最小化参数
integrator      = steep         ; 最速下降法最小化
emtol           = 1000.0        ; 当最大力 < 1000 kJ/mol/nm 时停止
emstep          = 0.01          ; 初始步长
nsteps          = 50000         ; 最大最小化步数
nstlist         = 1
cutoff-scheme   = Verlet
ns_type         = grid
coulombtype     = PME
rcoulomb        = 1.0
rvdw            = 1.0
pbc             = xyz
EOF

# 为能量最小化准备 TPR
gmx grompp \
    -f em.mdp \
    -c ionized.gro \
    -p topol.top \
    -o em.tpr

# 在 GPU 上运行能量最小化
gmx mdrun \
    -v \
    -deffnm em \
    -gpu_id 0 \
    -ntmpi 1 \
    -ntomp 8

# 检查能量是否收敛
gmx energy -f em.edr -o em_potential.xvg
# 选择 10（Potential），然后选择 0 退出

---

步骤 7 — NVT 平衡（温度）

cat > nvt.mdp << 'EOF'
; NVT 平衡
define              = -DPOSRES      ; 位置约束
integrator          = md            ; leap-frog 积分器
nsteps              = 50000         ; 100 ps（2 fs 步长）
dt                  = 0.002         ; 2 fs 步长
nstxout             = 500
nstvout             = 500
nstenergy           = 500
nstlog              = 500
continuation        = no
constraint_algorithm = lincs
constraints         = h-bonds
lincs_iter          = 1
lincs_order         = 4
cutoff-scheme       = Verlet
ns_type             = grid
nstlist             = 10
rcoulomb            = 1.0
rvdw                = 1.0
DispCorr            = EnerPres
coulombtype         = PME
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
EOF

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

gmx mdrun \
    -deffnm nvt \
    -gpu_id 0 \
    -ntmpi 1 \
    -ntomp 8 \
    -nb gpu \
    -bonded gpu \
    -pme gpu \
    -update gpu

---

步骤 8 — NPT 平衡（压力）

cat > npt.mdp << 'EOF'
; NPT 平衡
define              = -DPOSRES
integrator          = md
nsteps              = 50000
dt                  = 0.002
nstxout             = 500
nstvout             = 500
nstenergy           = 500
nstlog              = 500
continuation        = yes
constraint_algorithm = lincs
constraints         = h-bonds
cutoff-scheme       = Verlet
ns_type             = grid
nstlist             = 10
rcoulomb            = 1.0
rvdw                = 1.0
DispCorr            = EnerPres
coulombtype         = PME
tcoupl              = V-rescale
tc-grps             = Protein Non-Protein
tau_t               = 0.1 0.1
ref_t               = 300 300
pcoupl              = Parrinello-Rahman
pcoupltype          = isotropic
tau_p               = 2.0
ref_p               = 1.0
compressibility     = 4.5e-5
refcoord_scaling    = com
pbc                 = xyz
EOF

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

gmx mdrun \
    -deffnm npt \
    -gpu_id 0 \
    -ntmpi 1 \
    -ntomp 8 \
    -nb gpu \
    -bonded gpu \
    -pme gpu \
    -update gpu

---

步骤 9 — 生产阶段 MD 运行

cat > md.mdp << 'EOF'
; 生产阶段 MD 运行
integrator          = md
nsteps              = 5000000      ; 10 ns（2 fs 步长）
dt                  = 0.002
nstxout-compressed  = 5000        ; 每 10 ps 保存坐标
nstenergy           = 5000
nstlog              = 5000
continuation        = yes
constraint_algorithm = lincs
constraints         = h-bonds
cutoff-scheme       = Verlet
ns_type             = grid
nstlist             = 10
rcoulomb            = 1.0
rvdw                = 1.0
DispCorr            = EnerPres
coulombtype         = PME
tcoupl              = V-rescale
tc-grps             = Protein Non-Protein
tau_t               = 0.1 0.1
ref_t               = 300 300
pcoupl              = Parrinello-Rahman
pcoupltype          = isotropic
tau_p               = 2.0
ref_p               = 1.0
compressibility     = 4.5e-5
pbc                 = xyz
EOF

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

# 全 GPU 卸载的生产运行
gmx mdrun \
    -deffnm md \
    -gpu_id 0 \
    -ntmpi 1 \
    -ntomp 16 \
    -nb gpu \
    -bonded gpu \
    -pme gpu \
    -update gpu \
    -v

实时监控进度：

tail -f md.log | grep -E "(ns/day|Step|Time)"

---

步骤 10 — 分析

基本轨迹分析

# RMSD（随着时间的主链稳定性）
gmx rms \
    -s md.tpr \
    -f md.xtc \
    -o rmsd.xvg \
    -tu ns
# 对参考组和拟合组都选择 4（Backbone）

# RMSF（按残基的柔性）
gmx rmsf \
    -s md.tpr \
    -f md.xtc \
    -o rmsf.xvg \
    -res
# 选择 4（Backbone）

# 回转半径（紧凑性）
gmx gyrate \
    -s md.tpr \
    -f md.xtc \
    -o gyrate.xvg
# 选择 1（Protein）

# 氢键
gmx hbond \
    -s md.tpr \
    -f md.xtc \
    -num hbonds.xvg
# 对供体和受体都选择 1（Protein）

绘制 XVG 文件

import numpy as np
import matplotlib.pyplot as plt

# 加载 RMSD 数据
data = np.loadtxt('rmsd.xvg', comments=['@', '#'])
time = data[:, 0]      # 单位为 ns
rmsd = data[:, 1] * 10 # 将 nm 转换为埃

plt.figure(figsize=(10, 4))
plt.plot(time, rmsd)
plt.xlabel('Time (ns)')
plt.ylabel('RMSD (Å)')
plt.title('Backbone RMSD')
plt.grid(True, alpha=0.3)
plt.savefig('rmsd_plot.png', dpi=150, bbox_inches='tight')

传输结果

# 从本地机器：
rsync -avz -e "ssh -p <ssh-port>" \
    root@<server-ip>:/workspace/lysozyme/ \
    ./md_results/

---

多 GPU 模拟

对于非常大的体系，使用多块 GPU 配合域分解：

# 4-GPU 运行（调整 -ntmpi 以匹配 GPU 数量）
gmx mdrun \
    -deffnm md \
    -gpu_id 0123 \
    -ntmpi 4 \
    -ntomp 4 \
    -nb gpu \
    -bonded gpu \
    -pme gpu \
    -npme 1 \
    -update gpu \
    -v

多 GPU 效率： 超过 4 块 GPU 的扩展通常仅对 >1 百万原子系统有益。对于较小系统，在 Clore.ai 上使用单块高端 GPU 更具成本效益。

---

故障排除

致命错误：没有 GPU 卸载

# 检查 CUDA 是否工作正常
nvidia-smi
python3 -c "import ctypes; ctypes.CDLL('libcuda.so')"

# 强制回退到 CPU 以便测试
gmx mdrun -deffnm md -ntmpi 1 -ntomp 16

体系爆炸 / 负体积

这通常表示能量最小化存在问题：

# 用更小的步长进行更长时间的最小化
# 编辑 em.mdp：emstep = 0.001, emtol = 100

性能缓慢

# 在运行期间检查 GPU 利用率
watch -n 1 nvidia-smi

# 调整 GPU 卸载设置
gmx mdrun -deffnm md -nb gpu -pme gpu -bonded gpu -update gpu \
    -ntmpi 1 -ntomp $(nproc)

---

常用力场

力场	适合用于
amber99sb-ildn	蛋白质，通用用途
charmm36m	蛋白质 + 脂质膜
gromos54a7	类药物分子
oplsaa	有机分子，脂类

---

成本估算

模拟	系统规模	GPU	时间	成本
10 ns 蛋白质	25K 原子	RTX 3090	约 2 小时	~$0.60
100 ns 蛋白质	25K 原子	A100 40G	约 6 小时	~$4.50
100 ns 膜体系	200K 原子	A100 80G	约 12 小时	~$9
1 μs 蛋白质	25K 原子	A100 80G	约 3 天	~$55

---

附加资源

• GROMACS 文档

• GROMACS 教程（Justin Lemkul）

• NVIDIA HPC GROMACS 容器

• GROMACS GitHub

• Amber 力场参数

• CHARMM-GUI 膜构建器

---

在 Clore.ai 上运行 GROMACS 使研究人员能够以远低于 AWS 或 Azure 的价格使用 A100 和 RTX 4090 GPU——使长期 MD 模拟对学术实验室和个人研究人员在经济上可行。

---

Clore.ai 的 GPU 建议

在 Clore.ai 上的预估费用	开发/测试	RTX 3090（24GB）
~$0.12/每 GPU/每小时	生产	RTX 4090（24GB）
标准 MD 模拟	大规模	A100 80GB
大体系 / 长时程运行	💡 本指南中的所有示例均可部署在	Clore.ai

GPU 服务器上。浏览可用 GPU 并按小时租用 — 无需承诺，提供完整的 root 访问权限。 Clore.ai GPU 服务器。浏览可用 GPU 并按小时租用 — 无需承诺，提供完整的 root 访问权限。