The GROMACS .mdp File Explained: The Parameters That Actually Matter

The .mdp file is the run-parameter file GROMACS grompp reads to build a .tpr. Of its 100+ options, roughly ten decide whether your simulation is correct: integrator, dt, nsteps, constraints, cutoff-scheme, rcoulomb and rvdw, coulombtype, tcoupl, pcoupl, gen-vel and pbc. The rest can stay at their documented defaults.

Most students meet the .mdp file the same way. You download em.mdp, nvt.mdp, npt.mdp and md.mdp from a tutorial, they work, and you never look inside. Then an examiner asks why you used a 2 fs time step, or why your barostat is what it is, and there is no answer. This guide from the StemSkills Lab team (10+ years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling) walks the parameters that carry real physics, using the actual defaults from the official documentation, so you can write and defend your own files. It sits inside our pillar guide on learning molecular dynamics with GROMACS.

What is a .mdp file and what does grompp do with it?

The .mdp (molecular dynamics parameters) file is plain text. Each line sets one option. It is not read by gmx mdrun. It is read by gmx grompp, the pre-processor, which combines three inputs into one binary run file:

• -f your .mdp file (how to run)
• -c a structure file (where the atoms are)
• -p a topology (what the atoms are, and how they interact)

The output is a .tpr, and this is the detail that matters: the .tpr is frozen. Editing the .mdp after grompp changes nothing. You must re-run grompp. That single fact explains a large share of confused debugging, where a student edits md.mdp, re-runs mdrun, and sees the old behaviour.

One convention to know before reading any parameter list: in the GROMACS documentation, as the mdp options reference states, “Default values are given in parentheses, or listed first among choices.” Anything you omit from your file silently takes that default. Omission is a choice, and that is why the defaults below are worth knowing.

Which .mdp parameters actually matter?

These are the settings where the default is either wrong for your system or where the value is genuinely yours to justify. Every default in this table comes from the official mdp options reference.

Read the failure column again. Five of these defaults (nsteps, constraints, coulombtype, tcoupl, pcoupl) are wrong for an ordinary protein simulation. GROMACS is not being unhelpful. It refuses to guess your science for you.

Which integrator and time step should you use?

Use integrator = steep (steepest descent) for energy minimisation and integrator = md for dynamics. The md integrator is leap-frog, and the reference manual describes it plainly: “The default MD integrator in GROMACS is the so-called leap-frog algorithm… It produces trajectories that are identical to the Verlet algorithm.”

The time step is the parameter students most often cannot defend, so here is the reasoning. Your step must resolve the fastest motion in the system. That motion is bond vibration involving hydrogen, which has a period of roughly 10 fs, and a safe integrator needs several steps per period. That is why the default dt is 0.001 ps (1 fs).

To go faster you must remove the fast motion rather than ignore it. Setting constraints = h-bonds holds bonds to hydrogen at fixed length, so they no longer vibrate, and 2 fs becomes safe. The constraint is solved by LINCS (constraint-algorithm = LINCS by default), from Hess and co-workers’ paper LINCS: A linear constraint solver for molecular simulations.

So dt = 0.002 and constraints = h-bonds are a single decision, not two. Change one without the other and you have either a slow simulation or an unstable one. That is the answer to the examiner’s question.

How do you choose the cut-off and electrostatics settings?

Non-bonded interactions are split in two, and they are not symmetric.

Van der Waals decays fast, so truncating it at rvdw = 1 nm is physically reasonable. Electrostatics decays as 1/r, which is far too slow to cut off. Truncating it is the classic way to produce a plausible-looking but wrong trajectory, especially with charged residues or ions.

This is why the coulombtype = Cut-off default must be changed to coulombtype = PME. Particle Mesh Ewald computes long-range electrostatics under periodic boundary conditions in N·log(N) time, from Darden, York and Pedersen’s 1993 paper in the Journal of Chemical Physics (volume 98, pages 10089 to 10092). Its supporting parameters, fourierspacing (default 0.12 nm) and pme-order (default 4), are tuned defaults you should leave alone.

One rule that is not optional: match rcoulomb and rvdw to the force field you chose. Force fields are parameterised with a specific cut-off, and changing it means you are no longer running the force field you cited in your methods. If you are still choosing one, see our guide on the GROMACS workflow and force field selection.

Modern GROMACS uses cutoff-scheme = Verlet by default, which manages the neighbour list with a buffer set by verlet-buffer-tolerance (default 0.005 kJ mol⁻¹ ps⁻¹ per atom). Because that buffer is automatic, rlist (default 1 nm) is adjusted for you. Do not tune it by hand.

Which thermostat and barostat should you use?

Use tcoupl = v-rescale for temperature. Use pcoupl = C-rescale or Parrinello-Rahman for pressure in your NPT and production runs. Avoid the Berendsen versions of both.

The reason is a real physics distinction, and it is the most quotable thing in the manual. On the Berendsen thermostat, the GROMACS reference manual states it “does not generate a proper canonical ensemble, so rigorously, the sampling will be incorrect… we strongly recommend against using it for new simulations.” On velocity-rescale, the same manual says: “This thermostat produces a correct canonical ensemble and still has the advantage of the Berendsen thermostat: first order decay of temperature deviations and no oscillations.”

Berendsen suppresses the fluctuations that the canonical ensemble is supposed to have. Velocity-rescale adds a stochastic term that restores the correct kinetic energy distribution, from Bussi, Donadio and Parrinello’s Canonical sampling through velocity rescaling (Journal of Chemical Physics, volume 126, article 014101, 2007). Many tutorials still circulate with tcoupl = berendsen. They are out of date.

