How to Perform Energy Minimization in GROMACS: A Step-by-Step Tutorial

Energy minimization in GROMACS relaxes a starting structure so no atoms sit in high-energy overlaps before dynamics begin. You write an em.mdp file with the steepest-descent integrator, build a run input with gmx grompp, then run gmx mdrun. Minimization is done when the maximum force (Fmax) drops below your emtol threshold.

You have solvated your protein, added ions, and neutralised the box. The next command you type decides whether your simulation runs or crashes in its first few steps. Skip energy minimization and you carry hidden strain into dynamics: two atoms placed almost on top of each other, a clashing side chain, a water molecule wedged into the protein core. Those tiny geometric problems create enormous forces, and the integrator responds by flinging atoms across the box. This guide from the StemSkills Lab team, who have spent 10+ years running structural bioinformatics and production molecular dynamics for drug-discovery and materials projects, walks through the energy-minimization step with the exact commands and settings, shows you how to confirm it worked, and fixes the errors that trip students up. New to the engine itself? Start with our GROMACS molecular dynamics guide, then come back here for the minimization step.

What is energy minimization in GROMACS, and why run it before MD?

Energy minimization is a geometry-relaxation step that moves atoms downhill on the potential-energy surface until the forces on them are small. It is not molecular dynamics: there are no velocities, no temperature, and no time. The algorithm simply nudges each atom in the direction that lowers the energy, one step at a time, until the structure reaches a nearby local minimum.

The reason it matters is numerical stability. A structure straight from a PDB file or a solvation tool almost always contains close contacts. When two atoms sit too close, the Lennard-Jones repulsion between them is huge, so the force is huge. Molecular dynamics integrates those forces over a finite timestep (commonly 2 fs), and a huge force over even a small timestep displaces an atom by an unphysical distance. That is what people mean when a run “blows up” or throws LINCS warnings. Minimization removes the worst clashes first, so dynamics starts from a sane geometry. The GROMACS reference manual documents steepest descent as the standard minimizer, which is exactly why it is the default first pass: it stays stable even when the starting structure is badly strained.

What files do you need before energy minimization?

Minimization is a mid-workflow step, so it assumes you already have three inputs from the setup stage:

• A coordinate file (for example solv_ions.gro): your protein solvated in a water box with counter-ions added.
• A topology file (topol.top): the force-field definition of every atom, bond, angle, and interaction in the system.
• A parameter file (em.mdp): the run settings you are about to write, telling GROMACS to minimize rather than to integrate dynamics.

If you followed the canonical Lysozyme in Water tutorial by Justin A. Lemkul, these are the exact filenames produced by the gmx solvate and gmx genion steps. If your files are named differently, substitute your own names in the commands below.

How do you write an em.mdp file for energy minimization?

The .mdp file is a plain-text list of options. For minimization the key line is the integrator, which you set to a minimizer rather than a dynamics integrator like md. A standard beginner em.mdp looks like this:

; em.mdp - passed to grompp to generate em.tpr
integrator  = steep      ; steepest-descent energy minimization
emtol       = 1000.0     ; stop when max force < 1000 kJ/mol/nm
emstep      = 0.01       ; initial step size (nm)
nsteps      = 50000      ; maximum number of minimization steps

nstlist         = 1        ; update neighbor list every step
cutoff-scheme   = Verlet   ; buffered neighbor searching
coulombtype     = PME      ; particle-mesh Ewald for electrostatics
rcoulomb        = 1.0      ; short-range electrostatic cut-off (nm)
rvdw            = 1.0      ; short-range van der Waals cut-off (nm)
pbc             = xyz      ; periodic boundaries in all three dimensions

Three parameters control when minimization stops. emtol is the convergence target: GROMACS keeps going until the maximum force on any atom falls below this value, in kJ/mol/nm. The official mdp-options documentation lists the default emtol as 10.0 kJ/mol/nm, which is stricter than the 1000.0 used here. For pre-dynamics cleanup, 1000.0 is a common and perfectly adequate target. emstep is the initial step size in nm, and nsteps caps how many steps run before GROMACS gives up. Minimization stops at whichever comes first: Fmax below emtol, or nsteps reached.

How do you run grompp and mdrun for energy minimization?

GROMACS runs in two moves. First gmx grompp (the GROMACS preprocessor) combines your coordinates, topology, and parameters into a single binary run-input file, em.tpr. Then gmx mdrun executes that file.

gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr

The flags map directly to your files: -f is the parameter file, -c the coordinate file, -p the topology, and -o the output run-input file. Read the grompp output before continuing. If it prints notes or warnings, understand each one; grompp catching a problem here is far cheaper than mdrun failing later. Once em.tpr is written, run:

gmx mdrun -v -deffnm em

