Styles: the plug-in design

1 · Understand: sockets and plugs

Think of a wall socket. Any lamp, kettle or charger works in it, because the socket fixes the shape and every device is made to fit. In LAMMPS the socket is a base class such as Pair; each pair style is a plug: a class built from Pair that fills in its functions. The main loop only ever talks to the socket.

The socket is a list of functions in src/pair.h:158. Those ending in = 0 have no default: every style must write them.

158  virtual void compute(int, int) = 0;                     // all forces, every step
164  virtual double single(int, int, int, int, double, ...)   // one pair (default: "not available")
187  virtual void settings(int, char **) = 0;                 // the pair_style arguments
188  virtual void coeff(int, char **) = 0;                    // the pair_coeff arguments
192  virtual double init_one(int, int) { return 0.0; }        // set up one type pair, return its cutoff

The plug says which socket it fits and which functions it fills, in src/pair_lj_cut.h:27:

27  class PairLJCut : public Pair {          // "built from Pair"
31    void compute(int, int) override;        // "my own version of compute"
32    void settings(int, char **) override;
33    void coeff(int, char **) override;
35    double init_one(int, int) override;
42    double single(int, int, int, int, double, double, double, double &) override;

2 · Registration: how "lj/cut" finds PairLJCut

1. The header has one line: PairStyle(lj/cut,PairLJCut); (src/pair_lj_cut.h:16).
2. At build time CMake collects these lines from every header of the core and the switched-on packages into one generated file, style_pair.h. A package that is off contributes nothing: that is why lj/smooth/linear was missing on page B.1.
3. When LAMMPS starts, src/force.cpp:92 turns each line into an entry of a dictionary: name → "a function that does new PairLJCut(lmp)" (src/force.cpp:41).
4. pair_style lj/cut 2.5 looks the name up (src/force.cpp:270). Found: the class is created and its settings() receives 2.5. Not found: the "Unrecognized pair style" error.

So adding a new style means: write one class, add one PairStyle line, rebuild. No other file of LAMMPS changes. Part E does exactly that.

A surprise: lj/expand with shift yes

Choose lj/expand, tick shift yes and run. The energies do not match. The hand side subtracts the energy at the cutoff, which is what shift yes is meant to do (the manual: the energy goes to 0 at the cutoff). For lj/expand the cutoff is r = r_c + Δ, where r − Δ = r_c, so the offset should be 4ε[(σ/r_c)¹² − (σ/r_c)⁶].

The code computes something else (src/pair_lj_expand.cpp:239):

239    double ratio = sigma[i][j] / (cut[i][j] + shift[i][j]);     // σ/(rc + Δ), not σ/rc
240    offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio, 12.0) - pow(ratio, 6.0));

Tick "use the offset exactly as written" and the hand values match LAMMPS to 15 digits. With ε = σ = 1, Δ = 0.2, r_c = 2.3, the shifted energy just inside the cutoff (r = 2.4999) is −0.0105 instead of 0, while lj/cut goes to 0 as it should. The forces are not affected, only the reported energy. The same line is in the current development version. This is what a hand calculation is for: it finds the one line where the code and the formula disagree.

What the sockets show

• The last row is always 0 when r reaches the cutoff: LAMMPS calls single() only if r² < r_c² (src/pair.cpp:1943). At exactly r_c the pair does not interact.
• shift yes subtracts one constant, computed once in init_one(), from every energy. The force column does not change at all: a constant has no slope.
• The F column is the force itself. Inside compute() LAMMPS uses F/r (called fpair), so it can multiply by Δx, Δy, Δz without a square root; pair_write multiplies back by r.
• Three different plugs, one socket: pair_write, compute()'s callers and the rest of LAMMPS never needed to know which one was in.

Your turn