Protein–ligand MD in GROMACS 2026 with AMBER99SB-ILDN and GAFF2 (ACPYPE): HIV-1 protease + indinavir

This tutorial builds and runs a complete protein–ligand simulation: HIV-1 protease with its inhibitor indinavir (ligand code MK1, PDB 1HSG). The protein uses AMBER99SB-ILDN. The ligand uses GAFF2 with AM1-BCC charges, generated with ACPYPE. Every command below was run with GROMACS 2026.3 on a CPU-only machine (8 threads), and every number quoted is from that run.

It follows on from the Vina docking tutorial, which prepared the same protein and ligand. Here we start from the crystal pose. Starting MD from a docked pose works the same way, but it is a hypothesis test, not a refinement step.

Why AMBER + GAFF2 rather than CHARMM36 + CGenFF? Both are good choices. The AMBER route can be done entirely with open-source tools on your own machine (AmberTools + ACPYPE), while CGenFF parameters come from a web server. What matters is not mixing families: AMBER protein force fields go with GAFF/GAFF2 ligands, and CHARMM36 goes with CGenFF.

0Environment

micromamba create -n gmx -c conda-forge gromacs=2026.3 acpype ambertools gmx_mmpbsa
micromamba activate gmx

The conda-forge GROMACS build is fine for learning and small systems. For production on a cluster, use a GPU build compiled for your hardware.

1Split the complex and choose protonation states

wget https://files.rcsb.org/download/1HSG.pdb -O 1hsg.pdb
grep -E '^(ATOM|TER)' 1hsg.pdb > protein.pdb && echo END >> protein.pdb
grep '^HETATM.\{11\}MK1' 1hsg.pdb > mk1_xtal.pdb    # 45 heavy atoms, no H

Crystal waters are removed here, including the structural “flap water” (HOH 308). It will be replaced by bulk water during solvation, and in practice a water molecule often re-enters that site within nanoseconds. For a careful study, keep HOH 308 explicitly.

pdb2gmx and the catalytic aspartates

gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top -i posre.itp \
            -ff amber99sb-ildn -water tip3p -ignh

pdb2gmx picks histidine tautomers from the hydrogen-bond network (here HISE for His69 in both chains) and, by default, makes every Asp, Glu, Lys and Arg charged. The result has total charge +4.

This default is a real modelling decision for HIV protease. The catalytic dyad (Asp25/Asp25′) shares a proton when a hydroxyl-containing inhibitor such as indinavir is bound. Most MD studies of protease–inhibitor complexes protonate one of the two aspartates. We kept the default (both deprotonated) to keep this tutorial short. For publication-quality work, run pdb2gmx with -asp and choose the neutral form (ASH) for Asp25 of one chain, and justify the choice, for example with literature on your specific inhibitor or with constant-pH/pKa calculations.

2Ligand: hydrogens and bond orders first

ACPYPE (through antechamber) needs a ligand with all hydrogens and correct bond orders. A PDB HETATM record has neither. We reuse the RDKit procedure from the docking tutorial: assign bond orders from the indinavir SMILES, read stereochemistry from the crystal coordinates, and add hydrogens in place. That gives mk1_xtal_h.sdf with 92 atoms (45 heavy + 47 H), formal charge 0, in the crystal frame. Then convert it to MOL2:

obabel mk1_xtal_h.sdf -O mk1_h.mol2

Check the net charge now. Indinavir is neutral at physiological pH (its most basic nitrogen has a pKa of about 6). If your ligand is charged, pass that charge to ACPYPE with -n.

3GAFF2 / AM1-BCC parameters with ACPYPE

acpype -i mk1_h.mol2 -b MK1 -c bcc -n 0 -a gaff2 -o gmx

ACPYPE runs antechamber (atom typing and the AM1-BCC charge procedure, which calls the sqm semi-empirical program), parmchk2 (fills missing parameters) and tleap, and then converts the result to GROMACS format. For this 92-atom ligand it took 15 min 57 s on one CPU core. Almost all of that time was the AM1 geometry optimisation inside sqm. Small fragments take seconds, but large flexible drugs can take tens of minutes, so start this step early.

Output in MK1.acpype/:

The assigned GAFF2 types are worth a look: 17 ca (aromatic carbons), 2 ns/n amides, n3 for the piperazine nitrogens, nb for the pyridine nitrogen, and oh/ho for the two hydroxyls. A wrong type here, such as an aromatic ring typed as aliphatic, is the usual cause of a ligand that distorts during MD.

4Merge ligand and protein

GROMACS requires all [ atomtypes ] to come before the first [ moleculetype ]. The protein [ moleculetype ] sections are pulled in by the chain .itp includes, so the ligand’s atom types must be split into their own file and included right after the force field.

# split the ACPYPE itp into atomtypes + moleculetype
itp = open("MK1.acpype/MK1_GMX.itp").read().splitlines()
i = next(k for k, l in enumerate(itp) if l.startswith("[ moleculetype"))
open("MK1_atomtypes.itp", "w").write("\n".join(itp[:i]) + "\n")
open("MK1.itp", "w").write("\n".join(itp[i:]).replace(" UNL ", " MK1 ") + "\n")

