Linearized augmented-plane-wave method

At the core of density functional theory the Hohenberg-Kohn theorems state that every observable of an interacting many-electron system is a functional of its ground-state charge density and that this density minimizes the total energy of the system.[5] The theorems do not answer the question how to obtain such a ground-state density. A recipe for this is given by Walter Kohn and Lu Jeu Sham who introduce an auxiliary system of noninteracting particles constructed such that it shares the same ground-state density with the interacting particle system.[6] The Schrödinger-like equations describing this system are the Kohn-Sham equations. With these equations one can calculate the eigenstates of the system and with these the density. One contribution to the Kohn-Sham equations is the effective potential which itself depends on the density. As the ground-state density is not known before a Kohn-Sham DFT calculation and it is an input as well as an output of such a calculation, the Kohn-Sham equations are solved in an iterative procedure together with a recalculation of the density and the potential in every iteration. It starts with an initial guess for the density and after every iteration a new density is constructed as a mixture from the output density and previous densities. The calculation finishes as soon as a fixpoint of a self-consistent density is found, i.e., input and output density are identical. This is the ground-state density.

A method implementing Kohn-Sham DFT has to realize these different steps of the sketched iterative algorithm. The LAPW method is based on a partitioning of the material's unit cell into non-overlapping but nearly touching so-called muffin-tin (MT) spheres, centered at the atomic nuclei, and an interstitial region (IR) in between the spheres. The physical description and the representation of the Kohn-Sham orbitals, the charge density, and the potential is adapted to this partitioning. In the following this method design and the extraction of quantities from it are sketched in more detail. Variations and extensions are indicated.

The central aspect of practical DFT implementations is the question how to solve the Kohn-Sham equations

${\displaystyle \left[{\hat {T}}_{\text{s}}+V_{\text{eff}}(\mathbf {r} )\right]\left|\Psi _{j}^{\mathbf {k} }(\mathbf {r} )\right\rangle =\epsilon _{j}^{\mathbf {k} }\left|\Psi _{j}^{\mathbf {k} }(\mathbf {r} )\right\rangle }$

with the single-electron kinetic energy operator ${\displaystyle {\hat {T}}_{\text{s}}}$, the effective potential ${\displaystyle V_{\text{eff}}(\mathbf {r} )}$, Kohn-Sham states ${\displaystyle \Psi _{j}^{\mathbf {k} }(\mathbf {r} )}$, energy eigenvalues ${\displaystyle \epsilon _{j}^{\mathbf {k} }}$, and position and Bloch vectors ${\displaystyle \mathbf {r} }$ and ${\displaystyle \mathbf {k} }$. While in abstract evaluations of Kohn-Sham DFT the model for the exchange-correlation contribution to the effective potential is the only fundamental approximation, in practice solving the Kohn-Sham equations is accompanied by the introduction of many additional approximations. These include the incompleteness of the basis set used to represent the Kohn-Sham orbitals, the choice of whether to use the pseudopotential approximation or to consider all electrons in the DFT scheme, the treatment of relativistic effects, and possible shape approximations to the potential. Beyond the partitioning of the unit cell, for the LAPW method the central design aspect is the use of the LAPW basis set ${\displaystyle \left\lbrace \phi _{\mathbf {k} ,\mathbf {G} }(\mathbf {r} )\right\rbrace }$ to represent the valence electron orbitals as

${\displaystyle \left|\Psi _{j}^{\mathbf {k} }(\mathbf {r} )\right\rangle =\sum \limits _{\mathbf {G} }c_{j}^{\mathbf {k} ,\mathbf {G} }\left|\phi _{\mathbf {k} ,\mathbf {G} }(\mathbf {r} )\right\rangle ,}$

where ${\displaystyle c_{j}^{\mathbf {k} ,\mathbf {G} }}$ are the expansion coefficients. The LAPW basis is designed to enable a precise representation of the orbitals and an accurate modelling of the physics in each region of the unit cell.

