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

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 plane wave expansion of the density in the interstitial region and the vacuum regions is limited by the $G_\text{max}$ cutoff specified in calculationSetup/cutoffs/@Gmax. The spherical harmonics expansion in the MT spheres is limited by the $l_\text{max}^\alpha$ cutoffs specified for each species in atomSpecies/species/atomicCutoffs/@lmax.

The charge density is constructed by calculating

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

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.

In the Fleur input file the temperature of the Fermi distribution is typically specified in terms of the related energy in calculationSetup/bzIntegration/@fermiSmearingEnergy. The number of valence electrons is specified in calculationSetup/bzIntegration/@valenceElectrons.

After determining the Fermi energy the occupation numbers are calculated as

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.

The number of core electron states is specified for each species separately in atomSpecies/species/@coreStates. By default all core electron states are fully occupied. Of course, for the calculation of core holes these occupations can be changed by explicitely providing an electron configuration in the optional element atomSpecies/species/electronConfig.

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.

A rigorous treatment of the core electron density is enabled if the calculationSetup/coreElectrons/@ctail switch is set to T. Otherwise the charge stored in the tails of the core electron states is just evenly distributed over the interstitial region.

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

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

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

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

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

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

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}$.

In the Fleur input file the $l_\text{max}^\text{coretail}$ cutoff parameter is specified in the optional XML attribute calculationSetup/coreElectrons/@coretail_lmax. If for a certain atom species this value is larger than the respective $l_\text{max}^\alpha$ cutoff it is reduced to this cutoff for this atom species. If the $l_\text{max}^\text{coretail}$ parameter is not set in the input file its default value is 99.

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

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.