Lennard-Jones Gas 3D: Condensation, Pressure Gauge and Thermometer

Setting up the lab bench...

Presets

Live readouts

Graphs

Energy per atom against time: kinetic, potential and total

Temperature and pressure against time

Radial distribution function g(r)

Speed distribution against Maxwell-Boltzmann

Mean squared displacement (diffusion)

The pair potential U(r) and force F(r)

Numerical checks

running...

How to use

Pick a preset or build your own gas. Experiments: (1) Dilute gas and compare the pressure with the ideal value ρT in the readout (Z close to 1); (2) Condensing and watch the potential energy fall and clusters form while the box stays dilute; (3) compress the density slider and see P rise steeply once atoms touch; (4) Crystal then press Heat x1.6 several times to melt it and watch the g(r) peaks wash out and the MSD start to grow; (5) switch the thermostat off to see energy conservation. Drag to orbit, wheel to zoom. Space pauses, R resets, click and hold on the box to heat a spot.

The 3D page uses a Blender-built gas cell with glass walls, a pressure gauge, a thermometer and a heater plate whose parts react to the live simulation.

The science

Atoms of a noble gas attract each other weakly at distance (induced-dipole, or London, dispersion forces) and repel strongly when their electron clouds overlap. The Lennard-Jones pair potential captures both with two terms:

U(r) = 4 eps [ (sigma/r)¹² − (sigma/r)⁶ ] F(r) = −dU/dr = (24 eps / r) [ 2 (sigma/r)¹² − (sigma/r)⁶ ]

The minimum is at r = 2^1/6 sigma = 1.122 sigma, where U = −eps and F = 0. Everything is in reduced units (sigma = eps = atom mass = k_B = 1), so one run describes argon, krypton or xenon alike; for argon sigma = 0.34 nm, eps/k_B = 120 K and the time unit is 2.2 ps. Here the potential is truncated at the cut-off r_c and shifted by U(r_c) so that the energy is continuous.

Molecular dynamics. Newton's equations are integrated with the velocity-Verlet scheme, which is time-reversible and symplectic: it does not drift in total energy, it only fluctuates. With the thermostat off the run is a micro-canonical (NVE) system, and the readout shows the relative energy drift.

x(t+dt) = x + v dt + ½ a dt² v(t+dt) = v + ½ [a(t) + a(t+dt)] dt

Temperature and pressure. Temperature is defined from the kinetic energy by equipartition, T = 2 K / (d N − d) (d momentum constraints removed). A weak Berendsen thermostat rescales the velocities towards the target, v → v sqrt(1 + (dt/tau)(T0/T − 1)); it fixes the mean temperature but not the exact fluctuations of a canonical ensemble. The pressure comes from the virial theorem:

P = [ N k_B T + (1/d) Σ_i<j r_ij · F_ij ] / V ⇒ P = ρ T + (virial term)

For a dilute gas the leading correction to the ideal law is the second virial coefficient, P / (ρ T) = 1 + B₂(T) ρ + …, with B₂ = −(1/2) ∫ (e^−U/T − 1) d^dr. B₂ is negative at low temperature (attraction wins, the pressure is below ideal) and positive at high temperature (repulsion wins). The checks below compute B₂ by numerical quadrature and compare it with the pressure measured in the simulation.

Structure and dynamics. The radial distribution function g(r) counts neighbours at distance r relative to an ideal gas of the same density: a gas has g ≈ 1, a liquid has a first peak near sigma and decaying waves, a crystal has sharp peaks at the lattice distances. The mean squared displacement grows as ⟨Δr²⟩ = 2 d D_s t in a fluid (Einstein relation, D_s is the self-diffusion coefficient) and stays flat in a solid. The speed histogram converges to the Maxwell-Boltzmann law (2D: f = (v/T) e^−v²/2T, 3D: 4π v² (2πT)^−3/2 e^−v²/2T).

Honest simplifications. A few hundred atoms in a periodic box, so finite-size effects are visible near the liquid-gas transition; short runs, so the 'state' label is a heuristic from g(r), density and diffusion, not a thermodynamic phase diagram; the potential is pairwise and spherical (no three-body terms, no molecular rotation); the O(N²) force loop uses the minimum-image convention and r_c < L/2 is assumed. Quantitative critical values (3D: T_c ≈ 1.31, ρ_c ≈ 0.31 for the full potential) are not reproduced by such small systems.