How to Calculate a Binding Free Energy with Umbrella Sampling in GROMACS (Pulling, Windows and gmx wham)

Umbrella sampling computes a potential of mean force (PMF) along one reaction coordinate. In GROMACS it runs in four stages: pull the ligand away from the protein, extract frames as windows, run a restrained simulation in each window, then combine them with gmx wham. The binding free energy is the depth between the bound minimum and the dissociated plateau.

Most students meet free energy through MM-PBSA with gmx_MMPBSA, which reuses one equilibrium trajectory and returns a number the same afternoon. Then a reviewer, an examiner or a supervisor asks the awkward question: how do you know the ligand does not simply refuse to leave on your timescale, and what did your end-point estimate do about the entropy of unbinding? Umbrella sampling is the standard answer, and it is a different kind of commitment. It is not one simulation. It is a set of simulations arranged along a path, which is why it belongs in a project plan rather than in an afternoon.

This guide is written by the StemSkills Lab team from 10+ years in structural bioinformatics, drug design and multiscale molecular modeling. Every keyword and flag below is checked against the GROMACS mdp options reference and the gmx wham help page for GROMACS 2026.3, and against the official GROMACS umbrella sampling tutorial, which targets GROMACS 2021 or later. The .mdp accepts both pull-coord1-k and pull_coord1_k style spellings of the pull keywords. This article prints the hyphenated form throughout, and you should pick one style and stay with it.

When should you use umbrella sampling instead of MM-PBSA?

Use it when the question is about a free energy along a path, not about a single bound state. Three situations justify the cost.

• You need a defensible binding free energy for one system. A PMF that flattens at large separation gives you a well depth relative to a genuinely dissociated reference state, with sampling in between rather than assumed.
• The interesting physics is a barrier. Residence time arguments, a gate that must open, a ligand that must slide through a channel. Plain MD will not cross these, which is exactly why our post on how to build a free energy landscape in GROMACS names umbrella sampling as the method you move to when the barrier is never sampled.
• You are comparing two states of the same system. Wild type against mutant, or the same ligand in two pockets, along the same reaction coordinate.

Do not use it to rank a library. Every additional ligand is a whole new set of windows, so the cost scales in a way that screening cannot absorb. That is a job for docking, and then for an end-point method.

The three are not interchangeable and they do not answer the same question, so pick by question first and by budget second. If you are working out what to learn in what order, our computational biology skills roadmap puts free energy work after you can already run and analyse a stable trajectory, for a reason: every failure mode below is a sampling failure wearing a free energy costume.

What is the reaction coordinate, and why does choosing it badly ruin everything?

The reaction coordinate is the single geometric quantity you will restrain, scan and plot the free energy against. For protein and ligand work it is usually the distance between the centre of mass of the ligand and the centre of mass of the protein, or of a defined set of pocket residues.

Everything downstream inherits this choice. A PMF is a free energy projected onto your coordinate, with every other degree of freedom averaged over. If the real unbinding event requires a loop to open and your coordinate says nothing about that loop, the loop will not open during your windows, and your profile will describe a process that does not happen. Choose the coordinate before you start, write down why, and be ready to defend it. Three practical rules:

• Pull along a clear exit. Inspect the pocket and pick a direction where the ligand leaves without shearing through the protein. If a buried ligand has no straight exit, a distance coordinate is the wrong tool.
• Use centre of mass groups, not single atoms. The GROMACS pull code is built on group centres of mass, and a single-atom coordinate couples your free energy to one atom’s fluctuations.
• Keep the pocket definition fixed. If group 1 is a residue selection, use the same index group for the pull run, every window and the analysis.

How do you prepare the system?

Exactly as for any other production run. Build the topology, solvate, add ions, minimise, then run NVT and NPT equilibration. Nothing about umbrella sampling changes the setup, so follow the standard protein and ligand MD workflow in GROMACS and the NVT and NPT equilibration guide rather than repeating them here.

One setup decision is specific to this method and it is the one that quietly ruins runs: the box. You are going to drag the ligand a long way along one axis, and the ligand must never approach its own periodic image. Build the box elongated along the pull axis with enough solvent that the fully separated state still has a real distance to the nearest image, using gmx editconf with explicit box vectors rather than a cubic default. Deciding this after the pull run means running it again.

