Setting Up a Protein-Ligand MD Simulation in GROMACS: Step by Step

You have a docked pose from AutoDock Vina, a protein crystal structure, and a nagging question: does this ligand actually stay put when the system is allowed to move? A single docking score is a static guess. Molecular dynamics answers the dynamic question, and in practice that means learning the exact GROMACS workflow below. The steps rarely change between systems, but the details of ligand parameterization trip up nearly everyone the first time. This walkthrough takes you from a raw structure to a running production trajectory, explaining why each stage exists rather than just what to type.

The Workflow at a Glance

Every protein-ligand MD run in GROMACS follows the same skeleton: prepare structures, build topologies, define a box, add solvent and ions, minimize, equilibrate temperature and pressure, then run production dynamics. Each stage feeds the next through a small set of files, and understanding that flow keeps you from getting lost.

If this is your first exposure to MD, read what molecular dynamics simulation actually computes before diving into commands. The rest of this guide assumes you know what a force field and an integration timestep are.

Preparing the Protein and the Starting Pose

Start with a clean structure. Download your protein from the PDB, then strip crystallographic waters, ions, and any co-crystallized ligands you are not studying. Inspect the file in VMD or PyMOL and confirm there are no missing loops or alternate conformations lurking in the coordinates. Missing residues must be modeled back in before you proceed, because gmx pdb2gmx will happily build a topology around a broken chain and leave you debugging a mysterious blow-up later.

Your ligand pose usually comes from docking. If you have not generated one yet, the AutoDock Vina docking workflow produces exactly the bound conformation you need. Extract the top-ranked pose, then separate the protein and ligand into two files. The protein goes to pdb2gmx; the ligand takes a different route entirely, because standard force fields do not know your small molecule.

Building the Protein Topology

Run gmx pdb2gmx -f protein.pdb -o protein.gro -ter and select the CHARMM36 force field when prompted, plus a water model such as TIP3P. This assigns atom types, bonded parameters, and partial charges for every standard residue, writing them into topol.top.

gmx pdb2gmx -f protein_clean.pdb -o protein.gro -water tip3p -ff charmm36-jul2022

Watch the terminal output for warnings about protonation states and disulfide bonds. Histidine tautomers in particular can change hydrogen-bonding behavior in a binding site, so decide these deliberately rather than accepting defaults blindly.

Parameterizing the Ligand with CGenFF

This is the stage that separates a smooth run from a lost afternoon. Standard residues live inside CHARMM36, but your ligand does not, so you must generate its parameters with the CHARMM General Force Field (CGenFF). Save the ligand as a .mol2 file with correct bond orders and explicit hydrogens, then submit it to the CGenFF server (or a local silcsbio install) to receive a .str stream file.

Convert that stream file into GROMACS-readable topology with the community cgenff_charmm2gmx_py3.py script:

python cgenff_charmm2gmx_py3.py LIG ligand.mol2 ligand.str charmm36.ff

This emits ligand.itp, a .prm file of any new bonded parameters, and a fresh .gro. Pay close attention to the penalty scores in the .str output. Penalties above roughly 50 mean CGenFF is analogizing from poor chemical matches, and those parameters deserve validation or manual refinement before you trust any binding conclusion.

Merging Topologies and Building the System

Now stitch the two worlds together. Combine the protein .gro and ligand .gro coordinates into a single complex.gro, adjusting the atom count on the second line. Then edit topol.top to #include "ligand.itp" after the force-field include but before the water topology, add the ligand's .prm parameters near the top, and append the ligand as a molecule under the [ molecules ] directive with a count of 1.

Ordering matters here. The .prm parameters must load before any [ moleculetype ] that uses them, and the molecule list at the bottom must match the physical order of atoms in complex.gro exactly. A mismatch produces the infamous "number of coordinates does not match topology" error, which is almost always a bookkeeping mistake rather than a scientific one.

Solvation and Adding Ions

A protein in vacuum is physically meaningless, so build a solvent box. Define the box with gmx editconf, keeping at least 1.0 nm between the solute and every box edge to prevent a molecule from interacting with its own periodic image.

gmx editconf -f complex.gro -o box.gro -bt dodecahedron -d 1.0 -c
gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top

A rhombic dodecahedron box has only about 71% of the volume of a cubic box for the same minimum-image distance, so it needs roughly 30% fewer water molecules, which is a free speedup on every step of the simulation. After solvation, neutralize the system and set physiological salt. Preprocess with gmx grompp to create a run-input file, then let gmx genion swap random waters for ions:

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15

The -neutral flag cancels the system's net charge and -conc 0.15 adds 0.15 M NaCl, approximating cellular ionic strength.

Energy Minimization and Equilibration

Freshly built systems contain bad contacts: atoms placed slightly too close, awkward bond geometries from the merge, water molecules jammed into crevices. Energy minimization uses steepest descent to relax these clashes before you ever apply thermal energy, preventing the forces from launching an atom across the box on step one.

gmx grompp -f minim.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

Check the output. The maximum force Fmax should fall below the target in your .mdp file (commonly 1000 kJ/mol/nm), and the potential energy should be a large negative number that plateaus smoothly. A positive final energy signals an unresolved clash, usually traceable back to a topology or coordinate ordering slip.

Once minimization converges, you still cannot jump straight to production. The system needs to reach the target temperature and density gradually, in two phases, with the protein and ligand held in place by position restraints so the solvent equilibrates around a fixed solute.

NVT (constant number, volume, temperature) brings the system to 300 K using a thermostat such as V-rescale, typically over 100 ps. NPT (constant number, pressure, temperature) then adds a barostat like Parrinello-Rahman or C-rescale to settle the system to 1 bar and the correct density, again for 100 ps or more.

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt

Before moving on, plot temperature from the NVT run and pressure and density from the NPT run using gmx energy. Temperature should hover near 300 K, and density should converge to roughly 1000 kg/m³ for aqueous systems. If density is still drifting, extend NPT rather than pressing forward.

Production and What to Do With the Trajectory

With a properly equilibrated system, remove the position restraints and run unrestrained dynamics. A production run of at least 100 ns is a reasonable starting point for a binding study, though timescales depend entirely on your question.

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md

When the trajectory finishes, the analysis begins. Use gmx rms for backbone RMSD to judge stability, gmx rmsf for per-residue flexibility, gmx gyrate for compactness, and gmx hbond to count protein-ligand hydrogen bonds over time. A ligand that drifts out of the pocket, seen as a rising ligand RMSD, tells you the docked pose was not a stable binder no matter how good its Vina score looked. Load the .xtc into VMD to watch the binding event with your own eyes before trusting any number.

Practical Takeaways

• Ligand parameterization with CGenFF is the make-or-break step; check penalty scores and never trust parameters analogized from poor chemical matches.
• Keep topology bookkeeping meticulous: include order in topol.top and atom order in the .gro must agree exactly, or grompp will reject the system.
• Use a rhombic dodecahedron box with a 1.0 nm minimum distance to cut water count and speed up every step.
• Never skip equilibration; run NVT then NPT with position restraints and confirm temperature, pressure, and density have converged before production.
• Validate the run through backbone RMSD, ligand RMSD, and hydrogen-bond analysis rather than eyeballing a single frame.
• Treat MD as a hypothesis-testing tool: it refines docking predictions but still requires experimental confirmation.

If you have not installed the supporting toolkit yet, start with a clean AutoDock Vina installation for generating starting poses, and pair it with a dedicated conda environment so your CGenFF conversion scripts and their pinned dependencies never collide with the rest of your stack.