6. Basis sets

6.1. LAPW basis

The reader is supposed to be familiar with the LAPW basis set. This section is intended as a short recapitulation and to define the notation used in the manual.

The LAPW method relies on a partitioning of space into non-overlapping MT spheres centered at the atomic nuclei and the interstitial region. The basis functions are defined differently in the two regions, plane waves in the interstitial and numerical functions in the spheres, consisting of a radial part \({u_{lp}^a(r)}\) (\({a}\) is the atomic index) and the spherical harmonics \({Y_{lm}(\hat{\mathbf{r}})}\). \({p}\) is an index to distinguish different types of radial functions: \({u_{l0}^a=u_l^a}\), \({\,u_{l1}^a=\dot{u}_l^a=\partial u_l^a/\partial\epsilon}\), \({\,u_{lp}^a=u_l^{a,\mathrm{LO}},\,p\ge 2}\), where we have used the usual notation of LAPW. Augmented plane waves are formed by matching linear combinations of the \({u}\) and \({\dot{u}}\) to the interstitial plane waves, forming the standard LAPW basis set. Local orbitals \({u^\mathrm{LO}}\) can be used to extend the basis set, to enable the description of semicore and high-lying conduction states. The plane waves and the angular part of the MT functions can be converged straightforwardly with the reciprocal cutoff radius \({g_\mathrm{max}}\) and the maximal l quantum number \({l_\mathrm{max}}\), respectively, whereas the radial part of the MT functions is more difficult to improve systematically, but it is possible, see [Comput. Phys. Commun. 184, 2670 (2013)].

6.1.1. MTACCUR

(*) The accuracy of the radial MT basis can be analyzed with the keyword MTACCUR followed by two numbers \(e_1\) and \(e_2\). Spex then calculates the MT representation error [Phys. Rev. B 83, 081101] in the energy range between \(e_1\) and \(e_2\). (If unspecified, upper and lower bounds are chosen automatically.) The results are written to the output files spex.mt.t where t is the atom type index, or spex.mt.s.t with the spin index s(=1 or 2) for spin-polarized calculations. The files contain sets of data for all l quantum numbers, which can be plotted separately with “gnuplot” (e.g., plot "spex.mt.1" i 3 for \({l=3}\)).

MTACCUR -1 2 Calculate MT representation error between -1 and 2 Hartree.
MTACCUR Automatic energy range.

6.2. Mixed product basis

The methods implemented in Spex distinguish themselves from conventional DFT with local or semilocal functionals by the fact that individual two-particle scattering processes are described explicitly. In such processes, the state of an incoming particle and the state that it scatters into form products of wavefunctions. Many of the physical quantities, such as the Coulomb interaction, the polarization function, and the dielectric function, become matrices when represented in a suitable product basis. In the context of the FLAPW method, the mixed product basis is an optimal choice. It consists of two separate sets of functions that are defined only in one of the spatial regions, while being zero in the other, the first is defined by products of the LAPW MT functions and the second by interstitial plane waves. (The products of plane waves are plane waves again.) The corresponding parameters are defined in the section MBASIS of the input file.

6.2.1. GCUT (MBASIS)

If \({N}\) is the number of LAPW basis functions, one would naively expect the number of product functions to be roughly \({N^2}\). In the case of the interstitial plane waves, this is not so, since, with a cutoff \({g_\mathrm{max}}\), the maximum momentum of the product would be \({2g_\mathrm{max}}\), leading to \({8N}\) as the number of product plane waves. Fortunately, it turns out that the basis size can be much smaller in practice. Therefore, we introduce a reciprocal cutoff radius \({G_\mathrm{max}}\) for the interstitial plane waves and find that, instead of \({G_\mathrm{max}=2g_\mathrm{max}}\), good convergence is achieved already with \({G_\mathrm{max}=0.75g_\mathrm{max}}\), the default value. The parameter \({G_\mathrm{max}}\) can be set to a different value with the keyword GCUT in the section MBASIS of the input file.

GCUT 2.9 Use a reciprocal cutoff radius of 2.9 Bohr-1 for the product plane waves.

6.2.2. LCUT (MBASIS)

