Wannier Exercise 2: Fe

This exercise has been developed by Dongwook Go. For questions, please contact via email d.go@fz-juelich.de.

In the second session of the hands-on tutorial, we will cover wannierization of metallic systems. Obtaining WFs for metals is not so different from the case of insulators. The only difference is that we need to use more number of Bloch states than the number of WFs because we need to find the optimal dimensional subspace which is disentangled from the rest of the Bloch states. We will get WFs for bcc Fe, a good old example.

Let's first define the current directory as EXERCISE_Fe as an environment variable.

EXERCISE_Fe=$PWD

and check whether it is properly defined

echo $EXERCISE_Fe

As in the previous hour, we define a data directory $(DATA_Fe)

cd Fe
DATA_Fe=$PWD

and check again the address is right

echo $DATA_Fe

Okay then let's go back to $EXERCISE_Fe and start the tutorial.

cd $EXERCISE_Fe

1. Converging the charge density and analyzing the electronic structure

1.1 Self-consistent calculation

The first thing you need to do is to converge the charge density as usual FLEUR calculations. This time, we converge the charge density by ourselves since it doesn't take much time even in most labtops.

Inpgen

Let's make a directory for generating input files from the inpgen.

mkdir inpgen

and change the directory into inpgen.

cd inpgen

An inpgen file can be found in $DATA_Fe/inpgen.Fe. Let's copy this to our current directory.

cp $DATA_Fe/inpgen.Fe .

You may check how it looks like.

cat inpgen.Fe

Note that we have included the SOC. Now, let's generate input files.

inpgen -f inpgen.Fe

Self-consistent cycles

Then let's begin the self-consistent calculation. Let's make another directory for converging the charge density.

mkdir ../conv

and copy essential files to start the self-consistent calculation.

cp inp.xml kpts.xml sym.xml ../conv

and change directory to ../conv

cd ../conv

Open inp.xml and set itmax="25" or any large number. FLEUR will stop automatically once the charge density does not change significantly over the iterations. Now we're ready run FLEUR and just wait until it gets converged...

fleur_MPI

To make sure that our converged charge density is "okay", let's check the spin and orbital magnetic moment. If you get about , it's fine.

grep -A 5 'magnetic moment' out

Band structure

We make a band structure plot, which we will compare with the Wannier band structure later. To do so, we do the band structure calculation in a separate directory.

mkdir ../band

Let's copy essential files: inp.xml, cdn.hdf, kpts.xml, sym.xml.

cp inp.xml cdn.hdf kpts.xml sym.xml ../band

and we also change the directory

cd ../band

Open band/inp.xml and set band="T". Also, change kPointList <kPointListSelection listName="path-2"/>. Then run

fleur_MPI

To make a plot of the band structure, we use a python script plots_Fe.ipynb. Open the notebook file for the plot.

2. Preparing files necessary for Wannierization

To geneate files necessary for the wannierization from FLEUR, let's copy cdn.hdf, inp.xml, sym.xml, and kpts.xml from conv directory into wann.

mkdir ../wann
cd ../wann
cp ../conv/{cdn.hdf,inp.xml,sym.xml,kpts.xml} .

2.1. kpoints

Let's make an equally-distant kpoint mesh that spans the full Brillouin zone for the wannierization.

inpgen -inp.xml -kpt wannier#gamma@grid=8,8,8 -noKsym

When you open kpts.xml, you can find kPointList whose name is "wannier"

<kPointList name="wannier" count="512" nx="8" ny="8" nz="8" type="mesh">
<kPoint weight="     0.0019531250000">   0.00/8.00    0.00/8.00    0.00/8.00</kPoint>
<kPoint weight="     0.0019531250000">   0.00/8.00    0.00/8.00    1.00/8.00</kPoint>
<kPoint weight="     0.0019531250000">   0.00/8.00    0.00/8.00    2.00/8.00</kPoint>
<kPoint weight="     0.0019531250000">   0.00/8.00    0.00/8.00    3.00/8.00</kPoint>
<kPoint weight="     0.0019531250000">   0.00/8.00    0.00/8.00    4.00/8.00</kPoint>
...

Then, open inp.xml and set kPointList <kPointListSelection listName="wannier"/> so that FLEUR uses the kPointLlist "wannier".

2.2. Defining initial projections

Now let's define initial projections. For bcc Fe, it works fine to use "usual" spd projections. Let's make a new file with its name projgen_inp

touch projgen_inp

and copy and paste

Fe 0 0 0
Fe 1 0 0
Fe 2 0 0

This should define total 18 WFs (s, px, py, pz, dxy, dyz, dzx, dz2, dx2y2 = 9, multiplied by 2 for the spin). Note that in the third column, which is supposed to indicate "m" for the angular part of the WFs, is "0". Specifying "0" chooses all the orbital having a value of "l" For example,

Fe 1 0 0

is equilvalent to

Fe 1 1 0
Fe 1 2 0 
Fe 1 3 0

Next, we modify inp.xml. Open inp.xml and find a section

<output dos="F" band="F" slice="F">
...
</output>

Add wannier="T" in output and insert a section for wannier like this:

<output dos="F" band="F" slice="F" wannier="T">
    <wannier>
    <bandSelection minSpinUp="9" maxSpinUp="44"/>
    <jobList>  projgen stopopt </jobList>
    </wannier>
    ...
````

We defined ```bandSelection``` from 9 to 44 (36 Bloch states) because there are 8 semicore states far below the valence states (we defined 3s3p states as local orbitals). Also, set ```numbands="44"``` so that FLEUR calculates all the Bloch states we need.

If you finished modifying ```inp.xml```, run