How do you run the pull simulation?

The pull run is a normal MD run with the pull code switched on and a non-zero pull rate, so the reference distance moves and the ligand is dragged out. Its only job is to generate a continuous trajectory of configurations spanning the whole coordinate, from bound to dissociated. It is not a free energy calculation, and its work values are not your answer.

The pull block goes in a standard .mdp file next to your usual integrator, thermostat and barostat settings. If any of those are unfamiliar, our post on the GROMACS .mdp file explained covers the rest of the file. The pull section looks like this:

pull                     = yes
pull-ngroups             = 2
pull-ncoords             = 1
pull-group1-name         = Protein          ; index group, your pocket or whole protein
pull-group2-name         = LIG              ; index group, the ligand
pull-coord1-type         = umbrella
pull-coord1-geometry     = distance
pull-coord1-groups       = 1 2
pull-coord1-dim          = N N Y            ; pull along z only, in this example
pull-coord1-rate         = RATE             ; nm ps^-1, YOU choose this, see below
pull-coord1-k            = K                ; kJ mol^-1 nm^-2, YOU choose this
pull-nstxout             = 1000
pull-nstfout             = 1000

RATE and K are placeholders. GROMACS will not accept them, and no honest article can hand you a number here, because the right values depend on your pocket, your ligand and how far you need to travel. What the documentation does fix is the units: pull-coord1-k is in kJ mol^-1 nm^-2 and pull-coord1-rate is in nm ps^-1. Two decisions are yours:

• The rate. Slower is better and slower is more expensive. Too fast and you rip the ligand out through the protein instead of along the exit, leaving distorted structures that no amount of umbrella sampling repairs. Watch the pulled trajectory before you trust it.
• The force constant. Strong enough that the reference distance actually drags the ligand, not so strong that the spring dominates the dynamics.

Read the entry for pull-coord1-start in the mdp reference before you use it, because it changes how pull-coord1-init is interpreted relative to your starting structure. Also set pull-nstxout and pull-nstfout explicitly. These write the pullx and pullf files, and gmx wham reads exactly those files later.

The type stays umbrella for both stages. What distinguishes the pulling stage from the production windows is the rate, and this is the single point beginners get wrong. The pull run uses a rate greater than zero. Every window uses a rate of exactly 0.0, so the restraint sits still at its assigned distance.

The harmonic restraint itself is what the GROMACS reference manual describes for umbrella pulling: “A harmonic potential is applied between the centers of mass of two groups” (GROMACS reference manual, pull code). That is the entire physical content of the method, and it dates to Torrie and Valleau’s 1977 paper in the Journal of Computational Physics, volume 23, pages 187 to 199, doi:10.1016/0021-9991(77)90121-8, which introduced biasing potentials to force sampling of regions ordinary simulation ignores.

How do you choose the umbrella windows?

A window is a starting structure plus the reference distance its restraint will hold. You harvest both from the pull trajectory. The procedure has three steps and none of them involves guessing.

Step 1: measure what you actually have. Do not assume the pull produced evenly spaced configurations, because a constant reference velocity does not give a constant centre-of-mass separation. Measure the real separation frame by frame with gmx distance, using a selection built from the centre of mass of each group. The gmx distance help page documents the selection options, including the reference-position types that give you a centre of mass rather than an atom.

Step 2: pick target distances, then pull the matching frames. Walk your measured distance list from the bound minimum out to a separation where the ligand is clearly free, and choose the frames whose real distances give you the spacing you want. Extract each one with gmx trjconv and its -dump option, which extracts a frame at or near one specific time:

gmx trjconv -s pull.tpr -f pull.xtc -o conf_window12.gro -dump 187

Step 3: set each window’s restraint to its own distance. Each window gets its own .mdp with pull-coord1-init equal to that window’s measured separation and pull-coord1-rate = 0.0. The official tutorial’s production block, for its two-pyrimidine system, is a working example of the shape:

pull                     = yes
pull-ngroups             = 2
pull-group1-name         = pyrimidine_1
pull-group2-name         = pyrimidine_2
pull-ncoords             = 1
pull-coord1-groups       = 1 2
pull-coord1-geometry     = distance
pull-coord1-type         = umbrella
pull-coord1-k            = 5000.0
pull-coord1-rate         = 0.0
pull-coord1-init         = 0.834

