Молекулярная динамика GROMACS

Молекулярно-динамические симуляции с ускорением на GPU — от сворачивания белков до поиска лекарств

GROMACS (GROningen MAchine for Chemical Simulations) — самый широко используемый в мире пакет для молекулярно-динамического моделирования. Изначально разработанный в Университете Гронингена, сейчас он поддерживается глобальным сообществом и является рабочей лошадкой вычислительной химии и структурной биологии в лабораториях по всему миру.

Благодаря ускорению на GPU GROMACS может моделировать системы из миллионов атомов со скоростью, для которой на оборудовании только с CPU потребовались бы недели. Доступная аренда GPU от Clore.ai делает крупномасштабные MD-симуляции доступными для отдельных исследователей и небольших лабораторий.

---

Что можно моделировать?

• Сворачивание и динамика белков — наблюдайте конформационные изменения за наносекунды–микросекунды

• Связывание лекарств с белками — рассчитывайте свободные энергии связывания для поиска лекарств

• Симуляции мембран — липидные бислои, мембранные белки, транспорт ионов

• Взаимодействия белок-белок — изучайте образование комплексов и динамику интерфейсов

• Наука о материалах — полимеры, наночастицы, модели воды

• Расчёты свободной энергии — алхимические преобразования, PME

---

Предварительные требования

• Аккаунт Clore.ai с арендой GPU

• Базовые знания командной строки Linux

• Файлы молекулярной системы (топология + координаты) или используйте примеры систем

• Необязательно: GROMACS локально для визуализации (VMD, Pymol)

---

Зачем использовать GROMACS с ускорением на GPU?

GROMACS с выгрузкой вычислений на GPU обеспечивает значительное ускорение:

Размер системы	Только CPU (нс/день)	Одна A100 (нс/день)	Ускорение
25 тыс. атомов	~50	~800	~16x
100 тыс. атомов	~15	~400	~27x
500 тыс. атомов	~3	~150	~50x
1 млн атомов	~1	~80	~80x

Ускорение на GPU наиболее полезно для крупных систем (>100K атомов). Для небольших тестовых систем производительность CPU может быть сопоставимой из-за накладных расходов на передачу данных.

---

Шаг 1 — Арендуйте GPU на Clore.ai

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 использует многоядерность для невалентных взаимодействий)

• RAM: 32 ГБ+

• Диск: 50 ГБ+ (траектории могут быть большими)

---

Шаг 2 — разверните контейнер GROMACS

Используйте официальный HPC-образ GROMACS от NVIDIA — он оптимизирован для GPU NVIDIA с поддержкой CUDA:

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

Переменные окружения NVIDIA для GROMACS:

• GMX_GPU_DD_COMMS=true — включает коммуникации разложения по доменам на GPU

• GMX_GPU_PME_PP_COMMS=true — включает коммуникации PME-PP на GPU

• 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 нм от белка)
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.15 М 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 кДж/моль/нм
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 пс (шаг 2 фс)
dt                  = 0.002         ; Шаг 2 фс
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 нс (шаг 2 фс)
dt                  = 0.002
nstxout-compressed  = 5000        ; Сохранять координаты каждые 10 пс
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]      # в нс
rmsd = data[:, 1] * 10 # преобразовать нм в ангстремы

plt.figure(figsize=(10, 4))
plt.plot(time, rmsd)
plt.xlabel('Время (нс)')
plt.ylabel('RMSD (Å)')
plt.title('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 миллиона атомов. Для меньших систем одна высокопроизводительная GPU экономически выгоднее на Clore.ai.

---

Устранение неполадок

Фатальная ошибка: нет выгрузки на 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 нс	25 тыс. атомов	RTX 3090	~2h	~$0.60
Белок 100 нс	25 тыс. атомов	A100 40G	~6 ч	~$4.50
Мембрана 100 нс	200 тыс. атомов	A100 80G	~12 ч	~$9
Белок 1 мкс	25 тыс. атомов	A100 80G	~3 дня	~$55

---

Дополнительные ресурсы

• Документация GROMACS

• Учебные материалы по GROMACS (Джастин Лемкул)

• Контейнер NVIDIA HPC GROMACS

• GROMACS на GitHub

• Параметры силового поля Amber

• Конструктор мембран CHARMM-GUI

---

Запуск GROMACS на Clore.ai позволяет исследователям получать доступ к GPU A100 и RTX 4090 за долю стоимости AWS или Azure, делая длительные МД-симуляции экономически выгодными для академических лабораторий и отдельных исследователей.

---

Рекомендации по GPU для Clore.ai

Сценарий использования	Рекомендуемый GPU	Оценочная стоимость на Clore.ai
Разработка/тестирование	RTX 3090 (24 ГБ)	$0.07–0.21/гпу/ч
Стандартные МД-симуляции	RTX 4090 (24 ГБ)	$0.14–0.42/гпу/ч
Крупные системы / длительные запуски	A100 80 ГБ	голое железо

💡 Все примеры в этом руководстве можно развернуть на Clore.ai GPU-серверах. Просматривайте доступные GPU и арендуйте по часам — без обязательств, с полным root-доступом.