GROMACS Molecular Dynamics

GPU-accelerated molecular dynamics simulations — from protein folding to drug discovery

GROMACS (GROningen MAchine for Chemical Simulations) is the world's most widely used molecular dynamics simulation package. Originally developed at the University of Groningen, it's now maintained by a global community and is the workhorse of computational chemistry and structural biology labs worldwide.

With GPU acceleration, GROMACS can simulate systems of millions of atoms at speeds that would take weeks on CPU-only hardware. Clore.ai's affordable GPU rentals make large-scale MD simulations accessible to individual researchers and small labs.

---

What Can You Simulate?

• Protein folding and dynamics — observe conformational changes in nanoseconds to microseconds

• Drug-protein binding — calculate binding free energies for drug discovery

• Membrane simulations — lipid bilayers, membrane proteins, ion transport

• Protein-protein interactions — study complex formation and interface dynamics

• Materials science — polymers, nanoparticles, water models

• Free energy calculations — alchemical transformations, PME

---

Prerequisites

• Clore.ai account with GPU rental

• Basic Linux command line knowledge

• Molecular system files (topology + coordinates), or use example systems

• Optional: GROMACS locally for visualization (VMD, Pymol)

---

Why Use GPU-Accelerated GROMACS?

GROMACS with GPU offloading provides dramatic speedups:

System Size	CPU Only (ns/day)	Single A100 (ns/day)	Speedup
25K atoms	~50	~800	~16x
100K atoms	~15	~400	~27x
500K atoms	~3	~150	~50x
1M atoms	~1	~80	~80x

GPU acceleration is most beneficial for large systems (>100K atoms). For small test systems, CPU performance may be comparable due to data transfer overhead.

---

Step 1 — Rent a GPU on Clore.ai

1. Go to clore.ai → Marketplace

2. Filter by GPU: A100, RTX 4090, or RTX 3090 recommended

3. For large systems (>500K atoms): choose A100 40GB or 80GB

4. For standard simulations: RTX 4090 or RTX 3090 are excellent value

Recommended specs:

• GPU: A100 40GB or RTX 4090

• CPU: 16+ cores (GROMACS uses multi-core for non-bonded interactions)

• RAM: 32GB+

• Disk: 50GB+ (trajectories can be large)

---

Step 2 — Deploy GROMACS Container

Use NVIDIA's official HPC GROMACS image — it's optimized for NVIDIA GPUs with CUDA support:

Docker Image:

nvcr.io/hpc/gromacs:2023.2

Exposed Ports:

22

Environment Variables:

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 Environment Variables for GROMACS:

• GMX_GPU_DD_COMMS=true — enables GPU-based domain decomposition communications

• GMX_GPU_PME_PP_COMMS=true — enables GPU-based PME-PP communications

• GMX_FORCE_UPDATE_DEFAULT_GPU=true — forces GPU coordinate update (significant speedup)

---

Step 3 — Connect and Verify

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

# Check GROMACS version
gmx --version

# Check GPU availability
nvidia-smi

# Verify GROMACS can see GPU
gmx mdrun -h 2>&1 | grep -i gpu

Expected output from gmx --version should show:

GROMACS version: 2023.2
CUDA version: 11.x or 12.x
GPU support: CUDA

---

Step 4 — Prepare Your System

Using an Example System (Lysozyme in Water)

This is the classic GROMACS tutorial system — perfect for testing your setup:

# Create working directory
mkdir -p /workspace/lysozyme && cd /workspace/lysozyme

# Download lysozyme PDB structure
wget https://files.rcsb.org/download/1AKI.pdb -O 1AKI.pdb

# Remove water molecules from crystal structure
grep -v HOH 1AKI.pdb > 1AKI_clean.pdb

# Generate topology using AMBER99SB force field
gmx pdb2gmx \
    -f 1AKI_clean.pdb \
    -o processed.gro \
    -water spce \
    -ff amber99sb-ildn

When prompted for force field selection, choose amber99sb-ildn (option 6 typically).

---

Step 5 — Build the Simulation Box

# Define simulation box (dodecahedron, 1.0 nm from protein)
gmx editconf \
    -f processed.gro \
    -o boxed.gro \
    -c \
    -d 1.0 \
    -bt dodecahedron

# Solvate the box with water
gmx solvate \
    -cp boxed.gro \
    -cs spc216.gro \
    -o solvated.gro \
    -p topol.top

# Add ions to neutralize charge
# First create tpr file for ion addition
gmx grompp \
    -f /usr/local/gromacs/share/gromacs/top/em.mdp \
    -c solvated.gro \
    -p topol.top \
    -o ions.tpr

# Add Na+ and Cl- ions (0.15M NaCl)
gmx genion \
    -s ions.tpr \
    -o ionized.gro \
    -p topol.top \
    -pname NA \
    -nname CL \
    -neutral \
    -conc 0.15
# Select group 13 (SOL) when prompted

---

Step 6 — Energy Minimization

# Create energy minimization MDP file
cat > em.mdp << 'EOF'
; Energy minimization parameters
integrator      = steep         ; Steepest descent minimization
emtol           = 1000.0        ; Stop when max force < 1000 kJ/mol/nm
emstep          = 0.01          ; Initial step size
nsteps          = 50000         ; Maximum minimization steps
nstlist         = 1
cutoff-scheme   = Verlet
ns_type         = grid
coulombtype     = PME
rcoulomb        = 1.0
rvdw            = 1.0
pbc             = xyz
EOF

# Prepare TPR for energy minimization
gmx grompp \
    -f em.mdp \
    -c ionized.gro \
    -p topol.top \
    -o em.tpr

