# Calculate the current density at the fermi level from right to left going states
cd = CurrentDensity(
    configuration,
    energies = [0.0]*eV,
    energy_weights = Right,
    kpoints=MonkhorstPackGrid(1,11),
    )

# Calculate the volume element and x,y integration area
dX, dY, dZ = cd.volumeElement()
dAXY = numpy.linalg.norm( numpy.cross(dX,dY) )

# Calculate the current density along z integrated over x,y
shape = cd.shape()
# Take out spin up (index 0) z component (index 2) in every grid point.
cd_z = numpy.array( [ [ [ cd[i,j,k][0][2] for k in range(shape[2])] for j in range(shape[1])] for i in range(shape[0]) ] )

# Sum across x and y dimensions
cd_z_integrated = numpy.array([cd_z[:,:,k].sum() * dAXY for k in range(shape[2])])


# Plot the data
import pylab
pylab.figure()
dz = dZ.norm()
z = numpy.arange(shape[2])*dz.inUnitsOf(Bohr)
pylab.plot(z,cd_z_integrated)
pylab.show()