ACPYPE names the residue UNL in some files and MK1 in others. Renaming it consistently avoids confusion in index groups later. Edit topol.top so that the top reads:

#include "amber99sb-ildn.ff/forcefield.itp"

; GAFF2 atom types for the ligand (must come before any [ moleculetype ])
#include "MK1_atomtypes.itp"

#include "topol_Protein_chain_A.itp"
#include "topol_Protein_chain_B.itp"

; Ligand topology (GAFF2 / AM1-BCC via ACPYPE)
#include "MK1.itp"
#ifdef POSRES
#include "posre_MK1.itp"
#endif

#include "amber99sb-ildn.ff/tip3p.itp"

Add the ligand to the end of [ molecules ]:

[ molecules ]
Protein_chain_A     1
Protein_chain_B     1
MK1                 1

Putting posre_MK1.itp under the same POSRES define as the protein means define = -DPOSRES restrains both during equilibration.

Build the complex coordinate file by appending the ligand lines to the protein .gro file. The ligand lines go after the protein, in the same order as in [ molecules ], and the atom count on line 2 must be updated: 3,128 + 92 = 3,220.

P = open("protein.gro").read().splitlines()
L = open("MK1.acpype/MK1_GMX.gro").read().splitlines()
lig = [l.replace("UNL  ", "MK1  ") for l in L[2:-1]]
body = P[2:-1] + lig
lines = ["HIV-1 PR + MK1", str(len(body))] + body + [P[-1]]
open("complex.gro", "w").write("\n".join(lines) + "\n")

5Box, water and ions

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

A rhombic dodecahedron with 1.0 nm clearance gives a box volume of 380.6 nm³, about 29% smaller than the equivalent cube. That saves the same fraction of water and compute. solvate added 11,519 waters.

gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr -maxwarn 1
echo SOL | gmx genion -s ions.tpr -o ionized.gro -p topol.top \
                      -pname NA -nname CL -neutral -conc 0.15

-maxwarn 1 is used here for one specific, understood warning: “You are using Ewald electrostatics in a system with net charge”. That is precisely what genion is about to fix. This is the only step where we allow a warning. genion replaced waters with 34 Na⁺ and 38 Cl⁻ (neutralising the +4 and adding about 0.15 M NaCl). The final system has 37,633 atoms.

6Parameter files

The four .mdp files were generated with the Build tab of the mdp checker using the AMBER settings: PME, 1.0 nm cut-offs with potential-shift, DispCorr = EnerPres, a 2 fs time step with constraints = h-bonds, V-rescale at 300 K, and C-rescale at 1 bar. They use two coupling groups, Protein_MK1 and Water_and_ions. The reasoning behind each line is in Choosing GROMACS mdp parameters.

The coupling groups need an index file. The ligand goes into the protein group: a 92-atom ligand on its own is too small to be thermostatted as a separate group.

printf '1 | 13\nname 21 Protein_MK1\nq\n' | gmx make_ndx -f em.gro -o index.ndx

(Group 1 is Protein and 13 is MK1 in the default list. Check the numbers printed by make_ndx for your system.)

7Energy minimisation

gmx grompp -f em.mdp -c ionized.gro -p topol.top -o em.tpr      # -maxwarn 0
gmx mdrun -deffnm em

Steepest descent converged to Fmax < 1000 in 989 steps (11 s). Final potential energy: −5.93 × 10⁵ kJ/mol. Maximum force: 984 kJ mol⁻¹ nm⁻¹, on a protein atom. A ligand atom carrying the maximum force after EM would point to a clash or a bad topology.

8Equilibration: NVT, then NPT

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -deffnm nvt

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 -deffnm npt

Both stages ran for 100 ps with heavy atoms of the protein and ligand restrained (-DPOSRES; -r gives the reference coordinates). grompp accepted all four stages with -maxwarn 0. The only message was the harmless NOTE about removing centre-of-mass motion with restraints.

Check the equilibration with gmx energy. Use -b 20 to drop the first 20 ps:

echo Temperature | gmx energy -f nvt.edr -b 20 -o nvt_T.xvg
printf 'Temperature\nPressure\nDensity\n' | gmx energy -f npt.edr -b 20 -o npt.xvg

These are the numbers from this run (8 CPU threads, no GPU):

The pressure fluctuation of ±150 bar is not a problem. Instantaneous pressure in a box of about 37,000 atoms fluctuates by hundreds of bar. Judge the average against the 1 bar target, and only over a long enough window. A 2 bar average over 80 ps is consistent with 1 bar.

The density is a better indicator. It went from 1013 kg/m³ in the first 10 ps to a steady 1016 kg/m³, which is reasonable for TIP3P with a protein and 0.15 M NaCl at 300 K. If the density is still drifting at the end of NPT, extend NPT before going on to production.

9Production

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

Production has no restraints (no -r, no -DPOSRES). It continues from the NPT checkpoint, so velocities and box size carry over (continuation = yes, gen_vel = no). grompp passed with -maxwarn 0 and printed no notes beyond the informational rlist message.

