Back

Computational Biochemistry · Free Energy Methods

SMD & Umbrella Sampling Simulation

Steered MD pulling · Umbrella sampling windows · WHAM analysis · PMF & ΔGbind · GROMACS

GROMACS GROMOS96 53A6 Steered MD Umbrella Sampling WHAM · PMF · ΔG
System
Prep
🧬
Prepare Topology
Step 1 pdb2gmx
📦
Define Unit Cell
Step 2 editconf
💧
Solvation & Ions
Step 3 solvate
📉
EM & Equilibration
Step 4 mdrun
SMD ·
US · PMF
🏹
Steered MD
Step 5 pull code
☂️
Umbrella Sampling
Step 6 23 windows
📊
WHAM & PMF
Step 7 gmx wham
✅
Summary
Step 8 ΔGbind
1
System Preparation
🧬
Prepare the Topology pdb2gmx
Aβ42 protofibril (2BEG) · GROMOS96 53A6 · N-acetylated termini · Chain B position restraints
PDB 2BEG GROMOS96 53A6 SPC water POSRES_B COO⁻ termini
▼

System Overview

The system is the dissociation of a single peptide (chain A) from the growing end of an Aβ42 protofibril. The structure file is the wild-type Aβ42 protofibril (PDB: 2BEG), acetylated at the N-terminus of each chain. The reaction coordinate is the z-axis COM distance between chain A and chain B.

Generate Topology

gmx pdb2gmx -f 2BEG_model1_capped.pdb -ignh -ter -o complex.gro
# Select: GROMOS96 53A6 · SPC water · "None" N-termini · "COO-" C-termini

Add Chain B Position Restraints

Chain B will serve as the immobile reference during the pulling simulation. Add a special position restraint block to topol_Protein_chain_B.itp:

#ifdef POSRES_B
#include "posre_Protein_chain_B.itp"
#endif
Why restrain chain B? Without restraining chain B, the extensive non-covalent interactions between chains A and B would cause the entire complex to be towed through the simulation box rather than separating chain A from the protofibril.
Force fieldWater modelN-terminusC-terminus
GROMOS96 53A6SPCNone (acetylated)COO⁻
📦
Define the Unit Cell editconf
Elongated box along z-axis · 12.0 nm z-dimension · satisfies minimum image convention during pull
6.56 × 4.36 × 12.0 nm z-axis pull Min. image conv. 5.0 nm pull dist.
▼

Critical Requirement for Pulling Simulations

The minimum image convention must be satisfied throughout the entire pulling simulation. GROMACS calculates distances considering periodicity — if you pull over a distance greater than half the box length in the pulling direction, the periodic image distance becomes the reference and completely corrupts the results.

Rule: Pull distance must always be less than ½ the box length in the pulling direction. Here: pull distance = 5.0 nm, box z = 12.0 nm → 5.0 < 6.0 ✓

Place Protofibril in the Box

gmx editconf -f complex.gro -o newbox.gro -center 3.280 2.181 2.4775 -box 6.560 4.362 12

The protofibril is centered such that the peptide to be pulled has sufficient space to travel 5.0 nm along the z-axis without encountering periodic images.

Box x
6.560 nm
Box y
4.362 nm
Box z
12.000 nm
Max pull dist.
5.0 nm
💧
Solvation & Ions solvate + genion
SPC water · 100 mM NaCl physiological concentration · neutralizing counterions
SPC water 100 mM NaCl -neutral -conc 0.1 NA + CL
▼

Add Water

gmx solvate -cp newbox.gro -cs spc216.gro -o solv.gro -p topol.top

Add Ions at Physiological Concentration

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.1
# Select group 13 (SOL) to replace water molecules

Unlike the basic protein simulation, umbrella sampling uses physiological ionic strength (100 mM NaCl) on top of the neutralizing counterions. This is important for accurately capturing the binding free energy in a realistic environment.

Water model
SPC
Salt conc.
100 mM NaCl
NA⁺ (neutralize)
+ physiological
CL⁻
+ physiological
📉
Energy Minimization & Equilibration mdrun EM + NPT
Steepest descent EM · NPT equilibration with position restraints · prepare system for pulling
Steepest descent NPT equil. POSRES_B npt.cpt
▼

Energy Minimization

gmx grompp -f minim.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

NPT Equilibration

gmx grompp -f npt.mdp -c em.gro -p topol.top -r em.gro -o npt.tpr
gmx mdrun -deffnm npt

The -r em.gro flag applies position restraints to the protein heavy atoms. Chain B restraints (POSRES_B) ensure the reference fibril structure remains stable while the system equilibrates thermally and mechanically.

