Lab 4: Solids and Electronic Band Structures

Back to Course Overview

Reminder Don’t forget to copy the lab04 folder from /opt/Courses/MSE404/lab04 to your home directory.

This week we are going to start doing some calculations on solids, i.e. fully periodic systems. Many of the principles will be the same, but as you will see there are a few things that need to be done differently.

Diamond

As our first example of a solid we’re going to look at diamond. You can find an annotated input file in 01_carbon_diamond/C_diamond_detailed.in. Much of the input is similar to the molecular input files we have seen, just with a few key differences. These include the fact that we have specified an fcc lattice (using ibrav), we have listed the atomic positions in fractional coordinates (ATOMIC_POSITIONS crystal) and have specified the k-point sampling (K_POINTS automatic). ATOMIC_POSITIONS crystal means that the coordinates of the atoms are expressed in relative coordinates of the unit cell vectors as defined via the ibrav and A variables. Absolute Cartesian coordinates (x,y,z) and fractional coordinates (fx,fy,fz) are related by: (x,y,z) = fx*v1 + fy*v2 fz*v3, where v1,v2,v3 are the unit lattice vectors. For diamond, which has the same atomic structure as Zinc Blende, we set ibrav=2, i.e. face-centred cubic (fcc) Bravais lattice, and a basis of two carbon atoms, whose fractional coordinates are (0,0,0) and (1/4,1/4,1/4). Internally, Quantum ESPRESSO sets the the fcc lattice vectors as:
v1 = (A/2)(-1,0,1), v2 = (A/2)(0,1,1), v3 = (A/2)(-1,1,0). Hence, the absolute Cartesian coordinates for the two carbon atoms are given by : C1: (x,y,z) = 0*v1 + 0*v2 + 0*v3 and C2: (x,y,z) = 1/4*v1 + 1/4*v2 + 1/4*v3.

Task

  • Run the input file for diamond. There are a couple of extra things to notice in the output file.
  • The output lists the automatically generated k-points. We have 10 here since crystal symmetries have been taken into account after generating the 64 points on the 4x4x4 grid we requested.
  • The eigenvalues are output for each k-point.

Silicon

Silicon crystallises in the same structure as diamond. So in this case, all we need to do is change the lattice length to the silicon value, and change the atoms from carbon to silicon. And of course we need to specify a silicon pseudopotential rather than use the carbon one. If you use diff you can see the differences between the silicon input file relative to the carbon diamond input file. If you’re in the lab02 directory remember you can type diff 01_carbon_diamond/C_diamond.in 02_silicon/Si.in. You’ll be able to see that we’ve changed four lines in our input file and everything else is the same. So switching to a system with the same structure but different atoms can be quite simple. Note we’ve just used the measured lattice constant for both systems; in later labs we’ll see how to find the lattice constant predicted by DFT.

For periodic systems it is important to test the convergence with respect to k-point sampling. The folder 02_silicon contains a template file and script to check the convergence of the total energy with respect to the k-point sampling in silicon.

Task

  • Run the script and plot the convergence of total energy with respect to k-point sampling.

Optional Task

  • For every periodic system you simulate, you should converge both the cut-off energy and k-points. Try adapting one of the scripts from last week to also converge the energy of silicon with respect to the cut-off energy. How does the convergence behaviour of the two parameters compare?

The Electronic Band Structure

While the electronic density obtained from DFT is meaningful, the Kohn-Sham states are not strictly the electronic states of the system. Nonetheless, they are in practice often a good first approximation of the electronic states of a system, so can be useful in understanding the properties of a system.

We have now seen how to converge our calculations with respect to the sampled k-point grid density. And you’ll have seen in the calculations you have done that the calculated eigenvalues are a bit different at each calculated k-point. Examining how the state energies change from one k-point to the next can tell us useful things such as if a material is likely to have a direct or indirect optical gap for example. For this we need to visualize how the energies of the states vary with k-point. The usual way this is done is to plot the band energies along lines between the various high-symmetry k-points in the Brillouin zone. The details of how this can be done is beyond the scope of this course, but an outline is given at the end of this lab.

If you have a particular structure and you want to find out which are the important k-points, then this website is a useful tool.

The directory 03_bandstructure/01_diamond contains input files to calculate and plot the band structure of diamond.

