ScientifyStep 1
The ligand and the protein must belong to the same force field family. The three common routes are below. GROMACS 2026 ships amber14sb.ff and amber19sb.ff; the CHARMM36m port for GROMACS is downloaded from the MacKerell lab.
grompp does not complain when force fields are mixed
The CHARMM port's [ defaults ] uses fudgeLJ 1.0 and fudgeQQ 1.0; AMBER and GAFF use 0.5 and 0.8333. If a GAFF ligand is placed into a CHARMM system, grompp still writes a tpr, but the ligand's 1-4 interactions are scaled by CHARMM rules and its internal energy is wrong. Tested on this page (GROMACS 2026.3): with an ACPYPE GAFF2 ligand inside a charmm36-feb2026 system, grompp gives no WARNING at all; the log only shows Generating 1-4 interactions: fudge = 1.
The CGenFF version must match the port version
The feb2026 port comes as two packages, CGenFF 4.6 and 5.0, with 5.0 as the default. Download the port that matches the CGenFF program version the server used for your ligand. A ligand generated with 5.0 combined with the jul2022 port (4.6 only) gives Atomtype not found or missing parameters.
Always pass -n to ACPYPE
Without -n, ACPYPE guesses the net charge by summing antechamber Gasteiger charges. Charged ligands such as carboxylates or quaternary amines are easily guessed wrong, and the AM1-BCC calculation itself does not flag it.
AM1-BCC or ABCG2
ABCG2 lowers the GAFF2 hydration free energy RMSE on FreeSolv (642 molecules) to 0.99 kcal/mol. In a relative binding free energy test with 12 targets and 507 perturbations, however, the ΔΔG RMSE was 1.38 kcal/mol for ABCG2 and 1.31 kcal/mol for AM1-BCC, a non-significant difference. AM1-BCC remains a safe default for binding simulations.
| Route | Protein force field | Ligand tool and charges | Quality indicator | Water model | Typical use |
|---|---|---|---|---|---|
| CHARMM | CHARMM36m (charmm36-feb2026_cgenff-5.0.ff) | CGenFF server (cgenff.com) exports GROMACS format directly; charges assigned by analogy | Penalty: <10 usable as is, 10–50 basic validation, >50 needs reoptimization | CHARMM-modified TIP3P (tip3p in the port) | Following the Lemkul tutorial; membrane proteins or existing CHARMM systems |
| AMBER + GAFF2 | AMBER14SB or AMBER19SB (shipped with GROMACS 2026) | ACPYPE driving antechamber; AM1-BCC by default (-c bcc), ABCG2 available (-c abcg2) | Parameters filled in by parmchk2 (entries marked ATTN in the frcmod); no single numeric threshold | TIP3P for AMBER14SB; OPC or OPC3 for AMBER19SB | Quick setup after docking; later gmx_MMPBSA |
| AMBER + OpenFF | AMBER14SB | OpenFF Toolkit + Interchange; Sage 2.3.0 uses NAGL graph-neural-network charges (model file openff-gnn-am1bcc-1.0.0.pt in the offxml), Sage 2.2 and earlier use AM1-BCC | No penalty; coverage of elements and chemistries in Sage | TIP3P (Sage is co-optimized with TIP3P) | Drug-like molecules well covered by Sage training data; sharing parameters with OpenMM workflows |
Step 2
A Vina PDBQT has no bond orders and no hydrogens bonded to carbon. Converting the PDBQT to mol2 directly with OpenBabel infers bond orders from geometry; the example in the Meeko documentation shows a spurious double bond added to a ring.
If the ligand was prepared with Meeko before docking, the PDBQT REMARK lines store a SMILES string and an atom mapping, and mk_export.py uses them to restore bond orders and all hydrogens. Hydrogens on carbon are rebuilt by RDKit with simple geometric rules, and the later energy minimization corrects them.
The ligand protonation state is fixed before docking, and neither Meeko nor Vina changes it. Use the same protonation state in MD as in docking. If it has to change (for example a deprotonated carboxylic acid at pH 7.4), change the SMILES and redock, or at least rebuild the hydrogens from the new SMILES template.
When assembling the complex, take the ligand coordinates from the file written by the parameterization tool (LIG_GMX.gro from ACPYPE, the pdb returned by the CGenFF server). Their atom order and names match the topology one to one, and the coordinates come from the docked pose you supplied. Pasting coordinates from the docking output almost always breaks the atom order.
# Ligand prepared with Meeko before docking: restore bond orders and all H from the PDBQT REMARKs
mk_export.py vina_out.pdbqt -s lig_docked.sdf # older releases use -o; check mk_export.py --help
# Keep only the pose you want (first record = top-ranked pose), then write mol2 for CGenFF / ACPYPE
obabel lig_docked.sdf -O lig.mol2 -l 1
# obabel names the residue UNL1; rename it so ACPYPE writes residue LIG (make_ndx, gmx_MMPBSA and the analysis commands select "LIG")
sed -i.bak 's/UNL1/LIG1/' lig.mol2
# Check the net charge before parameterization
python3 -c "from rdkit import Chem; m=Chem.MolFromMolFile('lig_docked.sdf',removeHs=False); print(Chem.GetFormalCharge(m), m.GetNumAtoms())"# Ligand NOT prepared with Meeko: rebuild bond orders from a SMILES template, keep the docked coordinates
from rdkit import Chem
from rdkit.Chem import AllChem
pose = Chem.MolFromPDBFile("lig_pose.pdb", removeHs=True) # heavy atoms of the docked pose
template = Chem.MolFromSmiles("CCCc1ccccc1O") # same protonation state you want in MD
pose = AllChem.AssignBondOrdersFromTemplate(template, pose)
pose_h = Chem.AddHs(pose, addCoords=True) # H placed on the docked heavy atoms
Chem.MolToMolFile(pose_h, "lig_docked.sdf")
print("formal charge:", Chem.GetFormalCharge(pose_h))
# Then, as above: obabel lig_docked.sdf -O lig.mol2 && sed -i.bak 's/UNL1/LIG1/' lig.mol2# Protein only: remove ligand, crystal waters and additives first
grep -v -e HETATM -e CONECT complex_docked.pdb > protein.pdb
# CHARMM36m route (force field directory unpacked in the working directory)
# -ter: choose NH3+ and COO-. Without it an N-terminal MET gets MET1 from ethers.n.tdb and pdb2gmx stops
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top -ff charmm36-feb2026_cgenff-5.0 -water tip3p -ignh -his -ter
# AMBER route (amber14sb.ff ships with GROMACS 2026)
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top -ff amber14sb -water tip3p -ignh -his- By default pdb2gmx assigns histidine protonation from hydrogen-bond geometry. For His, Asp and Glu near the binding site, compute pKa values with a tool such as PROPKA and set the states manually with -his and the related interactive options.
- -ignh makes pdb2gmx ignore hydrogens in the input and rebuild them from the force field, which avoids hydrogen names from docking software that do not match the rtp.
- In GROMACS 2026, pdb2gmx -rtpres defaults to auto and renames residues when needed so that grompp assigns CMAP terms correctly.
- If the protein has missing residues or side-chain atoms, pdb2gmx reports Long bonds and/or missing atoms; complete the structure first.
Step 3
Merging follows three rules: the whole system has exactly one [ defaults ]; the ligand [ atomtypes ] go after forcefield.itp and before any [ moleculetype ]; a position restraint file follows directly after the moleculetype it belongs to.
The lig_gmx.top exported by the CGenFF server is a standalone system topology. Convert it to an itp as in the current Lemkul tutorial: remove the forcefield.itp include, paste in the content of lig_ffbonded.itp, delete the water and ion sections, and rename the moleculetype from Other to LIG. The ligand's new dihedral parameters sit in [ dihedraltypes ], so this itp also goes after forcefield.itp and before the protein moleculetype.
Atom numbers written by genrestr are relative to the moleculetype. Running it on the full complex gro gives system-wide numbers, and grompp reports Atom index (n) in position_restraints out of bounds. ACPYPE already writes posre_LIG.itp, which can be used directly.
With multiple chains, pdb2gmx writes the protein moleculetypes to topol_Protein_chain_A.itp and similar files, and topol.top only includes them. Put the ligand include after those includes and before the water model include.
The [ atomtypes ] written by ACPYPE have no atomic-number column, so grompp records atomnumber -1 for the ligand atoms. Forces are not affected, but gmx hbond (GROMACS 2024 and later) identifies donors and acceptors by element and stops on the ligand with Selection 'resname LIG' has no donors AND has no acceptors! Nothing to be done. (tested on this page with GROMACS 2026.3). The second awk in the AMBER route code fills in the atomic number from the mass; with it, the same command reports the hydrogen bond between the ligand and Gln102. The CGenFF port and OpenFF Interchange write atomtypes with atomic numbers, so they do not need this step.
# GAFF2 + AM1-BCC; always pass the net charge with -n
acpype -i lig.mol2 -b LIG -c bcc -a gaff2 -n 0
cd LIG.acpype
# LIG_GMX.itp starts with [ atomtypes ]; split it so atomtypes can go before every [ moleculetype ]
# ACPYPE atomtypes carry no atomic number; the second awk adds it (from the mass) so that
# gmx hbond (2024+) can find the ligand's N/O donors and acceptors
awk '/^\[ *atomtypes *\]/{f=1} /^\[ *moleculetype *\]/{f=0} f' LIG_GMX.itp | \
awk '!/^ *[;[]/ && NF>=7 {m=$3; z=(m<1.5)?1:(m<13)?6:(m<15)?7:(m<17)?8:(m<20)?9:(m<31.5)?15:(m<33)?16:(m<36)?17:(m<80)?35:53; $2=$2" "z} 1' \
> ../LIG_atomtypes.itp
awk '/^\[ *moleculetype *\]/{f=1} f' LIG_GMX.itp > ../LIG.itp
cp LIG_GMX.gro ../LIG.gro
cp posre_LIG.itp ../posre_LIG.itp # heavy-atom restraints, 1000 kJ/mol/nm^2
cd ..# OpenFF Sage ligand -> GROMACS files (Interchange)
from openff.toolkit import ForceField, Molecule
from openff.interchange import Interchange
lig = Molecule.from_file("lig_docked.sdf") # coordinates = docked pose, H included
lig.name = "LIG"
sage = ForceField("openff-2.3.0.offxml")
inter = Interchange.from_smirnoff(force_field=sage, topology=[lig], box=[5, 5, 5])
inter.to_top("LIG_openff.top")
inter.to_gro("LIG_openff.gro")
# Then split it in bash; [ defaults ] is dropped, and the moleculetype part stops at [ system ]:
# awk '/^\[ *atomtypes *\]/{f=1} /^\[ *moleculetype *\]/{f=0} f' LIG_openff.top > LIG_atomtypes.itp
# awk '/^\[ *moleculetype *\]/{f=1} /^\[ *system *\]/{f=0} f' LIG_openff.top > LIG.itp
# cp LIG_openff.gro LIG.gro # atom order = lig_docked.sdf; make posre_LIG.itp with genrestr as in the CGenFF route# CGenFF route: generate ligand heavy-atom restraints on the ligand-only coordinates
gmx editconf -f lig_gmx.pdb -o LIG.gro
printf '0 & ! a H*\nq\n' | gmx make_ndx -f LIG.gro -o index_LIG.ndx
echo 3 | gmx genrestr -f LIG.gro -n index_LIG.ndx -o posre_LIG.itp -fc 1000 1000 1000# Append ligand atoms to the protein coordinates (same frame, nm)
python3 - <<'EOF'
p = open("protein.gro").read().splitlines()
l = open("LIG.gro").read().splitlines()
atoms = p[2:-1] + l[2:-1]
open("complex.gro", "w").write("\n".join([p[0], f"{len(atoms):5d}", *atoms, p[-1]]) + "\n")
EOF; topol.top (AMBER + GAFF2 example; CHARMM route is identical in layout)
#include "amber14sb.ff/forcefield.itp"
#include "LIG_atomtypes.itp" ; ligand atomtypes: after forcefield.itp, before any [ moleculetype ]
[ moleculetype ]
; Name nrexcl
Protein_chain_A 3
... ; written by pdb2gmx
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Ligand topology and its restraints (must follow the ligand [ moleculetype ])
#include "LIG.itp"
#ifdef POSRES_LIG
#include "posre_LIG.itp"
#endif
; Include water topology
#include "amber14sb.ff/tip3p.itp"
...
[ molecules ]
; Compound #mols (same order as complex.gro)
Protein_chain_A 1
LIG 1Step 4
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 # OPC/TIP4P-type water: -cs tip4p.gro
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
echo SOL | gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15; ions.mdp (only used to build ions.tpr for genion)
integrator = steep
emtol = 1000.0
emstep = 0.01
nsteps = 50000
nstlist = 1
cutoff-scheme = Verlet
coulombtype = cutoff ; plain cutoff: no "Ewald with net charge" warning before neutralization
rcoulomb = 1.0
rvdw = 1.0
pbc = xyzBox size
-d 1.0 keeps the complex at least 2.0 nm from its periodic image, more than the 1.2 nm cutoff, so the complex does not interact directly with itself. The strict condition in the GROMACS manual is a box length of at least the molecule size plus twice the cutoff (equivalent to -d 1.2 with CHARMM settings); relaxing it to -d 1.0 to reduce the number of waters is common practice. A rhombic dodecahedron has 71% of the volume of a cube with the same image distance.
Water box file
Three-site models such as TIP3P and SPC use spc216.gro; OPC with AMBER19SB uses the four-site box tip4p.gro.
Ion names
Both the AMBER ports shipped with GROMACS and the CHARMM port define NA and CL; in the CHARMM port NA/SOD and CL/CLA use the same atom type. -conc 0.15 adds salt up to 0.15 mol/L on top of neutralization.
Wrong group in genion
When prompted for the group to replace, choose SOL. Choosing Water or System gives The solvent group ... is not continuous if the atoms in the group are not contiguous.
Step 5
The parameters below are written for CHARMM36m and match the CHARMM36 settings in the GROMACS manual and the current Lemkul tutorial. For the AMBER route, replace only the nonbonded cutoff block.
; em.mdp (CHARMM36m settings; see the AMBER note below)
integrator = steep
emtol = 1000.0 ; kJ/mol/nm
emstep = 0.01
nsteps = 50000
nstlist = 1
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = no
pbc = xyz; nvt.mdp 100 ps, protein and ligand heavy atoms restrained
define = -DPOSRES -DPOSRES_LIG
integrator = md
nsteps = 50000 ; 50000 * 2 fs = 100 ps
dt = 0.002
nstxout-compressed = 5000
nstenergy = 500
nstlog = 500
continuation = no
constraint-algorithm = lincs
constraints = h-bonds
lincs-iter = 1
lincs-order = 4
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.2
fourierspacing = 0.16
pme-order = 4
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = no
tcoupl = V-rescale
tc-grps = System
tau-t = 1.0
ref-t = 300
pcoupl = no
pbc = xyz
gen-vel = yes
gen-temp = 300
gen-seed = -1 ; new seed per replica; npt.mdp 100 ps, restraints kept, C-rescale barostat
define = -DPOSRES -DPOSRES_LIG
integrator = md
nsteps = 50000
dt = 0.002
nstxout-compressed = 5000
nstenergy = 500
nstlog = 500
continuation = yes
constraint-algorithm = lincs
constraints = h-bonds
lincs-iter = 1
lincs-order = 4
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.2
fourierspacing = 0.16
pme-order = 4
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = no
tcoupl = V-rescale
tc-grps = System
tau-t = 1.0
ref-t = 300
pcoupl = C-rescale
pcoupltype = isotropic
tau-p = 5.0
ref-p = 1.0
compressibility = 4.5e-5
refcoord-scaling = com ; required with position restraints + pressure coupling
pbc = xyz
gen-vel = no; md.mdp 100 ns production, no restraints
integrator = md
nsteps = 50000000 ; 50,000,000 * 2 fs = 100 ns
dt = 0.002
nstxout = 0
nstvout = 0
nstfout = 0
nstxout-compressed = 5000 ; one frame every 10 ps
compressed-x-grps = System
nstenergy = 5000
nstlog = 5000
continuation = yes
constraint-algorithm = lincs
constraints = h-bonds
lincs-iter = 1
lincs-order = 4
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.2
fourierspacing = 0.16
pme-order = 4
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = no
tcoupl = V-rescale
tc-grps = System
tau-t = 1.0
ref-t = 300
pcoupl = C-rescale
pcoupltype = isotropic
tau-p = 5.0
ref-p = 1.0
compressibility = 4.5e-5
pbc = xyz
gen-vel = no; AMBER14SB / AMBER19SB + GAFF2 or Sage: replace the nonbonded block with
rcoulomb = 1.0
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 1.0
DispCorr = EnerPres
; and delete rvdw-switch- There is no official cutoff recommendation for AMBER force fields in GROMACS. rvdw = rcoulomb = 1.0 nm, Potential-shift and DispCorr = EnerPres are the common community settings on the GROMACS forum for AMBER14SB and AMBER19SB.
- constraints = h-bonds goes with dt = 2 fs. GPU-resident update (-update gpu) requires constraints on h-bonds only.
- C-rescale has a default tau-p of 5 ps and can be used for both equilibration and production; switching to Parrinello-Rahman for production, as old tutorials do, is no longer necessary.
- refcoord-scaling = com is required in NPT with position restraints; otherwise grompp warns You are using pressure coupling with absolute position restraints.
- NVT and NPT of 100 ps each is the length used in the Lemkul tutorial. The criterion is that temperature, pressure and density are flat in the second half; check Temperature, Pressure and Density with gmx energy.
- A 100 ns xtc with one frame every 10 ps has about 10,000 frames. If only the complex matters, set compressed-x-grps to Protein_LIG (requires -n index.ndx) to make the trajectory much smaller.
Version differences
Most second-hand tutorials repeat the 2018 version of the Lemkul tutorial. The table lists where it differs from GROMACS 2026 and the current Lemkul tutorial.
The 2018 tutorial puts the ligand in the same thermostat group as the protein because groups with few atoms, such as a ligand or ions, have large kinetic energy fluctuations and make the thermostat unstable when coupled separately. With V-rescale, a single System group avoids the issue.
# Optional: two coupling groups (the 2018 tutorial layout)
printf '"Protein" | "LIG"\nq\n' | gmx make_ndx -f em.gro -o index.ndx
# in nvt/npt/md.mdp: tc-grps = Protein_LIG Water_and_ions tau-t = 1.0 1.0 ref-t = 300 300
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr| Old setting | Result in GROMACS 2026 | Current setting |
|---|---|---|
| pcoupl = Berendsen (NPT equilibration) | grompp issues a WARNING; -maxwarn is needed to continue | pcoupl = C-rescale, tau-p = 5.0 |
| tcoupl = Berendsen | Also issues a WARNING | tcoupl = V-rescale |
| tc-grps = Protein_JZ4 Water_and_ions | Works, but needs make_ndx and -n in grompp | tc-grps = System (current Lemkul tutorial); two groups are also correct, but never put the ligand or the ions in a group of their own |
| charmm36-jul2022.ff + cgenff_charmm2gmx.py | Contains only CGenFF 4.6 parameters; mismatches ligands generated with 5.0 | charmm36-feb2026_cgenff-5.0.ff, with GROMACS format exported by the server |
| ns_type = grid | Ignored, prints Ignoring obsolete mdp entry | Delete the line |
| Misspelled or removed mdp parameters | Unknown left-hand ... in parameter file, counted as one WARNING | Rename according to the current mdp documentation |
| gmx mdrun -nsteps to change the step count | Still works, but discouraged since 2019 | gmx convert-tpr -extend / -until / -nsteps |
| mdrun -deffnm | Still available in 2026, marked for removal since 2021 | Can still be used; a restart must use the same -deffnm as the first run |
Step 6
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt -nb gpu -pme gpu
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt -nb gpu -pme gpu
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 8 -nb gpu -pme gpu -bonded gpu -update gpu -pin on# Three independent replicas: same minimized structure, new velocities in NVT
for r in 1 2 3; do
mkdir -p rep$r && cd rep$r
gmx grompp -f ../nvt.mdp -c ../em.gro -r ../em.gro -p ../topol.top -o nvt.tpr
gmx mdrun -deffnm nvt -nb gpu -pme gpu
gmx grompp -f ../npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p ../topol.top -o npt.tpr
gmx mdrun -deffnm npt -nb gpu -pme gpu
gmx grompp -f ../md.mdp -c npt.gro -t npt.cpt -p ../topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -nb gpu -pme gpu -bonded gpu -update gpu
cd ..
done- Energy minimization is acceptable when the log shows converged to Fmax < 1000 and the potential energy is negative; in the Lemkul tutorial a simple protein in water is on the order of 10^5–10^6 kJ/mol, depending on system size and number of waters.
- -r em.gro provides the reference coordinates for position restraints and has been mandatory since GROMACS 2018.
- -update gpu is not set explicitly for NVT and NPT. mdrun defaults to -update auto and falls back to the CPU for unsupported features; an explicit -update gpu fails with an error in that case.
- On a single GPU use -ntmpi 1 and set -ntomp to the number of physical cores assigned to the job. With PME on the GPU, -npme can only be 1.
- In GPU-resident runs, frequent energy calculation and coupling slow things down; keep nstcalcenergy at the default of 100 or larger.
- Experience figures: for a protein-in-water system of about 60,000 atoms with all work offloaded (-nb -pme -bonded -update gpu), one author measured about 80, 150 and 250 ns/day on T4, V100 and A100, and about 20 ns/day on a 24-core CPU (Tencent Cloud experience post, described by the author as a rough test).
- Tested on this page (CPU, for reference only; a GPU is much faster): T4 lysozyme L99A/M102Q with 2-propylphenol (PDB 3HTB), AMBER14SB + GAFF2 (AM1-BCC), rhombic dodecahedron -d 1.0, 0.15 M NaCl, 33,379 atoms in total; conda-forge GROMACS 2026.3 CPU build, Apple M2 with 8 cores, 1 thread-MPI rank × 8 OpenMP threads, with other jobs running on the machine at the same time. ACPYPE took 31 s for the AM1-BCC charges; energy minimization reached Fmax < 1000 in 955 steps (potential energy −5.34×10^5 kJ/mol) in about 8 min; the 50 ps production segment ran at 4.8 ns/day under heavy load; at lower load (1-minute load average about 13), mdrun -nsteps 5000 -resethway gave 10.3 ns/day, so 100 ns would take about 10 days at this speed.
-maxwarn
-maxwarn N lets grompp write a tpr as long as there are no more than N WARNINGs. NOTEs are not counted. The grompp help states that the option is not for normal use and may generate unstable systems. Decide by reading each warning.
| Warning text (excerpt) | Meaning | What to do |
|---|---|---|
| You are using Ewald electrostatics in a system with net charge | PME is used before the system is neutralized | When it appears only while building ions.tpr, -maxwarn 1 is acceptable because genion neutralizes next; an ions.mdp with coulombtype = cutoff avoids it. If it still appears from EM onward, the system is not neutral or the ligand charges do not sum to an integer; fix the topology |
| N non-matching atom names ... atom names from topol.top will be used | Atom order or names in the coordinates do not match the topology | Never skip. Coordinates would receive the parameters of the wrong atoms. Reassemble with the ligand coordinates written by the parameterization tool |
| Atomtype X was defined previously ... and has now been defined again | The ligand itp redefines an atom type that already exists in the force field | Never skip. The later definition overrides the earlier one and may change protein force field parameters. Remove the duplicate or rename the ligand atom type |
| You are using pressure coupling with absolute position restraints | NPT with restraints but no refcoord-scaling | Do not skip; add refcoord-scaling = com |
| The Berendsen thermostat does not generate the correct kinetic energy distribution; The Berendsen barostat does not generate any strictly correct ensemble | Berendsen is used (tested on this page; each counts as one WARNING) | Switch to V-rescale / C-rescale |
| Unknown left-hand '...' in parameter file | Misspelled or removed mdp parameter | Fix the name; skipping means the parameter has no effect |
| System has non-zero total charge: 0.000999 | This is a NOTE and needs no -maxwarn | Check the ligand charge sum; go back to parameterization if it is more than 0.01 from an integer |
Restarts
# 1) Interrupted run: continue from the last checkpoint (written every 15 min by default)
gmx mdrun -deffnm md -cpi md.cpt -nb gpu -pme gpu -bonded gpu -update gpu
# 2) Finished 100 ns, extend by another 100 ns (-extend is in ps)
gmx convert-tpr -s md.tpr -extend 100000 -o md_200ns.tpr
gmx mdrun -deffnm md -s md_200ns.tpr -cpi md.cpt -nb gpu -pme gpu -bonded gpu -update gpu
# Or set the end time directly: run until 200 ns
gmx convert-tpr -s md.tpr -until 200000 -o md_200ns.tpr
# 3) Queue system with a wall-time limit: stop cleanly before 23.5 h and write a checkpoint
gmx mdrun -deffnm md -cpi md.cpt -maxh 23.5
# 4) Changing mdp options while continuing (e.g. output frequency)
gmx grompp -f md_new.mdp -c npt.gro -t md.cpt -p topol.top -o md_new.tpr
gmx mdrun -s md_new.tpr -cpi md.cpt -noappend| Error text (excerpt) | Cause | What to do |
|---|---|---|
| The input requested N steps, however the checkpoint file has already reached step M | The total step count in the tpr is smaller than the step already reached in the checkpoint, typically when the original md.tpr is used again after one extension. Tested on this page: restarting a run that has exactly finished with -cpi gives no error; the log shows continuing from step N and mdrun ends immediately without running a single step | Create a new tpr with convert-tpr -extend or -until and pass it with -s |
| Checksum wrong for 'md.log'. The file has been replaced or its contents have been modified | An output file was changed or replaced between runs | Restore the original file, or add -noappend to write new .part0002 files |
| Some output files listed in the checkpoint file md.cpt are not present or not named as the output files by the current program | File names differ from the first run, typically -deffnm in the first run but not in the restart | Use the same -deffnm or file names as the first run |
| Checkpoint file is for a system of X atoms, while the current system consists of Y atoms | The cpt and the tpr belong to different systems | Find the matching tpr and cpt |
| Cannot restart with appending because the previous simulation part used ... precision | Single- and double-precision builds of GROMACS were mixed | Use the same build, or add -noappend |
- -extend and -until are in ps. To add 100 ns, write -extend 100000; a user on the GROMACS forum wrote -extend 5000000 thinking in steps and ended up with 5 μs.
- A checkpoint is written every 15 minutes by default (-cpt 15), so an unexpected stop loses at most about 15 minutes.
- When changing mdp parameters during a restart, read full-precision coordinates and velocities with grompp -t md.cpt, and keep continuation = yes and gen-vel = no.
- Do not restart from the last gro frame instead of the cpt. A gro file holds coordinates with three decimals and no thermostat or barostat state.
- With -noappend, the number in the new file names is the simulation part: the first restart gives md.part0002.*, and it increases from there. Tested on this page: part 4 of the same system run with -noappend produced md.part0004.xtc, and gmx trjcat -f md.xtc md.part0004.xtc joined them directly.
Errors
| Error text (excerpt) | Actual cause | What to do |
|---|---|---|
| Residue 'LIG' not found in residue topology database | A pdb containing the ligand was given to pdb2gmx | Run pdb2gmx on the protein only and parameterize the ligand separately |
| Atom HB3 in residue XXX not found in rtp entry | Hydrogen names in the input do not match the force field rtp | Add -ignh so pdb2gmx rebuilds hydrogens |
| Invalid order for directive atomtypes | The ligand [ atomtypes ] appear after the protein moleculetype | Split the atomtypes into their own file and include it right after forcefield.itp |
| Found a second defaults directive | The ligand top/itp contains [ defaults ] | Delete [ defaults ] from the ligand file |
| Atomtype XXX not found | Ligand atomtypes not included, or CGenFF version does not match the port | Check the include order; use the CHARMM port that matches the CGenFF version |
| No such moleculetype LIG | The name in [ molecules ] differs from the ligand moleculetype name | Use the same name in both places |
| number of coordinates in coordinate file (solv.gro, N) does not match topology (topol.top, M) | The ligand is missing from [ molecules ], or counts or order differ from the gro | Check each molecule's count in the order of the gro file |
| Atom index (n) in position_restraints out of bounds (1-m) | Ligand restraints placed under the wrong moleculetype, or genrestr run on the complex | Include the file right after the ligand itp; rerun genrestr on the ligand-only gro |
| Steepest Descents converged to machine precision ... but did not reach the requested Fmax < 1000 | Overlapping atoms or abnormal ligand geometry | Read Maximum force = ... on atom N in the log to locate the atom; if the docked pose clashes with the protein or hydrogens are misplaced, minimize the ligand alone in vacuum first and inspect it |
| LINCS WARNING ... relative constraint deviation | The system blows up early in NVT, usually because of the ligand topology or insufficient minimization | Follow the diagnosis steps in the GROMACS documentation: find which atoms become unstable first; minimize the protein in water, the ligand in vacuum and the ligand in water separately; use gmx energy to find the bonded term that spikes |
| atom C1 not found in buiding block 1MET while combining tdb and rtp | ethers.n.tdb in the CHARMM36 port defines a terminus named MET1, and pdb2gmx selects it by default for a protein whose N-terminal residue is Met (tested on this page with both the jul2022 and feb2026 ports on GROMACS 2026.3) | Add -ter to pdb2gmx and choose NH3+ for the N terminus and COO- for the C terminus |
| step N: One or more water molecules can not be settled (early in energy minimization, with stepNb.pdb and stepNc.pdb written) | One steepest-descent step was too large for the rigid water constraints; the step is rejected and the algorithm reduces the step size | Ignore it if minimization ends with converged to Fmax < 1000, and delete step*.pdb; if it appears in NVT or production, treat it like a LINCS error |
Length and replicas
Number of replicas
The MD reliability checklist published in Communications Biology in 2023 asks for at least 3 independent simulations per condition with statistical analysis. Knapp et al. (JCTC 2018) tested 100 replicas, suggested 5–10 as a rule-of-thumb minimum, and found that several shorter replicas give more reliable conclusions than one long trajectory.
How replicas are generated
Start from the same minimized structure and regenerate velocities in NVT with gen-seed = -1; each replica goes through its own NVT, NPT and production run. Changing only the final production segment does not give independent replicas.
How to judge the length
Use the convergence of the quantity you will report: whether the ligand RMSD relative to the protein and key hydrogen-bond distances are flat in the second half, and whether the means of the three replicas agree. The 10 ns in the Lemkul tutorial is for demonstration only.
When the ligand drifts away
If the ligand RMSD keeps rising in one replica and the ligand leaves the pocket, that is part of the result; report it together with the other replicas and give the fraction of replicas in which the ligand stays in the pocket.
Analysis
Commands and interpretation for RMSD, RMSF, radius of gyration, hydrogen bonds and SASA are in the trajectory analysis guide. This section covers only the preprocessing required before MM/GBSA and the changes in gmx_MMPBSA 1.7.
# Make molecules whole and keep the complex centered
printf '"Protein" | "LIG"\nq\n' | gmx make_ndx -f em.gro -o index.ndx
printf 'Protein_LIG\nSystem\n' | gmx trjconv -s md.tpr -f md.xtc -n index.ndx \
-o md_center.xtc -center -pbc mol -ur compact
# gmx_MMPBSA 1.7.x, GB with the 1.7 defaults written out explicitly
cat > mmpbsa.in <<'EOF'
&general
sys_name="complex", startframe=5001, endframe=10000, interval=10,
PBRadii=4,
/
&gb
igb=8, saltcon=0.150,
/
EOF
gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ct md_center.xtc -ci index.ndx \
-cg Protein LIG -cp topol.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv- gmx_MMPBSA 1.7.0 was released on 2026-09-12. GROMACS calculations require -cp topol.top, and the ligand topology must already be in topol.top.
- 1.7.0 changed defaults: when not set, igb goes from 5 to 8, PBRadii from 3 to 4 (mbondi3), and exdi from 80 to 78.5. Write these values explicitly in the input before comparing with older results.
- qh_entropy = 1 is rejected in 1.7.0. Use interaction entropy or C2 entropy when an entropy correction is needed.
- Tested ranges: GROMACS 2022–2026, AmberTools ≥24.8 and <27, Python 3.11–3.12.
- Absolute MM/GBSA values usually deviate substantially from experimental binding free energies. The method is suited to comparing structurally similar ligands on the same target, using the same protocol and number of frames for every ligand.
- Neighboring frames are highly correlated, so doubling the number of frames does not halve the error. The example uses every 10th frame (100 ps) and estimates the error from the spread between the three replicas.
- Tested on this page: gmx_MMPBSA 1.7.0 (installed with pip, AmberTools 26, GROMACS 2026.3) ran on the 3HTB system above with -cg Protein LIG (-cg accepts group names or zero-based group numbers); a 6-frame GB calculation took about 70 s, and ΔTOTAL appears at the end of FINAL_RESULTS_MMPBSA.dat. This is only a connectivity test on a 50 ps trajectory, and the number has no physical meaning.
Practice in China
The AMBER route most often seen on the Keinsci (Computational Chemistry Commune) forum and on Sobereva's blog is: compute RESP or RESP2 charges with the scripts shipped with Multiwfn, then build a GAFF topology with Sobtop. The force field is the same as in the ACPYPE route above; the charge method and the tools differ. The content below comes from the Sobtop home page, Sobereva's posts and forum answers; commands were checked against Sobtop 2026.1.16 and the scripts shipped with Multiwfn.
# 1) Charges: scripts shipped with Multiwfn (examples/RESP/ in the Multiwfn directory), run on Linux
# Before running, set the ORCA= and orca_2mkl= paths, nprocs and maxcore inside the script
# Arguments: structure file, net charge, spin multiplicity; input is the docked ligand with H,
# the script first optimizes it at B97-3c and then fits the charges
./RESP2_ORCA.sh lig.mol2 0 1
# If the structure is already optimized, use the version without optimization
./RESP2_ORCA_noopt.sh lig.mol2 0 1
# With Gaussian, use RESP2.sh or RESP.sh the same way
# Output lig.chg: element and coordinates in the first four columns, charge in the last;
# atom order is the same as in the input file
# 2) Topology: start ./sobtop and enter in turn (GAFF route, example 2 on the Sobtop page)
# lig.mol2 docked pose, with H and correct bond orders
# lig.chg type the chg path at the main menu to load the RESP2 charges
# 2, [Enter] write the gro file; coordinates come from the mol2
# 1, 2, 4 write a GROMACS topology; assign GAFF atom types; take parameters from
# the library and let the program guess missing ones
# [Enter], [Enter] default top and itp pathsHow Sobtop output goes into topol.top
Sobtop writes lig.itp and lig.top. lig.top contains only [ defaults ] (fudgeLJ 0.5, fudgeQQ 0.8333, the same as AMBER) and an include, and is not used when merging. lig.itp starts with [ atomtypes ]; split it out by the rules in step 3 and place it after forcefield.itp. Sobtop atomtypes include an atomic-number column (at.num), so the atomic-number fix needed in the ACPYPE route is not required. The residue name defaults to MOL and the molecule name is taken from the file name; the commands on this page select LIG, so either rename it to LIG in the itp and gro or use MOL in the commands. Checked against the Methyl_benzoate example files in the Sobtop 2026.1.16 package.
Without a charge column in the mol2, all charges are 0
Sobtop writes the atomic charges recorded in the mol2 into the itp; if the mol2 has no charges it writes 0 for every atom and reports no error. A mol2 exported by docking software may carry Gasteiger charges, and these are used as is. After generating the topology, check that the charge column in [ atoms ] really comes from the chg file and sums to the net charge.
Option 4 guesses missing parameters
With option 4, bonds and angles missing from the parameter library take the current structure as the equilibrium value with approximate force constants, and missing dihedrals get a rotational barrier of 0. Ordinary organic ligands are rarely missing parameters; when the screen lists missing terms, read them one by one, and handle missing rotatable dihedrals separately. For more accurate bond and angle parameters, choose option 7: bonds and angles are derived from the Hessian in a Gaussian fchk or ORCA hess file, while dihedrals still come from GAFF (example 3 on the Sobtop page).
Level of theory for RESP
RESP.sh by default optimizes at B3LYP-D3(BJ)/def2-SVP and runs the single point at B3LYP-D3(BJ)/def2-TZVP with IEFPCM water as the default solvent; the ORCA scripts optimize at B97-3c, and differences of a few hundredths between the two are normal. Sobereva noted in a forum answer that 6-311G** is fine for optimization and frequencies but on the low side for RESP fitting in Multiwfn, and that diffuse functions on hydrogen are unnecessary. For elements 18 and beyond (for example Br and I) Gaussian has no built-in fitting radii, so the scripts cannot finish automatically and the charges must be computed by hand as described in sobereva.com/441. For simulations in solution, Sobereva recommends RESP2(0.5).
Which conformation to use for ligand charges
Sobereva's advice (sobereva.com/441, appendix 2): a docked ligand conformation may be unreasonable and should not be used directly for RESP. Use the docked conformation as the starting guess for a regular geometry optimization, compute RESP, then run the complex MD. If the dominant ligand conformation in the MD differs markedly from the optimized one, cluster the trajectory, take the representative structure, minimize it with the force field, extract the ligand and compute a single-point RESP (no further quantum-chemical optimization), then use the new charges for the production simulations. For a ligand from a high-resolution crystal structure, optimize only the hydrogen positions before computing RESP.
AM1-BCC, ABCG2 and RESP: two views on the forum
Sobereva's view is that ABCG2 only replaces AM1-BCC and cannot replace RESP or RESP2. A member of the group that developed ABCG2 replied in the same thread that ABCG2 was tuned specifically for GAFF2, outperforms RESP on solvation free energies, and that both GAFF2-ABCG2 and GAFF2-RESP(2) should be regarded as valid combinations. Another thread notes that for molecules with P=O (phosphate esters, FAD and the like), when antechamber calls sqm for AM1-BCC a hydrogen may migrate onto the P=O during optimization and the charges come out wrong; for such molecules compute RESP with Multiwfn instead.
gmx_MMPBSA, gmx_mmpbsa and g_mmpbsa are three different tools
Chinese materials often mix them up. gmx_MMPBSA, used in the analysis section of this page, is the Python program by Valdés-Tresanco et al. built on MMPBSA.py from AmberTools. gmx_mmpbsa is a bash script by Jerkwin (Jicun Li) that reads parameters from the tpr with gmx dump and calls APBS for the PB term; the 2019 version had PB only, without GB or entropy, and a 2021 update added screening effects and the entropy contribution. g_mmpbsa supports only specific versions of GROMACS and APBS; it and GMXPBSAtool are what section 9 of Jerkwin's Chinese GROMACS tutorial describes. Before copying commands, check which tool a tutorial uses; their inputs and outputs are not interchangeable.
Hand it to an agent
Example instruction: "Take the top-ranked pose in docking/vina_out.pdbqt and receptor.pdb, build the system with AMBER14SB + GAFF2 (AM1-BCC, ligand net charge 0) in 0.15 M NaCl, run 3 replicas of 100 ns, run gmx_MMPBSA on the last 50 ns, and check whether the ligand stays in the pocket."
The agent exports the protonated ligand, runs ACPYPE and pdb2gmx, merges the topology, solvates and adds ions, minimizes and equilibrates, rents a GPU when needed for production, restarts from checkpoints after interruptions, and then runs PBC treatment, RMSD and hydrogen-bond analysis and gmx_MMPBSA. The outputs include topol.top and all itp files, the four mdp files, tpr/xtc/edr/log for each replica, analysis plots and FINAL_RESULTS_MMPBSA.dat. The Scientify cases "Lysozyme molecular dynamics simulation" and "Small-molecule hydration free energy" use the same environment.
You still need to judge whether the ligand protonation state and net charge are correct, whether key residues such as His are protonated sensibly, whether the docked pose itself is credible, and how the MM/GBSA numbers should be interpreted in a paper.
References
- GROMACS 2026 mdp options — Official description of Berendsen, C-rescale, constraints and refcoord-scaling
- GROMACS 2026 new features: AMBER14SB/AMBER19SB and pdb2gmx -rtpres — Shipped AMBER force fields and validation-pending status
- GROMACS force fields — Recommended CHARMM36 settings; OPC/OPC3 with AMBER19SB
- GROMACS source readir.cpp, grompp.cpp, toppush.cpp, topio.cpp (release-2026) — Exact text of each grompp warning and whether it is a WARNING or a NOTE
- GROMACS: Managing long simulations — -cpi, appending and checksums, extending with convert-tpr
- gmx convert-tpr help — -extend/-until are in ps
- GROMACS mdrun performance guide — GPU-resident runs and the -npme limit
- GROMACS reference manual: periodic boundary conditions — Rhombic dodecahedron volume and box size conditions
- Lemkul: Lysozyme in Water, energy minimization — EM potential energy magnitude and Fmax check
- Genheden & Ryde: review of MM/PBSA and MM/GBSA — Scope and limitations of MM/GBSA
- GROMACS common errors — Official explanations of topology errors
- GROMACS terminology: blowing up and diagnosing unstable systems — Steps for investigating LINCS warnings
- Lemkul: Protein-Ligand Complex tutorial (2025.x) — Current workflow, mdp files and penalty thresholds
- Lemkul: Protein-Ligand Complex tutorial (2018) — Old settings for comparison
- MacKerell lab: CHARMM36 port for GROMACS — feb2026 port and CGenFF version matching
- ACPYPE (GitHub) — Command options, ABCG2 option, net charge guessing
- He et al., ABCG2, JCTC 2025 — ABCG2 hydration free energy data
- Behera, Gapsys, de Groot, JCIM 2025 — ABCG2 versus AM1-BCC for binding free energies
- OpenFF force field releases — Sage 2.3.0 uses AshGC charges and is co-optimized with TIP3P
- OpenFF Interchange construction and export docs — from_smirnoff and GROMACS export
- Meeko: exporting docking results — PDBQT bond-order problem and mk_export.py
- gmx_MMPBSA 1.7.0 release notes — Changed defaults and the -cp requirement
- gmx_MMPBSA compatibility and migration — Supported GROMACS, AmberTools and Python versions
- Reliability and reproducibility checklist for MD simulations, Commun. Biol. 2023 — At least 3 independent simulations per condition
- Knapp, Ospina-Forero, Deane, JCTC 2018 — Rule-of-thumb minimum number of replicas
- GROMACS forum: nonbonded mdp parameters for AMBER force fields — Experience post: common community cutoff settings for AMBER
- GROMACS forum: Extend option in gmx convert-tpr — Experience post: -extend given in steps extends far too much
- gmx-users: v-rescale fatal error (ions coupled separately) — Experience post: Lemkul explains that ions should not be a separate thermostat group
- gmx-users: CGenFF validation / optimization — Experience post: interpreting and validating CGenFF penalties
- Tencent Cloud developer community: GROMACS GPU build and speed test — Experience post: rough speeds for a 60,000-atom system on T4/V100/A100
- Sobtop home page (2026.1.16, examples 2 and 3) — Sobtop menu flow, reading charges from mol2, missing-parameter handling and Hessian-derived parameters
- Sobereva: a super lazy script to calculate RESP atomic charges — Default levels, arguments and noopt variant of RESP.sh; limitation for elements 18 and beyond
- Sobereva: lazy scripts for RESP, RESP2 and 1.2*CM5 charges with ORCA and Multiwfn — RESP2_ORCA.sh settings (ORCA path, nprocs, maxcore) and B97-3c optimization
- Sobereva: principle of RESP fitting and its calculation in Multiwfn (appendix 2) — Which conformation to use for ligand RESP charges in protein-ligand systems
- Keinsci forum: Gaussian keywords for RESP charges and Sobtop parameters of a small molecule — Experience post: Sobereva on the basis set for the RESP single point and diffuse functions
- Keinsci forum: is ABCG2 with GAFF2 better than RESP/RESP2? — Experience post: two views from Sobereva and a member of the ABCG2 developer group
- Keinsci forum: a small question on computing ligand RESP charges — Experience post: AM1-BCC via sqm fails for molecules with P=O
- Jerkwin: gmx_mmpbsa user guide (Chinese) — The APBS-based gmx_mmpbsa script, similar in name to gmx_MMPBSA but a different tool