LAMMPS from scratch · 25 min

Reading C++ with school maths

C++ looks foreign, but the parts of LAMMPS that do physics are ordinary maths: named numbers, tables, sums and cases. You learn the handful of ideas and symbols needed, then step line by line through the real LAMMPS force loop for three atoms, watching every number, and match the result to LAMMPS.

The six C++ ideas behind every physics routine, each matched to school maths
The symbols you will meet in LAMMPS and what they mean in words
Step through the real force loop of pair lj/cut for three atoms and match LAMMPS

1 · Understand: six ideas, all of them maths

C++	In words	In maths
`double r2inv;`	make a box called r2inv that holds a decimal number	let r2inv be a real number
`r2inv = 1.0 / rsq;`	put the value of 1/rsq into the box (a command, not an equation: `n = n + 1` is fine)	r2inv := 1/r²
`x[i][0]`	row i, column 0 of the table x (counting starts at 0)	x-coordinate of atom i
`for (j = 0; j < n; j++) { … }`	repeat the lines in { } for j = 0, 1, …, n−1	Σ over j
`if (rsq < cutsq) { … }`	do the lines in { } only when the test is true	a case of a piecewise function
`double init_one(int i, int j)`	a named recipe that takes i and j and returns a number	f(i, j)

A seventh idea, the class, is a form with blanks. LAMMPS has a form called Pair with a blank named compute; every force style fills in that blank its own way (section 4).

2 · The symbols you will meet

Symbol	Read it as	Example from LAMMPS
`;`	end of one instruction	`r6inv = r2inv * r2inv * r2inv;`
`{ }`	a group of lines that belong together	the body of a loop
`//`	a comment for humans; the compiler ignores the rest of the line	`// loop over neighbors of my atoms`
`*`	multiply (in declarations such as `double **x` it marks a table: ignore it)	`delx * fpair`
`+=`, `-=`	add to, subtract from	`f[i][0] += delx * fpair;`
`++`	add one	`ii++`
`[ ]`	pick an entry from a table	`lj1[itype][jtype]`
`->`	the part called … inside …	`atom->x`: the positions kept by the atom manager
`::`	belongs to	`PairLJCut::compute`: the compute of PairLJCut
`&`	here: keep only some bits of a number (a technical trick; skip it)	`j &= NEIGHMASK;`
`\|\|`, `&&`	or, and	`if (newton_pair \|\| j < nlocal)`
`#`	an instruction to the compiler, not to the computer	`#include "atom.h"`

3 · Modify, observe, match: step through the real loop

These are the lines of pair_lj_cut.cpp, lines 93–133 that do the work, for three atoms on a line (ε = σ = 1, cutoff 2.5). Choose the positions, then press Step and watch every box change.

x₁ x₂ x₃

What the loop is, in one line of maths

Line 119 adds each pair's push to atom i, and line 123 the opposite push to atom j. Over the whole loop that is

F_i = Σ_j F(r_ij) × (x_i − x_j) / r_ij with F(r) = 24(2r⁻¹³ − r⁻⁷) E = Σ_pairs 4(r⁻¹² − r⁻⁶)

For x = 0, 1.1, 2.3 the three pairs are at r = 1.1, 2.3 and 1.2:

pair	delx	fpair = F/r	evdwl
1–2	−1.1	1.4437231	−0.98337245
1–3	−2.3	−0.030232992	−0.026837948
2–3	−1.2	−1.8430778	−0.89096529

f₁ = (−1.1)(1.4437231) + (−2.3)(−0.030232992) = −1.5185595 f₂ = −(−1.1)(1.4437231) + (−1.2)(−1.8430778) = +3.7997887 f₃ = −(−2.3)(−0.030232992) − (−1.2)(−1.8430778) = −2.2812292 E = −0.98337245 − 0.026837948 − 0.89096529 = −1.9011757 check: f₁ + f₂ + f₃ = 0 (Newton's third law)

4 · The class: a form with a blank

In pair.h, line 158, the general form for every force style says:

virtual void compute(int, int) = 0;

Read it as: "every pair style must have a recipe called compute; this form leaves it blank" (= 0 means blank). In pair_lj_cut.h the Lennard-Jones style fills it in:

class PairLJCut : public Pair { PairLJCut is a filled-in Pair form public: void compute(int, int) override; "I fill in the compute blank" double init_one(int, int) override; "and the init_one blank" …

When your input says pair_style lj/cut, LAMMPS picks this filled-in form, and every timestep it calls its compute: the loop you just stepped through. Part B explains how the name lj/cut is connected to PairLJCut.

Your turn

1. Translate into maths: forcelj = r6inv * (lj1 * r6inv - lj2); with lj1 = 48 and lj2 = 24.

2. Move atom 3 to x = 3.0 in the stepper. Which pairs are skipped at line 113, and what is f₃?

3. Why does LAMMPS write f[j][0] -= delx * fpair instead of looping over the pair a second time from atom j's side?