Wannier Functions for Interpolation in Reciprocal Space

by Sebastian Tillack for exciting oxygen

Purpose: In this tutorial, you will learn how to compute and visualize Maximally Localized Wannier Functions (MLWFs) from a calculation using hybrid functionals. Furthermore, we show how to use MLWFs in order to achieve the electronic band-structure and density of states.

---

0. Theoretical background

---

☛ Theory of Wannier functions

Wannier functions form an alternative representation of electronic energy-bands compared to the usual Bloch functions. First, we consider a set of $J$ isolated energy bands, i.e., a set of bands that has a finite gap above and below, such as the four valence bands in diamond. Bands within such an isolated group form a closed subspace of the full space of solutions to the Schrödinger equation. Thus, the single-particle wavefunctions $\psi_{n,{\bf k}\!}({\bf r})$, $n=1,\dots,J$ form a basis of this subspace. We have the freedom to mix the bands within this subspace according to some unitary transformation $U^{\bf k}_{mn}$ obtaining a new but equally valid basis $\phi_{n,{\bf k}\!}({\bf r})$ of the subspace under consideration

(1)

\begin{align} \phi_{n,{\bf k}\!}({\bf r}) = \sum\limits_{m=1}^J U^{\bf k}_{mn}\; \psi_{m,{\bf k}\!}({\bf r}). \end{align}

The task is to find a set of unitary matrices $U^{\bf k}$ such that the Wannier functions

(2)

\begin{align} w_{n,{\bf R}\!}({\bf r}) = \frac{1}{N_{\bf k}} \sum\limits_{\bf k} {\rm e}^{-{\rm i}{\bf k}\cdot{\bf R}}\; \phi_{n,{\bf k}\!}({\bf r}) \end{align}

are as localized as possible according to some localization criterion. To this purpose, we employ the Marzari-Vanderbilt functional

(3)

\begin{align} \Omega = \sum\limits_{n = 1}^J \left(\langle{w_{n,{\bf 0}}|{\bf r}^2|w_{n,{\bf 0}}}\rangle - \langle{w_{n,{\bf 0}}|{\bf r}|w_{n,{\bf 0}}}\rangle^2\right), \end{align}

which measures the spread of the Wannier functions and which is the quantity to be minimized.

In the case of entangled energy bands, i.e., bands that are crossing each other and that are not separated from all other bands by an energy gap, first, an optimal $J$-dimensional subspace of all states within a given energy window $[E_1;E_2]$ has to be disentangled. The optimal subspace is described by a set of semi-unitary rectangular $\mathcal{J}_{\bf k}\times J$ matrices $\mathcal{U}^{\bf k}$ where $\mathcal{J}_{\bf k}$ is the number of bands, that fall within the given energy window at the point k. Once the optimal subspace

(4)

\begin{align} \widetilde{\psi}_{n,{\bf k}\!}({\bf r}) = \sum\limits_{m=1}^{\mathcal{J}_{\bf k}} \mathcal{U}^{\bf k}_{mn}\, \psi_{m,{\bf k}\!}({\bf r}) \end{align}

has been disentangled, again, we search for a set of square unitary matrices $U^{\bf k}$ to mix the states within the optimal subspace according to Eq. (1). In addition, a second (inner) energy window $[e_1;e_2]$ can be introduced within which the original states remain unchanged in the disentangled subspace, i.e., $\mathcal{U}^{\bf k}_{mn} = \delta_{mn}$ for all bands that fall inside the inner energy window.

Once we have obtained a set of maximally localized Wannier functions, they provide an optimal tight-binding basis which allows for an efficient interpolation of the wavefunction and the eigenenergies in reciprocal space. This is the so called Wannier interpolation scheme. In order to obtain the wavefunctions and the corresponding eigenenergies at any arbitrary point q in reciprocal space we set up a basis of Bloch functions

(5)

\begin{align} \phi_{n,{\bf q}\!}({\bf r}) = \sum\limits_{\bf R} {\rm e}^{{\rm i}{\bf q}\cdot{\bf R}}\, w_{n,{\bf R}\!}({\bf r}) \end{align}

in which we will expand the interpolated wavefunctions

