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