Molecular Dynamics Simulation - 3D Ball-and-Stick Reaction Chamber Lab (Water, Argon, CO2, Methane)

Setting up the lab bench...

Presets

Live readouts

Graphs

Temperature against time (green line = thermostat target)

Energy per molecule: kinetic, potential, total

Radial distribution function g(r) of the heavy atoms

Mean squared displacement of the molecules (slope = 6D)

Numerical checks

running...

How to use

Pick a preset (liquid water, argon crystal, steam at 6000 K...), or choose a substance, number of molecules, density and target temperature yourself.
Drag any atom on the canvas to pull it with a spring; the neighbours are dragged along and the thermostat takes the heat away.
Heat pulse and Quench change the speeds instantly; the thermostat then pulls the temperature back to the target. Set the thermostat to None to watch an isolated system at constant energy.
Argon crystal: raise the target temperature above about 85 K and watch the sharp g(r) peaks wash out as it melts, while the diffusion coefficient jumps.
Steam breaks apart: Morse bonds on and 6000 K, then count the broken bonds in the readout.
Turn Partial charges to zero in liquid water to see what hydrogen bonding does to g(r), the density and the energy.
Drag to orbit, wheel to zoom. Space pauses, R resets. Time warp sets how many steps run per frame; the page lowers it automatically if your device cannot keep up.

The 3D page is a Blender-built laboratory: baked-PBR bench, glass reaction chamber with a heater and a cooler plate, pressure gauge, thermometer, model-kit shelf and glassware. The molecules inside the chamber are the very same simulation as in the 2D page, drawn as ball-and-stick models.

The science

This is a classical force-field molecular dynamics (MD) engine written from scratch in JavaScript. Every atom is a point mass; the force on each one is the negative gradient of a potential energy made of bonded and non-bonded terms, and Newton's second law is integrated step by step. The same engine runs in the 2D and the 3D page, so both give identical numbers.

The potential energy

U = Σ_bonds ½k_b(r − r₀)² + Σ_angles ½k_θ(θ − θ₀)² + Σ_i<j 4ε_ij[(σ_ij/r)¹² − (σ_ij/r)⁶] + Σ_i<j q_iq_je²/(4πε₀r)

The Lennard-Jones term gives short-range repulsion (Pauli exclusion) and long-range dispersion attraction; unlike pairs use the Lorentz-Berthelot rules σ_ij = (σ_i+σ_j)/2 and ε_ij = √(ε_iε_j). The Coulomb term between partial charges is what creates hydrogen bonds in water. Atoms of the same molecule do not feel LJ or Coulomb forces from each other (they are held by bonds and angles). With the Morse bonds switch each bond becomes U = D(1 − e^{−a(r−r₀)})² with a = √(k/2D), so a bond can stretch and rupture when it is given about D of energy.

Integration and thermostats

v(t+½Δt) = v(t) + F(t)Δt/2m; r(t+Δt) = r(t) + v(t+½Δt)Δt; v(t+Δt) = v(t+½Δt) + F(t+Δt)Δt/2m

This velocity-Verlet scheme is time-reversible and symplectic, so with the thermostat off (NVE) the total energy stays constant to a few parts per 10⁵. The Langevin thermostat adds friction and random kicks, v → c₁v + √(1−c₁²)√(k_BT/m) ξ with c₁ = e^−γΔt, and samples the canonical ensemble; Berendsen rescales velocities by λ = √(1 + (Δt/τ)(T₀/T − 1)), which relaxes the temperature smoothly but does not give the exact fluctuations. The temperature is T = 2K/(3N−3)k_B.

Observables

P = (2K_cm/3 + W_mol/3)/V, W_mol = Σ_{molecule pairs} R_AB·F_AB g(r) = ⟨n(r, r+dr)⟩ / (4πr²dr ρ) D = lim_t→∞ ⟨|Δr(t)|²⟩ / 6t

The pressure is the molecular pressure: centre-of-mass kinetic energy plus the virial of the forces between molecules, which stays correct with multiple time steps. The radial distribution function g(r) uses heavy atoms (all atoms except hydrogens); a sharp first peak followed by oscillations means a liquid, a split series of peaks means a crystal, g ≈ 1 means a gas. The self-diffusion coefficient comes from the slope of the molecules' mean squared displacement (Einstein relation). Real values to compare: liquid water at 300 K D ≈ 2.3–2.5 × 10⁻⁹ m²/s, liquid argon at 90 K ≈ 2 × 10⁻⁹ m²/s (this small, cut-off system gives somewhat different values).

Units, parameters and honest limits

Internal units: Å, amu, kJ/mol; one internal time unit is 0.1 ps (100 fs). Water, CO₂, methane and nitrogen use parameters of the published SPC/Fw, EPM2, OPLS-AA and TraPPE families, simplified and rounded; argon uses σ = 3.405 Å, ε = 0.996 kJ/mol.
Electrostatics use a molecule-based cutoff with a reaction-field correction (dielectric constant 78), the usual cheap replacement for Ewald summation. LJ is force-shifted to zero at the cutoff. Systems of 27–216 molecules with periodic boundaries have finite-size effects: densities and pressures are only qualitatively right and the pressure fluctuates by hundreds of bar.
Bond topology is fixed: with Morse bonds a bond can break and the atoms fly apart, but atoms never form new partners, so this is thermal dissociation, not a reactive force field and not quantum chemistry. Classical MD also has no zero-point energy: high-frequency vibrations (O–H stretch) carry k_BT instead of the much smaller quantum energy.
Multiple time steps (r-RESPA): bonds and angles are integrated with inner steps of at most 0.5 fs, while the slow non-bonded forces are applied once per outer step (default 2 fs). This makes flexible water about four times cheaper; above about 3 fs the energy drift grows. The Verlet scheme still overestimates stiff-vibration amplitudes by 1/(1 − (ωδt)²/4), so inner steps must stay small.
Temperatures above about 400 K in liquid water without a barostat boil the box into a gas only if you also lower the density; the density slider is a piston that compresses or expands the periodic box smoothly (0.1% per step).

Built-in numerical checks

The panel at the bottom runs nine checks against analytic or independent results on small, separate simulations: the force-free LJ pair distance (the root of the shifted-force equation, which differs from 2^1/6σ by 0.24%), the N₂ bond-vibration period against 2π√(μ/k), energy and momentum conservation in NVE, a finite-difference test that the atomic virial equals −dU/dV for 64 flexible water molecules (it validates every force term, the switching function and the reaction field), the Langevin mean temperature, the Maxwell kurtosis of 3, and the ideal-gas law P V = N k_B T for dilute N₂.