Normalization condition in evaluating the radial distribution function for a hard-spheres system

I am trying to simulate a 3D hard-spheres in a box system in Python. The idea behind this is that I have to create a system in a 1x1x1 box with periodic boundary conditions and with $$N = 500$$ particles of diameter $$d\approx0.09$$ (approximately) since packing fraction $$\phi = 0.2$$, which cannot overlap.

My algorithm:

1. Initiate a random box by picking a particle $$(x, y, z)$$ from a uniform distribution.
2. Place another particle in the box, by picking $$(x, y, z)$$ from a uniform distribution. Make sure it is not overlapping with the first one.
3. Place a third one, make sure it does not overlap with the first $$2$$.
4. Repeat this process for all particles up to $$N = 500$$.

I just have to make sure that I am imposing my periodic boundary conditions correctly.

What I have to do next is evaluate the radial distribution function. Which simply means that I have to pick a particle, evaluate its distance from every other particle, and based on the distance, I place the other particle in a bin corresponding to its distance. My bins are essentially intervals like $$[0,0.05), [0.05,0.1), [0.1, 0.15), ... [0.95, 1]$$, so a particle at a distance $$r$$ from my reference particle will go in to the bin which has $$r$$ in its interval. If $$r = 0.82$$, it will go to bin corresponding to $$[0.8, 0.85)$$. The thickness of my shell in this example is $$\Delta r = 0.05$$.

If $$n_i(r)$$ is the number of particles in the shell at a distance $$r$$ from particle $$i$$, I define $$n(r) = \sum _{i=1} ^N n_i(r)\tag{1}$$ Then I define $$\tag{2}g(r) = \frac{\frac{n(r)}{N}}{4\pi r^2 \Delta r \rho }$$ where $$\rho=N/V$$ is the number density. In my code, I am only considering the closest image of my particles.

However, when I plot $$g(r)$$ against $$r$$, I get the following plot

I kind of believe that it goes to zero because I have periodic boundary conditions, and I am only considering the closest image of my molecules. I don't think this is correct though for multiple reasons: the graph should be plateauing around 1, and it should not be going to zero. Looking at the plots here, I definitely believe I am going wrong somewhere.

Is my normalization incomplete for my radial distribution function?

Do periodic boundary conditions lead to $$g(r)$$ going to zero for hard spheres?

• As I understand, the numerator in the g(r) talks about the average number of particles in a shell. So shouldn't it be divided by the number of shells instead of total number of particles? Adding up the particles and dividing by the total number of particles seems to cancel out. Maybe there's a better way to calculate the average number of particles in a shell at a distance r. Dec 2 '20 at 6:38
• @KavyaMrudula, thanks for your answer. What do you mean by dividing by the number of shells? you mean the number of bins I have? Dec 2 '20 at 15:06
• sorry, I'm not exactly sure. But this source seems to be explaining it well. Take a look at it. google.com/amp/s/physicspython.wordpress.com/2019/07/31/… Dec 2 '20 at 16:14

As you suggesting, your computed $$g(r)$$ goes to zero at large $$r$$ because of your handling of the periodic boundary conditions. You should include every particle at a distance of $$r$$ whether it's in the main cell or not, and not only the closest particle.
If it's not clear for you, the following picture may help you. Imagine you're computing $$g(r)$$ for the red particle of the bottom-left corner, then getting only the closest particle would mean that you only take $$r_1$$ into account, while $$r_2$$, $$r_3$$ and $$r_4$$ are also smaller than $$1$$ and should be considered when computing $$g(r)$$.
• @megamence this is quite interesting. Why do you divide $n(r)$ by $N$ in $g(r)$? Dec 2 '20 at 20:21
• @megamence, can you try and see what you get when you take PBC into account, the way Hebo suggests, and also by removing $N$? I think the $1$ represents the fact that at larger distances the subject will see $4\pi r^2\Delta r\rho$ number of particles. Dec 2 '20 at 23:09
• @megamence the best way to get convinced about whether to remove the $N$ in your normalisation factor may be to look at your limit at large $r$ and check whether it's around $1$