Module 4: Molecular Mechanics & Force Fields
Every method that follows — docking scores, MD trajectories, conformational searches — is ultimately calling the same function: a force field that turns a set of atomic coordinates into a single number, the energy. This module is that engine.
1. Atoms and Springs: The Central Approximation
Molecular mechanics (MM) treats a molecule as a set of point masses (atoms) connected by springs (bonds), governed entirely by classical physics. It throws away the electrons.
That sounds like an outrageous simplification — and it is. Quantum mechanics gives the true energy, but its cost scales steeply with system size (formally N⁴ or worse for correlated methods), which puts a solvated protein permanently out of reach. MM scales roughly as N², so it can evaluate a 100,000-atom system in milliseconds. The entire field of computational drug discovery exists inside that trade-off.
What you buy
Speed — enough to score millions of docking poses (Module 8) or integrate a trajectory for microseconds (Module 10). Without MM, none of it is possible.
What you pay
No electrons means no bond breaking or forming, no charge transfer, and no polarization. A classical force field can never model a chemical reaction — which is why covalent docking and enzymology need QM/MM.
2. The Force Field Energy Decomposition
A force field is two things: a set of functional forms for the energy, and the parameters that go in them. The total energy is a simple sum of independent terms:
The split matters enormously. Bonded terms involve atoms 1–4 bonds apart and describe a molecule's internal geometry. Non-bonded terms act between all remaining atom pairs — including atoms on entirely different molecules. Non-bonded terms are where drug binding lives: when a ligand meets a protein, no bonds change, so the entire interaction is vdW plus electrostatics. That is also why they dominate the cost — there are N² of them.
The Bonded Terms
Bond stretching
E = k_b (r − r₀)²A harmonic (Hooke's law) spring. r₀ is the equilibrium bond length, k_b the stiffness. Because it is harmonic, the bond can never break — stretch it far enough and the energy just rises forever, quadratically.
Angle bending
E = k_θ (θ − θ₀)²Also harmonic, penalizing deviation of a three-atom angle from its equilibrium value θ₀ (e.g. ~109.5° for sp³ carbon).
Torsion (dihedral)
E = Σ (Vₙ/2)[1 + cos(nφ − γ)]Periodic, not harmonic — because rotating a bond by 360° must return you to where you started. n is the periodicity (n = 3 for the three-fold barrier of an sp³–sp³ bond). This is the softest term and therefore the one that governs conformation.
Improper / out-of-plane
E = k_ω (ω − ω₀)²A correction term with no direct physical bond behind it: it exists to keep groups planar that ought to be planar — aromatic rings, sp² centres, amide bonds — which the other three terms would otherwise allow to pucker.
The Non-Bonded Terms
Van der Waals — Lennard-Jones 12-6
E = 4ε[(σ/r)¹² − (σ/r)⁶]You have already met this curve in Module 3: an r⁻⁶ attractive tail (London dispersion) against a brutally steep r⁻¹² repulsive wall (Pauli exclusion). The r⁻¹² term is why steric clashes are so catastrophically expensive — halve the distance and the penalty rises 4,096-fold.
Electrostatics — Coulomb
E = q_i·q_j / (4πε₀ε_r · r)Fixed partial charges on each atom, decaying only as 1/r — long-range, which is why it is the expensive term to compute properly. Because the charges are fixed, a standard force field cannot let a molecule polarize in response to its environment.
Interactive Playground: The Butane Torsional Profile
The single most important curve in conformational analysis. Rotate the central C–C–C–Cdihedral of butane and watch the Newman projection turn while the torsion term reports the energy. The numbers are not illustrative — this is a real OPLS Fourier torsion, and it reproduces butane's textbook profile.
dark = front carbon · grey = rear carbon
0.00
Syn — fully eclipsed CH₃/CH₃
- Rotatable bonds are cheap to rotate. The barriers here are only ~3–5 kcal/mol — comparable to thermal energy — so at 298 K butane is constantly interconverting. That is precisely why a flexible ligand has so many conformers to search (Module 6) and why each rotatable bond costs entropy on binding (Module 3).
- The bioactive conformation is rarely the global minimum. Anti dominates in solution (~68%), but a receptor can pay a kcal or two to bind a gauche-like shape. A conformer search that only keeps the global minimum will miss it.
- Note the two-to-one degeneracy. Gauche sits 0.86 kcal/mol above anti yet still accounts for ~32% of molecules — because there are two gauche wells and only one anti. Energy alone never determines population; you must count states.
3. Force Field Families
The functional forms above are nearly universal. What distinguishes force fields is their parameters and what they were fitted to reproduce.
| Force field | Built for | Notes |
|---|---|---|
| AMBER | Proteins, nucleic acids | Charges fitted to reproduce the QM electrostatic potential (RESP). ff19SB is the current protein standard. |
| CHARMM | Proteins, lipids, membranes | CGenFF extends it to drug-like molecules. Strong lipid/membrane parameters. |
| OPLS | Organic liquids, small molecules | Parameterized against experimental liquid densities and heats of vaporization — the torsion you just used is OPLS-style. |
| GROMOS | Biomolecular simulation | Historically united-atom (hydrogens absorbed into heavy atoms) for speed. |
| GAFF | Arbitrary drug-like ligands | The general AMBER force field — the usual answer to "my ligand has no parameters". |
Force fields assume an atom type behaves the same everywhere: a carbonyl carbon in one molecule is parameterized identically to a carbonyl carbon in another. Proteins are made of 20 repeating amino acids, so this works beautifully. Your ligand is the problem.It is a novel molecule nobody parameterized, so it gets assigned by analogy via GAFF/CGenFF — and unusual chemotypes get poor parameters. A meaningful fraction of "the docking score was wrong" is really "the ligand parameters were wrong".
4. Energy Minimization & the Multiple-Minimum Problem
A force field gives energy as a function of coordinates — the potential energy surface (PES). Minimization walks downhill on that surface to find a stable geometry, using the gradient (the forces):
Steepest descent
Step directly downhill along the negative gradient. Robust and forgiving of terrible starting structures, but converges slowly near the minimum. The standard first pass for relieving clashes after adding hydrogens or building a homology model.
Conjugate gradient / L-BFGS
Use information from previous steps to choose smarter directions. Much faster convergence, but less tolerant of a bad start. The usual second pass.
This is the defining limitation. Minimization finds the nearest local minimum — the bottom of whatever valley you happened to start in — and can never cross a barrier to find a better one. In the butane playground above, start at 150° and minimization will slide you to anti; start at 90° and it will trap you in gauche. It has no idea the other well exists.
A drug-like molecule with 5 rotatable bonds has hundreds of local minima. Escaping them requires sampling, not minimization: systematic or stochastic conformer searches (Module 7), Monte Carlo and genetic algorithms (Module 6), or molecular dynamics with enough thermal energy to hop barriers (Module 10). Minimization is how you polish a structure; it is never how you find one.
5. Where Classical Force Fields Break
Knowing the failure modes is what separates using a force field from trusting it blindly.
No reactivity
Harmonic bonds cannot break. Any covalent inhibitor, any enzyme mechanism, any proton transfer is outside the model. Fix: QM/MM — treat the reacting centre quantum-mechanically and the rest classically (Module 10).
Fixed charges
A real molecule's electron density rearranges when it enters a charged pocket. Fixed point charges cannot. Fix: polarizable force fields (AMOEBA, Drude) — more accurate, considerably more expensive.
Parameter quality
The model is only as good as the numbers in it, and novel chemotypes get parameters by analogy. Fix: machine-learned force fields, which learn the PES from QM data directly instead of assuming a functional form — covered in Module 10.
Each of these limitations is the reason a later technique exists. Keep the master equation in mind for the rest of the course: when a docking score misleads you (Module 6) or a trajectory drifts (Module 10), the explanation is almost always hiding in one of these six terms.