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:

v[i] = calc_volume(r[i + 1], d) - calc_volume(r[i], d)
view raw gistfile1.py hosted with ❤ by GitHub
The relevant radial distribution function calculation now looks like this:
for i in range(len(bin_edges) - 1) :
v = calc_volume(bin_edges[i+1], box_size[0] / 2.) - calc_volume(bin_edges[i], box_size[0] / 2.)
ideal = v * bulk_density
actual = hist[i]
r[i] = (bin_edges[i + 1] + bin_edges[i]) / 2.
n[i] = actual
gofr[i] = actual / ideal
view raw gistfile1.py hosted with ❤ by GitHub
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.

#!/usr/bin/python
import numpy as np
from scipy.integrate import quad
from math import pi
def read_run_params(fname):
params = {}
with open(fname, 'r') as f:
for line in f.readlines():
sline = line.split()
if(len(sline) == 2):
params[sline[0]] = sline[1]
return params
def pair_distance(x, y, ndims, img, PBC):
dsum = 0
dist = 0
for j in range(ndims):
dist = (x[j] - y[j])
if(PBC):
dist = dist - round(dist / img[j]) * img[j]
dsum += dist ** 2
return np.sqrt(dsum)
def bin_atom(i, distances, positions, bin_edges, ndims, img, PBC):
for j in range(i+1, len(positions)):
distances[j] = pair_distance(positions[i], positions[j], ndims, img, PBC)
return np.histogram(distances[i+1:], bin_edges)
#square regions
def region_1(rho, d):
return 0.5 * (d**2 * np.sqrt(rho**2 - 2 * d ** 2))
def region_1_2_theta(rho, d):
return np.arccos(d / np.sqrt(rho ** 2 - d ** 2))
#curved square regions
def region_2_integrand(theta, rho, d):
return np.sqrt(np.cos(theta) ** 2 - (d / rho)**2) / (np.cos(theta) ** 3)
def region_2(rho, d):
i4 = d**3 / 6. * (rho**2 / d **2 - 1) * (pi / 4 - region_1_2_theta(rho, d))
i3 = d ** 2 * rho / 3. * quad(region_2_integrand, region_1_2_theta(rho,d), pi / 4, args=(rho, d))[0]
return i4 + i3
#spherical region
def region_3_integrand(theta, rho, d):
return np.sqrt(np.cos(theta) ** 2 - (d / rho)**2) / np.cos(theta)
def region_3(rho, d):
return rho ** 3 / 3. * (d / rho * (pi / 4 - region_1_2_theta(rho, d)) - quad(region_3_integrand, region_1_2_theta(rho, d), pi / 4, args=(rho, d))[0])
def calc_volume(rho, d):
alpha = rho / d
if(alpha <= 1):
return 4./3 * pi * (rho) ** 3
if(alpha <= np.sqrt(2)):
return 4. / 3 * pi * (rho) ** 3 - 6. * pi * (2 * (rho)**3 - 3 * d * (rho)**2 + d**3) / 3.
if(alpha < np.sqrt(3)):
return 16. * (region_1(rho,d) + region_2(rho, d) + region_3(rho,d))
return 8. * d ** 3
def calc_gofr(xyz_file, box_size, bin_width=0.1, ndims=3, PBC=True):
if(ndims != 3):
raise Exception("Have not implemented n-dimensional g(r)")
#if we don't have a cube, I don't yet have the g(r) cube-sphere intersection implemented
if(not(box_size[0] == box_size[1] and box_size[1] == box_size[2])):
bin_edges = np.arange(0, min(box_size) / 2, bin_width)
else:
bin_edges = np.arange(0, np.sqrt(3) * box_size[0] / 2, bin_width)
(hist, bin_edges) = np.histogram([-1], bin_edges)
with open(xyz_file, 'r') as f:
lines_read = 2
frames_read = 0
natoms = 0
try:
natoms = int(f.readline())
except ValueError:
raise Exception("Doesn't look like %s in an xyz file" % xyz_file)
if(natoms == 0):
return None
positions = [0 for x in range(natoms)]
distances = [0 for x in range(natoms)]
print "Processing frame...0",
while(1):
try:
#skip over frames and number of atoms
line = f.readline()
while(line != "" and line.split(":")[0] != "Frame"):
line = f.readline()
for i in range(natoms):
tokens = f.readline().split()
if(len(tokens) != ndims + 1):
raise EOFError()
positions[i] = [float(x) for x in tokens[1:]]
for i in range(natoms):
(newhist, bin_edges) = bin_atom(i,distances, positions, bin_edges, ndims, box_size, PBC)
hist = np.add(hist, newhist)
frames_read += 1
print "\rProcessing frame...%d" % frames_read,
except EOFError:
break
print "\nCalculating g(r)..."
volume = 1.
for x in box_size:
volume = volume * x
bulk_density = sum(hist) / volume
r = np.empty(len(bin_edges) - 1)
n = np.empty(len(bin_edges) - 1)
gofr = np.empty(len(bin_edges) - 1)
for i in range(len(bin_edges) - 1) :
v = calc_volume(bin_edges[i+1], box_size[0] / 2.) - calc_volume(bin_edges[i], box_size[0] / 2.)
ideal = v * bulk_density
actual = hist[i]
r[i] = (bin_edges[i + 1] + bin_edges[i]) / 2.
n[i] = actual
gofr[i] = actual / ideal
return (r,n,gofr)
if __name__ == "__main__":
import sys
print "Note: This program assumes your box is a cube. It does not check"
if(len(sys.argv) != 5):
print "usage: [gofr.py] [input xyz file] [output file] [run_params_file] [bin_width]"
exit()
bin_width = float(sys.argv[4])
params = read_run_params(sys.argv[3])
ndims = int(params['n_dims'])
(r, n, gofr) = calc_gofr(sys.argv[1],
box_size=[int(params['box_%d_size' % x]) for x in range(1,ndims+1)],
bin_width=bin_width,
ndims=ndims)
with open(sys.argv[2], 'w') as f:
for (ri, ni, gofri) in zip(r,n,gofr):
f.write("%g %g %g\n" % (ri, ni, gofri))
view raw gistfile1.py hosted with ❤ by GitHub