The band structure of bulk silicon

Table of Contents

Introduction

In this example, you will learn how to

  • set up periodic structures, exemplified by a silicon crystal,

  • calculate eigenstates and band structure of silicon

using NanoLanguage. We will also go into some of the numerical details specific for periodic systems, since the understanding of these is important in order to calculate the systems correctly with ATK.

Prerequisites

It will be assumed that you are familiar with the most basic operations in NanoLanguage, such as

  • specifying chemical elements,

  • setting up atom positions,

  • saving samples to VNL files,

  • using some of the functionality in the ATK.KohnSham module.

If these concepts are new to you, we recommend that you consult the tutorial Geometry of a water molecule before you continue any further.

Introduction to periodic systems

A crystal is characterized by being periodic. If we take the unit cell associated with the crystal, and repeat this in all three spatial dimensions, this operation will generate the actual crystal. Put differently, we say that the crystal is invariant to translational symmetry operations. Our first task will be to see how such a periodic structure is constructed using ATK.

We will do this step-by-step by creating a NanoLanguage script which

  1. Sets up a bulk silicon crystal

  2. Prints out the crystal configuration

  3. Saves the crystal sample to a VNLFile for later use.

Creating and examining a periodic configuration

Our first step is to set up the Bravais lattice associated with the silicon crystal. The following script shows how this is done using NanoLanguage:

from ATK.KohnSham import *

# Create Bravais lattice 
lattice_constant = 5.43*Angstrom
fcc_lattice = FaceCenteredCubic(lattice_constant)

The first line (after the import statement) defines a keyword a, which is set to the lattice constant of silicon. After that, we define the lattice for the silicon crystal as a FaceCenteredCubic object, which we then assign to the variable fcc_lattice. NanoLanguage has native support for all 14 Bravais lattices, so we never have to to specify a Bravais lattice in terms of its unit cell vectors.

So far, we have just created an empty lattice without any atoms. The atoms of the crystal are specified using a BulkConfiguration object in the following way:

# Specify bulk configuration with 2 silicon atoms
# as basis in a FCC Bravais lattice
si_fcc = BulkConfiguration(
    bravais_lattice = fcc_lattice,
    elements = 2*[Silicon],
    fractional_coordinates = [
        [0.00, 0.00, 0.00], 
        [0.25, 0.25, 0.25]
        ]
    )

The arguments to BulkConfiguration are:

  1. the variable fcc_lattice, which we used above to define the lattice type of the crystal

  2. a list of elements specifying the atoms in the primitive basis

  3. a list of positions for the corresponding atoms where each set of coordinates in the list also is given as a Python list.

There are two NanoLanguage methods for specifying atomic positions, namely

  • absolute Cartesian coordinates

  • relative coordinates expressed in the three unit vectors

We generally recommend the latter approach, as relative coordinates are independent of lattice constants, but cases exists where Cartesian coordinates will be easier to use. See the reference documentation for BulkConfiguration for more specific details.

Now that we have defined the bulk crystal, it is a good idea to make sure that it has been set up properly. For convenience, let us define a general Python function that takes one single argument, namely the bulk configuration, and prints atomic symbols for all atoms, the corresponding positions, as well as the unit cell vectors. Below is a Python function printBulkConfig() that takes care of the task

from ATK.KohnSham import *

TAB = "\t"

def printAtomsInBulk(bulkSample):
    for (atom,pos) in zip(bulkSample.elements(), bulkSample.cartesianCoordinates()): 
        print atom.symbol(), TAB,
        for i in range( len(pos) ):
            print pos[i].inUnitsOf(Ang), TAB,
        print

def printLatticeVectors(lattice):
    for vector in lattice.primitiveVectors():
        for i in range( len(vector) ):
            print vector[i].inUnitsOf(Ang), TAB,
        print

def printBulkConfig(bulkSample):
    print "\nPrimitive vectors (in Angstrom)"
    printLatticeVectors(bulkSample.bravaisLattice())
    print "\nAtom list (positions in Angstrom)"
    printAtomsInBulk(bulkSample)
    print

printBulkConfig.py

Download and save the above script under the name printBulkConfig.py. The function printBulkConfig() can then be imported and used subsequently by adding the two lines

