Skip to content

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 DFT-HSE06 effective single-particle 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 DFT-HSE06 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.14-2.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 so-called optical gap of the system, that is, to the lowest-energy optically active neutral excitations of the system. In DFT-HSE06 (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 single-particle 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/0953-8984/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 2-3 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 single-crystal 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.0-2.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 temperature-dependent 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 single-particle 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) non-selfconsistent or self-consistent computation of the GW self-energy, (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 GW-predicted fundamental gaps for Fe\(_2\)O\(_3\), between 1.44 eV and 5.05 eV. Piccinin eventually settles for a non-GW treatment as the empirically "best" approach: a parameterized hybrid density functional (PBEh with an exchange mixing parameter of 1/7).

  • The presence of d-electron correlation in Fe adds some significant complexity to all of this, a discussion not pursued here, as does the absence of electron-phonon 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 room-temperature 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 direction-less spin arrangements (called "up" or "down" but not coupled to the lattice). We would need a better treatment (spin-orbit 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 DFT-HSE06 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.5-3.0 eV) here attributed to Fe\(_2\)O\(_3\), an non-materials 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 materials-specific fit. On the other hand, the presence of the transition metal Fe introduces aspects of the electronic structure (d-electron 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).

  1. Create a new folder bands_and_dos.
  2. Preparation of geometry.in: Take the geometry.in.next.step from the final HSE06 relaxation of Part 2, copy it to bands_and_dos and make sure to re-add the proper spin initialization, which was not printed in geometry.in.next.step.
  3. 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.
  4. We will need to choose a particular k-space path to plot the energy band structure, which we will set up using GIMS or CLIMS. Both utilities rely on the standardized k-space 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 k-space 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 to 1e-2 in the settings (upper right corner) and, then, press the button Primitive. Alternatively, you can do it with clims by the following command: clims-unit-cell-info --primitive --symthresh 1e-2.
  5. Preparation of control.in: If you did not use GIMS, you may use clims-prepare-run --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.
  6. 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 k-point 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 k-points, this strategy is decidedly useful. For the present example, a restart may not be essential.
  7. After the calculation of band structure and density of states is complete, visualize the results either with clims-aimsplot 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 real-space periodic cell ("Born-von Karman cell") corresponding to the k_grid used during the s.c.f. cycle. This is a process that works exactly it the real-space k_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 the exx_band_structure_version keyword to change to a different, but far more expensive reciprocal-space 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 k-space 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 FHI-aims - typically, the Fermi level in periodic FHI-aims 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 in control.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 atoms-in-molecules (Mulliken) decomposition of qualitative atomic partial charges and spin states in the aims.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 k-point grid and not necessarily exact. | For converged numbers, create a DOS and/or band structure plot on a denser k-grid.

Highest occupied state (VBM) at -8.57347524 eV (relative to internal zero) | Occupation number: 1.00000000 | K-point: 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 | K-point: 373 at 0.500000 0.200000 0.900000 (in units of recip. lattice) | Spin channel: 1

ESTIMATED overall HOMO-LUMO gap: 3.31545773 eV between HOMO at k-point 276 and LUMO at k-point 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 k-point grid, this system is most likely an insulator or a semiconductor.

So, the s.c.f. k-space 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, first-principles 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 k-points 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. k-space 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 k-points 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 k-point, etc.) will NOT be written to separate file 'Mulliken.out'. This file can be requested by stating 'output mulliken' explicitly. Summary of the per-atom 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 per-atom 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 spin-up and spin-down 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 FHI-aims 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.