gmx rdf Tutorial: Radial Distribution Function in GROMACS

To calculate a radial distribution function in GROMACS, run gmx rdf -f traj.xtc -s topol.tpr -ref '<reference>' -sel '<selection>' -o rdf.xvg. It histograms distances from the reference positions to the selected positions over the trajectory and normalises them to g(r). Add -cn rdf_cn.xvg to get the coordination number, read at the first minimum after the first peak.

This guide is written by the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling). It covers three RDFs a student actually needs: water around water as a sanity check, water around a bound ligand, and Na⁺/Cl^– around a protein. Every option below is taken from the official gmx rdf documentation (GROMACS 2026.4 at the time of writing). If you have not run a production simulation yet, start with our GROMACS molecular dynamics learning path.

What is a radial distribution function, in plain terms?

A radial distribution function, g(r), tells you how likely you are to find a particle at distance r from a reference particle, compared with a completely uniform (ideal gas) arrangement at the same density. A value of 1 means “no structure, just average density”. A value above 1 means particles pile up at that distance; below 1 means they avoid it.

For liquid water the curve has a sharp first peak (the first hydration shell), a dip (the first minimum), a weaker second peak, and then it flattens to 1. The experimental benchmark is useful because it gives you a number to check your own output against. Skinner and co-workers combined four X-ray diffraction data sets and reported the first O-O peak of ambient water at 2.80 Å (0.280 nm) with a height of 2.57 (Skinner et al., J. Chem. Phys. 138, 074506, 2013). Water models used in biomolecular force fields land close to that position; the exact height depends on the model.

Mathematically, the coordination number up to a cutoff R is n(R) = 4πρ ∫ g(r) r² dr from 0 to R, where ρ is the number density of the selected particles. Bulk water has about 33.4 molecules per nm³ (55.5 mol/L multiplied by Avogadro’s number), which is the ρ that turns the water O-O curve into “about four neighbours”. You do not have to integrate by hand, because -cn does it for you.

What do I need before running gmx rdf?

• A trajectory (.xtc or .trr) from a production run, not equilibration.
• The run input file (.tpr) passed to -s. Use the .tpr, not a .gro: the default -rmpbc yes makes molecules whole using bond information, and -excl reads exclusions, and both need the topology that only the .tpr carries.
• Atom and residue names for your system. In most GROMACS force-field ports water is residue SOL with oxygen OW, and gmx genion names ions NA and CL by default. Your ligand name is whatever you gave it during parameterisation (we use LIG below). Check with grep on your .gro file before you start.

gmx rdf takes selections directly on the command line, so you rarely need an index file. If you prefer named groups, our gmx make_ndx tutorial shows how to build them and pass them with -n index.ndx.

How do I run gmx rdf for water around water?

Start with the water O-O RDF. It takes a minute, and if it does not peak near 0.28 nm, something is wrong with your system or your command, not with the science.

gmx rdf -f md.xtc -s md.tpr \
        -ref 'resname SOL and name OW' \
        -sel 'resname SOL and name OW' \
        -b 10000 -bin 0.002 -rmax 1.0 \
        -o rdf_ow_ow.xvg -cn cn_ow_ow.xvg

What each option does:

• -ref: the positions you measure from.
• -sel: the positions you measure to. You can give several selections and get one curve per selection in the same file.
• -b 10000: skip the first 10,000 ps (10 ns). -b and -e are in ps by default; use only the equilibrated part of the run.
• -bin 0.002: bin width in nm. 0.002 nm is the default and is fine for thousands of waters.
• -rmax 1.0: largest distance in nm. The default (0) means half the box length with PBC, which is usually more than you need and slower.

Because -ref and -sel here are the same atoms, each oxygen can pair with itself at r = 0. If you get a spike in the very first bin, add -cut 0.05, which clears the RDF below 0.05 nm. The manual calls -cut “a rougher alternative to exclude intra-molecular peaks”, compared with -excl, which we cover in the troubleshooting section.

How do I calculate the RDF of water around a bound ligand?

This is the RDF most docking and MD students want: is the binding site dry or wet, and does the ligand keep a hydration shell when bound? Measure from the ligand to water oxygens:

gmx rdf -f md.xtc -s md.tpr \
        -ref 'resname LIG' \
        -sel 'resname SOL and name OW' \
        -b 10000 -bin 0.01 -rmax 1.2 \
        -o rdf_lig_ow.xvg -cn cn_lig_ow.xvg

Two things change compared with water-water. First, the bin width goes up to 0.01 nm. With one ligand as reference you have far fewer pairs per frame than with thousands of waters, so a 0.002 nm bin gives a noisy, spiky curve. Second, you must decide what “the ligand” means:

• 'resname LIG' uses every ligand atom as a reference and averages over them. Good for “how much water sits near the ligand as a whole”.
• 'com of resname LIG' uses one point, the ligand centre of mass. Good for small, compact ligands; misleading for long, flexible ones, because the centre of mass can sit in empty space.
• 'resname LIG and name O1' uses one specific atom (replace O1 with a real atom name from your topology). Best for asking whether a particular hydroxyl or carbonyl stays hydrated.

The -seltype option (default atom) changes the output positions of all selections at once, for example to mol_com. Writing com of inside one selection is safer because it only changes the selection you mean. -selrpos (also default atom) sets the reference positions the selection engine uses when it evaluates position-based keywords; leave it at the default unless you know you need it.

How do I measure ions around a protein with gmx rdf?

