Linearized augmented-plane-wave method

The linearized augmented-plane-wave method (LAPW) is an implementation of Kohn-Sham density functional theory (DFT) adapted to periodic materials. It allows to consider the full potential without any shape approximations and to accurately treat valence as well as core electrons within the DFT framework. It does not rely on the pseudopotential approximation and employs a systematically extendable basis set. With these characteristics it belongs to the most precise implementations of DFT, applicable to all crystalline materials, irrespective of their chemical composition. It can be used as a reference to evaluate other approaches.

Introduction

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 is sketched in more detail and variations and extensions are indicated.

Solving the Kohn-Sham equations

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

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

with the single-electron kinetic energy operator ${\hat {T}}_{s}$, the effective potential $V_{\text{eff}}(\mathbf {r} )$, Kohn-Sham states $\Psi _{j}^{\mathbf {k} }(\mathbf {r} )$, energy eigenvalues $E_{j}$, and position and Bloch vectors $\mathbf {r}$ and $\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, solving the Kohn-Sham equations in practice is accompanied by the introduction of many additional approximations. This includes the incompleteness of the basis set used to represent the Kohn-Sham orbitals, the choice whether the pseudopotential approximation is used or all electrons are considered in the DFT scheme, the treatment of relativistic effects, and possible shape approximations to the potential. These different method design decisions are intertwined. Beyond the partitioning of the unit cell, for the all-electron full-potential LAPW method the central design aspect is the usage of the LAPW basis set $\left\lbrace \phi _{\mathbf {k} ,\mathbf {G} }(\mathbf {r} )\right\rbrace$ to represent the valence electron orbitals as

$\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 $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 $\Omega$ covering atoms $\alpha$ at positions $\mathbf {\tau } _{\alpha }$, an LAPW basis function is characterized by a reciprocal lattice vector $\mathbf {G}$ and the considered Bloch vector $\mathbf {k}$. It is given as

$\phi _{\mathbf {k} ,\mathbf {G} }(\mathbf {r} )=\left\lbrace {\begin{array}{l l}{\frac {1}{\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}({\hat {r}}_{\alpha })&{\text{for }}{\vec {r}}{\text{ in MT}}_{\alpha }\end{array}}\right.,$

where $\mathbf {r} _{\alpha }=\mathbf {r} -\mathbf {\tau } _{\alpha }$ is the position vector relative to atom nucleus $\alpha$. An LAPW basis function is thus a plane wave in the IR and a linear combination of the radial functions $u_{l,\alpha }(r_{\alpha },E_{l,\alpha })$ and ${\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })$ multiplied by spherical harmonics $Y_{l,m}$ in each MT sphere. The radial function $u_{l,\alpha }(r_{\alpha },E_{l,\alpha })$ hereby is the solution to the spherically averaged problem for the predetermined energy parameter $E_{l,\alpha }$ with regular behavior at the nucleus. Together with its energy derivative ${\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 $a_{l,m}^{\mathbf {k} ,\mathbf {G} ,\alpha }$ and $b_{l,m}^{\mathbf {k} ,\mathbf {G} ,\alpha }$ are automatically determined by enforcing the basis function to be continuously differentiable in the respective $(l,m)$ channel. The set of LAPW basis functions is defined by specifying a cutoff parameter $K_{\text{max}}=|{\vec {k}}+{\vec {G}}|_{\text{max}}$. In each MT sphere the spherical harmonics expansion is limited by an angular momentum cutoff $l_{{\text{max}},\alpha }$.

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. 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). These are tailored to provide a precise representation of the semicore states.

The considered degree of relativistic physics differs for core and valence electrons. The strong localization of core electrons is connected to large kinetic energy contributions and thus a fully relativistic treatment is desirable and common. For the determination of the radial functions $u_{l,\alpha }(r_{\alpha },E_{l,\alpha })$ and ${\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 (SRA) or similar approaches. The dominant effect neglected by these approximations is spin-orbit coupling. The construction of the Hamiltonian matrix within the respective approximation is trivial because the energy parameters are strongly connected to the application of the Hamiltonian to the radial functions. Spin-orbit coupling can additionally be included, though this leads to a more complex Hamiltonian, connected to increased computational demands. In the interstitial region it is reasonable and common to describe the valence electrons without considering relativistic effects.

Beyond the treatment of spin-orbit coupling also the inclusion of nonspherical parts of the potential is not directly covered by the construction of the radial functions. For the full-potential description (FLAPW) this has to be additionally included for each MT sphere in the setup of the Hamiltonian matrix.

Representation of the charge density and the potential

After calculating the Kohn-Sham eigenfunctions the next step is the construction of the charge density by occupying the energetically lowest eigenstates up to the Fermi level with electrons. Most commonly the occupation hereby follows the Fermi-Dirac distribution. The Fermi level itself is determined in this process by keeping the number of electrons fixed. The resulting charge density $\rho (\mathbf {r} )$ then has a region-specific form

$\rho (\mathbf {r} )=\left\{{\begin{array}{l l}\sum \limits _{\mathbf {G} }\rho _{\mathbf {G} }e^{i{\vec {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}({\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 of 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. A notable challenge in the construction of the potential is the divergence of individual potential contributions like the Hartree potential or the external potential in infinite periodic crystals. To overcome this challenge a common approach to calculate the full potential in a reliable way without the appearance of divergent sums is employing Weinert's method for solving the Poisson equation.

Variants and extensions of the LAPW method

• APW: The augmented-plane-wave method is the predecessor of LAPW. It only uses the solution to the spherically averaged problem for the augmentation in the MT spheres. The energy derivative of this 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. In comparison to LAPW this is a more complex problem to solve.
• Local orbitals extensions: The LAPW basis can be extended by local orbitals (LOs). These are additional basis functions having nonvanishing values only in a single MT sphere. They are composed of the radial functions $u_{l,\alpha }(r_{\alpha },E_{l,\alpha })$, ${\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 or the elimination of the linearization error for the valence states.
• APW+lo: In the APW+lo method the augmentation in the MT spheres only consists of the function $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 ${\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 approach the APW+lo method leads to less stiff basis functions. 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 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 $u_{l,\alpha }(r_{\alpha },E_{l,\alpha })$ and ${\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha })$. This yields basis functions continuously differentiable also in the $(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 pairs of local orbitals introducing the functions $u_{l,\alpha }(r_{\alpha },E_{l,\alpha }^{\text{lo}})$ and ${\dot {u}}_{l,\alpha }(r_{\alpha },E_{l,\alpha }^{\text{lo}})$ are added to the LAPW basis. The energy parameters $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 the augmentation in the MT spheres additionally includes the second energy derivative ${\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. In comparison to a pure LAPW basis this approach can precisely represent Kohn-Sham orbitals in a boader 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 specify vacuum regions with own augmentations of the plane waves. This enables efficient calculations for lower-dimensional systems. To describe thin films and surfaces an extension of the LAPW method to two-dimensional systems has been implemented. For the treatment of atomic chains an extension to one-dimensional setups has been formulated.

Software implementations

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

• Elk
• Exciting
• Fleur
• Wien2k