The MT functions are constructed from the products \({u_{lp}^a(r)Y_{lm}(\hat{\mathbf{r}})u_{l'p'}^a(r)Y_{l'm'}(\hat{\mathbf{r}})}\). (Obviously, the atom index \({a}\) must be identical.) The product of spherical harmonics expands into a linear combination of \({Y_{LM}(\hat{\mathbf{r}})}\) with \({L=|l-l'|,...,l+l'}\) and \({|M|\le L}\). In principle, this leads to a cutoff of \({2l_\mathrm{max}}\) for the products, but, as with the plane waves, one can afford to use a much smaller cutoff parameter in practice. We use \({L_\mathrm{max}=l_\mathrm{max}/2}\) as the default. The corresponding keyword in the input file is called LCUT.

LCUT 5 4 Use l cutoffs of 5 and 4 for two atom types in the system.


In the LAPW basis, the matching conditions at the MT boundaries require large l cutoff values, typically around \({l=10}\), giving a correspondingly large number of radial functions. Not all of these functions must be considered in the product basis set. The keyword SELECT enables a specific choice of the radial functions. When specified, an entry for each atom type is required. The best way to describe the syntax of SELECT is with examples:

SELECT 2;3 Use products of \({u_{lp} u_{l'p'}}\) with \({l\le 2}\) and \({l'\le 3}\) for \({p=0}\) (so-called \({u}\)) and no function of \({p=1}\) (so-called \({\dot{u}}\)).
SELECT 2,1;3,2 Same as before but also include the \({\dot{u}}\) functions with \({l\le 1}\) and \({l'\le 2}\).
SELECT 2,,1100;3,,1111 Same as first line but also include the local orbitals \({p\ge 2}\), which are selected (deselected) by “1” (“0”): here, the first two and all four LOs, respectively. The default behavior is to include semicore LOs but to exclude the ones at higher energies.
SELECT 2,1,1100;3,2,1111 Same as second line with the LOs.

The LOs are selected by series of 1s and 0s. If there are many LOs, a definition like 0000111 can be abbreviated by 4031 (four times 0, three times 1).

One of the most important quantities to be calculated in GW calculations (and related methods) are the “projections” \({\langle M_{\mathbf{k}\mu} \phi_{\mathbf{q}n} | \phi_{\mathbf{k+q}n'} \rangle}\), which describe the coupling of a propagating particle to an interaction line. Often, the two states with \(n\) and \(n'\) are occupied and unoccupied, respectively. In this sense, the first set of SELECT parameters before the semicolon refers to the occupied states, the set after the semicolon refers to the unoccupied states. For \(l=2\) with \(2l\) or \(2l+1\) being the l cutoff specified by LCUT, the default would be 2;3.

6.2.4. TOL (MBASIS)

(*) The set of MT products selected by SELECT can still be highly linearly dependent. Therefore, in a subsequent optimization step one diagonalizes the MT overlap matrix and retains only those eigenfunctions whose eigenvalues exceed a predefined tolerance value. This tolerance is 0.0001 by default and can be changed with the keyword TOL in the input file.

TOL 0.00001 Remove linear dependencies that fall below a tolerance of 0.00001 (see text).


Despite the favorable convergence behavior of the mixed product basis, it can still be quite large. In the calculation of the screened interaction, each matrix element, when represented in the basis of Coulomb eigenfunctions, is multiplied by \({\sqrt{v_\mu v_\nu}}\) with the Coulomb eigenvalues \({\{v_\mu\}}\). This gives an opportunity for reducing the basis-set size further by introducing a Coulomb eigenvalue threshold \({v_\mathrm{min}}\). The reduced basis set is then used for the polarization function, the dielectric function, and the screened interaction. The parameter \({v_\mathrm{min}}\) can be specified after the keyword OPTIMIZE MB in three ways: first, as a “pseudo” reciprocal cutoff radius \({\sqrt{4\pi/v_\mathrm{min}}}\) (which derives from the plane-wave Coulomb eigenvalues \({v_\mathbf{G}=4\pi/G^2}\)), second, directly as the parameter \({v_\mathrm{min}}\) by using a negative real number, and, finally, as the number of basis functions that should be retained when given as an integer. The so-defined basis functions are mathematically close to plane waves. For testing purposes, one can also enforce the usage of plane waves (or rather projections onto plane waves) with the keyword OPTIMIZE PW, in which case the Coulomb matrix is known analytically. No optimization of the basis is applied, if OPTIMIZE is omitted.

OPTIMIZE MB 4.0 Optimize the mixed product basis by removing eigenfunctions with eigenvalues below \(4\pi/4.0^2\).
OPTIMIZE MB -0.05 Optimize the mixed product basis by removing eigenfunctions with eigenvalues below 0.05.
OPTIMIZE MB 80 Retain only the eigenfunctions with the 80 largest eigenvalues.
OPTIMIZE PW 4.5 Use projected plane waves with the cutoff 4.5 Bohr-1 (for testing only, can be quite slow).

In summary, there are a number of parameters that influence the accuracy of the basis set. Whenever a new physical system is investigated, it is recommendable to converge the basis set for that system. The parameters to consider in this respect are GCUT, LCUT, SELECT, and OPTIMIZE.

Finally, we show a section MBASIS for a system with three atom types as an example

  GCUT 3.0
  LCUT 5 4 4
  SELECT 3,2;4;2 2;3 2;3
  TOL 0.0001


(*) By default, Spex constructs augmented plane waves from the mixed product basis in a similar way as in the LAPW approach. In this way, the unphysical “step” at the MT sphere boundaries is avoided. It also leads to a further reduction of the basis set without compromising accuracy. This APW construction can be avoided with the keyword NOAPW.


(*) This keyword augments the mixed product basis by the radial \({u_{lp}}\) functions, the \({l}\) and \({p}\) parameters are according to the first part of the SELECT definition. This improves the basis convergence for COHSEX calculations.


(*) Checks the accuracy of the MT product basis.


(**) Removes wavefunction coefficients belonging to radial functions that are not used to construct the MT product basis. (Currently the coefficients are actually not removed from memory.)

6.2.10. FFT (WFPROD)

When the interaction potential is represented in the mixed product basis, the coupling to the single-particle states involve projections of the form \({\langle M_{\mathbf{k}\mu} \phi_{\mathbf{q}n} | \phi_{\mathbf{k+q}n'} \rangle\,.}\) The calculation of these projections can be quite expensive. Therefore, there are a number of keywords that can be used for acceleration. Most of them are, by now, somewhat obsolete. An important keyword, though, is FFT in the section WFPROD of the input file. When used, the interstitial terms are evaluated using Fast Fourier Transformations (FFTs), i.e., by transforming into real space (where the convolutions turn into products), instead of by explicit convolutions in reciprocal space. For small systems the latter is faster, but for large systems it is recommendable to use FFTs because of a better scaling with system size. A run with FFTs can be made to yield results identical to the explicit summation. This requires an FFT reciprocal cutoff radius of \({2g_\mathrm{max}+G_\mathrm{max}}\), which can be achieved by setting FFT EXACT, but such a calculation is quite costly. It is, therefore, advisable to use smaller cutoff radii, thereby sacrificing a bit of accuracy but speeding up the computations a lot. If given without an argument, Spex will use 2/3 of the above exact cutoff. One can also specify a cutoff by a real-valued argument explicitly, good compromises between accuracy and speed are values between 6 and 8 Bohr-1.

FFT 6 Use FFTs with the cutoff 6 Bohr-1.
FFT Use the default cutoff \({\frac{2}{3}(2g_\mathrm{max}+G_\mathrm{max})}\).
FFT EXACT Use the exact cutoff \({2g_\mathrm{max}+G_\mathrm{max}}\).


(**) Since version 02.03 the interstitial part of the projections is calculated exactly, leading to favorable convergence with respect to GCUT. However, the exact evaluation takes considerably more time. If APPROXPW is set, the old and approximate way of calculating the products is used instead. (Only included for testing.)

6.2.12. LCUT and MTTHR (WFPROD)

(**) The MT part of the projections can be accelerated as well using the the keywords LCUT and MTTHR in the section WFPROD. (The keyword LCUT should not be confused with the keyword of the same name in the section MBASIS, Section 6.2.2.) The first can be used to restrict the l cutoff of the wavefunctions (only for the projections), one l cutoff for each atom type. The second is a threshold value below which wave-function coefficients are neglected. By default, the two options are unused.


(**) The keyword MINCPW modifies the plane-wave wave-function coefficients such that their absolute values become as small as possible. This reduces the error in the evaluation of the wave-function products arising from the finite G cutoff and leads to better convergence with respect to GCUT and smoother band structures, when APPROXPW is used. The coefficients can be modified, because the set of IPWs is overcomplete. If set, the eigenvectors of the overlap matrix with the n smallest eigenvalues are used for this modification. Note that it only makes sense to use this together with APPROXPW.

6.3. Wannier orbitals


For Wannier functions to be available, the Spex configure script for compilation has to be run with the option --with-wan. See Section 1.

For many features of the Spex code, Wannier functions are required: Wannier interpolation, spin-wave calculations, GT self-energy, Hubbard U parameter. Wannier functions are constructed by a linear combination of single-particle states. They can be viewed as Fourier transforms of the Bloch states. While Bloch states have a definite k vector and are delocalized over real space, the Wannier functions are localized at specific atomic sites and delocalized in k space. The mathematical definition is

(1)\[\displaystyle w_{\mathbf{R}n}(\mathbf{r}) = \frac{1}{N} \sum_\mathbf{k} e^{-i\mathbf{kR}} \sum_m U_{\mathbf{k}m,n} \phi_{\mathbf{k}m}(\mathbf{r})\]

with the lattice vector \(\mathbf{R}\) (the Wannier orbital is localized in the respective unit cell), the orbital index \(n\), the number of k vectors \(N\), and the transformation matrix \(U\). (The latter is unitary if there are as many Wannier orbitals as there are Bloch bands included in the construction. There can be more Bloch bands but not less.)

Wannier functions are defined in the section WANNIER:

  ORBITALS 1 4 (sp3)


This line is the main definition of Wannier orbitals giving the orbital characters for the first guess orbitals and the energy window. The energy window is defined by the first two arguments, the indices of the lower and upper band [summation index \(m\) in Eq. (1)]. Note that if, in the above example, the fourth state is degenerate with the fifth at some k point, the fifth state will automatically be included as well. (This is the default behavior. It can be changed to excluding the fourth band or to ignoring degeneracies altogether by setting preprocessor macros, see the source file “wannier.f”.) After the energy window follows the definition of the orbital characters of the projection functions used to calculate the first-guess Wannier functions. Although strictly not guaranteed, the orbital character is usually preserved not only after projection in the first-guess functions but even after the maximal localization procedure (see Section 6.3.2). The orbital character is defined for each atom (round brackets) or each atom type (square brackets). Empty brackets, such as () or [], are allowed. For example, (p) () (d) gives three p orbitals in the first atom, five d orbitals in the third, and none in the second. The number of bands in the energy window must not be smaller than the number of Wannier orbitals (but can be larger).

ORBITALS 1 18 (s,p,d) Define s, p, and d Wannier orbitals in the first (and perhaps only) atom.
ORBITALS 6 11 [t2g] Define t2g Wannier orbitals in the first atom type (here two atoms) from the sixth through the 11th band.
ORBITALS 1 9 () (d) Define d Wannier orbitals in the second atom.

The above example for the WANNIER section is for bulk silicon, which has two atoms in the unit cell. Four sp3 hybrid orbitals are defined for the first and none for the second. (Trailing empty brackets as in (sp3) () can be omitted.) In this case, the four bands of the energy window (1 4) are all formed by bonding sp3 orbitals. Therefore, by projection, the resulting Wannier orbitals are evenly distributed over the two atoms. This is an effect of the system’s symmetry and not explicitly caused by the particular Wannier definition. If one wants to describe the antibonding orbitals as well, one would have to define, for example, (sp3) (sp3) or just [sp3], and increase the energy window. Since the empty states of silicon are entangled, simply doubling the energy window (1 8) does not suffice. We will discuss the case of entangled bands below.

The orbital characters are defined in the same way as in Wannier90. The following orbital labels are available: s, p, d, f, px, py, pz, dxy, dyz, dz2, dxz, dx2y2, t2g, eg, fz3, fxz2, fyz2, fzxy, fxyz, fxxy, fyxy, sp, sp2, sp3, sp3d, and sp3d2. Consult the Wannier90 manual for details. In addition, there are capitalized labels (S, P, D), in which case complex spherical harmonics are used instead of real spherical harmonics. (Obviously, S is equivalent to s and included only for completeness. In practice, real and complex spherical harmonics behave identically in the present context.) Wannier functions (except for s orbitals) have an angular dependence. Sometimes it is necessary to specify a particular spatial orientation of Wannier orbitals. For this purpose, one can give Euler angles (in degrees), for example, (s,px,30,60,-90,d), through which the Wannier orbital (of px symmetry in this case) is rotated.

The so-generated Wannier function are first-guess Wannier functions. To be more precise, the expansion coefficients \(U_{\mathbf{k}n,m}\) are obtained from projecting the Bloch functions \(\phi_{\mathbf{k}n}\) onto localized functions with the chosen orbital character. Then, an orthonormalization procedure is applied to make the Wannier set orthonormal. For many purposes, the resulting set of Wannier functions is already usable. However, it is not yet maximally localized. The following keyword enforces maximal localization with Wannier90 (using the Wannier90 library), yielding the set of maximally localized Wannier functions (MLWFs).


The maximal localization procedure of Wannier90 is applied starting from the set of first-guess Wannier functions. The input file for Wannier90 (“wannier.win”) is generated automatically by Spex if it is not present yet. The user can still provide his own “wannier.win” and, in this way, fine-tune the localization procedure. However, some of the parameters (num_wan and num_bands) might still be changed by Spex if they are inconsistent with the Spex input (in ORBITALS). The output of Wannier90 is written to the file “wannier.wout”. This file should be checked for convergence of the maximal localization. For further details, consult the Wannier90 manual.


In the case of entangled bands, it is often necessary to add an extra step in the construction of MLWFs, the disentanglement, which is also an iterative optimization procedure. To this end, one defines a frozen energy window, which lies inside the outer energy window defined by ORBITALS. The frozen energy window is defined by two energies, the lower and the upper bound. The lower bound is assumed to coincide with the minimum energy of the lowest band of the outer window, i.e., the first argument of ORBITALS. If FROZEN is given without an argument, the upper bound is guessed by Spex. This guess can be overridden by giving the upper energy as an argument to FROZEN. Alternatively, one can specify the Wannier90 parameters dis_froz_min and dis_froz_max in the file “wannier.win”. Again, you should check the convergence of the disentanglement procedure in the file “wannier.wout”. For further details, consult the Wannier90 manual.

FROZEN Run disentanglement procedure. Frozen window boundaries are guessed by Spex.
FROZEN 5eV Run disentanglement procedure and use 5 eV as upper boundary.


(*) The set of Wannier functions should reflect the symmetry of the system. This can be used to define matrix representations that determine how the Wannier functions transform when a given symmetry operation is applied. These matrix representations can be used to speed up the evaluation of certain Wannier-expanded quantities, for example, in spin-wave and GT calculations. (The name IRREP is a bit misleading because the matrices are not irreducible in general.) Unfortunately, Wannier90 tends to break the symmetry of the Wannier set, so IRREP can, in most cases, not be applied to MLWFs. However, unless the ORBITALS definition itself does not break the symmetry, the first-guess Wannier functions reflect the full symmetry, and IRREP can be used.

6.3.5. RESTART

If RESTART is set, Spex writes the Wannier functions (the U matrix) to the file “spex.uwan” if it does not exist. Otherwise, the Wannier functions are read from “spex.uwan”. In the same way the file “spex.kolap” is written or read containing the overlap of wavefunctions at neighboring k points. These overlaps are needed for the maximalization procedure. (The file “spex.kolap” can be reused if the orbitals in ORBITALS are changed but not the outer energy window, i.e., the first two arguments.)


(**) This keyword makes Spex read Wannier functions generated by Fleur. (Works only with old version of Fleur. Currently unmaintained.)


The keyword INTERPOL has already been discussed in Section 5.1.17. It is not related to the Wannier construction, but we explain it here briefly for completeness. When given, Spex performs a Wannier interpolation for the reference mean-field system (output to the file “bands0”) and of the calculated GW (HF, COHSEX, et cetera) quasiparticle energies (output to the file “bands1”).