GROMACS Molecular Dynamics Simulation Protocol: Step-by-Step Protein-Ligand Simulation

Molecular dynamics (MD) simulations have emerged as a premier computational methodology in computer-aided drug design (CADD), structural pharmacology, and biophysics. While static crystallographic or cryo-EM structures provide atomic snapshots of biological macromolecules, biomolecules execute their physiological roles through continuous thermal fluctuations, induced-fit conformational transitions, and transient solvent reorganizations.

This publication-grade protocol guides researchers through setting up, parameterizing, simulating, and analyzing an unrestrained all-atom protein-ligand complex using GROMACS 2024 and the AMBER force field suite. Whether investigating inhibitor residence times, validating virtual screening docking poses, or calculating relative binding affinities, following a standardized protocol ensures numerical stability and reproducible trajectories.

For researchers and commercial drug discovery teams seeking turnkey execution, high-throughput cluster runs, or advanced free-energy perturbation (FEP/TI), consult the Bioinformatics Daily MD Simulation Services for dedicated computational assistance.

1. Introduction & Real-World Biological Context

In modern medicinal chemistry, static docking poses generated by scoring functions (e.g., AutoDock Vina, Glide) often score favorably in terms of predicted geometry yet fail downstream when subjected to dynamic solvated environments. Molecular docking relies on rigid or semi-flexible receptor models and approximate scoring functions that drastically simplify entropic penalties and solvent displacement thermodynamics.

Molecular dynamics bridges this divide by solving Newton’s classical equations of motion:

F_i = m_i * (d^2 r_i / dt^2) = -grad_i V(r_1, r_2, ..., r_N)

Where:

• m_i represents the mass of atom i,
• r_i defines the Cartesian spatial coordinate vector of atom i,
• V(r_1, ..., r_N) denotes the empirical force field potential energy function, encompassing bonded interactions (bond stretching, angle bending, torsional dihedral rotations) and non-bonded potentials (Coulombic electrostatics and Lennard-Jones 12-6 van der Waals dispersion).

Key Questions Answered by Protein-Ligand MD Simulations

1. Pose Stability & Residence Verification: Does the docked ligand remain firmly anchored within the orthosteric or allosteric pocket, or does it dislodge into the solvent bulk within 20–50 ns?
2. Induced Fit & Cryptic Pocket Opening: How do dynamic loop movements (e.g., the Gly-rich loop in protein kinases or the active-site flap in HIV-1 protease) flex to accommodate or lock the small molecule? Read our deep-dive on Protein-Ligand Interaction Mechanics.
3. Solvation Dynamics & Water Bridges: Are structural water molecules trapped between polar ligand pharmacophores and backbone amide groups, stabilizing the bound complex?
4. Binding Energetics: What is the approximate free energy of binding (ΔG_bind) when computed via Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) or Generalized Born (MM-GBSA) models?

2. Prerequisites & Environment Setup

To guarantee reproducibility and eliminate library mismatches, execute all workflows inside a dedicated conda environment containing GROMACS, Open Babel, AmberTools (for ligand charge assignment), and ACPYPE (AnteChamber PYthon Parser InterfacE).

2.1 Conda Environment Construction

# Create isolated computational biology environment
conda create -n gromacs-cadd python=3.10 -y
conda activate gromacs-cadd

# Install GROMACS (GPU-enabled build where CUDA driver is available)
conda install -c conda-forge gromacs=2024.1 -y

# Install AmberTools for antechamber and parmchk2
conda install -c conda-forge ambertools=23.3 -y

# Install ACPYPE and Open Babel for automated ligand parameterization
conda install -c conda-forge acpype openbabel -y

# Install Python scientific stack for trajectory post-processing
pip install MDAnalysis MDAnalysisTests numpy pandas matplotlib seaborn gmx_MMPBSA

2.2 Hardware Verification & Acceleration