This a three step process:

  1. Calculate a converged density with a standard self-consistent field (SCF) calculation. This means solving the Kohn-Sham (KS) equations on a uniform k-grid self-consistently. At the ith-iteration we compute the electron density n^(i) from the KS states and use this density to generate the KS potential for the next iteration V_KS^(i+1) and solve KS equations with this new potential. The loop is terminated when |n^(i+1) - n^(i)| < thr, where thr is the required convergence threshold.

  2. Use that density to perform a non self-consistent (NSCF) calculation for k-points along chosen high-symmetry lines. In a NSCF calculation, the energy is not minimised with respect to the electron density, instead to obtain the KS energies and KS states for a particular k-point, one generates the KS potential from the converged electron density (obtained in the SCF calculation) and diagonalises the resulting KS Hamiltonian at that specific k-point.

  3. Extract the energies from this calculation and convert it to a dataset we can plot.

  • We can do step 1 above exactly as previously, as you can see in 01_C_diamond_scf.in but for step 2, we’ll need to pick some suitable high-symmetry points.
  • The diamond lattice is FCC. In 02_C_diamond_nscf.in we have set up an input file with the path already set as Γ-K-X-Γ'-L-X-W-L where Γ' indicates the gamma point in a different Brillouin zone. You could use this same path for any other FCC system.
  • The other inputs that are changed in this file with respect to the scf calculation.
    • We set the calculation = bands for a band structure plot.
    • We choose the crystal_b option for the K_POINTS section. This means that we will enter just the high-symmetry points and points in between will automatically be generated. Then we give the number of high symmetry points we’ll enter with the coordinates of each in reciprocal lattice coordinates followed by the number of points to generate between it and the next point. Note, we have cut down the number of points for segments that are shorter.
  • For the third step we’ll be using the bands.x tool from the espresso package. The input for this is 03_C_diamond_scf.in. It has a single section: BANDS, and we can just accept all defaults. As usual full details of the inputs and options are given in the INPUT_BANDS.txt file in the quantum espresso documentation folder.

To perform the calculation, we run pw.x with the first two input files in turn. And then run bands.x with the third. There’s a simple script provided as run_all.sh that will explicitly run the three steps. The third step will produce several files, with different output formats. The one we’ll use is bands.out.gnu which can be used easily with gnuplot.

We have also set up a script file for gnuplot as plotbands.gplt:

set encoding utf8 # This lets us use the Gamma symbol directly

# The locations of the tics are given in the output of the bands.x calculation
set xtics ("Γ" 0.0000, "K" 1.0607,"X" 1.4142, "Γ" 2.4142, "L" 3.2802, "X" 4.1463, "W" 4.6463, "L" 5.3534)

# This gives us a full vertical line
set grid xtics lt -1 lw 1.0

# We don't need a legend
unset key

# set a label and a title
set ylabel "Energy (eV)"
set title "Carbon Diamond Electronic Band Structure"

# This tells gnuplot to plot all the points from this file connected with lines
plot "bands.out.gnu" with lines

# And if you run this script directly as an argument to gnuplot, rather than
# by loading it within gnuplot, you can uncomment the following to keep the
# plot window open until clicked. You can save to a file from here.
pause mouse

The bands.x calculation sets the x-axis scale so that the distance between high symmetry points corresponds to the actual Cartesian distance between the points. To find the x-coordinates for the high-symmetry points you can check the output of the bands.x calculation. Note in this case it doesn’t detect K as a high-symmetry point so we need to add that in ourselves. If you examine the points given in the second input, you’ll see the chosen gamma to K to X path is actually a single straight line in reciprocal space, so we can find where K should fall along this path.

You can bring up the gnuplot graph with gnuplot plotbands.gplt. You’ll notice we still need to shift the vertical axis to put the valence band maximum at 0. There’s a second gnuplot script called plotbands_shifted.gplt where you can see how this could be done with gnuplot. As we could see the valence band max was at gamma (the first point on our path), we could read the value of the energy at this point from one of the other output files, bands.out.

Task

  • Follow the above steps to generate a plot of the electronic band structure of diamond.

Optional Task

  • Generate a plot of the electronic band structure of silicon. Note this has an fcc Bravais lattice as did carbon diamond, so you can use the same set of high symmetry k-points for your plot.

Density of States

