A tour of the source tree

Every command works the same way. The first word is looked up in src/input.cpp. If the command chooses a style (pair_style, fix, compute …), the style name is then looked up in a list of registered styles, and the matching class does the work. Its file name follows a rule:

lj/expand is LJ with the whole curve moved outwards by a distance Δ: E(r) = 4ε[(σ/(r−Δ))¹² − (σ/(r−Δ))⁶], with the cutoff also moved to r_c + Δ. It is used for particles of different sizes. Δ is the 5th word of pair_coeff. With Δ = 0 it must give exactly the lj/cut numbers.

Its file (section 3 rule: src/pair_lj_expand.cpp) is 90 % the same as pair_lj_cut.cpp. The maths inside the double loop over atoms i and their neighbours j:

106  if (rsq < cutsq[itype][jtype]) {          // r² below (rc + Δ)²
107    r = sqrt(rsq);
108    rshift = r - shift[itype][jtype];       // r − Δ   ← the only new idea
109    rshiftsq = rshift * rshift;
110    r2inv = 1.0 / rshiftsq;
111    r6inv = r2inv * r2inv * r2inv;          // 1/(r−Δ)⁶
112    forcelj = r6inv * (lj1 * r6inv - lj2);  // lj1 = 48εσ¹², lj2 = 24εσ⁶
113    fpair = factor_lj * forcelj / rshift / r;
       …
125    evdwl = r6inv * (lj3 * r6inv - lj4) - offset;   // lj3 = 4εσ¹², lj4 = 4εσ⁶

src/pair_lj_expand.cpp:106. Lines 112–113 are the force F = 24ε[2σ¹²/(r−Δ)¹³ − σ⁶/(r−Δ)⁷], divided by r so it can be multiplied by Δx, Δy, Δz. (The [itype][jtype] after lj1 … lj4 is dropped here for space.)