Ensure that your host machine recognizes your NVIDIA GPU for PME (Particle Mesh Ewald) electrostatics and non-bonded force offloading:

# Verify NVIDIA GPU driver and CUDA availability
nvidia-smi

# Check GROMACS acceleration capabilities
gmx -version | grep -E "(Precision|GPU support|SIMD)"

A properly configured build will report Precision: single, GPU support: CUDA, and active AVX2/AVX-512 SIMD vectorization.

3. Input Data Format & Preprocessing

The primary input files for a protein-ligand MD setup consist of:

1. Protein Receptor Model: High-resolution X-ray or cryo-EM structure file (.pdb) devoid of non-crystallographic water molecules, buffer additives, and alternate conformations. Alternatively, an AlphaFold2 predicted model can be employed if experimental coordinates are unavailable.
2. Small Molecule Ligand: Clean 3D coordinate file with correct protonation states and bond orders (.sdf or .mol2).

3.1 Raw PDB Structure Cleaning

Experimental PDB files contain crystallographic artifacts (e.g., glycerol, sulfate ions, PEG, azide) and missing loop regions.

# 1. Download target PDB (e.g., Human CDK2 in complex with an inhibitor, PDB ID: 1H1Q)
wget https://files.rcsb.org/download/1H1Q.pdb -O raw_complex.pdb

# 2. Extract protein coordinates, stripping heteroatoms and crystallographic waters
grep -E "^ATOM" raw_complex.pdb > protein_clean.pdb

# 3. Extract the ligand coordinates into a distinct file (e.g., residue name 'LBD')
grep -E "^HETATM.*LBD" raw_complex.pdb > ligand_raw.pdb

# 4. Check for missing atoms or incomplete residues using pdb4amber
pdb4amber -i protein_clean.pdb -o protein_amber.pdb --dry

3.2 Ligand Valency and Protonation

Small molecule force fields (e.g., GAFF2) require explicit hydrogens and net formal charges corresponding to physiological pH (7.4):

# Convert PDB ligand to SDF with proper bond orders and explicit hydrogens
obabel ligand_raw.pdb -O ligand_h.sdf -p 7.4 --gen3d

Review the resulting structure in PyMOL or UCSF Chimera to ensure nitrogen atom hybridization, carboxylic acid deprotonation, and chiral center orientations match experimental chemistry.

4. The Step-by-Step Computational Workflow

[Protein PDB] ---> [gmx pdb2gmx] -------------\
                                              +--> [gmx solvate] -> [gmx genion] -> [EM] -> [NVT] -> [NPT] -> [Production MD]
[Ligand SDF]  ---> [Antechamber/ACPYPE] ------/

Step 1: Generate Protein Topology (pdb2gmx)

GROMACS converts Cartesian coordinates into mathematical topologies (masses, charges, bond types, dihedral angles) using force field definition files:

# Generate topology with AMBER99SB-ILDN and TIP3P explicit water model
gmx pdb2gmx -f protein_amber.pdb \
            -o protein_processed.gro \
            -p topol.top \
            -i posre.itp \
            -ff amber99sb-ildn \
            -water tip3p \
            -ignh

Flag Explanations:

• -ff amber99sb-ildn: Assigns the modern AMBER99SB force field with updated isoleucine, leucine, and aspartate torsional potentials.
• -water tip3p: Selects the 3-site transferable intermolecular potential water model.
• -ignh: Strips existing hydrogen atoms and recalculates them according to standard residue ionization states at pH 7.0.

Step 2: Parameterize the Small Molecule Ligand with ACPYPE

Because standard macromolecular force fields do not contain parameters for arbitrary synthetic drug scaffolds, we parameterize the ligand using the General AMBER Force Field (GAFF2) with Semi-Empirical AM1-BCC partial atomic charges:

# Run ACPYPE to compute partial charges and produce GROMACS .itp & .gro
acpype -i ligand_h.sdf \
       -b LIG \
       -c bcc \
       -n 0 \
       -a gaff2