Now let’s analyse the electronic states by computing the density of states (DOS). This is a little easier to visualise and shows how many electronic states (in fact Kohn-Sham states for our DFT calculation) are at a given energy. More precisely, the DOS tells us how many electronic states, for a system of volume V, can be occupied in a small (infinitesimal) energy range near a specific energy. The DOS is directly related to the band structure as high DOS at a specific energy level means that a large number of electronic states can be occupied. Hence, bands with large energy dispersion in the Brillouin zone result in low DOS, whereas less dispersive (more flat) bands result in high DOS. In insulators and semiconductors the DOS is zero inside the band gap, as there are no available states in that energy range. In a similar way to the electronic band structure, we produce the density of states plot in three steps.

  1. Perform a self consistent calculation as before, producing a converged charge density.
  2. Take the density calculated calculated in the previous step and use it to perform a non-self-consistent calculation on a more dense grid of k-points. We want a good representation of how the state energies vary as we move around the Brillouin zone so we use a much denser grid here than we need to obtain a converged density in the previous step.
  3. Convert the state energies calculated on this dense k-point grid to a density of states.
    • Since we still only have a representation of the state energies at a fixed set of points in reciprocal space, we need to interpolate between them in some sensible way if we turn this into a count of the total number of states at an arbitrary energy.
    • The most common way this is done is to use some energy broadening scheme. Then the number of states at some arbitrary energy is given as the sum over energies at the calculated k-points times some weighting given by the broadening. In practice this is quite fast and straight-forward, although you’ll need to tune the broadening energy so that your calculated density of states is smooth in the correct way.
      • If you use too large a broadening, you may smear out important features.
      • If you use too small a broadening you may introduce spurious features and your density of states plot will look very bumpy.
      • In principle you would want the smearing to be comparable to the typical change in energy of a state from a k-point to its neighbours. In practice though it’s easiest to just try different values until it looks right.
    • The other way to interpolate is to use the so-called tetrahedron method. Essentially this corresponds to doing a three dimensional linear interpolation from a regular grid of values. This calculation can be noticeably slower than using a broadening but there is no need to to worry about using the correct smearing. The density of states will simply become more finely featured as you increase the density of the k-point grid in the non-self-consistent calculation.
    • It’s important to note that in a real measurement of the density of states of a system the there is an implicit broadening that comes from
      1. Electron-phonon coupling: the states are not simply at a fixed energy, but will have some distribution as the atoms vibrate.
      2. Any measurement probe will have a finite energy width associated with it, which will limit how finely it can resolve density of states features.
    • So while tetrahedron may seem the more accurate approach, you shouldn’t necessarily think of it as a more correct representation of a real system.

The directory 04_densityofstates/01_diamond contains an example set of input files to calculate the density of states for diamond following these three steps, and using a Gaussian broadening for the density of states calculation.

  1. 01_C_diamond_scf.in is a usual scf input file for diamond as you have seen before.
  2. 02_C_diamond_nscf.in is the input file for the non-self-consistent calculation with pw.x. Compared to the scf calculation we’ve changed the following:
    • We need to change the calculation type to nscf with calculation = nscf
    • We’ve increased the number of bands here to 8 (pw.x will just include enough bands for the occupied states in an insulator or semiconductor by default, and we’d like to see some of the conduction band density of states too).
    • We’ve increased the k-point sampling to a 20x20x20 grid, and we have removed the shift. Many systems have a valence band maximum or conduction band minimum at the gamma point, so it is good to ensure it’s explicitly included in the grid.
  3. 03_C_diamond_dos.in is the input file for dos.x. This code input file requires just a DOS section. Here we’ve specified the Gaussian broadening to use with degauss in Rydberg. We also specify DeltaE which is the spacing between points in the output file, in eV. Note - we’ve picked values for these of similar magnitude despite their different units. In fact if degauss is not specified, and no broadening scheme is used in the DFT calculation, degauss will take the value of DeltaE by default. You can check the documentation file INPUT_DOS.txt for more details.

Now we need to run all three inputs, the first two with pw.x and the third with dos.x. There’s a simple script to do these three steps explicitly in run_all.sh.

The final step produces a file named pwscf.dos by default. This is a simple text file you can plot in whatever software you like. It has three columns:

  1. Energy (eV)
  2. Density of States (states/eV)
  3. Integrated Density of States (states)

