Lab: a molecule by hand

The setup

Our molecule is n-butane, CH₃–CH₂–CH₂–CH₃, in the united-atom picture: each carbon together with its hydrogens is one bead. Four beads, three bonds, two angles, one dihedral. Parameters are those of the TraPPE force field. From here on we use the units of GROMACS, so the numbers can be compared with real input files.

(TraPPE keeps the bonds rigid; for the flexible runs we give them a typical C–C spring.) Every energy below is compared with k_B T = 2.494 kJ/mol: anything much larger is almost never visited at 300 K, anything smaller is visited all the time.

Step 1

Why a molecule needs bonded terms, and exclusions

Why: If LJ were enough, we could just put four beads together and let them stick. Let us check what LJ says about two carbons at the length of a real C–C bond.

What the numbers say: LJ describes how separate molecules touch, at about 0.4 nm. A chemical bond is a much shorter, much stronger thing: at 0.154 nm LJ gives 123 607 kJ/mol, 50 000 times k_B T, and would blow the molecule apart in one step. Two consequences. First, the shape of a molecule is held by its own terms: bonds, angles, dihedrals. Second, atoms 1 bond apart (1–2) and 2 bonds apart (1–3) must be excluded from LJ, otherwise 230 kJ/mol of fake repulsion fights the angle term. Every MD code keeps such an exclusion list.

Step 2

The bond: a stiff spring that sets the time step

Why: The bond keeps two beads near 0.154 nm. The simplest description is a spring. We need its force for the dynamics, and its vibration period, because the time step has to follow the fastest motion in the system (Lab 1).

What the numbers say: At 300 K the bond only wiggles by 2 %, so for most purposes it is effectively rigid, yet it is by far the fastest motion: one period every 35 fs. With about 15–20 steps per period (Lab 1), this alone limits us to about 2 fs. In an all-atom model the C–H bond vibrates every 11 fs, forcing 0.5–1 fs. We spend most of the computer time following a motion we do not care about. Step 7 shows the fix.

Step 3

The angle: dot product, energy, force

Why: Bonds alone would let the chain fold flat onto itself. The angle term keeps C–C–C near 114°. To use it we must get the angle from coordinates, then turn the energy into forces on the atoms.

What the numbers say: Turning an angle into forces needs nothing beyond the derivative and a lever arm: the torque dU/dθ divided by the bond length. At 300 K the angle wanders by about 4°, much more than the bond (2 %), and it vibrates 2.5 times more slowly (91 fs vs 37 fs). That ratio is exactly what we gain in step 7 when the bonds are frozen and the angle becomes the fastest motion.

Step 4

The dihedral: where the chemistry is

Why: With bonds and angles fixed, butane can still twist around its central bond. This twist decides whether the chain is straight (trans) or kinked (gauche), and in polymers and lipids it controls the whole shape. The dihedral angle φ needs four atoms and two cross products.

How fast does the molecule switch? An escape over a barrier happens at a rate of (vibration frequency) × e^(−barrier/k_B T):

Trans fraction in a Langevin run of 200 molecules at 300 K; one jump every 23 ps per molecule.

What the numbers say: The dihedral needs only two cross products, and its three-term formula gives the whole story of butane: trans is lowest, gauche is 3.6 kJ/mol higher (1.4 k_B T, so still common: about one third of the molecules), cis is almost never visited. The full sum (0.663) differs a little from the three-state count (0.677) because the wells have different widths. Limitation, by numbers: the molecule switches only every few tens of picoseconds. Starting from all-trans or all-gauche, 200 molecules need about 100 ps before both runs agree on 0.65–0.68. A longer chain, or a higher barrier, makes this nanoseconds to microseconds: sampling dihedrals is often the real bottleneck of a simulation (Chapter 7, metadynamics).

Step 5

1–4 pairs: exclude, scale or keep?

Why: Atoms 1 and 4 are three bonds apart. They are not held at a fixed distance by bonds and angles: their distance changes with φ. Should they feel LJ? Different force fields answer differently, and the numbers show this is not a detail.

What the numbers say: In gauche the two CH₃ ends are only 0.313 nm apart, inside σ = 0.375, so their LJ is strongly repulsive (+19 kJ/mol). TraPPE fitted its dihedral with this pair excluded; OPLS fitted its own dihedral with this pair at half strength. Mixing the two, TraPPE's dihedral with OPLS's rule, turns butane 96 % trans, almost no gauche left: the chains come out far too straight. The lesson: a force field is a package. Dihedral parameters, 1–4 scaling and exclusions were fitted together, and must be used together. In GROMACS this is gen-pairs, fudgeLJ and fudgeQQ; in LAMMPS special_bonds.

Step 6

Combining rules: CH₃ meets CH₂

Why: We have LJ parameters for CH₃–CH₃ and CH₂–CH₂, but a simulation also needs CH₃–CH₂, and with many atom types there are far too many pairs to fit each one. Combining rules make the cross parameters from the pure ones.

What the numbers say: For similar atoms like CH₃ and CH₂, the two rules agree to 0.0001 nm and the choice does not matter. For very different sizes, such as a small cation next to a large anion, they disagree by 12 % in σ: at 0.25 nm one rule says strong repulsion, the other mild attraction, and the contact distance moves by 0.034 nm. That directly shifts the radial distribution function and the binding. Like 1–4 scaling, the combining rule belongs to the force field (GROMACS comb-rule, LAMMPS pair_modify mix). Ion–ion and ion–solvent pairs are often fitted separately for exactly this reason.

Step 7

Constraints: SHAKE by hand, and the time step it buys

Why: The bond is stiff (2 % wiggle) and fast (35 fs). If we freeze it at 0.154 nm, the fastest motion becomes the angle (91 fs) and the time step can grow. SHAKE does this: after each normal step, it moves atoms along their old bond directions until every bond length is right again.

Energy fluctuation in NVE (kJ/mol, 20 molecules × 10 ps): smaller means a more accurate trajectory.

What the numbers say: The error falls very fast (1.5×10⁻³ → 7.5×10⁻⁶ → 1.8×10⁻¹⁰ nm), so two or three iterations are enough for one bond; coupled bonds in a chain need about 10, because fixing one bond disturbs its neighbours. The payoff: rigid bonds at 6 fs are as accurate as flexible bonds at 2 fs (0.039 vs 0.048), three times fewer steps for the same physics, and flexible bonds blow up at 10 fs while rigid ones still run at 12 fs. This matches the period ratio of step 3 (91/37 ≈ 2.5). In all-atom models, constraining only the bonds to H already doubles the time step (1 -> 2 fs). Limitations: the iterations cost time and are hard to parallelise (GROMACS prefers LINCS, and SETTLE for water); constraining angles in long chains converges badly; and every constraint removes a degree of freedom, which must be subtracted when computing the temperature.

Summary