WIEN2k: An APW+lo program for calculating the properties of solids

Quantum mechanical calculations play a central role in understanding the properties of materials and, increasingly, predicting the properties of new materials. While in the early days, the emphasis was mainly on understanding the energy, atom positions, and band structure, modern codes now calculate a large number of different properties ranging from piezoelectric response to nuclear magnetic resonance (NMR) shielding, examples of which will be given later. With the advent of increasingly sophisticated methods and the ever increasing speed of computers over the last decades, in some cases the accuracy of quantum mechanical calculations rivals or even surpasses the accuracy of experimental measurements.

There are many different methods of theoretically modeling the behavior of electrons and atoms in materials. While earlier approaches focused on dealing with the electrons via wave functions,¹ many current methods use density functional theory (DFT),² which has significant speed advantages. Following the method outlined by Kohn and Sham³ (KS), the interacting many-body system of electrons is mapped onto a non-interacting system of quasiparticles, characterized by KS orbitals with a specific KS energy. They have many of the properties of the true electron wave functions and of particular importance is that one can fill up these KS orbitals as a function of their KS energy yielding the true electron density. The KS approach needs an exchange-correlation (XC) functional and the corresponding XC potential. However, the exact functional is unknown and approximations are needed (see Sec. II B).

A second split in terms of methods is how the atomic positions are considered, and there are two main methods: cluster calculations for a finite number of atoms, which focus on the local properties of some atomic arrangement, and those which are designed to exploit the periodic nature of most solids; the WIEN2k code is an example of the latter. We represent the solid by a unit cell, which is repeated in all three directions, corresponding to periodic boundary conditions. This assumes that the solid is perfect, ordered, and infinite; however, a real crystal differs from this ideal situation, since it is finite, may contain defects or impurities, and may deviate from its ideal stoichiometry. For these important aspects and how to handle them using supercells, see Chap. 8.2 of Ref. 4.

There are many computational methods for solving the KS equations, for instance, linear combination of atomic orbitals (LCAO), numerical basis sets, pseudopotential schemes, or space partitioning methods. A recent comparison⁵ of these methods showed that especially all-electron codes predict essentially identical results, demonstrating a high reproducibility, whereas some pseudopotential codes lead to large deviations. One of the most accurate codes is our WIEN2k code,⁶ which is the focus of this paper and is based on the augmented plane wave (APW) method. Detailed descriptions including many conceptual and mathematical details are given in Ref. 7. The term all-electron (see Chap. 8.4 of Ref. 4) means that all electrons from the core (starting from the 1s shell) to the valence states are included.

Turning to some historical specifics of our approach (see also Ref. 8), Slater⁹ proposed the original APW method. Unfortunately, the original formulation leads to a nonlinear eigenvalue problem due to the energy-dependent radial basis functions, which is computationally expensive. An important improvement came from Andersen,¹⁰ who introduced a linearization of this energy dependency, and Koelling and Arbman¹¹ made the linearized-APW (LAPW) method a practical computational scheme using the muffin-tin (MT) approximation (see Sec. II). This was taken a step further by Freeman and collaborators who made the LAPW method a full-potential all-electron total energy method.^12,13

This LAPW method formed the basis for the original WIEN code.¹⁴ However, the LAPW method had the drawback that only one principal quantum number per angular momentum ℓ could be described and thus failed to give reliable results for all elements on the left of the periodic table because these atoms require a proper description of shallow core states (semi-core) and valence states at the same time (e.g., 1s and 2s in Li or 3sp and 4sp in Ti). This problem was solved by Singh,¹⁵ who introduced local orbitals (LOs) for the description of semi-core states. He also noted that the LAPW method needed a larger plane-wave basis set than the APW method. To overcome this problem, he suggested the augmented plane wave plus local orbitals (APW+lo) method,^16,17 where the linearization of the energy-dependent radial wave function was facilitated by an extra local orbital (lo, different from an LO, see Sec. II), which has a superior plane-wave convergence compared to LAPW. Last but not least, the linearization of the energy dependency can introduce some inaccuracy in high precision calculations. This problem was finally solved by introducing additional higher (second) derivative LOs (HDLOs).^18,19 These latest developments form the basis of the present WIEN2k_19 code,⁶ while previous versions have been described in several reviews.^7,8,20–23 The method of our choice could be named (L)APW+lo+LO+HDLO, but we use a shorter acronym APW+lo. It is described in detail in Sec. II A.

In the APW-based methods, the unit cell is decomposed into spheres centered at the nuclear sites and an interstitial region,⁷ as shown in Fig. 1. These atomic spheres with radii R_MT must not overlap, but should be chosen for computational efficiency as large as possible with the additional constraint that R_MT for d-elements should be chosen to be about 10%–20% bigger than for sp-elements, while f-elements should get even larger spheres because for identical sphere sizes the number of plane-waves (PWs) to reach convergence is largest for the localized 4f (5f) electrons, medium for 3d (4d, 5d)-electrons, and much smaller for sp-states. An exception is the H atom, whose sphere with short C–H or O–H bonds should be chosen approximately half the size of R_MT(C) or R_MT(O). In WIEN2k, these sphere radii can be set automatically in an optimal way using the setrmt utility. Note that non-optimal sphere sizes may lead to poor convergence (eventually only for one particular atom) and significantly longer computing time or suffer from truncated Fourier or spherical-harmonic expansions. In the worst case, they can even produce “ghost-states” (unphysical eigenvalues in the occupied spectrum) if the R_MT of an sp-element is much larger than that of the other atoms.

The electron density ρ and KS potential v^KS (defined in Sec. II B) are expanded as a Fourier series in the interstitial (I) region (K denotes a reciprocal lattice vector in units of inverse bohr) and as lattice harmonics (symmetry adapted combinations of spherical harmonics Z_LM)²⁴ times radial functions ρ_LM(r) [$vLMKS(r)$ for the potential] inside the spheres (S_t, where t is the atom index),

By default, the Fourier expansion runs up to |K| = 12 for large R_MT (16 for R_MT < 1.2; 20 for R_MT < 0.7 bohr), while the angular momentum expansion truncates at L = 6. Note that the old “MT” approximation uses a constant value in the interstitial (i.e., only K = 0) and a spherically symmetric density/potential inside the spheres (i.e., only L = 0).

This space decomposition plays a crucial role in the definition of core and valence electrons, which are treated differently in APW-type methods. Core states are defined as having wave functions (densities) completely confined inside the atomic spheres. Thus, we do not use the standard definitions of core and valence, but, e.g., in 3d transition metals (TMs), the 3s and 3p states are also considered as valence, since a couple of percent of their charge leaks out of the atomic sphere. To distinguish them from the conventional definition, we call them semi-core states. Typically, these states are treated using LOs (see below) and their energies are less than 6 Ry below the Fermi energy, but in special situations (small spheres due to short nearest neighbor distances or high pressure), even lower lying states (such as Al-2p) have to be included. Relativistic effects are important for the core states, and thus, they are calculated by numerically solving the radial Dirac equation in the spherical symmetric part of the potential v^KS. Core states are constrained to be localized and not hybridized with states at the neighboring atoms, but we use a thawed core (no frozen core approximation), i.e., the core states are recalculated in each self-consistent field cycle.⁷ The semi-core and valence electrons are commonly treated scalar relativistically, i.e., including mass velocity and Darwin s-shift corrections, but neglecting spin–orbit (SO) interactions.^7,25 The SO effects can later on be included in a second variational step using the scalar-relativistic orbitals as a basis.^7,26 Since p_1/2 radial wave functions differ considerably from scalar relativistic (or p_3/2) orbitals, one can also enrich the basis set with additional p_1/2 local orbitals, specifically, an LO (see below) with a p_1/2 radial wave function, which is added in the second-variational SO calculation.²⁷ 

The basis functions for the valence electrons consist of APWs, which are plane waves in the interstitial region augmented with radial wave functions u_tℓ(r, E_tℓ) defined at a fixed energy E_tℓ, and lo.^7,16,17 An APW is given by

where k is a point in the first Brillouin zone (BZ), $Yℓm(r^)$ are spherical harmonics, and u_tℓ are solutions of the scalar-relativistic radial KS equation⁷ inside the sphere S_t. Note that these radial functions u_tℓ are recalculated in each self-consistent-field (SCF) cycle, allowing for an expansion/contraction corresponding to the given charge state (ionicity) of the atom. These adaptive basis functions are part of the reason for the high accuracy of APW-based methods. The coefficients $Atℓmk+K$ are chosen such that the interstitial and sphere parts of the APW match at the sphere boundary. However, these APWs allow no variations of the radial functions for eigenvalues different than E_tℓ and thus would be a poor basis. To overcome this constraint, the energy dependency is handled by a lo, which is nonzero only inside a MT sphere, and given by

where $u̇tℓ$ is the first energy derivative of u_tℓ. The coefficients $Atℓmlo$ and $Btℓmlo$ are chosen such that $ϕtℓmlo$ is zero at R_MT and normalized.

The APW+lo basis set has the advantage of a superior PW convergence as compared to the standard LAPW method,^7,17 reducing the number of PWs by almost 50%, but it needs additional lo basis functions. Thus, in WIEN2k, the default is to restrict the ℓ_max to the chemical ℓ-values sp(d, f), which are hard to converge, but use a standard LAPW basis set inside the spheres for the higher ℓ values (by default up to ℓ_max = 10),

where the coefficients $Atℓmk+K$ and $Btℓmk+K$ are chosen such that $ϕk+KLAPW$ and its first derivative are continuous at the sphere boundary.

As mentioned before, semi-core states (or also high-lying empty states) cannot be described accurately by APW+lo. For these states, the basis set has to be improved, and this can be done by adding another type of local orbitals, the LOs, containing radial functions u_tℓ calculated at the appropriate (e.g., semi-core) energy $EtℓLO,i$⁠,

For instance, for TiO₂, one would use the Ti-3p energy to calculate $uTi,1(r,ETi,1LO,i)$ and then add some Ti-4p radial function u_Ti,1(r, E_Ti,1), choosing the coefficients $AtℓmLO,i$ and $CtℓmLO,i$ such that the LO is zero at the R_MT and is normalized. By adding such LOs (representing Ti-3p states), a consistent and accurate description of both the Ti-3p semi-core and Ti-4p valence states is possible, retaining orthogonality, which is not assured when the multiple-window approach is used (see Ref. 7). Cases where this improvement is essential is the electric field gradient (EFG) calculation of rutile TiO₂²⁸ or lattice parameter calculations of compounds with such elements. Note how lo [Eq. (4)] and LO [Eq. (6)] differ in their respective second terms.

A clever choice of energy parameters E_tℓ in Eqs. (3)–(6) is essential for accurate results, and WIEN2k has several automatic ways to make an optimal choice in most cases.⁶ E_tℓ of semi-core states (actually, of all states whose energy in the free atom is more than 0.5 Ry below the highest occupied atomic orbital, e.g., also C-2s or Ar-3s states) are determined by taking the average of the two energies E_bottom and E_top, where the corresponding u_tℓ(R_MT) is zero or has zero slope. For localized d or f valence electrons, the same procedure is used, but E_top is searched only 0.5 Ry above E_F to ensure that the energy parameters are set below E_F. The energy parameters of all other valence states are set to 0.2 Ry below E_F (0.2 Ry above E_F if there is a high lying semi-core LO). Thus, all our energy parameters are dynamically updated during the SCF cycle and not fixed by input.

