MD-Sim

Custom Hybrid CPU/GPU Molecular Dynamics Engine

A high-performance Molecular Dynamics (MD) engine written in C++ and CUDA. It uses a hybrid architecture where bond/angle forces are calculated on the CPU (OpenMP) and non-bonded interactions (Lennard-Jones + Coulomb) are accelerated on the GPU.

Features:

Forces are computed independently on CPU/GPU and summed into a shared force array before the integration step.

GPU Neighbor Searching (Cell-Linked List):

To avoid the $O(N^2)$ cost of checking every atom pair, I used a spatial hashing algorithm entirely on the GPU.

1. The simulation box is divided into cubic cells of side length $d\ge r_{cut}+\text{skin}$. Each atom is hashed to a 1D cell index $k=x_{grid}+y_{grid}\cdot D+z_{grid}\cdot D^2$.
2. Atoms are sorted by cell ID using CUB Radix Sort.
3. A lookup table (cell_start, cell_end) is generated to map cell IDs to start/end indices in the sorted atom array.
4. For every atom, the kernel searches only the 27 adjacent cells (3x3x3 block). Neighbors within $r<r_{cut}+\text{skin}$ are stored in a fixed-width adjacency list per atom.

Equations:

1. Equations of Motion - Langevin Dynamics

The engine uses Langevin Dynamics instead of standard Newtonian physics ($F=ma$). This simulates a system connected to a heat bath (thermostat) by adding friction and random thermal noise. $m_i\frac{d{\mathbf{v}}_i}{dt}=\underset{\text{Forces from Atoms}}{\underbrace{{\mathbf{F}}_{conservative}({\mathbf{r}}_i)}}-\underset{\text{Friction}}{\underbrace{\gamma {\mathbf{v}}_i}}+\underset{\text{Random Noise}}{\underbrace{{\mathbf{R}}_i(t)}}$

2. Non-Bonded Interactions (GPU)

Calculated between all pairs of atoms (excluding bonded neighbors) that fall within the cutoff radius ($r<r_{cut}$).

A. Lennard-Jones Potential

$V_{LJ}(r_{ij})=4{ϵ}_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right]$

B. Coulomb Potential (Electrostatics)

$V_{Coulomb}(r_{ij})=\frac{1}{4\pi {ϵ}_0}\frac{q_i{q}_j}{r_{ij}}$

3. Bonded Interactions (CPU)

A. Harmonic Bonds (2-Body)

$V_{bond}(r)=\frac{1}{2}{k}_b(r-r_0{)}^2$

B. Harmonic Angles (3-Body)

$V_{angle}(\theta )=\frac{1}{2}{k}_{\theta }(\theta -{\theta }_0{)}^2$

4. Integration Scheme (Velocity Verlet)

The engine solves these differential equations using a discrete time step ($\Delta t$).

1. Half-Kick: $\mathbf{v}(t+\frac{1}{2}\Delta t)=\mathbf{v}(t)+\frac{1}{2}\mathbf{a}(t)\Delta t$
2. Drift: $\mathbf{r}(t+\Delta t)=\mathbf{r}(t)+\mathbf{v}(t+\frac{1}{2}\Delta t)\Delta t$
3. Compute Forces: $\mathbf{a}(t+\Delta t)=\mathbf{F}(\mathbf{r}(t+\Delta t))/m$
4. Half-Kick: $\mathbf{v}(t+\Delta t)=\mathbf{v}(t+\frac{1}{2}\Delta t)+\frac{1}{2}\mathbf{a}(t+\Delta t)\Delta t$

Resources:

• Used this as a starting point: https://youtu.be/ChQbBqndwIA?si=14Vq-WGgh8yWk7DQ
• C++: https://www.learncpp.com/
• Langevin Dynamics: https://en.wikipedia.org/wiki/Langevin_dynamics
• CUDA basics: https://youtu.be/xwbD6fL5qC8?si=ivaTzz3BVbgJdFsF

Github