Potential Energy — Energy Minimization
Aβ42 Protofibril Complex · Steepest Descent · GROMOS96 53A6
Potential energy converges — system ready for steered MD pulling
2
Steered MD & Umbrella Sampling
🏹
Generating Configurations via Steered MD pull code
500 ps SMD pull · 0.01 nm/ps rate · extract 501 frames · measure COM distance · select ~23 windows
z-axis pull 500 ps k = 1000 kJ/mol/nm² 0.2 nm spacing pullf.xvg
▼

Umbrella Sampling Schematic

The reaction coordinate ζ is the z-axis COM distance between chain A and chain B. A steered MD run generates a continuous trajectory of configurations spanning the full reaction coordinate, from which individual window starting configurations are extracted.

Create Index Groups

gmx make_ndx -f npt.gro
> r 1-27
> name 19 Chain_A
> r 28-54
> name 20 Chain_B
> q

Pull Code .mdp Settings

pull = yes pull_ncoords = 1 ; one reaction coordinate pull_ngroups = 2 ; two groups define the coordinate pull_group1_name = Chain_A pull_group2_name = Chain_B pull_coord1_type = umbrella ; harmonic biasing potential pull_coord1_geometry = distance ; COM distance pull_coord1_dim = N N Y ; pull along z only pull_coord1_groups = 1 2 pull_coord1_start = yes ; initial COM distance = reference pull_coord1_rate = 0.01 ; 0.01 nm/ps = 10 nm/ns pull_coord1_k = 1000 ; kJ mol⁻¹ nm⁻²
  • pull_coord1_rate = 0.01 — constant velocity pulling; force builds until critical non-covalent interactions with chain B are broken
  • pull_coord1_dim = N N Y — pulling restricted to z-axis only; x and y motion is unrestricted
  • pull_coord1_k = 1000 — spring force constant; must not be so large as to deform system elements
  • POSRES_B — restraining chain B prevents dragging the whole fibril instead of separating chain A

Run Steered MD

gmx grompp -f md_pull.mdp -c npt.gro -p topol.top -r npt.gro -n index.ndx -t npt.cpt -o pull.tpr
gmx mdrun -deffnm pull -pf pullf.xvg -px pullx.xvg
Pull Force vs Time — Steered MD
Chain A dissociation from Aβ42 protofibril · 500 ps
Force builds then drops as critical contacts break — chain A begins dissociating ~halfway through

Extract Frames & Measure COM Distances

gmx trjconv -s pull.tpr -f pull.xtc -o conf.gro -sep
# Produces conf0.gro ... conf500.gro (501 frames at 1 ps intervals)
bash get_distances.sh # outputs summary_distances.dat

Review summary_distances.dat and identify frames with ~0.2 nm COM distance spacing for umbrella windows. Example:

6 0.500 ← use conf6.gro as window 1 ... 160 0.704 ← use conf160.gro as window 2 ... 449 4.900 ← use conf449.gro as window 23
NOT umbrella sampling yet: This SMD step only generates starting configurations. The actual umbrella sampling (restrained at fixed COM distance) is the next step. The pull rate here must be validated — too fast causes system deformation.
☂️
Umbrella Sampling Simulations md_umbrella
23 windows · 0.2 nm spacing · 0.5–5.0 nm range · NPT equilibration per window · rate = 0
23 windows 0.5–5.0 nm rate = 0 NPT/window umbrella*.tpr
▼

Step 1 — NPT Equilibration in Each Window

Each selected starting configuration must be briefly equilibrated in its window before data collection. This stabilizes the local configuration at the target COM distance.

gmx grompp -f npt_umbrella.mdp -c conf6.gro -p topol.top -r conf6.gro -n index.ndx -o npt0.tpr
gmx grompp -f npt_umbrella.mdp -c conf449.gro -p topol.top -r conf449.gro -n index.ndx -o npt22.tpr
...
gmx mdrun -deffnm npt0
gmx mdrun -deffnm npt22

Key Difference from SMD — pull_coord1_rate = 0

In umbrella sampling, pull_coord1_rate is set to zero. The spring does not move. The system is harmonically restrained at the initial COM distance of each window. Setting pull_coord1_start = yes means the reference distance is automatically read from each starting configuration — no manual specification needed.

Step 2 — Run Umbrella Sampling MD

gmx grompp -f md_umbrella.mdp -c npt0.gro -t npt0.cpt -p topol.top -r npt0.gro -n index.ndx -o umbrella0.tpr
...
gmx grompp -f md_umbrella.mdp -c npt22.gro -t npt22.cpt -p topol.top -r npt22.gro -n index.ndx -o umbrella22.tpr

