Have I missed something? RDF (radial distribution) script cannot capture the correct .xyz

Question

I am trying to learn how RDF is working. To learn how RDF is working in practice I am testing the script (below) by printing all of the variables using this .xyz file example:

I passed cartesian coordination data correctly, and the boundary condition and the printed variables are correct which I checked with manual calculations.

However, the resulting plot only shows a flat line at y = 0. Probably I missed something small thing but I cannot catch it. Have I missed something here?

      16
Energy           -200.0
A            1.977502779         1.825612486        -1.078815994
A            0.073484389        -2.915354734        -1.169129839
A           -1.682844787        -1.543503043        -2.245494959
A            0.226989000         0.103121000        -0.417822000
B            -0.693058883         0.156772052         1.151824239
B            -1.448474661        -3.517890885        -2.298992143
B            -1.458396055        -1.997135497        -0.344566446
B             0.985126104        -4.427396897        -0.775735938
B            -3.121800014        -1.219516661        -3.292662828
B             2.017385825         0.679529254         0.430803534
B             0.212637914        -1.726148783        -2.725852021
B             3.217540502         3.061671270        -1.526834132
B             1.456477430         0.338098844        -2.124519369
B            -1.187423538         0.191670365        -1.675742064
B             0.143729055         2.143022931        -0.717464213
B             1.083320805        -1.581985916        -0.284118283

the printed variables are

edges: [-12.    -11.999 -11.998 ...   1.199   1.2     1.201]
num_increments: 13201
x is [ 0.648635 -0.712093  1.175089  2.536351  0.431639  0.68227   1.776624
2.48202  -2.049614  2.48957  -1.061626  0.443473  0.890669 -0.407062
1.391732  3.873991]
len(x) is 16
S is 12
numberDensity is 0.009259259259259259
d is [4.16134966 5.64980835 4.16167321 1.64793118 5.16051574 4.84601737
2.79507052 3.09205521 7.29774691 2.7948553  5.45595986 5.45590615
3.0918364  4.84603638 5.16043525 2.4       ]
g[p, :]: [0. 0. 0. ... 0. 0. 0.]
result is [0 0 0 ... 0 0 0]
numberDensity is 0.009259259259259259

def pairCorrelationFunction_3D(x, y, z, S, rMax, dr):
"""Compute the three-dimensional pair correlation function for a set of
spherical particles contained in a cube with side length S.  This simple
function finds reference particles such that a sphere of radius rMax drawn
around the particle will fit entirely within the cube, eliminating the need
to compensate for edge effects.  If no such particles exist, an error is
returned.  Try a smaller rMax...or write some code to handle edge effects! ;)
Arguments:
    x               an array of x positions of centers of particles
    y               an array of y positions of centers of particles
    z               an array of z positions of centers of particles
    S               length of each side of the cube in space
    rMax            outer diameter of largest spherical shell
    dr              increment for increasing radius of spherical shell
Returns a tuple: (g, radii, interior_indices)
    g(r)            a numpy array containing the correlation function g(r)
    radii           a numpy array containing the radii of the
                    spherical shells used to compute g(r)
    reference_indices   indices of reference particles
"""
from numpy import zeros, sqrt, where, pi, mean, arange, histogram

# Find particles which are close enough to the cube center that a sphere of radius
# rMax will not cross any face of the cube
bools1 = x > rMax
bools2 = x < (S - rMax)
bools3 = y > rMax
bools4 = y < (S - rMax)
bools5 = z > rMax
bools6 = z < (S - rMax)

interior_indices, = where(bools1 * bools2 * bools3 * bools4 * bools5 * bools6)
num_interior_particles = len(interior_indices)

if num_interior_particles < 1:
    raise  RuntimeError ("No particles found for which a sphere of radius rMax\
            will lie entirely within a cube of side length S.  Decrease rMax\
            or increase the size of the cube.")

edges = arange(-S, rMax + 1.1 * dr, dr)
num_increments = len(edges) - 1
g = zeros([num_interior_particles, num_increments])
radii = zeros(num_increments)
numberDensity = len(x) / S**3

# Compute pairwise correlation for each interior particle
for p in range(num_interior_particles):
    index = interior_indices[p]
    d = sqrt((x[index] - x)**2 + (y[index] - y)**2 + (z[index] - z)**2)
    d[index] = 2 * rMax

    (result, bins) = histogram(d, bins=edges, normed=False)
    g[p,:] = result / numberDensity

# Average g(r) for all interior particles and compute radii
g_average = zeros(num_increments)
for i in range(num_increments):
    radii[i] = (edges[i] + edges[i+1]) / 2.
    rOuter = edges[i + 1]
    rInner = edges[i]
    g_average[i] = mean(g[:, i]) / (4.0 / 3.0 * pi * (rOuter**3 - rInner**3))

return (g_average, radii, interior_indices)
# Number of particles in shell/total number of particles/volume of shell/number density
# shell volume = 4/3*pi(r_outer**3-r_inner**3)



# preprocess the structure file (struc)
a_file = open(struc)
lines = a_file.readlines()
a_file.close()

# del first two lines
del lines[0]
del lines[0]

df = pd.read_fwf(struc)
df.to_csv('struc_file.csv')

df.dropna(inplace = True)

column_label = ["ID", "type", "b", "c"]
df = pd.read_csv('struc_file.csv', names=column_label)

df = df.drop([0, 1])    # first and second row
df = df.drop(columns = ["ID"])
new = df["b"].str.split(" ", n = 1, expand = True)

df["x"] = new[0]
df["y"] = new[1]
df["z"] = df["c"]
df = df.drop(columns = ["b", "c"])
df = df.reset_index(drop=True)


# Calculation setup
domain_size = 12
num_particles = 10

dr = 0.001
particle_radius = 0.1
rMax = domain_size / 10

g_r, r, reference_indeces = pairCorrelationFunction_3D(x_particle, y_particle, z_particle, domain_size, rMax, dr)

plt.figure()
plt.plot(r, g_r, color='black')
plt.xlabel('r')
plt.ylabel('g(r)')
plt.xlim( (-rMax, rMax) )
plt.ylim( (0, 1.05 * g_r.max()) )
plt.show()
#The script is from https://github.com/cfinch/Shocksolution_Examples/blob/master/PairCorrelation/example_3D.py

Answer 1

Ok, I think I figured out why your code is giving a straight line. This is too big to write as a comment, so I am putting this as answer.

At the start of the function, you are determining the edges of the bins of histograms:

edges = arange(-S, rMax + 1.1 * dr, dr)

So, your edges range between -S and rMax+1.1*dr. From the data you fed to the function, the edges range from -12 to 1.2 (edges: [-12. -11.999 -11.998 ... 1.199 1.2 1.201]).

Your edges are mainly in the negative region, with only a small part at the end of the range extending into the positive region.

Now, the array d holds the distances between atom pairs, all of which are obviously positive numbers. And the values of d you show:

d is [4.16134966 5.64980835 4.16167321 1.64793118 5.16051574 4.84601737
2.79507052 3.09205521 7.29774691 2.7948553  5.45595986 5.45590615
3.0918364  4.84603638 5.16043525 2.4       ]

all seem to be above 1.201. The function is calculating histogram in the range you gave, and since there are no values in that range, your histogram is $0$ all over.

It seems that there is some error in the formula that you used to calculate edges. I can't say what the correct formula should be without studying the code in more detail. As it's your code, you would be able to fix it much faster.

Personally, I would suggest not calculating edges explicitly, but rather using the range=() argument and just specifying the number of bins in numpy.histogram.