For counter-ions around a protein, a plain atom-based RDF from every protein atom is hard to read. The -surf option measures from each selected ion to the closest atom of each reference molecule or residue instead:

gmx rdf -f md.xtc -s md.tpr \
        -ref 'group Protein' -surf mol \
        -sel 'name NA' 'name CL' \
        -b 10000 -bin 0.01 -rmax 1.5 \
        -o rdf_prot_ions.xvg

This writes one curve for Na⁺ and one for Cl^–. Two rules from the manual matter here. With -surf, the -ref selection must select atoms (centres of mass are not supported). And with -surf no volume normalisation is applied, because the bins around an irregular protein surface do not have a computable volume. So read this output as a relative count profile (are ions concentrated within about half a nanometre of the surface, or not?), not as a true g(r) that goes to 1.

If you want a properly normalised ion RDF, pick a specific site instead, for example the carboxylate oxygens of one residue: -ref 'resid 45 and name OE1 OE2' with -sel 'name NA'. This pairs well with our salt bridge analysis guide when you suspect an ion is competing with a lysine or arginine.

Which -norm option should I use?

The -norm option decides what the y-axis means. The default is right for almost every student project.

How do I read the first peak and the first minimum?

Open the .xvg in xmgrace (xmgrace rdf_ow_ow.xvg) or plot it in Python with our xvg plotting tutorial. Then read three numbers:

1. First peak position. The most probable contact distance. For water O-O it should sit near 0.28 nm. For a ligand polar atom to water oxygen, a peak in the same 0.27 to 0.30 nm range suggests directly hydrogen-bonded water.
2. First peak height. How strongly structured that shell is. A ligand atom buried in a hydrophobic pocket may show no peak at all, which is a real result: the site is dry.
3. First minimum position. The outer edge of the first shell. This is the distance you use for the coordination number.

Write the numbers with the method, not just the plot. “The first minimum of the ligand O3 to water oxygen RDF was at 0.34 nm” is a sentence a reviewer can check; “water was present” is not. Our guide to writing the results section of a docking and MD study shows how to phrase this.

How do I get a coordination number from -cn?

The manual says “Option -cn produces the cumulative number RDF, i.e. the average number of particles within a distance r.” In practice:

1. Find the first minimum in rdf.xvg, say at r = 0.34 nm.
2. Open rdf_cn.xvg and read the y value at that same r.
3. That value is the average number of selected particles within 0.34 nm of each reference position: the coordination number.

A quick way to read it without a plot (it skips the xvg header lines and prints the first row at or beyond 0.34 nm):

grep -v '^[#@]' cn_ow_ow.xvg | awk '$1 >= 0.34 {print; exit}'

For water O-O this gives a value close to four in common rigid water models, matching the tetrahedral hydrogen-bond picture. For Na⁺ to water oxygen the first-shell count is commonly reported as five to six, depending on the ion parameters. If your number is far outside these ranges, check the troubleshooting list below before you interpret it.

Why does my gmx rdf output look wrong? (Troubleshooting)

If the error is not in this table, our list of common GROMACS errors and fixes covers file and topology problems that affect every analysis tool.

How long a trajectory do I need for a good RDF?

There is no single number, but the rule is simple: you need enough independent pairs per bin. For water around water, a few nanoseconds already gives a smooth curve because thousands of waters contribute every frame. For one ligand atom to water, or one site to a sparse ion, you are averaging a handful of pairs per frame, so use the full equilibrated production run and compare halves of the trajectory (-b and -e). If the first-peak height changes noticeably between the first and second half, the RDF is not converged. Running independent replicas, as described in our guide on how many replicate MD simulations to run, is the more defensible way to report error bars.

Where does gmx rdf fit in a full MD analysis?

An RDF answers a solvation or packing question. It sits alongside stability checks such as RMSD and RMSF and interaction checks such as protein-ligand hydrogen bonds. A typical protein-ligand chapter reports stability first, then contacts, then (if the binding site involves water) the ligand-water RDF and coordination number. For the method itself, cite the GROMACS paper (Abraham et al., SoftwareX 1-2, 19-25, 2015) and state the version you ran. To see where trajectory analysis sits among the other skills, use our computational biology skills roadmap.

Frequently asked questions

What is the difference between -ref and -sel in gmx rdf?

-ref is the set of positions distances are measured from, and -sel is the set they are measured to. g(r) is averaged over all reference positions. You can pass several -sel selections and get one curve for each in the same output file.

Why does my RDF not go to 1 at long distances?

Either -rmax is too short, the reference is inside a region without bulk solvent (a buried pocket or a membrane core), or you used -surf or -norm none, which do not produce a normalised g(r). Check which case applies before changing anything.

How do I get the coordination number in GROMACS?

Add -cn rdf_cn.xvg to your gmx rdf command. Find the first minimum of g(r) in rdf.xvg, then read the value of rdf_cn.xvg at that distance. That value is the average number of selected particles in the first shell around each reference position.

What bin width should I use with gmx rdf?

The default -bin 0.002 nm works for abundant pairs such as water-water. For one ligand or a few ions as reference, use about 0.01 nm to cut noise. If the peak shape changes when you halve the bin, you do not have enough data for that resolution.

Can gmx rdf use an index file instead of selections?

Yes. Pass -n index.ndx and refer to groups by name, for example -ref 'group LIG_O3'. Selections typed directly on the command line do the same job without an extra file and are easier to record in a methods section.

If you want to learn trajectory analysis on a real protein-ligand system with feedback on every step, the cohort below is built for that.