计算三维粒子对关联函数时遇numpy.ndarray与float比较错误

尝试从给定FDF文件计算三维球形粒子的对关联函数，该函数筛选完全处于LatticeVectors构成立方体内的参考粒子以消除边缘效应，返回元组(g, radii, interior_indices)。运行代码时，在bools1 = x > rmax行触发错误：

TypeError: '>' not supported between instances of 'numpy.ndarray' and 'float'

import numpy as np
import glob
import os

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

核心错误原因 #

报错的直接原因是**x数组的元素类型不是数值型**：比如coords(fdfpath)返回的x是包含非数值元素的列表，转成numpy数组后类型为object，导致无法和float类型的rmax进行比较运算。

修复步骤 #

1. 强制指定数组数值类型
   将数组转换代码修改为：

   x = np.array(x, dtype=np.float64)
   y = np.array(y, dtype=np.float64)
   z = np.array(z, dtype=np.float64)
   

   确保数组元素为浮点型，支持和float的比较操作。

2. 验证rmax的有效性
   添加代码检查rmax的类型和值，确保它是合法的浮点数值：

   rmax = max(element_diameter(specie) for specie in species_)
   assert isinstance(rmax, (int, float)), f"rmax must be a number, got {type(rmax)}"
   

   同时检查element_diameter函数的返回值，确保它返回的是数值类型而非字符串或其他非数值类型。

3. 优化内部粒子筛选逻辑
   用逻辑与运算符&代替布尔数组乘法，代码更清晰且效率更高：

   mask = (x > rmax) & (x < sx - rmax) & (y > rmax) & (y < sy - rmax) & (z > rmax) & (z < sz - rmax)
   interior_indices = np.where(mask)[0]
   

4. 修复数密度计算错误
   原代码中numberdensity = len(x) / sx * sy * sz因运算符优先级问题，实际计算的是(len(x)/sx)*sy*sz，正确的数密度应该是粒子数除以体积：

   numberdensity = len(x) / (sx * sy * sz)
   

import numpy as np
import glob
import os

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, dtype=np.float64)
    y = np.array(y, dtype=np.float64)
    z = np.array(z, dtype=np.float64)
    
    sx, sy, sz = lattice_vectors_mag(fdfpath)
    _, _, species_ = species(fdfpath)
    
    rmax = max(element_diameter(specie) for specie in species_)
    # 验证rmax是数值类型
    assert isinstance(rmax, (int, float)), f"rmax must be a number, got {type(rmax)}"
    
    # 使用逻辑与构建筛选掩码
    mask = (x > rmax) & (x < sx - rmax) & \
           (y > rmax) & (y < sy - rmax) & \
           (z > rmax) & (z < sz - rmax)
    interior_indices = np.where(mask)[0]
    
    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

内容的提问来源于stack exchange，提问作者Eftal Gezer