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



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>


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 wn,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.
wannier_function.png
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.

bs-PBE0+Wannier.png
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

DOS-k4-WA-k40.png

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.

DOS-k4-vs-k4-WA-k40.png

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)

 

 

http://exciting.wikidot.com/oxygen-wannier-functions