我正在努力学习rdf是如何工作的。为了了解rdf在实践中是如何工作的，我正在通过使用以下.xyz文件示例打印所有变量来测试脚本（如下所示）：
我正确地传递了笛卡尔坐标数据，边界条件和打印的变量都是正确的，我通过手动计算进行了检查。
然而，结果图仅显示y=0处的一条平线。可能我错过了一些小东西，但我无法捕捉到它。我错过什么了吗？

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

打印的变量是

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