Two supporting values you will be asked about. tau-t is the temperature coupling time constant, commonly 0.1 ps, and it must be given once per group in tc-grps. Couple the solute and solvent separately (a typical pair is Protein Non-Protein) so that a hot solute cannot hide inside a cold bath. For pressure, tau-p defaults to 5 ps and compressibility is set to the compressibility of your solvent, which for water at 300 K and 1 bar is 4.5e-5 bar⁻¹.

When should gen-vel be yes?

In nvt.mdp only. This is the parameter that quietly ruins equilibration protocols.

gen-vel = yes draws fresh velocities from a Maxwell distribution at gen-temp (default 300 K), with gen-seed (default -1, meaning a pseudo-random seed). You want that exactly once, when you first give your minimised, motionless structure some heat.

After that, your NPT and production runs must continue from the velocities you already equilibrated. So they set gen-vel = no and you carry the state forward by passing the checkpoint file to grompp:

• gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr

The -t nvt.cpt flag is doing the real work. Set gen-vel = yes in npt.mdp and you throw away the thermal equilibration you just spent hours on, then report it as an equilibrated system.

Which parameters should you leave at their defaults?

Not touching a parameter is a defensible answer when the default encodes tuned expertise. Leave these alone unless you have a specific reason:

• cutoff-scheme = Verlet, verlet-buffer-tolerance = 0.005, and rlist = 1: the buffer is set automatically
• constraint-algorithm = LINCS, lincs-order = 4, lincs-iter = 1: fine for a 2 fs step
• fourierspacing = 0.12 and pme-order = 4: standard PME accuracy
• nstlist = 10: mdrun tunes this at run time for performance anyway
• nsttcouple = -1: lets GROMACS pick an interval for accurate integration
• pbc = xyz: correct for any solvated box

The output-frequency options (nstxout-compressed, nstenergy, nstlog) are the harmless ones. They change file size and time resolution, not physics. Write too often on a long run and you will fill your disk.

What do the four files look like in practice?

A standard protein workflow uses four .mdp files, and they differ in only a few lines:

Note emtol, the minimisation convergence criterion. Its documented default is 10.0 kJ mol⁻¹ nm⁻¹, which is strict. Most tutorials relax it to 1000.0 because minimisation only needs to remove bad contacts before dynamics, not find a true minimum. That is a deliberate override of a default, and it is the kind of choice worth stating in your methods section.

For the order these files run in, and where MD sits relative to structure preparation and docking, see the computational biology skills roadmap.

Troubleshooting: .mdp errors and what they mean

These messages come from grompp or mdrun and point back to the .mdp file.

• “Unknown left-hand ‘XXX’ in parameter file”: a typo, or an option removed in your GROMACS version. Check spelling against the mdp options page for your version, not the latest.
• “The cut-off length is longer than half the shortest box vector”: not an .mdp bug. Your box is too small for your rcoulomb. Rebuild it with more padding in gmx editconf, or an atom would see its own periodic image.
• “Number of grps in tc-grps does not match the number of tau-t / ref-t values”: tc-grps, tau-t and ref-t must have the same count. Two groups need two values each.
• “You are generating velocities and reading a checkpoint”: gen-vel = yes while passing -t. Pick one. In npt.mdp the answer is gen-vel = no.
• The run finishes instantly with no trajectory: nsteps is at its default of 0. You never set the length.
• Repeated LINCS warnings: usually dt = 0.002 without constraints = h-bonds, or a system that was never properly minimised. Fix the pair, then check em.mdp converged.

For the wider set of messages beyond the parameter file, see our guide to common GROMACS errors.

Frequently asked questions

Can I use the same .mdp file for every system?
Mostly yes, with two exceptions. rcoulomb and rvdw must follow your force field, and tc-grps must name groups that exist in your system. A membrane system also needs pcoupltype = semiisotropic rather than the default isotropic.

Why does GROMACS ignore my edits to the .mdp file?
Because mdrun reads the .tpr, not the .mdp. Any change requires re-running gmx grompp to rebuild the .tpr. This is the single most common source of confusion for beginners.

Is a 4 fs time step possible?
Only if you remove more fast motion, typically with virtual sites for hydrogens. Simply setting dt = 0.004 with constraints = h-bonds will make your system unstable. Angle vibrations involving hydrogen are still active at that step.

Should I use C-rescale or Parrinello-Rahman?
Both give correct volume fluctuations. C-rescale is stochastic and tolerant of a system that is still settling, which makes it well behaved during NPT equilibration when the density is still moving. Parrinello-Rahman is a common production choice once the density is already stable. Either is defensible if you state which you used.

What does the .mdp file have to do with my methods section?
It is most of it. Integrator, time step, constraints, cut-offs, electrostatics, thermostat, barostat and run length are exactly what a reviewer expects to see reported. If you can annotate your .mdp, your methods section is already written.

Do I need position restraints in equilibration?
For a standard protein protocol, yes. define = -DPOSRES in nvt.mdp and npt.mdp holds heavy atoms near their starting positions while solvent relaxes around them, then you release them for production.

Where this fits in your workflow

The .mdp file is where your simulation stops being a tutorial and starts being your science. Every value in it is a claim about physics that you will be asked to justify, in a viva, in a methods section, or by a reviewer. The good news is the list is short. Ten parameters carry the argument, the documented defaults handle the rest, and the reasoning behind each one (fast motion sets the time step, 1/r decay demands PME, correct ensembles rule out Berendsen) is the same reasoning you will reuse in every simulation you ever run.