Implicit in this approximation is a linearization of the energy dependency of the radial wave functions. Since the true u_tℓ(r, E_tℓ = ε_i) varies most for more localized states (e.g., 3d or 4f) at larger distances from the nucleus, large spheres and a large valence bandwidth could cause a poor description of the variations of u_tℓ(r, ε_i) with energy, leading to a significant dependency of the results on R_MT (where smaller R_MT yield more correct results but with a larger computational effort). This can be solved by adding an LO, which involves the second energy derivative of u_tℓ, called an HDLO,^18,19

Figure 2 illustrates the effect of adding HDLOs for the lattice parameter of fcc-La as a function of R_MT. The lattice parameter becomes independent of R_MT when both d and f-HDLOs are added to the basis set, while the standard APW+lo+LO basis (LOs for 5s, 5p states) produces an error of 0.04 bohr for the largest R_MT.

The KS orbitals are expanded using the combined basis set described above (n is the band index),

and the coefficients $cnki$ are determined by the Rayleigh–Ritz variational principle. The number of APW (or LAPW) basis functions Eq. (3) [or Eq. (5)] is determined by the cutoff value K_max for the reciprocal lattice vectors K such that $k+K≤Kmax$ and depends on the smallest of the atomic radii $RMTmin$ and the type of atom. Typically, the necessary $RMTminKmax$ values range from 3 (for small H-spheres) to 7 for sp-elements, 8 for TM-d elements, and 9 for 4f lanthanides. These values can be reduced by 0.5–1 for low quality screening calculations and increased by 0.5–2 for highest precision. It should be mentioned that the efficiency of the APW+lo method depends crucially on the possible R_MT values. For instance, the O-2p states converge well with $RMTOKmax$ = 7. In MgO, one can use $RMTO$ = 2 bohrs, leading to a very small PW cutoff energy of 170 eV. However, in Mg(OH)₂, one has to use $RMTO$ = 1.1 bohrs due to the short O–H distances, leading to a PW cutoff of 550 eV, i.e., an order of magnitude larger effort.

The specific setup of all these basis functions can, of course, be selected manually by experts, but one of the great strengths of WIEN2k is that the default input usually works quite well and is fairly robust.

In DFT, the total energy of the system is given by³ 

The terms on the right-hand side represent the noninteracting kinetic, electron–nucleus, Hartree, XC, and nucleus–nucleus energies, respectively. The variational principle leads to the KS [or generalized KS²⁹ (gKS)] equations (in this section, the orbital index i is a shorthand notation for valence and core orbitals),

where v^KS is the KS potential,

which is the sum of the electron–nucleus, Hartree, and XC potentials. Choosing an appropriate functional E_xc in Eq. (9) [and potential v_xc in Eq. (11)] for the XC term is crucial in order to obtain reliable results for the problem at hand.^30–32 Several hundred^33,34 different functionals are available in the literature; some of them were proposed as general-purpose functionals, while others were devised for a specific property (e.g., bandgap) or types of systems (e.g., van der Waals). Numerous functionals have been implemented in the WIEN2k code, and below, we provide a brief overview of the different families of functionals. Note that for XC functionals that depend explicitly on the electron density ρ, e.g., the local density approximation (LDA) or the generalized gradient approximation (GGA), the XC potential v_xc is multiplicative, while for functionals that depend implicitly on ρ, e.g., meta-GGA (MGGA) or hybrids, the XC potential is non-multiplicative when implemented in the gKS scheme.

The LDA, GGA, and MGGA represent the first three rungs of Jacob’s ladder of XC functionals.³⁵ These approximations are semilocal, since E_xc is defined as

and the XC energy density ε_xc depends only locally on some properties of the system. In the LDA, ε_xc depends on the electron density $ρ=∑i=1Nψi2$⁠, while in the GGA, ε_xc depends also on the first derivative ∇ρ. At the MGGA level, the functionals depend additionally on the Laplacian of the electron density ∇²ρ and/or the kinetic-energy density $t=1/2∑i=1N∇ψi*⋅∇ψi$⁠. Semilocal functionals are the most commonly used methods in the solid-state community for the calculation of properties depending on the total energy such as the geometry, cohesive energy, or the adsorption energy of a molecule on a surface. The main reason is that they are faster than all other types of approximations and therefore allow calculations of larger systems.

There is a huge literature on the performance of semilocal functionals, concerning the geometry and cohesive energy of solids. Extensive benchmark studies have been conducted by us^36–39 and others (see, e.g., Refs. 40 and 41). The results of these works showed that among the GGA functionals, those with a small enhancement factor such as AM05,⁴² PBEsol,⁴³ or a few others^44–46 are the most accurate for the lattice constant and bulk modulus, while the standard PBE⁴⁷ is the best choice for the cohesive energy. At the MGGA level, the SCAN functional⁴⁸ is becoming increasingly popular and has been shown to be simultaneously as good as the best GGAs for the geometry (e.g., PBEsol) and the cohesive energy (PBE).^39,41 However, it should be mentioned that SCAN can be quite problematic for iterant magnetic systems^49,50 or alkali metals.⁵¹ 

Many semilocal functionals have been implemented directly in the WIEN2k code, but basically all existing semilocal functionals can be used because WIEN2k is interfaced to the Libxc^33,34 library of XC functionals. One current limitation is that the MGGA functionals are not yet implemented self-consistently (by default, the GGA PBE potential is used for generating the orbitals although the user can choose another potential).

Beginning in the 21st century, hybrid functionals,⁵² which belong to the fourth rung of Jacob’s ladder, started to be extensively used for calculations of solids.^53–55 The one that is currently the most popular is HSE06,^56–58 which is a screened version of the other well-known PBE0.^59,60 In (screened) hybrid functionals, the exchange energy is a linear combination of a semilocal (SL) functional and the Hartree–Fock (HF) expression,

where

In Eq. (14), v is either the bare Coulomb potential $v=1/r−r′$ for unscreened hybrids or a potential that is screened at short or long range for screened hybrids. For solids, it is computationally advantageous to use a potential that is short range, for instance, the Yukawa potential $v=e−λr−r′/r−r′$ (Ref. 61) or $v=erfc(μr−r′)/r−r′$ (Ref. 56), where erfc is the complementary error function. Although hybrid functionals are also used for the total energy (geometry optimization, and cohesive energy), they are particularly interesting for properties derived from the electronic band structure such as the bandgap, for which they significantly improve upon standard GGA functionals such as PBE (see Refs. 62–65 for recent extensive benchmarking).

In WIEN2k, unscreened and screened hybrid functionals are implemented⁶⁶ according to the scheme of Massidda et al.,⁶⁷ which is based on the pseudo-charge method for calculating the Coulomb potential.¹² The treatment of the Coulomb singularity is done by multiplying the Coulomb potential by a step function,⁶⁸ which is very efficient compared to other methods.⁶⁹ The screened hybrid functionals use the Yukawa potential, and in Ref. 66, it was shown that the results obtained with the PBE-based hybrid YS-PBE0 for the bandgap are almost identical to those obtained with HSE06 (which uses the erfc screened potential), provided that the screening parameter is chosen appropriately [$λ=3/2μ$⁠, see Ref. 70]. Because of the double integral and summations over orbitals in the HF exchange [Eq. (14)], the calculations are much more expensive (between 10 and 1000 times) than semilocal methods; however, there are a couple of ways to speed-up such calculations significantly. For instance, one can first use a rather crude k-mesh and later on improve the k convergence in a few additional iterations continuing the previous calculations. Furthermore, a reduced k-mesh for the internal loop in the HF potential is possible,⁷¹ and finally, often a one-shot procedure^72,73 is sufficient, where the hybrid orbitals and eigenvalues are calculated perturbatively on top of a calculation with the semilocal functional on which the hybrid functional is based.

Calculations using hybrid functionals in WIEN2k can be found in Refs. 74–77 for applications and in Refs. 39, 64, and 78–80 for various benchmark studies.

The high computational cost of hybrid methods discussed in Sec. II B 2 limits the size of the systems that can be treated. Alternatively, one can use an on-site method, namely, DFT+U,⁸¹ exact exchange for correlated electrons (EECE),⁸² or on-site hybrids,⁸³ which can be viewed as approximate but cheap versions of the hybrid or HF methods. In these methods, a hybrid/HF treatment is applied only to the electrons of a particular angular momentum belonging to a selected atom. However, using such an on-site scheme only makes sense when the considered electrons are well localized around the atom, which is, in general, the case for strongly correlated electrons. The on-site methods are mostly applied to open 3d-, 4f-, or 5f-shells in strongly correlated materials in order to improve the description of the electronic and magnetic properties. For such systems, the standard GGA methods provide results that are often even qualitatively inaccurate.⁸⁴ 

Different versions of DFT+U exist in the literature, and those available in WIEN2k are the following: (1) the original version,⁸¹ called Hubbard in mean field (HMF) in WIEN2k, (2) the fully localized limit version,^85,86 called self-interaction correction (SIC), and (3) the around mean-field (AMF) version.⁸⁶ The details of the implementation of DFT+U in the LAPW method can be found in the work of Shick et al.,⁸⁷ while a very good summary and discussion of the DFT+U flavors is given in Ref. 88. Note that since the on-site term is applied only inside the sphere surrounding the atom of interest,⁸⁷ the results may depend on the radius R_MT of this sphere (see Refs. 89 and 90 for illustrations), which is a drawback of the on-site methods. Among our works reporting DFT+U calculations, we mention Refs. 77 and 91–93.

For many technical aspects, the EECE and on-site hybrid methods are quite similar to DFT+U; however, there are two conceptual differences. The first one concerns the double-counting term. While in DFT+U, the double-counting term is derived using concepts from the Hubbard model (see Ref. 88 for a summary of the various expressions), in EECE and on-site hybrids, it is given by the semilocal expression of E_xc (e.g., PBE) evaluated with the density of the strongly correlated electrons.⁸² The second difference is the calculation of the Slater integrals in the Coulomb and Hartree–Fock terms. In DFT+U, they are parameterized with screened intra-atomic Coulomb (U) and exchange (J) interactions, which are usually chosen empirically. However, in EECE and on-site hybrids, the Slater integrals are calculated explicitly using the orbitals of the strongly correlated electrons.⁸² The results obtained with the on-site hybrid methods can be found in Refs. 83 and 94–98.

The results obtained with DFT+U and on-site hybrids should be qualitatively similar in many (but not necessary all) cases.^83,99–101 Actually, both methods contain empirical parameters: U and J (or only U_eff = U − J) in DFT+U and α_x [Eq. (13)] in on-site hybrids. In both cases, the results will depend crucially on the value of the parameters (U, J) or α_x. For applications, the EECE method⁸² is less interesting since it consists of 100% of unscreened Hartree–Fock exchange applied to correlated electrons, which usually is not accurate.

Technically, we mention that calculations with the on-site methods can only be done in a spin-polarized mode, i.e., with runsp_lapw. However, it is possible to apply on-site methods to non-magnetic systems by using the script runsp_c_lapw, which constrains the system to have no spin polarization.

It is well known that the GGA functionals that are commonly used for total-energy calculations, such as PBE, provide bandgaps that are much smaller than experiment.⁵⁷ Thus, one has to resort to other methods to get reliable results for the bandgap. Hybrid functionals and the GW method^102,103 (see Sec. III I 4) provide much more accurate values; however, they are also significantly more expensive than semilocal methods and cannot be applied easily to very large systems. Therefore, fast semilocal methods have been proposed specifically intended for bandgap calculations, and those which are available in WIEN2k are discussed below. Note that a more detailed discussion of the DFT methods for bandgaps is provided in Ref. 104.

