Dynamique moléculaire GROMACS

Simulations de dynamique moléculaire accélérées par GPU — du repliement des protéines à la découverte de médicaments

GROMACS (GROningen MAchine pour les simulations chimiques) est le logiciel de simulation de dynamique moléculaire le plus utilisé au monde. Initialement développé à l’Université de Groningue, il est désormais maintenu par une communauté mondiale et constitue l’outil de travail des laboratoires de chimie computationnelle et de biologie structurale du monde entier.

Avec l’accélération GPU, GROMACS peut simuler des systèmes de millions d’atomes à des vitesses qui prendraient des semaines sur du matériel uniquement CPU. Les locations de GPU abordables de Clore.ai rendent les simulations MD à grande échelle accessibles aux chercheurs individuels et aux petits laboratoires.

---

Que pouvez-vous simuler ?

• Repliement des protéines et dynamique — observer les changements conformationnels en nanosecondes à microsecondes

• Liaison médicament-protéine — calculer les énergies libres de liaison pour la découverte de médicaments

• Simulations de membranes — bicouches lipidiques, protéines membranaires, transport ionique

• Interactions protéine-protéine — étudier la formation de complexes et la dynamique des interfaces

• Science des matériaux — polymères, nanoparticules, modèles d’eau

• Calculs d’énergie libre — transformations alchimiques, PME

---

Prérequis

• Compte Clore.ai avec location de GPU

• Connaissances de base de la ligne de commande Linux

• Fichiers du système moléculaire (topologie + coordonnées), ou utilisez des systèmes d’exemple

• Optionnel : GROMACS en local pour la visualisation (VMD, PyMOL)

---

Pourquoi utiliser GROMACS accéléré par GPU ?

GROMACS avec déchargement sur GPU offre des accélérations spectaculaires :

Taille du système	CPU seul (ns/jour)	Un seul A100 (ns/jour)	Accélération
25 k atomes	~50	~800	~16x
100 k atomes	~15	~400	~27x
500 k atomes	~3	~150	~50x
1 M atomes	~1	~80	~80x

L’accélération GPU est surtout bénéfique pour les grands systèmes (>100 k atomes). Pour les petits systèmes de test, les performances CPU peuvent être comparables en raison du surcoût de transfert de données.

---

Étape 1 — Louer un GPU sur Clore.ai

1. Accédez à clore.ai → Place de marché

2. Filtrer par GPU : A100, RTX 4090 ou RTX 3090 recommandé

3. Pour les grands systèmes (>500 k atomes) : choisissez A100 40 Go ou 80 Go

4. Pour les simulations standard : les RTX 4090 ou RTX 3090 offrent un excellent rapport qualité-prix

Spécifications recommandées :

• GPU : A100 40 Go ou RTX 4090

• CPU : 16 cœurs ou plus (GROMACS utilise plusieurs cœurs pour les interactions non liantes)

• RAM : 32 Go ou plus

• Disque : 50 Go ou plus (les trajectoires peuvent être volumineuses)

---

Étape 2 — Déployer le conteneur GROMACS

Utilisez l’image HPC GROMACS officielle de NVIDIA — elle est optimisée pour les GPU NVIDIA avec prise en charge CUDA :

Image Docker :

nvcr.io/hpc/gromacs:2023.2

Ports exposés :

22

Variables d'environnement :

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

Variables d’environnement NVIDIA pour GROMACS :

• GMX_GPU_DD_COMMS=true — active les communications de décomposition en domaines basées sur le GPU

• GMX_GPU_PME_PP_COMMS=true — active les communications PME-PP basées sur le GPU

• GMX_FORCE_UPDATE_DEFAULT_GPU=true — force la mise à jour des coordonnées sur GPU (accélération significative)

---

Étape 3 — Se connecter et vérifier

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

# Vérifier la version de GROMACS
gmx --version

# Vérifier la disponibilité du GPU
nvidia-smi

# Vérifier que GROMACS détecte le GPU
gmx mdrun -h 2>&1 | grep -i gpu

La sortie attendue de gmx --version devrait afficher :

Version de GROMACS : 2023.2
Version de CUDA : 11.x ou 12.x
Prise en charge GPU : CUDA

---

Étape 4 — Préparez votre système

Utilisation d’un système d’exemple (lysozyme dans l’eau)

C’est le système classique du tutoriel GROMACS — parfait pour tester votre installation :

# Créer le répertoire de travail
mkdir -p /workspace/lysozyme && cd /workspace/lysozyme

# Télécharger la structure PDB du lysozyme
wget https://files.rcsb.org/download/1AKI.pdb -O 1AKI.pdb

# Supprimer les molécules d’eau de la structure cristalline
grep -v HOH 1AKI.pdb > 1AKI_clean.pdb

