LAMMPS from scratch · 20 min

How atoms are stored

Inside LAMMPS an atom is not an object: it is one row number in a set of long arrays (x, v, f, type, tag, image …). You see which arrays an atom_style creates, why every processor also stores copies of other atoms (ghosts) and count them by hand, watch LAMMPS reorder its atoms in memory, and check by hand how image flags rebuild the true path of an atom that crossed the periodic box.

Per-atom arrays, the local index i and the atom ID (tag)
Owned and ghost atoms: count the ghosts of the melt crystal by hand
Why LAMMPS sorts atoms in memory, and what that does to the order you see
Image flags: unwrapped x = x + ix·L, checked against LAMMPS

1 · Understand: an atom is a row number

The class Atom (src/atom.h:72) holds one array per property. Atom number i is simply row i of every array:

46  int nlocal, nghost;      // atoms this processor owns, and copies of other atoms
47  int nmax;                // length the arrays are allocated for
72  tagint *tag;             // the atom ID you see in data files and dumps
73  int *type, *mask;        // type 1, 2 …; mask = which groups it belongs to (one bit per group)
74  imageint *image;         // how many times it crossed the periodic box (section 5)
75  double **x, **v, **f;    // position, velocity, force: x[i][0], x[i][1], x[i][2]

Rows 0 … nlocal−1 are the atoms this processor owns and moves. Rows nlocal … nlocal+nghost−1 are ghosts: copies of atoms that sit just outside its region, kept so every owned atom sees all its neighbours (section 3). The row number i is not the atom ID: tag[i] is. They start equal and drift apart (section 4).

Which arrays exist is decided by atom_style. Every style has the defaults (src/atom_vec.cpp:33); each adds its own list (fields_grow):

atom_style	arrays per atom	bytes per atom (approx.)
`atomic`	tag, type, mask, image, x, v, f	4 × 4 + 9 × 8 = 88
`charge`	+ q (charge)	88 + 8 = 96
`full`	+ q, molecule, and lists of the bonds, angles, dihedrals, impropers and special neighbours of each atom	a few hundred, depending on the per-atom limits

src/atom_vec_charge.cpp:33 · src/MOLECULE/atom_vec_full.cpp:35. Integers are 4 bytes, doubles 8, in the default build.

2 · Predict: how many ghost atoms?

In serial, the box is the processor's region, and its neighbours across a periodic face are copies of its own atoms. LAMMPS makes the copies one direction at a time (src/comm_brick.cpp:775). For each face it copies every atom within the communication cutoff of the opposite face (page B.2: pair cutoff + skin):

429  slabhi[iswap] = sublo[dim] + cutghost[dim];   // atoms near the low face → ghosts beyond the high face
456  slablo[iswap] = subhi[dim] - cutghost[dim];   // atoms near the high face → ghosts beyond the low face
842  if (x[i][dim] >= lo && x[i][dim] <= hi)       // inside the slab: send it

Because x, then y, then z are done in turn, ghosts of ghosts fill the edges and corners: the result is the whole box grown by the cutoff on every side. At step 0 the melt atoms sit on FCC planes every a/2 = 0.8397981 apart (page A.5), 20 planes per direction, so the count can be done by hand.

pair cutoff skin

3 · Observe: the count moves in steps

Try skin 0.3, 0.5 and 0.9 with cutoff 2.5. The first two give the same 5841 ghosts; the third jumps. Ghosts come in whole planes, so nothing changes until the cutoff passes the next plane at a multiple of 0.8397981 (3 × 0.8397981 = 2.5193943, 4 × 0.8397981 = 3.3591924). After melting there are no planes, and the count changes smoothly: 5506 at the end of melt (page A.9).

Ghosts cost memory and time: here there are more ghosts (5841) than owned atoms (4000). A small box, or many processors each holding a small region, makes this worse (page B.7).

4 · Observe: LAMMPS reorders its atoms

Every 1000 steps by default (src/atom.cpp:103), Atom::sort() (src/atom.cpp:2379) reorders all arrays so atoms close in space sit close in memory, which makes the force loop faster. A dump without dump_modify sort id writes atoms in memory order, so the IDs come out shuffled. Here sorting is set to every 100 steps:

5 · Match: image flags by hand

When an atom leaves the box through a face, LAMMPS puts it back on the other side, so x always stays between 0 and L = 16.795962, and it adds or subtracts 1 from the atom's image flag ix. The true, unwrapped position is then (src/dump_custom.cpp:2665):

2665  buf[n] = x[j][0] + ((image[j] & IMGMASK) - IMGMAX) * xprd;    // xu = x + ix · L

Diffusion (the mean squared displacement lesson) must use the unwrapped xu: an atom that steps across a face would otherwise appear to jump by a whole box length. Below, the atoms in the step-200 dump with ix ≠ 0: predict xu from x and ix.

Run section 4 first.

Your turn

1. Why are there 4 ghost planes beyond the high face but only 3 beyond the low face, with cutoff 2.8?

2. About how much memory do the per-atom arrays of melt need at step 0?

3. An atom has x = 0.4 and ix = −2 in a box of length 20. Where has it really gone?