LAMMPS from scratch · 20 min

The conductor: the LAMMPS class

The whole program is one object of class LAMMPS. It owns fourteen members (atom, update, domain, force, neighbor …), and every line of an input script writes into one of them. You see which member each line of in.melt changes, then predict and check the bookkeeping that the members keep: step number, timestep, simulation time, ghost cutoff, and what clear throws away.

main() is three lines: create LAMMPS, read the script, delete LAMMPS
The fourteen members, what each one keeps, and the order they are created in
Which member each line of a script writes into
Predict step, dt, simulation time and the communication cutoff by hand, and what clear resets

1 · Understand: the whole program in three lines

When you type lmp -in in.melt, the function main() in src/main.cpp:69 does:

69  auto *lammps = new LAMMPS(argc, argv, lammps_comm);   // build the conductor
70  lammps->input->file();                                // read the script, line by line
71  delete lammps;                                        // clean up and finish

Everything else happens inside that one object. Its constructor builds the members in a fixed order (src/lammps.cpp:844), and the comments in the code say why the order matters:

852  comm = new CommBrick(this);       // Comm before Atom, so threads are known
855  neighbor = new Neighbor(this);
861  domain = new Domain(this);
865  atom = new Atom(this);
873  group = new Group(this);
875  force = new Force(this);          // after group, to create temperature
878  modify = new Modify(this);
880  output = new Output(this);        // after group, so "all" exists
882  update = new Update(this);        // after output, force, neighbor
883  timer = new Timer(this);

memory, error, universe and input are made just before. Every class in LAMMPS inherits from Pointers (src/pointers.h:110), which gives it a short name for each member. That is why code anywhere can simply write update->dt or atom->x.

Member	Class	Keeps
`memory`	Memory	allocates and frees all big arrays
`error`	Error	prints ERROR and WARNING messages and stops
`universe`	Universe	all processors, and partitions for multi-replica runs
`input`	Input	reads the script line by line and calls the right code
`atom`	Atom	every per-atom array: x, v, f, type, mass, charge …
`update`	Update	timestep, step number, units, the integrator (run)
`neighbor`	Neighbor	neighbour lists and the skin
`comm`	Comm	sends atoms between processors; ghost atoms
`domain`	Domain	the box, periodic boundaries, lattice, regions
`force`	Force	pair, bond, angle, dihedral, improper, kspace styles; unit constants
`modify`	Modify	the lists of fixes and computes
`group`	Group	named groups of atoms ("all" exists from the start)
`output`	Output	thermo, dumps and restart files
`timer`	Timer	stopwatches for the timing table (page A.9)

2 · Which member does each line change?

Tap a line of in.melt. The members it writes into light up.

memory

error

universe

input

atom

update

neighbor

comm

domain

force

modify

group

output

timer

3 · Modify: the bookkeeping, predicted by hand

The members keep running totals. Two of them are easy to get wrong. The simulation time is not step × dt: update adds each stretch of steps with the timestep that was in force, every time the timestep changes and at the end of every run (src/update.cpp:552):

552  void Update::update_time()
554    atime += (ntimestep - atimestep) * dt;   // time of the steps since the last update
555    atimestep = ntimestep;

The communication cutoff is how far beyond its own region a processor keeps copies (ghosts) of atoms: the pair cutoff plus the skin, because the neighbour list reaches that far (src/neighbor.cpp:379: cut = cutoff + delta, delta = skin).

first run then timestep second run skin pair cutoff

4 · Observe: clear rebuilds the conductor

clear is three lines of src/input.cpp:886: lmp->destroy(); lmp->create(); lmp->post_create();. Every member is thrown away and built fresh, as if LAMMPS had just started. (This page runs clear before every script for exactly that reason.) The new update sets the defaults of the chosen units (src/update.cpp:149):

units	default dt	default skin	code
`lj`	0.005	0.3	src/update.cpp:167
`real`	1	2	src/update.cpp:188
`metal`	0.001	2	src/update.cpp:209
`si`	1e-8	0.001	src/update.cpp:230
`cgs`	1e-8	0.1	src/update.cpp:251
`electron`	0.001	2	src/update.cpp:272
`micro`	2	0.1	src/update.cpp:293
`nano`	0.00045	0.1	src/update.cpp:314

after clear, units

What the conductor shows

An input script is a list of instructions to fourteen members. Knowing which member a command talks to tells you where its code is (page B.1) and what it can change.
Members remember: step count and simulation time carry over from one run to the next, until reset_timestep or clear.
Some commands check other members first. units asks domain whether a box exists and refuses if it does (src/input.cpp:2011): the box was measured in the old units.

Your turn

1. run 200 (dt 0.005), then timestep 0.002, run 500. What are step and time at the end?

2. After that, reset_timestep 0. What is the time now?

0. When the new step is smaller than the old one, Update::reset_timestep sets atime = 0 (src/update.cpp:508). Moving the step forward (e.g. reset_timestep 1000) would instead add the missing steps at the current dt.

3. Which member answers thermo_style custom vol, and which one answers temp?

vol is domain->xprd × yprd × zprd. temp is a compute (compute temp, made by output for thermo), kept in modify, which reads atom->v and atom->mass.