MD Simulation: Reference for the Neighbor's List Method

Question

With a rather basic knowledge in C++, I have written my own MD simulation code. Currently, I calculate forces in the most naive way: I go through all the atoms and account for their interactions. This of course is an $\mathcal{O}(N^2)$ algorithm, and very inefficient. It wouldn't be an issue, except for my research purposes, I need to simulate 1000 or more atoms, and currently this would take a very very long time.

I wanted to learn how to implement the neighbor's list method. I have read some of the standard MD books, but they are often in Fortran or some language I cannot decipher. I'm looking for a rather modern and easy-to-follow guide on how to implement this algorithm. It's not an important part of my research so I would like to spend as little time on it as possible. I have to admit, I'm basically asking if there's a quick and dirty way I can learn and implement this algorithm! Does anyone know of such a guide? Or has any experience with implementing the neighbor's list method?

Answer 1

I'd recommend "The Art of Molecular Dynamics Simulation" by D. C. Rapaport. The code samples are written in C. I'm not a huge fan of the programming style of the book, but at least it's not FORTRAN. Having said that, my advice would be to take any book where neighbour lists are explained (for instance the Frenkel & Smit book, which I guess it's what you are using now) and just implement the algorithm, which is usually written in pseudo-code.

If you want a TL;DR, there are several types of neighbour lists. All of these rely on your interaction potential being 0 outside a certain range $r_c$.

• Verlet lists ($\mathcal{O}(N^2)$ - $\mathcal{O}(N^{3/2})$ but with a small-ish prefactor):
  1. do an $N^2$ sweep where for each particle you fill an array of neighbours that are closer than $r_c + r_v$, where $r_v$ is a parameter and save the current position of each particle (let's call it $\vec{r}_i(t_0)$.
  2. Evolve your simulation by using the list of neighbours to compute the forces.
  3. After each integration step check whether, for any particle $i$, $|\vec{r}_i(t) - \vec{r}_i(t_0)| \geq r_v / 2$. If it is then update all the lists (see step 1).
• Cell lists ($\mathcal{O}(N)$):

1. Divide your simulation box in cells of linear size $l > r_c$.
2. Use linked lists to assign particles to the cells.
3. For each particle $i$, the list of neighbours is given by all the particles that are in $i$'s cells and in each of the (8 in 2D and 26 in 3D) neighbouring cells. Calculate the forces between $i$ and its neighbours.
4. After each integration step check whether particles have crossed cell boundaries and update the data structure accordingly

• Verlet lists built with cell lists (usually the best choice, $\mathcal{O}(N)$ with a smaller prefactor than just using cells):

1. Step 1 of the "Verlet lists" section is carried out by using cell lists. The linked lists that store the cell data structures are kept updated throughout the simulation

The reason why using Verlet lists instead of cell lists is (in MD simulations) better is that the average number of neighbours is smaller for the former than the latter. The difference is due to the fact that, for each particle, cell lists give you a list of neighbours that are in a volume $(3l)^3$, which, if $r_v$ is chosen wisely, is much larger than the volume of the Verlet sphere ($4/3 \pi (r_c + r_v)^3$). Therefore, on average, you compute much fewer distances with Verlet lists.