Some papers change everything. Aneesur Rahman’s 1964 work, “Correlations in the Motion of Atoms in Liquid Argon”, is one of them. Just eight pages in Physical Review, but it fundamentally invented molecular dynamics simulation.
Rahman’s idea was radical for its time: use a computer to solve Newton’s equations for 864 argon atoms and observe what happens. This was 1964—computers filled entire rooms, and simulating individual atoms seemed like science fiction. But it worked, brilliantly.
For the first time, researchers could “see” how atoms move in a liquid. The results weren’t just computational curiosities—they revealed fundamental physics. Rahman discovered that atoms in liquids don’t simply diffuse randomly; they get trapped by their neighbors and bounce back, creating the “cage effect” that remains central to our understanding of liquid dynamics.
I wanted to test whether Rahman’s results would hold up using modern computational tools. This project, available on GitHub, combines Python for analysis with LAMMPS for simulation. The question was straightforward: run the exact same physical conditions as 1964 and see what emerges.
The Setup: Modern Tools, Classic Physics
I preserved Rahman’s original conditions while leveraging modern computational advantages. The core parameters in my LAMMPS input script match his exactly:
- System: 864 argon atoms
- Interaction: Lennard-Jones potential with $\sigma = 3.4$ Å and $\epsilon/k_B = 120$ K
- Thermodynamic state: Target temperature of 94.4 K and density of 1.374 g/cm³
However, I incorporated modern best practices unavailable in 1964:
- Energy minimization before dynamics to eliminate poor contacts from the initial crystal lattice
- Proper equilibration with 500 ps in the NVT ensemble, allowing the crystal to melt into a realistic liquid
- Improved integration using Velocity Verlet with a 2 fs timestep instead of Rahman’s custom 10 fs method
The production run lasted 10 ps in the NVE ensemble, generating 5,001 frames. Temperature remained within 1% of target with an RMS fluctuation of just 0.0165—substantially more stable than Rahman could achieve with 1960s hardware.
The Results: 60 Years Later
My Python analysis toolkit processes the LAMMPS trajectory to recreate each figure from the original paper. The agreement is remarkable.
System Check: Temperature and Velocities
First, the basic checks. Figure 1a shows temperature fluctuating around our 94.4 K target—exactly what we want for proper thermal equilibrium. Mean temperature was 94.73 K (just 0.33 K off target) with a standard deviation of 1.56 K. That’s significantly better stability than Rahman achieved with 1960s computational resources.
Figure 1b shows the velocity distribution from 12.9 million velocity components. The result is a perfect Maxwell-Boltzmann distribution, as expected for thermal equilibrium. The distribution widths at various heights match Rahman’s results remarkably well: 1.77, 2.48, and 3.56 compared to his 1.77, 2.52, and 3.52.
Liquid Structure
How are atoms arranged in a liquid? The Radial Distribution Function $g(r)$ (left panel) provides the answer. It shows the probability of finding another atom at distance $r$ from any reference atom.
The sharp first peak at 3.82 Å reveals well-defined nearest neighbors, but subsequent peaks rapidly dampen—the characteristic signature of liquid structure. Short-range order persists, but long-range order disappears. Rahman found the first peak at 3.7 Å; our 0.1 Å difference is negligible given the 60-year gap in computational methods.
The Static Structure Factor $S(k)$ (right panel) is the Fourier transform of $g(r)$. This quantity directly corresponds to what you’d measure in X-ray or neutron scattering experiments. Peak locations at 6.71, 12.6, 18.2, and 24.94 (in units of $k\sigma$) match Rahman’s expected values of 6.8, 12.5, 18.5, and 24.8 with remarkable precision.
Diffusion: How Atoms Wander
Track how far atoms move over time and you obtain the Mean Square Displacement (MSD). Figure 3 shows the classic transition from ballistic motion (short times, where atoms move in straight lines) to diffusive motion (long times, where the linear regime emerges).
The slope of the linear region gives the diffusion coefficient: $D = 2.47 \times 10^{-5} \text{ cm}^2/\text{s}$ from my simulation versus Rahman’s $2.43 \times 10^{-5} \text{ cm}^2/\text{s}$. That represents 2% agreement after 60 years of technological advancement—remarkable consistency for completely different hardware and software implementations.
The Cage Effect: Rahman’s Big Discovery
Here’s where Rahman made his most profound discovery. The Velocity Autocorrelation Function (VACF) in Figure 4 measures how long an atom “remembers” its initial velocity direction.
In a gas, this function would decay exponentially to zero. In Rahman’s liquid, something fundamentally different occurs: the VACF goes negative around 0.3 ps before slowly returning to zero. This is the famous “cage effect”—atoms become trapped by their neighbors and bounce backward, temporarily moving opposite to their initial direction.
My simulation captures this phenomenon perfectly. The VACF reaches -0.020 at 1 ps, with a minimum of -0.083. This was Rahman’s key insight about liquid dynamics, and it remains fundamental to our understanding of how liquids behave at the atomic scale.
The Fourier transform of the VACF (Figure 5) reveals something even more intriguing. Instead of the simple curve predicted by classical diffusion theory, there’s a distinct peak around $\beta = 0.25$ (where $\beta = \hbar\omega/k_B T$). This peak represents atoms “rattling” in their cages—behavior more reminiscent of a solid than a gas. Rahman observed this feature too, and it became a key signature of liquid dynamics that bridges the gap between crystalline and gaseous behavior.
Advanced Analysis: Van Hove Functions
Rahman didn’t limit himself to basic properties. The Van Hove correlation function $G(r,t)$ describes how liquid structure evolves over time. Figure 6 shows the distinct part—how neighbor coordination shells “melt” as time progresses.
At 1.0 ps, the structure remains well-defined with clear shells. By 2.5 ps, it becomes increasingly diffuse. Rahman compared this evolution to theoretical predictions (the Vineyard approximation) and found that theory predicted overly rapid structural decay. My results confirm this finding—the Vineyard approximation underestimates how persistent local structure remains in liquids.
Figure 7 demonstrates that atomic motion deviates significantly from simple random diffusion. The non-Gaussian parameters $\alpha_n(t)$ measure deviations from Gaussian displacement behavior. If atoms moved via simple random walks, these parameters would be zero.
Instead, they’re positive, with $\alpha_2$ reaching 0.115. This indicates that some atoms make larger jumps while others remain more confined than simple diffusion theory predicts. Rahman was the first to identify this non-Gaussian character in liquid dynamics, revealing the heterogeneous nature of atomic motion in liquids.
Figure 8 tests Rahman’s “delayed convolution approximation”—his proposed improvement over existing theory. By comparing actual structure evolution with theoretical predictions, this shows how simulation can validate and improve theoretical models. The delayed approach works better than the standard Vineyard method, demonstrating how computation and theory can work together.
Implementation Notes
Beyond reproducing the physics, this project represents an exercise in modern scientific software engineering:
- Clean architecture: Analysis code organized into a Python package (
src/argon_sim) with dedicated modules for thermodynamics, dynamics, and structure calculations - Intelligent caching: Expensive computations like MSD and VACF are cached to prevent redundant calculations
- Memory optimization: Large datasets processed in chunks to efficiently handle the 864 × 5,001 atom-frame matrices
- Full reproducibility: A
Makefileenables one-command installation of dependencies, simulation execution, and figure generation
What required months of computation time on Rahman’s IBM mainframe now completes in less than an hour on a modern laptop, with superior numerical accuracy and statistical precision.
Why This Still Matters
Replicating Rahman’s work transcends mere historical curiosity. The quantitative agreement is striking: diffusion coefficients within 2%, structural peaks within 0.1 Å, velocity distributions matching to three significant figures.
This level of reproducibility, achieved with completely different hardware and software, validates something fundamental. Rahman’s physical model was remarkably sound—the Lennard-Jones potential captures liquid argon physics with extraordinary fidelity. His computational methodology, while constrained by 1960s technology, was scientifically rigorous.
Most importantly, Rahman didn’t simply simulate a liquid—he established the conceptual foundation for modern computational materials science. Molecular dynamics has evolved from 864 atoms to millions, but the core insight endures: Newton’s equations plus statistical mechanics provides a direct window into the atomic world.
This replication demonstrates that some physics is truly universal. The cage effect, velocity correlations, and structural evolution aren’t computational artifacts from 1960s hardware. They’re fundamental characteristics of how matter behaves at the atomic scale, as relevant today as they were six decades ago.