The Tran–Blaha modified Becke–Johnson (TB-mBJ) potential¹⁰⁵ consists of a modified version of the BJ potential¹⁰⁶ for exchange (which reproduces the exact KS potential of atoms very well^106,107) and LDA¹⁰⁸ for correlation. The exchange part, which is a MGGA since it depends on the kinetic-energy density t, is

where $vxBR$ is the Becke–Roussel (BR) potential¹⁰⁹ and

with

being the average of $∇ρ/ρ$ in the unit cell. The parameters in Eq. (16) are α = −0.012, β = 1.023 bohrs^1/2, and p = 1/2 and were determined by minimizing the mean absolute error of the bandgap for a set of solids.¹⁰⁵ As shown in benchmark studies,^64,65,80 the TB-mBJ potential is currently the most accurate semilocal method for bandgap prediction. Other parameterizations of Eq. (16) were proposed in Refs. 110 and 111 and are also available in WIEN2k.

The GLLB-SC potential^112,113 is given by

where $exPBEsol$ is the PBEsol exchange-energy density per electron, $vcPBEsol=δEcPBEsol/δρ$ is the PBEsol correlation potential, and $KxLDA=82/3π2$⁠. Since the GLLB-SC potential depends on the orbital energies (ε_H is the one at the valence band maximum), a non-zero derivative discontinuity^114,115 can be calculated and added to the KS bandgap for comparison with the experimental value.^113,116 Similar to TB-mBJ, the GLLB-SC potential is significantly more accurate than traditional GGA functionals, as shown in Refs. 80, 104, 113, and 117. On the other hand, we note that these potentials are not obtained as the functional derivative of an energy functional.

Among other DFT methods, which have been shown to provide bandgaps more accurately than PBE, are the GGAs EV93PW91,^118,119 AK13,^120,121 and HLE16,¹²² as well as the LDA-type functional Sloc.¹²³ As mentioned in Sec. II B 1, the potential of MGGA energy functionals is not implemented in WIEN2k; however, it is still possible to calculate bandgaps non-self-consistently using the total energy (see Ref. 124 for details). Such MGGAs that are particularly interesting for bandgaps are HLE17¹²⁵ and TASK.¹²⁶ 

Figure 3 shows results for the bandgap of 76 solids, which we considered in our previous works.^64,80 Compared to the standard PBE functional, the results are much improved when TB-mBJ, GLLB-SC, or HSE06 is used, since the mean absolute error (MAE) drops from 1.99 eV with PBE to 0.47 eV, 0.64 eV, or 0.82 eV, respectively, for the other methods. Among all methods considered in Refs. 64 and 80, TB-mBJ leads not only to the smallest MAE but also to a slope (b = 0.97) of the linear fit that is closest to 1 (the AK13¹²⁰ functional also leads to a slope of 0.97).

Finally, we also mention that the Slater¹²⁷ and Krieger–Li–Iafrate¹²⁸ potentials have been implemented in the WIEN2k code.^129,130 However, these ab initio potentials, which are as expensive as the HF/hybrid methods, are not really intended for bandgap calculations, but may be interesting for more fundamental studies or as a better starting point for approximating the exact KS exchange potential.

The semilocal and hybrid functionals are, in general, quite inaccurate for describing weak interactions.^131,132 This is mainly due to the London dispersion forces that are not included properly in these approximations. Nevertheless, much better results can be obtained by adding to the semilocal/hybrid functional a correlation term (E_c,disp) accounting for the dispersion forces. There are essentially two types of dispersion corrections. The first one is of the atom-pairwise (at-pw) type,

where $CnAB$ are the dispersion coefficients for atom pair A–B separated by the distance R_AB and $fndamp$ is a damping function. Among the at-pw dispersion methods, Grimme’s DFT-D2¹³³ and DFT-D3¹³⁴ can be used with WIEN2k through the separate software DFT-D3¹³⁵ that supports periodic boundary conditions.¹³⁶ A feature of the DFT-D3 method is that the dispersion coefficients are not precomputed and fixed but depend on the coordination number of the system.

The second type of dispersion methods available in WIEN2k are the nonlocal van der Waals (NL-vdW) functionals¹³⁷ 

where the kernel Φ is a function of ρ, ∇ρ, and the interelectronic distance $r−r′$⁠. Compared to the single integral for semilocal functionals [Eq. (12)], evaluating such a double integral is clearly more involved computationally. In order to make the calculations affordable, Román-Pérez and Soler¹³⁸ proposed a very efficient method based on fast Fourier transformation (FFT) for evaluating Eq. (20). Their method is used for the implementation of NL-vdW functionals in WIEN2k. However, to apply the FFT-based method efficiently in an all-electron method, it is necessary to make ρ smoother by removing the high-density region close to the nucleus. This is done with the following formula:

where A = 1 bohr³ and ρ_c is the density cutoff that determines how smooth the new density ρ_s should be. More details can be found in Ref. 139, where it is shown that converged benchmark results can easily be obtained at a relatively modest cost.

In our previous works,^39,140 a plethora of DFT-D3 and NL-vdW functionals were assessed on solids. The test set consists of strongly bound solids and van der Waals solids such as rare gases or layered compounds. The results showed that among the at-pw methods, PBE-D3/D3(BJ)^134,141 seems to be a pretty good choice, while among the nonlocal methods, rev-vdW-DF2¹⁴² is the most balanced and actually more or less the best among all tested functionals.

The total energy E_tot of a periodic solid (with frozen nuclear positions) is given by Eq. (9). The individual terms are of opposite sign and, in an all-electron method, very large. In order to cancel the Coulomb singularity, we follow the algorithm of Weinert et al.,¹⁴³ where the kinetic and potential-energy terms are combined and a numerically stable method is obtained.

The force exerted on an atom t residing at position R^t, defined as

is calculated from the Hellman–Feynman theorem and includes Pulay corrections,¹⁴⁴ thereby taking into account that parts of the basis set used in WIEN2k depend upon the position of atoms (see Refs. 17 and 145–147 for the derivation specific to APW-based methods).

WIEN2k exploits the self-consistency of the KS equations, running through a sequence of calculations where the target is a density (and other parameters), which when passed through these SCF calculations yield the same density. This is equivalent to finding the solution to a set of simultaneous equations and is, in the most general case, referred to mathematically as a “fixed-point problem,” although it goes under the different name of “mixing” in the DFT literature. The general method used for just the electron density ρ and other density-like variables (for instance, the density matrix) is discussed in Ref. 148, while the extension to include atomic positions is described in Ref. 149. At any given iteration n, the relevant variables can be written as a vector $(ρn,vnorb)$⁠, including the density ρ at the start of an iteration as a function of the Cartesian coordinates r and an orbital potential v^orb (if used as in DFT+U or on-site hybrids, see Sec. II B 3), as well as other relevant variables. After running through the SCF sequence, a new density [symbolized by the SCF mapping $KS(ρn,vnorb)$] is produced. The fixed-point for just the density is when the two are equal, i.e., for all of the variables, the set of simultaneous equations

where $D(ρn,vnorb)$ is the density residual. Since the total energy with respect to the position of the atoms also has to be minimal, when these are allowed to vary, this can be expanded to include the forces $Fn(ρn,Rnt,vnorb)=−∇RntE(ρn,Rnt,vnorb)$⁠, i.e., solve the larger problem

This is equivalent to solving the SCF problem for an extended KS-equation, finding the variational minimum of both the atomic positions and densities.

The Hellman–Feynman forces due to the input density ρ_n are calculated within lapw0, while the Pulay corrections are calculated in lapw2 using the new density $KS(ρn,Rnt,vnorb)$⁠, so the forces above are, in general, not true derivatives of the energy, rather pseudo-forces that converge to them as the density converges. The general method is to expand to first order, i.e., write for the next value of the variables,

where H_n is an approximation to the inverse Jacobian, which is constructed as a Simplex gradient, i.e., a multi-dimensional numerical derivative using some number of the prior steps. The approximate inverse Jacobian has two components:

1. A predicted component where the density, positions, and other variables have changed in a way that maps onto the prior steps in the SCF iterations, so some information is already available.

2. An unpredicted component where the changes in the variables are new, so no prior information is available.

Equation (25) ignores higher-order terms in the expansion, which can break down far from the fixed point. In addition, because it is generated by a type of numerical differentiation, it can have limited accuracy if the step sizes are inappropriate. In addition, the unpredicted component can lead to instabilities. The approach taken is to control the algorithm greed for the unpredicted step and also use trust regions for both the predicted and unpredicted steps. The general idea of trust regions is that Eq. (25) is only reasonably accurate for changes of the variables, which are smaller than some value, which is called the trust radius, as illustrated in Fig. 4. Starting from some initial defaults, at each iteration, the algorithm will check to see if the step the mixer proposes to use is small enough; if it is too large, then the step is reduced. In the next iteration, if the step used led to an adequate improvement when bounded by the trust radius, then this radius is increased; if the step was not good, then the trust radius is decreased. In addition to changing the trust radii based upon improvement (or not), the most recent algorithm¹⁵⁰ also looks at the last step to see how large it should have been for both the predicted and unpredicted parts. This approach is significantly more stable and often leads to much smaller steps than earlier versions of the mixer used.

By design, the algorithm requires minimal user input beyond an estimate for the initial step to take for the unpredicted step; the most recent version of the algorithm¹⁵⁰ automatically controls all the internal parameters. No algorithm is perfect, and the convergence of the mixing depends significantly upon the nature of the physical problem being considered. The better the description of the underlying quantum mechanical problem, primarily the XC potential, the more rapidly it will converge. A very badly posed problem, in contrast, may converge only very slowly or not at all.

Two main algorithms are used within WIEN2k. The first is MSEC3, which is an updated version of a multisecant Broyden method¹⁴⁸ with trust region controls. This is a conservative algorithm, which uses the least greedy approach at every iteration. It is recommended for problems that converge badly. The more pushy MSR1¹⁴⁹ uses a more aggressive algorithm, which is significantly better for problems with soft modes that may converge only very slowly with MSEC3.

One unique feature of WIEN2k is that it can simultaneously converge both the density and atomic positions by solving the fixed-point problem of Eq. (24) using a multisecant approach.¹⁴⁹ This is often considerably faster than converging them independently as done in many other DFT codes and different from molecular dynamics approaches such as Car–Parrinello.¹⁵¹ The convergence rate depends upon the number and width of the eigenvalue clusters¹⁴⁹ of the combined electron and atomic position Jacobian. This approach does not follow the Born–Oppenheimer surface, which is the energy surface when the density is converged, rather some other surface, which is a balance between having converged densities and pseudo-forces as illustrated in Fig. 5. As such it can be somewhat confusing to the user, particularly as the pseudo-forces can vary in a strange fashion. This mode can be used with both the more conservative MSEC3 and the more aggressive MSR1 algorithm.

For the optimization of lattice parameters, WIEN2k offers a couple of workflows and utilities to generate structures with different lattice parameters, running the corresponding SCF calculations and analyzing the results. The optimized lattice parameters, however, are found only from the lowest total energy since there is no stress tensor in WIEN2k yet. This makes the optimization tedious for low symmetry cases and practically impossible for triclinic lattices.