The -v flag is verbose mode, which prints the step number, potential energy, and Fmax as minimization proceeds so you can watch it converge. -deffnm em sets the default filename for every output, so you get em.gro (the minimized structure), em.edr (energies), and em.log (the full log) without naming each one. On a solvated protein this step usually finishes in seconds to a couple of minutes on a laptop.

How do you read the potential energy and Fmax convergence output?

When mdrun finishes, it prints a short summary. A converged run reports something in this form (numbers depend entirely on your system, so treat these as placeholders):

Steepest Descents converged to Fmax < 1000 in NNN steps
Potential Energy  = -X.XXXe+05
Maximum force     =  X.XXe+02 on atom ...

Two signals tell you the step succeeded. First, the line should say the run converged to Fmax below your emtol, not that it “did not converge”. Second, the potential energy should be large and negative for a solvated biomolecular system, because attractive interactions dominate once clashes are gone. A positive potential energy after minimization almost always means a serious setup problem upstream.

To inspect the energy trajectory rather than just the final value, use gmx energy on the .edr file:

gmx energy -f em.edr -o potential.xvg

At the interactive prompt, type Potential and press Enter, then type 0 and press Enter to finish. GROMACS writes potential.xvg, a plot of potential energy against step number. A healthy minimization curve drops steeply at first, then flattens as the structure settles into its local minimum. If the curve is still falling sharply at the last step, minimization stopped early because it hit nsteps, not because it converged.

Steepest descent vs conjugate gradient vs l-bfgs: which minimizer should you use?

GROMACS ships three minimization algorithms, set through the integrator keyword. Steepest descent is the default recommendation for the first pass because it is stable even when the starting structure is bad. Conjugate gradient and l-bfgs converge more efficiently near a minimum but are less forgiving of large initial forces.

For a standard protein-in-water workflow, steepest descent alone is enough before equilibration. A common two-stage habit is steepest descent followed by conjugate gradient when you need a tighter minimum, for example before a normal-mode analysis.

What are the common energy-minimization errors, and how do you fix them?

Most minimization failures come from the setup step, not from the em.mdp file. Here are the ones students hit most often.

• “number of coordinates in coordinate file does not match topology”: your .gro and .top disagree on atom count. This usually means you edited one file but not the other, or added ions to the structure without updating the [ molecules ] section of the topology. Check that every molecule count in topol.top matches the actual contents of solv_ions.gro.
• “System has non-zero total charge” as a grompp note: you have not fully neutralised the box. Go back to the gmx genion step and add the correct number of counter-ions before minimizing.
• “Steepest Descents did not converge to Fmax < … in … steps”: minimization ran out of steps. If the potential energy is already low and stable, the structure is often usable anyway. If not, raise nsteps, or check for a real geometry problem such as overlapping atoms from a bad insertion.
• mdrun stops with atoms moving unusually far, or the energy is positive: this points to a severe clash the minimizer cannot escape, frequently a misplaced ligand or a water molecule generated inside the protein. Inspect the starting structure in a viewer such as VMD or PyMOL before minimizing again.
• Wrong integrator: if you accidentally leave integrator = md in em.mdp, GROMACS runs dynamics, not minimization, and a clashing structure will explode. Confirm the integrator line reads steep, cg, or l-bfgs.

When minimization converges and the potential energy is low and negative, you are ready for the next stage. Minimization only settles the geometry; it does not give atoms velocities or bring the system to temperature and pressure. That is the job of NVT and NPT equilibration, which come next in the workflow. To see where this fits in the full skill path, follow our computational biology skills roadmap.

Frequently asked questions

Do I always need to run energy minimization before MD?

Yes, for essentially every biomolecular system. Any structure built from a PDB file and solvated will contain close contacts, and running dynamics on those produces enormous forces that destabilise the integrator. Minimization is the standard first step after solvation and ion addition.

What does emtol actually control?

emtol is the convergence threshold for the maximum force on any single atom, in kJ/mol/nm. GROMACS stops minimizing once Fmax falls below emtol. The documented default is 10.0, while many tutorials use 1000.0 for a quicker pre-dynamics cleanup. A larger emtol converges faster but leaves a slightly less relaxed structure.

How many steps should minimization take?

There is no fixed answer, because it depends on how strained the starting structure is. Set nsteps as a generous cap (50000 is common) and let convergence stop the run early. Watch the verbose output or the potential.xvg curve to confirm it flattened rather than hitting the step limit.

Why is my potential energy positive after minimization?

A positive potential energy for a solvated system almost always signals an upstream setup error, such as an unneutralised box, a mismatched topology, or overlapping atoms. Recheck the solvation and ion steps and inspect the structure in a molecular viewer before trying again.

Can I use conjugate gradient instead of steepest descent?

You can, but steepest descent is the safer first pass because it tolerates the large forces in a fresh structure. Conjugate gradient (cg) converges more efficiently near a minimum, so it is best used as a second, refining minimization after steepest descent.