GROMACS Position Restraints: gmx genrestr and POSRES
Position restraints in GROMACS hold selected atoms near reference coordinates with a harmonic force while the solvent equilibrates. pdb2gmx writes posre.itp for the protein only, so for a ligand you build an index group of its heavy atoms with gmx make_ndx, run gmx genrestr on the ligand’s own coordinate file, include the output inside the ligand’s [ moleculetype ] under #ifdef POSRES, and switch it on with define = -DPOSRES.
Most first protein-ligand simulations go wrong at the same point. The protein stays put during NVT and NPT because pdb2gmx gave it a restraint file, but the ligand has none, so it drifts out of the pocket before production even starts. Or the student does generate a ligand restraint file, pastes the include line in the wrong place in topol.top, and grompp stops with an “out of bounds” error that makes no sense on first reading. This guide from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modelling) walks through the whole mechanism: what posre.itp contains, how gmx genrestr numbers atoms, where the #ifdef block must sit, how to pick a force constant, and how to read the errors grompp prints when any of this is wrong.
It assumes you already have a protein-ligand complex built and solvated, as in our protein-ligand MD simulation tutorial, and a ligand topology from a parameterisation server (see how to parameterise a ligand for GROMACS). For the full sequence of GROMACS steps, start from our molecular dynamics with GROMACS learning path.
What is a position restraint in GROMACS?
A position restraint is a harmonic potential that pulls an atom back toward a fixed reference position. The GROMACS reference manual section on restraints gives the form as V = ½ k |r − R|², where r is the atom’s current position, R is its reference position and k is the force constant in kJ mol−1 nm−2. The manual also lets you give three separate force constants for x, y and z, so an atom can be held to a point, a line or a plane.
The atoms are not frozen. Justin Lemkul’s GROMACS lysozyme tutorial puts it plainly: “Movement is permitted, but only after overcoming a substantial energy penalty.” That is the whole point during equilibration. Water and ions reorganise around the solute and the box reaches the right density, while the protein and ligand stay close to the geometry you prepared and checked.
The reference positions R do not come from the restraint file. They come from the coordinate file you pass to gmx grompp with -r. The restraint file only says which atoms are restrained and how strongly.
What does posre.itp contain, and why is the ligand missing from it?
When you run gmx pdb2gmx, it writes three files: the processed structure, topol.top, and posre.itp. The gmx pdb2gmx documentation lists -i (default posre.itp) as the include file it writes and -posrefc (default 1000) as its force constant. The file restrains the protein’s heavy atoms, which means every atom that is not a hydrogen. Each line looks like this:
[ position_restraints ]
; atom type fx fy fz
1 1 1000 1000 1000
5 1 1000 1000 1000The columns are the atom index, the function type (1 is the ordinary harmonic restraint) and the three force constants. pdb2gmx also writes this block into topol.top, directly after the protein’s bonded terms:
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endifpdb2gmx only processes residues its force field knows, so a ligand, a cofactor or a drug molecule gets no restraint file. That is the gap gmx genrestr fills. If your protein has several chains, pdb2gmx writes one topology and one restraint file per chain (for example topol_Protein_chain_A.itp and posre_Protein_chain_A.itp), each already wired into the correct molecule.
How do you generate a ligand restraint file with gmx genrestr?
The workflow below follows Lemkul’s protein-ligand complex tutorial, which uses T4 lysozyme L99A (PDB 3HTB) with the ligand JZ4. The JZ4 coordinate file has 22 atoms, 10 of them heavy atoms, and the assembled complex has 2636 atoms before solvation. Replace jz4 with your own ligand’s name.
Step 1: Make an index group of the ligand’s heavy atoms
Run gmx make_ndx on the ligand’s own coordinate file, not on the complex:
gmx make_ndx -f jz4.gro -o index_jz4.ndx
> 0 & ! a H*
> qThe selection 0 & ! a H* means “group 0 (System) and not any atom whose name starts with H”. It creates a new group containing only heavy atoms. In the tutorial this new group is group 3. Our gmx make_ndx index groups tutorial explains the selection syntax in more detail.
Step 2: Run gmx genrestr on that group
gmx genrestr -f jz4.gro -n index_jz4.ndx -o posre_jz4.itp -fc 1000 1000 1000When prompted, pick the heavy-atom group you just made. The flags, as listed in the gmx genrestr manual page:
-f: the structure file (defaultconf.gro). Use the ligand alone.-n: the index file (optional). Without it you choose from the default groups.-o: the output include file (defaultposre.itp). Always rename it, or you will overwrite the protein’s file.-fc: force constants in kJ mol−1 nm−2 (default 1000 1000 1000). A single value applies the same constant in all three directions.
Step 3: Check the atom numbers
Open posre_jz4.itp. The indices must run from 1 up to the number of atoms in the ligand (22 for JZ4), skipping hydrogens. The genrestr manual carries an explicit warning that atom numbers in every [ moleculetype ] start at 1, so the tool “will only produce a useful file for the first molecule” in its input. Run it on the whole complex and you get system numbering, where the ligand’s atoms come after every protein atom, so the indices belong to no single molecule and grompp rejects them.
Want the guided, hands-on version?
Our live Molecular Modeling & MD Simulations cohort bootcamp takes you from zero to running real docking and MD workflows, with a portfolio project for your grad-school applications.
Where does the #ifdef POSRES block go in topol.top?
This is the step that produces most of the errors. A [ position_restraints ] block belongs to whichever [ moleculetype ] grompp is reading when it meets the block. So the ligand’s restraint include has to sit after the ligand’s [ atoms ] and bonded sections, and before the next [ moleculetype ] begins. The GROMACS common errors page shows the pattern for several molecules. Each restraint include follows its own molecule:
#include "topol_A.itp"
#ifdef POSRES
#include "posre_A.itp"
#endif
#include "topol_B.itp"
#ifdef POSRES
#include "posre_B.itp"
#endif
#include "ligand.itp"
#ifdef POSRES
#include "ligand_posre.itp"
#endifGrouping all three restraint includes together after the last molecule is the layout the same page labels WRONG, because every restraint would then be read as part of the ligand.
There is a second ordering rule. A ligand topology from CGenFF, ACPYPE or a similar tool usually brings new atom types or bonded parameters, and those must appear before any [ moleculetype ]. In Lemkul’s tutorial the ligand file goes right after the force field include, at the top of topol.top:
; Include forcefield parameters
#include "./charmm36-feb2026_cgenff-5.0.ff/forcefield.itp"
; Include ligand parameters
#include "jz4_gmx.itp"Because jz4_gmx.itp holds the JZ4 [ moleculetype ] itself, the restraint block goes inside that file, at the end of the ligand’s sections:
; Include Position restraint file
#ifdef POSRES
#include "posre_jz4.itp"
#endifFinally, the [ molecules ] directive at the bottom of topol.top must list the ligand after the protein, in the same order as the coordinate file. Our GROMACS topology file explained post walks through every directive in order.
How do you switch the restraints on and off in the .mdp file?
The #ifdef POSRES lines are preprocessor conditions. The block is included only when POSRES is defined, and you define it in the .mdp file. The GROMACS mdp options reference describes define as “defines to pass to the preprocessor” and notes that -DPOSRES triggers inclusion of posre.itp. In Lemkul’s complex tutorial both nvt.mdp and npt.mdp carry:
define = -DPOSRES ; position restrain the protein and ligandProduction md.mdp has no define line, so the restraints disappear. Two other settings travel with the restraints:
- The
-rflag on grompp. Since GROMACS 2018 you must supply the reference coordinates explicitly:gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr, then-c nvt.gro -t nvt.cpt -r nvt.grofor NPT. refcoord_scalingunder pressure coupling. The tutorial’s npt.mdp setsrefcoord_scaling = com, which scales the centre of mass of the reference coordinates with the box. The mdp reference warns that with the defaultnothe virial and pressure “might be ill defined”.
To check that a restraint really was included, run grompp with -pp processed.top. The gmx grompp manual describes this as writing out the pre-processed topology, so you can search it for your ligand’s [ position_restraints ] block. Our NVT and NPT equilibration guide covers the rest of these two runs.
How do you restrain the protein and ligand separately?
Use a different macro name for the ligand. The complex tutorial shows this alternative:
; Ligand position restraints
#ifdef POSRES_LIG
#include "posre_jz4.itp"
#endifNow define = -DPOSRES restrains only the protein, define = -DPOSRES_LIG only the ligand, and define = -DPOSRES -DPOSRES_LIG both. The tutorial adds that restraining both together is probably enough for a standard equilibration.
The same mechanism lets you change the force constant without regenerating files. The GROMACS preprocessor accepts a value with a define, so you can write a macro name in place of the number in the restraint file:
[ position_restraints ]
1 1 POSRES_FC POSRES_FC POSRES_FCand set it per run with define = -DPOSRES -DPOSRES_FC=500. CHARMM-GUI’s GROMACS inputs use this pattern, with separate macros for backbone and side chains (POSRES_FC_BB and POSRES_FC_SC), to step restraints down over several equilibration stages. If you take this route, every run must define the macro, or grompp will find a word where it expects a number.
What force constant should you use?
1000 kJ mol−1 nm−2 is the default in both pdb2gmx and genrestr, and it is the value the Lemkul tutorials use. The numbers show why it works. With k = 1000, an atom displaced 0.1 nm (1 Å) from its reference costs ½ × 1000 × 0.1² = 5 kJ/mol. Thermal energy at 300 K is RT = 0.008314 × 300 ≈ 2.49 kJ/mol. By equipartition, the mean-square displacement along each axis is kT/k, so the typical thermal wobble is about √(2.49/1000) ≈ 0.05 nm, or 0.5 Å per axis. The atoms move, but stay near where you put them.
| Force constant (kJ mol−1 nm−2) | Typical RMS wobble per axis at 300 K | When students use it |
|---|---|---|
| 1000 | about 0.05 nm (0.5 Å) | Default NVT and NPT for protein and ligand heavy atoms |
| 100 | about 0.16 nm (1.6 Å) | A softer later stage before production |
| 10 | about 0.5 nm (5 Å) | A last, very weak stage; barely holds a small ligand |
The wobble values are calculated from √(kT/k) at 300 K. They show scale, not a rule. A stronger constant than 1000 rarely helps and can make the restrained atoms fight the thermostat. A ligand that leaves its pocket the moment restraints are released usually points to a problem with the pose or the ligand topology, which no force constant will fix.
Which restraint method should you choose?
gmx genrestr and the GROMACS topology format support several ways to hold atoms. They behave differently, and the reference manual stresses that restraints and constraints are separate algorithms.
| Method | How you set it up | What it does | Typical use |
|---|---|---|---|
| Harmonic position restraint (function type 1) | posre.itp from pdb2gmx or gmx genrestr, plus define = -DPOSRES | Pulls atoms toward reference positions in x, y and z | Standard NVT and NPT equilibration |
| Flat-bottomed position restraint (function type 2) | Written by hand in a [ position_restraints ] block | No force inside a radius r, harmonic outside it | Keeping a ligand inside a region without fixing its position |
| Distance restraints | gmx genrestr -disre | Holds pairs of atoms near their current distance, not a fixed point in space | Keeping a protein’s shape while letting it move |
| Freeze groups | freezegrps and freezedim in the .mdp | Atoms do not move at all in the frozen dimensions | Special cases such as fixed surfaces; not for normal equilibration |
How do you fix the common position restraint errors?
These are the messages grompp prints, quoted from the GROMACS source and error documentation.
“Atom index (n) in position_restraints out of bounds (1-N)”
Here n is the offending index and N is the number of atoms in the molecule grompp was reading (22 for JZ4). The full message continues: “This probably means that you have inserted topology section “position_restraints” in a part belonging to a different molecule than you intended to.” Two causes account for almost every case:
- The include sits under the wrong molecule. The protein’s restraint file ended up after the ligand’s include, so grompp tries to apply protein atom numbers to a 22-atom ligand. Move each #ifdef block directly below its own molecule.
- genrestr was run on the complex. The indices are system numbers in the thousands instead of 1 to 22. Rerun genrestr on the ligand’s own .gro file.
“Invalid order for directive atomtypes”
The ligand file adds atom types after a [ moleculetype ] has already started, usually because it was included after the protein. The errors page says new types and parameters “must appear before any [moleculetype] directive”. Move the ligand include up, directly below the force field include.
“Cannot find position restraint file restraint.gro (option -r)”
POSRES is defined but you did not pass -r. The message itself explains that since GROMACS 2018 the restraint coordinates must be given explicitly, and you can reuse the file you give to -c. Add -r em.gro for NVT or -r nvt.gro for NPT.
The ligand still drifts during NVT, with no error
The restraint was never included. Common reasons: the #ifdef uses POSRES_LIG but the .mdp only defines POSRES, the define line is missing from nvt.mdp, or CGenFF’s file still says #include "posre.itp" instead of your ligand file (the tutorial lists renaming this line as one of the required edits). Run grompp with -pp and search the processed topology to confirm.
Do you need position restraints in the production run?
No, for a standard protein-ligand study. Production should be unrestrained, because restraints bias exactly the motions you want to measure: ligand stability, pocket flexibility and RMSD. Remove the define line or leave it empty, keep the same topology, and drop -r from grompp. If a reviewer could ask “was the ligand held in place?”, your methods section should state the restraint stages and force constants used in equilibration; our MD methods section guide shows how to word this.
Frequently asked questions
Does gmx genrestr restrain hydrogen atoms?
Only if they are in the group you select. The usual practice, matching pdb2gmx, is to restrain heavy atoms only by building an index group with 0 & ! a H*. Hydrogens then follow their heavy atoms through bonds or constraints.
Can I use gmx genrestr on a .pdb file instead of .gro?
Yes. The manual lists gro, g96, pdb, brk, ent, esp and tpr as accepted -f formats. What matters is that the file contains only the molecule you are restraining, with atoms in the same order as its topology.
What is the difference between -DPOSRES and -DPOSRES_LIG?
They are just macro names. -DPOSRES is the name pdb2gmx writes for the protein; POSRES_LIG is a second name you choose for the ligand so you can switch the two restraint sets on and off independently in the .mdp file.
Should I restrain water during equilibration?
Normally no; water is what you are equilibrating. Force field water files include an optional #ifdef POSRES_WATER block for special cases, and it stays off unless you define it.
Why set refcoord_scaling when restraints are on during NPT?
With pressure coupling the box changes size, but restraint reference positions do not unless refcoord_scaling tells GROMACS to scale them. Setting refcoord_scaling = com in npt.mdp, as the Lemkul tutorial does, is the common choice.
Once the ligand stays in its pocket through NVT and NPT, the next stages are production MD and analysis. See where restraints, equilibration and trajectory analysis fit among the other skills in our computational biology skills roadmap, and keep our common GROMACS errors guide open for the next grompp message.
Want the guided, hands-on version?
Our live Molecular Modeling & MD Simulations cohort bootcamp takes you from zero to running real docking and MD workflows, with a portfolio project for your grad-school applications.