(6)

\begin{align} \psi_{n,{\bf q}\!}({\bf r}) = \sum\limits_{m=1}^J V^{\bf q}_{mn}\, \phi_{n,{\bf q}\!}({\bf r}). \end{align}

The expansion coefficients $V^{\bf q}$ and the eigenenergies corresponding to the wavefunctions $\psi_{n,{\bf q}}$ are obtained by the diagonalization of the small $J \times J$ Hamiltonian $H^{\bf q}_{mn} = \langle \phi_{m,{\bf q}}|\hat{\bf H}|\phi_{n,\bf q}\rangle$. In practice, the interpolation points q can be points along high-symmetry paths through the Brillouin zone as we need it in the calculation of a band structure or a very dense grid of points as it is needed in the calculation of the density of states.

---

1. Prerequisites

The system we will consider throughout this tutorial is diamond.

The interpolation scheme based on the Wannier functions that is presented in this tutorial is of general validity for any system. However, it is of fundamental importance for all those systems in which, due to the complexity of the structure and/or of the used methodology, the calculation of electronic properties on a fine grid in reciprocal space is very expensive or even prohibitive. In particular, this applies to the case of calculations performed using GW and/or hybrid exchange-correlation functionals. Indeed, for this latter case the Wannier interpolation is the only way of plotting the electronic band-structure implemented in exciting. As an example, the application to hybrid functionals is presented in this tutorial.

We assume that you have successfully completed the tutorial on Hybrid-functional calculations and that the PBE, PBE0, and HSE folders already exist as quoted in that tutorial. All the calculations which are described in the present (Wannier) tutorial can be performed in the corresponding directories already considered in Hybrid-functional calculations.

In the next we assume that the calculation performed in Hybrid-functional calculations are stored (according to the general notation used in all tutorials) in the starting folder /home/tutorials/diamond-hybrids. This should be adapted if you used a different notation.

Then, to start the tutorial move to the reference directory

$ cd /home/tutorials/diamond-hybrids

First, we will focus on the PBE0 exchange-correlation functional. The same approach can be used for HSE, but it will not be treated explicitly in this tutorial.

---

2. Calculation of Wannier functions, electronic band-structure, and DOS for PBE0

First, move to the PBE0 directory

$ cd PBE0

Then, modify the input file input.xml present in the current directory. The changed file looks like this

<input>
 
   <title>Diamond PBE0 Wannier</title>
 
   <structure speciespath="$EXCITINGROOT/species" >
      <crystal scale="6.7425">
         <basevect> 0.0     0.5     0.5 </basevect>
         <basevect> 0.5     0.0     0.5 </basevect>
         <basevect> 0.5     0.5     0.0 </basevect>
      </crystal>
      <species speciesfile="C.xml">
         <atom coord="0.00 0.00 0.00" />
         <atom coord="0.25 0.25 0.25" />
      </species>
   </structure>
 
   <groundstate 
      do="skip"
      ngridk="4 4 4"
      rgkmax="5.0"
      nempty="20"
      xctype="HYB_PBE0">
      <Hybrid excoeff="0.25" />
   </groundstate>
 
   <properties>
 
      <!-- CALCULATION OF THE WANNIER FUNCTIONS -->
      <wannier input="hybrid">
         <group
            method="opfmax"
            fst="1"
            lst="4"/>
         <group
            method="disFull"
            outerwindow="0.0 2.5"
            innerwindow="0.0 1.5"
            nwf="13"/>
      </wannier>
 
   </properties>
 
</input>

So what have we changed?

• First, note that we have set the attribute do="skip" in the groundstate element in order not to redo the relatively expensive ground-state calculation.
• Second, we have added the element properties which in turn contains the subelement wannier which will be discussed in details in the next subsection.

Before continuing, do not forget to update the speciespath.

$ SETUP-excitingroot.sh

---

i) Calculation of the Wannier functions

The calculation of the Wannier functions is the crucial part in the Wannier interpolation scheme. The construction of the Wannier functions can be very time demanding in complex systems.

