# Setup special points in the Brillouin zone
bcc_symmetry_points = {
   "G"  : (0.00, 0.00, 0.00),
   "H"  : (0.50,-0.50, 0.50),
   "P"  : (0.25, 0.25, 0.25),
   "N"  : (0.00, 0.00, 0.50),
   "N'" : (0.50, 0.50, 0.00),
   "H'" : (0.50, 0.50, 0.50)
  }

# Specify route in the Brillouin zone
bcc_route = ["G","H","N","G","P","N'","P","H'"]

# Call auxillary function for plotting spin-polarized band structure
from bandstructure import calculateAndPlotSpinBandstructure
calculateAndPlotSpinBandstructure(
    ncfilename = nc_file,
    datafilename = "dft_Febcc.dat",
    plotfilename = "plotBandStructure.py",
    imagename = 'Fe_bcc.png',
    symmetry_points = bcc_symmetry_points,
    route = bcc_route,
    number_of_kpoints = 25

The two figures below show the band structure of bcc iron for a spin-polarized (Figure 22) and unpolarized (Figure 23) calculation. We see that the spin-up states (red lines) in general are located lower in energy than the spin-down states. This explains why more spin-up states than spin-down states are occupied.

In this example, we have used matplotlib to generate our band structure plot. The auto-generated plotting script is shown below. This script will fail to run if you do not have the right Python tools installed. Please see the homepage of matplotlib for specific details as to what should be installed on your particular platform if you want to use the script below.

import matplotlib
matplotlib.use('Agg')
from pylab import *

import pickle
f = open("dft_Febcc.dat","rb")
X    = pickle.load(f)
Y_up = pickle.load(f)
Y_dn = pickle.load(f)
f.close()

for i in range(len(Y_up)):
    plot(X,Y_up[i],'r-')
    plot(X,Y_dn[i],'b-')
  
ylabel('Energy (eV)')
xticks ([0.0, 0.8660254037844386, 1.5731321849709863, 2.0731321849709863, 
  2.5061448868632055, 2.9391575887554247, 3.3721702906476438, 3.805182992539863],
  ['G', 'H', 'N', 'G', 'P', "N'", 'P', "H'"])
grid(True)
v1 = axis()
v2 = [v1[0],3.80899,v1[2],v1[3]]
axis(v2)

This module imports other modules, which are either built into NanoLanguage or described below. The following sections describe each function contained in the bandstructure.py module. This module is generalized and can be used, in combination with the other modules listed below, to calculate the band structure of a bulk system along user specified directions in the Brillouin zone. For more details, see the code comments in the scripts.

from ATK.KohnSham import *
from numpy import array
from separations import *

def writeBandstructureDataFile (filename,symmetry_points,route,
                                k_points,bandstructure_up,bandstructure_dn):
    """
    Create a Pickle data file with the bandstructure
    Three objects are dumped: kx (1D values of the k-points
    for the horizontal axis), and energies for up and down spin
    # The last symmetry point is of course not included
    # except as endpoint for the last segment;
    # hence the -1 in the range commands
    """
    # Generate 1D k-points for the horizontal axis
    kx_points = []
    separations,separationsums = generateSeparations(symmetry_points,route)
    for point in range(len(route)-1):
        numKPoints = len(k_points[point])
        kx_points += convertVectorPointToX(numKPoints,point,separations,separationsums)
    # Make empty lists of the right size for the bandstructure data
    # Each list has as many elements as there are bands
    # Each element is a list, the length of which matches the total number of k-points
    energies_up = []
    energies_dn = []
    number_of_bands = bandstructure_up[0].numberOfBands()
    for band in range(number_of_bands):
        energies_up.append([])
        energies_dn.append([])
    # Extract the bands, strip off the unit, and store in the lists
    for point in range(len(route)-1):
        for band in range(number_of_bands):
            tmp = bandstructure_up[point].band(band)
            energies_up[band] = energies_up[band] + [E.inUnitsOf(eV) for E in tmp]
            tmp = bandstructure_dn[point].band(band)
            energies_dn[band] = energies_dn[band] + [E.inUnitsOf(eV) for E in tmp]
    # Store data in pickle file
    import pickle
    f = open(filename,"wb")
    pickle.dump(kx_points,f)
    pickle.dump(energies_up,f)
    pickle.dump(energies_dn,f)
    f.close()

def calculateBandstructure(calculation,symmetry_points,route,n_points,spin):
    """
    Calculate the bandstructure for the "calculation" object,
    using the provided symmetry points and the route that
    connects them. This routines defines the k-points,
    calculate the bandstructure, and return both. In principle
    each segment of the bandstructure plot can contain a
    different number of k-points, but for simplicity we choose
    an equal number for now.
    """
    band_k_points = []
    bandstructure = []
    # The last symmetry point is of course not included
    # except as endpoint for the last segment
    for point in range(len(route)-1):
        start_point = symmetry_points[route[point]]
        end_point = symmetry_points[route[point+1]]
        band_k_points.append(generateLineBetweenVectorsPoints(start_point,end_point,n_points))
        bandstructure.append(calculateEnergyBands(calculation,band_k_points[point],spin))
    return band_k_points, bandstructure

def writePyLabScript(plotfilename,datafilename,imagename,symmetry_points,route):
    # Create pyLab script to plot the data
    f = open(plotfilename,"w")
    separations,separationsums = generateSeparations(symmetry_points,route)
    # The text block to be written to the script
    s = """
import matplotlib
matplotlib.use('Agg')
from pylab import *

import pickle
f = open("%s","rb")
X    = pickle.load(f)
Y_up = pickle.load(f)
Y_dn = pickle.load(f)
f.close()

for i in range(len(Y_up)):
    plot(X,Y_up[i],'r-')
    plot(X,Y_dn[i],'b-')
  
ylabel('Energy (eV)')
xticks (%s,%s)
grid(True)
v1 = axis()
v2 = [v1[0],%g,v1[2],v1[3]]
axis(v2)
savefig('%s',dpi=100)
    """ % (datafilename,separationsums,\
    route,separationsums[-1]*1.001,imagename)
    f.write(s)

def calculateAndPlotSpinBandstructure(
        ncfilename,datafilename,plotfilename,imagename,
        symmetry_points,route,number_of_kpoints=25):

    calculation = restoreSelfConsistentCalculation(ncfilename)
    band_k_points,bandstructure_up = calculateBandstructure(
        calculation, symmetry_points, route,
        number_of_kpoints, Spin.Up)

    band_k_points,bandstructure_dn = calculateBandstructure(
        calculation, symmetry_points, route,
        number_of_kpoints, Spin.Down)

    writeBandstructureDataFile(datafilename,symmetry_points,route,
                               band_k_points,bandstructure_up,bandstructure_dn)

    writePyLabScript(plotfilename,datafilename,imagename,symmetry_points,route)

bandstructure.py

The functions contained in this module, could just as well have been included in the bandstructure.py module. However, they could be usable in other situations, where you might not be interested in the energy bands. Therefore, they have been placed in a separate module. It is generally a good idea to modularize your scripts once they grow in size and complexity. Furthermore, this allows you to reuse methods, thus saving yourself the time of rewriting a simple library function every time.

This function helps us to project a continuous path between a sequence of points in 3D space onto a one-dimensional axis. In our case the points will be k-points in the Brillouin zone, but the script is general, and makes no assumption regarding the nature of the points (only, that they are dimensionless).

The separation between two points on the one-dimensional axis will be proportional to the vector distance between two sequential points. We calculate this using the dot function from NumPy on the vector delta (a NumPy array) between two points on the route. These distances are stored in the variable separations, which is a simple list; the elements in this list will be the one-dimensional values corresponding to the points on the path. We also store the accumulated sum of these values in separationsums, for later use (see the function convertVectorPointToX() below).

This function generates the corresponding one-dimensional values for N points, evenly distributed between two points in space. It is necessary to have generated the separations using the function generateSeparations() (see above).

Given two points in space, this function generates a sequence of N evenly distributed points between them, and returns these in a list.

from numpy import array, dot
import math 

def generateSeparations (symmetry_points, route):
    """
    Given a dict of 3D vector points and the route which
    connects them, generate the so-called separations and
    separation sums.
    """
    separations = []
    # The first point is always placed at x=0
    separationsums = [0.0]
    # The last symmetry point is of course not included
    # except as endpoint for the last segment
    for p in range(len(route)-1):
        delta = array(symmetry_points[route[p+1]]) - array(symmetry_points[route[p]])
        # Calculate and store the vector distance between subsequent points
        separations.append(math.sqrt(dot(delta,delta)))
        # Cumulative sum of the distances
        separationsums.append(separationsums[-1] + separations[-1])
    return separations, separationsums

def convertVectorPointToX (N, id, separations, separationsums):
    """
    Use the separations and separation sums to map 3D vector
    points onto a one-dimensional dummy coordinate for the
    horizontal axis, e.g. to be used in a bandstructure plot.
    id is the segment number.
    """
    x = []
    for i in range(N):
        # If changes are made to the definition of the vector points
        # in generateLineBetweenVectorsPoints below,
        # the denominator here must be changed accordingly
        x.append(float(i)/(N-1)*separations[id] + separationsums[id])
    return x

def generateLineBetweenVectorsPoints (start_point,end_point,N):
    """
    Given two points (given as lists), generate a sequence of
    N points between them, evenly distributed and return in a list.
    """
    b = []
    for i in range(N):
        # The end-point is included; to exclude it, use N in the denominator instead
        # If this is changed, the change should be mirrored in convertVectorPointToX above
        a = (float(i)/(N-1),)*3*(array(end_point)-array(start_point)) + start_point
        b.append(a)
    return b

separations.py