thomson.py 文件源码

def elec_potential_xyz(x0):
    x0 = x0.reshape(3, x0.size/3)
    x = x0[0, :]
    y = x0[1, :]
    z = x0[2, :]
    dsq = 0.

    r = np.sqrt(x**2 + y**2 + z**2)
    x = x/r
    y = y/r
    z = z/r
    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
        d = (coord_i[indices]-coord_j[indices])**2
        dsq += d

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