Effective way to build the neighbor's list in MD

Question

I'm trying to implement the following form of the cell/neighbor list method in my MD code. I have divided my simulation box into a fixed number of cells, and according to its positions, I have correctly assigned each atom to the appropriate cell. I build the neighbor list for each individual cell $i$. When performing calculations, I only consider the atoms within the neighboring cells; furthermore, on this subset, I perform the Verlet search where I only consider the atoms within radius $r_{cutoff}$ of atom $k$. This even smaller subset of atoms will be interacting. The algorithm I want to implement is found in this paper, algorithm C 1. Algorithm 1 is the pseudocode for it.

My questions are regarding the building of the neighbors list. Pictured above is the simulation box represented in 3 layers (picture them on top of each other to form a 3D box).

1. I know that we should not consider all the neighboring cells because that would lead to double counting the interactions. Is it correct to only consider the L shaped cells I have drawn? In the figure, we are considering the atoms in cell 14 and are interested in its neighbors.
2. Currently, I can make a neighbors list if I use A LOT of if/else statements to get all the boundary conditions correctly. I have been thinking about a more efficient/concise algorithm by detecting patterns but I couldn't come up with anything. Is there an algorithm that is more efficient than doing a bunch of if statements? Perhaps one that automatically accounts for the periodic boundary conditions too?

Essentially, I'm asking how I can build the neighbors list. For this part, I'm looking to have something like int neighbors[N_c][26], where N_c is the number of cells, and 26 is how many neighbors each cell $i$ has. Data type is int because cells are identified by an int. neighbors[1] returns the array of neighbors of cell 1.

Answer 1

There are open-source MD codes, and also text books, that describe how this can be done in general. I don't want to recommend any particular one, but you can search "molecular simulation" or something similar. It might save you some time in your program design.

In your diagram you are double counting cells 5 and 23. The rule is, you need to look at neighbouring cells with relative positions $(\Delta n_x,\Delta n_y,\Delta n_z)$ where the relative cell indices take all possible combinations of the values $+1$, $0$, $-1$; except that if you look at $(\Delta n_x,\Delta n_y,\Delta n_z)$, you don't need to look at $(-\Delta n_x,-\Delta n_y,-\Delta n_z)$, so as to avoid double counting. This gives 13 neighbouring cells that need considering, and you can draw up a permanent table of the necessary relative index triplets in the computer code. And, as you realize, you need to look at the central cell, 14 in your diagram, for which $(\Delta n_x,\Delta n_y,\Delta n_z)=(0,0,0)$: when you start looking at the coordinates of pairs of atoms, it needs slightly different treatment to avoid double counting.

This leads to my main suggestion: give the cells a 3D integer index $(n_x,n_y,n_z)$, and store the cell contents (e.g. positions, or a pointer to a list of positions) in an array (or similar) using these indices. Then, identifying the neighbouring cells is a simple matter of adding $(\Delta n_x,\Delta n_y,\Delta n_z)$ to $(n_x,n_y,n_z)$ to examine each neighbouring cell in turn; the modulo or mod function will handle the wraparound of the indices associated with periodic boundary conditions. In other words, calculate modulo(nx+delta_nx,nc) (or whatever the syntax is in your preferred language) where nx runs from zero to nc-1 and there are nc cells along each coordinate direction; similarly for ny and nz.

When it comes to applying the periodic boundary conditions to the vectors between pairs, it is possible to avoid if statements by using arithmetic statements such as xij = xij - rint(xij/box)*box where xij is an uncorrected difference in $x$-coordinates between atoms i and j, box is the box length and rint is a function that gives the nearest integer to its argument (again, this may be different for different programming languages).

Answer 2

If this is a file that VMD recognizes, you can create the cell neighbor list with a tcl script that simply defines the total volume of the box you are working with and the cell side-length. Then in 3 nested loops (for x, y, z axis), shift along the whole volume of the box (working in layers) and select all the atoms within the loop-defined cell coordinates and writing them to an external list.