I am trying to calculate the pair correlation function from a given FDF file. It computes the three-dimensional pair correlation function for a set of spherical particles in the LatticeVectors. 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.

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

My code is:

def pair_correlation_function(label=None, path="i*", dr=0.1, plot_=True):
    fdfpath = glob.glob(f"{path}{os.sep}{label}.fdf")[0]
    x, y, z = coords(fdfpath)
    x = np.array(x)
    y = np.array(y)
    z = np.array(z)
    sx, sy, sz = lattice_vectors_mag(fdfpath)
    _, _, species_ = species(fdfpath)
    rmax = max(element_diameter(specie) for specie in species_)
    bools1 = x > rmax
    bools2 = x < sx - rmax
    bools3 = y > rmax
    bools4 = y < sy - rmax
    bools5 = z > rmax
    bools6 = z < sz - rmax
    (interior_indices,) = np.where(bools1 * bools2 * bools3 * bools4 * bools5 * bools6)
    num_interior_particles = len(interior_indices)
    if num_interior_particles < 1:
        raise RuntimeError("No particles found. Increase the LatticeVectors.")
    edges = np.arange(0.0, rmax   1.1 * dr, dr)
    num_increments = len(edges) - 1
    g = np.zeros([num_interior_particles, num_increments])
    radii = np.zeros(num_increments)
    numberdensity = len(x) / sx * sy * sz
    for p in range(num_interior_particles):
        index = interior_indices[p]
        d = np.sqrt((x[index] - x)**2   (y[index] - y)**2   (z[index] - z)**2)
        d[index] = 2 * rmax
        result, bins = np.histogram(d, bins=edges, density=False)
        g[p, :] = result / numberdensity
    g_average = np.zeros(num_increments)
    for i in range(num_increments):
        radii[i] = (edges[i]   edges[i   1]) / 2.0
        router = edges[i   1]
        rinner = edges[i]
        g_average[i] = np.mean(g[:, i]) / (4.0 / 3.0 * np.pi * (router**3 - rinner**3))
    return g_average, radii, interior_indices

I am getting TypeError: '>' not supported between instances of 'numpy.ndarray' and 'float' error for the line bools1 = x > rmax.

CodePudding user response：

You have to test the condition for every element of the list. Try smtg like:

a = [1,0,3,5,6,0]

b = [i>0 for i in a]

Output:

[True, False, True, True, True, False]

CodePudding user response：

Check the elements in numpy.ndarray,maybe they are 'str' or other types instead of 'int' or 'float'.