WIEN2k consists of a large set of individual programs (mostly written in Fortran 90), which are linked together via tcsh-shell scripts representing a particular workflow. With this modular structure, WIEN2k is, on the one hand, very flexible and one can run a dedicated program for a particular task. On the other hand, there is not just one program and the specific task will be determined by directives in the input file, but a user has to know which program performs this specific task.

WIEN2k can be driven either from the command line or using a web-based graphical user interface (GUI), called w2web, which can be accessed by any web browser. Most likely, an experienced user will use the command line and explore all advanced features of WIEN2k, but for the beginner, the web-based GUI provides a very good starting point and it also teaches the user the corresponding command line.

The first task of every calculation is to define the structural data. As an example, the StructGen@w2web page is shown in Fig. 6 for the case of TiCoSb. The necessary basic input consists of the following:

• the lattice type (P, B, F, and H for primitive, body centered, face centered, and hexagonal, respectively) or, if already known, one of the 230 space groups; for the Heusler compound TiCoSb, we can select F lattice or space group 216 (⁠$F4¯3m$⁠);

• the lattice parameters a, b, and c (in bohr or Å) and the angles α, β, and γ (in degree);

• the atoms and their positions; if the space group is given, only one of the equivalent atoms has to be specified.

When the new structure is saved, the setrmt utility determines the nearest neighbor distances and automatically sets optimized atomic sphere sizes R_MT for this structure. The choice of R_MT has nothing to do with ionic radii but depends on the convergence properties of the atoms as discussed before (see the end of Sec. II A). It is important to note that if one wants to compare total energies for a series of calculations (e.g., for volume optimization), the R_MT should be kept constant.

An alternative on the command line is the makestruct utility, which works analogous to StructGen@w2web. More complex structures can be converted from cif or xyz files using the cif2struct or xyz2struct utilities. The generated structures can be conveniently visualized using XCrysDen¹⁵² or VESTA.¹⁵³ 

Starting from a basic (simple) structure, WIEN2k has powerful tools to generate supercells and manipulate them. supercell generates quickly h × k × l supercells (with/without B- or F-centering so that the supercell size can be increased by factors of two) and can add vacuum for surface slab generation. The structeditor, a collection of GNU Octave scripts, is even more powerful, since it can create arbitrary supercells (e.g., $3×3×l$⁠), rotate or merge structures, and delete or add atoms.

As mentioned above, WIEN2k consists of many individual programs and most of them have their own input file. Although this sounds very tedious at first, there are default inputs for all programs and several tools for changing the most important parameters on the fly. In w2web, the next step would be to check the symmetry of the newly generated structure and generate the input files for the SCF calculation (initialize@w2web). The user can provide a couple of parameters (only needed if one wants to change the defaults, see below) and run the following steps in batch mode or step by step:

• nn: Determines the distances between all atoms up to twice the nearest neighbor distance. In addition, it checks for overlapping spheres and will issue an error message if the spheres overlap. It also checks if identical elements have the same environment and eventually regroups them into equivalent sets.

• sgroup: Checks the structure and determines the space-group. It will group the atoms into sets of equivalent ones according to the Wyckoff positions of the corresponding space-group. In addition, it will check and determine the smallest possible (primitive) cell and create the corresponding structure file if necessary. For instance, if one enters the NaCl structure as a primitive cubic structure with four Na and four Cl atoms, it will automatically create a primitive FCC cell with only one Na and Cl atom.

• symmetry: Finds the symmetry operations of the space group as well as the point group symmetry of each atom and the corresponding LM expansion for the density/potential [Eqs. (1) and (2)].

• lstart: Solves numerically the radial Dirac equation for free atoms and creates atomic densities. Using the eigenvalues (or the localization within the atomic spheres) of all atomic states, it groups them into core and valence states. It selects automatically LOs for semi-core states and writes the starting energy parameters E_ℓ to case.in1 (during the SCF cycle, they are searched and adapted automatically to ensure best possible settings in all cases).

• kgen: Generates a shifted or non-shifted equidistant k-mesh with a user specified density in the irreducible part of the BZ.

• dstart: Superposes the atomic densities and creates the starting density for the SCF cycle.

On the command line, a corresponding script is called, which optionally allows us to specify various parameters (the most important ones are given below with their default values for reference):

init_lapw [-b -vxc PBE -ecut -6.0 -rkmax 7.0 -numk 1000]

The switches are described as follows: -b indicates batch mode (instead of step by step), -vxc selects the DFT functional, -ecut gives the core–valence separation energy (in Ry), -rkmax determines the plane wave cut-off parameter $RMTminKmax$⁠, and -numk determines the total number of k-points in the full BZ.

The most critical parameter is $RMTminKmax$⁠, which determines not only the quality but also the required computing time. The type of atom with the smallest R_MT determines this value because the R_MT for other atoms are set such that when the smallest atom is converged with the number of PWs, all others are also converged with the number of PWs. If the smallest sphere is a H atom (for instance, in short O–H bonds), $RMTminKmax=3$ is sufficient, most sp/d/f-elements converge with $RMTminKmax=7$/8/9. For lower (higher) precision, one can decrease (increase) these values by 10%–20%.

Of similar importance is the selection of a k-mesh. Generally speaking, small unit cells and metallic character require a large number of k-points (typical starting values for the SCF cycle would be a 10 × 10 × 10 mesh), while large cells (100 atoms) and insulators can be started with only one k-point. In any case, after the first SCF cycle, the k-mesh should be increased and the results (e.g., the forces on the atoms) should be checked. Certain properties (DOS, optics, and NMR) may need an even denser mesh, which in WIEN2k can be easily done and is fairly cheap, since it is only used for the property of interest.

The SCF cycle consists in WIEN2k of a complex workflow using several different programs. The main steps are as follows:

• lapw0: Calculates the Coulomb and XC potential from the density.

• lapw1: Calculates the valence and semi-core eigenvalues and eigenvectors at all requested k-points

• lapw2: Calculates the valence electron density

• lcore: Calculates the core eigenvalues and the core density

• mixer: Adds up the core and valence densities and mixes the total density with densities from previous iterations. In addition, it may update the atomic positions according to the calculated forces (see Sec. II C) and also the density matrices or orbital potentials when DFT+U or on-site hybrid methods are used.

Additional programs may be called depending on the requested options to include SO coupling or one of the specialized functionals discussed in Sec. II B (DFT-D3, NL-vdW, DFT+U, on-site hybrid/EECE, or hybrid-DFT/HF).

In w2web, the SCF cycle can be started by clicking on SCF@w2web. In this interface, one can then specify several parameters such as convergence criteria, parallelization, simultaneous optimization of internal atomic positions, or adding SO coupling.

The most important parameters for the corresponding command line script are as follows:

run_lapw [-ec 0.0001 -cc 0.0001 -fc 1.0 -p -so -min]

The SCF cycle will stop when the (optional) convergence criteria -ec (energy in Ry), -cc (charge in e⁻), and -fc (forces in mRy/bohr) are fulfilled three times in a row. SO coupling (only possible after a previous init_so_lapw step) is switched on using -so, -min relaxes the atomic positions simultaneously with the electron density (Sec. II C), and -p switches on parallelization (Sec. II E).

The basic summary of the SCF cycle is written into the case.scf file and all relevant quantities are labeled :LABEL: and can be searched/monitored using analyse@Utils@w2web or the Linux grep command. If the desired convergence has been reached, it is advisable to save all relevant input/output files using either save@Utils@w2web or the save_lapw utility. One can now either check results using more k-points (kgen@single_prog@Execution@w2web) or modify other inputs (input files@Files@w2web) such as $RMTminKmax$ or the XC functional and then continue with the SCF cycle. Later on, it is always possible to come back to a previously saved calculation using restore@Utils@w2web (restore_lapw).

Once these steps have been finished, one could, for instance, optimize the lattice parameters (optimize@w2web) or perform various other tasks (Tasks@w2web) such as Bandstructure@Tasks@w2web, DOS@Tasks@w2web, ElectronDensities@Tasks@w2web, XSPEC@Tasks@w2web, TELNES@Tasks@w2web, or OPTIC@Tasks@w2web (see Fig. 6). Each of these tasks consists of a guided workflow and let the user prepare the necessary inputs, executes various small programs, and visualizes the results.

WIEN2k runs on any Linux platform and also on Mac. It is written mainly in Fortran 90 (a few programs are written in C), and the workflows are managed by tcsh scripts. Most of the time-critical parts use libraries such as BLAS and LAPACK, and efficient libraries are therefore mandatory. There is direct installation support for the standard Linux tools GFortran+OpenBLAS (at least gcc 6.x) and Intel IFORT+MKL. The latter still gives the best performance. w2web is a specialized web server written in perl and listens on a user-defined high port. Its access is, of course, password protected and can be limited to specific IP addresses.

For the optional installation of the MPI-parallel version (useful only on clusters with InfiniBand network or larger shared memory workstations with at least 16 cores), one needs obviously MPI (e.g., Open MPI or Intel MPI) and also ScaLAPACK (included in the MKL), FFTW,¹⁵⁴ and, optional but highly recommended, ELPA.^155,156

The following Linux tools are necessary (not all of them are always installed by default): tcsh, Perl 5, Ghostscript, gnuplot, GNU Octave, and Python 2.7.x+NumPy.

Optional, but highly recommended, programs for certain tasks include the following:

• XCrysDen¹⁵² and VESTA¹⁵³ for structure and electron density visualization but also generation of band structure k-meshes or plotting Fermi-surfaces.

• Libxc^33,34 for XC functionals not directly implemented in WIEN2k.

• DFT-D3¹³⁵ for DFT+D3 calculations of van der Waals systems.

• Wannier90^157,158 for constructing Wannier functions using the wien2wannier utility.

• phonopy,^159,160 Phonon,^161,162 or PHON^163,164 for phonon calculations.

• BoltzTraP2¹⁶⁵ for transport calculations (see Sec. III I 1).

• fold2bloch¹⁶⁶ to fold supercell band structures back to the primitive BZ.

• SKEAF¹⁶⁷ to extract de Haas-van Alphen frequencies from WIEN2k.

• Critic2^168,169 is an alternative program to the WIEN2k program aim to analyze 3D scalar fields such as the electron density by using the “atoms in molecules” (AIM) theory of Bader.^170,171

WIEN2k is a highly accurate all-electron code based on the APW method and thus certainly not as fast as some other (pseudopotential or minimal basis set) codes. However, it takes advantage of inversion symmetry and when present it will automatically use the “real” instead of the “complex” version of the code, thus saving half of the memory and running almost four times as fast. In addition, it is highly optimized and efficiently parallelized at three different levels, which can be optimally chosen depending on the size of the problem and the available hardware.

Except for OpenMP parallelization (see below), parallelization is activated by a -p switch in our scripts and needs a .machines file as listed and described below:

• OpenMP parallelization:

The main (time consuming) programs are all parallelized using OpenMP and can use the corresponding threaded BLAS, LAPACK, and FFTW libraries. It is activated by either setting the OMP_NUM_THREADS variable globally or using omp_prog:N directives in .machines. While many parts of the code scale very well with the number of parallel threads on a multi-core shared memory machine, unfortunately, the scaling of the matrix diagonalization is, at present, limited to 2–4 cores due to performance bottlenecks in the corresponding OpenBLAS or MKL libraries.

• k-point parallelization:

