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
result
andbins
and see if they have numbers in them. It looks like something is going wrong in that loop where you assigng[p,:]
, because that's where everything is zero. $\endgroup$d
contains value andbins = edges
butg[p, :]
contains 1D zeros array. I have checkednumpy.histogram
but that didn't help to spot the problem. $\endgroup$result
andbins
i.e. the output of numpy.histogram. Are they both zeros ? Also, the numpy documentation saysnormed=False
is deprecated and gives wrong results. You may want to try thedensity
argument instead. 9https://numpy.org/doc/stable/reference/generated/numpy.histogram.html) $\endgroup$result
is zeros,bins
is not which is the same asedges
, and yes I trieddensity=False
but that didn't work. $\endgroup$