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

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)
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])
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))

# 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)
a_file.close()

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

df.to_csv('struc_file.csv')

df.dropna(inplace = True)

column_label = ["ID", "type", "b", "c"]

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
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

• Try to print result and bins and see if they have numbers in them. It looks like something is going wrong in that loop where you assign g[p,:], because that's where everything is zero. Jul 27 at 17:00
• @SRMaiti I thought that was the problem! d contains value and bins = edges but g[p, :] contains 1D zeros array. I have checked numpy.histogram but that didn't help to spot the problem. Jul 27 at 17:26
• Print out the value of result and bins i.e. the output of numpy.histogram. Are they both zeros ? Also, the numpy documentation says normed=False is deprecated and gives wrong results. You may want to try the density argument instead. 9https://numpy.org/doc/stable/reference/generated/numpy.histogram.html) Jul 28 at 5:28
• @SRMaiti result is zeros, bins is not which is the same as edges, and yes I tried density=False but that didn't work. Jul 28 at 10:01

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.