Lesson 4.2 · 15 min

Dihedrals

A dihedral (torsion) angle describes the twist around a bond, defined by four atoms. Its energy, a short sum of cosines, controls the shapes molecules adopt: trans or gauche in butane, the backbone of proteins, the flatness of rings.

The idea

Bonds (two atoms) and angles (three atoms) fix distances and bends. The last internal coordinate is the twist around a bond. Take four bonded atoms 1–2–3–4. Atoms 1, 2, 3 define one plane, and 2, 3, 4 another. The dihedral angle φ is the angle between these planes. Looking straight down the 2–3 bond (a Newman projection), φ is the angle between atom 1 in front and atom 4 behind. By convention φ = 180° is trans (anti), ±60° gauche, and 0° cis.

Torsion energy. Rotating around a single bond is periodic, so the energy is a short Fourier series: U(φ) = Σ kn[1 + cos(nφ − δn)], with multiplicities n = 1, 2, 3 (sometimes up to 6). An sp³–sp³ bond like ethane's has a 3-fold term; a double bond or amide a 2-fold term that keeps it planar. GROMACS also offers the equivalent Ryckaert–Bellemans form, a polynomial in cos ψ.

Butane. The trans state is lowest. Gauche is about 3 kJ/mol higher (≈ 1 kBT at 300 K), so both are populated: about 65% trans and 35% gauche. The barrier between them is about 12 kJ/mol (≈ 5 kBT), so changes are rare events, a few per 100 ps. Torsions with higher barriers (peptide bonds, about 80 kJ/mol) essentially never flip in a normal simulation.

Impropers. To keep planar groups flat (aromatic rings, sp² carbons, amides) or to stop a chiral centre inverting, force fields add an improper dihedral: a harmonic or cosine term on the angle defined by four atoms that are not in a chain.

1–4 interactions. Atoms 1 and 4 of a dihedral are only three bonds apart. Their Lennard-Jones and Coulomb interactions are kept, but scaled (OPLS: ×0.5 for both; AMBER: LJ ×0.5, Coulomb ×1/1.2), because the torsion term was fitted together with them. 1–2 and 1–3 pairs are excluded entirely.

Key idea. Dihedrals control conformations. The torsion term, the scaled 1–4 non-bonded forces, and impropers together decide which shapes a molecule adopts and how often it changes between them.

The algorithm

Computing φ and the forces:

b1 = r2 - r1;  b2 = r3 - r2;  b3 = r4 - r3
n1 = b1 × b2;  n2 = b2 × b3              # normals of the two planes
m1 = n1 × (b2 / |b2|)
phi = atan2(m1 · n2, n1 · n2)            # signed, -180..180, trans = 180

dU = dU/dphi                             # e.g. -Σ n k_n sin(n phi - delta_n)
F1 = -dU * |b2| / |n1|² * n1             # (sign depends on the phi convention)
F4 = +dU * |b2| / |n2|² * n2
F2, F3 from the projections of b1, b3 on b2, with F1 + F2 + F3 + F4 = 0
  1. Use atan2, not arccos: arccos loses the sign of φ and is inaccurate near 0° and 180°.
  2. Sign conventions differ between references. Always check analytic forces against a numerical derivative (the code tab does).
  3. A rotatable bond usually has several dihedrals (e.g. 9 H–C–C–H style terms in ethane); their sum is what matters.
  4. GROMACS: [ dihedrals ] function 9 (multiple cosines), 3 (Ryckaert–Bellemans), 4 or 2 (impropers); fudgeLJ, fudgeQQ set the 1–4 scaling. LAMMPS: dihedral_style charmm / opls / multi/harmonic, special_bonds for 1–4 scaling.

See it inside LAMMPS

The same idea in the real code, run live in your browser and checked by hand.