Chapter 6.22. WIEN2k: an augmented plane wave plus local orbital package for the electronic structure of solids

An efficient and accurate scheme for solving the many-electron problem of a crystal (with nuclei at fixed positions) is the Kohn–Sham theory within DFT (Hohenberg & Kohn, 1964; Kohn & Sham, 1995). The key quantity therein is the spin density ρ_σ(r), in terms of which the total energy is is the kinetic energy (of non-interacting particles), and E_NN are the repulsive electron–electron and nuclear–nuclear Coulomb energies, respectively, is the nuclear–electron attraction and is the exchange–correlation energy. While the equation above is formally still exact and we can nowadays calculate almost all of these terms virtually exactly by numerical methods, we do not know the exact form of . This exchange–correlation energy should account for the approximations made in the kinetic energy (of non-interacting particles) and cancel the self-interaction of an electron with itself in the classical term. The variation of E_tot with respect to the density leads to the Kohn–Sham equations, which must be solved self-consistently in an iterative process,where are single-particle orbitals which are used to define the total density ρ_σ = . Note that the corresponding energies are formally just Lagrange parameters. Many different methods exist to solve the Kohn–Sham equations using different basis sets and under various approximations. One of the most accurate schemes is based on the augmented plane-wave (APW) method, as implemented in the WIEN2k code and discussed in the following section. Besides the accurate solution of the Kohn–Sham equations, the approximations for and determine the quality of the results. Perdew & Schmidt (2001) have introduced a `Jacob's Ladder into DFT heaven', which classifies the various approximations according to accuracy but also complexity.

WIEN2k natively supports the local density approximation (LDA), several common versions of the generalized gradient approximation (GGA; see, for example, Perdew et al., 1996, 2008) and also various meta-GGAs (Sun et al., 2015), into which the kinetic energy density t(r) also enters. In addition, there is an interface to the libxc library (Marques et al., 2012), which includes an even larger variety of different functionals. For insulators, various versions of hybrid functionals can be used, which mix a certain fraction of Hartree–Fock (HF) exchange with DFT. These functionals are computationally rather expensive but became quite popular, in particular in quantum chemistry (B3LYP; Becke, 1993) and for insulating solids (HSE06; Krukau et al., 2006; Tran & Blaha, 2011). However, because of limited correlation they fail badly for metals, predicting palladium or platinum to be ferromagnetic metals (Tran et al., 2012).

Nowadays, most systems are simulated by GGA calculations with acceptable but not perfect accuracy. Exceptions are systems with van der Waals bonds, where no bonding at all may occur, and correlated 3d or 4f electrons, where GGAs may break down and give an incorrect ground state (metal instead of insulator, nonmagnetic instead of magnetic solution). For these systems, one often applies a semi-empirical DFT+U approximation (Anisimov et al., 1993) or the `exact exchange for correlated electrons' (EECE) method (Novák et al., 2006), where a certain fraction of HF exchange is applied selectively to the correlated 3d or 4f electrons. These very cheap approximations are used a lot for 3d transition-metal or lanthanide compounds. For van der Waals systems one can use the latest meta-GGAs (Sun et al., 2015) or explicitly add long-range interactions as in DFT-D3 (Grimme et al., 2010), DFT-D4 (Caldeweyher et al., 2020) or van der Waals-DFT (Tran et al., 2017, 2019). If the energy gap is the main quantity of interest, one can use the modified Becke–Johnson (mBJ) potential (Tran & Blaha, 2009), which predicts gaps with high accuracy at virtually no cost. Alternatively, WIEN2k also serves as a basis for many-body perturbation theory, such as GW (Jiang et al., 2013; Jiang & Blaha, 2016), BSE (Laskowski & Blaha, 2010) or dynamical mean-field theory (DMFT; Held, 2007; Kuneš et al., 2010; Haule, 2017).