# Display Cartesian coordinates for silicon crystal
from printBulkConfig import printBulkConfig
printBulkConfig(si_fcc)

to the main script.

Finally, we store the silicon crystal sample in a VNLFile, so we can reuse it later in a different script. This is accomplished using the following two lines:

# Store crystal structure in a VNL file for later use
vnl_file=VNLFile("si.vnl")
vnl_file.addToSample(si_fcc,"si_fcc")

The entire script for setting up the fcc silicon crystal now reads

from ATK.KohnSham import *

# Create Bravais lattice 
lattice_constant = 5.43*Angstrom
fcc_lattice = FaceCenteredCubic(lattice_constant)

# Specify bulk configuration with 2 silicon atoms
# as basis in a FCC Bravais lattice
si_fcc = BulkConfiguration(
    bravais_lattice = fcc_lattice,
    elements = 2*[Silicon],
    fractional_coordinates = [
        [0.00, 0.00, 0.00], 
        [0.25, 0.25, 0.25]
        ]
    )

# Display Cartesian coordinates for silicon crystal
from printBulkConfig import printBulkConfig
printBulkConfig(si_fcc)

# Store crystal structure in a VNL file for later use
vnl_file=VNLFile("si.vnl")
vnl_file.addToSample(si_fcc,"si_fcc")

bulk-setup.py

Download and save the bulk-setup.py script to your disk, and process it with ATK by issuing the command

atk bulk-setup.py

This will print the following lines of output:

Primitive vectors (in Angstrom)
0.0     2.715   2.715   
2.715   0.0     2.715   
2.715   2.715   0.0     

Atom list (positions in Angstrom)
Si      0.0     0.0     0.0     
Si      1.3575  1.3575  1.3575

and generate the VNLFile si.vnl, which will now be located on your disk. The file si.vnl can be imported into Virtual NanoLab if you would like to study a 3D rendering of the silicon crystal. A 3D rendering of the silicon structure is in Figure 18.

The structure of the silicon fcc crystal, plotted using Virtual NanoLab. The unit cell of the crystal is indicated with a wire box.

Figure 18: The structure of the silicon fcc crystal, plotted using Virtual NanoLab. The unit cell of the crystal is indicated with a wire box.


k-point sampling

The fact that the periodicity is an inherent feature of a crystal structure, has very important consequences for the numerical treatment of such structures; these must be considered when we calculate periodic systems in ATK.

The fact that a crystal is invariant in certain directions with respect to translations by a unit cell length implies that the wave vector is a conserved quantity (note that, you can have 1D and 2D periodic systems also). As a result, the eigenfunctions obey Bloch's theorem and are therefore labeled by the wave vector k (in 1D, k is just a scalar, the wave number k). This also means that the Schrödinger equation can be solved independently for each possible value of k.

When constructing the total electron density \rho({r}), it is, however, necessary to know the solution in all points of the first Brillouin zone, since

\displaystyle

	\rho(r) = \sum_{n} \int d{\mathbf{k}}\ f_{n\mathbf{k}}\ |\psi_{n\mathbf{k}}(r)|^2,

where f_{n\mathbf{k}}\, is the occupation of the state k in the n'th band and \psi_{n\mathbf{k}}(r)\, is the corresponding Bloch function.

In practice, the integral is approximated by a sum, which is then evaluated over a finite number of k-point samples in the Brillouin zone. Clearly, as the number of k-points increases, the sum approximates the 'true' value of the above integral to a greater accuracy.

In the present version of ATK, the k-point sampling is specified in terms of how many points that should be used in the three directions of the unit cell vectors. These directions need not always coincide with the Cartesian axes.

The k-point sampling is set up by calling the function brillouinZoneIntegrationParameters() as

brillouinZoneIntegrationParameters([n1,n2,n3])

with a list argument containing three positive, non-zero integers, where each integer corresponds to the number of k-points for the respective axis directions.

This sets up a mesh with n_1 \times n_2 \times n_3 k-points uniformly spaced in the first Brillouin zone (i.e. in the reciprocal lattice) following the Monkhorst-Pack scheme [7].

