Home Research Project Teaching Seminar Notes CV Blog

FCC Molecular Dynamics Simulator

A molecular dynamics simulator built in Python and C++ to explore numerical methods in atomistic simulation and scientific software development.

I independently developed this project in order to connect the continuum modeling and numerical simulation in my liquid crystals research with modeling and simulation at the atomistic scale. I also wanted experience building scientific simulation software. This simulator initializes a chosen face-centered cubic crystal lattice, evolves particle trajectories using Lennard-Jones pair potentials and Newton's laws of motion, and computes structural, thermodynamic, and dynamical properties of the material. I implemented the core simulation algorithms, analysis routines, automated reports and plots, and a Streamlit web interface. I then accelerated the computational kernels using C++ and pybind11.

GitHub · Interactive Streamlit Demo · Documentation

FCC molecular dynamics Streamlit interface showing simulation controls and a nickel crystal. — Simulation controls and visualization of initial 4 × 4 × 4 Ni FCC crystal.

Technical Highlights

Velocity Verlet integration for NVE molecular dynamics
Berendsen thermostat for NVT simulations
Periodic boundary conditions with minimum-image convention
Verlet neighbor lists for efficient force calculations and time integration
Python and C++ simulation backends connected using pybind11
Automated structural, thermodynamic, and transport analysis
100+ automated tests for numerical and physical validation

Numerical Methods and Implementation

Particle trajectories are evolved using the velocity Verlet algorithm with periodic boundary conditions and the minimum-image convention. Pair interactions are modeled using a 12-6 Lennard-Jones potential, while Verlet neighbor lists reduce the cost of identifying interacting particle pairs.

The original implementation was written in Python. I subsequently rewrote the computationally expensive force, integration, and neighbor-list kernels in C++ and exposed them to Python using pybind11. This preserves the Python analysis and visualization workflow while moving the performance-critical numerical kernels to compiled code.

Verification and Validation

I developed an automated validation suite to test both the numerical implementation and the expected physical behavior of the simulator. Tests include:

Energy conservation in NVE simulations
Temperature stability and equipartition in NVT simulations
Linear momentum conservation
Timestep refinement and convergence studies
Recovery of the expected FCC coordination number of 12
Agreement between independent Python and C++ implementations
Automated unit and regression testing with pytest

Energy conservation validation — **Energy conservation.** Relative energy drift remains well below the prescribed tolerance during a 10 ps NVE simulation.

Time step convergence validation — **Time step convergence.** The observed convergence order of 1.85 is consistent with the expected second-order accuracy of velocity-Verlet.

Simulation Analysis

Post-processing tools compute structural, thermodynamic, and dynamical properties from particle trajectories, including the radial distribution function, structure factor, mean squared displacement, velocity autocorrelation function, coordination number, pressure, heat capacity, and diffusion coefficients.

Radial distribution function for nickel at 300 K. — **Radial distribution function.** The pronounced coordination-shell peaks reflect the ordered structure of the FCC crystal.

Temperature versus time for an NVT nickel simulation. — **Temperature.** The NVT simulation fluctuates around the target temperature of 300 K.

Mean squared displacement for nickel at 300 K. — **Mean squared displacement.** Atomic displacements remain bounded at room temperature, as expected for atoms vibrating about lattice sites in the solid phase.

Velocity autocorrelation function for nickel at 300 K. — **Velocity autocorrelation function.** The decaying oscillations characterize atomic vibrations and loss of velocity correlation.

Performance

To explore performance-oriented scientific computing, I implemented a second simulation backend in C++ and connected it to Python using pybind11. The force evaluation, neighbor-list, and integration kernels were moved to compiled code while preserving the Python-based analysis and visualization workflow.

Benchmarking showed large improvements in both kernel throughput and end-to-end simulation runtime compared with the original non-vectorized Python implementation.

Lennard-Jones force throughput comparison between Python and C++ backends. — **Force-kernel throughput.** The C++ backend evaluates substantially more neighbor-pair interactions per second than the original Python implementation across all tested system sizes.

Runtime comparison for Python and C++ NVE simulations. — **End-to-end runtime.** For 500 NVE steps, the C++ backend achieves substantially lower runtimes across all tested system sizes than the Python backend.

Technologies

Python · C++ · NumPy · pybind11 · pytest · Matplotlib · Streamlit