# Import an electron density
density = nlread('results.hdf5', ElectronDensity)[0]

# Get the shape of the data.
shape = density.shape()

# Find the volume elements.
dX, dY, dZ = density.volumeElement()

# Calculate the unit area in the y-z plane.
dAYZ = numpy.cross(dY,dZ)[0]

# calculate density along x integrated over y,z
n_x = [ density[i,:,:].sum() * dAYZ for i in range(shape[0]) ]

# Integrate the density along x to get the number of electrons.
dx = dX.norm()
number_of_electrons = 0

for i in range(shape[0]):
    number_of_electrons += n_x[i]*dx

print('Total number of electrons = {:.5f}'.format(number_of_electrons))

 # Alternative to calculating the integral
dX, dY, dZ = density.volumeElement()

# Volume of a single grid point
dV = numpy.dot(numpy.cross(dX, dY), dZ)

density_unit = density.unit()
density_array = density.toArray() * density_unit
integral = density_array.sum() * dV

print('Alternative method:')
print('Total number of electrons= {:.5f}'.format(integral))