thomson.py 文件源码

def xyz2U(x, y, z):
    """
    compute the potential
    """
    dsq = 0.

    indices = np.triu_indices(x.size, k=1)

    for coord in [x, y, z]:
        coord_i = np.tile(coord, (coord.size, 1))
        coord_j = coord_i.T
        dsq += (coord_i[indices]-coord_j[indices])**2

    d = np.sqrt(dsq)
    U = np.sum(1./d)
    return U