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

Question

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?

Answer 1

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