In exciting this task is performed inside the element properties by the block defined by the new element wannier. Here, is an XML snippet of this block:

      <!-- CALCULATION OF THE WANNIER FUNCTIONS -->
      <wannier input="hybrid">
         <group
            method="opfmax"
            fst="1"
            lst="4"/>
         <group
            method="disFull"
            outerwindow="0.0 2.5"
            innerwindow="0.0 1.5"
            nwf="13"/>
      </wannier>

---

☛ The element "wannier": Attributes and subelements

• With input="hybrid" we define from which kind of calculation the used eigenvectors and eigenenergies come from. The default is "gs", which stands for ground state calculations with local/semi-local functionals.
• Within the wannier element we have declared two groups of bands.
  • The first group describes the isolated set of the four valence bands in diamond specified by the first state fst="1" and the last state lst="4". lst is always the last occupied state, while the optimal first valence state fst needs to be determined.
  • In order to compute the band structure of the unoccupied states as well we further have to construct Wannier functions describing these states. This is accounted for by the second group. Since the conduction bands do not form an isolated group of bands we do not specify a band range but energy windows instead. The outer window is given by the attribute outerwindow="0.0 2.5". By default, the values are relative to the Fermi level, i.e., in this case we would like to disentangle the subspace from a 2.5 Hartree window above the Fermi energy which corresponds to the lowest conduction bands. Optionally, one can define a second energy window within which the states are described exactly. We do this by setting the attribute innerwindow="0.0 1.5". The inner window has to be fully contained within the outer window. With nwf="13" we specify the dimension of the subspace or, equivalently, the number of Wannier functions to compute from the given windows. In general, setting up the energy windows and the number of Wannier functions for entangled bands can be a delicate task. For a detailed description and analysis, we refer to paper [1] describing the present implementation.
• The method being employed for both groups differ. For the first group of isolated bands, first, an initial guess to the Wannier functions is obtained using so called optimized projection functions (OPF). In a second step, maximally localized Wannier functions are computed by minimizing the spread functional.
• For the second group of entangled bands four steps in total are carried out. First, again, we create an initial guess using the OPF method. In the second step, starting from the initial guess, the optimal disentangled subspace at each k-point is computed. In the third step, we again employ OPFs (but now restricted to the disentangled subspace) to generate an initial guess for the final minimization of the spread functional (step four). For this last step, two different approaches are implemented. The first approach is the standard introduced by Souza, Marzari, and Vanderbilt. Here, the disentangled subspace is kept fixed during the minimization of the spread functional and we minimize only w.r.t. square unitary matrices. This behaviour is triggered by method, which is set here to the default value "disFull". This choice allows to vary the disentangled subspace during the minimization of the spread, i.e., we minimize w.r.t. the full rectangular semi-unitary matrices.

---

So let's run the calculation with the command

$ time exciting_smp

Once the calculation is finished, the newly created file WANNIER_INFO.OUT informs us about the run and the result. Information about the obtained Wannier functions is printed at the end of the file. To look at the result in our example type the following command.

$ tail -n 33 WANNIER_INFO.OUT 

********************************************************************************
* Wannier functions                                                            *
********************************************************************************

   #             localization center (lattice)           Omega      Omega_I      Omega_D     Omega_OD
================================================================================
   1        0.125000     0.125000    -0.375000        2.338824     2.061300     0.000000     0.277524
   2        0.125000     0.625000     0.125000        2.338816     2.061292    -0.000000     0.277524
   3        0.125000     0.125000     0.125000        2.338837     2.061311     0.000000     0.277526
   4        0.625000     0.125000     0.125000        2.338822     2.061297     0.000000     0.277525
--------------------------------------------------------------------------------
                                           total:     9.355299     8.245200    -0.000000     1.110099