Exactly how many k-points to choose depends to a large extent on the geometry and dimensionality of the system. For example, the bigger the unit cell is, the fewer points are generally needed. Please note, that it is only necessary to sample the directions in which the crystal is periodic. A carbon nanotube, for example, only requires a large k-point sampling along the tube axis (and a single point, of course, in the two transverse directions).

Let us use the script below to check how many k-points are needed in order to obtain an acceptable value of the total energy in a silicon crystal. By acceptable we mean within a certain degree of accuracy.

from ATK.KohnSham import *

# Restore previously created silicon structure
si_file = VNLFile("si.vnl")
si_fcc  = si_file.readAtomicConfigurations()["si_fcc"]

# Define basis set and mesh cut-off
basis = basisSetParameters(type=SingleZeta)
mesh  = electronDensityParameters(mesh_cutoff=75.0*Rydberg)

print '\nConverging bulk fcc Si in number of k-points'
print 'kpts\t\tTotal energy (eV)\t(rel. diff)'
print '---------------------------------------------'

convergence_criterion = 1.0e-5

kpts = 0
E = Eprev = 1.0*eV

# Loop determine number of k-points for convergence
while ( (kpts<1) or (abs( (Eprev-E)/E ) > convergence_criterion) ):
    # Increase number of k-points for each iteration
    kpts = kpts + 1
    Eprev = E
    kgrid = brillouinZoneIntegrationParameters(
        monkhorst_pack_parameters=3*[kpts]
        )
    dft_method = KohnShamMethod(
        basis_set_parameters=basis,
        electron_density_parameters=mesh,
        brillouin_zone_integration_parameters=kgrid
        )
    # Perform actual calculation
    if kpts > 1:
        # Use the previous scf as initial calculation
        scf = executeSelfConsistentCalculation(
            atomic_configuration=si_fcc,
            method=dft_method,
            initial_calculation=scf,
            )
    else:
        # For (1,1,1) there is no initial calculation to use
        scf = executeSelfConsistentCalculation(
            atomic_configuration=si_fcc,
            method=dft_method
            )
    E = calculateTotalEnergy(scf)
    if kpts > 1:
        print 3*(kpts,), '\t', E.inUnitsOf(eV), '\t\t', abs((Eprev-E)/E)
    else:
        print 3*(kpts,), '\t', E.inUnitsOf(eV)
        
print '\nConverged with',kpts,'k-points in each direction\n'

bulk-kpoint.py

Before we actually run the script, a few remarks are in place regarding its structure

  1. First, we load the silicon sample from our previous run and store this configuration in the variable si_fcc.

  2. The basis set and the mesh cut-off are stored in the variables basis and mesh using the respective functions basisSetParameters() and electronDensityParameters(). Note, to speed up the calculation a bit, we have chosen both a small basis set and mesh cut-off.

  3. A Python while loop is constructed. In each iteration step, the number of k-points is incremented by 1 in all directions, after which the total system energy E is calculated. When the relative difference between the current energy E and the energy Eprev of the previous step is lower than the variable convergence_criterion, the while-loop is terminated.

  4. In each iteration step, the method for the self-consistent calculation is set using a KohnShamMethod object, which is initialized with the chosen basis set, mesh cut-off, and the number of k-points as specified by the respective variables basis, mesh, and kgrid. The returned method object for solving the Kohn-Sham equations is assigned to the variable method.

  5. A self-consistent calculation is now performed for the silicon system. This is done in two steps:

    1. For the first set of k-points (1,1,1), we calculate the electron density from scratch and change the Boolean value of use_old from False to True.

    2. In every subsequent calculation we reuse the density calculated in the previous calculation by passing the variable old_scf to executeSelfConsistentCalculation() using the optional keyword initial_calculation. This reduces the number of steps required for the calculation to converge.

    In every step we pass the system si_fcc as the argument to executeSelfConsistentCalculation() along with the parameters previously specified in our dft_method.

  6. The total energy is evaluated using the function calculateTotalEnergy(), which is called with the object calculation returned by the executeSelfConsistentCalculation() function.

Observe, that the above script explicitly uses a single zeta basis set. In this case, this simple basis set is sufficient for determining the optimal k-point sampling.

Now, download and save the script, and process it through ATK. You will get the following output

