This repository contains an implementation of a simple Molecular Dynamics simulation of a Lennard Jones potential. It follows the standard MD workflow, beginning with random particle initialization in a cubic simulation box with periodic boundary conditions. The system advances in time using the Velocity Verlet integration algorithm, which requires force computation at each timestep.
The force computation employs a cell list algorithm, a spatial decomposition technique that reduces computational complexity from O(N²) to O(N) for short-range interactions. The simulation domain is divided into cubic cells of size equal to the interaction cutoff radius, where each particle interacts only with particles in its own cell and 13 neighboring cells using an L-shaped optimization pattern, dramatically reducing the number of pair evaluations compared to a naive O(N²) approach.
The force computation kernel implements the standard 12-6 Lennard-Jones potential given by:
\[V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]\]where $\epsilon$ represents the potential well depth, $\sigma$ the zero-potential distance, and $r$ the inter-particle separation. The force is computed as the negative gradient of this potential, with a configurable cutoff radius of 2.5σ. Periodic boundary conditions are enforced using the minimum image convention, which wraps interactions across domain boundaries. The implementation uses arithmetic masks instead of branches for cutoff enforcement, enabling branch-free vectorization.
Additionally, the data structure employs a Structure of Arrays (SoA) layout, where separate contiguous arrays store x, y, z positions, velocities, and forces.
You may enable/disable OpenMP and energy checks in the Makefile. Set your compiler with appropriate flags. Then simple run make.
$ make
For a debug build run
$ make debug
This will generate a binary called md_cell and md_cell_debug respectively.
The number of timesteps and the interval for cell list rebuilding can be set in the “simulation.params” file. The physical problem parameters can be changed in the header include/particle.h.
Please execute the binary as follows to set active thread count and affinity.
$ OMP_NUM_THREADS=xx OMP_SCHEDULE=xx OMP_PLACES=cores OMP_PROC_BIND=close ./md_cell
Currently only the main force computation loop uses schedule(runtime). Both OMP_SCHEDULE=guided,1 and OMP_SCHEDULE=dynamic seem to solve the load balance issue but the former seems to have less overhead.