Treat 5000.0 and 0.834 as that tutorial’s values for that system, not as defaults. The -init value is per window by definition. The force constant of 5000 kJ mol^-1 nm^-2 is a reasonable order of magnitude to start from for a distance coordinate, and whether it suits your system is decided by one test, not by argument.

What is the acceptance criterion for a set of windows?

Histogram overlap. Every window produces a distribution of sampled distances, and neighbouring distributions must overlap. This is not a stylistic preference. WHAM reconstructs one continuous profile by stitching windows together through their shared regions, so where there is no overlap, there is nothing to stitch, and the free energy difference across that gap is unconstrained.

The test is visual and immediate. Run gmx wham and open the histogram file it writes. If the bells sit shoulder to shoulder with visible overlap all the way from bound to plateau, the window set is acceptable. If any pair barely touches, add a window between them. Two failure shapes tell you which knob to turn:

• Spikes with no overlap: your force constant is too strong. Each window is pinned so tightly it explores almost nothing.
• Broad, heavily merged bells: your force constant is too weak, and neighbouring windows are sampling the same region, which wastes compute without extending coverage.

Fixing a gap by raising the force constant everywhere is the wrong move and makes the overlap worse. Add windows in the gap.

How do you run the windows?

Each window is an independent simulation: grompp it, then mdrun it. Equilibrate each window briefly under its own restraint before collecting production data, so the ligand settles at its assigned distance rather than relaxing during the data you keep.

The useful property here is that the windows do not talk to each other. They are embarrassingly parallel, so on a shared cluster you submit them as an array job and they queue independently instead of blocking each other. This is what makes umbrella sampling feasible on modest resources: the wall-clock cost is one window, repeated, not the sum of a serial chain. Script the loop over windows rather than typing commands by hand, because you will run this more than once.

How do you run gmx wham and read profile.xvg?

gmx wham is GROMACS’s implementation of the Weighted Histogram Analysis Method, the method published by Kumar, Bouzida, Swendsen, Kollman and Rosenberg in the Journal of Computational Chemistry, 1992, volume 13, issue 8, pages 1011 to 1021, doi:10.1002/jcc.540130812. The tool itself is described in Hub, de Groot and van der Spoel, “g_wham, a free weighted histogram analysis implementation including robust error and autocorrelation estimates”, Journal of Chemical Theory and Computation, 2010, volume 6, issue 12, pages 3713 to 3720, doi:10.1021/ct100494z. Cite both if you publish a PMF from GROMACS.

It takes lists of files, not the files themselves. You write two plain text files, one listing every window’s .tpr with -it and one listing every window’s pull output with -if for forces or -ix for coordinates:

gmx wham -it tpr-files.dat -if pullf-files.dat -o profile.xvg -hist histo.xvg -unit kJ

Three details matter. The two lists must be in the same order, line for line, because gmx wham pairs them by position. The default output names are profile.xvg for the PMF and histo.xvg for the histograms. The default energy unit is kJ, changed with -unit (kJ, kCal or kT), and the default analysis temperature is 298 K, changed with -temp. If your simulations ran at 310 K and you leave -temp alone, the numbers are wrong in a way nothing else will flag.

profile.xvg is two columns: the reaction coordinate and the free energy. Plot it as you would any other GROMACS output, for example with the approach in our guide to plotting .xvg files in Python. Read it as follows: the minimum near the bound distance is the bound state, the flat region at large separation is the dissociated reference, and the vertical difference between them is your binding free energy. If the profile has not flattened at your largest distance, you did not pull far enough, and the depth you would quote is a number relative to a state that is still interacting.

How do you put error bars on the PMF?

By bootstrapping, which gmx wham does for you. Add -nBootstrap with the number of bootstrap iterations and select the method with -bs-method, whose allowed values are b-hist (the default), hist, traj and traj-gauss. The results are written to bsResult.xvg, the bootstrapped average with standard deviations, and bsProfs.xvg, every individual bootstrapped profile. The error estimation described in the Hub 2010 paper is precisely what this machinery is for, so a PMF published without it invites the first question at your viva.

How do you know the PMF is converged?

