Wannier Exercise 2: Fe

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 $N_\mathrm{wann}$ 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 $2.2\mu_\mathrm{B}$, 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 $A_{mn}(\mathbf{k})$ as well as the overlap matrix $M_{mn}^{(\mathbf{k},\mathbf{b})}$. bkpts contains vectors connecting neighboring kpoints, which are necessary for evaluating finite-difference expression for $M_{mn}^{(\mathbf{k},\mathbf{b})}$. 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. $A_{mn}(\mathbf{k})$ and $M_{mn}^{(\mathbf{k},\mathbf{b})}$

Now let's calculate $A_{mn}(\mathbf{k})$ (initial projection) and $M_{mn}^{(\mathbf{k},\mathbf{b})}$ (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=wannierise 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 restart=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.