Details on the Hamiltonian and Overlap matrix setup

Due to the partitioning of the unit cell into different regions it is natural to express the Hamiltonian and Overlap matrices in terms of sums of contributions originating from each of these regions. Also the construction of the radial functions in the MT spheres on the basis of the spherical part of the effective potential leads to a seperataion of spherical and nonspherical contributions in the spheres. In consequence the Hamiltonian is expressed as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} H_{\vec{G},\vec{G'}}^{\vec{k}} = H_{\vec{G},\vec{G'}}^{\vec{k},\text{IR}} + \left[\sum_\alpha H_{\vec{G},\vec{G'}}^{\vec{k},\alpha,\text{sphr}} + H_{\vec{G},\vec{G'}}^{\vec{k},\alpha,\text{nsphr}} \right] + \left[\sum_\text{vac} H_{\vec{G},\vec{G'}}^{\vec{k},\text{vac}} \right] \end{equation}$$

and the Overlap matrix as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} S_{\vec{G},\vec{G'}}^{\vec{k}} = S_{\vec{G},\vec{G'}}^{\vec{k},\text{IR}} + \left[\sum_\alpha S_{\vec{G},\vec{G'}}^{\vec{k},\alpha} \right] + \left[\sum_\text{vac} S_{\vec{G},\vec{G'}}^{\vec{k},\text{vac}} \right]. \end{equation}$$

In the following the basic setup of these matrices is sketched. Since this is mainly intended to introduce cutoff parameters and other quantities that are important for the usage of the FLAPW method some contributions to the matrices will be neglected. For example the contributions related to local orbitals are not discussed here. These additional basis functions would increase the matrix size. Also the vacuum contributions will be ignored.

Interstitial contributions

For the construction of the interstitial contributions to the matrices one defines the step function

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \Theta(\vec{r}) = \left\{ \begin{array}{l l} 1 & \text{for } \vec{r} \in \text{IR} \\0 & \text{for } \vec{r} \in \text{MT}^\alpha\\ \end{array}\right.. \end{equation}$$

It is analytically constructed in reciprocal space as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \Theta(\vec{G}) = \delta_{\vec{G},\vec{0}} - \sum\limits_\alpha e^{-i\vec{G}\vec{\tau}_\alpha} \frac{4\pi (R_{\text{MT}^\alpha})^3}{\Omega} \frac{j_1(G R_{\text{MT}^\alpha})}{G R_{\text{MT}^\alpha}} \end{equation}$$

up to the reciprocal cutoff parameter $G_\text{max}$. $j_1$ is the spherical Bessel function for $l=1$ and $\Omega$ is the unit cell volume. The step function is the indicator function for the interstitial region but due to the finite reciprocal cutoff parameter its real space representation obtained by a Fourier transformation is only approximately correct. Fortunately only the correctness in reciprocal space up to the cutoff parameter is needed and this is fulfilled.

With the definition of the step function the interstitial contributions to the matrices become

$$\begin{align} \renewcommand{\vec}[1]{\mathbf{#1}} H_{\vec{G}\vec{G}'}^{\vec{k},\text{IR}} & = \left\langle \phi_{\vec{k}\vec{G}} \middle\vert \hat{H} \middle\vert \phi_{\vec{k}\vec{G}'} \right\rangle_\text{IR} \nonumber \\ & = \left\langle \tilde{\phi}_{\vec{k}\vec{G}} \middle\vert V(\vec{r})\Theta(\vec{r}) -\frac{1}{2}\Theta(\vec{r}) \nabla^2 \middle\vert \tilde{\phi}_{\vec{k}\vec{G}'} \right\rangle_\Omega \nonumber \\ & = \frac{1}{\Omega} \int\limits_\Omega e^{-i(\vec{G}-\vec{G}')\vec{r}} V(\vec{r}) \Theta(\vec{r}) d^3r + \frac{\left( \vec{k} + \vec{G}'\right)^2}{2\Omega} \int\limits_\Omega e^{-i(\vec{G}-\vec{G}')\vec{r}} \Theta(\vec{r}) d^3r \\ & = \left(V\Theta\right)(\vec{G}-\vec{G}') + \frac{\left( \vec{k} + \vec{G}'\right)^2}{2} \Theta(\vec{G}-\vec{G}') \end{align}$$

and

$$\begin{align} \renewcommand{\vec}[1]{\mathbf{#1}} S_{\vec{G}\vec{G}'}^{\vec{k},\text{IR}} & = \left\langle \tilde{\phi}_{\vec{k}\vec{G}} \middle\vert \Theta(\vec{r}) \middle\vert \tilde{\phi}_{\vec{k}\vec{G}'} \right\rangle_\Omega \nonumber = \frac{1}{\Omega} \int\limits_\Omega e^{-i(\vec{G}-\vec{G}')\vec{r}} \Theta(\vec{r}) d^3r \\ & = \Theta(\vec{G}-\vec{G}') \end{align}$$

The reciprocal plane wave cutoff for the step function $G_\text{max}$ is also used for the function $\left(V\Theta\right)(\vec{G})$.

Note that the expression for the interstitial contributions to the Hamiltonian does not lead to a Hermitian matrix. Since this property is required the actual expression to calculate this contribution is altered to guarantee it. It becomes

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} H_{\vec{G}\vec{G}'}^{\vec{k},\text{IR}} = \left(V\Theta\right)(\vec{G}-\vec{G}') + \frac{\left( \vec{k} + \vec{G}\right)^2 + \left( \vec{k} + \vec{G}'\right)^2}{4} \Theta(\vec{G}-\vec{G}'). \end{equation}$$

Spherical MT contributions

Since the radial functions are already constructed for the spherical potential in each MT sphere the respective contributions to the Hamiltonian matrix can easily be expressed as

$$\begin{align} \renewcommand{\vec}[1]{\mathbf{#1}} H_{\vec{G}\vec{G}'}^{\vec{k},\alpha,\text{sphr}} = \sum\limits_{l=0}^{l_\text{max}^\alpha} \sum\limits_{m=-l}^{l} E_l \left[ \left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast a_{\vec{k}\vec{G}'}^{L\alpha} + \langle \dot{u}_l^\alpha \vert \dot{u}_l^\alpha \rangle \left(b_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast b_{\vec{k}\vec{G}'}^{L\alpha} \right] + \left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast b_{\vec{k}\vec{G}'}^{L\alpha}. \end{align}$$

Accordingly the MT contribution to the Overlap matrix becomes

$$\begin{align} \renewcommand{\vec}[1]{\mathbf{#1}} S_{\vec{G}\vec{G}'}^{\vec{k},\alpha} = \sum\limits_{l=0}^{l_\text{max}^\alpha} \sum\limits_{m=-l}^{l} \left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast a_{\vec{k}\vec{G}'}^{L\alpha} + \langle \dot{u}_l^\alpha \vert \dot{u}_l^\alpha \rangle \left(b_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast b_{\vec{k}\vec{G}'}^{L\alpha}. \end{align}$$

Note that the expression for the Hamiltonian contributions again does not lead to a Hermitian matrix. Therefore it is altered to

$$\begin{align} \renewcommand{\vec}[1]{\mathbf{#1}} H_{\vec{G}\vec{G}'}^{\vec{k},\alpha,\text{sphr}} = \sum\limits_{l=0}^{l_\text{max}^\alpha} \sum\limits_{m=-l}^{l} E_l \left[ \left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast a_{\vec{k}\vec{G}'}^{L\alpha} + \langle \dot{u}_l^\alpha \vert \dot{u}_l^\alpha \rangle \left(b_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast b_{\vec{k}\vec{G}'}^{L\alpha} \right] + \frac{1}{2}\left(\left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast b_{\vec{k}\vec{G}'}^{L\alpha} + \left(b_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast a_{\vec{k}\vec{G}'}^{L\alpha} \right). \end{align}$$

The summation over the magnetic quantum numbers $m$ can actually be performed analytically. This makes the calculation of the spherical MT contributions efficient.

Nonspherical contributions in the MT spheres

For the calculation of the contribution due to the nonspherical part of the potential in a MT sphere one first defines a local Hamiltonian $t_{LL'}^{\alpha[\circ][\diamond]}$ in the basis of the radial functions $u_{l}^\alpha(r_\alpha,E_{l}^{\alpha})$ and $\dot{u}_{l}^\alpha(r_\alpha,E_{l}^{\alpha})$ as

$$\begin{align} \renewcommand{\vec}[1]{\mathbf{#1}} t_{LL'}^{\alpha[\circ][\diamond]} & = \int\limits_{\text{MT}^\alpha} \sum\limits_{L''=1}^{l_\text{max}^\alpha} \overset{[\circ]}{u}{}_{l}^\alpha(r_\alpha,E_l^\alpha) V_{L''}^\alpha(r_\alpha) \overset{[\diamond]}{u}{}_{l'}^\alpha(r_\alpha,E_{l'}^\alpha) Y_L^\ast(\vec{\hat{r}_\alpha}) Y_{L''}(\vec{\hat{r}}_\alpha) Y_{L'}(\vec{\hat{r}_\alpha}) d^3r \nonumber \\ & = \int\limits_0^{R_{\text{MT}^\alpha}} \sum\limits_{L''=1}^{l_\text{max}^\alpha} r_\alpha^2 \overset{[\circ]}{u}{}_{l}^\alpha(r_\alpha,E_l^\alpha) V_{L''}^\alpha(r_\alpha) \overset{[\diamond]}{u}{}_{l'}^\alpha(r_\alpha,E_{l'}^\alpha) G_{ll''l'}^{mm''m'} dr_\alpha. \end{align}$$

In this expression $\circ$ and $\diamond$ can each stand for a $u_{l}^\alpha(r_\alpha,E_{l}^{\alpha})$ or a $\dot{u}_{l}^\alpha(r_\alpha,E_{l}^{\alpha})$ and $V_{L''}^\alpha(r_\alpha)$ is the potential in the sphere. The Gaunt coefficients $G_{ll''l'}^{mm''m'}$ are defined as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} G_{ll''l'}^{mm''m'} = \int\limits_0^{2\pi} \int\limits_0^{\pi} Y_L^\ast(\theta,\phi) Y_{L''}(\theta,\phi) Y_{L'}(\theta,\phi) \mathrm{sin}(\theta) d\theta d\varphi. \end{equation}$$

The local Hamiltonian then has to be incorporated into the global Hamiltonian for the whole unit cell by applying the matching coefficients. One obtains

$$\begin{align*} \renewcommand{\vec}[1]{\mathbf{#1}} H_{\vec{G}\vec{G}'}^{\vec{k},\alpha,\text{nsphr}} & = \sum\limits_L \sum\limits_{L'} \left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast t_{LL'}^{\alpha[\phantom{\cdot}][\phantom{\cdot}]} a_{\vec{k}\vec{G'}}^{L'\alpha} + \left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast t_{LL'}^{\alpha[\phantom{\cdot}][\cdot]} b_{\vec{k}\vec{G'}}^{L'\alpha} \nonumber \\ & \phantom{=} + \left(b_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast t_{LL'}^{\alpha[\cdot][\phantom{\cdot}]} a_{\vec{k}\vec{G'}}^{L'\alpha} + \left(b_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast t_{LL'}^{\alpha[\cdot][\cdot]} b_{\vec{k}\vec{G'}}^{L'\alpha} \\ & = \sum\limits_L \left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast \left(\sum\limits_{L'}t_{LL'}^{\alpha[\phantom{\cdot}][\phantom{\cdot}]} a_{\vec{k}\vec{G'}}^{L'\alpha}\right) + \left(a_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast \left(\sum\limits_{L'}t_{LL'}^{\alpha[\phantom{\cdot}][\cdot]} b_{\vec{k}\vec{G'}}^{L'\alpha}\right) \nonumber \\ & \phantom{=} + \left(b_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast \left(\sum\limits_{L'}t_{LL'}^{\alpha[\cdot][\phantom{\cdot}]} a_{\vec{k}\vec{G'}}^{L'\alpha}\right) + \left(b_{\vec{k}\vec{G}}^{L\alpha}\right)^\ast \left(\sum\limits_{L'}t_{LL'}^{\alpha[\cdot][\cdot]} b_{\vec{k}\vec{G'}}^{L'\alpha}\right). \end{align*}$$

In these expressions the sums over $L$ and $L'$ are actually limited by the reduced $l$ cutoff parameter $l_\text{nonsph}^\alpha$. While in the calculation for the spherical contributions very high $l$ cutoff parameters were required to obtain reasonable kinetic energies and Overlap matrices, the sums involved for the nonspherical contributions are related to the potential only. In practice this implies that smaller cutoffs can be used here. Note that the reduction of the $l$ cutoffs reduces the computational demands significantly.