One the one hand this is fairly straightforward Python programming. At the same time there are some challenging moments.
Assuming you have performed a spin-polarized calculation of a nanoribbon, the following rather self-explanatory lines will sum up the electron density, defined as the Bloch function squared, of the HOMO (or LUMO) band from Gamma to Z:
bulk_configuration = nlread("gnr_z14.nc", BulkConfiguration)[0]
conf = nlread("gnr_z14.nc", BulkConfiguration, read_state=False)[0]
# Some magic tricks to automatically determine the LUMO/HOMO
dm = bulk_configuration.calculator()._densityMatrixCalculator()
number_of_electrons = int(dm.fermiDistribution().numberOfElectrons())
LUMO = number_of_electrons/2
HOMO = LUMO-1
# Which state to compute?
spin = Spin.Up
state = HOMO
filename = "bloch_z14_HOMO_up.nc"
# Number of points to sample from G to Z
Nk = 100
# Here we go!
density = None
for kz in numpy.linspace(0,0.5,Nk):
b = BlochState(bulk_configuration, quantum_number=state, k_point=[0,0,kz], spin=spin)
psi = b.toArray()
if density == None:
density = (psi*psi.conj()).real
else:
density += (psi*psi.conj()).real
Now we have the density - but how to visualize it? This gets a bit trickier, but here's the code for it.
import NLEngine
def makeGridValues(datagrid, cell, origin, unit):
n0, n1, n2 = datagrid.shape
gA, gB, gC = cell
u0 = NLEngine.Cartesian3D(gA[0],gA[1],gA[2])
u1 = NLEngine.Cartesian3D(gB[0],gB[1],gB[2])
u2 = NLEngine.Cartesian3D(gC[0],gC[1],gC[2])
cell_origin = NLEngine.Cartesian3D(origin[0],origin[1],origin[2])
grid_descriptor = NLEngine.GridDescriptor(n0,n1,n2,
NLEngine.UnitCell(u0,u1,u2,cell_origin))
grid3d = NLEngine.RealGrid3D(grid_descriptor,
NLEngine.doubleSequenceToRealVector(datagrid.flatten()),True)
return GridValues(grid3d,unit)
cell = conf.bravaisLattice().primitiveVectors().inUnitsOf(Bohr)
origin = conf.bravaisLattice().origin().inUnitsOf(Bohr)
grid = makeGridValues(density,cell,origin,unit = eV/eV)
nlsave(filename, grid)
nlsave(filename, conf)
The order of the blocks got a bit messed up here - you can use the attached script instead. The picture for a zigzag ribbon of width=14 for HOMO/down is attached as example, but actually the picture is virtually identical for up/down (graphene is after all not magnetic) - but also for HOMO/LUMO!
(For a bit better picture it might be a good idea to do the calculation for a supercell repeated 2x in Z.)
The results here can be compared to the tutorial on Bloch states (http://quantumwise.com/documents/tutorials/latest/GrapheneBloch), where we see that at the Z point the Bloch states are separated on two sides of the ribbon - but when we sum up the entire band, this is not the case!