This together with OpenMP is a very simple and highly efficient parallelization, which works even on a loosely coupled cluster of simple PCs with a slow network for small to medium sized cases (up to 100 atoms/cell) where the eigenvalue problem needs to be solved for several k-points. It requires a common (NFS) filesystem on all machines and password-less ssh (private/public keys). N lines speed:hostname in .machines will split the list of k-points into N junks, and N jobs will be started in parallel on the corresponding hosts, followed by a summation step of the partial densities. Since WIEN2k can use temporary (local) storage for the eigenvectors, we are not limited in the number of k-points and our personal record is a NMR chemical shift calculation for fcc Al with 10⁶k-points.

• MPI-parallelization:

With a sufficiently powerful hardware (at least 16 cores or a cluster with InfiniBand network) and for medium to large sized problems (more than 50 atoms/cell), it is possible, and actually necessary, to parallelize further using MPI. Besides a tremendous speedup that can be achieved by parallelization over atoms and, in particular, over basis functions, this version will distribute the necessary memory on all requested computers, thus allowing calculations for unit cells with more than 1000 atoms.^172–174 Such cells require basis sets of about 10⁵ APWs, and the resulting Hamiltonian, overlap, and eigenvector matrices may need about 500 GB of memory, which are distributed over the nodes in the standard ScaLAPACK block-cyclic distribution. For such large systems, the solution of the general eigenvalue problem can become the time-limiting (N³) step (depending on the number, type, and R_MT of the atoms, the setup of the complicated matrix elements can take a comparable fraction of the total time), but the ELPA library provides a highly efficient and scalable (1000 cores) diagonalization. In cases with fewer atoms in large cells (isolated molecules in a large box or surface slabs with sufficient vacuum), we can use an iterative diagonalization¹⁷⁵ using the previous eigenvectors as start. Depending on the requested number of eigenvalues, this method may be up to 100 times faster than full diagonalization and still scales very well with the number of cores. The MPI version of the code is used when lines with speed:hostname:N (or lines with more than one hostname) in .machines are specified. Of course, coupling of k- and MPI-parallelization (and/or OpenMP) is possible.

Once a self-consistent solution for a chosen atomic structure is done, one can focus on the electronic structure. The energy eigenvalues as a function of the k-vector obtained at the end of the KS calculation define the band structure. The k-path along high symmetry lines in the irreducible BZ (see the Bilbao Crystallographic Server¹⁷⁶) can be either obtained from WIEN2k default templates or generated graphically using XCrysDen.¹⁵² In WIEN2k, one can plot the energy bands by using the program spaghetti, indicating that its interpretation is difficult. However, there are some tools to help. A first tool is a symmetry analysis, which determines the irreducible representation (of the corresponding point group) for each KS eigenvalue. With this knowledge, one can connect the KS eigenstates to bands by using compatibility relations and satisfying the non-crossing rule. The chemical bonding information of state nk is contained in the corresponding wave function ψ_nk, which is complex and three dimensional. However, when computing the square of its modulus, one obtains an electron density, which is a real function and easy to visualize. Integrating this (normalized) electron density, one obtains a charge q, which can be decomposed into contributions from the interstitial region I and the atomic spheres S_t (labeled by the atom number t and the quantum number ℓ according to the atomic-like basis set),

This allows us to compress the detailed information contained in the wave function of a single eigenstate state ψ_nk to a few numbers that can be stored and analyzed. In addition, WIEN2k decomposes the q_tℓ according to the symmetry of the corresponding point group. For example, the five d-orbitals of a TM atom surrounded by ligands in octahedral symmetry are split into the t_2g and e_g manifold (crystal field splitting), while for lower symmetry, a splitting into five different d-orbitals is obtained. A review paper⁸ (Sec. 6.2) illustrates these tools for TiC, a refractory metal (crystallizing in the sodium chloride structure) that is almost as hard as diamond but has metallic, covalent, and ionic bonding contributions. These data are very useful to analyze the electronic structure, and we illustrate this for the band structure of the Heusler compound TiCoSb (see Ref. 177). If one wants to know which atomic states (e.g., Co-d, Ti-d, or Sb-p states) contribute most to a certain band, one can show the character of the bands, sometimes called “fat bands.” For each eigenvalue ε_nk, the size of the circle represents the weight of the chosen character (e.g., a particular q_tℓm). Figure 7 shows which band states originate mainly from Co-d, Ti-d, or Sb-p states, giving the band structure a chemical interpretation.

From the KS eigenvalues calculated on a sufficiently fine k-grid in the irreducible BZ, one can obtain the density of states (DOS), usually by means of the (modified) tetrahedron method.¹⁷⁸ By using the partial charges [Eq. (26)], one can decompose the total DOS into partial DOS (PDOS), which are useful for understanding chemical bonding and interpreting various spectroscopic data. Figure 8(a) shows how much each region (atomic spheres of Co, Ti, Sb, and the interstitial) contributes to the total DOS. We show the dominating valence contributions from the Co-3d and Ti-3d electrons in Fig. 8(b) and those from Sb-5s/5p in Fig. 8(c). Unfortunately, the interstitial PDOS cannot be decomposed into atomic and ℓ-like contributions uniquely. However, by analyzing the atomic orbitals in the free atom, we see (Table I) that only a fraction of the related electron density resides inside the corresponding atomic sphere, e.g., 81% for the Co-3d but only 15% for the Co-4s orbital. Therefore, a significant part of the density lies outside the atomic sphere, leading to a non-negligible PDOS from the interstitial region [Fig. 8(a)]. A simple scheme to eliminate this interstitial part is to renormalize the partial DOS with a factor $qtℓren$ (determined by a least squares fit) such that the sum of the renormalized PDOS contributions yields the total DOS,

This sum runs only over the “chemical” ℓ, which are the main contributions. In Table I, we see that the $qtℓren$ are close to the free atom situation for more localized orbitals but differs significantly, for example, for Co-4s, which is more localized in the solid than in the free atom. For the importance of this effect, see Sec. III B.

The fundamental variable in DFT is the electron density ρ, which can be compared to experimental data. The total ρ, which is obtained by summing over all occupied states, can be decomposed into its contributions coming from the core, semi-core, and valence states. A variety of tools (such as XCrysDen¹⁵² or VESTA¹⁵³) allows one to visualize the density ρ along a line, in a plane (2D), or in the unit cell (3D). One can easily compute the density corresponding to a selected energy window of electronic states in order to visualize their bonding character. By taking the difference between the crystalline density and a superposition of atomic densities (placed at the atomic position of the crystal), one obtains a difference density Δρ = ρ_crystal − ρ_atoms, which shows chemical bonding effects much more clearly than the total or valence density. Figure 9 provides an illustration of Δρ for TiCoSb within the (110) plane, where we can observe the strong asphericities in the electron density around the Ti and Co atoms originating from different occupations of the five 3d orbitals, as well as the charge transfer (discussed below using Bader charges).

When we want to compare the computed electron densities or the related X-ray structure factors (computed using the lapw3 module) with experimental data, we must take into account the motion of the nuclei. In DFT calculations, we assume that the nuclei are at rest, whereas in an experiment, this motion must be considered, for example, by means of the Debye–Waller factors, which can also be calculated by phonon calculations.

It is a strength of theory to allow various decompositions (of the DOS or electron density, for instance), which are often useful for interpreting properties, but these may depend on the basis set used in a calculation, for example, when deriving atomic charges. In a LCAO scheme, one takes the weights of all atomic orbitals centered at a given atom to determine how much charge corresponds to that atom (Mulliken’s population analysis). In an APW scheme, the charge inside the related atomic sphere would give an atomic charge, but this value clearly depends on the chosen atomic radius and lacks the interstitial contribution. However, the renormalized partial charges $Qtren$ obtained from an integral over the renormalized PDOS gives a meaningful measure of charge transfer, as shown in Table II.

A basis-set independent alternative is the AIM procedure proposed by Bader,^170,171 which is based on a topological analysis of the density. It uniquely defines volumes (called “atomic basins”) that contain exactly one nucleus by enforcing a zero-flux boundary: $∇ρ⋅n^=0$⁠. Inside such an atomic basin, this scheme uniquely defines the Bader charge for a given density independently of the basis-set method that was used to calculate the electron density.^170,171 An example of the application of the AIM method is given in Table II, which shows the charge inside the atomic basins of TiCoSb determined with the aim module of WIEN2k. $Qt,crystalBader$ is the nuclear charge Z_t minus the number of electrons (a positive value indicates a depletion of electrons) using the SCF density, and $Qt,superBader$ is the same quantity but using a density from a superposition of the free neutral atoms. In this crystal structure, according to a Bader analysis, even the superposition of neutral densities leads to a significant charge transfer from Ti to Co and Sb, which is enhanced for Co and Ti but reduced for Sb during the SCF cycle. These Bader charges can be compared to the $Qtren$⁠. As we can see, in both methods, there is a transfer of electrons from the Ti to the Co atom, but the specific amount and, in particular, the charge state of Sb differs significantly depending on the way it is calculated. An inspection of the difference density (Fig. 9) shows a negative Δρ around Sb and thus indicates a positive charge of the Sb atom in contrast to the Bader charges, which seems to pick up a lot of charges in the interstitial region leading to a negative Sb charge. In essence, one should be careful with quantitative charge state assignments.

An alternative to aim is the Critic2 package,^168,169 which determines Bader charges using a pre-calculated 3D mesh of densities. It is very fast; however, the integration of total charges on such a crude mesh is inaccurate and one should restrict its usage for magnetic moments (integrating the spin densities) or valence charges densities (be careful with the 3D mesh).

Experimental valence band photoelectron spectra (PES) are often just compared to the total DOS. Such a comparison, however, can at best reproduce certain peak positions, but usually not the experimental intensities. This is even more true with modern synchrotron-based hard X-ray PES (HAXPES), where the spectra differ considerably depending on the excitation energy. This is because the cross sections of different atomic orbitals change dramatically as a function of excitation energy and this effect should be taken into account. The pes module¹⁷⁹ of WIEN2k uses the partial DOS (PDOS_tℓ) and multiplies it with the corresponding energy-dependent atomic orbital cross sections σ_tℓ, taken from various tables,^180,181

In addition, pes can use the renormalized PDOS (see Sec. III A), so that the contributions from the less localized orbitals (whose wave functions are mainly in the interstitial region) are also properly taken into account. This module allows us to specify the X-ray energy and can handle unpolarized and linearly polarized light as well as linear dichroism in angular distribution (LDAD). It was successfully applied for various examples¹⁷⁹ (SiO₂, PbO₂, CeVO₄, In₂O₃, and ZnO).

Here, we compare in Fig. 10 the experimental HAXPES spectrum¹⁷⁷ of TiCoSb at 6 keV with the theoretical calculation. The theoretical spectrum reproduces the experimental intensities very well, but the bandwidth is too small so that the Sb-s peak has about 1 eV less binding energy. This is a well-known DFT problem and concerns all states at lower energy. The decomposition of the total spectrum allows us to analyze the contributions to the different peaks in the spectrum. The low energy feature is almost exclusively from Sb-s states, the double peak at −5 eV is from Sb-p states, and the double peak at −2 eV is from Co-d states. However, Co-s also contributes significantly to the feature at −6 eV, and Sb-p contributions are even larger than Co-3d for the lowest binding energy peak at −1.5 eV. All Ti contributions are very small and thus not shown in Fig. 10.

