Choosing GROMACS .mdp parameters: a stage-by-stage guide for GROMACS 2026

Most .mdp files in circulation are copies of copies, and most have picked up a stale setting somewhere: a Berendsen barostat, a group cut-off scheme, or rcoulomb = 1.0 next to rvdw = 1.2. This guide works through the four stages of a standard biomolecular simulation (energy minimisation, NVT equilibration, NPT equilibration, production) and explains each parameter group and the reason for its value.

Every recommendation here was checked against GROMACS 2026.3. The four templates produced by the Build tab of the mdp checker pass gmx grompp -maxwarn 0 and run with gmx mdrun on a solvated test system (Trp-cage, PDB 1L2Y, 7,654 atoms). The AMBER-flavoured versions were also used for the protein–ligand tutorial (37,633 atoms). You can generate the files yourself in the tool and read this page alongside them.

The principle: change as little as possible

GROMACS defaults have improved a lot. In 2026 you can leave these unset, or set them explicitly to their defaults:

cutoff-scheme = Verlet (the only scheme since 2020)
verlet-buffer-tolerance = 0.005 (the pair-list buffer is sized automatically, and mdrun may raise nstlist for performance)
nstcalcenergy, nsttcouple and nstpcouple (-1 lets grompp choose values consistent with your coupling constants)
constraint-algorithm = lincs, comm-mode = Linear

The parameters you must decide are the integrator and time step, constraints, the non-bonded treatment (dictated by the force field), temperature and pressure coupling, and output frequency. These are covered below.

Stage 1: energy minimisation

MDP

integrator  = steep
emtol       = 1000.0
emstep      = 0.01
nsteps      = 50000

Steepest descent is robust for removing clashes after solvation and ion placement. emtol = 1000 kJ mol⁻¹ nm⁻¹ is the conventional target before starting dynamics with a 2 fs step. You do not need a lower value unless you are doing normal-mode analysis. Our 37,633-atom protease–inhibitor system converged in 989 steps (11 s on 8 CPU threads), finishing with Fmax = 984 kJ mol⁻¹ nm⁻¹.

What to check afterwards: the log must say converged to Fmax < 1000. If it says “converged to machine precision” with a larger Fmax, the minimiser got stuck. The log names the atom with the largest force. Inspect it, because it usually points to a structural problem rather than a parameter one.

Minimisation uses the same non-bonded settings as dynamics. Leave constraints off (constraints = none). Note that constraint-algorithm = SHAKE gets past grompp but mdrun refuses it with “Can not do energy minimization with Shake, use Lincs”.

Non-bonded settings: set by the force field, not by taste

Cut-offs and van der Waals modifiers are part of how a force field was parametrised. Use what its developers validated.

CHARMM36 (GROMACS manual recommendation):

MDP

coulombtype   = PME
rcoulomb      = 1.2
vdwtype       = Cut-off
vdw-modifier  = Force-switch
rvdw-switch   = 1.0
rvdw          = 1.2
DispCorr      = no

AMBER (ff14SB, ff99SB-ILDN, GAFF/GAFF2 ligands):

MDP

coulombtype   = PME
rcoulomb      = 1.0
vdwtype       = Cut-off
vdw-modifier  = Potential-shift
rvdw          = 1.0
DispCorr      = EnerPres

AMBER force fields were developed with an 8–10 Å cut-off and a long-range dispersion correction. Values of 0.9–1.0 nm are typical in GROMACS ports. Two hard rules apply with the Verlet scheme: rcoulomb must equal rvdw (or be larger, with PME), otherwise grompp stops; and the cut-off must be less than half the shortest box vector. fourierspacing = 0.12 and pme-order = 4 are the defaults and give good PME accuracy. Do not loosen them to speed things up. Let mdrun’s PME tuning shift work between real and reciprocal space instead.

Time step and constraints

MDP

dt           = 0.002
constraints  = h-bonds

A 2 fs step with bonds to hydrogen constrained (LINCS for the solute, SETTLE for water, which is automatic) is the standard. constraints = all-bonds is not needed with modern force fields, and CHARMM36 was parametrised with h-bonds.

4 fs with hydrogen mass repartitioning (HMR). Since GROMACS 2022 you do not need to edit topologies: add mass-repartition-factor = 3 and grompp moves mass from heavy atoms onto their bonded hydrogens. With dt = 0.004 this roughly doubles throughput. On our Trp-cage test system it went from about 85 to 131 ns/day on 8 CPU threads. HMR changes kinetics slightly (dynamics are a little slower) but not thermodynamics, so it is widely used for sampling. If you report rates or diffusion, say that HMR was used. Equilibrate at 2 fs and switch to 4 fs for production.

Without constraints, keep dt at or below 1 fs. grompp will warn about the “estimated oscillational period” of the fastest bonds if it is too large.

Temperature coupling

MDP

tcoupl  = V-rescale
tc-grps = Protein Non-Protein   ; or simply: System
tau-t   = 0.1 0.1
ref-t   = 300 300

V-rescale (Bussi–Donadio–Parrinello) samples the canonical ensemble correctly and is stable from the first step, so you can use it for every stage. grompp now warns that Berendsen “does not generate the correct kinetic energy distribution”. Nosé–Hoover is also correct, but it oscillates if the system starts far from the target temperature.

