Wannier Functions for Interpolation in Reciprocal Space

by Sebastian Tillack for exciting nitrogen

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 very accurate 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. 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

We assume that you have successfully completed the tutorial on Hybrid-Functional calculations. If you have followed all the instructions change to the directory that contains the PBE0 calculation.

$ cd /home/tutorials/diamond-pbe0/PBE0

---

2. Calculation of Wannier Functions

First, we need do edit the input file input.xml 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"
      nempty="50"
      xctype="HYB_PBE0">
      <Hybrid 
         exchangetype="HF" 
         excoeff="0.25" />
   </groundstate>
 
   <properties>
      <wannier input="hybrid">
         <group
            method="opfmax"
            fst="1"
            lst="4"
            />
         <group
            method="disentangle"
            outerwindow="0.0 2.5"
            innerwindow="0.0 1.5"
            nwf="13"
            />
      </wannier>
 
      <!--bandstructure>
         <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-->
 
   </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 groundstate calculation. Second, we have added a new element wannier inside the properties element. With input="hybrid" we define from which kind of calculation the used eigenvectors and eigenenergies come from. 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". The corresponding Wannier functions are going to be calculated using the method="opfmax" which means that we will compute maximally localized Wannier functions starting from an initial guess obtained using the optimized projection functions (OPF) method. This is the recommended option for isolated groups of bands.

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 first have to disentangle an optimal subspace by choosing method="disentangle". Therefore we have to specify an energy window from which we want to disentangle the subspace. This 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. Once the optimal subspace has been disentangled, maximally localized Wannier functions describing the states within this subspace are constructed using the OPF method.

Lastly, for the moment, we have commented out the bandstructure element.

So let's run the calculation with the commands

$ SETUP-excitingroot.sh 
$ time excitingser

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. In this case the result looks like this

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

   #             localization center (lattice)           Omega      Omega_I      Omega_D     Omega_OD
================================================================================
   1       -0.375000     0.125000     0.125000        2.336738     2.058068     0.000000     0.278671
   2        0.125000    -0.375000     0.125000        2.336739     2.058069     0.000000     0.278670
   3        0.125000     0.125000    -0.375000        2.336739     2.058069     0.000000     0.278670
   4        0.125000     0.125000     0.125000        2.336738     2.058068     0.000000     0.278670
--------------------------------------------------------------------------------
   5        0.328563    -0.251972    -0.254088        4.375840     3.218150     0.012214     1.145476
   6        0.126959     0.125980     0.621102        2.709870     2.158849     0.000103     0.550918
   7       -0.037920    -0.183492    -0.176403        5.314959     3.450804     0.037860     1.826295
   8       -0.272731     0.269726     0.268635        4.320478     3.293540     0.016872     1.010066
   9        0.620562     0.126433     0.126612        2.726057     2.169360     0.000043     0.556654
  10        0.127146     0.620604     0.126417        2.768403     2.188094     0.000110     0.580200
  11       -0.275585     0.271960    -0.268321        4.266192     3.270256     0.020738     0.975199
  12        0.127264     0.125657     0.126312        2.802881     2.204766     0.000149     0.597965
  13        0.335265    -0.251735     0.168452        4.371655     3.209141     0.011206     1.151308
  14       -0.040184     0.400780    -0.175191        5.220403     3.405928     0.036545     1.777930
  15       -0.040406    -0.183562     0.409204        5.093304     3.353998     0.043103     1.696203
  16        0.329522     0.176198    -0.253867        4.381938     3.219398     0.012284     1.150256
  17       -0.273801    -0.266006     0.269343        4.313935     3.290692     0.017796     1.005447
================================================================================
                                           total:    62.012872    46.665250     0.209025    15.138598
                                         average:     3.647816     2.745015     0.012296     0.890506

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.

---

3. 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.

   ...
   <properties>
      <wannier input="hybrid" do="fromfile">
         ...
      </wannier>
 
      <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>
      ...
   </properties>
   ...

