Over the summer of 2012, I worked at the Naval Research Laboratory in Washington, DC, just as I did in 2011. In fact, I even had the same mentor. However, unlike last time, my project didn’t really have to do with the high performance computation resources the Naval Research Laboratory has. Instead, my project focused on molecular dynamics – a branch of computation physics that treats atoms as particles behaving under the Lennard-Jones potential and simulates their movement in discrete time steps.

Introduction

In general, there are a couple things we want with a molecular dynamics simulation. Obviously, we want it to be accurate. However, we also need it to be fast. Each particle interacts with every other particle at once, meaning that for a simulation with $N$ particles, we could expect to worry about $O(N^2)$ interactions. Thus, efficiency is what my project focused on; particularly, what algorithm leads to the most efficient molecular dynamics simulations (while not sacrificing much accuracy)?

Algorithms

The general way these algorithms work is that there are a bunch of particles, which start off in some lattice configuration and presumably some random velocity to represent temperature. The simulation is divided into discrete time steps, where, for each timestep, the force between every pair of particles is computed and the net force for each particle is summed up to determine the particle’s current acceleration. The acceleration of course will affect the velocity and thus position, but these are first and second order effects that appear only after many time steps. Over time, the particles should bounce around near each other, without ever getting too close.

Lennard-Jones Potential. As is standard, the cutoff is chosen at r = 2.5σ.

As mentioned earlier, the largest inefficiency comes from having to do $O(N^2)$ calculations of force for each time step. Considering that we only start to observe interesting behavior when $N$ is sufficiently large, it makes sense that we would want to improve the efficiency of this step as much as possible. It turns out we can do this by assuming the LJ potential only acts over a small distance (which is a good approximation, since the end behavior scales as $r^{-12}$). In a sense, this reduces the problem to $O(N)$ assuming $r_c$ is much smaller than the dimensions of the simulation. Note that $r_c \approx 2.5 \sigma$ is the cutoff radius of the LJ potential.

While in theory, this is a good solution, the main problem is figuring out which particles are close to a given one – that is, which particles are within $r_c$ of the particle of interest. My goal was to experiment with different integration step techniques, including the ones shown below. These weren’t the only ones, but I risk revealing too much information if I go into all the algorithms I tried.

Brute Force (all-pairs)

Ok, this one is really simple. For each pair of particles, check to see if they’re within $r_c$ of each other and if they are, sum the forces as usual. Basically this one is still $O(N^2)$ because we do an $O(N^2)$ check to find the particles within $r_c$ of the particle of interest. Way too slow.

Linked-Cell

Here we get a little more interesting. We create a secondary structure for the particles, namely a grid of cells (cubes) of side length $r_c$. Then, for each particle, we only check if particles within the current or one of the $26$ neighboring cells is within $r_c$ of the particle of interest for doing MD force calculations. When we move each particle, we check to make sure the particle hasn’t escaped its current cell, and if it has, we make sure the linked-cell structure knows about it.

Monotonic Lagrangian Grid

This method is really interesting, though I didn’t actually get much success with it, if I recall correctly (either it didn’t work or was too slow). So essentially, you pretend that the particles are always in a 3D grid-like pattern and try to enumerate them as if that was the case. Of course, the particles are in reality, way off from the perfect grid-like structure. However, the MLG method allows us to sort of fit an enumeration to the particles so that particles that are close to each other in enumeration are also close to each other in real space. The supposed benefits of this method are that it’s light-weight and can easily be used to retrieve the nearest neighbors (supposedly).

Results

Well, I can’t actually say what the experimental results were due to security restrictions (I don’t really see how this would be of any threat, but I guess I’ll play it safe). I can show a very simple demo of my application working on one of the algorithms described above though, since all of this stuff has been done before.

NRL Pictures

Finally, since this was my last summer at NRL, I thought it would be fitting to add some pictures of the place.