Parameters Explained:

• -b LIG: Re-names the ligand residue to LIG.
• -c bcc: Executes the Austin Model 1 Bond Charge Correction (AM1-BCC) method.
• -n 0: Specifies the net molecular charge (adjust to -1, +1, or +2 depending on functional groups).
• -a gaff2: Employs GAFF2 parameter tables.

ACPYPE outputs a folder named LIG.acpype/. The essential files are:

• LIG_GMX.gro: Cartesian coordinates of the ligand in GROMACS format.
• LIG_GMX.itp: Molecular topology containing atom types, charges, bonded pairs, and dihedrals.
• LIG_GMX_posre.itp: Positional restraint parameters for heavy atoms during equilibration.

Step 3: Complex Assembly and Master Topology Integration

We merge the protein and ligand coordinate files into a single complex coordinate file and include the ligand topology inside topol.top.

# 1. Combine coordinates: copy protein coordinates and append ligand coordinates
grep -v "ENDMDL" protein_processed.gro | head -n -1 > complex.gro
tail -n +3 LIG.acpype/LIG_GMX.gro | head -n -1 >> complex.gro

# 2. Update atom count in line 2 of complex.gro
protein_atoms=$(sed -n '2p' protein_processed.gro | tr -d ' ')
ligand_atoms=$(sed -n '2p' LIG.acpype/LIG_GMX.gro | tr -d ' ')
total_atoms=$((protein_atoms + ligand_atoms))
sed -i "2s/.*/ ${total_atoms}/" complex.gro

# 3. Append original simulation box dimensions from protein_processed.gro to bottom
tail -n 1 protein_processed.gro >> complex.gro

Now update topol.top using a text editor to include the ligand parameters:

; In topol.top:
; Include forcefield parameters
#include "amber99sb-ildn.ff/forcefield.itp"

; INCLUDE LIGAND PARAMETERS HERE (BEFORE [ moleculetype ])
#include "LIG.acpype/LIG_GMX.itp"

; [ moleculetype ] definition for Protein follows...
...
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; INCLUDE LIGAND POSITIONAL RESTRAINTS
#ifdef POSRES
#include "LIG.acpype/LIG_GMX_posre.itp"
#endif

; Include water topology
#include "amber99sb-ildn.ff/tip3p.itp"

; Include topology for ions
#include "amber99sb-ildn.ff/ions.itp"

[ molecules ]
; Compound        #mols
Protein_chain_A     1
LIG                 1

Step 4: Define Periodic Simulation Unit Cell

We define a rhombic dodecahedron box. Dodecahedral geometries save ~29% solvent volume compared to standard cubic boxes, drastically accelerating computation:

# Center complex in a rhombic dodecahedron with 1.2 nm minimum buffer to boundary
gmx editconf -f complex.gro \
             -o complex_box.gro \
             -bt dodecahedron \
             -d 1.2 \
             -c

Step 5: Explicit Solvation

Fill the unit cell with TIP3P water molecules:

# Solvate the simulation box
gmx solvate -cp complex_box.gro \
            -cs spc216.gro \
            -o complex_solv.gro \
            -p topol.top

gmx solvate automatically updates the [ molecules ] directive at the foot of topol.top with the added SOL count.

Step 6: Ion Addition and System Neutralization

Biological proteins function at physiological ionic strength (~150 mM NaCl). We neutralize residual system net charges and add counter-ions using gmx genion:

Create a preprocessing parameter file ions.mdp:

; ions.mdp - Preprocessor parameter file for genion
integrator  = steep
nsteps      = 500
emtol       = 1000.0
emstep      = 0.01
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb    = 1.0
rvdw        = 1.0
pbc         = xyz

Compile and introduce ions:

# Compile portable binary run input (.tpr)
gmx grompp -f ions.mdp -c complex_solv.gro -p topol.top -o ions.tpr -maxwarn 1