Converging bulk fcc Si in number of k-points
kpts            Total energy (eV)       (rel. diff)
---------------------------------------------
(1, 1, 1)       -337.095933988
(2, 2, 2)       -356.154925831          0.0535132058008
(3, 3, 3)       -355.723973423          0.00121147979962
(4, 4, 4)       -356.327501856          0.00169374642619
(5, 5, 5)       -356.273593339          0.000151312132773
(6, 6, 6)       -356.329983582          0.000158252871966
(7, 7, 7)       -356.32241598           2.1238075445e-05
(8, 8, 8)       -356.330060267          2.14528268082e-05
(9, 9, 9)       -356.328781041          3.59001521208e-06

Converged with 9 k-points in each direction

From the results, we see, that a (9,9,9)-grid gives excellent convergence for the silicon crystal. Inspecting the numbers a bit closer reveals that the accuracy is acceptable even for a (4,4,4)-grid. It is noteworthy how small the change in the total energy is when going from a (4,4,4) to a (5,5,5) grid, or more generally from an even number in each direction to an odd number. Clearly the (5,5,5) calculation takes longer time, but does not improve the accuracy much. This is a typical feature when using a small number of k-points in all directions.

Reusing the old density can be used to save considerable time. Say, you have an old calculation which was converged using a given set of k-points but now you want to use an even larger mesh to obtain a better accuracy. Using the initial_calculation keyword this can be done in a few iterations instead of several steps if ATK had to perform the calculation from scratch.

The band structure of bulk Si

In this example, we will demonstrate some general steps involved in ATK based band structure calculations. To illustrate the basic principles, it is sufficient to calculate and plot the band structure for a small portion of the Brillouin zone. The power of NanoLanguage, however, easily enables you to create a script that generates and plots the band structure for the entire zone, for example in order to produce a plot similar to the one shown for silicon in Figure 19.

The band structure of silicon, plotted in the entire Brillouin zone along the traditional symmetry directions for the fcc lattice. The data and plot script for the 2D plotting program gnuplot were generated entirely with NanoLanguage.

Figure 19: The band structure of silicon, plotted in the entire Brillouin zone along the traditional symmetry directions for the fcc lattice. The data and plot script for the 2D plotting program gnuplot were generated entirely with NanoLanguage.


The self-consistent calculation

The first step is to perform the self-consistent calculation of the silicon crystal. We know that a k-point mesh of (4,4,4) is sufficient to obtain converged total energies; in the tutorial Calculating molecular properties, we examined the effect of choosing a different basis set and exchange-correlation functional. For improved accuracy, we will here use the GGA-PBE- approximation for the exchange-correlation functional, as well as the DoubleZetaPolarized basis set. The following script starts the self-consistent field calculation and stores the results:

from ATK.KohnSham import *
import ATK

# Restore previously defined Silicon structure
si_fcc = VNLFile("si.vnl").readAtomicConfigurations()["si_fcc"]

#Specify NetCDF file for storing results
ATK.setCheckpointFilename('si_gga_dzp.nc')

# Define new KohnShamMethod using a non-default
# exchange-correlation potential (GGA)
method = KohnShamMethod(
    exchange_correlation_type = GGA.PBE,
    basis_set_parameters = basisSetParameters(type=AtomicOrbitals.DZP),
    brillouin_zone_integration_parameters = \
        brillouinZoneIntegrationParameters(monkhorst_pack_parameters = [4,4,4])
    )

# Perform the actual calculation
gga_calculation = method.apply(atomic_configuration = si_fcc)

bulk-gga_calculation.py

The results of the calculation are stored in the NetCDF file si_gga_dzp.nc. Note, that the checkpoint file is overwritten without warning if it already exists! More information can be found in the reference manual under setCheckpointFilename().

After the calculation has finished, we can use the file si_gga_dzp.nc as the starting point for calculating the band structure.

Calculating and visualizing Bloch states

One of the analysis functions for studying bulk samples is the calculateEigenstates() function. By specifying the quantum number for specific states at a desired k-vector, it is possible to visualize the eigenstates of silicon using VNL.

If we restore our DFT calculation from above (stored in the si_gga_dzp.nc file), the eigenstates can be calculated using a minimum of effort.

from ATK.KohnSham import *

