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$ 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
    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$ 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
    Dec 2, 2020 at 15:49
  • $\begingroup$ @megamence this is quite interesting. Why do you divide $n(r)$ by $N$ in $g(r)$? $\endgroup$ Dec 2, 2020 at 20:21
  • 1
    $\begingroup$ to find the average number of particles in each bin @HitanshuSachania $\endgroup$
    – megamence
    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$ 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
    Dec 3, 2020 at 9:02

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.