# Replace solvent molecules with ions to reach 0.15 M NaCl concentration
# Select group 'SOL' (typically group 13 or 15) when prompted
echo "SOL" | gmx genion -s ions.tpr \
                        -o complex_solv_ions.gro \
                        -p topol.top \
                        -pname NA \
                        -nname CL \
                        -neutral \
                        -conc 0.15

Step 7: Energy Minimization (EM)

Steepest descent energy minimization relieves steric clashes and abnormal bond lengths introduced during crystal packing and solvation:

Create em.mdp:

; em.mdp - Steepest Descent Minimization
integrator  = steep
emtol       = 100.0        ; Stop when max force < 100.0 kJ/(mol*nm)
emstep      = 0.01         ; Energy step size (nm)
nsteps      = 50000        ; Maximum number of steps

nstenergy   = 500
cutoff-scheme = Verlet
ns_type     = grid
coulombtype = PME
rcoulomb    = 1.0
rvdw        = 1.0
pbc         = xyz

Execute energy minimization:

# Compile and run energy minimization
gmx grompp -f em.mdp -c complex_solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em -nb gpu

Verify convergence: the maximum force (Fmax) must fall below 1000 kJ/(mol*nm), and potential energy must remain negative (typically -10^5 to -10^6 kJ/mol).

Step 8: Two-Phase Restrained Equilibration (NVT & NPT)

Equilibration occurs in two distinct thermal stages to prevent the ligand and protein backbone from collapsing while solvent molecules establish thermal motion.

Phase 1: NVT Ensemble (Constant Volume and Temperature, 100 ps)

Create nvt.mdp:

; nvt.mdp - NVT Ensemble Equilibration
define      = -DPOSRES     ; Enable position restraints on protein & ligand
integrator  = md
dt          = 0.002        ; 2 femtoseconds
nsteps      = 50000        ; 100 picoseconds

nstxout-compressed = 5000
nstenergy   = 1000
nstlog      = 1000

continuation = no
constraint_algorithm = lincs
constraints = h-bonds
lincs_iter  = 1
lincs_order = 4

cutoff-scheme = Verlet
ns_type     = grid
coulombtype = PME
rcoulomb    = 1.0
rvdw        = 1.0

; Temperature coupling: Velocity-rescaling thermostat
tcoupl      = V-rescale
tc-grps     = Protein_LIG Water_and_ions
tau_t       = 0.1   0.1
ref_t       = 300   300

; Velocity generation
gen_vel     = yes
gen_temp    = 300
gen_seed    = -1

pbc         = xyz

To prevent temperature coupling errors between heterogeneous groups, build an index file clustering Protein and LIG into one thermal group:

# Create unified temperature coupling group
gmx make_ndx -f em.gro -o index.ndx <<EOF
1 | 13
name 14 Protein_LIG
q
EOF

# Run NVT equilibration
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -v -deffnm nvt -nb gpu

Phase 2: NPT Ensemble (Constant Pressure and Temperature, 100 ps)

Create npt.mdp:

; npt.mdp - NPT Ensemble Equilibration
define      = -DPOSRES
integrator  = md
dt          = 0.002
nsteps      = 50000        ; 100 ps

nstxout-compressed = 5000
nstenergy   = 1000
nstlog      = 1000

continuation = yes
constraint_algorithm = lincs
constraints = h-bonds

cutoff-scheme = Verlet
coulombtype = PME
rcoulomb    = 1.0
rvdw        = 1.0

tcoupl      = V-rescale
tc-grps     = Protein_LIG Water_and_ions
tau_t       = 0.1   0.1
ref_t       = 300   300

; Pressure coupling: Parrinello-Rahman barostat
pcoupl           = Parrinello-Rahman
pcoupltype       = isotropic
tau_p            = 2.0
ref_p            = 1.0
compressibility  = 4.5e-5
refcoord_scaling = com

pbc         = xyz

