Extending Radial Distributions

Radial distribution functions are one of the most important quantities to calculate when running a molecular dynamics simulation. They give good insight into the phase and pair potentials between molecules. Practically, they are calculated according to the following equation:

where v_i is a small volume increment that corresponds to the radius increment r_i, n_i is the average number of particles found in the volume v_i, and rho_bulk is the bulk density of the entire simulation. Essentially, at each incremental volume we're calculating the ratio of particles found to what we would find if the particles in that volume behaved like the bulk density. The incremental volume in which we search for particles may be seen in the image below:

However, we notice there is a problem as r_i grows: the concentric circles (or spheres in 3D) begin to intersect with the simulation box that contains the particles. This is is shown below.

Generally this problem is approached by stopping a radial distribution calculation when r_i is equal to half the side length of the simulation box. That actually means we lose a significant amount of statistics, pi / 8 in 3 dimensions to be exact. This in fact may fixed by using the results of my previous post about cube-sphere intersections. In fact, I wrote it specifically about this issue.

The way to allow radial distribution functions to be calculated all the way out until all particles are counted is to calculate the true v_i, which is the intersection volume between the simulation box and the concentric spheres that make up v_i . Using the code I wrote, v_i may be calculated according to:

The relevant radial distribution function calculation now looks like this:
Finally, here are the results on Lennard-Jones fluid calculation:

The simulation box ends at the end of the line. As you can see, there is a huge amount of space lost in traditional radial distribution function calculations. Another cool thing about correcting the volume is that the radial distribution does not have to be normalized to converge to 1. Often in these calculations people struggle with double counting and incorrect volumes so that the value is just forced to converge to 1 at the edge of the plot. That is the value which it should converge to. Here's the full radial distribution calculator, which is part of a simulation engine I'm writing.