Lab 2.7 · 30 min

Lab: building a system by hand

Before the first time step, a simulation needs units, starting positions, starting velocities and a periodic box, plus rules for which pairs to compute. We build a 9-atom system by hand and let the numbers show why each of these choices is made, and what goes wrong without it.

Convert argon to reduced units, and place 9 atoms on a lattice instead of at random
Prepare velocities: remove the drift, count degrees of freedom, scale to the target temperature
Use periodic boundaries, the minimum image, a cutoff and a neighbour list, and see what each one fixes

The setup

To keep every number checkable, we work in two dimensions with 9 argon atoms in a square box of side L = 3.6 (reduced units, Lesson 2.1). The steps are exactly those of a real 3D simulation; only the counting is smaller.

atoms N = 9 box L = 3.6 × 3.6, periodic in x and y target T = 0.72 (86 K for argon) time step dt = 0.005 cutoff r_c = 1.8 (the largest this box allows, see step 6)

Step 1

Units: why we do not use joules and metres

Why: Argon quantities in SI units are tiny numbers like 10⁻²¹ J. Writing them out at every step invites mistakes, and the same simulation would have to be repeated for every gas. Reduced units fix both problems.

argon: σ = 3.405×10⁻¹⁰ m ε/k_B = 119.8 K mass = 39.948 g/mol energy unit ε = 119.8 × 1.3806×10⁻²³ J/K = 1.6540×10⁻²¹ J mass unit m = 39.948 × 1.6605×10⁻²⁷ kg = 6.6335×10⁻²⁶ kg time unit τ = σ × √(m/ε) = 3.405×10⁻¹⁰ × √(4.0106×10⁻⁵) = 2.156×10⁻¹² s = 2.156 ps T* = 0.72 -> 0.72 × 119.8 = 86.3 K dt* = 0.005 -> 0.005 × 2.156 ps = 10.8 fs

What the numbers say: In reduced units every number is close to 1, and one simulation describes all Lennard-Jones gases: to get krypton, convert back with krypton's σ, ε and mass. Limitation: you must convert back at the end, and a mixture needs one species chosen as the reference.

Step 2

Starting positions: why a lattice and not random

Why: Atoms must start somewhere. Dropping them at random looks natural, but two atoms almost always land on top of each other, and the LJ energy explodes at short distance.

Random (9 positions drawn by the script, rounded): atom 1 at (3.40, 1.84) and atom 6 at (3.25, 1.72):

dx = −0.15, dy = −0.12 r = √(0.0225 + 0.0144) = √0.0369 = 0.192 U = 4 × (0.192⁻¹² − 0.192⁻⁶) = 1.6 × 10⁹ F = 9.9 × 10¹⁰ velocity change in ONE step = F × dt = 4.9 × 10⁸ (a typical thermal speed is about 1)

Lattice: a 3 × 3 square grid with spacing a = L/3 = 1.2, positions x = (i + ½)·a, y = (j + ½)·a:

(0.6,0.6) (1.8,0.6) (3.0,0.6) (0.6,1.8) (1.8,1.8) (3.0,1.8) (0.6,3.0) (1.8,3.0) (3.0,3.0) every atom has 4 neighbours at 1.2 and 4 at 1.2×√2 = 1.697 U(1.2) = −0.8910 U(1.697): 1.697² = 2.88, 2.88³ = 23.89, r⁻⁶ = 0.04186, r⁻¹² = 0.001752 U = 4 × (0.001752 − 0.04186) = −0.1604 energy per atom = ½ × (4 × −0.8910 + 4 × −0.1604) = −2.1028 (the ½: each pair energy is shared by two atoms) total = 9 × −2.1028 = −18.925

What the numbers say: The random start puts 1.6 billion ε into one pair; the first step would fire the atoms apart at absurd speeds. The lattice has no overlaps and every atom has the same environment. Limitations: the start is an ordered crystal, so a liquid must first be melted and equilibrated; N must fit the lattice (n² in 2D, 4n³ for fcc in 3D); molecules need packing tools such as Packmol.

Step 3

Starting velocities: drift, degrees of freedom, scaling

Why: Velocities set the temperature. We draw them from a bell-shaped (Gaussian) distribution, as in a real gas (Maxwell–Boltzmann), but a few random numbers never average exactly to zero, and temperature must be counted correctly.

9 random (vx, vy), two decimals: (−0.33,−0.88) (−0.66,−0.67) (+0.38,−0.11) (+1.48,−1.83) (−0.00,−0.89) (+0.78,−2.12) (−0.34,+0.21) (−1.48,+0.99) (+0.18,+1.01) 1. drift (average velocity): vx: (−0.33 −0.66 +0.38 +1.48 −0.00 +0.78 −0.34 −1.48 +0.18)/9 = 0.01/9 = +0.0011 vy: (−0.88 −0.67 −0.11 −1.83 −0.89 −2.12 +0.21 +0.99 +1.01)/9 = −4.29/9 = −0.4767 kinetic energy of the drift = ½ × 9 × 0.4767² = 1.02 (of 8.87 in total) -> subtract the drift from every atom 2. temperature, 2D: T = Σ m v² / N_f N_f = 2N − 2 = 16 (2 directions per atom, minus 2 for the fixed total momentum) Σ v² after removing the drift = 15.696 T = 15.696 / 16 = 0.9810 3. scale every velocity by λ = √(0.72 / 0.9810) = 0.8567 -> T = 0.7200 wrong count N_f = 2N = 18: T would read 0.64, 11 % too low

What the numbers say: Without step 1 the whole box would slide downwards, and 11 % of the "thermal" energy would really be motion of the box. With only 9 atoms, using 2N instead of 2N − 2 gives an 11 % error in T; for 10 000 atoms it is 0.01 %. Limitations: 9 random numbers make a very rough bell curve, but any shape relaxes to Maxwell–Boltzmann within a few collisions; what must be right is the total momentum (zero) and T.