It is customary to shift the x-axis in the plot such that the Fermi energy or valence band max is at 0. While a value for the Fermi level is given in the file header of the generated pwscf.dos, this is determined in a simple way from the integrated density of states. It may be worth obtaining this from a separate calculation using a relatively small broadening if you’re looking a metallic system, while for semiconductors and insulators you could find the maximum valence band state energy manually. If you’re plotting in gnuplot you can shift the x-axis origin within the plot command:

plot "pwscf.dos" using ($1-13.180):2 with lines

Where we have used 13.180 as the value of the Fermi energy. If you want to plot the integrated DOS using the right hand axis you can do the following:

set ytics nomirror
set y2tics
set xlabel "Energy (eV)"
set ylabel "Density of states"
set y2label "Integrated density of states"
set key top left
plot "pwscf.dos" using ($1-13.180):2 with lines title "Density of States", \
     "pwscf.dos" using ($1-13.180):3 axes x1y2 with lines title \
     "Integrated density of states"

Task

  • Calculate and plot the density of states and integrated density of states for diamond.

Optional Task

  • Calculate and plot the density of states and integrated density of states for silicon, making sure to shift the 0 of the x-axis to be the highest energy in the valence bands. Save it as a png graphic.
    • How does this compare to the density of states for carbon diamond?

Summary

  • In this lab we looked at how to calculate:
    • k-point convergence in solids
    • the electronic band structure of a solid
    • the electronic density of states of a solid
  • We have seen how several calculations may be chained together where the output of one is used as an input for a subsequent calculation.
  • We have used the bands.x and dos.x codes from the Quantum Espresso package.
  • We have done some more plotting in gnuplot.
  • We should always keep in mind that the Kohn-Sham eigenvalues as obtained from a DFT calculation do not correspond to the real interacting electron energy levels, but are often useful as a first approximation.

Extra: High Symmetry Points

Finding appropriate high symmetry points and their labels for a band structure plot is beyond the scope of this course, but generally you need to find the Brillouin zone for the system you’re interested in, along with the names of the high symmetry points, and how these should be represented in terms of your reciprocal lattice vectors.

These can often be looked up in a table for a given structure, while keeping in mind that even for the same structure you may come across papers where some zone boundary high symmetry points will be labelled differently. You can access these tables online at e.g. the Library of Crystallographic Prototypes by typing the name of the mineral into the search box. This will give you the space group number, which you could then use with e.g. the Bilbao Crystallographic Server to find the reciprocal space coordinates of the high symmetry points and their labels. Another good reference is TU Graz which has nice interactive visualizations of the most common types along with their labels.

For the diamond lattice example, we might do this as follows:

  • The diamond lattice is FCC. We can find find the space group number from http://aflow.org/CrystalDatabase/A_cF8_227_a.html and see that it is number 227.
  • We can enter this number at http://www.cryst.ehu.es/cryst/get_kvec.html and find several high symmetry points. Often, you’ll want to pick the same path that was chosen in some previous work to ensure you can reproduce it correctly. In 02_C_diamond_nscf.in we have set up an input file for the path Γ-K-X-Γ'-L-X-W-L where Γ' indicates the gamma point in a different Brillouin zone. Note - these labels don’t exactly match points given in the table. Many points have a number of equivalent positions on the Brillouin zone surface, and often different conventions can be used for different materials with the same structure.

Extra: Electronic states of molecules

We have already seen the occupied energy levels for molecules, but it is also possible to calculate unoccupied states, just as with solids, except we only have a single set of energy levels. So we don’t need to calculate anything like an electronic band structure or a density of states. If you recall our methane example previously, we just calculated the states at the gamma point. To find the energies of these states we can just read the output files from a pw.x scf calculation directly.

The directory 05_state_energies/01_methane has an example input file for methane. There are no new inputs used here; we perform the calculations as we’ve done previously except that we add two additional states to the calculation (the default would give us just the four doubly occupied orbitals) using the keyword nbnd.

Run the calculation as usual, and you’ll see in the output just before the total energy, the highest occupied and lowest unoccupied level energies are given in eV. And the energies of each band are given just above this. We can use this to estimate the HOMO-LUMO gap as the difference between these two values.

Last updated on 2021-11-08