LAMMPS from scratch · 25 min

The families of styles

Pair, bond, angle, dihedral, improper, kspace, fix and compute styles each have their own moment in a timestep and their own number in the output. One united-atom butane molecule, built from bond lengths, angles and a dihedral you choose, lets every family speak once: you predict each energy and each geometric value by hand and LAMMPS reports the same. On the way, a hand calculation shows that fix restrain dihedral holds the molecule 180° away from the angle the manual promises.

The families, the base class of each, and where in a timestep each one is called
Forces come from pair, bonded and kspace styles; fixes change the motion; computes only measure
Build a molecule from internal coordinates and predict E_bond, E_angle, E_dihed, E_vdwl and the fix energy by hand
Special bonds: why the 1–4 LJ weight matters

1 · Understand: who speaks when in one timestep

Page B.3 showed one family (pair) in detail. Every family works the same way: a base class with sockets, many plugs. What differs is when the main loop calls them. The order inside one step of Verlet::run() (src/verlet.cpp:258):

258  modify->initial_integrate(vflag);        // FIX: first half of velocity Verlet (nve, nvt …)
306  force_clear();                           // all forces f = 0
316  force->pair->compute(eflag,vflag);       // PAIR
321  if (force->bond) force->bond->compute(eflag,vflag);          // BOND
322  if (force->angle) force->angle->compute(eflag,vflag);        // ANGLE
323  if (force->dihedral) force->dihedral->compute(eflag,vflag);  // DIHEDRAL
324  if (force->improper) force->improper->compute(eflag,vflag);  // IMPROPER
329  force->kspace->compute(eflag,vflag);     // KSPACE (long-range Coulomb)
347  if (n_post_force_any) modify->post_force(vflag);   // FIX: add or change forces (restrain, addforce, langevin …)
348  modify->final_integrate();               // FIX: second half of velocity Verlet
349  if (n_end_of_step) modify->end_of_step();          // FIX: averages, printing, thermostats …
356  output->write(ntimestep);                // asks COMPUTEs for the numbers it prints

Family	Base class	Its job	Thermo keyword	In this browser
pair	Pair	forces between nearby atoms	`evdwl`, `ecoul`	lj/cut, morse, buck …
bond	Bond	force along each listed bond	`ebond`	harmonic, fene, morse …
angle	Angle	force on each listed angle	`eangle`	harmonic, cosine, charmm …
dihedral	Dihedral	twist around the middle bond of 4 atoms	`edihed`	opls, harmonic, charmm …
improper	Improper	keeps 4 atoms flat or a centre chiral	`eimp`	harmonic, cvff, umbrella
kspace	KSpace	the long-range part of Coulomb	`elong`	none here (page D.6)
fix	Fix	anything that changes the motion: integrate, thermostat, add forces	`f_ID`	79 styles
compute	Compute	measures, never changes anything	`c_ID`	68 styles

A useful rule: forces come from pair, bonded and kspace styles; fixes move atoms and may add forces; computes only read.

2 · Modify: build a butane molecule

Four united atoms, CH₃–CH₂–CH₂–CH₃, with TraPPE-like parameters. The page places the atoms from three bond lengths (1.55, 1.53, 1.56 Å), two angles (110°, 116°) and the dihedral φ you choose, writes the data file (page A.8), and gives it to LAMMPS.

dihedral φ (°) 1–4 LJ weight add fix restrain restrain K restrain φ₀ (°)

butane.data and the script

3 · Predict by hand, then run

4 · Observe: what each family said

bond, angle, dihedral: each one sums over its own list from the data file (3 bonds, 2 angles, 1 dihedral). LAMMPS's harmonic forms have no ½ (page A.8). The OPLS dihedral is ½K₁(1 + cos φ) + ½K₂(1 − cos 2φ) + ½K₃(1 + cos 3φ) + ½K₄(1 − cos 4φ) (src/MOLECULE/dihedral_opls.cpp:177).
pair: LJ is computed for every pair within the cutoff, then special_bonds multiplies 1–2, 1–3 and 1–4 pairs by its three weights. In butane only atoms 1 and 4 are a 1–4 pair, and with weight 0 (TraPPE) E_vdwl is exactly 0. At φ = 65° with weight 1 they add over 4 kcal/mol, more than all bonded terms together: the dihedral parameters were fitted assuming that weight, so the weight is part of the force field.
fix: restrain adds a force in post_force and, with fix_modify energy yes, its energy to PotEng.
compute: bond/local, angle/local and dihedral/local measured the geometry and changed nothing. Their numbers equal the hand geometry.

5 · A surprise: where does fix restrain hold the dihedral?

The manual for fix restrain says E = K[1 + cos(nφ − d)] with d = φ₀ + 180°, so that the energy is lowest at φ = φ₀. The code (src/fix_restrain.cpp:135 and src/fix_restrain.cpp:642) uses d = φ₀:

135  target[nrestrain] = DEG2RAD * utils::numeric(FLERR,arg[iarg+7],false,lmp);   // φ₀, no + 180
642  p = p*cos_target[m] + df1*sin_target[m];      // cos(φ − φ₀)
645  p += 1.0;
647  edihed += k * p;                              // E = K[1 + cos(φ − φ₀)]: HIGHEST at φ = φ₀

Section 3 already showed it: the hand value with the code's formula matches f_r; the manual's formula does not. To see what it does to a molecule, minimize with the restraint as the only dihedral term and look where φ ends:

The same lines are in the current development version and in versions back to at least 2018. If you use fix restrain dihedral (for example to scan a torsion profile), check φ after the run, or give φ₀ − 180°.

Your turn

1. Which φ gives E_dihed = 0? Predict from the OPLS formula, then set it.

2. Set φ = 0 (cis) with 1–4 weight 1. Why does E_vdwl explode?

3. Which family would you use to keep atom 1 from moving, and why not a compute?

A fix: fix 2 one setforce 0 0 0 (with v = 0) or fix spring/self. A compute can only measure; it is never called in a way that changes x, v or f.