Execute NPT:

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr
gmx mdrun -v -deffnm npt -nb gpu

Step 9: Unrestrained Production Simulation (100 ns)

Create md.mdp:

; md.mdp - 100 ns Production MD Run
integrator  = md
dt          = 0.002        ; 2 fs
nsteps      = 50000000     ; 50,000,000 steps = 100,000 ps = 100 ns

nstxout-compressed = 50000 ; Write coordinates every 100 ps (1000 frames total)
compressed-x-grps  = System
nstenergy   = 5000
nstlog      = 5000

continuation = yes
constraint_algorithm = lincs
constraints = h-bonds

cutoff-scheme = Verlet
coulombtype = PME
rcoulomb    = 1.0
rvdw        = 1.0

tcoupl      = V-rescale
tc-grps     = Protein_LIG Water_and_ions
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

Launch production run:

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_100ns.tpr
gmx mdrun -v -deffnm md_100ns -nb gpu -pme gpu -bonded gpu

5. Data Visualization & Result Interpretation

Before analyzing coordinates, you must post-process the trajectory to eliminate periodic boundary condition (PBC) jumping, where the protein appears split across unit cell faces.

# 1. Cluster complex and center in unit cell
echo "Protein_LIG System" | gmx trjconv -s md_100ns.tpr \
                                       -f md_100ns.xtc \
                                       -o md_center.xtc \
                                       -pbc mol \
                                       -center \
                                       -n index.ndx

# 2. Extract Protein Backbone RMSD
echo "Backbone Backbone" | gmx rms -s md_100ns.tpr -f md_center.xtc -o rmsd_protein.xvg

# 3. Extract Ligand RMSD fitted to Protein Backbone
echo "Backbone LIG" | gmx rms -s md_100ns.tpr -f md_center.xtc -o rmsd_ligand.xvg -n index.ndx

# 4. Extract per-residue RMSF (Root Mean Square Fluctuation)
echo "C-alpha" | gmx rmsf -s md_100ns.tpr -f md_center.xtc -o rmsf.xvg -res

# 5. Extract Hydrogen Bond Count between Protein and Ligand
echo "Protein LIG" | gmx hbond -s md_100ns.tpr -f md_center.xtc -num hbond.xvg -n index.ndx

5.1 Publication-Quality R Visualization Script

Use this copy-pasteable R script using ggplot2 and patchwork to visualize RMSD, RMSF, and hydrogen bonds:

library(ggplot2)
library(dplyr)
library(patchwork)

# Function to read GROMACS XVG files safely
read_xvg <- function(filepath, col_names) {
  lines <- readLines(filepath)
  data_lines <- grep("^[^#@]", lines, value = TRUE)
  df <- read.table(text = data_lines, col.names = col_names)
  return(df)
}

# 1. Load Data
rmsd_prot <- read_xvg("rmsd_protein.xvg", c("Time_ps", "RMSD_nm")) %>%
  mutate(Time_ns = Time_ps / 1000, Type = "Receptor Backbone")

rmsd_lig <- read_xvg("rmsd_ligand.xvg", c("Time_ps", "RMSD_nm")) %>%
  mutate(Time_ns = Time_ps / 1000, Type = "Ligand (Pocket-Fitted)")

rmsd_comb <- rbind(rmsd_prot, rmsd_lig)
rmsf_data <- read_xvg("rmsf.xvg", c("Residue", "RMSF_nm"))
hbond_data <- read_xvg("hbond.xvg", c("Time_ps", "HBonds", "Pairs")) %>%
  mutate(Time_ns = Time_ps / 1000)

# Plot A: Dual RMSD Trajectory
p1 <- ggplot(rmsd_comb, aes(x = Time_ns, y = RMSD_nm, color = Type)) +
  geom_line(alpha = 0.85, linewidth = 0.7) +
  scale_color_manual(values = c("Receptor Backbone" = "#2563eb", "Ligand (Pocket-Fitted)" = "#e11d48")) +
  theme_minimal(base_size = 13) +
  labs(title = "Structural Stability (RMSD)",
       x = "Simulation Time (ns)",
       y = "RMSD (nm)",
       color = NULL) +
  theme(legend.position = "top",
        panel.grid.minor = element_blank())

