Header Ads Widget

Ab initio Calculations Using Q-Espresso Code

Last Posts

10/recent/ticker-posts

Lab 4: The Electronic Band Structure and Density of States

 

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:

  1. Calculate a converged density with a standard scf calculation.
  2. 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).
  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 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.

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

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

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.

 

https://eamonnmurray.gitlab.io/modelling_materials/lab04/ 

Post a Comment

0 Comments