# import an electron difference density
difference_density = nlread('results.nc', ElectronDifferenceDensity)[0]

# Get the shape of the data.
ni, nj, nk = difference_density.shape()
# Find the volume elements.
dX, dY, dZ = difference_density.volumeElement()
length_unit = dX.unit()
# Calculate the volume of the volume element.
dV = numpy.dot(dX, numpy.cross(dY,dZ)) * length_unit**3
# Get a reference for the data of the difference density.
grid_data = difference_density[:,:,:]*Units.e

# Calculate the absoute density sum.
density_abs_sum = abs(grid_data).sum()

# Calculate center of mass
Mi = sum([ i*abs(grid_data[i,:,:]).sum() for i in range(ni)])/density_abs_sum
Mj = sum([ j*abs(grid_data[:,j,:]).sum() for j in range(nj)])/density_abs_sum
Mk = sum([ k*abs(grid_data[:,:,k]).sum() for k in range(nk)])/density_abs_sum
center_of_mass = Mi*dX + Mj*dY+ Mk*dZ

# Calculate monopole
density_sum = grid_data.sum()
m0 = density_sum*dV

#calculate dipole
dens_i = sum([ (i-Mi)*grid_data[i,:,:].sum() for i in range(ni)])
dens_j = sum([ (j-Mj)*grid_data[:,j,:].sum() for j in range(nj)])
dens_k = sum([ (k-Mk)*grid_data[:,:,k].sum() for k in range(nk)])
m1 = (dX*dens_i + dY*dens_j + dZ*dens_k)*dV

# Calculate quadropoles (3 x_i *x_j-r*r delta_ij)*dens
dens_ii = sum([ (i-Mi)**2 * grid_data[i,:,:].sum() for i in range(ni)])
dens_jj = sum([ (j-Mj)**2 * grid_data[:,j,:].sum() for j in range(nj)])
dens_kk = sum([ (k-Mk)**2 * grid_data[:,:,k].sum() for k in range(nk)])
dens_ij = sum([ (i-Mi)*(j-Mj)*grid_data[i,j,:].sum()
	          for i in range(ni) for j in range(nj)])
dens_ik = sum([ (i-Mi)*(k-Mk)*grid_data[i,:,k].sum()
	          for i in range(ni) for k in range(nk)])
dens_jk = sum([ (j-Mj)*(k-Mk)*grid_data[:,j,k].sum()
	          for j in range(nj) for k in range(nk)])

m2 = (dX.outer(dX)*dens_ii + dY.outer(dY)*dens_jj + dZ.outer(dZ)*dens_kk\
   + (dX.outer(dY)+dY.outer(dX))*dens_ij + (dX.outer(dZ)+dZ.outer(dX))*dens_ik\
   + (dY.outer(dZ)+dZ.outer(dY))*dens_jk) * dV

r2 = sum(m2)
m2 = 3*m2 - numpy.identity(3)*r2

print "center of mass (bohr)  = ", center_of_mass

print "monopole   (e)         = ", m0
print "dipole     (e*bohr)    = ", m1
print "quadropole (e*bohr**2) = "
print m2