Even within a given crystal structure, structural parameters such as lattice parameters or the free internal coordinates of some atomic sites (Wyckoff positions) may change the calculated properties (such as an X-ray absorption spectrum) significantly. While for most crystalline materials the experimental lattice parameters have been determined with high precision (usually much more accurately than theoretical equilibrium lattice parameters, which often have errors of 1–2%), the exact atomic positions may not be known accurately in cases where single crystals are not available or light elements are placed next to heavy ones. In particular, for `non-ideal' solids, for example surfaces, interfaces, alloys, impurities or vacancies, the exact atomic positions are usually not known experimentally and one must optimize them theoretically. This can be performed in a two-step procedure: one can simply use the accurate experimental lattice parameters or optimize the lattice parameters theoretically via the total energy, i.e. varying them until the total energy reaches a minimum. In any case, for a given set of lattice parameters one should always optimize the atomic positions. We utilize the calculated forces acting on the atoms and an efficient algorithm (Marks, 2013) to simultaneously find the equilibrium positions and the self-consistent electron density, i.e. we move the atoms during the self-consistent field (SCF) cycles until the forces are close to zero and the density is self-consistent.

In particular, it should be noted that lattice parameters depend crucially on the selected DFT functional. As a rule of thumb, one can expect that the LDA will overbind systems (implying lattice parameters that are too short) except for heavy elements (5d series), while PBE (Perdew et al., 1996) will underbind such heavy elements (implying lattice parameters that are too long) but gives pretty good parameters for compounds with atoms from the 3d series. Functionals such as PBESOL (Perdew et al., 2008) are often a good compromise (Tran et al., 2016). Once the structure has been relaxed, one can check distances and angles and compare the structural parameters with available experiments.

For the fully relaxed structure one can calculate phonons via a frozen-phonon approach using tools such as PHONON (Parlinski et al., 1997) or PHONOPY (Togo & Tanaka, 2015), in which long-range Coulomb effects can also be taken into account via Born effective charges calculated by BerryPI (Ahmed et al., 2013).

Once the structure has been relaxed, the next step is to obtain the ground-state electronic structure. WIEN2k allows the energy band structure, densities of states (DOS) and electron densities to be calculated in a convenient way. By a decomposition of the total DOS into partial DOS according to the contributions of different atoms and angular momenta, one can nicely explore the detailed electronic structure and unravel the chemical bonding. These partial DOS also serve as a first guess for the interpretation of XAS spectra. Key quantities for many materials such as the energy band gap or their optical properties can be calculated with good predictive power due to the TB-mBJ method (Tran & Blaha, 2009). Electron densities of the valence states, or sometimes of a specific energy region only, let one visualize the electronic structure in a more intuitive way and directly reveal ionic or covalent interactions of bonding or antibonding nature.

Since WIEN2k is an all-electron program, it is very well suited to calculate various hyperfine interaction parameters for NMR and Mössbauer spectroscopy. The electric field gradient (quadrupole interactions), as well as the Mössbauer isomer shift and hyperfine fields, can be obtained directly from the ground-state electron density near or at the nucleus (Blaha, 2010; Body et al., 2007). Much more effort is required for the calculation of NMR chemical shifts because the application of an external magnetic field breaks translational symmetry. A perturbative approach has been implemented in WIEN2k (Laskowski & Blaha, 2014, 2015a), and recently we have also extended these calculations to Knight shifts in metallic systems (Laskowski & Blaha, 2015b; Kalantari et al., 2017).

Both near-edge XAS (NEXAFS or XANES) and electron energy-loss (EELS) spectroscopy (Hébert, 2007) represent a process in which a core electron of one particular atom is excited by a photon (or a fast electron) to an empty conduction-band state. When we want to simulate such processes, we usually rely on eigenvalues ɛ_n,k obtained from DFT band-structure calculations, neglecting the fact that this is not justified in principle. In order to calculate the spectrum, we need matrix elements between the initial and the final state including the electromagnetic radiation. This is typically performed within the dipole approximation and the intensity can be calculated from the squared dipole (momentum) matrix elements between the initial (core, I) and final (conduction band, F) state,

This leads to the well known dipole-selection rule, and XAS reveals transitions from core states on atom X with angular momentum l into unoccupied orbitals with Δl ± 1 on the same atom. In essence, the spectrum is calculated from the corresponding partial DOS times the squared radial matrix element. In the case of polarized light and oriented samples one can further fine-tune the information and obtain orientation-dependent spectra, which can be understood by substituting the l-like partial DOS with an appropriate lm-like DOS; for example, for a K-edge spectrum of a hexagonal system one might replace the total p-DOS by a p_z and p_x + p_y DOS.

Using such a scheme usually leads to quite poor agreement with experiment because the excited electron and the core hole will interact with each other, leading to strong excitonic effects and, according to the final-state rule (von Barth & Grossmann, 1979), the electronic structure of the final state determines the spectrum. To model these effects one should produce a supercell (as large as possible, but depending on the specific system and your computer) of about 32–128 atoms, and on one of the atoms remove one of the core electrons and add it to the valence electrons (or as a constant background charge) to keep the system neutral. In the following SCF cycle the valence orbitals on the atom with the core hole will be attracted much more by the less screened nucleus than orbitals on other atoms and will be lowered in energy. In addition, the supercell allows also some (static) screening by electrons of the surrounding atoms. The resulting XAS spectrum will indicate much stronger localization (at low energies) and in general should agree significantly better with experiment. However, the screening in such static supercell calculations will not be perfect and exchange effects are treated only at a GGA level. Thus, in particular for metals (Luitz et al., 2001) or small-gap semiconductors (Mizoguchi et al., 2004) one sometimes uses a smaller core hole (a `half core hole') to obtain better agreement with experiment. Since energy band gaps are often very badly described by standard GGA, one can use the TB-mBJ (Tran & Blaha, 2009) potential, which can give band gaps close to experiment and improve the XANES simulation considerably in certain cases (Hetaba et al., 2012). In particular, the K edges of light elements are very well reproduced and explained by such simulations. Recent applications using WIEN2k include C K spectra in fullerenes (Erbahar et al., 2016), O K and Nb M₃ edges in Nb₃O₇(OH) (Khan et al., 2016) and B and N edges in h-BN monolayers on Ni, Rh and Pt(111) surfaces (Laskowski et al., 2009). On the other hand, the single-particle approach fails to some extent for the L_2,3 edges of 3d transition-metal (TM) compounds (Laskowski & Blaha, 2010). For the early TMs the branching ratio between the L₂ and L₃ spectra is not 1:2, as would be expected naively from the occupation numbers of the 2p_1/2 and 2p_3/2 states, while for the later TMs strong correlation effects are not fully treated in the GGA. A partial remedy is to treat the excitonic effects in a more rigorous way using the Bethe–Salpether approach (BSE) as described below (Laskowski & Blaha, 2010).

