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 bandstructure and density of states.
Table of Contents

0. Theoretical background
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 exchangecorrelation functionals. Indeed, for this latter case the Wannier interpolation is the only way of plotting the electronic bandstructure 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 Hybridfunctional 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 Hybridfunctional calculations.
In the next we assume that the calculation performed in Hybridfunctional calculations are stored (according to the general notation used in all tutorials) in the starting folder /home/tutorials/diamondhybrids. This should be adapted if you used a different notation.
Then, to start the tutorial move to the reference directory
$ cd /home/tutorials/diamondhybrids
First, we will focus on the PBE0 exchangecorrelation 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 bandstructure, 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 groundstate 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.
$ SETUPexcitingroot.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>
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 realspace lattice vector R labeling the Wannier function w_{n,R}(r). The realspace 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
$ PLOTwannierfunction.sh
ii) Electronic bandstructure 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 BANDSTRUCTURE > <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 bandstructure calculated using the PBE0 functional and the Wannier (WA) interpolation scheme, you can execute the script PLOTbandstructure.py (for a detailed description of the script arguments see The python script "PLOTbandstructure.py") as follows:
$ PLOTbandstructure.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 bandstructure, 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 BANDSTRUCTURE" 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 PLOTdos.py (for a detailed description of the script arguments see The python script "PLOTdos.py") as follows:
$ PLOTdos.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
$ PLOTdos.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
 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)