X-ray photoelectron spectroscopy (XPS) determines the binding energy (BE) of core states. These BEs are specific to certain atoms, but the possible small changes of BEs (core level shifts) provide important additional information about the chemical environment and, in particular, the oxidation state of that element. Since WIEN2k is an all-electron method, it has the self-consistent core eigenvalues available and one could calculate their BE as energy difference with respect to the Fermi level. However, according to Janak’s theorem,¹⁸² DFT eigenvalues represent the partial derivative of the total energy with respect to the orbital occupancy and are therefore not necessarily good approximations of experimental excitation energies. Such BEs are typically 10%–20% too small. In fact, even BE differences (core level shifts) from ground-state calculations might not be reliable because screening effects of the final state are not included. Much better approximations to experimental BEs can be obtained according to Slater’s transition state theory,¹⁸³ where half of a core electron is removed.^173,184,185 The corresponding SCF eigenvalue, which represents the slope of the total energy vs occupation at half occupation, is a much better approximation to the actual energy difference for n and n − 1 occupation, and typical BE errors are reduced to a few percent. For solids, such calculations should employ large supercells where only one of the atoms gets excited. This method also allows for some possible screening due to the valence electrons but suffers from the fact that, in solids, a neutral unit cell is required. The missing half electron can be compensated by adding a negative background charge, by increasing the number of valence electrons by one half, or by playing slightly with the nuclear charge¹⁸⁶ according to a virtual crystal approximation. However, it is not always clear which of these methods should be preferred. Successful applications include, for instance, the N-1s shifts of h-BN covered Pt, Rh, and Ru(111) surfaces (with a unit cell of the “nanomesh” containing more than 1000 atoms)¹⁷³ or to Pb-5d and Ta-4f shifts in the misfit layer compound (PbS)_1.14TaS₂.¹⁸⁵ 

The work function can also be obtained from surface slab calculations as the difference of the Fermi level and the Coulomb potential in the middle of the vacuum region. An example can be found for free and h-BN covered Ni and Rh(111) surfaces in Refs. 187 and 188. One has to carefully check the convergence of the work function with respect to the size of the vacuum region.

Experimental techniques such as X-ray emission (XES), near-edge X-ray absorption (XAS, NEXAFS, and XANES), and electron energy loss spectroscopy (EELS) represent an electronic transition between a core state and a corresponding valence/conduction band state, which leads to the measurement of emitted/absorbed X-rays or the energy loss of transmitted electrons. The intensity of such a spectrum is given by Fermi’s golden rule according to dipole transitions between an initial (Ψ_I) and a final (Ψ_F) state,

The dipole selection rule is valid when the X-ray energy is not too large and limits transitions between a core state on atom X and angular momentum ℓ into/from a conduction/valence band state with Δℓ ± 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, the orientation-dependent spectra can be obtained by substituting the ℓ-like partial DOS by an appropriate ℓm-like DOS, e.g., for a K-spectrum of a tetragonal/hexagonal system by replacing the total p-DOS by p_z and p_x + p_y-DOS.¹⁸⁹ 

Such a scheme leads to very good results for XES spectra, where the final state has a filled core-hole and the valence-hole is usually well screened. For XANES and EELS spectra, however, the final state¹⁹⁰ determines the spectrum. The final state has a core hole and an excited electron in the conduction band, and they will interact with each other leading to strong excitonic effects. In order to describe this effect in a DFT-based band-structure code, one has to create a supercell (as large as possible, but depending on the specific system and the hardware resources) of about 32–256 atoms and remove a core electron from one of the atoms. This electron should be added either to the valence electrons (if there are proper states in the conduction band, e.g., in B–K edges of BN) or to the constant background charge (if the lowest conduction band states are of completely wrong character, e.g., for the O–K edge of a TM-oxide) to keep the system neutral. In the following SCF cycle, the valence states on the atom with the core hole will get a lower energy and localize, but the surrounding electrons are allowed to partially contribute to the screening.

As an example, the Cs-L₃ spectrum of CsK₂Sb is shown in Fig. 11(a) and compared with state-of-the-art calculations¹⁹¹ using the Bethe–Salpeter equation (BSE) (for further details on BSE, see Sec. III I 3). Obviously, the spectrum calculated with the ground state electronic structure is very different from core hole supercell or BSE calculations, where the spectral weight is redistributed into the first (excitonic) peak. In particular, the Cs-5d states come down in energy and partially screen the core hole [Fig. 11(b)]. The Cs-6s contributions (not shown) are 1–2 orders of magnitude smaller because both their dipole matrix elements and the 6s-PDOS are smaller than the corresponding Cs-5d quantities. On the other hand, the BSE and core hole supercell calculations agree quite well. Note that the BSE data have a larger broadening. In this example, the size of the supercell is easy to converge, but it should be noted that for details of the spectrum, a rather good k-mesh (8 × 8 × 8 for the 128 atom supercell) is necessary.

While the core hole approach works generally quite well, it also has clear limitations or needs extensions:

• For metals or small bandgap semiconductors, a full core hole is sometimes too much because the static screening in the supercell might not be enough. Better results can be obtained using a “partial” hole,^192,193 although adjusting the size of the hole until the resulting spectra match the experiment is not fully ab initio anymore.

• The used DFT approximation may not be accurate. This can concern the ground state, for instance, the O–K edge of NiO could be greatly improved using the TB-mBJ approximation¹⁹⁴ instead of PBE or PBE+U. Even more problematic are excited-state effects due to the additional d electron in strongly correlated materials (for instance, 3d in TM oxides), where very poor L_2,3 edges are obtained in single particle approaches and sophisticated methods such as dynamical mean field theory (DMFT, see Sec. III I 5) or configuration interactions¹⁹⁵ are needed.

• For early (“d⁰”) TM compounds, the L_2,3 edges are influenced strongly by interactions and interference effects between the 2p_1/2 and 2p_3/2 states, which are split only by a few eV. This can be accounted for using fully relativistic BSE calculations (see Sec. III I 3), where both the 2p_1/2 and 2p_3/2 states are taken into account simultaneously.¹⁹⁶ 

• The B-K edge of hexagonal BN has been investigated many times in the literature.^189,197 While the first strong excitonic peak originating from antibonding B-p_z (π^*) states is well described by supercell calculations and in full agreement with BSE calculations,¹⁹⁷ the experimental double peak at around 7 eV above the π^* peak originating from σ^* (B-p_x,y) states shows up in the calculations only as a single peak. This can only be fixed by taking electron–phonon interactions into account. The approach is based on statistical averages over all vibrational eigenmodes of the system.¹⁹⁷ Thus, one calculates first the vibrational modes of h-BN and then, using a supercell of e.g., 128 atoms, all atoms are displaced according to the vibrational eigenmodes with amplitudes determined by the Bose–Einstein occupations. Even at T = 0 K, the zero-point motion is enough to split the degenerate p_x and p_y states, since, on average, each B will have one N neighbor at a smaller (larger) distance than the two others. This leads to the desired splitting of the single σ^* peak into a double peak.¹⁹⁷ 

EELS has fairly similar basic principles as XAS but differs slightly because of the finite momentum of the electrons.¹⁹⁸ The telnes3 module of WIEN2k calculates the double differential scattering cross section on a grid of energy loss values and impulse transfer vectors. This double differential cross section is integrated over a certain momentum transfer q to yield a differential cross section, which can be compared to the experiment. This formalism allows the calculation of relativistic EELS including transitions of arbitrary order (i.e., non-dipole transitions), and it can take into account the relative orientation between the sample and beam.¹⁹⁹ Practical aspects on how to perform EELS calculations have been given by Hébert,²⁰⁰ and some examples can be found in Refs. 201–204.

The optics module of WIEN2k uses the independent-particle approximation (IPA) and calculates the direct transitions (conserving k) between occupied nk and unoccupied n′k states, where for both states, KS eigenvalues are used.²⁰⁵ The joint density of states is modified by transition probabilities given by the square of the momentum matrix elements M = ⟨n′k|A·p|nk⟩ between these states, which determine the intensity of optical spectra using dipole selection rules and clearly distinguish between optically allowed and forbidden transitions. From the resulting imaginary part ε₂ of the dielectric function, its real part ε₁ can be obtained by the Kramers–Kronig transformation and then additional optical functions such as conductivity, reflectivity, absorption, or the loss function can also be calculated. In a metallic solid, an additional Drude term accounts for the free-electron intraband contribution. For insulators and semiconductors, where the DFT gap is often too small when compared to the experiment, one can use a “scissor operator.” This sounds complicated but is nothing else than a rigid shift of the unoccupied DFT bands to adjust the (too small) DFT bandgap, either using the experimental gap or, more ab initio, using the gap calculated with TB-mBJ. Note that TB-mBJ usually gives very good bandgaps, but the bandwidth of both the valence and conduction bands are too small, and thus, the optical properties with TB-mBJ might not be very accurate, but still more accurate than standard GGA (see Sec. III B 2 in Ref. 104 for a brief summary of literature results with TB-mBJ). Alternatively, hybrid-DFT functionals⁶⁶ can be used, which give quite good bandgaps for semiconductors, but one should be aware that the optical properties usually require a quite dense k-mesh, which makes hybrid calculations fairly expensive.

As an example we present in Fig. 12 the imaginary part of the dielectric function ε₂ for CsK₂Sb using various approximations. First, we note that very dense k-meshes are necessary for converged results, which makes the application of more expensive many-body perturbation methods such as GW even more difficult. As expected, PBE calculations yield the smallest bandgap of 1.06 eV, while hybrid YS-PBE0 gives 1.68 eV and TB-mBJ 2.08 eV. This can be compared to G₀W₀ results¹⁹¹ of 1.62 eV or early experimental estimates of 1.0–1.2 eV, which, however, have been criticized. It should be noted that the GW result¹⁹¹ for ε₂ probably suffers from an under converged k-mesh and even with their large smearing (note the large tail below 1.62 eV in the GW results shown in Fig. 12), a distinct peak structure emerges, which is not present in the k-converged results.

Experimental results for optical conductivity, reflectivity, or absorption as well as the low energy valence electron energy loss spectrum (VEELS) can often be successfully interpreted in the IPA.^203,204 However, sometimes (in particular, for wide gap insulators) the frequency-dependent dielectric function ε in the IPA may have little in common with the experimental situation. This is because excitations are two-particle processes and the missing electron–hole interaction (i.e., the excitonic effect already mentioned in Sec. III C) can significantly affect the calculated optical response of a material when they are strong. In order to overcome this problem, one needs to include the electron–hole correlation explicitly by solving the BSE (see Sec. III I 3).

When magnetism occurs in a solid, it may come from localized electrons (e.g., from f electrons of rare-earth atoms) or itinerant (delocalized) electrons (e.g., in Fe, Co, or Ni). In any case, magnetism comes mainly from exchange splitting causing a partial occupation of states, which differ between the spin-up (N_↑) and spin-down (N_↓) electrons. The corresponding magnetic moment M is defined as the difference between these occupation numbers (M = N_↑ − N_↓). For such cases, one must perform spin-polarized calculations (runsp_lapw) and needs the spin density in addition to the total electron density. The default is collinear magnetic order as found in ferromagnets, for example, in Fe, Co, Ni, or antiferromagnets, for example, in Cr. In addition to collinear magnets mentioned here, one can also handle non-collinear magnetism (for example, systems with canted magnetic moments or spin spirals), as described in Sec. III I 2.