gmx mdrun -deffnm umbrella0
...
gmx mdrun -deffnm umbrella22
Windows
23
COM range
0.5 – 5.0 nm
Spacing
~0.2 nm
pull rate
0 (fixed)
Umbrella Sampling Windows — COM Distance Histograms
23 windows · 0.5–5.0 nm · Neighboring windows must overlap for valid PMF
Overlapping Gaussian distributions confirm adequate sampling — gaps indicate windows needing additional simulation
Overlap requirement: Neighboring histograms must overlap for WHAM to reconstruct a continuous PMF. If a gap is found (e.g., at ζ = 0.8 nm), add an extra window centered there and re-run WHAM — the other windows do not need to be repeated.
3
WHAM Analysis & Results
📊
WHAM Analysis & PMF Extraction gmx wham
Weighted Histogram Analysis Method · reconstruct PMF curve · ΔG = 46 kcal/mol · kCal or kJ/mol output
gmx wham PMF curve ΔGbind = 46 kcal/mol pullf-files.dat tpr-files.dat
▼

Prepare Input File Lists

# tpr-files.dat:
umbrella0.tpr
umbrella1.tpr
...
umbrella22.tpr

# pullf-files.dat (files must have unique names):
umbrella0_pullf.xvg
umbrella1_pullf.xvg
...
umbrella22_pullf.xvg

Run WHAM

gmx wham -it tpr-files.dat -if pullf-files.dat -o -hist -unit kCal

WHAM opens each umbrella*.tpr and umbrella*_pullf.xvg sequentially and reconstructs the unbiased PMF from the biased sampling in each window. Output units can be kCal, kJ, or kT.

PMF Curve & ΔGbind

Potential of Mean Force (PMF)
Aβ42 Chain A Dissociation · WHAM reconstruction · ΔG = 46 kcal mol⁻¹
ΔGbind = difference between PMF maximum and plateau at large ζ. Reported value: 46 kcal mol⁻¹ (published: 50.5 kcal mol⁻¹)
ΔGbind (tutorial)
46 kcal/mol
ΔGbind (published)
50.5 kcal/mol
PMF plateau
> 3.5 nm
Energy units
kcal mol⁻¹
PMF defect at ζ ≈ 0.8 nm: A defect in the PMF reflects insufficient sampling at a high-energy transition state region. Fix: add one extra umbrella window centered at ζ = 0.8 nm and re-run WHAM — no other windows need repeating. This is the power of independent umbrella windows.
✅
Summary & Best Practices ΔGbind
Complete umbrella sampling workflow · key considerations · common pitfalls · validation
PMF validation Convergence check Window overlap Bootstrap error
▼

Complete Workflow Summary

StepToolKey Output
1. Topologypdb2gmxtopol.top · complex.gro
2. Unit celleditconfnewbox.gro (12 nm z)
3. Solvationsolvate + genionsolv_ions.gro (100 mM NaCl)
4. EM + NPTmdrunnpt.gro · npt.cpt
5. Steered MDpull codeconf0-500.gro · pullf.xvg
6. US windowsmd_umbrellaumbrella0-22.tpr · pullf.xvg
7. WHAMgmx whamPMF · histograms · ΔGbind

Critical Settings & Pitfalls

  • Box length rule: pull distance < ½ box length in pulling direction — always
  • Pull rate validation: too fast deforms the system; validate using published literature or trial simulations
  • Window overlap: histograms must overlap in neighboring windows; gaps produce PMF artefacts
  • Rate = 0 in US: pull_coord1_rate = 0 is mandatory during umbrella sampling — the window must not drift
  • POSRES_B reference: restraining the immobile group prevents rigid-body drift of the whole complex
  • Thermostat coupling: do not couple small groups (JZ4, CL) separately — unstable kinetic energy fluctuations
  • Bootstrap error: run gmx wham -bs to compute statistical error on the PMF via bootstrapping

PMF Interpretation

ΔGbind = PMF maximum − PMF plateau at large ζ (fully dissociated state). This represents the free energy required to dissociate the peptide from the fibril. The PMF must converge to a stable flat value at large separations for the result to be physically meaningful.
⚙
Simulation Parameters
Force field
GROMOS96 53A6
Water model
SPC
Salt conc.
100 mM NaCl
Box z-length
12.0 nm
Pull rate (SMD)
0.01 nm/ps
Force constant
1000 kJ/mol/nm²
US windows
23
Window spacing
~0.2 nm
📈
Key Results
ΔGbind
46 kcal mol⁻¹
binding free energy
PMF
Potential of mean force
along z-axis
23
Umbrella windows
0.5–5.0 nm
WHAM
Weighted histogram
analysis method
pullf.xvg
Force vs time
from each window
Histograms
Window overlap
convergence check
System: Aβ42 Protofibril (2BEG)  ·  Force field: GROMOS96 53A6  ·  Method: SMD + Umbrella Sampling + WHAM
Reference: J. Lemkul & D. Bevan, J. Phys. Chem. B (2010) · Virginia Tech · mdtutorials.com/gmx/umbrella