# Plot B: Protein C-alpha RMSF Fluctuations
p2 <- ggplot(rmsf_data, aes(x = Residue, y = RMSF_nm)) +
  geom_line(color = "#059669", linewidth = 0.8) +
  geom_area(fill = "#10b981", alpha = 0.2) +
  theme_minimal(base_size = 13) +
  labs(title = "Local Flexibility (RMSF)",
       x = "Residue Number",
       y = "RMSF (nm)") +
  theme(panel.grid.minor = element_blank())

# Plot C: Intermolecular Hydrogen Bond Persistence
p3 <- ggplot(hbond_data, aes(x = Time_ns, y = HBonds)) +
  geom_line(color = "#7c3aed", alpha = 0.6) +
  geom_smooth(method = "loess", color = "#4c1d95", span = 0.2, se = FALSE) +
  theme_minimal(base_size = 13) +
  labs(title = "Protein-Ligand Hydrogen Bonds",
       x = "Simulation Time (ns)",
       y = "Number of H-Bonds") +
  theme(panel.grid.minor = element_blank())

# Compose publication multi-panel
final_figure <- p1 / (p2 | p3) + plot_annotation(tag_levels = 'A')
ggsave("Figure_MD_Publication_Panel.pdf", final_figure, width = 11, height = 8, dpi = 300)

5.2 How to Interpret the Output Figures

• Backbone RMSD plateauing below 0.25–0.30 nm (2.5–3.0 Å): Confirms the global tertiary fold converged without global denaturation.
• Ligand RMSD plateauing with fluctuations under 0.15 nm: Indicates the ligand remains locked in its active site. Large jumps (> 0.4 nm) signify unbinding or binding-pose transitions.
• RMSF Peaks (> 0.2 nm): Naturally correspond to exposed loops and terminal residues. The catalytic pocket residues should exhibit low RMSF (under 0.1 nm), reflecting rigid binding cavity walls.
• Continuous H-Bonds (1–3 steady bonds): Verifies the key pharmacophoric hydrogen bonds remain intact across solvent exposure.

6. Common Errors & Troubleshooting

7. Frequently Asked Questions (FAQ)

What hardware is required for a 100 ns GROMACS protein-ligand simulation?

An NVIDIA RTX 3080 or A100 GPU processes ~100–350 ns per day for a typical 60,000-atom solvated system, whereas CPU-only execution achieves fewer than 10 ns per day.

Can GROMACS simulate covalent protein-ligand complexes?

Yes, but covalent linkages require defining custom chemical bonds, angles, and dihedrals directly between the reactive residue (e.g., Cys, Ser) and the inhibitor within the master topology file rather than treating the ligand as an isolated molecule.

How does MM-PBSA calculate binding free energy?

The gmx_MMPBSA tool evaluates the thermodynamic cycle: Delta G_bind = mean(G_complex) - mean(G_receptor) - mean(G_ligand), calculating gas-phase molecular mechanics energies plus polar (Poisson-Boltzmann) and non-polar solvation terms across trajectory frames.

Why should NVT and NPT equilibrations use positional restraints?

Positional restraints (-DPOSRES) hold protein backbone and ligand heavy atoms fixed while allowing random water and ion coordinates to relax around macromolecular crevices, preventing system blowup.

Continue Your Training & Professional MD Services

Ready to take your structural biology workflows further?

• Explore our foundational guide on Protein-Ligand Interactions & Computational Docking.
• Learn how to predict novel target geometries using AlphaFold2 & ColabFold.
• Need high-throughput production simulations, binding free-energy profiles, or dedicated computational power? Submit your target to our certified MD Simulation Services.