4.4.1. How to obtain XAS spectra using WIEN2k: a tutorial

In the following, a B K-edge calculation using WIEN2k is described for hexagonal BN. As input for such a calculation one has to provide the crystal structure (space group, lattice parameters and atomic positions): for hexagonal BN, space group P63/mmc (No. 193); a = 2.504, c = 6.661 Å; B (1/3, 2/3, 3/4), N (1/3, 2/3, 1/4).

Basically, three commands suffice to obtain a self-consistent ground-state calculation.

(i) makestruct_lapw: generates the structure. Enter the data above (or import a cif file).

(ii) init_lapw -b: use the default options to generate all other input files.

(iii) run_lapw: run SCF cycle with defaults.

As mentioned above, for a core-hole calculation one should create a supercell using the following.

(i) x supercell: for example a 4 × 4 × 1 or 5 × 5 × 2 supercell and mark one of the B atoms as `B 1'.

(ii) init_lapw: accept the symmetry changes made by the nn and sgroup programs.

(iii) Edit the input files for core and valence, removing one (half) of the B 1s electrons and adding it to the valence.

(iv) run_lapw: after the SCF cycle remove the extra valence electron from the input.

(v) x lapw2 -qtl: calculate partial charges.

(vi) x xspec: calculate the spectrum after specification of the B K edge in the input.

(vii) x broadening: apply Gaussian, Lorentzian and E-dependent broadening.

The obtained spectra are shown in Fig. 1, where the theoretical and experimental spectra have been aligned at the first peak. Obviously, without the core hole the first peak in the theoretical spectrum is far too broad and is missing the clear sign of excitonic effects. Introduction of a full (half) core hole sharpens this feature and the rather good agreement between experiment and theory is obvious. Closer inspection, however, reveals that the theory has problems with the exact peak position at higher energies (203 and 215 eV) and also that the shape of the peak at 198 eV does not fully match the experiment. The experimental double-peak structure at 198 eV can also be explained when one takes electron–phonon inter­actions and zero-point vibrations into account (Karsai et al., 2018). It should be mentioned that broadening has a large influence on the theoretical spectrum and one should be careful when trying to match a particular experiment, because broadening that is too large may obscure some features that may show up in an improved resolution experiment in the future.

Figure 1 Experimental B K edge of hexagonal BN (taken from Franke et al., 1997) together with theoretical simulations. The theoretical curves are broadened by Gaussian and Lorentzian (0.1 eV) and a linear E-dependent broadening (γ = ΔE/10).

Measurements on well defined surfaces or single crystals allows orientation-dependent spectra to be obtained, which can be interpreted employing the B-p_z and B-p_x,y partial DOS in the theoretical simulations (Fig. 2). Obviously, the first peak stems solely from B-π* bands, while the antibonding σ* bands contribute at higher energies.

Figure 2 Angle-dependence of the theoretical B K spectra (full core hole).

An analysis using the B-p partial DOS for the ground-state and core-hole calculations shows that the B-p_z DOS of the antibonding p* states at the bottom of the conduction bands form a very narrow, excitonic-like peak within the gap upon introduction of the core hole (Fig. 3). All other features in the partial DOS of the atom with a core hole also sharpen and move down in energy owing to the reduced screening of the nucleus by the core electrons.

Figure 3 B-p DOS in h-BN in the ground state and with a full core hole.