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.

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

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