**by****Sebastian Tillack****for****exciting****nitrogen**

** Purpose:** In this tutorial, you will learn how to compute and visualize

`(`

**Maximally Localized Wannier Functions**`) 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.`

**MLWFs**### 0. Theoretical background

### 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

`calculation.`

**PBE0**`$ 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`=

`in the`

**"skip"**`element in order not to redo the relatively expensive groundstate calculation. Second, we have added a new element`

**groundstate**`inside the`

**wannier**`element. With`

**properties**`=`

**input**`we define from which kind of calculation the used eigenvectors and eigenenergies come from. Within the`

**"hybrid"**`element we have declared two groups of bands.`

**wannier**The first ` group` describes the isolated set of the four valence bands in diamond specified by the first state

`=`

**fst**`and the last state`

**"1"**`=`

**lst**`. The corresponding Wannier functions are going to be calculated using the`

**"4"**`=`

**method**`which means that we will compute maximally localized Wannier functions starting from an initial guess obtained using the optimized projection functions (`

**"opfmax"**`) method. This is the recommended option for isolated groups of bands.`

**OPF**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`=

`. Therefore we have to specify an energy window from which we want to disentangle the subspace. This window is given by the attribute`

**"disentangle"**`=`

**outerwindow**`. By default, the values are relative to the Fermi level,`

**"0.0 2.5"***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**`. The inner window has to be fully contained within the outer window. With`

**"0.0 1.5"**`=`

**nwf**`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`

**"13"**`method.`

**OPF**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`=

`to the`

**"fromfile"**`element we tell the program to read the already computed Wannier functions from the file`

**wannier**`. Within the new element`

**WANNIER_TRANSFORM.OUT**`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`

**wannierplot**`=`

**cell**`we define in which unit cell we want the functions to be localized. This corresponds to the real-space lattice vector`

**"1 1 1"****labeling the Wannier function $w_{n,{\rm R}\!}({\rm r})$. The real-space grid on which the functions are evaluated is specified by**

`R``=`

**grid**`. The three`

**"60 60 60"**`elements define a box (given by cartesian vectors) in which the functions are computed. The box is always centered around the`

**point**`which is given relative to the localization center of the corresponding Wannier function,`

**origin***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

`by executing the script`

**xcrysden**`$ 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

**for real part. Once**

`'r'``has opened you have to chose an isovalue to plot. A value around`

**xcrysden****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

`element.`

**wannierplot**... <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`=

`in the`

**"true"**`element.`

**bandstructure**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

`directory (`

**PBE**`) and the`

**PBE**`directory (`

**PBE0**`). The generated plot`

**PBE0**`compares the Wannier (`

**PLOT.png**`) interpolated band structure to the one obtained from a simple Fourier (`

**W**`) interpolation.`

**F**For a converged calculation with ` ngridk`=

`and`

**"8 8 8"**`=`

**nempty**`the result will be the following.`

**"100"**### 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**`-points. Therefore, go back into the`

**k**`directory and open the input file. Comment out again the`

**PBE0**`element and add a new element`

**bandstructure**`with the following attributes`

**dos**`$ 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`=

`defines the energy window for which the`

**"-1 1"**`is calculated and`

**DOS**`=`

**nwdos**`gives the number of energy subdivisions within this window. Run the calculation as usual with`

**"1000"**`$ 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

`element`

**dos**... <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

`which will look like this:`

**PLOT.png**In contrast to the ` DOS` obtained from the coarse grid, the Wannier-interpolated

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

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

### 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

`that results from interpolating the eigenenergies from a coarse grid on a dense grid of`

**DOS**

**80$\times$80$\times$80**`-points with the exact`

**k**`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`

**DOS**`xc-functional on a coarse`

**GGA_PBE**`-grid (`

**k**`) 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`

**4$\times$4$\times$4**`).`

**GGA PBE**The following plot compares the ` DOS` that was obtained by performing a groundstate calculation on a coarse grid (

`) and then interpolating the eigenenergies on a dense grid (`

**4$\times$4$\times$4**`) using the Wannier interpolation scheme (in red) with the`

**80$\times$80$\times$80**`obtained by directly performing the groundstate calculation on the dense grid (in blue).`

**DOS**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

`does not change a lot. The impact on the conduction band region, however, is slightly bigger.`

**DOS**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.

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