How to calculate the electron density for a specific band

[jerry]

Dear,

How to calculate the electron density of the HOMO and LUMO band and a certain spin direction? i.e. the spin up electron density of the HOMO band of a graphene nanoribbon.

Thanks a lot!

[zh]

Currently, there is no straightforward way to do it. One may call "BlochState()" to calculate the Bloch state of a specified eigenvalue, and then estimate the modul of such state.

[jerry]

Thank you very much!
Is there any small programs using python in atk to calculate the electron density?

[zh]

Not available. You have to do such programming by yourself.

[Anders Blom]

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:

Code: python

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.

Code: python

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, 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!

[kaypu]

Dear professor Anders, i got a problem with the script.
***************************************
from NanoLanguage import *
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)

bulk_configuration = nlread("D:/DATA/ATK/band/zig8-1.nc", BulkConfiguration)[0]
conf = nlread("D:/DATA/ATK/band/zig8-1.nc", BulkConfiguration, read_state=False)[0]

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_z8_HOMO_up.nc"

# Number of points to sample from G to Z
Nk = 20

density = None
for kz in numpy.linspace(0.42,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

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 wrong message shows as follows:

Traceback (most recent call last):
  File "c:\users\kaypu\appdata\local\temp\0851755004935871.py", line 50, in <module>
    grid = makeGridValues(density,cell,origin,unit = eV/eV)
  File "c:\users\kaypu\appdata\local\temp\0851755004935871.py", line 20, in makeGridValues
    return GridValues(grid3d,unit)
NameError: global name 'GridValues' is not defined

what's wrong with it? i use atk 11.8.2

[Anders Blom]

12.8 is required for this to work.

[kaypu]

thank you  professor