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
评论列表
文章目录