Computational Biochemistry · Protein–Ligand MD
Protein–Ligand Complex MD Simulation
Step-by-step workflow for simulating a protein with a small-molecule ligand · GROMACS · CHARMM36 · CGenFF
Prep
Prod.
Download & Separate the PDB
Download PDB code 3HTB (T4 lysozyme L99A/M102Q with ligand JZ4 — 2-propylphenol). Remove crystal waters and save the JZ4 ligand to a separate file — it cannot be handled by pdb2gmx since it is not in any GROMACS force field topology file.
# Then delete JZ4 lines from 3HTB_clean.pdb
Install CHARMM36 Force Field
# Produces charmm36-jul2022.ff/ subdirectory in working directory
Generate Protein Topology
# Select option 1: CHARMM36 (July 2022) from current directory
| Output file | Contents |
|---|---|
| 3HTB_processed.gro | Protein structure with H atoms added by force field |
| topol.top | Protein topology — bonds, angles, dihedrals, charges |
| posre.itp | Position restraint file for heavy atoms |
Ligand Parameterization Tools
| Force Field | Tool | Description |
|---|---|---|
| AMBER | Antechamber / acpype | GAFF parameters via Antechamber |
| CHARMM | CGenFF | Official CHARMM General Force Field server |
| GROMOS96 | ATB | Automated Topology Builder for GROMOS 54A7 |
| OPLS-AA | LigParGen | Jorgensen group server for OPLS topologies |
Step 1 — Add H Atoms with Avogadro
Open jz4.pdb in Avogadro → Build → Add Hydrogens. Save as jz4.mol2 (Sybyl Mol2 format). CHARMM36 is all-atom, so all H must be explicitly included.
Step 2 — Fix the .mol2 File
# Set all residue names and numbers to: 1 JZ4
perl sort_mol2_bonds.pl jz4.mol2 jz4_fix.mol2
Step 3 — Upload to CGenFF Server
Upload jz4_fix.mol2 at the CGenFF website. Download the GROMACS format ZIP. CGenFF provides penalty scores per parameter — values below 10 are acceptable, 10–50 need validation, above 50 need reparametrization.
Step 4 — Adapt jz4_gmx.itp
- Remove
#include "charmm36.ff/forcefield.itp" - Replace
#include "charmm36.ff/jz4_ffbonded.itp"with actual file content - Remove all lines from
; Include water topologyto end of file - Rename
[moleculetype]fromOthertoJZ4 - Change
#include "posre.itp"to#include "posre_jz4.itp"
Step 5 — Build Complex Coordinate File
cp 3HTB_processed.gro complex.gro
# Paste jz4.gro atom lines into complex.gro before box vectors
# Increment atom count on line 2 of complex.gro (+22 atoms → 2636 total)
Step 6 — Update topol.top
#include "jz4_gmx.itp" <-- MUST appear before any [moleculetype]
[molecules]
Protein_chain_A 1
JZ4 1
#include "jz4_gmx.itp" line MUST appear after the force field include but BEFORE any [moleculetype] section. Ligand parameters define new atom types that must be known before molecules are built.Define Simulation Box
A rhombic dodecahedron is used instead of a cubic box for globular proteins. Its volume is ~71% of an equivalent cubic box — requiring far fewer water molecules and saving substantial computation time.
Solvate the System
Note: GROMACS always stores coordinates in the most numerically efficient triclinic representation. The dodecahedron shape is recovered after generating a .tpr file using trjconv -ur compact.
Assemble & Neutralize
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral
# Select group 15 "SOL" when prompted
The 3HTB protein has a net charge of +6e. Six Cl⁻ ions replace water molecules to achieve electrical neutrality.
Final [molecules] directive in topol.top:
JZ4 1
SOL 10228
CL 6
gmx mdrun -v -deffnm em
The complex system converged quickly due to the high quality of both topologies. The potential energy decreases monotonically and Fmax drops well below the 1000 kJ mol⁻¹ nm⁻¹ threshold.
Restrain the Ligand
> 0 & ! a H*
> q
gmx genrestr -f jz4.gro -n index_jz4.ndx -o posre_jz4.itp -fc 1000 1000 1000
Non-hydrogen JZ4 atoms are position-restrained during equilibration. This keeps the ligand in its binding pose while water and ions equilibrate around the complex.
Thermostat Coupling — Important
Do NOT couple each molecule separately. Coupling small groups (JZ4, CL) individually causes instability due to kinetic energy fluctuations. Use a combined group:
# NOT: Protein JZ4 SOL CL ← this will blow up
Run NVT
gmx mdrun -deffnm nvt
gmx mdrun -deffnm npt
The -r nvt.gro flag provides the reference coordinates for position restraints, continuing to restrain protein and ligand while pressure equilibrates.
Run Production MD
gmx mdrun -deffnm md_0_10
Recenter & Rewrap Trajectory
# Select Protein for centering, System for output
For protein-ligand complexes, centering is more complex than a single protein — for long simulations it may be necessary to create a custom index group combining protein + ligand for centering.
Protein Backbone RMSD
JZ4 Ligand Pose RMSD
> 13 & ! a H*
> name 23 JZ4_Heavy
gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg
# Backbone for fit, JZ4_Heavy for RMSD
Protein–Ligand H-bond Distance
The JZ4 hydroxyl (O10) can hydrogen bond with the Gln102 carbonyl OE1. Average distance over 10 ns: 0.43 ± 0.17 nm. Over last 5 ns: 0.30 ± 0.03 nm — consistent with a stable H-bond (<0.36 nm threshold).
Protein–Ligand Interaction Energy
gmx mdrun -deffnm ie -rerun md_0_10.xtc -nb cpu
gmx energy -f ie.edr -o interaction_energy.xvg
# Select Coul-SR:Protein-jz4 and LJ-SR:Protein-jz4
Unlike standard amino acids, small-molecule ligands must be parameterized using external tools before they can be included in an MD simulation. Each tool is designed to be consistent with a specific force field family. Always validate penalty scores before using a topology.
| Penalty Score | Interpretation | Action |
|---|---|---|
| < 10 | Excellent — well-characterized | Use directly |
| 10 – 50 | Moderate — assigned by analogy | Validate topology |
| > 50 | Poor — unreliable | Manual reparametrization required |