Part 3: Bulk Properties with the HSE06 Hybrid Density Functional: Energy Band Structure and Density of States
Objective
We now have properly relaxed the atomic and electronic structure of hematite Fe\(_2\)O\(_3\) and can analyze the electronic structure and spin state in more detail. Specifically, we will focus on the electron band structure, local spin moment, partial and full densities of states as calculated from the DFTHSE06 effective singleparticle eigenstates, using "light" settings. As we noted before, a better set of species_defaults ("intermediate" or better) should be used for such a calculation.
However, before we launch into the discussion, we need to clarify a few fundamental issues. While we can calculate a band structure and density of states using the DFTHSE06 density functional, it is no a priori clear if and how this data will correspond to any experimentally measurable observable. We need to understand the physics of what we are doing before we can move on to the technical aspects of the calculation.
Current Experiment and Theory for the Fundamental Band Gap of Fe\(_2\)O\(_3\)
A recent (2019) paper by Piccinin (https://doi.org/10.1039/C8CP07132B) provides a summary of the state of knowledge for bulk hematite Fe\(_2\)O\(_3\), including state of the art theory and its connection to experimental data.
The crux is that there are several different experimental observables that we might like to simulate:

Neutral, including optical, electronic excitations of the system. A subset of these excitations will be measurable in optical spectroscopies, for instance, absorption  reflected most effectively by the color of the material. Piccinin (whose work is computational) settles for an experimental optical absorption onset of approximately 2 eV. The actual quoted value is somewhat higher (2.142.20 eV) but in this author's view, the raw experimental data (reproduced in Figure 10 of Piccinin in 2019) appear to leave some room for discussion. However, a critical observation is that the absorption onset corresponds to the socalled optical gap of the system, that is, to the lowestenergy optically active neutral excitations of the system. In DFTHSE06 (the subject of this tutorial), we do not approximate the optical gap. At best, we approximate the fundamental gap of the system (ionization potential minus electron affinity), see next point.

The fundamental gap of a molecule or solid is defined as the difference between the ionization potential (IP) and the electron affinity (EA) of the system  i.e., the difference between (definition of IP:) the lowest energy required to remove an electron from an initially neutral crystal and (definition of EA:) the highest energy that can be gained by attaching an electron to an initially electronically neutral crystal. Piccinin quotes an experimental value of approximately 2.6 eV for the experimentally determined, fundamental gap. (Table 1)
For reasons given below, one might hope to simulate an approximation to the fundamental gap of Fe\(_2\)O\(_3\) using the singleparticle eigenvalues of the HSE06 density functional. However, let us first assess what is the experimental situation regarding the fundamental gap of hematite Fe\(_2\)O\(_3\).

The experimental fundamental gap value of ~2.6 eV as quoted by Piccinin is taken from a 1999 paper, Zimmermann et al. https://doi.org/10.1088/09538984/11/7/002.

In Figure 15 of Zimmermann et al., photoemission data (IP) and inverse photoemission, also known as Bremsstrahlung isochromat spectroscopy (BIS) data are shown, substantiating (visually, without the deeper analysis presumably made in that work) an approximate difference of 23 eV between the IP and the EA. This data is, in this authors' view, as good as it gets when trying to estimate fundamental gaps from electron spectroscopies. One remaining issue is, however, that these measurements are surface measurements, obtained by charging an otherwise insulating singlecrystal sample. The impact of the surface structure after preparation of the single crystal in vacuum is not clear to this author, but it could be substantial, particularly for the experimentally derived EA.

The experimental fundamental gap of Fe\(_2\)O\(_3\) of Zimmermann et al. is further substantiated by a range of literature values, 2.02.7 eV, attributed to Ref. 80 of Zimmermann et al. This Ref. 80 is a 1986 photoemission study by Fujimori et al. https://doi.org/10.1103/PhysRevB.34.7318. However, crucially, the measured data in this paper pertain only to photoemission (i.e., to the IP), but no experimental data is provided for electron attachment (EA). Their quoted fundamental gap range is traced backwards to two temperaturedependent conductivity measurements, Ref. 31 of Fujimori et al.

Piccinin's own paper does report computational results for the fundamental gap of Fe\(_2\)O\(_3\) ontained by the GW approximation in various flavors. In short, "GW" is an approximation to the actual quasiparticle equation, i.e., it is based on an equation that incorporates the right physics for the quasiparticle gap, which we are hoping to approximate. While "GW" typically the "best one can do" for singleparticle like quasiparticle excitations in electronic structure theory of solids, the method does depend on several further detailed approximations, such as: (i) The underlying density functional used to calculate the Green's function and the screened Coulomb interaction, (ii) nonselfconsistent or selfconsistent computation of the GW selfenergy, (iii) the absence of a "vertex function", which covers the local correlation and which is already approximated to be unity within "GW". Each of these approximations could be critical and Piccinin shows that different, ostensibly choices for (i)(ii) can lead to a significant variation of GWpredicted fundamental gaps for Fe\(_2\)O\(_3\), between 1.44 eV and 5.05 eV. Piccinin eventually settles for a nonGW treatment as the empirically "best" approach: a parameterized hybrid density functional (PBEh with an exchange mixing parameter of 1/7).

The presence of delectron correlation in Fe adds some significant complexity to all of this, a discussion not pursued here, as does the absence of electronphonon renormalization (i.e., impact of nuclear motion and relaxation on the experimental quasiparticle properties).

Our spin state is the expected zero Kelvin spin state. It turns out that hematite has a spin phase transition at around 250 K, which could affect the observed experimental gap. Within scalar relativistic DFT, however, we cannot simulate this phase transition, since the roomtemperature spin arrangement is "canted", i.e., the local spin moments found on the Fe ions are not collinear. Scalar relativistic DFT unfortunately, only allows for directionless spin arrangements (called "up" or "down" but not coupled to the lattice). We would need a better treatment (spinorbit coupled) to obtain an idea of the canted spin arrangement.
This is the backdrop before which the band structure predictions attempted here must be considered.
Why DFTHSE06 for Band Structure Predictions of Fe\(_2\)O\(_3\)?
In general, DFT eigenvalues (including those of the HSE06 hybrid density functional, used below) do not have a formal justification  except as approximate sources of information on the ionization potential (electron removal energy) and the electron affinity (electron attachment). Both statements follow from Janak's theorem (which holds for exact density functional theory, not approximate density functionals).
In practice, an approximate density functional must have a mathematically discontinuous derivative as a function of particle number at an integer electron number in order to be able to predict both IP and EA in the form of highest occupied and lowest unoccupied eigenvalues of the same, neutral DFT calculation. The discontinuity is not the only condition (the functional should also be linear as a function of particle number, between different integer particle numbers, and it should have the right slope).
In brief, the HSE06 functional (and any hybrid density functional) has a derivative discontinuity at integer particle numbers by construction. This sets hybrid density functionals apart from semilocal functionals, such as PBE. Thus, we can in principle use an approximate hybrid density functional as a starting point of a qualitative prediction of the fundamental gap and perhaps of other properties of the electronic quasiparticle band structure. Piccinin's work ultimately follows this approach (using PBEh, not HSE06) and simultaneously shows that one must nevertheless pay attention to whether an actual density functional parameterization can be expected to yield reasonable values for the predicted fundamental gap.
Indeed, it is this authors' (VB) viewpoint that, for materials with a rough fundamental band gap range (1.53.0 eV) here attributed to Fe\(_2\)O\(_3\), an nonmaterials specific calculation by the HSE06 density functional with fixed parameters (screening parameter 0.11 (Bohr radii)\(^{1}\), exchange parameter 0.25) has a reasonable chance of predicting the fundamental gap within a range of several tenths of an eV.
In fact, in some more detailed work on sulfides and on halide perovskites in VB's group, the fundamental band gap is typically somewhat underestimated, not overestimated, by the HSE06 density functional with the above, fixed parameters (see, e.g., Refs. https://dx.doi.org/10.1021/acs.chemmater.8b03380, http://dx.doi.org/10.1021/acs.chemmater.7b02638, https://dx.doi.org/10.1021/acs.chemmater.8b03380, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.121.146401, https://pubs.acs.org/doi/abs/10.1021/jacs.9b02909, https://doi.org/10.1021/acs.chemmater.9b05107 and references therein for such examples). As one advantage, the choice of the HSE06 functional with fixed parameters avoids a materialsspecific fit. On the other hand, the presence of the transition metal Fe introduces aspects of the electronic structure (delectron correlation) that may not follow the trends observed in electronically simpler semiconductors and it is well known that the redox levels of transition metal ions can pose significantly larger problems.
Tasks
In order to calculate an approximate HSE06 energy band structure, here using "light" settings for the tutorial, we proceed as follows (see solutions
folder for illustrations of the expected output).
 Create a new folder
bands_and_dos
.  Preparation of
geometry.in
: Take thegeometry.in.next.step
from the final HSE06 relaxation of Part 2, copy it tobands_and_dos
and make sure to readd the proper spin initialization, which was not printed ingeometry.in.next.step
.  In the construction of
control.in
for the band structure and DOS calculation, several options are available, including editing at the command line, graphical setup in GIMS, or command line setup in CLIMS. We recommend to try the "Band Structure workflow" in GIMS.  We will need to choose a particular kspace path to plot the energy band structure, which we will set up using GIMS or CLIMS. Both utilities rely on the standardized kspace paths for different space groups as defined by Setyawan and Curtarolo (https://doi.org/10.1016/j.commatsci.2010.05.010). In other words, the choice of kspace path for the band structure calculation depends on the crystal having the right space group. However, if you check the symmetry and the Bravais lattice of the relaxed structure from Part 2, you will find that the geometry relaxation has slightly changed the structure and, thus, slightly altered the space group. To get back the space group of the hematite structure, we must ensure that GIMS and CLIMS recognize it correctly and discard small deviations. You can either do this using GIMS, by changing the
symmetry threshold
to1e2
in the settings (upper right corner) and, then, press the button Primitive. Alternatively, you can do it with clims by the following command:climsunitcellinfo primitive symthresh 1e2
.  Preparation of
control.in
: If you did not use GIMS, you may useclimspreparerun hse06 bands dos
. Again, to save computer time for this tutorial, we only use light species defaults. For publishing results we recommend to use intermediate settings at least.  As mentioned in Part 2 of this tutorial, if you want to just restart an earlier HSE06 calculation from its converged electronic structure, you could store density matrices for each kpoint in the earlier step (many of them for this small structure), then copy these stored density matrices over to the present working directory and add the appropriate restart keyword to
control.in
also here:elsi_restart read_and_write 100
. For larger hybrid DFT calculations with fewer kpoints, this strategy is decidedly useful. For the present example, a restart may not be essential.  After the calculation of band structure and density of states is complete, visualize the results either with
climsaimsplot
or upload them to GIMS.
As a consistency check, the control.in
format for this step may end up looking as follows:
xc hse06 0.11
hse_unit bohr
hybrid_xc_coeff 0.25
exx_band_structure_version 1
relativistic atomic_zora scalar
k_grid 10 10 10
output dos_tetrahedron 20 10 3001
output species_proj_dos_tetrahedron 20 10 3001
output band 0.00000 0.00000 0.00000 0.50000 0.00000 0.00000 23 G L
output band 0.50000 0.00000 0.00000 0.50000 0.23420 0.23420 18 L B1
output band 0.76580 0.50000 0.23420 0.50000 0.50000 0.50000 20 B Z
output band 0.50000 0.50000 0.50000 0.00000 0.00000 0.00000 21 Z G
output band 0.00000 0.00000 0.00000 0.36710 0.00000 0.36710 27 G X
output band 0.63290 0.36710 0.00000 0.50000 0.50000 0.00000 10 Q F
output band 0.50000 0.50000 0.00000 0.63290 0.63290 0.23420 8 F P1
output band 0.63290 0.63290 0.23420 0.50000 0.50000 0.50000 17 P1 Z
output band 0.50000 0.00000 0.00000 0.76580 0.36710 0.36710 14 L P
spin collinear
[....]
Important new keywords:

exx_band_structure_version 1
is a specific keyword that instructs the program to create the hybrid density functional theory energy band structure from a lattice Fourier transform based installation of the electronic structure determined using the realspace periodic cell ("Bornvon Karman cell") corresponding to thek_grid
used during the s.c.f. cycle. This is a process that works exactly it the realspacek_grid
is sufficiently dense. It is also competitively fast. However, in cases in which the s.c.f.k_grid
is too small, the lattice Fourier transform can introduce very significant aliasing artifacts into the band structure (one cannot miss these artifacts). In these cases, one could use theexx_band_structure_version
keyword to change to a different, but far more expensive reciprocalspace computation of the energy band structure. In practice, however, the much better way forward is to rerun the calculation with a better converged s.c.f.k_grid
choice. 
A single
output band [...]
line requests energy band structure output along a path between two points in reciprocal space. The full kspace path is set up by GIMS and can also be visualized there. 
output dos_tetrahedron 20 10 3001
writes the density of states, integrated using the tetrahedron method, in a particular energy interval: here, between 20 and +10 eV in internal energy units of FHIaims  typically, the Fermi level in periodic FHIaims calculations is a few eV away from this internal zero. The number of energy output points is 3001, i.e., one point for each 0.01 eV of the energy axis. Be sure to choose sufficiently many output energy points to obtain a reasonably finely resolved density of states. 
output species_proj_dos_tetrahedron 20 10 3001
creates a projected density of states, resolved into different species (typically, chemical elements) as defined incontrol.in
, as well as spin and angular momentum channels. The keyword follows the same logic as the density of states. However, additionally, this keyword will write an atomsinmolecules (Mulliken) decomposition of qualitative atomic partial charges and spin states in theaims.out
file  see there.
Results
The results for the corresponding band structure and DOS can be again found in the solutions
folder.
In addition to the band structure and DOS, we learn the following important predicted quantities (towards the end of aims.out
file, last s.c.f. iteration):
What follows are estimated values for band gap, HOMO, LUMO, etc.  They are estimated on a discrete kpoint grid and not necessarily exact.  For converged numbers, create a DOS and/or band structure plot on a denser kgrid.
Highest occupied state (VBM) at 8.57347524 eV (relative to internal zero)  Occupation number: 1.00000000  Kpoint: 276 at 0.300000 0.900000 0.700000 (in units of recip. lattice)  Spin channel: 1
Lowest unoccupied state (CBM) at 5.25801751 eV (relative to internal zero)  Occupation number: 0.00000000  Kpoint: 373 at 0.500000 0.200000 0.900000 (in units of recip. lattice)  Spin channel: 1
ESTIMATED overall HOMOLUMO gap: 3.31545773 eV between HOMO at kpoint 276 and LUMO at kpoint 373  This appears to be an indirect band gap.  Smallest direct gap : 3.32916001 eV for k_point 81 at 0.100000 0.300000 0.700000 (in units of recip. lattice)  Direct gap HOMO spin channel : 1  Direct gap LUMO spin channel : 2 The gap value is above 0.2 eV. Unless you are using a very sparse kpoint grid, this system is most likely an insulator or a semiconductor.
So, the s.c.f. kspace grid provides us with an estimate of the band gap of 3.33 eV  somewhat higher than the experimentally conjectured value, but, in this author's opinion, actually not unreasonably high. The reference by Piccinin shows that it will be hard to arrive at a single, firstprinciples predicted value for hematite Fe\(_2\)O\(_3\) with standard methods, without resorting to experimental data for calibration.
There is a second band gap estimate in control.in
, this one derived from only the kpoints on the requested band structure paths. If the valence band maximum (VBM) and conduction band minimum (CBM) were both part of the plotted band structure, this band gap value would be authoritative and the band gap would have to be smaller than the gap estimated on the s.c.f. kspace grid. However, in this case the band gap estimate from the band structure output is larger:
"Band gap" of total set of bands:
 Lowest unoccupied state: 5.24794573 eV
 Highest occupied state : 8.63662801 eV
 Energy difference : 3.38868228 eV
It would appear that our band structure did not include the kpoints at which the exact VBM and CBM are located. The estimate from the s.c.f. k_grid
is better.
The charge and spin analysis of the atoms in the structure looks as follows:
Performing Mulliken charge analysis on all atoms. Full analysis (per state, per kpoint, etc.) will NOT be written to separate file 'Mulliken.out'. This file can be requested by stating 'output mulliken' explicitly. Summary of the peratom charge analysis:   atom electrons charge l=0 l=1 l=2 l=3  1 24.628499 1.371501 6.153311 12.387560 6.031940 0.055687  2 24.628499 1.371501 6.153310 12.387560 6.031941 0.055687  3 24.628499 1.371501 6.153311 12.387560 6.031941 0.055687  4 24.628499 1.371501 6.153311 12.387560 6.031941 0.055687  5 8.914810 0.914810 3.853794 5.049802 0.011214  6 8.914096 0.914096 3.853512 5.049376 0.011208  7 8.914097 0.914097 3.853512 5.049377 0.011208  8 8.914810 0.914810 3.853794 5.049802 0.011214  9 8.914096 0.914096 3.853512 5.049377 0.011208  10 8.914096 0.914096 3.853512 5.049377 0.011208   Total 152.000000 0.000000
Summary of the peratom spin analysis:   atom spin l=0 l=1 l=2 l=3  1 4.208603 0.017182 0.031514 4.161340 0.001433  2 4.208603 0.017182 0.031514 4.161340 0.001433  3 4.208602 0.017181 0.031514 4.161340 0.001433  4 4.208602 0.017181 0.031514 4.161340 0.001433  5 0.000000 0.000000 0.000000 0.000000  6 0.000031 0.000004 0.000027 0.000000  7 0.000031 0.000004 0.000027 0.000000  8 0.000000 0.000000 0.000000 0.000000  9 0.000031 0.000004 0.000028 0.000000  10 0.000031 0.000004 0.000027 0.000000   Total 0.000000
Atomic charges from a Mulliken analysis are not expected to be equal to chemical charges. The fact that the sign and trend is right is good.
However, the spin state and spin moment is remarkably close to literature values (4.21, denoting the difference between spinup and spindown electron count attributed to a given atom). We did retain the accepted spin state.
Finally, here are the calculated energy band structure and densities of states. By the way, the FHIaims band structure and DOS output files are all column formatted ASCII files, so the data can also be read into any 2D plotting program for further analysis there.
The choice of the energy zero in the band structure plot for a system with a gap is arbitrary, by the way. It has no physical significance  all it does is place the zero somewhere in the gap (between occupied and unoccupied states). Specifically, the choice of this plotted energy zero between different structures with a gap must NEVER be interpreted as a relative shift of band structures. For proper aligmnent of the energy levels of different materials (if that can be done), it is necessary to identify and analyze a well defined internal reference level of the structure.