On groups: the old advice to couple solute and solvent separately came from the “hot solvent / cold solute” problem with weak coupling. With V-rescale, a single System group is fine and avoids index-file bookkeeping. If you do split into groups, each must be large enough to have meaningful kinetic energy. Never couple a single small ligand on its own; merge it with the protein (Protein_LIG) using make_ndx. Whatever you choose, tau-t and ref-t must have exactly one value per group.

Pressure coupling

MDP

pcoupl          = C-rescale
pcoupltype      = isotropic
tau-p           = 1.0     ; equilibration; 2.0–5.0 is also fine
ref-p           = 1.0
compressibility = 4.5e-5
refcoord-scaling = com    ; only with position restraints

C-rescale (stochastic cell rescaling, Bernetti & Bussi 2020) is the modern default. It relaxes the box as quickly as Berendsen but gives correct volume fluctuations, so it works for both equilibration and production. Parrinello-Rahman is also correct for production, but it can oscillate badly if the box is far from equilibrium. Use it only after a C-rescale equilibration.

Constraints that grompp enforces:

Steps per coupling period. tau-p must be at least 25 × nstpcouple × dt for C-rescale and Berendsen, and 200× for Parrinello-Rahman. With nstpcouple = -1 grompp picks a valid value.
Number of values. ref-p and compressibility need one value for isotropic, two for semiisotropic (xy, then z; used for membranes), and six for anisotropic.
Compressibility. 4.5 × 10⁻⁵ bar⁻¹ is water at about 300 K. For membranes with semiisotropic coupling, use 4.5e-5 4.5e-5.

When position restraints and pressure coupling are on together, set refcoord-scaling = com. Otherwise grompp warns that the absolute restraints “will give artifacts”.

The equilibration stages

NVT (about 100 ps)

MDP

define      = -DPOSRES
gen-vel     = yes
gen-temp    = 300
gen-seed    = -1
continuation = no
pcoupl      = no

Velocities are drawn from a Maxwell distribution and the box volume is fixed, while heavy atoms of the solute are restrained to their minimised positions (-r em.gro on the grompp command line). Check that the temperature reaches the target and stays there: gmx energy -f nvt.edr, select Temperature.

NPT (about 100 ps or more)

MDP

define       = -DPOSRES
gen-vel      = no
continuation = yes
pcoupl       = C-rescale
refcoord-scaling = com

Continue from the NVT checkpoint (-t nvt.cpt) and let the density settle. For a protein in TIP3P water at 300 K and 1 bar, expect around 1,000–1,020 kg m⁻³ depending on solute fraction. Pressure fluctuates by hundreds of bar in a system of this size. That is normal. Judge the running average, not individual frames.

Membrane systems and systems with large ligands, cofactors or metal sites often need a longer equilibration with gradually weaker restraints. CHARMM-GUI’s six-step protocol is a good model.

Production

MDP

continuation = yes
gen-vel      = no
; no define = -DPOSRES
pcoupl       = C-rescale   ; or Parrinello-Rahman
tau-p        = 2.0

Remove the restraints, continue from npt.cpt, and set nsteps to the time you need: nsteps = time_in_ps / dt, so 100 ns at 2 fs is 50,000,000 steps. Write the trajectory with nstxout-compressed every 10–100 ps. Analysis rarely needs more, and a 100-ns run of a 40,000-atom system at 10 ps intervals already produces roughly 1–2 GB of .xtc. Set nstxout, nstvout and nstfout to 0. Checkpoints hold the full-precision state needed to restart.

Other parameters you may meet

comm-mode = Linear, nstcomm = 100. This removes the drift of the system’s centre of mass. Use Angular only for an isolated molecule in vacuum.
nstenergy and nstlog. Every 5,000–10,000 steps is plenty for monitoring.
compressed-x-grps. Write only the protein and ligand to the .xtc file to save disk space. Remember that later analysis then needs a matching .tpr or index.
free-energy, pull, awh. These are separate topics with their own rules. Do not copy them from another project without understanding them.

A pre-submission check

Before queuing a long run, paste the file into the mdp checker. It applies all of the rules on this page offline (cut-off consistency, coupling-constant ratios, group counts, continuation logic) and also shows total simulated time and approximate output size. Then run grompp with -maxwarn 0 and do a short test with mdrun. If either complains, the error guide quotes the exact messages.

References

GROMACS 2026 documentation: Molecular dynamics parameters (.mdp options) and Force fields in GROMACS.
Bussi G, Donadio D, Parrinello M. Canonical sampling through velocity rescaling. J Chem Phys 2007;126:014101. doi:10.1063/1.2408420
Bernetti M, Bussi G. Pressure control using stochastic cell rescaling. J Chem Phys 2020;153:114107. doi:10.1063/5.0020514
Hopkins CW, Le Grand S, Walker RC, Roitberg AE. Long-time-step molecular dynamics through hydrogen mass repartitioning. J Chem Theory Comput 2015;11:1864–1874. doi:10.1021/ct5010406
Huang J, MacKerell AD Jr. CHARMM36 all-atom additive protein force field: validation based on comparison to NMR data. J Comput Chem 2013;34:2135–2145. doi:10.1002/jcc.23354
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
Hess B. P-LINCS: A parallel linear constraint solver for molecular simulation. J Chem Theory Comput 2008;4:116–122. doi:10.1021/ct700200b