Periodic boundary conditions

The idea

Even a million atoms is a very small piece of matter. In a cube of 1000 atoms, about half of them touch a wall or surface. Walls change everything: density, structure and dynamics would all be wrong for a bulk liquid.

Periodic boundary conditions (PBC) remove the walls. Imagine the box copied in every direction to fill all space. Every copy moves exactly like the original, so we only simulate one box, but its atoms feel neighbours on all sides.

• Wrapping: an atom that leaves through the right face re-enters through the left face. The number of atoms stays constant.
• Minimum image: when computing forces, each atom interacts with the nearest copy of every other atom. Two atoms at opposite edges are really neighbours.
• Cutoff rule: the cutoff must be at most L/2. Otherwise an atom could interact with two copies of the same neighbour, or with its own image.
• Unwrapped coordinates: for diffusion (Lesson 6.2) you need the true distance travelled. Store image counters, or the "unwrapped" positions, alongside the wrapped ones.

PBC is not perfect. The system is still artificially periodic, so very long-range effects such as large-scale fluctuations, slow hydrodynamic flows and some electrostatics feel the box size. Always check important results with a bigger box.

The algorithm

Two one-line functions do almost all the work in a cubic box of side L:

def wrap(x):                 # keep positions in [0, L)
    return x - L * floor(x / L)

def min_image(d):            # shortest vector to the nearest copy
    return d - L * round(d / L)

1. After each move, wrap the positions. Optionally count crossings in an integer "image" array, so the unwrapped position is wrap(x) + image·L.
2. In the force loop, compute d = r_i − r_j and apply min_image before taking the length.
3. Check that r_c ≤ L/2 before running. LAMMPS and GROMACS stop with an error if it is not.
4. For analysis, use unwrapped coordinates for diffusion, and minimum-image distances for structure (RDF).

Non-cubic (triclinic) boxes need a slightly more general version, using the box matrix. The idea is identical.