📑 In This Article
Overview
GROMACS (GROningen MAchine for Chemical Simulations) is one of the world's most widely used and highly optimized molecular dynamics (MD) simulation packages. Originally developed for biomolecular simulations, it is now used across chemistry, materials science, and drug discovery. Best of all — it's completely free and open-source. This guide takes you from zero to running your first GROMACS simulation of a protein in water, covering every command you need along the way.
✅ What You'll Be Able to Do
- Prepare and clean a PDB protein structure
- Generate topology with pdb2gmx
- Solvate in TIP3P water + add ions
- Run energy minimization
- Perform NVT + NPT equilibration
- Launch a 10 ns production run
- Analyse RMSD, RMSF & radius of gyration
- Apply GPU acceleration for faster runs
🌊 What is Molecular Dynamics Simulation?
Molecular Dynamics simulates the physical movements of atoms and molecules over time, governed by Newton's equations of motion. Each atom is treated as a particle with mass, position, and velocity. Forces between atoms are calculated from a force field, and positions are updated at each timestep (typically 1–2 femtoseconds).
Timescale
Typical simulations run nanoseconds to microseconds of biological time
Force Fields
AMBER, CHARMM, GROMOS, OPLS-AA describe bonded and non-bonded interactions
Explicit Solvent
Protein is surrounded by explicit water molecules (TIP3P, SPC/E models)
Analysis
RMSD, RMSF, Rg, hydrogen bonds, free energy — extract dynamic properties
🔄 Workflow Overview
💻 Installation
GROMACS can be installed on Linux, macOS, or Windows (WSL). For best performance, compile from source with GPU (CUDA) support.
Ubuntu / Debian
sudo apt-get update
sudo apt-get install gromacs
# Verify installation
gmx --versionCompile from Source (GPU Support)
wget https://ftp.gromacs.org/gromacs/gromacs-2024.tar.gz
tar xvzf gromacs-2024.tar.gz
cd gromacs-2024
mkdir build && cd build
cmake .. -DGMX_BUILD_OWN_FFTW=ON -DGMX_GPU=CUDA
make -j$(nproc)
sudo make install
source /usr/local/gromacs/bin/GMXRC📥 Step 1 — Prepare Your Protein
Download your target protein from RCSB PDB. For this tutorial, we use lysozyme (PDB: 1AKI) — a classic GROMACS tutorial protein.
grep -v "HOH" 1AKI.pdb > 1AKI_clean.pdbOpen your PDB file in PyMOL and check the header for "MISSING" residues. Incomplete structures can cause simulation instability. Use MODELLER or Swiss-Model to rebuild missing loops.
⚙️ Step 2 — Generate Topology with pdb2gmx
The pdb2gmx
tool converts your PDB to GROMACS format and assigns force field parameters.
gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water tip3p
# When prompted, select force field:
# 6: AMBER99SB-ILDN protein, nucleic AMBER94This generates: topol.top (topology),
posre.itp (position restraints), and
1AKI_processed.gro (coordinates).
📦 Step 3 — Define Box and Solvate
Define simulation box
Create a dodecahedral box with 1.0 nm distance from protein to box edge.
Fill with water
Solvate the protein with explicit TIP3P water molecules.
Add ions
Neutralize system charge and add 0.15 M NaCl to mimic physiological salt concentration.
# Define simulation box
gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt dodecahedron
# Add water
gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top
# Prepare ion addition
gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr
# Add Na+ and Cl- ions
gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral⬇️ Step 4 — Energy Minimization
Before any simulation, you must remove steric clashes and bad contacts through energy minimization. This relaxes the structure to a local energy minimum.
integrator = steep ; Steepest descent algorithm
emtol = 1000.0 ; Stop when max force < 1000 kJ/mol/nm
emstep = 0.01 ; Initial step size
nsteps = 50000 ; Max steps
nstlist = 1
cutoff-scheme = Verlet
ns_type = grid
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
pbc = xyzgmx grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
# Check convergence
gmx energy -f em.edr -o potential.xvg
# Select: 10 (Potential Energy), then 0 to quitPlot potential.xvg — the potential energy should steadily decrease and converge to a stable negative value (typically −500,000 to −1,000,000 kJ/mol for a small protein in water).
🌡️ Step 5 — Equilibration (NVT + NPT)
Equilibration is performed in two phases to stabilize the system at target temperature and pressure before the production run.
| Phase | Ensemble | Duration | Purpose |
|---|---|---|---|
| NVT | Constant N, V, T | 100 ps | Stabilize temperature (300 K) |
| NPT | Constant N, P, T | 100 ps | Stabilize pressure (1 bar) and density |
# NVT Equilibration (temperature)
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -v -deffnm nvt
# NPT Equilibration (pressure + density)
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -v -deffnm npt🚀 Step 6 — Production MD Run
After equilibration, run the production simulation. A typical starting point is 10 nanoseconds (5,000,000 steps at 2 fs timestep).
integrator = md
nsteps = 5000000 ; 10 ns total (5M × 2fs)
dt = 0.002 ; 2 fs timestep
nstxout = 5000 ; Save coordinates every 10 ps
nstvout = 5000
nstenergy = 5000
nstlog = 5000
continuation = yes
constraint_algorithm = lincs
constraints = h-bonds
tcoupl = V-rescale
tau_t = 0.1
ref_t = 300
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 2.0
ref_p = 1.0
compressibility = 4.5e-5gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -v -deffnm md -ntmpi 1 -ntomp 4 # 4 CPU threads
# For GPU acceleration:
gmx mdrun -v -deffnm md -ntmpi 1 -ntomp 4 -gpu_id 0📊 Step 7 — Trajectory Analysis
Once your simulation completes, extract key structural and dynamic properties:
RMSD — Structural Stability
backbonegmx rms -s md.tpr -f md.xtc -o rmsd.xvg -tu ns
# Select: 4 (Backbone) for both reference and fit groupRMSF — Per-Residue Flexibility
C-alphagmx rmsf -s md.tpr -f md.xtc -o rmsf.xvg -res
# Select: 4 (C-alpha atoms)Radius of Gyration — Compactness
proteingmx gyrate -s md.tpr -f md.xtc -o gyrate.xvg
# Select: 1 (Protein)Representative MD trajectory analysis — RMSD stabilisation indicates an equilibrated simulation. A plateau in RMSD (typically below 0.3 nm for a globular protein) confirms structural stability.
⚡ Performance Tips
Use GPU
GROMACS on a modern GPU (e.g., RTX 3090) runs 50–100× faster than CPU-only.
Thread Tuning
Run gmx mdrun -ntmpi 1 -ntomp N where N = number of physical CPU cores.
Checkpoint Files
Use -cpi md.cpt to resume interrupted simulations without starting over.
Trajectory Compression
Use XTC format (lossy) for coordinate saving to significantly reduce file sizes.