# How to calculate $g_r$ when using the cell method algorithm

I have a self-written MD code in C++. I am simulating atoms using the Langevin dynamics model. I use the linked-cell method to speed up the simulation; as a result, atomic interactions are only happening among the nearest neighboring atoms. That means, when calculating the distances between atom $i$ and the rest of the system, only those atoms within the neighboring cells are considered. This naturally affects the calculation of the radial distribution function $g_r(r)$. For example, $g_r(r>r_c)=0$ uniquely for some cutoff distance $r_c$ (due to a lack of interactions). This isn't good since you want to simulate all the way up to larger distances to ensure $g_r \rightarrow 1$.

Is there a way to use the cell method and calculate $g_r(r)$ correctly? If not, am I stuck with using the traditional $\mathcal{O}(N^2)$ algorithm for this calculation, in which all atomic interactions are considered? (I simulate about 1000 atoms).

I believe that you can make a small improvement, but not a dramatic one. Assume that your force loop consists of an outer loop over cells, and an inner one over nearest-neighbour cells (and then, more loops over the atoms within those cells). For the pair distribution function calculation, the inner loop can go over not just the nearest-neighbour cells, but over all cells covering a roughly spherical region out to the desired distance, relative to the cell being considered in the outer loop. You can construct a list of these cells (with indices relative to the central cell) at the start of the program. How you do this depends on how you handle the cell indexing, as you mentioned in your earlier question.

If the maximum distance of interest is not as large as half the box length $L$, this can make a reasonable saving of time over the all-pairs method. However, if you want to calculate out to half the box length, the maximum saving is the ratio of the spherical volume you are looking at, $\frac{4}{3}\pi (L/2)^3$ to the volume $L^3$ that you would look at with all pairs: around 50%. Essentially, you are just avoiding looking in the "corners" of the periodic box surrounding a given cell. For a system of size $N=1000$, I would think that it's not worth it, and you would do best to stick to all-pairs (discarding any for which the separation is larger than $L/2$).

The silver lining here is that the statistics of the calculation improve along with the extra work done. In a sense, you get what you pay for: you look at more pairs, for the higher values of $r$, so the cost goes up, but the precision of the calculation is better for large $r$. If you wish, you can offset this by doing the calculation less often (certainly not every step!).

There is another approach, which should give you the long-range part of $g(r)$ more cheaply, but again it is probably only worth doing if you have much larger systems. You can create a 3D histogram of particle density, $\rho(\mathbf{r})$ from a snapshot of the configuration, by counting the atoms in a cubic grid of cells, but on a smaller length scale than for the force calculation. The cell size should be a fraction of the particle diameter, typically. Then perform a 3D discrete Fourier transform, using an FFT algorithm, giving the instantaneous $\hat{\rho}(\mathbf{k})$ and take the square modulus. This function of wave-vector $|\hat{\rho}(\mathbf{k})|^2$ can be averaged over snapshots taken at intervals, and eventually gives the structure factor $S(\mathbf{k})$. Inverse-Fourier-transforming this would give you $g(r)$. Again, however, I think this would be too much trouble for $N=1000$ atoms.

• Thank you for the insight and suggestions. I guess I'll stick with the $\mathcal{O}(N^2)$ operation. It will take ages though! – Ptheguy Jul 26 '18 at 21:38
• In that case, the take-home message from my answer is the simplest one: don't do this calculation every step! You lose very little by allowing the configuration to evolve over many steps between calculations of configurational properties: enough for the correlations in time to decay away. – user28077 Jul 26 '18 at 21:43
• Yeah for sure not every step. I actually do it every 100 time step. Still, given my computing powers (surface pro 3) it takes a very very long time. – Ptheguy Jul 26 '18 at 21:53

That is not how it should work. Simulations using linked-list cells (and/or Verlet lists) or $\mathcal{O}(N^2)$ calculations should (must!) give the exact same results. The only side effect should be a speed-up. You need to decouple the operations, using linked-list cells to evaluate forces and a $\mathcal{O}(N^2)$ loop to calculate the radial distribution function. This should not slow down the simulation too much, since you don't really need to update the $g(r)$ every time step. Depending on the system, doing it every $10^3$ or even $10^4$ time steps shouldn't worsen the quality of the $g(r)$ too much.

• I disagree with you on "calculations must give the exact results". (1) the cell method is an approximation, so some disagreement is inherent to the method; (2) I have outputted the potential and force calculations using both methods and they agree; (3) plots such as energy, velocity, and other parameters are in agreement with theory, as well as the results obtained via $\mathcal{O}(N^2)$ method. And as for $g_r$, the larger I make the cutoff distance the better the results are. That means the more interactions I consider the better. So I'm pretty sure it has to do with the approximation. – Ptheguy Jul 27 '18 at 17:44
• @Ptheguy Yeah I probably took for granted that you are simulating a system of particles interacting through an interaction potential that goes to 0 at a finite system $r_c$. If that's the case if you take the size of the cells to be $> r_c$ then the $\mathcal{O}(N^2)$ and the linked-list algorithms must give the same results. If that's not the case then I agree that the cell method is an approximation. – lr1985 Jul 30 '18 at 21:08