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 the electronic states of a system, so can be useful in understanding the properties of a system.

## The Electronic Band Structure

In the previous lab we saw 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. Hopefully, you have a better understanding of what this k-point sampling represents from your lectures. 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.

## Electronic states of molecules

Of course for isolated molecules, this is all much easier. In this case, 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 `01_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).

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.

## The Electronic Band Structure of Diamond

The directory `02_bandstructure/01_diamond`

contains input files to calculate and plot the band structure of diamond.
This a three step process:

- Calculate a converged density with a standard scf calculation.
- Use that density to perform a non self-consistent (nscf) calculation for k-points along chosen high-symmetry lines (c.f. we performed a nscf calculation on a dense grid in the same way for the density of states).
- 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.

- We set the
- 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 an 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*

- 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. 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. In a similar way to the electronic band structure, we produce the density of states plot in three steps.

- Perform a self consistent calculation as before, producing a converged charge density.
- 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.
- 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 comparible to the typical change in energy of a state from a k-point to its neighbors. 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
- Electron-phonon coupling: the states are not simply at a fixed energy, but will have some distribution as the atoms vibrate.
- 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 `03_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.

`01_C_diamond_scf.in`

is a usual scf input file for diamond as you have seen before.`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.

- We need to change the calculation type to nscf with
`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:

- Energy (eV)
- Density of States (states/eV)
- 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 semi-conductors 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 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?

### Density of states using the tetrahedron method

As mentioned, the other way to generate the density of states is using the
tetrahedron method. To do this within quantum espresso, you must use
the tetrahedron occupation scheme for the DFT calculation, and then don’t
set any `degauss`

value in the `dos.x`

input file. This is outlined in the
`INPUT_DOS.txt`

help file.

An example where we have modified the previous diamond input files to
calculate the density of states using tetrahedron integration is given in the
`03_densityofstates/02_diamond_tet`

directory. You can check the `diff`

between the input files to see what we
changed.

If you run the three steps again here, you’ll once again produce a `pwscf.dos`

file as before. If you check the output file from the `dos.x`

calculation
you’ll also see the note `Tetrahedra used`

, whereas the previous output notes
details of the broadening.

Additionally if you look at the header of the `pwscf.dos`

file you’ll see a
different value for the Fermi energy than we obtained in the previous
calculation. This would mainly come from a small underestimation that happens
in the tetrahedron method. You can also see this in the integrated DOS in
comparison to the broadening case. For this reason a separate calculation of
the Fermi energy may often be useful. Note, you’ll get more well defined
features and a smoother curve systematically as you increase the k-point
sampling density. The calculation will take a lot longer; if you use e.g. a
60x60x60 grid in the nscf calculation it takes about 15 mins run in serial,
and the DOS calculation takes about 20 minutes, though this may be worth it if
producing an image for a paper or report. We’ve tried to keep run time for
example calculations provided down to a couple of minutes or less.

*Task*

- Plot the density of states and compare it that obtained using Gaussian broadening.

## Summary

- In this lab we looked at how to calculate:
- the energy levels of a molecule and hence the HOMO-LUMO gap
- 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 advanced 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 Tasks*

Check the convergence of the HOMO-LUMO gap of a methane molecule with respect to plane-wave energy cut-off.

- It’s all right just to do this manually for a few values of the energy cut-off.
The output file will list the highest occupied and lowest unoccupied orbital energy before the final total energy output.

If you do want to try to write a script to parse out the HOMO-LUMO gap from an output file along with the energy cut-off, based on the scripts you used last week, you could use the following awk command to parse the output from a set of output files:

```
awk '/kinetic-energy/{ecut=$4} /highest/{print ecut, $8-$7}' *out > gap_v_ecut.dat
```

This tells awk to save the fourth word on the line to the variable “ecut” when it finds a line containing the text “kinetic-energy”, and when it finds the text “highest” it outputs the value stored in this variable along with the difference between the 7th and 8th fields.

There are more details of the awk command in the summary of Linux commands given in Lab 1. You can also check the man page by typing

`man awk`

in a terminal.

Generate a plot of the electronic band structure of aluminium. Even though this has a different number of atoms per cell, and is a metal, the Bravais lattice is again fcc, so you can use the same high-symmetry k-points as before again. (Aluminium was covered in lab 2).

Using tetrahedra, calculate and plot the density of states and integrated density of states for silicon, and aluminium.

- In both cases, shift the 0 of the x-axis to the Fermi energy in your
plots, and be sure to label them correctly. Save them as
`png`

graphics.

- In both cases, shift the 0 of the x-axis to the Fermi energy in your
plots, and be sure to label them correctly. Save them as

### 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.