```bash
fleur_MPI

It should be done in a few second. You wil find that now proj file is generated.

cat proj

2.3. Bloch states and other miscellaneous files

Next, we need to prepare a file containing the Bloch states (eig.hdf) as well as other miscellaneous files (bkpts, WF1.win). eig.hdf is required to construct the initial projection as well as the overlap matrix . bkpts contains vectors connecting neighboring kpoints, which are necessary for evaluating finite-difference expression for . WF1.win is the input file for the WANNIER90 program.

These files can be geneated by specifying prepwan90 in inp.xml. Open imp.xml. Let's remove projgen and stopopt that we specified in the previous step, and insert prepwan90 switch:

<wannier>
<bandSelection minSpinUp="9" maxSpinUp="44"/>
<jobList>  prepwan90 </jobList>
</wannier>
````

Then run FLEUR with ```-eig hdf flag```. If you specify ```-eig hdf``` flag, the program generates ```eig.hdf``` file.



```bash
fleur_MPI -eig hdf

2.4. and

Now let's calculate (initial projection) and (overlap matrix). Also, don't forget to specify eig66="T" so that FLEUR keep using the same eig.hdf file for the consistency of the gauge.

Our inp.xml should have the following:

<output dos="F" band="F" slice="F" wannier="T" eig66="T">
    <wannier>
    <bandSelection minSpinUp="9" maxSpinUp="44"/>
    <jobList>  matrixamn matrixmmn </jobList>
    </wannier>

    ...

</output>
fleur_MPI -eig hdf

This will also generate WF1.eig that contains energy eigenvalues of the Bloch states.

3. Running WANNIER90

Now, we are finally ready to run wannier90.x. For metals, "disentanglement" is done before iteratively improving the unitary matrix. We set the maximum of the frozen window as +4 eV above the Fermi energy. Let's check the Fermi energy first.

grep fermi ../conv/out

You will see something like this:

 fermi energy and band-weighting factors:
     -->  new fermi energy            :   0.408182 htr
 fermi energy and band-weighting factors:
     -->  new fermi energy            :   0.408182 htr
 fermi energy and band-weighting factors:
     -->  new fermi energy            :   0.408182 htr
 fermi energy and band-weighting factors:
     -->  new fermi energy            :   0.408182 htr
  2   determination of fermi energy               4.59sec=   0h  0min  4sec ->   2.3%  (17calls:    0.232sec -    0.367sec)
 ```

Grep the last number and multiply by 27.2114 to get the Fermi energy in eV unit. For example, in my case, 0.408182*27.2114 = 11.1072. Open ```WF1.win``` and set ```dis_froz_max``` to be +4 eV with respect to the Fermi energy. For example, I set ``` dis_froz_max=15.1072```. Also, set ```num_iter=300```. Then run WANNIER90.


```bash
wannier90.x WF1

If you open WF1.wout the file shows WF centers and spreads after each iteration. If "Delta" becomes sufficiently smaller than the original value of the spread, then you can stop the calculation. Otherwise, we continue more iterations. This can be done by setting restart=wannierize in WF1.win. Then run again

wannier90.x WF1

I hope by now the value of the spread does not change significantly and WFs have converged. In WF1.wout file, you can find a message that the disentanglement is finished, like this:

     411      34.39115571      34.39115571      -1.134E-10      8.88    <-- DIS
     412      34.39115571      34.39115571      -1.087E-10      8.90    <-- DIS
     413      34.39115571      34.39115571      -1.041E-10      8.92    <-- DIS
     414      34.39115571      34.39115571      -9.979E-11      8.94    <-- DIS
     415      34.39115571      34.39115571      -9.562E-11      8.96    <-- DIS
     416      34.39115571      34.39115571      -9.161E-11      8.98    <-- DIS

             <<<      Delta < 1.000E-10  over  3 iterations     >>>
             <<< Disentanglement convergence criteria satisfied >>>

4. Post-processing: Band interpolation

Finally, let's interpolate the band structure. This is one of the most straightforward way to check if you WFs are okay. Open WF1.win and uncomment wannierize=plot and bands_plot=true. Also, set kpoint path,

 begin kpoint_path
 G  0.0   0.0   0.0   H -0.5   0.5   0.5
 H -0.5   0.5   0.5   N  0.0   0.5.  0.0
 N  0.0   0.5   0.0   P  0.25  0.25. 0.25
 P  0.25  0.25. 0.25  G  0.0   0.0   0.0
 G  0.0   0.0   0.0   N  0.0   0.5   0.0
 end kpoint_path

WF1.win should look like this

 length_unit=Bohr
 num_wann=          18
 num_iter=         100
 num_bands=          36

 search_shells=         200
 !optional parameters for wannierization
 !num_cg_steps=
 !trial_step=
 !fixed_step=
 !restart=wannierise

 ! optional parameters for disentangling
 !dis_win_min=
 !dis_win_max=
 !dis_froz_min=
 dis_froz_max=15.1072
 dis_num_iter=10000
 !dis_mix_ratio=
 !dis_conv_tol=
 !dis_conv_window=

 ! optional parameters for plotting
 spin=up
 restart=plot
 !wannier_plot=true
 !wannier_plot_supercell=3
 bands_plot=true
 !fermi_surface_plot=true

 ...

  begin kpoint_path
 G  0.0   0.0   0.0   H -0.5   0.5   0.5
 H -0.5   0.5   0.5   N  0.0   0.5   0.0
 N  0.0   0.5   0.0   P  0.25  0.25  0.25
 P  0.25  0.25  0.25  G  0.0   0.0   0.0
 G  0.0   0.0   0.0   N  0.0   0.5   0.0
 end kpoint_path

Then run

wannier90.x WF1

You will be able to see WF1_band.dat. We will use plots_Fe.ipynb for plotting.