# Run energy minimization on GPU
gmx mdrun \
    -v \
    -deffnm em \
    -gpu_id 0 \
    -ntmpi 1 \
    -ntomp 8

# Check that energy converged
gmx energy -f em.edr -o em_potential.xvg
# Select 10 (Potential) then 0 to exit

---

Step 7 — NVT Equilibration (Temperature)

cat > nvt.mdp << 'EOF'
; NVT Equilibration
define              = -DPOSRES      ; Position restraints
integrator          = md            ; Leap-frog integrator
nsteps              = 50000         ; 100 ps (2 fs timestep)
dt                  = 0.002         ; 2 fs timestep
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

---

Step 8 — NPT Equilibration (Pressure)

cat > npt.mdp << 'EOF'
; NPT Equilibration
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

---

Step 9 — Production MD Run

cat > md.mdp << 'EOF'
; Production MD Run
integrator          = md
nsteps              = 5000000      ; 10 ns (2 fs timestep)
dt                  = 0.002
nstxout-compressed  = 5000        ; Save coordinates every 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

# Full GPU offload production run
gmx mdrun \
    -deffnm md \
    -gpu_id 0 \
    -ntmpi 1 \
    -ntomp 16 \
    -nb gpu \
    -bonded gpu \
    -pme gpu \
    -update gpu \
    -v

Monitor progress in real time:

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

---

Step 10 — Analysis

Basic Trajectory Analysis

# RMSD (backbone stability over time)
gmx rms \
    -s md.tpr \
    -f md.xtc \
    -o rmsd.xvg \
    -tu ns
# Select 4 (Backbone) for both reference and fitted groups

# RMSF (per-residue flexibility)
gmx rmsf \
    -s md.tpr \
    -f md.xtc \
    -o rmsf.xvg \
    -res
# Select 4 (Backbone)

# Radius of gyration (compactness)
gmx gyrate \
    -s md.tpr \
    -f md.xtc \
    -o gyrate.xvg
# Select 1 (Protein)

# Hydrogen bonds
gmx hbond \
    -s md.tpr \
    -f md.xtc \
    -num hbonds.xvg
# Select 1 (Protein) for both donor and acceptor

Plot XVG Files

import numpy as np
import matplotlib.pyplot as plt

# Load RMSD data
data = np.loadtxt('rmsd.xvg', comments=['@', '#'])
time = data[:, 0]      # in ns
rmsd = data[:, 1] * 10 # convert nm to Angstroms

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

Transfer Results

# From your local machine:
rsync -avz -e "ssh -p <ssh-port>" \
    root@<server-ip>:/workspace/lysozyme/ \
    ./md_results/

---

Multi-GPU Simulations

For very large systems, use multiple GPUs with domain decomposition:

# 4-GPU run (adjust -ntmpi to match GPU count)
gmx mdrun \
    -deffnm md \
    -gpu_id 0123 \
    -ntmpi 4 \
    -ntomp 4 \
    -nb gpu \
    -bonded gpu \
    -pme gpu \
    -npme 1 \
    -update gpu \
    -v

Multi-GPU efficiency: Scaling beyond 4 GPUs is typically only beneficial for systems >1 million atoms. For smaller systems, a single high-end GPU is more cost-effective on Clore.ai.

---

Troubleshooting

Fatal Error: No GPU Offload

# Check CUDA is working
nvidia-smi
python3 -c "import ctypes; ctypes.CDLL('libcuda.so')"

# Force CPU fallback for testing
gmx mdrun -deffnm md -ntmpi 1 -ntomp 16

System Explodes / Negative Volumes

This usually indicates a problem with energy minimization:

# Run longer minimization with smaller step size
# Edit em.mdp: emstep = 0.001, emtol = 100

Slow Performance

# Check GPU utilization during run
watch -n 1 nvidia-smi

# Tune GPU offload settings
gmx mdrun -deffnm md -nb gpu -pme gpu -bonded gpu -update gpu \
    -ntmpi 1 -ntomp $(nproc)

---

Common Force Fields

Force Field	Best For
amber99sb-ildn	Proteins, general use
charmm36m	Proteins + lipid membranes
gromos54a7	Drug-like molecules
oplsaa	Organic molecules, lipids

---

Cost Estimation

Simulation	System Size	GPU	Time	Cost
10 ns protein	25K atoms	RTX 3090	~2h	~$0.60
100 ns protein	25K atoms	A100 40G	~6h	~$4.50
100 ns membrane	200K atoms	A100 80G	~12h	~$9
1 μs protein	25K atoms	A100 80G	~3 days	~$55

---

Additional Resources

• GROMACS Documentation

• GROMACS Tutorials (Justin Lemkul)

• NVIDIA HPC GROMACS Container

• GROMACS GitHub

• Amber Force Field Parameters

• CHARMM-GUI Membrane Builder

---

Running GROMACS on Clore.ai allows researchers to access A100 and RTX 4090 GPUs at a fraction of AWS or Azure pricing — making long MD simulations economically viable for academic labs and individual researchers.

---

Clore.ai GPU Recommendations

Use Case	Recommended GPU	Est. Cost on Clore.ai
Development/Testing	RTX 3090 (24GB)	~$0.12/gpu/hr
Standard MD Simulations	RTX 4090 (24GB)	~$0.70/gpu/hr
Large Systems / Long Runs	A100 80GB	~$1.20/gpu/hr

💡 All examples in this guide can be deployed on Clore.ai GPU servers. Browse available GPUs and rent by the hour — no commitments, full root access.