Construction of the charge density

The partitioning of the unit cell into different regions does not only affect the LAPW basis but also the representation of the charge density which also becomes region-dependent. In detail the charge density is given by

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \rho(\vec{r})=\left\lbrace \begin{array}{l l l} \sum\limits_{\vec{G}} \rho_{\vec{G}} e^{i\vec{G}\vec{r}} & \text{for } \vec{r} \in \text{IR} \\ \sum\limits_{L} \rho_L^\alpha(r_\alpha) Y_L(\vec{\hat{r}_\alpha}) & \text{for } \vec{r} \in \text{MT}^\alpha\\ \sum\limits_{\vec{G}_\parallel} \rho_{\vec{G}_\parallel}^\text{vac}(z) e^{i\vec{G}_\parallel \vec{r}_\parallel} & \text{for } \vec{r} \in \text{VR}^\text{vac} \\ \end{array} \right.. \end{equation}$$

In the interstitial region it is given as a plane-wave expansion, in the MT spheres as a spherical harmonics expansion, and in the vacuum region as one-dimensional functions times plane waves parallel to the respective thin film.

The charge density is constructed by calculating

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \rho(\vec{r}) = 2 \sum\limits_{\vec{k}} \sum\limits_{\nu} \omega_{\vec{k}}^{\nu} \left\vert \psi_{\vec{k}}^{\nu}(\vec{r}) \right\vert^2, \end{equation}$$

where $\psi_{\vec{k}}^{\nu}(\vec{r})$ is the Kohn-Sham wave function for eigenvalue $\nu$ at k point $\vec{k}$ and $\omega_{\vec{k}}^{\nu}$ is the occupation of the respective state times the weight of the k point. The factor 2 is due to the spin-degeneracy in non-spinpolarized systems. For magnetic calculations it is replaced by a sum over the spins.

Only parts of the charge density construction are discussed below. For the valence electrons the determination of the occupation numbers is discussed. It is important to users since it depends on input parameters. In the construction of the core electron density contribution several minor approximations come into play. These are also discussed below.

Occupation numbers for the valence electrons

The occupation numbers $\omega_{\vec{k}}^{\nu}$ for the valence electrons are calculated by occupying the states according to a Fermi distribution. In a first step this approach is used together with the predefined number of valence electrons in the unit cell $N_\text{elec}$ to determine the Fermi energy $E_\text{F}$ by enforcing

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} N_\text{elec} = 2 \sum\limits_{\vec{k}} \sum\limits_{\nu} \omega(\vec{k}) \frac{1}{e^{(\epsilon_{\vec{k}}^{\nu}-E_\text{F})/k_\text{B}T} + 1}, \end{equation}$$

where the factor 2 is again due to spin-degeneracy, $\omega(\vec{k})$ is the weight of the respective k point, $\epsilon_{\vec{k}}^{\nu}$ is the eigenenergy of the respective state, and $T$ is the temperature used for the definition of the Fermi distribution.

After determining the Fermi energy the occupation numbers are calculated as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \omega_{\vec{k}}^{\nu} = \omega(\vec{k}) \frac{1}{e^{(\epsilon_{\vec{k}}^{\nu}-E_\text{F})/k_\text{B}T} + 1}. \end{equation}$$

Density contributions from the core electrons

Besides the density due to the valence electrons there are also density contributions due to the core electrons. It is constructed by occupying the core electron states and summing the densities related to each state up for each of the atoms.

One should be aware of some details of the core electron density construction. Approximately the core electrons are confined within the MT spheres, but for a rigorous treatment one also has to consider the tails of the core electron states reaching out of the MT spheres.

In an exact treatment the core electron density can be represented similarly to the overall charge density as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \rho_\text{c}(\vec{r})=\left\lbrace \begin{array}{l l l} \sum\limits_{\vec{G}} \rho_{\text{c}\vec{G}} e^{i\vec{G}\vec{r}} & \text{for } \vec{r} \in \text{IR} \\ \sum\limits_{L} \rho_{\text{c}L}^\alpha(r_\alpha) Y_L(\vec{\hat{r}_\alpha}) & \text{for } \vec{r} \in \text{MT}^\alpha\\ \sum\limits_{\vec{G}_\parallel} \rho_{\text{c}\vec{G}_\parallel}^\text{vac}(z) e^{i\vec{G}_\parallel \vec{r}_\parallel} & \text{for } \vec{r} \in \text{VR}^\text{vac} \\ \end{array} \right.. \end{equation}$$

