7
$\begingroup$

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 1-by-1-by-1 box with periodic boundary conditions and with $N = 500$ particles of diameter $d\approx0.09$ (approx) 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 dont 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 https://demonstrations.wolfram.com/RadialDistributionFunctionForHardSpheres/, 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?

$\endgroup$
3
  • 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$ – Kavya Mrudula Dec 2 '20 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 '20 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$ – Kavya Mrudula Dec 2 '20 at 16:14
7
$\begingroup$

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

pbc_and_rdf

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