# Restore old NetCDF results
gga_calculation = restoreSelfConsistentCalculation('si_gga_dzp.nc')

# Calculate eigenstates
eigenstates = calculateEigenstates(
    self_consistent_calculation = gga_calculation,
    quantum_numbers = ( (0,1), (0.375,0.375,0.75) )
)

# Store eigenstates in a VNL file
vnl_file=VNLFile("si.vnl")
for eigenstate in eigenstates:
    vnl_file.addToSample(eigenstate, "si_fcc")

bulk-eigenstates.py

The quantum numbers given to the VNLFile function are given as a tuple sequence of

  • the band index (a non-negative integer, less than the total number of bands), and

  • k-points, in units of the reciprocal lattice vectors.

As we see in the example, the band index can be a sequence in itself, to get the Bloch functions of several bands at the same k-point. Comparing with the list of symmetry points below, we see that the selected point is the 'K' point. The eigenstates are stored in the same VNLFile using the same sample name si_fcc as used above to associate the eigenstates with the silicon structure. If there are several band indices, the returned object eigenstates is a list of Bloch functions, which is why we loop over eigenstate when saving the plots to the VNL file.

The generated VNLFile can now be imported into Virtual NanoLab, which then can be used to generate a 3D rendering of the eigenstates similar to the plot shown in Figure 20.

Eigenstates visualized using VNL.

Figure 20: Eigenstates visualized using VNL.


Setting up symmetry points and lines

Our next task is to calculate the band structure of silicon along a specific direction in the Brillouin zone.

Because the complete script for calculating and plotting the band structure of the silicon crystal is somewhat lengthy, we will go through the key elements of the script step-by-step. We start by creating a dictionary containing common symmetry points in the Brillouin zone for the fcc lattice

# Define common symmetry points in the Brillouin zone of silicon
fcc_symmetry_points = {
    'G' : (0.000, 0.000, 0.000),
    'X' : (0.500, 0.000, 0.500), 
    'W' : (0.500, 0.250, 0.750), 
    'L' : (0.500, 0.500, 0.500), 
    'K' : (0.375, 0.375, 0.750), 
    'U' : (0.625, 0.250, 0.625) 
}

The points are given in lattice coordinates independent of the lattice constant. ATK will automatically convert these points into Cartesian coordinates, using the reciprocal lattice vectors of the fcc lattice. We have used the label 'G' to denote the Γ-point of the Brillouin zone.

Now, we want to create a "path" of points in the Brillouin zone connecting two points in the above dictionary. For simplicity, we only consider the path between the Γ-point and and the point 'X'.

#Create line in Brillouin zone along which the
#band structure is calculated
start_point = fcc_symmetry_points['G']
end_point = fcc_symmetry_points['X']
band_k_points = lineOfVectorPoints(start_point,end_point,25)

Here, we have used a function lineOfVectorPoints() (defined below), which takes three arguments

  • The first two arguments are the two points that should be connected.

  • The third argument is an optional integer N, that specifies the number of points to include in the path between the two points

The function lineOfVectorPoints() returns a list of N points positioned on a straight line between the two given points. The definition of lineOfVectorPoints() is simple, namely

def lineOfVectorPoints(start_point,end_point,N=20):
    '''
    @return: A list containing (default: 20) k-points between two
    specified points in the Brillouin zone.
    '''
    from numpy import array
    points = []
    for i in range(N):
        new_point = (float(i)/(N-1),)*3*(array(end_point) - array(start_point))+start_point
        points.append(new_point)
    return points

The statement

    from numpy import array

is necessary in order to convert the tuples of points in the fcc_symmetry_points dictionary to numeric arrays; otherwise, we would not be able to calculate the difference between the two points. The returned list includes the end-point; to exclude it (which is sometimes desirable), simply use N instead of N-1 in the denominator.

Calculating the band structure

Here are the lines of the script is the part where the actual band structure gets calculated

#Restore old SCF calculation
calculation = restoreSelfConsistentCalculation('si_gga_dzp.nc')
band_structure = calculateEnergyBands(
    self_consistent_calculation = calculation,
    k_point_list = band_k_points)

The input to the function calculateEnergyBands() should be prepared in the form of a list of points, given as tuples (k_1, k_2, k_3); this input corresponds exactly to the list returned by the user defined function lineOfVectorPoints() described above, which sets up the sequence of points between the symmetry points.