Split and compare. Re-run gmx wham using only the first half of each window’s data, then again using only the second half, and overlay the two profiles on the finished one. If the well depth and the shape move appreciably between halves, the windows are not converged and the answer is more sampling per window, not more windows. This is the honest test, it costs nothing but analysis time, and it is the one examiners ask about.

What breaks a PMF, and how do you fix it?

• The PMF has an unphysical jump. Look at histo.xvg: there is a gap between two histograms. Add windows inside the gap. Do not raise the force constant globally.
• Every histogram is a narrow spike. The force constant is too high. Lower it and re-run the affected windows.
• Neighbouring histograms are nearly identical. The force constant is too low, so windows are not held apart. Raise it, or space the windows further and let overlap decide.
• The profile never flattens. The pull did not travel far enough for a dissociated plateau to exist. Extend the pull and add windows at the far end. If the box cannot accommodate that, rebuild the box first.
• The ligand meets its own periodic image. The box was too short along the pull axis. Rebuild with a longer vector along that axis and start again. Nothing salvages this after the fact.
• The protein rotates or drifts during the pull. Your distance coordinate stops meaning what you intended. Check the centre-of-mass motion removal setting, and consider restraining the protein, using a pocket-residue pull group instead of the whole protein, or a geometry that fixes the pull direction rather than a free distance.
• gmx wham runs but the profile is nonsense. Check that the .tpr list and the pull-file list are in the same order. A mismatched pairing produces output without producing an error.
• gmx wham cannot find the pull data. Your windows ran without pull-nstxout or pull-nstfout set, so there is nothing to analyse. Re-run the windows with the output frequencies set.
• Something fails at grompp instead. Pull group names must match index groups exactly, and the same index file must be passed to every window. Our list of common GROMACS errors and how to fix them covers the general grompp failures.

How do you report the result in a thesis or paper?

State, in the methods, every one of the following: the reaction coordinate and the two pull groups by definition; the pull rate and force constant used for the pulling stage; how windows were selected and how many there were; the umbrella force constant and the simulation time per window; the total sampling summed over windows; the WHAM implementation with its citation; the analysis temperature; and the bootstrap method and number of iterations. In the results, give the well depth with its bootstrap uncertainty, show the PMF with the histogram overlap plot next to it, and state the convergence test you ran.

Reviewers reject PMFs for missing exactly these items, so write the methods paragraph while the runs are still on disk. Our guide to writing a molecular dynamics methods section for GROMACS has the structure, and the same discipline applies here with the pull and WHAM parameters added.

Frequently asked questions

Is umbrella sampling more accurate than MM-PBSA?

It answers a different question. Umbrella sampling gives a free energy profile along a coordinate you chose, with explicit solvent and with the entropy along that coordinate sampled. MM-PBSA gives an end-point estimate from bound-state snapshots. A PMF from a badly chosen coordinate or from unconverged windows is not better than anything. The comparison table above compares them on what they compute, not on a ranking.

How many windows do I need?

As many as it takes to get overlapping histograms from the bound minimum to a flat plateau. There is no correct count independent of your system, your force constant and how far you pulled. Choose spacing from real measured distances, run gmx wham, look at histo.xvg, and add windows where the overlap fails.

Should I use pullf or pullx files for gmx wham?

Either. Pass force files with -if or coordinate files with -ix. Whichever you choose, the file list must be in the same order as the .tpr list you pass to -it.

Can I run umbrella sampling on a laptop or on Colab?

The windows are independent, so you can technically run them one at a time anywhere GROMACS runs, including our GROMACS on Google Colab setup. Whether it finishes in project time is a different question, because the total is the cost of one window multiplied by the number of windows. Plan it as cluster work.

Why does my binding free energy come out with the wrong sign?

Check which end of the profile you subtracted from which. The PMF is conventionally plotted with free energy against separation, so the bound state is the minimum and the dissociated plateau is the reference. The binding free energy is the plateau minus the minimum, and reversing the two flips the sign.

Where should I learn the surrounding workflow?

Start from the pillar guide, learn molecular dynamics with GROMACS, which sequences setup, equilibration, production and analysis, then come back here once your unbiased trajectories are already stable. The official GROMACS umbrella sampling tutorial is the reference implementation to practise on before you attempt your own system.