Purpose: In this tutorial, you will learn what is a species file and what is the physics behind it.
Table of Contents
exciting, as a code based on plane-wave augmentation, distinguishes between the interstitial and muffin-tin regions. In the interstitial region, various quantities, such as wavefunctions, potential, and density are described by means of plane waves. The quality (completeness) of this basis is determined by the maximum wavevector, defined via either gmaxvr or rgkmax. Basis functions in the muffin-tin region are much more cumbersome to define and, thus, require a detailed description. Furthermore, each muffin tin is associated with an atom with a specific nuclear charge, mass, etc. All this information is summarized in species files.
Consider a standard species file for carbon.
<?xml version="1.0" encoding="UTF-8"?>
<spdb xsi:noNamespaceSchemaLocation="../../xml/species.xsd"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<sp chemicalSymbol="C" name="carbon" z="-6.00000" mass="21894.16673">
<muffinTin rmin="0.100000E-04" radius="1.4500" rinf="21.0932" radialmeshPoints="250"/>
<atomicState n="1" l="0" kappa="1" occ="2.00000" core="true"/>
<atomicState n="2" l="0" kappa="1" occ="2.00000" core="false"/>
<atomicState n="2" l="1" kappa="1" occ="1.00000" core="false"/>
<atomicState n="2" l="1" kappa="2" occ="1.00000" core="false"/>
<basis>
<default type="lapw" trialEnergy="0.1500" searchE="false"/>
<custom l="0" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
<custom l="1" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
</basis>
</sp>
</spdb>
The information contained in the first two lines is not relevant for the physics of the system. They are only used to give advanced XML instructions to define the XML language, style sheet, and schema. For further details on this topic, see XML and exciting data format. Considering only the physical background of this file and dividing it into segments by adding irrelevant blank lines, the file shown above can be simplified as follows.
<spdb>
<sp chemicalSymbol="C" name="carbon" z="-6.00000" mass="21894.16673">
<muffinTin rmin="0.100000E-04" radius="1.4500" rinf="21.0932" radialmeshPoints="250"/>
<atomicState n="1" l="0" kappa="1" occ="2.00000" core="true"/>
<atomicState n="2" l="0" kappa="1" occ="2.00000" core="false"/>
<atomicState n="2" l="1" kappa="1" occ="1.00000" core="false"/>
<atomicState n="2" l="1" kappa="2" occ="1.00000" core="false"/>
<basis>
<default type="lapw" trialEnergy="0.1500" searchE="false"/>
<custom l="0" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
<custom l="1" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
</basis>
</sp>
</spdb>
This way, we can identify four logical segments.
<sp chemicalSymbol="C" name="carbon" z="-6.00000" mass="21894.16673">
rmin stands for the radial coordinate of the innermost grid point. radius represents the default value of the muffin-tin radius. Keep in mind that it may be, and frequently is, overridden in the input file. rinf is the coordinate of the effective infinity. This parameter is used in calculations of atoms during the initialization and every time core states are calculated during the self-consistency steps. Finally, radialmeshPoints defines the number of grid points in the muffin tin. In practice, there will be more points between the muffin-tin boundary and the effective infinity, but their number is calculated automatically. Even though these attributes are not further discussed here, you should be aware that the default choice is not necessarily optimal. The choice of the muffin-tin radius clearly depends on inter-atomic distances and is system-dependent. Additionally, rmin and radialmeshPoints are parameters that you might want to explore occasionally.<muffinTin rmin="0.100000E-04" radius="1.4500" rinf="21.0932" radialmeshPoints="250"/>
kappa in exciting differs from that in the textbooks of quantum mechanics where it is defined as $\kappa = (l - j)(2j + 1)$. Here, kappa adopts the absolute value of the "textbook $\kappa$". According to the standard species file, the shell $1s^2$ is assigned to the core, while $2s^22p^2$ are considered as the valence shells.<atomicState n="1" l="0" kappa="1" occ="2.00000" core="true"/>
<atomicState n="2" l="0" kappa="1" occ="2.00000" core="false"/>
<atomicState n="2" l="1" kappa="1" occ="1.00000" core="false"/>
<atomicState n="2" l="1" kappa="2" occ="1.00000" core="false"/>
trialEnergy defines initial values of the energy parameters which are kept constant during all self consistent steps. It is possible to adapt them during the calculations by setting the parameter searchE to true. However, this is not recommended since it can lead to unstable calculations. The energy parameters are always considered on the absolute scale, and not as the distance to the Fermi energy.<basis>
<default type="lapw" trialEnergy="0.1500" searchE="false"/>
<custom l="0" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
<custom l="1" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
</basis>
The muffin-tin part of the augmented-plane-wave (APW) basis is
\begin{equation} \phi_{\textbf{G+k}}(\textbf{r})=\sum_{lm}A^{\textbf{G+k}}_{lm}u_l(r;\varepsilon_l)Y_{lm}(\hat{\textbf{r}}) \; . \end{equation}In order to make APW accurate, the energy parameter $ε_l$ must be a very good approximation to the eigenenergies of states with angular momentum $l$. It is difficult or even impossible to arrange this in practice. To overcome this problem, the energy dependence of the basis is linearized. The radial function with the exact eigenenergy as the energy parameter ϵ is expressed as
\begin{equation} u_l(r;\epsilon)=u_l(r;\varepsilon_l)+ (\epsilon-\varepsilon_l)\dot{u_l}(r;\varepsilon_l) \; , \end{equation}where $\dot{u_l}=\partial{u_l}/\partial{\epsilon}$. This expression cannot be applied directly, however, it shows that the improved muffin-tin basis should contain $u_l(r;\varepsilon_l)$ and $\dot{u_l}(r;\varepsilon_l)$.
One possibility how to implement the linearization are linearized augmented plane waves (LAPW). According to LAPW, the basis functions in muffin tins are defined as
\begin{equation} \phi_{\textbf{G+k}}(\textbf{r})=\sum_{lm}\left[ A^{\textbf{G+k}}_{lm}u_l(r;\varepsilon_l) + B^{\textbf{G+k}}_{lm}\dot{u_l}(r;\varepsilon_l) \right]Y_{lm}(\hat{\textbf{r}}) \;, \end{equation}where the prefactors $A^{\textbf{G+k}}_{lm}$ and $B^{\textbf{G+k}}$ are chosen in a such way that $\phi_{\textbf{G+k}}(\textbf{r})$ is smooth and continuous at the muffin-tin boundary. To define the smooth (LAPW) augmentation in species, set the attribute type to "lapw". In the example below, the species file for carbon is transformed in such a way that only the augmentation type LAPW is used.
<basis>
<default type="lapw" trialEnergy="0.15" searchE="false"/>
<custom l="0" type="lapw" trialEnergy="-0.20" searchE="false"/>
<custom l="1" type="lapw" trialEnergy="0.15" searchE="false"/>
</basis>
Here, the energy parameters used for LAPW are $\varepsilon_0=−0.20$, $\varepsilon_1=0.15$ and $\varepsilon_{l>1}=0.15$.
The alternative approach to implant the linearization is to introduce a different kind of basis functions - local orbitals. They are defined as
\begin{equation} \phi_{\nu}(\textbf{r})=\left[ a_{nu}u_l(r;\varepsilon_l) + b_{nu}\dot{u_l}(r;\varepsilon_l) \right]Y_{lm}(\hat{\textbf{r}}) \;, \end{equation}in one particular muffin-tin and as 0 everywhere else. The example below shows how to transform your basis into APW+lo.
<basis>
<default type="apw" trialEnergy="0.1500" searchE="false"/>
<custom l="0" type="apw+lo" trialEnergy="-0.20" searchE="false"/>
<custom l="1" type="apw+lo" trialEnergy="0.15" searchE="false"/>
</basis>
This specification of the basis involves APWs for all $l$-channels. However, for Here, APW+lo is used for $l=0$ and $l=1$ channels. Higher $l$ are considered using LAPW. there are local orbitals added to the basis. They are requested by specifying type = "apw+lo". These local orbitals have exactly the same form as introduced above, namely, they are linear combinations of the "primitive" functions $u_l(r;\varepsilon_l)$ and $\dot{u_l}(r;\varepsilon_l)$. The energy parameters are the same as in the LAPW case. In theory LOs can also consist of more than two "primitive" functions which can also be higher order derivatives. These are used to describe semi-core states which are explained in the section on semi-core states.
As you have seen in the example of the first section, APW+lo and LAPW can be combined.
<basis>
<default type="lapw" trialEnergy="0.1500" searchE="false"/>
<custom l="0" type="apw+lo" trialEnergy="-0.20" searchE="false"/>
<custom l="1" type="apw+lo" trialEnergy="0.15" searchE="false"/>
</basis>
Here, APW+lo is used for $l=0$ and $l=1$ channels. Higher $l$ are considered using LAPW.
In general, both APW+lo and LAPW significantly improve upon APW. While, in practice, APW+lo is more accurate than LAPW, the former methods introduces additional basis functions. Therefore, APW+lo is not practical for replacing APW for high-$l$ channels. For this reason, we recommend to use APW+lo for low $l$, and LAPW for high $l$. This strategy is implemented in the standard species files.
The linearization procedure described above is useful for improving the description of valence states. However, local orbitals can be used also for describing semi-core states.
Consider the lithium atom. It has the shell structure $1s^22s^1$, which contains just one valence electron. The $1s$ state can be attributed to the core, but it is neither sufficiently localized within Li muffin tins nor well separated from valence states in terms of energy. Therefore, it makes sense to consider the $1s$ state on the same footing as valence states. This idea is implemented in a species file for Li which you can find below.
<spdb>
<sp chemicalSymbol="Li" name="lithium" z="-3.00000" mass="12652.66897">
<muffinTin rmin="0.100000E-04" radius="1.7000" rinf="29.9495" radialmeshPoints="250"/>
<atomicState n="1" l="0" kappa="1" occ="2.00000" core="false"/>
<atomicState n="2" l="0" kappa="1" occ="1.00000" core="false"/>
<basis>
<default type="lapw" trialEnergy="0.1500" searchE="false"/>
<custom l="0" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
<lo l="0">
<wf matchingOrder="0" trialEnergy="0.1500" searchE="true"/>
<wf matchingOrder="1" trialEnergy="0.1500" searchE="true"/>
<wf matchingOrder="0" trialEnergy="-1.8784" searchE="true"/>
</lo>
</basis>
</sp>
</spdb>
The attribute core = "false" explicitly defines that $1s$ and $2s$ are supposed to be treated as valence states. Another important ingredient is the element that describes a local orbital.
<lo l="0">
<wf matchingOrder="0" trialEnergy="0.1500" searchE="true"/>
<wf matchingOrder="1" trialEnergy="0.1500" searchE="true"/>
<wf matchingOrder="0" trialEnergy="-1.8784" searchE="true"/>
</lo>
his element consists of three subelements that define three "primitive" functions that satisfy the radial Schrödinger equation. The first subelement contains the attributes matchingOrder = "0" and trialEnergy = "0.1500". They mean that the first "primitive" function corresponds to the zeroth-order energy-derivative and the energy parameter is $0.15$. In other words, the first subelement corresponds to $u_0(r;0.15)$. The second subelement contains matchingOrder = "1". It means that the second "primitive" function is the first energy-derivative $\dot{u_0}(r;0.15)$. Finally, the third function is $u_0(r;−1.8784)$. In theory also higher order derivatives can be used to define a local orbital.
The local orbital is a linear combination of these "primitive" functions:
\begin{equation} \phi_{\nu}(\textbf{r})=\left[ a_{nu}u_l(r;\varepsilon_l) + b_{nu}\dot{u_l}(r;\varepsilon_l) + c_{nu}u_l(r;\varepsilon_{lo}) \right]Y_{lm}(\hat{\textbf{r}}) \;, \end{equation}where $l$, $\varepsilon_l$, and $\varepsilon_{lo}$ are of the values listed above. The prefactors $a_{\nu}$, $b_{\nu}$, and $c_{\nu}$ are automatically chosen in a such way that the local orbital normalizes to 1 and turns to 0 at the muffin-tin boundary. These two conditions are essential and are employed for all local orbitals. In order to be able to determine $a_{\nu}$, $b_{\nu}$, and $c_{\nu}$ uniquely, the third condition is that the radial derivative of the local orbital turns to 0 at the muffin-tin boundary.
Choosing appropriate linearization energies for local orbitals can be challenging. To avoid this, linearization energies can be computed automatically. This is done by using the Wigner-Seitz rules. They provide an energy range in which the ideal linearization energy can be found.To get the upper and lower bound energy, the radial Schrödinger equation is solved for the corresponding principal quantum number ($n$), angular momentum ($l$), and initial potential for different energies.The upper bound energy is then equal to the energy for which the wave function (the solution of SE) is zero at the muffin tin radius.
\begin{equation} u_l(r; E_{\text{max}})\bigg | _{r=R_{\text{MT}}} = 0 \; , \end{equation}The lower bound energy is given by the energy for which the first radial derivative of the wave function vanishes at the muffin tin radius.
\begin{equation} \frac{\partial u_{l}(r, E_{\text{min}})}{\partial r} \bigg | _{r=R_{\text{MT}}} = 0 \; , \end{equation}In the end, the mean of these two energies is taken as the lineratization energy $\epsilon_{l} = (E_{\text{max}}+E_{\text{min}})/2$.
To utilize this algorithm, the species file must specify the principal quantum number rather than the linearization energy. For the example given above, the trialEnergy = "0.1500" corresponds to the $2s$ state and therefore needs to be replaced by n = "2". The lower energy trialEnergy = "-1.8784" corresponds to the $1s$ state and needs to be replaced by n = "1". The whole local orbital subelement then looks like the following.
<lo l="0">
<wf matchingOrder="0" n="2" searchE="false"/>
<wf matchingOrder="1" n="2" searchE="false"/>
<wf matchingOrder="0" n="1" searchE="false"/>
</lo>
The calculated linearization energies can be checked in the LINENGY.OUT file.
The standard local orbitals in exciting are scalar relativistic basis functions. While this is a good approximation for most materials, those with strong spin-orbit coupling (SOC) effects are not properly described using a scalar relativistic (sr) basis. Therefore, we introduce so called Dirac-type local orbitals.
The "primitive" functions of these Dirac-type local orbitals are calculated with the radial Dirac equation instead of the Schrödinger equation and, thus, depend on the total angular momentum $j$. To define them, we use the relativistic quantum number $\kappa$.
Dirac-type orbitals can be added to the basis, by including an lo subelement with the attribute kappa to the basis subelement as follows.
<lo l="1">
<wf matchingOrder="0" trialEnergy="-3.74" searchE="false" kappa="1"/>
<wf matchingOrder="1" trialEnergy="-3.74" searchE="false" kappa="1"/>
</lo>
Note that here, in contrast to the atomicState part of the species file, the actual value of $\kappa$ and not the absolute value is used.
In the above shown code block, we add a $p$-orbital with total angular momentum $j=1/2$ to the basis. This orbital is especially important since it differs strongly from the sr $p$-orbital. Unlike the sr one, it has a non-zero value at $r=0$ and therefore, cannot be represented in an sr basis.
Dirac type local orbitals should always be used in combination with second variation with local orbitals which can be enabled by setting the input parameter svlo to true.