The following three lines of code appends the energy bands of our silicon crystal to a list. This may then be used for exporting the band structure to other applications.

#Extract calculated bands into a list
energy_bands = []
for band_index in range(band_structure.numberOfBands()):
    energy_bands.append(band_structure.band(band_index))

The number of bands in the band structure depends on the crystal lattice basis (i.e. the atoms in the unit cell) as well as on the size of the basis set (in our case DoubleZetaDoublePolarized). The latter quantity is difficult to know beforehand had we not written the script ourselves, and we therefore use the method numberOfBands() to query for the number of bands contained in the band structure. Once we know the number of bands, we can easily iterate through them one by one, and extract the band structure for each band. To accomplish this, we use the method band(), which takes the band index as argument. The band index ranges from 0 to n-1, where n is the total number of bands.

Exporting the band structure

We now combine the small pieces of code described above into the following complete script for calculating the energy bands of silicon.

from ATK.KohnSham import *

def lineOfVectorPoints(start_point,end_point,N=20):
    '''
    @return: A list containing (default: 20) k-points between two
    specified points in the Brillouin zone.
    '''
    from numpy import array
    points = []
    for i in range(N):
        new_point = (float(i)/(N-1),)*3*(array(end_point) - array(start_point))+start_point
        points.append(new_point)
    return points

# Define common symmetry points in the Brillouin zone of silicon
fcc_symmetry_points = {
    'G' : (0.000, 0.000, 0.000),
    'X' : (0.500, 0.000, 0.500), 
    'W' : (0.500, 0.250, 0.750), 
    'L' : (0.500, 0.500, 0.500), 
    'K' : (0.375, 0.375, 0.750), 
    'U' : (0.625, 0.250, 0.625) 
}

#Create line in Brillouin zone along which the
#band structure is calculated
start_point = fcc_symmetry_points['G']
end_point = fcc_symmetry_points['X']
band_k_points = lineOfVectorPoints(start_point,end_point,25)

#Restore old SCF calculation
calculation = restoreSelfConsistentCalculation('si_gga_dzp.nc')
band_structure = calculateEnergyBands(
    self_consistent_calculation = calculation,
    k_point_list = band_k_points)

#Extract calculated bands into a list
energy_bands = []
for band_index in range(band_structure.numberOfBands()):
    energy_bands.append(band_structure.band(band_index))

#print band energies in units of eV
for i in range(len(band_k_points)):
    print band_k_points[i][0],
    for band in range(band_structure.numberOfBands()):
        print energy_bands[band][i].inUnitsOf(eV),
    print

bulk-calculateBands.py

The script will print the energy bands of silicon to your terminal. If you want to store the information in a file, you should redirect the script output to a file in the following manner:

atk bulk-calculateBands.py > si_bands.dat

The first data column contains the value of k_x (band_k_points[i][0]), followed by the energy eigenvalues for each band associated with the chosen k-point in each column.

The obtained band structure of silicon is shown in Figure 21. The plot indicates that silicon would have an indirect band gap of approximately 2 eV. Experimentally, however, silicon is known to have an indirect band gap of 1.13 eV.

The GGA band structure of silicon between the Γ and X point, generated from the data produced by the presented NanoLanguage script.

Figure 21: The GGA band structure of silicon between the Γ and X point, generated from the data produced by the presented NanoLanguage script.


The reason for this discrepancy is, in this case, that the SingleZeta basis set (which for silicon contains the sp3 valence electrons) is not versatile enough to describe a silicon crystal. In particular, it does not have enough angular degrees of freedom to correctly form the strong tetragonal covalent bonds present in the diamond structure of silicon. Going to DoubleZeta, one would only add radial degrees of freedom, and thus it is necessary as a minimum to have a SingleZetaPolarized (but ideally DoubleZetaPolarized) to get accurate results for silicon.

If you re-run the above calculation with the DoubleZetaPolarized basis set (remember to change the name of the checkpoint file!), you should find a band gap below 0.6 eV, i.e. the band gap is severely underestimated. This is a well known artifact of DFT calculations using either the LDA or GGA-PBE exchange correlation functional.