Force–Extension Analysis in Steered Molecular Dynamics

Simple concept

In constant-velocity pulling, a virtual spring moves at a prescribed speed. The molecular system may lag behind that moving reference because of binding contacts, steric barriers, solvent resistance, protein deformation, membrane resistance, or internal rearrangement. The spring extension produces the recorded pulling force.

Reaction coordinate ξ may be a distance, projection along a direction, angle, dihedral, or another supported collective variable. In ligand unbinding it is often the separation between ligand and protein reference groups.

Why it is often more mechanistically useful

Force vs distance reveals where along the mechanical extension the system resisted pulling. It is especially useful for ligand unbinding, domain separation, membrane permeation, peptide detachment, protein unfolding, and rupture of macromolecular contacts.

• Peak position identifies the extension at which the largest resistance occurred.
• Area under the curve gives mechanical work along the sampled pulling path.
• Multiple peaks can represent sequential rupture of contacts or structural intermediates.
• Comparison between systems is only meaningful under matched protocols.

Meaning

Force vs position maps the force directly onto the collective variable ξ. For a straight directional pull, the x-, y-, or z-projected coordinate may be used. For a distance geometry, the reaction coordinate may simply be group separation.

This representation is valuable for choosing umbrella-sampling windows because it shows coordinate regions where resistance changes, structural transitions occur, or denser window spacing may be necessary.

Conceptual equation

W(ξ) = ∫ F(ξ) dξ

The signed area under a force–coordinate curve gives the mechanical work performed by the pulling force along that path. In discrete data, work is often estimated numerically with a trapezoidal integration.

• Work depends on the exact pulling path and protocol.
• Fast pulling usually produces larger dissipative work.
• A single nonequilibrium trajectory does not normally provide an equilibrium free energy.
• Repeated nonequilibrium work trajectories may be analyzed with specialized estimators, but convergence can be demanding.

Illustrative pull-code settings

pull                     = yes
pull-ngroups             = 2
pull-ncoords             = 1
pull-group1-name         = Protein
pull-group2-name         = Ligand
pull-coord1-type         = umbrella
pull-coord1-geometry     = direction
pull-coord1-dim          = N N Y
pull-coord1-groups       = 1 2
pull-coord1-vec          = 0 0 1
pull-coord1-start        = yes
pull-coord1-rate         = 0.001
pull-coord1-k            = 1000
pull-nstxout             = 100
pull-nstfout             = 100

Values are examples only. Geometry, dimensions, rate, force constant, reference groups, PBC handling, and direction must be designed for the system.

Run SMD and write force/coordinate files

gmx grompp -f md_pull.mdp -c npt.gro -r npt.gro -p topol.top -n index.ndx -o pull.tpr
gmx mdrun -deffnm pull -pf pullf.xvg -px pullx.xvg

Umbrella sampling and PMF with WHAM

# Prepare text files listing window inputs
# tpr-files.dat  -> one umbrella .tpr path per line
# pullf-files.dat or pullx-files.dat -> matching output files

gmx wham -it tpr-files.dat -if pullf-files.dat -o pmf.xvg -hist histo.xvg -bsres bsResult.xvg -bsprof bsProfs.xvg

Exact options depend on GROMACS version and input mode. Inspect histogram overlap, discard equilibration appropriately, and use bootstrap or block analysis for uncertainty.

Example reporting sentence

“Constant-velocity steered molecular dynamics was performed along a defined protein–ligand separation coordinate using a harmonic pulling restraint. Pulling force and coordinate were recorded throughout the trajectory, and Force–Time and Force–Distance profiles were analyzed across independent replicates. Representative configurations were subsequently used to initialize umbrella-sampling windows for WHAM-based PMF reconstruction.”