--------------------------------------------------------------------------------
   5        0.120930     0.128501     0.127383        2.941415     2.337524     0.001128     0.602763
   6       -0.289975     0.424800    -0.204073        3.918751     3.199509     0.022979     0.696264
   7       -0.036571    -0.466187    -0.023396        4.116418     3.522015     0.022880     0.571522
   8        0.123271     0.625753     0.126672        2.699352     2.174199     0.000953     0.524199
   9        0.623196     0.127398     0.128491        2.940993     2.337185     0.001136     0.602671
  10       -0.483137    -0.017098    -0.016904        4.167311     3.349846     0.056859     0.760605
  11        0.124324     0.126654     0.625750        2.699603     2.174401     0.000956     0.524247
  12        0.269960     0.276855    -0.270852        4.219640     3.504343     0.025629     0.689668
  13       -0.284412     0.284418     0.284226        4.132022     3.432774     0.037851     0.661398
  14        0.069336    -0.203440     0.424872        3.919580     3.200111     0.023047     0.696423
  15        0.283943    -0.283960    -0.283717        3.913075     3.310386     0.016546     0.586143
  16       -0.473647    -0.023365     0.533685        4.116625     3.522252     0.022834     0.571539
  17       -0.276292    -0.270570     0.276717        4.217628     3.502779     0.025547     0.689302
--------------------------------------------------------------------------------
                                           total:    48.002413    39.567324     0.258345     8.176744
================================================================================
                                           total:    57.357712    47.812524     0.258345     9.286843
                                         average:     3.373983     2.812501     0.015197     0.546285

We see the Wannier functions corresponding to the two groups that we have specified. The first four Wannier functions correspond to the first group and the valence bands. The remaining 13 Wannier functions correspond to the second group. The second to fourth column display where in the crystal the corresponding Wannier functions are localized. The last four columns give us information about the individual spreads of the Wannier functions.

---

ii) Visualization of Wannier functions

Now we would like to take a closer look on how the actual functions we have computed in the previous step look like. From the result printed above we already see that the four Wannier functions corresponding to the four valence bands are perfectly centered in the middle of the four tetrahedral bonds of a carbon atom. In order to visualize that we have to make a few changes to the input file.

• Modify the line containing the element wannier as follows:

      <wannier input="hybrid" do="fromfile">

• Add the following XML snippet inside the element properties

      <!-- VISUALIZATION OF THE WANNIER FUNCTIONS -->
      <wannierplot fst="1" lst="4" cell="1 1 1">
         <plot3d>
            <box grid="60 60 60">
               <origin coord="0.0 0.0 0.0"/>
               <point  coord="1.0 0.0 0.0"/>
               <point  coord="0.0 1.0 0.0"/>
               <point  coord="0.0 0.0 1.0"/>
            </box>
         </plot3d>
      </wannierplot>

By adding do="fromfile" to the wannier element we tell exciting to read the already computed Wannier functions from the file WANNIER_TRANSFORM.OUT. Within the new element wannierplot we have do define a range of Wannier functions we want to visualize labeled by the first column in the above output. Here, we are going the visualize the four Wannier functions corresponding to the valence bands. With cell="1 1 1" we define in which unit cell we want the functions to be localized. This corresponds to the real-space lattice vector R labeling the Wannier function w_n,R(r). The real-space grid on which the functions are evaluated is specified by grid="60 60 60". The three point elements define a box (given by cartesian vectors) in which the functions are computed. The box is always centered around the origin which is given relative to the localization center of the corresponding Wannier function, i.e., setting the origin to zero already selects the region in real space where the Wannier function is localized.

Again we start the calculation with the command

$ time exciting_smp

After the calculation is finished, you can visualize the result in xcrysden by executing the script

$ PLOT-wannier-function.sh

You are asked which function you want to visualize (type in a number from 1 to 4) and whether you want to see the real part, the imaginary part, or the square modulus. Note, that in this case the corresponding Wannier functions are fully real valued so type 'r' for real part. Once xcrysden has opened you have to chose an isovalue to plot. A value around 0.025 is a good choice here. For the third function the result will look like this.

---

ii) Electronic band-structure using the Wannier interpolation scheme

Now that we have successfully constructed a set of maximally localized Wannier functions we can make use of it to calculate an accurate band structure. In order to perform the calculation using Wannier interpolation, do the following changes in the input file.

• Check if the line containing the element wannier is written as follows:

      <wannier input="hybrid" do="fromfile">