By adding do="fromfile" to the wannier element we tell the program 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,{\rm R}\!}({\rm 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 excitingser

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 fourth function the result will look like this.

---

4. Band structure using Wannier interpolation

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 do so we reactivate the bandstructure element in the input file and comment out the wannierplot element.

   ...
   <properties>
      <wannier input="hybrid" do="fromfile">
         ...
      </wannier>
 
      <!--wannierplot fst="1" lst="4" cell="1 1 1">
         ...
      </wannierplot-->
 
      <bandstructure wannier="true">
         ...
      </bandstructure>
   </properties>
   ...

The only change we have to make is to set the attribute wannier="true" in the bandstructure element.

Run the calculation with

$ time excitingser

In order to display the result change to the parent directory and run the visualization script

$ cd ../
$ PLOT-compare-bands.py -25 25 wannier

 Comparison for WANNIER calculations

################################################

 Enter the names of the 2 directories to compare

------------------------------------------------

 Directory 1 ==> PBE
 Directory 2 ==> PBE0

################################################

Enter the name of the PBE directory (PBE) and the PBE0 directory (PBE0). The generated plot PLOT.png compares the Wannier (W) interpolated band structure to the one obtained from a simple Fourier (F) interpolation.

For a converged calculation with ngridk="8 8 8" and nempty="100" the result will be the following.

[converged]

---

5. Density of states using Wannier interpolation

The calculation of the density of states (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 order to reduce the error of this approximation a very dense sampling is needed which makes the Wannier interpolation scheme a desirable tool for this purpose.

First, however, let us compute the DOS from the solutions on the original grid of 4$\times$4$\times$4 k-points. Therefore, go back into the PBE0 directory and open the input file. Comment out again the bandstructure element and add a new element dos with the following attributes

$ cd PBE0/

   ...
   <properties>
      <wannier input="hybrid" do="fromfile">
         ...
      </wannier>
 
      <!--wannierplot fst="1" lst="4" cell="1 1 1">
         ...
      </wannierplot-->
 
      <!--bandstructure wannier="true">
         ...
      </bandstructure-->
 
      <dos
         winddos="-1 1"
         nwdos="1000"
         ngrdos="500"
         newint="true"
          />
   </properties>
   ...

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. Run the calculation as usual with

$ time excitingser

Now, we would like to use the Wannier interpolation scheme to interpolate to a much denser grid. In order to do so make the following changes to the dos element

      ...
      <dos
         wannier="true"
         ngridkint="80 80 80"
         winddos="-1 1"
         nwdos="1000"
         ngrdos="500"
         newint="true"
          />
      ...

This specifies that we want to use Wannier interpolation and that we are going to interpolate the energies on a grid of 80$\times$80$\times$80 k-points. Again, run

$ time excitingser

To visualize the result run the script

$ PLOT-compare-dos.py wannier

This will generate a plot called PLOT.png which will look like this:

In contrast to the DOS obtained from the coarse grid, the Wannier-interpolated DOS shows the correct square-root-like behaviour on the valence band edges which one would expect for nearly parabolic bands in a bulk material such as in diamond.

The DOS obtained for the settings of the converged hybrid functional calculation can be seen below (note different energy scale).

[converged]

---

6. Validation of the Wannier interpolation quality

In this last section we would like do discuss the quality of the Wannier interpolation. We do this at the example of the DOS. What we have to do is to compare the DOS that results from interpolating the eigenenergies from a coarse grid on a dense grid of 80$\times$80$\times$80 k-points with the exact DOS that we would obtain by performing the original calculation on the same dense grid. However, it is practically unfeasible to perform the expensive hybrid functional calculation on such a dense grid. What we can do instead is to compare the DOS obtained from Wannier interpolation based on a calculation using the GGA_PBE xc-functional on a coarse k-grid (4$\times$4$\times$4) with the exact result which can be obtained by performing the calculation on the dense grid (which is doable with (semi) local xc-functionals such as GGA PBE).

The following plot compares the DOS that was obtained by performing a groundstate calculation on a coarse grid (4$\times$4$\times$4) and then interpolating the eigenenergies on a dense grid (80$\times$80$\times$80) using the Wannier interpolation scheme (in red) with the DOS obtained by directly performing the groundstate calculation on the dense grid (in blue).

Again, note the smaller gap obtained using GGA PBE compared to the hybrid functional calculation. Besides that the change in the xc-functional only has a minor effect on the valence bands and the valence DOS does not change a lot. The impact on the conduction band region, however, is slightly bigger.

In order to further increase the quality of the Wannier interpolation we have to increase the k-point density of the supporting coarse grid. The same plot as above but with a coarse grid of 8$\times$8$\times$8 points can be found below.

[converged]

This plot shows that in this case a supporting grid of 8$\times$8$\times$8 k-points is enough for the Wannier interpolation to perfectly reproduce the exact converged result of the DOS. In turn, it is very likely that the Wannier interpolated DOS for the hybrid functional calculation on the same grid (converged plot in section 5) can be considered to be exact as well.

• wannierbugfix.tar.gz

---

Literature

• 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)