Considering a unit cell of volume ${\displaystyle \Omega }$ covering atoms ${\displaystyle \alpha }$ at positions ${\displaystyle \mathbf {\tau } _{\alpha }}$, an LAPW basis function is characterized by a reciprocal lattice vector ${\displaystyle \mathbf {G} }$ and the considered Bloch vector ${\displaystyle \mathbf {k} }$. It is given as

${\displaystyle \phi _{\mathbf {k} ,\mathbf {G} }(\mathbf {r} )=\left\lbrace {\begin{array}{l l}{\frac {1}{\sqrt {\Omega }}}e^{i(\mathbf {k} +\mathbf {G} )\mathbf {r} }&{\text{for }}\mathbf {r} {\text{ in IR}}\\\sum \limits _{l=0}^{l_{{\text{max}},\alpha }}\sum \limits _{m=-l}^{l}\left[a_{l,m}^{\mathbf {k} ,\mathbf {G} ,\alpha }u_{l,\alpha }(r_{\alpha },E_{l,\alpha })+b_{l,m}^{\mathbf {k} ,\mathbf {G} ,\alpha }{\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })\right]Y_{l,m}(\mathbf {\hat {r}} _{\alpha })&{\text{for }}\mathbf {r} {\text{ in MT}}_{\alpha }\end{array}}\right.,}$