Step 4

Periodic boundaries: why the box has no walls

Why: We want to simulate a piece of a bulk liquid, not a tiny droplet. In a small group of atoms almost all of them sit on the surface, where their environment is completely different.

system (no periodic boundaries)	atoms on the edge	fraction
3 × 3 (our 2D box)	8 of 9	89 %
10 × 10	36 of 100	36 %
100 × 100	396 of 10 000	4 %
10 × 10 × 10 (3D)	488 of 1000	49 %

With periodic boundaries the box is surrounded by copies of itself. An atom that leaves on one side comes back on the other:

x_wrapped = x − L × floor(x / L) x = 3.75 -> 3.75 − 3.6 × 1 = 0.15 x = −0.30 -> −0.30 − 3.6 × (−1) = 3.30

What the numbers say: Even 1000 atoms in 3D would be half surface without periodic boundaries. With them, every atom is surrounded by atoms and there is no surface at all. Limitations: the system is artificially periodic, so the box must be larger than any structure or correlation you study, and some properties need size corrections (for example diffusion, Lesson 6.2).

Step 5

The minimum image: neighbours through the wall

Why: With periodic copies, every atom has infinitely many images. For short-range forces we only want the nearest copy of each other atom. Measuring distances inside the box alone misses neighbours that sit just across a wall.

atom 1 at (0.6, 0.6), atom 3 at (3.0, 0.6) raw dx = 3.0 − 0.6 = 2.4 (looks far away: beyond the cutoff 1.8) minimum image: dx − L × round(dx / L) = 2.4 − 3.6 × 1 = −1.2 -> atom 3 is a NEAREST neighbour of atom 1, just across the left wall

	neighbours of atom 1	energy of atom 1	total energy
with minimum image	4 at 1.2, 4 at 1.697	−2.1028	−18.925
without (raw distances)	2 at 1.2, 1 at 1.697	−0.9712	−11.975

What the numbers say: Forgetting the minimum image throws away 54 % of the energy of a corner atom and 37 % of the total: atoms near the walls would feel almost no neighbours and drift to the middle. The rule needs one line of code and is applied to every pair, every step.

Step 6

Cutoffs: why we ignore far pairs, and the price

Why: The number of pairs grows as N², but the LJ energy of a far pair is tiny. A cutoff r_c skips pairs farther than r_c. It saves time, but it changes the energy, and the minimum image only works if r_c ≤ L/2.

cutoff r_c	energy per atom (infinite lattice)	error
1.8 (our box)	−2.1028	4.8 %
2.5	−2.1445	2.9 %
4.0	−2.2015	0.3 %
∞	−2.2089	0

jump at the cutoff: U(1.8) = 4 × (1.8⁻¹² − 1.8⁻⁶) = 4 × (0.00087 − 0.02940) = −0.1141 -> when a pair crosses r = 1.8 the energy jumps by 0.114: energy is not conserved fix: shifted potential U(r) − U(r_c), which is exactly 0 at the cutoff box size rule: r_c ≤ L/2, otherwise an atom could interact with two images of the same atom our box: L/2 = 1.8, so r_c = 2.5 is impossible with 9 atoms r_c = 2.5 needs L ≥ 5.0: 2D at this density N ≥ (5.0/1.2)² ≈ 17 3D liquid (ρ = 0.84): N ≥ 0.84 × 5³ ≈ 105

What the numbers say: A cutoff of 2.5 already loses about 3 % of the energy in this dense lattice; codes add a "tail correction" for the missing part. The cutoff also creates a small jump in energy unless the potential is shifted. And it is why real simulations need at least ~100 atoms: the box must be twice the cutoff.

Step 7

Neighbour lists: computing only the pairs that matter

Why: Even with a cutoff, a code that checks all N(N−1)/2 pairs every step to find the close ones wastes most of its time. A neighbour list stores the pairs within r_c + skin, and is rebuilt only when atoms have moved far enough (Lesson 2.6).

N	all pairs N(N−1)/2
9	36
1 000	499 500
100 000	4 999 950 000

3D liquid, ρ = 0.84, r_c = 2.5, skin = 0.3, list radius 2.8: neighbours per atom = (4/3) π × 2.8³ × 0.84 = 77.2 N = 1000: 1000 × 77.2 / 2 = 38 620 pairs (13× fewer than 499 500) when to rebuild: when some atom has moved more than skin/2 = 0.15 speed ≈ 1.5, dt = 0.005 -> 0.15 / (1.5 × 0.005) ≈ every 20 steps skin 0.6 instead: 104.8 neighbours (+36 % pairs), rebuild every ~40 steps

What the numbers say: For our 9 atoms the list holds all 36 pairs, so it saves nothing: lists pay off only for large N, where they turn an N² cost into a cost proportional to N. The skin is a trade-off: larger skin, fewer rebuilds but more pairs to check each step.

Summary

choice	what it fixes	price / limitation
reduced units	awkward 10⁻²¹ numbers; one run for all LJ gases	convert back at the end
lattice start	overlaps (U = 1.6×10⁹ for one random pair)	must melt and equilibrate
remove drift, N_f = 2N−2	sliding box; 11 % error in T for N = 9	none
periodic boundaries	89 % surface atoms in a 3×3 group	artificial periodicity, size effects
minimum image	54 % of a corner atom's energy missing	requires r_c ≤ L/2
cutoff (shifted)	N² far pairs	~3 % energy at r_c = 2.5 (tail correction)
neighbour list	499 500 pair checks for 1000 atoms	rebuilds; memory

See it inside LAMMPS

The same idea in the real code, run live in your browser and checked by hand.

A.5 · Building the system