Lab: finite-size effects by hand

The setup

Two kinds of finite-size effects appear:

Systems: an ideal gas (exact answers), the 2D LJ liquid of earlier labs in boxes of 9, 25 and 100 atoms, and a 3D LJ liquid (ρ = 0.8, T = 1.0) with 32 to 500 atoms. For real liquids we use water (η = 0.85 mPa s) and an ionic liquid (η = 50 mPa s).

Step 1

The g(r) tail is 1 − 1/N, not 1

Why: g(r) compares the number of pairs at distance r with an ideal gas of the same density. With a fixed number of particles, each one sees only N − 1 others, while the usual normalization uses the density N/V. That leaves a small offset at large r, even for an ideal gas, where the exact g is 1.

What the numbers say: The ideal gas confirms the formula exactly: the tail sits at 1 − 1/N everywhere. In a dense liquid the offset is smaller, about S(0)/N, where S(0) = ρk_BTκ_T is small for a stiff liquid. The first peak and minimum are what usually matter, and they hardly change from 25 to 100 atoms (2.92 vs 2.95). The 9-atom box is too small: its half-width 1.8 cuts right through the first minimum, so its structure is distorted (peak 2.62, minimum 0.985). Rule of thumb: the box must be several times larger than the structure you want to see.

Step 2

Integrals of g(r): Kirkwood–Buff and the Ganguly–van der Vegt correction

Why: Kirkwood–Buff integrals G = ∫(g − 1) dV over all space connect g(r) to compressibilities and partial molar volumes, and running coordination numbers are integrals too. A tiny offset in g(r), multiplied by a large volume, becomes a large error.

What the numbers say: The running integral drifts away linearly with the area (or volume) instead of levelling off, and the drift is pure artefact. The Ganguly–van der Vegt factor accounts for the particles that are "missing" because they are already counted within R; for the ideal gas it restores g = 1 exactly. In 3D with 500 molecules, the error in a Kirkwood–Buff integral is about 9 cm³/mol, comparable to the integrals themselves (often tens of cm³/mol). Limitations: the correction assumes the box is large compared with the range of correlations, and integrals still need R well beyond the last peak. Coordination numbers of the first shell are not affected noticeably.

Step 3

Diffusion: the Yeh–Hummer correction by hand

Why: A moving particle drags the surrounding liquid along. In a periodic box, the flow it creates also pushes on its own images, and the periodic flow field slows it down. Hydrodynamics gives the size of the effect: D in the box is lower than the true D by ξk_BT/(6πηL), where η is the viscosity and ξ = 2.837297 for a cubic box.

What the numbers say: The correction does not depend on the particle: every species gets the same absolute shift. It falls only as 1/L, so doubling the box (8 times the atoms) only halves it: for the ionic liquid a 20 nm box, about a million atoms, still underestimates D by 6 %. It matters most where it is least expected: the absolute correction is smaller in a viscous liquid, but D is smaller too, so the relative error is larger (25 % vs 11 %). Limitation: it needs the viscosity, ideally computed with the same force field (next lab), and it assumes a cubic box and a simple liquid.

Step 4

Measuring the 1/L law: four box sizes

Why: Is the correction real, and is it the right size? Run the same 3D LJ liquid in boxes of 32, 108, 256 and 500 atoms, measure D in each with the MSD (several independent runs, more for the noisy small boxes), and plot D against 1/L.

What the numbers say: D grows steadily with box size; even the 500-atom box is 11 % below the infinite-box value 0.0688. The three larger boxes fall on a straight line in 1/L whose slope gives η = 2.43, the same as the viscosity computed directly in the next lab (2.42), and the Yeh–Hummer correction brings all three to 0.0687–0.0688. The 32-atom box (L = 3.4, only about three molecular diameters) does not follow the law: including it in the fit gives a wrong viscosity (1.50) and a wrong extrapolation (0.0743). The hydrodynamic picture needs a box several molecules wide. Practical rule: correct one reasonably sized box with a viscosity computed independently, or extrapolate from several boxes that are all large enough, never from a tiny one.

Step 5

Consequences for conductivity and the Haven ratio

Why: Transport numbers built from self-diffusion inherit the Yeh–Hummer error. The Nernst–Einstein conductivity σ_NE uses D₊ and D₋ directly, and the Haven ratio σ_NE/σ compares it with the collective conductivity.

What the numbers say: A 14 % change in σ_NE moves the Haven ratio from 1.50 to 1.71, larger than many differences reported between force fields. The collective conductivity is expected to need a much smaller correction, because the hydrodynamic flow carries cations and anions together and largely cancels in the charge current; but this depends on the system (masses, sizes, how the box momentum is handled), so check it by running two box sizes. Either way: state the box size with any D, σ_NE or Haven ratio, and say whether D was corrected.

Step 6

Choosing the box: a checklist with numbers

Why: Box size is a trade-off between cost (N grows as L³) and the size-dependent errors above. Putting the main ones side by side helps choose L before a project starts.

What the numbers say: Some limits are hard: the cutoff (and any structure you want in g(r)) must stay below L/2, and no wave longer than L exists, so S(k) below 2π/L and slow collective modes are missing. Others are corrections you can calculate: D (Yeh–Hummer), charged boxes (Lab 5.5), long integrals of g(r) (Ganguly–van der Vegt). And some need checking: conductivity, viscosity and anything involving clusters or interfaces. A practical habit: run the key observable in two box sizes; if it changes by more than its error bar, either correct it or go bigger.

Summary