• Inside the element properties, substitute the XML snippet appearing under "VISUALIZATION OF THE WANNIER FUNCTIONS" with the following one:

      <!-- CALCULATION OF THE ELECTRONIC BAND-STRUCTURE -->
      <bandstructure wannier="true">
         <plot1d>
            <path steps="200">
               <point coord="1.0     0.0     0.0" label="Gamma"/>
               <point coord="0.625   0.375   0.0" label="K"/>
               <point coord="0.5     0.5     0.0" label="X"/>
               <point coord="0.0     0.0     0.0" label="Gamma"/>
               <point coord="0.5     0.0     0.0" label="L"/> 
            </path>
         </plot1d>
      </bandstructure>

Notice that in comparison with the standard case, the element bandstructure must have the explicit attribute wannier = "true".

Run the calculation with

$ time exciting_smp

For looking explicitly at the electronic band-structure calculated using the PBE0 functional and the Wannier (WA) interpolation scheme, you can execute the script PLOT-band-structure.py (for a detailed description of the script arguments see The python script "PLOT-band-structure.py") as follows:

$ PLOT-band-structure.py -a WA  -e -30 30

The generated plot PLOT.png shows the PBE0 band structure obtained using the Wannier interpolation scheme.

---

ii) Density of states using the Wannier interpolation scheme

After the calculation of the electronic band-structure, we show how to obtain the density of states (DOS) with and without the use of the Wannier interpolation scheme. In order to perform the calculation using Wannier interpolation, do the following changes in the input file.

• Check if the line containing the element wannier is written as follows:

      <wannier input="hybrid" do="fromfile">

• Inside the element properties, substitute the XML snippet appearing under "CALCULATION OF THE ELECTRONIC BAND-STRUCTURE" with the following one:

      <!-- CALCULATION OF THE DENSITY OF STATES -->
      <dos
         wannier="true"
         ngridkint="40 40 40"
         winddos="-1 1"
         nwdos="1000"
         inttype="tetra"/>

The calculation of the DOS involves an integration over the Brillouin zone. In practice, this is approximated by a summation over a finite grid of reciprocal space points sampling the Brillouin zone.

In the snippet above, the attribute winddos="-1 1" defines the energy window for which the DOS is calculated and nwdos="1000" gives the number of energy subdivisions within this window. With inttype="tetra" we specify the integration method (tetrahedron integration in this case).

Notice that in order to use the Wannier interpolation scheme, the element dos must contain the explicit attribute wannier = "true". Thus, the initial mesh specified by ngridk = "4 4 4" inside groundstate is "expanded" by Wannier interpolation to the mesh specified by ngridkint = "40 40 40".

Run the calculation as usual with

$ time exciting_smp

For looking explicitly at the DOS calculated using the PBE0 functional and the Wannier (WA) interpolation scheme, you can execute the script PLOT-dos.py (for a detailed description of the script arguments see The python script "PLOT-dos.py") as follows:

$ PLOT-dos.py -a WA  -t "ngridk = '4 4 4' + WA('40 40 40')"

The generated plot PLOT.png will look like the following

For the sake of comparison, we show in the next the calculation of the DOS using the PBE0 functional and the original mesh (KS) specified by ngridk = "4 4 4" without Wannier interpolation. To this purpose, modify the attributes of the element dos as follows:

      <!-- CALCULATION OF THE DENSITY OF STATES --> 
      <dos
         winddos="-1 1"
         nwdos="1000"
         inttype="tetra"/>

Run the calculation with

$ time exciting_smp

Now, we can compare the DOS calculated using the original mesh ("4 4 4") and the previus one obtained by Wannier (WA) interpolation of the original mesh to a "40 40 40" grid in reciprocal space. In order to visualize the results, execute

$ PLOT-dos.py -a WA KS  -ll "ngridk = '4 4 4' + WA('40 40 40')" "ngridk = '4 4 4'"  -rp

The resulting PLOT.png file should look like the following.

---

Literature

1. S. Tillack, A. Gulans, and C. Draxl, Phys. Rev. B 101, 235102 (2020)

• Some further reading:
  • N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997)
  • I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001)
  • J. I. Mustafa, S. Coh, M. L. Cohen, and S. G. Louie, Phys. Rev. B 92, 165134 (2015)