For spin-polarized calculations of a specific complex antiferromagnetic structure, in most cases, it is essential to specify proper (antiferromagnetic) atomic spin-moments as an input for the SCF cycle and the WIEN2k tool instgen allows us to set this easily. If one is interested in the orientation of the magnetic moments with respect to the crystal structure (easy or hard axis) or the magneto-crystalline anisotropy,²⁰⁶ the SO interaction must be included. For heavy elements or when orbital moments become important, one needs this full relativistic treatment by including the SO interaction.

Nowadays, one can study complicated systems, for example, BaFe₂O₅, an oxygen-deficient perovskite-like structure, which shows a Verwey transition. At low temperature, this system has a charge-ordered state (with Fe²⁺ and Fe³⁺ at different Fe-sites), but above the Verwey transition temperature (at about 309 K), a valence mixed state with the formal oxidation state Fe^2.5+ appears. DFT calculations made it possible to interpret this complicated situation (for details see Ref. 207 and Sec. 7.4.1 of Ref. 208). Another complex system is PrO₂, which has a Jahn–Teller-distorted CaF₂ structure.²⁰⁹ It contains Pr-4f electrons, which form a localized band (lower Hubbard band) for one 4f electron, but the others hybridize with the valence electrons forming a situation between Pr³⁺ and Pr⁴⁺. This compound is an antiferromagnetic insulator that requires a relativistic treatment. PBE+U calculations⁹³ provide results that are consistent with all experimental data for the bandgap, magnetic moment, and structural distortion.

In the 1980s, a numerical problem occurred in connection with several studies of the Fe₆₅Ni₃₅ INVAR alloy, which has a vanishing thermal expansion around room temperature. This is one of the systems for which the magnetization shows a hysteresis when a magnetic field is applied. The hysteresis causes numerical difficulties, since for a certain magnetic field, there are three solutions (magnetic moments) with very similar total energies causing difficulties in the convergence of a conventional SCF procedure. In order to solve this problem, the fixed spin moment (FSM) method was proposed.^210,211 It is a computational trick interchanging dependent and independent variables. Physically speaking, one applies a magnetic field and obtains a moment, but computationally one chooses the moment (as input) and calculates the field afterwards. In a conventional spin-polarized calculation, the Fermi energy must be the same for the spin-up and spin-down electrons. The magnetic moment M is an output. In the FSM scheme, one does several constrained calculations, where the moment M is an input, but allowing different Fermi energies for the two spin states. One can interpret the difference in the Fermi energies as a magnetic field. Although one needs to perform several calculations (instead of a single conventional one), they converge rather rapidly. The FSM method allows expanding the usual total-energy vs volume curve to an energy surface E_tot(V, M) as a function of volume V and moment M, which also provides new insights.

All aspects of nucleus–electron interactions, which go beyond the electric point-charge model for a nucleus, define the hyperfine interactions. Nuclei with a nuclear quantum number I ⩾ 1 have an electric quadrupole moment Q. The nuclear quadrupole interaction (NQI) stems from the interaction of such a moment and the EFG, the second derivative of the Coulomb potential at the corresponding nuclear site. One can measure the EFG with Mössbauer spectroscopy, NMR, nuclear quadrupole resonance (NQR), or perturbed angular correlation (PAC). The NQI determines the product of Q and the EFG, a traceless tensor. The latter has a principal component and an asymmetry parameter η. The EFG is a ground-state property that can be determined experimentally (measuring NQI), provided the nuclear quadrupole moment is known. In early studies, the EFG was interpreted as a simple point charge model with additional corrections (Sternheimer factor²¹²). However, later it was shown that the EFG can be calculated with DFT, as was illustrated²¹³ for LiN₃. The EFG is sensitive to the asymmetric charge distribution around a given nucleus and thus is a local probe, which often helps in clarifying the local atomic arrangement. The reader can find a short description of several EFG calculations for selected examples in Chap. 6.4 of Ref. 8, and here, we describe some important aspects and examples below.

Computationally, it is important to treat both valence and semi-core states very accurately because, due to the 1/r³ factor in the EFG expression, even small asphericities near the nucleus lead to important contributions. This was demonstrated for the first time for TiO₂ in the rutile structure, where the radial functions of the fully occupied Ti-3p semi-core orbitals are slightly different for p_x, p_y, and p_z states and thus contribute significantly to the EFG.²⁸ 

The mapping of the two Cu-EFGs in the high-temperature superconductor YBa₂Cu₃O₇ to the “plane” and “chain” Cu sites provides insight into which Cu atom is responsible for superconductivity, and the analysis of the EFG on all other sites helped us to interpret chemical bonding.²¹⁴ It should be mentioned that the EFG at the (superconducting) Cu-plane site comes out quite inaccurate and GGA+U calculations are necessary, leading to a redistribution of 0.07 e⁻ from a $3dx2−y2$ orbital to a $3dz2$ orbital and an EFG in agreement with the experiment.⁹² 

Next, we briefly discuss the study on 16 fluoroaluminates from Ref. 215, for which experimental NMR data were compared to DFT results for the ²⁷Al EFG. In all of these compounds, the aluminum atoms occur in Al$F63−$ octahedra, which have a wide diversity of connectivity and distortions. One of these structures is shown in the inset of Fig. 13 for Ba₃Al₂F₁₂. These fluoroaluminates illustrate how sensitive the EFG is to the exact position of neighboring atoms. A perfect octahedral symmetry would have a vanishing EFG, but small distortions cause an EFG. The calculations were first done using the less accurate powder diffraction data for the atomic positions [Fig. 13(a)], and the correlation between experiment and theory is not very good. Then, a DFT structure optimization was done leading to an almost perfect correlation between experimental and theoretical EFGs [Fig. 13(b)]. This structure optimization has an even more pronounced effect on the asymmetry parameter, as shown in Ref. 215. We should mention that there can also be a sensitivity to DFT functionals for EFG calculations as described in Refs. 66 and 80.

Last but not least, we demonstrate how one can determine the nuclear quadrupole moment Q from a combination of theoretical EFG calculations and experimental measurements of the quadrupole splitting, which is proportional to the product of EFG and Q. From the slope of a linear regression for the EFG of several Fe compounds, we could deduce the nuclear quadrupole moment Q of ⁵⁷Fe, the most important Mössbauer isotope. It was found to be about twice as large (Q = 0.16 b) as the previous literature value (Q = 0.082 b), suggesting to revise this nuclear property using electronic structure calculations.²¹⁶ 

The magnetic hyperfine field (HFF) at a nucleus originates from a Zeeman interaction between the magnetic moment I of this nucleus and the magnetic field at this site produced by the spin-polarized electrons in a ferromagnet. The HFF has contributions from the Fermi-contact term (the spin density at the nucleus), an orbital, and a spin dipolar contribution. Here, we skip the details but mention that an all-electron treatment is crucial, especially for the Fermi-contact term, since one needs accurate values of the spin-density close to the nucleus, for which the basis set used in WIEN2k is extremely useful. We recalculate the densities of all electrons, including the core, in each cycle of the SCF scheme, in contrast to the frozen-core approximation. The resulting core polarization can often be the main contribution to the HFF. The high quality of such calculations was demonstrated for the double perovskite BaFe₂O₅, for which the DFT calculations of EFG and HFF provided new insights (for details, see Ref. 207).

The NMR shielding $σ⃡$ tensor is defined as a constant between an induced magnetic field B_ind at the nucleus at site R and the external uniform field B_ext,

Its value is usually in the range of ppm (part per million). Since the magnetic field cannot be controlled with such a precision, the tensor is measured only with respect to some reference,

and often only its isotropic part $σ(R)=tr[σ⃡(R)]$ is known.

The external magnetic field is a relatively weak perturbation compared to the typical energy scale of the electronic structure; therefore, its effect on the spin and orbit of an electron can be separated in the theoretical calculations. Here, we only outline the main features that are specific to the APW+lo method and the WIEN2k code and are vital for achieving high accuracy of the computed NMR tensor. A more in-depth discussion can be found in the original publications. The formalism for computing the orbital part of the shielding (chemical shift) has been described in Refs. 217 and 218, and the spin part (Knight shift) of the response has been described in Ref. 219. The formalism has been applied for computing shielding in various insulating^78,220,221 and metallic systems.^222–224 

As will be explained below, this approach can reach the basis set limit for NMR calculations, and benchmark calculations for small molecular systems have proven that standard quantum chemistry methods can only reach this precision with very large uncontracted quintuple-zeta basis sets and only for light atoms.²²⁵ 

The orbital part of the shielding, i.e., the orbital component of the induced field B_ind, is obtained directly from the Biot–Savart law (in atomic units, with c as speed of light),

where j(r) is the induced orbital current, evaluated as expectation value of the current operator

WIEN2k separates the calculation of valence and core states. The core state contribution to the induced current is computed using the spherically symmetric core density only,

The method for computing the valence contribution to j(r) is based on a linear response approach^226–228 originally developed by Mauri, Pfrommer, and Louie.²²⁶ The expression for the induced current involves only the first-order terms with respect to the external field B_ext,

where $ψo(0)$ is an unperturbed KS occupied orbital, J⁰(r) is the paramagnetic part of the current operator [the first term in Eq. (33)], and J¹(r) is the diamagnetic component of the current operator [the second term in Eq. (33)]. $ψo(1)$ is the first-order perturbation of $ψo(0)$⁠, given by the standard formula involving a Green’s function,

and involves a sum over empty states (first term) as well as a sum over core states (second term) because the core states have been calculated before [Eq. (34)] and Eq. (36) is only correct if all states of a system are included. Note that core wave functions appear in Eq. (36) as if they were unoccupied states. H⁽¹⁾ is the perturbation due to the external magnetic field in symmetric gauge,

As discussed in Sec. II A, the basis functions of the APW+lo method are highly tuned to describe the occupied Bloch states everywhere in the unit cell (in particular, also close to the nucleus), but they are not a complete basis set. Therefore, depending on the perturbation of the Hamiltonian, they may not be well suited to expand the perturbations of wave functions. In fact, magnetic fields and NMR are such cases where the perturbed wave functions, in particular, near the nucleus, are very different, which means that the sum in Eq. (36) cannot be converged with the available set of orbitals. In order to remedy this issue, we had to enhance the original WIEN2k basis set. First, the standard set of local orbitals is extended significantly, both in the number of LOs per ℓ (typically to 5–8 LOs) and also in ℓ (typically to ℓ + 1), where ℓ refers to the maximal “chemical” ℓ of this atom.²¹⁷ These extra local orbitals are referred to as NMR-LOs, and the energy at which the radial functions of those NMR-LOs are computed is chosen such that each of the NMR-LO radial functions has zero value at the sphere boundary and the number of nodes inside the sphere of subsequent LOs increases by one corresponding to the next principal quantum number.²¹⁷ However, these NMR-LOs alone cannot completely improve the variational flexibility close to the nuclei.

The perturbation of the Hamiltonian due to the external magnetic field is proportional to a product of position and momentum operators. As a result, the perturbation of the radial wave function u_tℓ contains components proportional to u_tℓ±1 and their radial derivative $r∂∂rutℓ±1$⁠. A direct introduction of basis functions based on $r∂∂ru$ is not convenient within the APW formalism because such functions are not eigenstates of the radial Schrödinger equation. Therefore, we have proposed to add the desired term directly to the Green’s function present in the formula for the first-order perturbation of the valence state wave function.²¹⁸ It is referred to as the “$∂∂ru$ correction” (DUC).

