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 enter image description here

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?

  • 2
    $\begingroup$ 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. $\endgroup$ Commented Dec 2, 2020 at 6:38
  • $\begingroup$ @KavyaMrudula, thanks for your answer. What do you mean by dividing by the number of shells? you mean the number of bins I have? $\endgroup$
    – megamence
    Commented Dec 2, 2020 at 15:06
  • $\begingroup$ 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/… $\endgroup$ Commented Dec 2, 2020 at 16:14

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


  • $\begingroup$ thanks for your answer! Do you have any suggestions on how to normalize g(r)? $\endgroup$
    – megamence
    Commented Dec 2, 2020 at 15:49
  • $\begingroup$ @megamence this is quite interesting. Why do you divide $n(r)$ by $N$ in $g(r)$? $\endgroup$ Commented Dec 2, 2020 at 20:21
  • 1
    $\begingroup$ to find the average number of particles in each bin @HitanshuSachania $\endgroup$
    – megamence
    Commented Dec 2, 2020 at 21:29
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Dec 2, 2020 at 23:09
  • 1
    $\begingroup$ @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$ $\endgroup$
    – Hebo
    Commented Dec 3, 2020 at 9:02

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .