Bonds and angles

The idea

So far our atoms only interacted through Lennard-Jones forces. In a molecule, some atoms are chemically bonded. A force field describes this with simple terms that depend only on the internal geometry (Lesson 4.3 puts them all together):

• Bonds (two atoms): the energy rises when the distance r moves away from its equilibrium value r₀. Near r₀ a parabola, U = ½k_b(r − r₀)², is accurate. Real bonds can break; the Morse potential captures that, but most force fields keep the harmonic form because bonds do not break at normal temperatures.
• Angles (three atoms, i–j–k): the energy rises when the angle θ at the central atom j moves from θ₀: U = ½k_θ(θ − θ₀)².

Forces. The force is minus the gradient of the energy with respect to each atom's position. For a bond the force acts along the bond and pulls the atoms back towards r₀. For an angle the forces on the outer atoms act perpendicular to their bonds, opening or closing the angle; the central atom gets minus their sum. In every case the forces add up to zero (Newton's third law). Always test analytic forces against a numerical derivative of the energy; the code tab does this.

Fast vibrations. Bonds are stiff and hydrogen is light, so an O–H or C–H bond vibrates with a period of only about 10 fs. With 10–20 steps per vibration needed (Lesson 1.4), flexible bonds to hydrogen force Δt ≈ 0.5–1 fs. This is why most simulations constrain these bonds (Lesson 4.4) and use 2 fs.

The algorithm

Forces from a harmonic bond (i–j) and angle (i–j–k):

# bond
d = r_i - r_j;  r = |d|
F_i = -k_b (r - r0) d / r;   F_j = -F_i

# angle at j
a = r_i - r_j;  b = r_k - r_j
cos_t = a·b / (|a||b|);  theta = arccos(cos_t)
g = k_theta (theta - theta0) / sin(theta)
F_i = g * ( b/(|a||b|) - cos_t a/|a|² )
F_k = g * ( a/(|a||b|) - cos_t b/|b|² )
F_j = -F_i - F_k

1. Bonded atoms (1–2) and atoms sharing an angle (1–3) are usually excluded from the Lennard-Jones and Coulomb interactions. Their interaction is already in the bonded terms.
2. Watch the convention! GROMACS writes ½k(r − r₀)². AMBER and LAMMPS write K(r − r₀)², so K = k/2. Units also differ: kJ mol⁻¹ nm⁻² vs kcal mol⁻¹ Å⁻².
3. θ near 180° makes sin θ → 0 in the formula. Linear molecules need special treatment (or a cosine-based angle potential).
4. Check: finite-difference forces must match the analytic ones to many digits, and all forces must sum to zero.