# Générer la topologie à l’aide du champ de force AMBER99SB
gmx pdb2gmx \
    -f 1AKI_clean.pdb \
    -o processed.gro \
    -water spce \
    -ff amber99sb-ildn

Lorsque vous êtes invité à choisir le champ de force, sélectionnez amber99sb-ildn (option 6 en général).

---

Étape 5 — Construire la boîte de simulation

# Définir la boîte de simulation (dodécaèdre, à 1,0 nm de la protéine)
gmx editconf \
    -f processed.gro \
    -o boxed.gro \
    -c \
    -d 1.0 \
    -bt dodecahedron

# Solvater la boîte avec de l’eau
gmx solvate \
    -cp boxed.gro \
    -cs spc216.gro \
    -o solvated.gro \
    -p topol.top

# Ajouter des ions pour neutraliser la charge
# Créer d’abord le fichier tpr pour l’ajout d’ions
gmx grompp \
    -f /usr/local/gromacs/share/gromacs/top/em.mdp \
    -c solvated.gro \
    -p topol.top \
    -o ions.tpr

# Ajouter des ions Na+ et Cl- (NaCl 0,15 M)
gmx genion \
    -s ions.tpr \
    -o ionized.gro \
    -p topol.top \
    -pname NA \
    -nname CL \
    -neutral \
    -conc 0.15
# Sélectionner le groupe 13 (SOL) lorsque vous y êtes invité

---

Étape 6 — Minimisation d’énergie

# Créer le fichier MDP de minimisation d’énergie
cat > em.mdp << 'EOF'
; Paramètres de minimisation d’énergie
integrator      = steep         ; Minimisation par descente la plus raide
emtol           = 1000.0        ; Arrêter lorsque la force maximale < 1000 kJ/mol/nm
emstep          = 0.01          ; Taille de pas initiale
nsteps          = 50000         ; Nombre maximal d’étapes de minimisation
nstlist         = 1
cutoff-scheme   = Verlet
ns_type         = grid
coulombtype     = PME
rcoulomb        = 1.0
rvdw            = 1.0
pbc             = xyz
EOF

# Préparer le TPR pour la minimisation d’énergie
gmx grompp \
    -f em.mdp \
    -c ionized.gro \
    -p topol.top \
    -o em.tpr

# Exécuter la minimisation d’énergie sur GPU
gmx mdrun \
    -v \
    -deffnm em \
    -gpu_id 0 \
    -ntmpi 1 \
    -ntomp 8

# Vérifier que l’énergie a convergé
gmx energy -f em.edr -o em_potential.xvg
# Sélectionner 10 (Potentiel) puis 0 pour quitter

---

Étape 7 — Équilibration NVT (température)

cat > nvt.mdp << 'EOF'
; Équilibration NVT
define              = -DPOSRES      ; Contraintes de position
integrator          = md            ; Intégrateur leap-frog
nsteps              = 50000         ; 100 ps (pas de temps de 2 fs)
dt                  = 0.002         ; Pas de temps de 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             = Protéine Non-protéine
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

---

Étape 8 — Équilibration NPT (pression)

cat > npt.mdp << 'EOF'
; Équilibration 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             = Protéine Non-protéine
tau_t               = 0.1 0.1
ref_t               = 300 300
pcoupl              = Parrinello-Rahman
pcoupltype          = isotrope
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

---

Étape 9 — Exécution MD de production

cat > md.mdp << 'EOF'
; Exécution MD de production
integrator          = md
nsteps              = 5000000      ; 10 ns (pas de temps de 2 fs)
dt                  = 0.002
nstxout-compressed  = 5000        ; Sauvegarder les coordonnées toutes les 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             = Protéine Non-protéine
tau_t               = 0.1 0.1
ref_t               = 300 300
pcoupl              = Parrinello-Rahman
pcoupltype          = isotrope
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

# Exécution de production avec déchargement complet sur GPU
gmx mdrun \
    -deffnm md \
    -gpu_id 0 \
    -ntmpi 1 \
    -ntomp 16 \
    -nb gpu \
    -bonded gpu \
    -pme gpu \
    -update gpu \
    -v

Surveillez la progression en temps réel :

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

---

Étape 10 — Analyse

Analyse de base de la trajectoire

# RMSD (stabilité de la chaîne principale au fil du temps)
gmx rms \
    -s md.tpr \
    -f md.xtc \
    -o rmsd.xvg \
    -tu ns
# Sélectionner 4 (chaîne principale) pour les groupes de référence et ajusté

# RMSF (flexibilité par résidu)
gmx rmsf \
    -s md.tpr \
    -f md.xtc \
    -o rmsf.xvg \
    -res
