# 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 倍 |

{% hint style="success" %}
**GPU 加速对大型系统（>100K 原子）最为有益。** 对于小型测试系统，由于数据传输开销，CPU 性能可能具有可比性。
{% endhint %}

***

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

1. 前往 [clore.ai](https://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
```

{% hint style="info" %}
**GROMACS 的 NVIDIA 环境变量：**

* `GMX_GPU_DD_COMMS=true` —— 启用基于 GPU 的域分解通信
* `GMX_GPU_PME_PP_COMMS=true` —— 启用基于 GPU 的 PME-PP 通信
* `GMX_FORCE_UPDATE_DEFAULT_GPU=true` —— 强制 GPU 坐标更新（显著加速）
  {% endhint %}

***

## 第 3 步 — 连接并验证

```bash
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 教程体系——非常适合测试你的环境：

```bash
# 创建工作目录
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 — 构建模拟盒

```bash
# 定义模拟盒（十二面体，距蛋白 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 — 能量最小化

```bash
# 创建能量最小化 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 平衡（温度）

```bash
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 平衡（压力）

```bash
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 运行

```bash
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
```

{% hint style="info" %}
**实时监控进度：**

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

{% endhint %}

***

## 步骤 10 — 分析

### 基本轨迹分析

```bash
# 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 文件

```python
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')
```

### 传输结果

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

***

## 多 GPU 模拟

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

```bash
# 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
```

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

***

## 故障排除

### 致命错误：没有 GPU 卸载

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

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

### 体系爆炸 / 负体积

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

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

### 性能缓慢

```bash
# 在运行期间检查 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 文档](https://manual.gromacs.org/)
* [GROMACS 教程（Justin Lemkul）](http://www.mdtutorials.com/gmx/)
* [NVIDIA HPC GROMACS 容器](https://catalog.ngc.nvidia.com/orgs/hpc/containers/gromacs)
* [GROMACS GitHub](https://github.com/gromacs/gromacs)
* [Amber 力场参数](https://ambermd.org/)
* [CHARMM-GUI 膜构建器](https://www.charmm-gui.org/)

***

*在 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](https://clore.ai/marketplace) GPU 服务器。浏览可用 GPU 并按小时租用 — 无需承诺，提供完整的 root 访问权限。