It is constructed by first considering the (spherical) internal core electron density originating from a certain atom $\rho_\text{c}^{\alpha,\text{ int}}(r_\alpha)$. From this one constructs the respective interstitial contribution by replacing the MT part of $\rho_\text{c}^{\alpha,\text{ int}}(r_\alpha)$ by a smooth pseudodensity and then performing a plane-wave expansion of the resulting function. In detail one defines the pseudodensity as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \tilde{\rho}_\text{c}^\alpha(r_\alpha)=\left\lbrace \begin{array}{l l l} a_\alpha e^{-b_\alpha r_\alpha^2} & \text{for } r_\alpha < R_{\text{MT}^\alpha} \\ \rho_\text{c}^{\alpha,\text{ int}}(r_\alpha) & \text{for } r_\alpha \ge R_{\text{MT}^\alpha} \\ \end{array} \right., \end{equation}$$

where $a_\alpha$ and $b_\alpha$ are determined by enforcing a continuous value and slope of the pseudo charge density. The respective plane-wave expansion coefficients are then calculated as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \rho_{\text{c}\vec{G}}^\alpha = \sqrt{\frac{4\pi}{\Omega}} e^{i\vec{G}\vec{\tau}_\alpha} \int\limits_0^\infty \tilde{\rho}_\text{c}^\alpha(r_\alpha) j_0(Gr_\alpha) r_\alpha^2 dr_\alpha \end{equation}$$

and the overall core density interstitial plane-wave expansion coefficients become

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \rho_{\text{c}\vec{G}} = \sum\limits_\alpha \rho_{\text{c}\vec{G}}^\alpha. \end{equation}$$

Of course, the interstitial plane-wave expansion stretches over the whole unit cell and core electron tails originating from a certain atom may reach into the MT sphere of another atom. These contributions also have to be added. To do so one re-expands the interstitial core electron density in terms of spherical harmonics in each MT sphere. This defines the external core electron density

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \rho_{\text{c}L}^{\alpha,\text{ ext}}(r_\alpha) = \frac{4\pi}{\sqrt{\Omega}} \sum\limits_{\vec{G}} \rho_{\text{c}\vec{G}} i^l e^{i\vec{G}\vec{\tau}_\alpha} j_l(Gr_\alpha) Y_L^*({\hat{\vec{G}}}). \end{equation}$$

The overall core electron density in the MT sphere of atom $\alpha$ then becomes

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \sum_L \rho_{\text{c}L}^\alpha(r_\alpha) Y_L(\hat{\vec{r}}_\alpha) = \left[\rho_\text{c}^{\alpha,\text{ int}}(r_\alpha) - a_\alpha e^{-b_\alpha r_\alpha^2}\right] Y_{l=0,m=0}(\hat{\vec{r}}_\alpha) + \sum\limits_L \rho_{\text{c}L}^{\alpha,\text{ ext}}(r_\alpha) Y_L(\hat{\vec{r}}_\alpha). \end{equation}$$

Note that the re-expansion of the interstitial core electron density in terms of spherical harmonics is computationally expensive and it is questionable whether a large computational effort should be accepted since only a very small charge is affected by this calculation. Therefore one limits the $L$ for this expansion by a new cutoff parameter $l_\text{max}^\text{coretail}$.

For the vacuum regions one assumes that the core electron density in these regions originates from a narrow energy range which lies far below the vacuum potential at infinity. Under this assumption it is resonable to approximate the density in the vacuum by an exponentially decaying function. In detail one constructs the vacuum density expansion coefficients as

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} \rho_{\text{c}\vec{G}_\parallel}^\text{vac}(z) = \rho_{\text{c}\vec{G}_\parallel}^\text{vac}(z^\text{vac}) e^{-\frac{z^\text{vac}}{\vert z^\text{vac}\vert} \alpha_{\vec{G}_\parallel}^\text{vac} (z-z^\text{vac})}. \end{equation}$$

In this equation $z^\text{vac}$ is the z coordinate of the vacuum boundary. The coefficients $\rho_{\text{c}\vec{G}_\parallel}^\text{vac}(z^\text{vac})$ and $\alpha_{\vec{G}_\parallel}^\text{vac}$ are again determined by enforcing continuity of value and slope of the density.