# Sélectionner 4 (chaîne principale)

# Rayon de giration (compacité)
gmx gyrate \
    -s md.tpr \
    -f md.xtc \
    -o gyrate.xvg
# Sélectionner 1 (Protéine)

# Liaisons hydrogène
gmx hbond \
    -s md.tpr \
    -f md.xtc \
    -num hbonds.xvg
# Sélectionner 1 (Protéine) pour le donneur et l’accepteur

Tracer les fichiers XVG

import numpy as np
import matplotlib.pyplot as plt

# Charger les données RMSD
data = np.loadtxt('rmsd.xvg', comments=['@', '#'])
time = data[:, 0]      # en ns
rmsd = data[:, 1] * 10 # convertir les nm en angstroms

plt.figure(figsize=(10, 4))
plt.plot(time, rmsd)
plt.xlabel('Temps (ns)')
plt.ylabel('RMSD (Å)')
plt.title('RMSD de la chaîne principale')
plt.grid(True, alpha=0.3)
plt.savefig('rmsd_plot.png', dpi=150, bbox_inches='tight')

Transférer les résultats

# Depuis votre machine locale :
rsync -avz -e "ssh -p <ssh-port>" \
    root@<server-ip>:/workspace/lysozyme/ \
    ./md_results/

---

Simulations multi-GPU

Pour les systèmes très grands, utilisez plusieurs GPU avec décomposition en domaines :

# Exécution sur 4 GPU (ajustez -ntmpi pour correspondre au nombre de GPU)
gmx mdrun \
    -deffnm md \
    -gpu_id 0123 \
    -ntmpi 4 \
    -ntomp 4 \
    -nb gpu \
    -bonded gpu \
    -pme gpu \
    -npme 1 \
    -update gpu \
    -v

Efficacité multi-GPU : Le passage à plus de 4 GPU n’est généralement bénéfique que pour les systèmes de plus d’un million d’atomes. Pour les systèmes plus petits, un seul GPU haut de gamme est plus rentable sur Clore.ai.

---

Dépannage

Erreur fatale : aucun déchargement GPU

# Vérifier que CUDA fonctionne
nvidia-smi
python3 -c "import ctypes; ctypes.CDLL('libcuda.so')"

# Forcer le repli sur le CPU pour les tests
gmx mdrun -deffnm md -ntmpi 1 -ntomp 16

Le système explose / Volumes négatifs

Cela indique généralement un problème de minimisation d’énergie :

# Exécuter une minimisation plus longue avec une taille de pas plus petite
# Modifier em.mdp : emstep = 0.001, emtol = 100

Performances lentes

# Vérifier l’utilisation du GPU pendant l’exécution
watch -n 1 nvidia-smi

# Ajuster les paramètres de déchargement GPU
gmx mdrun -deffnm md -nb gpu -pme gpu -bonded gpu -update gpu \
    -ntmpi 1 -ntomp $(nproc)

---

Champs de force courants

Champ de force	Idéal pour
amber99sb-ildn	Protéines, usage général
charmm36m	Protéines + membranes lipidiques
gromos54a7	Molécules de type médicament
oplsaa	Molécules organiques, lipides

---

Estimation des coûts

Simulation	Taille du système	GPU	Temps	Coût
Protéine de 10 ns	25 k atomes	RTX 3090	~2h	~$0.60
Protéine de 100 ns	25 k atomes	A100 40 Go	~6 h	~$4.50
membrane de 100 ns	200 000 atomes	A100 80 Go	~12 h	~$9
protéine de 1 μs	25 k atomes	A100 80 Go	~3 jours	~$55

---

Ressources supplémentaires

• Documentation GROMACS

• Tutoriels GROMACS (Justin Lemkul)

• Conteneur NVIDIA HPC GROMACS

• GitHub GROMACS

• Paramètres du champ de force Amber

• Constructeur de membrane CHARMM-GUI

---

Exécuter GROMACS sur Clore.ai permet aux chercheurs d’accéder aux GPU A100 et RTX 4090 à une fraction du prix d’AWS ou d’Azure — rendant les longues simulations MD économiquement viables pour les laboratoires universitaires et les chercheurs indépendants.

---

Recommandations GPU Clore.ai

Cas d’utilisation	GPU recommandé	Coût estimé sur Clore.ai
Développement/Test	RTX 3090 (24 Go)	0,07–0,21 $/gpu/h
Simulations MD standard	RTX 4090 (24 Go)	0,14–0,42 $/gpu/h
Grands systèmes / longues simulations	A100 80 Go	bare metal

💡 Tous les exemples de ce guide peuvent être déployés sur Clore.ai des serveurs GPU. Parcourez les GPU disponibles et louez à l'heure — sans engagement, accès root complet.