4.4 Plotting the band structure and density of states of a solid

An example for plotting band structures and the density of states is contained in the example for bulk silicon (diamond lattice) in the testcases directory. Likewise, the fcc Al test case contains the necessary lines for control.in (just commented out).

This section contains a walk-through of the necessary steps and keywords to obtain this band structure and density of states. First, let us look at the relevant lines in control.in as given in the test case.

# output DOS
  output dos_tetrahedron -18.  0.  200

# output band structure
  output band 0.5   0.5   0.5    0.0   0.0   0.0    50  L     Gamma
  output band 0.0   0.0   0.0    0.0   0.5   0.5    50  Gamma X
  output band 0.0   0.5   0.5    0.25  0.5   0.75   50  X     W
  output band 0.25  0.5   0.75   0.375 0.375 0.375  50  W     K

To get an output of the band structure, use the keyword output band. Each ’band’ line covers a single line in reciprocal space, from a starting point to an end point and with a given number of points inbetween.

The total density of states of the system is obtained by the ’output dos_tetrahedron’ line, from a starting single-particle energy to an ending single-particle energy, with a specified number of points inbetween to create a smooth line (see section 3.52 for details). The tetrahedron here indicates the implementation of tetrahedron method to achieve better integration accuracy.

The density of states of a solid is a continuous function, but originates from a Brillouin zone integral,

g(ϵ)=nBZd3kδ(ϵϵn(𝐤)). (4.3)

The index n runs over all bands, and 𝐤 is the crystal momentum.

As it stands, this expression for the density of states creates a problem. We obtain it from a discrete set of k-points, rather than as a continuous integral. The density of states from a discrete k-mesh would thus be a sequence of δ-function peaks.

With the lines for band structure output, DOS output_tetrahedron in control.in, we can run the example calculation in the testcases directory.

Once this is done, a number of new files have appeared (in addition to the usual standard output stream): The files band1XXX.out, KS_DOS_total_tetrahedron.dat, and KS_DOS_total_raw_tetrahedron.dat.

Refer to caption
Refer to caption
Figure 4.8: Band structure and density of states of the bulk Si test case as generated by the aimsplot.py script, which is located in the “utilities” directory. A schematic visualization of the Brillouin zone of the fcc Bravais lattice including some high-symmetry k-points is also shown (taken from Wikipedia).

These files contain the necessary information in column formats that can be manipulated for further output with standard plotting tools. The result of such a manipulation is shown in Fig. 4.8.

However, it would be nice to not have to personally massage the output data every time to get this kind of output. The script ’aimsplot.py’, located in the utilities directory, should do this for you (type “aimsplot.py --help” for more information regarding available options). Simply run it in the output directory, and a version of the band structure and DOS should automatically appear as a png file.

Note: You will need python3 (not python2) to run aimsplot.py and you will need a compatible version of matplotlib (a common python library) installed, or else the aimsplot.py script will exit with a verbose but (to the uninitiated) cryptic error message.

As there are various conventions for plotting the density of states, each possible DOS plot produces two different output files. The first one, ending in “_raw.dat”, contains the exact density of states, using exactly the same eigenvalue energies and energy zero as internally in FHI-aims (in cluster systems, the energy zero is simply the vacuum level; in periodic systems, the energy zero is set by the 𝐆=0 component of the long-range electrostatic potential).

The file KS_DOS_total_tetrahedron.dat contains the same information, but the energy zero is shifted to the Fermi-level of the system. Note that, for metals, the Fermi level is usually well defined. For insulators, it is not well defined, and may end up anywhere in the gap (at some location where FHI-aims finds the normalization criterion for the number of electrons to be sufficiently fulfilled).

If, for a semiconductor, you wish to have the energy zero of a DOS shifted to the top of the valence band, you will have to do this as an extra step. Identify the location of the valence band edge (highest occupied level) either from the FHI-aims output or from the band structure, and shift the x-axis of KS_DOS_total_tetrahedron.dat and other files by the appropriate amount.

Except for the total density of states, there are also a few other options to obtain densities of states with some local, spatial resolution. In particular:

  • an atom-projected density of states, using output atom_proj_dos_tetrahedron

  • an angular-momentum component projected density of states averaged over all atoms of a single species, using the keyword output species_proj_dos_tetrahedron. This is useful because different atoms of the same element may be grouped together as a separate species in control.in, e.g., all atoms in the first layer of a slab. The species-projected density of states (but not currently the atom-projected version) will also be automatically plotted by aimsplot.py.

Note that there are two useful projections of the DOS, both of which are based on a Mulliken analysis. output species_proj_dos_tetrahedron can be used to get a species-dependent angular momentum projection, while output atom_proj_dos_tetrahedron resolves the states into the various atomic angular-momentum contributions.

As mentioned above, the species_proj_dos_tetrahedron keyword can also be used to get the DOS contribution from a whole group of atoms at once. Just duplicate the species defaults for the element(s) in question in control.in, give the species a new name, and then use that new species (say) for all surface atoms in geometry.in. This will get you the density of states contribution from all surface atoms.

Finally, an element-decomposed version of the band structure can also be obtained, using the outputband_mulliken keyword as an alternative. Two separate scripts for visualization are available in the utilities directory, band_mlk.py and band_mlk_soc.py, for non-SOC and SOC, respectively. Like aimsplot.py, these scripts are intended as starting points for plotting purposes and should be customized further, if necessary. The python matplotlib library is used to generate these plots and documentation of this package is available online.