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.

  • $\begingroup$ What kind of neighbour list do you want to construct? Verlet, cell,...? $\endgroup$ – Anton Menshov Jul 13 '18 at 21:39
  • $\begingroup$ It's both combined. I create the cell and then do verlet within the neighboring cells. Also, I'm simulating 1000-2000 atoms in Langevin Dynamics. I don't believe I'm doing anything fancy and also I'm not too advanced so I would like to stay with basics if possible. $\endgroup$ – Ptheguy Jul 13 '18 at 22:21

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).

  • $\begingroup$ I find your suggestions helpful. But, since the posting of this question, I have been working on the algorithm myself and now I have a very efficient way to automatically assign the neighbors of each cell. The thing is that I assign all the neighbors--total of 26 (excluding the cell we're focusing on). There is double counting in all the cells. Can I just divide my potential and force calculations by 2? I really like the way I wrote the algorithm and prefer not to have to rethink it. If not, I could just maybe work in your suggestions. $\endgroup$ – Ptheguy Jul 16 '18 at 18:34
  • 1
    $\begingroup$ You can certainly proceed as you suggest. The program is likely to be slower, up to a factor 2, because generally most computational work lies in the calculation of potentials and forces, rather than in the construction of lists. But this may not be critical to your application. $\endgroup$ – user28077 Jul 16 '18 at 19:37

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.

  • $\begingroup$ This is a c++ code that I have written from scratch. Thus, I need to build the cell/verlet list myself as well. $\endgroup$ – Ptheguy Jul 14 '18 at 19:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.