The convergence test with respect to the number of NMR-LO and DUC corrections is presented in Fig. 14. The induced current and shielding calculated within our linear response formalism is compared to the exact value for an isolated Ar atom. The current and shielding for a spherically symmetric atom can be computed exactly using only its density and the same formula [Eq. (34)] as for the core states (diamagnetic current).

Both DUC and several additional NMR-LOs are needed to reproduce the shape of the exact diamagnetic current in a region within 0.5 bohr from the nucleus. It appears in Fig. 14(a) that any error in the representation of the current in this region results in substantial errors of the computed shielding values [Fig. 14(b)]. This method can therefore reach the basis set limit.

The all-electron nature and the modular concept of WIEN2k makes it very easy to perform NMR calculations with wave functions including SO interactions for heavy nuclei or using different DFT approximations, including DFT+U (Sec. II B 3), but, in particular, also hybrid functionals (Sec. II B 2). Thus, we can compare the theoretical shielding with the experimental chemical shifts for several different compounds, and from the correlation and slope of the linear regression curve, the quality of a particular approximation to the XC effects can be evaluated. Ideally, the slope of this linear regression line should be −1, but typically for ionic compounds,⁷⁸ the slope with PBE is too large (−1.2), while with the hybrid functional YS-PBE0,⁶⁶ it is too small (−0.8). This is quite in contrast to organic molecules, where hybrid functionals perform much better than the GGAs. Surprisingly, the BJ potential¹⁰⁶ performs quite well for ionic oxides or halides and yields slopes close to −1.

Most importantly, theory should not only reproduce measured experimental values but also provide insights. WIEN2k allows us to analyze and identify the contributions to the NMR shielding, and for instance, for the F-shielding in the alkali-fluoride series (Fig. 15), the following observations can be made:²²¹ (i) Basically, all contributions to the F-shielding come from a region inside the F atomic sphere (the current in the rest of the unit cell contributes negligibly). (ii) The large (diamagnetic) shielding comes from the constant F-1s core contribution. (iii) The contributions from the F-2s bands are still diamagnetic, but much smaller and again constant within the series. (iv) The diamagnetic metal-p semi-core contributions (Na-2p to Cs-5p) increase within the series. This can be explained by the fact that for heavier elements the metal p-states increase in energy and come closer to the F-2p band. This leads to an increased bonding (in the metal-d band) and anti-bonding (in the F-p band) metal-p–F-p interaction giving slightly different (non-canceling) diamagnetic and paramagnetic contributions. (v) The trend of σ within the series comes mainly from the F-2p valence band. The most important ingredient, which determines the size of the (mostly paramagnetic) F-2p contribution, is the position of the unoccupied metal-d band. The perturbation due to the magnetic field couples the occupied F-p states with Δℓ ± 1 to unoccupied d states, and due to the energy denominator in Eq. (36), the Cs-d states give the largest contribution because they are the closest in energy to the valence bands. We can even artificially apply a (large) U value to the empty Cs-5d states shifting them further up. In this way, we do not alter the occupied states but still can increase the F shielding in CsF to reproduce the LiF or NaF shifts.

By a similar analysis, we could explain why the ³³S magnetic shielding decreases with the metal nuclear charge Z in the ionic alkali/alkali-earth sulfides but increases in TM sulfides.²²⁰ 

In order to compute the induced spin density and spin part of the NMR shielding tensor, we use a direct approach,²¹⁹ instead of applying the linear response formalism proposed, for instance, in Ref. 229. This is possible because the interaction of the spin with the external magnetic field does not break the periodicity. Therefore, we perform self-consistent spin polarized calculations with a finite external magnetic field B_ext acting on the electron spin only. The interaction with B_ext is cast into a spin-dependent potential leading to a spin splitting of eigenstates and a finite spin magnetization. The induced magnetic field at a given nucleus is computed using an expression for the magnetic hyperfine field,²³⁰ 

The first term (B_c) is the Fermi contact term, where m_av is the average of the spin density in a region near the nucleus with a diameter equal to the Thomson radius. The second term (B_sd) captures the spin-dipolar contribution to the hyperfine field, where Φ₁ is the large component of the wave function, S is the reciprocal relativistic mass enhancement, and μ is the magnetic moment operator of the electron. B_sd comes almost entirely from within the atomic sphere, which simplifies its calculation. The spin contribution (σ_s), i.e., the Knight shift to the shielding, is therefore given by two terms,

In order to obtain a sizable response and evaluate the NMR shielding with a numerical precision at the level of 1 ppm, we apply in our calculations an external magnetic field of 100 T, which induces a spin-splitting of approximately 1 mRy. These small changes require an extremely fine k-point sampling (for fcc Al, 10⁶k-points are needed), and this must always be carefully checked.

Further details and results of our approach can be found in Refs. 219, 222, and 223, but, here, we summarize the main findings: (i) The previously accepted point of view, namely, that Knight shifts are proportional to the partial s-DOS at the Fermi energy and the orbital contribution σ_orb is identical to that in the (ionic) reference compound, is only true for simple sp metals. (ii) In TM d-elements or metallic compounds, the orbital part σ_orb can be as important as the spin part σ_s. (iii) The s-DOS at E_F is always important, but an induced TM-d magnetic moment (proportional to the partial d-DOS at E_F) polarizes the core states in the opposite direction so that the valence and core polarizations can partly cancel. (iv) The dipolar contribution σ_sd is usually small, but, in anisotropic materials, a large dominance of one particular orbital at E_F can eventually lead to a very large contribution. We have found this in BaGa₂, where the p_z-DOS has a large and sharp peak at E_F leading to an aspherical magnetization density and a large dipolar contribution.²²² 

A single particle state of a periodic system is conventionally represented as a Bloch state ψ_nk(r), which is labeled by a band index n and a vector k inside the first BZ. It satisfies Bloch’s theorem

where u_nk(r) = u_nk(r + R) is a lattice periodic function and R is a Bravais lattice vector. Alternatively, one can define Wannier functions (WFs) w_nR(r) in unit cell R for a set of J bands as^231,232

where V is the unit cell volume and U^k are unitary transformation matrices that mix Bloch states at a given k. Because of the arbitrary phase of ψ_nk(r), the resulting WFs are usually not localized. Maximally localized WFs can be obtained by choosing the U^k such that the spread Ω is minimized,

This involves overlap integrals $Mmnk,b=⟨umk|unk+b⟩$ between the periodic part of the wave functions on a uniform grid in the BZ.

The maximally localized WFs are calculated by Wannier90,²³¹ and the transformation matrices U^k are provided by the wien2wannier²³² module of WIEN2k. The resulting WF can be used for various tasks. They can be visualized and are useful for the interpretation of chemical bonding to generate tight-binding models or interpolations to very fine k-meshes for properties that require fine k-meshes. Such properties can be transport, anomalous Hall conductivity, linear and non-linear optics, Berry curvatures and topology, or electron–phonon interactions. In particular, for the improved description of electronic correlations via the DMFT approximation,^233,234 WFs provide a realistic starting point.

In connection with wien2wannier, there is also the BerryPi²³⁵ module in WIEN2k, which calculates the polarization of solids using the Berry phase approach.²³⁶ BerryPi can calculate the change of polarization ΔP in response to an external perturbation to study ferroelectricity, the Born effective charges, pyroelectric coefficients, or the piezoelectric tensor. In addition, one can define Wilson loops and calculate Chern numbers to study topological properties.²³⁷ 

WIEN2k is interfaced with the BoltzTraP2 program¹⁶⁵ for calculating transport coefficients within the relaxation time approximation. The calculation is based on evaluating the transport distribution function,

using a fine mesh in k-space. To obtain the group velocities, v_nk, and also quasi-particle energies on a fine mesh or effective masses, BoltzTrap2 relies on interpolating the eigenvalues, ε_nk, and possibly also the relaxation times, τ_nk, using Fourier sums.

The interpolation is performed so that the calculated eigenvalue energies are reproduced exactly. Within KS theory, the multiplicative potential [Eq. (10)] means that it is often computationally very efficient to calculate a fine mesh of eigenvalues, which can then be interpolated further to evaluate Eq. (43). This argument no longer holds when hybrid functionals (Sec. II B 2), or the GW method (Sec. III I 4), are used to obtain the band structure. Therefore, BoltzTraP2 can also include the k-space derivatives in the interpolation scheme, as can be obtained from the momentum matrix elements introduced in Sec. III D. This allows a more efficient interpolation and the use of a coarser k-mesh in the actual DFT calculation.

Once the transport distribution has been obtained, the temperature and chemical potential dependent transport coefficients

can be obtained by a simple numerical integration. Figure 16 shows the highest valence bands of TiCoSb calculated using the interpolation scheme of BoltzTraP2. Compared to Fig. 7, a very good agreement with the band structure obtained using DFT is found. Since the Fourier interpolation is done bandwise (marked by color in Fig. 16), the band crossing along the Γ − X direction is not reproduced. However, the error is hardly visible by eye and, in accordance with the intention of the original algorithm,^238,239 no Fourier ripples are seen, which means that errors in the derivatives are isolated to the points where the crossing occurs. The calculated thermoelectric power factor using a constant relaxation time is also shown. The power factor peaks at a high value close to the band edge, where a steep transport distribution can be expected. The high power factor can be attributed to the complex constant energy surface shown in the inset, which is typical for half-Heusler compounds with favorable p-type thermoelectric performance.²⁴⁰ 

BoltzTraP2 is written mainly in PYTHON3 and can be used as PYTHON library. The interpolation is handled by a single PYTHON call, fitde3D. Once the Fourier coefficients have been obtained, the interpolation of the bands onto the direction needed for plotting the band structure (getBands), or the fine mesh needed for obtaining the transport distribution (getBTPbands), can be performed. The use of BoltzTraP2 as a library gives a reproducible and flexible work flow. The analysis associated with Fig. 16 can thus be performed with a single PYTHON script, which is included in the most recent distribution.

WIEN2k can only compute the electronic structure of magnetic systems with a collinear spin arrangement. For performance benefits, WIEN2k assumes that the spin density matrix along some direction is diagonal for each eigenstate. When SO interactions are taken into account, this condition is sometimes not satisfied. In such cases, the off-diagonal terms of the spin density matrix are simply ignored during the SCF procedure and only the z component of the spin density is properly converged. Generally, this is not a big issue for cases with collinear spin arrangement and is, in fact, common practice. However, such an approximation cannot be applied for systems with a non-collinear spin arrangement for which the full spin density matrix has to be considered. For that purpose, we have written a non-collinear spin version of WIEN2k, referred to as WIENNCM.

Our implementation is based on a mixed spinor basis set approach.^241,242 In the interstitial region, the basis functions are pure spinors given in a global (g) spin coordinate frame,

where $χ↑g=10$ and $χ↓g=01$⁠. Inside the atomic spheres, the basis functions are a combination of both up and down spinors, which are set in a local spin coordinate frame with a quantization axis pointing along the direction of the average magnetization of the given atomic sphere. This direction does not have to be the same for each sphere, and the basis functions are [for a LAPW basis set, where we drop the (r, E)-dependency in the radial u functions]

where $χσt$ is a spinor given in a local coordinate frame. This choice of the spin coordinate frame allows us to use spin-polarized radial functions with the quantization axis along the direction of the average magnetization. The matching of the $φK+k,σLAPW$ basis to the plane waves at r = R_MT is done for up and down plane waves in a global spin coordinate frame,