The demonstration run is 0.5 ns (250,000 steps of 2 fs), writing compressed coordinates every 10 ps (51 frames). It ran at 27.3 ns/day on 8 CPU threads, which is about 26 minutes. Over the run, the average temperature was 300.0 K, the average pressure 7 ± 5 bar (error estimate from gmx energy, with fluctuations of ±143 bar), and the density 1014.7 kg/m³.

0.5 ns is enough to check that the setup works. It is not enough to study the binding. For real work, plan at least 100 ns per replica and three or more replicas started with different velocities. On a single modern GPU this system runs at hundreds of ns/day, so that is a matter of days, not months.

10Is the complex stable? First analysis

Make the trajectory whole and centred before analysing or viewing it:

printf '1\n0\n' | gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -pbc mol -center

-pbc mol -center with group 1 (Protein) for centring and 0 (System) for output keeps molecules whole and puts the protein in the middle of the box. Next, make an index group with the ligand heavy atoms only (hydrogens add rotational noise to the RMSD):

printf '13 & ! a H*\nname 22 MK1_heavy\nq\n' | \
    gmx make_ndx -f md.tpr -n index.ndx -o ana.ndx

# protein backbone RMSD (fit and calculate on group 4, Backbone)
printf '4\n4\n'  | gmx rms -s md.tpr -f md_center.xtc -n ana.ndx -o rmsd_bb.xvg  -tu ns
# ligand RMSD *in the binding site*: fit on the backbone (4), compute on MK1_heavy (22)
printf '4\n22\n' | gmx rms -s md.tpr -f md_center.xtc -n ana.ndx -o rmsd_lig.xvg -tu ns

The second command is the one that answers “does the ligand stay in its pose?”. Fitting on the protein and measuring the ligand without a second fit picks up translation and rotation of the ligand inside the pocket. If you fit on the ligand itself, you only see its internal flexibility, and a ligand that leaves the pocket can still look fine.

Results over the 0.5 ns (the reference is the md.tpr structure, i.e. the end of NPT):

The ligand stays in its crystal binding mode. Its RMSD rises slightly in the second half and stays well below the ~0.2–0.3 nm level where you would start to look for a change of pose. That is all a 0.5 ns run can tell you.

A GROMACS 2026 pitfall with ACPYPE ligands: the new gmx hbond (the trajectory-analysis-framework version) stopped with

Inconsistency in user input:
Selection MK1' has no donors AND has no acceptors! Nothing to be done.

That happens even though indinavir has several OH and NH groups. Most likely, the new tool did not recognise the GAFF2 atom types in the ACPYPE topology as donors or acceptors. We did not investigate further. The legacy implementation finds them, so the H-bond count above comes from:

printf '1\n13\n' | gmx hbond-legacy -s md.tpr -f md_center.xtc -n ana.ndx -num hbnum.xvg

If you use the new tool with a GAFF ligand, check the donor and acceptor counts it prints before you trust a zero.

To estimate binding energies from this trajectory, continue with MM-GBSA/MM-PBSA with gmx_MMPBSA. That tutorial also explains why the absolute numbers should not be read as binding free energies.

Checklist

• The ligand has correct bond orders, every hydrogen, and the intended net charge before ACPYPE.
• sqm.out ends with “Calculation Completed”. You have read frcmod for estimated parameters.
• The ligand’s [ atomtypes ] are included before any [ moleculetype ].
• The ligand is listed in [ molecules ] in the same order as in the .gro file.
• Protonation states are a deliberate choice, especially for catalytic residues.
• grompp passes with -maxwarn 0 at every stage except the ion-placement .tpr.
• EM converged, and the maximum force is not on the ligand.
• Temperature, pressure and density reach their targets during equilibration.
• Ligand RMSD in the binding site has been checked, not only protein RMSD.

References

1. Abraham MJ, et al. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015;1–2:19–25. doi:10.1016/j.softx.2015.06.001
2. Sousa da Silva AW, Vranken WF. ACPYPE – AnteChamber PYthon Parser interfacE. BMC Res Notes 2012;5:367. doi:10.1186/1756-0500-5-367
3. Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. Development and testing of a general amber force field. J Comput Chem 2004;25:1157–1174. doi:10.1002/jcc.20035
4. Jakalian A, Jack DB, Bayly CI. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J Comput Chem 2002;23:1623–1641. doi:10.1002/jcc.10128
5. Lindorff-Larsen K, et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins 2010;78:1950–1958. doi:10.1002/prot.22711
6. Chen Z, Li Y, Chen E, et al. Crystal structure at 1.9-Å resolution of human immunodeficiency virus (HIV) II protease complexed with L-735,524, an orally bioavailable inhibitor of the HIV proteases. J Biol Chem 1994;269:26344–26348. (PDB 1HSG; the deposited molecule is HIV-1 protease.)
7. Case DA, et al. AmberTools. J Chem Inf Model 2023;63:6183–6191. doi:10.1021/acs.jcim.3c01153