where ${\displaystyle \mathbf {r} _{\alpha }=\mathbf {r} -\mathbf {\tau } _{\alpha }}$ is the position vector relative to the position of atom nucleus ${\displaystyle \alpha }$. An LAPW basis function is thus a plane wave in the IR and a linear combination of the radial functions ${\displaystyle u_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$ and ${\displaystyle {\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$ multiplied by spherical harmonics ${\displaystyle Y_{l,m}}$ in each MT sphere. The radial function ${\displaystyle u_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$ is hereby the solution of the Kohn-Sham Hamiltonian for the spherically averaged potential with regular behavior at the nucleus for the given energy parameter ${\displaystyle E_{l,\alpha }}$. Together with its energy derivative ${\displaystyle {\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$ these augmentations of the plane wave in each MT sphere enable a representation of the Kohn-Sham orbitals at arbitrary eigenenergies linearized around the energy parameters. The coefficients ${\displaystyle a_{l,m}^{\mathbf {k} ,\mathbf {G} ,\alpha }}$ and ${\displaystyle b_{l,m}^{\mathbf {k} ,\mathbf {G} ,\alpha }}$ are automatically determined by enforcing the basis function to be continuously differentiable for the respective ${\displaystyle (l,m)}$ channel. The set of LAPW basis functions is defined by specifying a cutoff parameter ${\displaystyle K_{\text{max}}=|\mathbf {k} +\mathbf {G} |_{\text{max}}}$. In each MT sphere, the expansion into spherical harmonics is limited to a maximum number of angular momenta ${\displaystyle l_{{\text{max}},\alpha }\approx K_{\text{max}}R_{{\text{MT}}_{\alpha }}}$, where ${\displaystyle R_{{\text{MT}}_{\alpha }}}$ is the muffin-tin radius of atom ${\displaystyle \alpha }$. The choice of this cutoff is connected to the decay of expansion coefficients for growing ${\displaystyle l}$ in the Rayleigh expansion of plane waves into spherical harmonics.

While the LAPW basis functions are used to represent the valence states, core electron states, which are completely confined within a MT sphere, are calculated for the spherically averaged potential on radial grids, for each atom separately applying atomic boundary conditions. Semicore states, which are still localized but slightly extended beyond the MT sphere boundary, may either be treated as core electron states or as valence electron states. For the latter choice the linearized representation is not sufficient because the related eigenenergy is typically far away from the energy parameters. To resolve this problem the LAPW basis can be extended by additional basis functions in the respective MT sphere, so called local orbitals (LOs).[7] These are tailored to provide a precise representation of the semicore states.

The plane-wave form of the basis functions in the interstitial region makes setting up the Hamiltonian matrix

${\displaystyle H_{\mathbf {G'} ,\mathbf {G} }^{\mathbf {k} }=\left\langle \phi _{\mathbf {k} ,\mathbf {G'} }{\Big |}{\hat {H}}{\Big |}\phi _{\mathbf {k} ,\mathbf {G} }\right\rangle =\left\langle \phi _{\mathbf {k} ,\mathbf {G'} }{\Big |}{\hat {T}}_{\text{s}}+V_{\text{eff}}(\mathbf {r} ){\Big |}\phi _{\mathbf {k} ,\mathbf {G} }\right\rangle }$

for that region simple. In the MT spheres this setup is also simple and computationally inexpensive for the kinetic energy and the spherically averaged potential, e.g., in the muffin-tin approximation. The simplicity hereby stems from the connection of the radial functions to the spherical Hamiltonian in the spheres ${\displaystyle {\hat {H}}_{\text{sphr}}^{\alpha }}$, i.e., ${\displaystyle {\hat {H}}_{\text{sphr}}^{\alpha }\left|u_{l,\alpha }(r_{\alpha },E_{l,\alpha })\right\rangle =E_{l,\alpha }\left|u_{l,\alpha }(r_{\alpha },E_{l,\alpha })\right\rangle }$ and ${\displaystyle {\hat {H}}_{\text{sphr}}^{\alpha }\left|{\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })\right\rangle =E_{l,\alpha }\left|{\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })\right\rangle +\left|u_{l,\alpha }(r_{\alpha },E_{l,\alpha })\right\rangle }$. In comparison to the MT approximation, for the full-potential description (FLAPW) contributions from the non-spherical part of the potential are added to the Hamiltonian matrix in the MT spheres and in the IR contributions related to deviations from the constant potential.

After the Hamiltonian matrix ${\displaystyle H_{\mathbf {G'} ,\mathbf {G} }^{\mathbf {k} }}$ together with the overlap matrix ${\displaystyle S_{\mathbf {G'} ,\mathbf {G} }^{\mathbf {k} }=\left\langle \phi _{\mathbf {k} ,\mathbf {G'} }{\Big |}\phi _{\mathbf {k} ,\mathbf {G} }\right\rangle }$ is set up, the Kohn-Sham orbitals are obtained as eigenfunctions from the algebraic generalized dense Hermitian eigenvalue problem

${\displaystyle \sum \limits _{\mathbf {G} }H_{\mathbf {G'} ,\mathbf {G} }^{\mathbf {k} }c_{j}^{\mathbf {k} ,\mathbf {G} }=\epsilon _{j}^{\mathbf {k} }\sum \limits _{\mathbf {G} }S_{\mathbf {G'} ,\mathbf {G} }^{\mathbf {k} }c_{j}^{\mathbf {k} ,\mathbf {G} }~~,}$

where ${\displaystyle \epsilon _{j}^{\mathbf {k} }}$ is the energy eigenvalue of the j-th Kohn-Sham state at Bloch vector ${\displaystyle {\mathbf {k} }}$ and the state is given as indicated above by the expansion coefficients ${\displaystyle c_{j}^{\mathbf {k} ,\mathbf {G} }}$.

The considered degree of relativistic physics differs for core and valence electrons. The strong localization of core electrons due to the singularity of the effective potential at the atomic nucleus is connected to large kinetic energy contributions and thus a fully relativistic treatment is desirable and common. For the determination of the radial functions ${\displaystyle u_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$ and ${\displaystyle {\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$ the common approach is to make an approximation to the fully relativistic description. This may be the scalar-relativistic approximation[8][9] (SRA) or similar approaches.[10][11] The dominant effect neglected by these approximations is the spin-orbit coupling. As indicated above the construction of the Hamiltonian matrix within such an approximation is trivial. Spin-orbit coupling can additionally be included, though this leads to a more complex Hamiltonian matrix setup or a second variation scheme,[12][13] connected to increased computational demands. In the interstitial region it is reasonable and common to describe the valence electrons without considering relativistic effects.

After calculating the Kohn-Sham eigenfunctions, the next step is to construct the electron charge density by occupying the lowest energy eigenstates up to the Fermi level with electrons. The Fermi level itself is determined in this process by keeping charge neutrality in the unit cell. The resulting charge density ${\displaystyle \rho (\mathbf {r} )}$ then has a region-specific form

${\displaystyle \rho (\mathbf {r} )=\left\{{\begin{array}{l l}\sum \limits _{\mathbf {G} }\rho _{\mathbf {G} }e^{i\mathbf {G} \mathbf {r} }&{\text{for }}\mathbf {r} {\text{ in IR}}\\\sum \limits _{\alpha }\sum \limits _{l=0}^{l_{{\text{max}},\alpha }}\sum \limits _{m=-l}^{l}\rho _{l,m}^{\alpha }(r_{\alpha })Y_{l,m}(\mathbf {\hat {r}} _{\alpha })&{\text{for }}\mathbf {r} {\text{ in MT}}_{\alpha }\end{array}}\right.,}$

i.e., it is given as a plane-wave expansion in the interstitial region and as an expansion into radial functions times spherical harmonics in each MT sphere. The radial functions hereby are numerically given on a mesh.

The representation of the effective potential follows the same scheme. In its construction a common approach is to employ Weinert's method for solving the Poisson equation.[14] It efficiently and accurately provides a solution of the Poisson equation without shape approximation for an arbitrary periodic charge density based on the concept of multipole potentials and the boundary value problem for a sphere.

Because they are based on the same theoretical framework, different DFT implementations offer access to very similar sets of material properties. However, the variations in the implementations result in differences in the ease of extracting certain quantities and also in differences in their interpretation. In the following, these circumstances are sketched for some examples.

The most basic quantity provided by DFT is the ground-state total energy of an investigated system. To avoid the calculation of derivatives of the eigenfunctions in its evaluation, the common implementation[15] replaces the expectation value of the kinetic energy operator by the sum of the band energies of occupied Kohn-Sham states minus the energy due to the effective potential. The force exerted on an atom, which is given by the change of the total energy due to an infinitesimal displacement, has two major contributions. The first contribution is due to the displacement of the potential. It is known as Hellmann-Feynman force. The other, computationally more elaborate contribution, is due to the related change in the atom-position-dependent basis functions. It is often called Pulay force[16] and requires a method-specific implementation.[17][18]

Independent of the actual size of an atom, evaluating atom-dependent quantities in LAPW is often interpreted as calculating the quantity in the respective MT sphere. This applies to quantities like charges at atoms, magnetic moments, or projections of the density of states or the band structure onto a certain orbital character at a given atom. Deviating interpretations of such quantities from experiments or other DFT implementations may lead to differences when comparing results. On a side note also some atom-specific LAPW inputs relate directly to the respective MT region. For example in the DFT+U approach the Hubbard U only affects the MT sphere.[19]

A strength of the LAPW approach is the inclusion of all electrons in the DFT calculation, which is crucial for the evaluation of certain quantities. One of which are hyperfine interaction parameters like electric field gradients whose calculation involves the evaluation of the curvature of the all-electron Coulomb potential near the nuclei. The prediction of such quantities with LAPW is very accurate.[20]

Kohn-Sham DFT does not give direct access to all quantities one may be interested in. For example most energy eigenvalues of the Kohn-Sham states are not directly related to the real interacting many-electron system. For the prediction of optical properties one therefore often uses DFT codes in combination with software implementing the GW approximation (GWA) to many-body perturbation theory and optionally the Bethe-Salpeter equation (BSE) to describe excitons. Such software has to be adapted to the representation used in the DFT implementation. Both the GWA and the BSE have been formulated in the LAPW context and several implementations of such tools are in use.[21][22][23] In other postprocessing situations it may be useful to project Kohn-Sham states onto Wannier functions. For the LAPW method such projections have also been implemented[24][25] and are in common use.

• APW: The augmented-plane-wave method[26] is the predecessor of LAPW. It uses the radial solution to the spherically averaged potential for the augmentation in the MT spheres. The energy derivative of this radial function is not involved. This missing linearization implies that the augmentation has to be adapted to each Kohn-Sham state individually, i.e., it depends on the Bloch vector and the band index, which subsequently leads to a non-linear, energy-dependent eigenvalue problem. In comparison to LAPW this is a more complex problem to solve. A relativistic generalization of this approach, RAPW, has also been formulated.[27]
• Local orbitals extensions: The LAPW basis can be extended by local orbitals[7] (LOs). These are additional basis functions having nonvanishing values only in a single MT sphere. They are composed of the radial functions ${\displaystyle u_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$, ${\displaystyle {\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$, and a third radial function tailored to describe use-case-specific physics. LOs have originally been proposed for the representation of semicore states. Other uses involve the representation of unoccupied states[28][29] or the elimination of the linearization error for the valence states.[30]
• APW+lo: In the APW+lo method[31][32] the augmentation in the MT spheres only consists of the function ${\displaystyle u_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$. It is matched to the plane wave in the interstitial region only in value. As an alternative implementation of the linearization the function ${\displaystyle {\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$ is included in the basis set as an additional local orbital. While the matching conditions result in an unphysical kink of the basis functions at the MT sphere boundaries, a careful consideration of the kink in the construction of the Hamiltonian matrix suppresses it in the Kohn-Sham eigenfunctions. In comparison to the classical LAPW method the APW+lo approach leads to a less stiff basis set. The outcome is a faster convergence of the DFT calculations with respect to the basis set size.
• Soler-Williams formulation of LAPW: In the Soler-Williams formulation of LAPW[33] the plane waves cover the whole unit cell. In the MT spheres the augmentation is implemented by replacing up to the angular momentum cutoff the plane waves by the functions ${\displaystyle u_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$ and ${\displaystyle {\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$. This yields basis functions continuously differentiable also in the ${\displaystyle (l,m)}$ channels above the angular momentum cutoff. As a consequence the Soler-Williams approach has reduced angular momentum cutoff requirements in comparison to the classical LAPW formulation.
• ELAPW: In the extended LAPW method[34][35] pairs of local orbitals introducing the functions ${\displaystyle u_{l,\alpha }(r_{\alpha },E_{l,\alpha }^{\text{lo}})}$ and ${\displaystyle {\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha }^{\text{lo}})}$ are added to the LAPW basis. The energy parameters ${\displaystyle E_{l,\alpha }^{\text{lo}}}$ are chosen to systematically extend the energy region in which Kohn-Sham states are accurately described by the linearization in LAPW.
• QAPW: In the quadratic APW method[36][37] the augmentation in the MT spheres additionally includes the second energy derivative ${\displaystyle {\ddot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })}$. The matching at the MT sphere boundaries is performed by enforcing continuity of the basis functions in value, slope, and curvature. This is similar to the super-linearized APW (SLAPW) method in which radial functions ${\displaystyle u_{l,\alpha }}$ and/or their derivatives ${\displaystyle {\dot {u}}_{l,\alpha }}$ at more than one energy parameter are used for the augmentation.[7] In comparison to a pure LAPW basis these approaches can precisely represent Kohn-Sham orbitals in a broader energy window around the energy parameters. The drawback is that the stricter matching conditions lead to a stiffer basis set.
• Lower-dimensional systems: The partitioning of the unit cell can be extended to explicitly include semi-infinite vacuum regions with their own augmentations of the plane waves.[38] This enables efficient calculations for lower-dimensional systems such as surfaces and thin films. For the treatment of atomic chains an extension to one-dimensional setups has been formulated.[39]

There are various software projects implementing the LAPW method and/or its variants. Examples for such codes are

• Elk
• Exciting
• Flair[40]
• FLEUR
• HiLAPW
• Wien2k