Wannier Exercise 1: PbTiO3
This exercise has been developed by Dongwook Go based on the previous documentation written by Frank Freimuth. For questions, please contact via email d.go@fz-juelich.de.
In this tutorial, we will learn how to obtain Wannier functions (WFs) for PbTiO3, one of the famous insulating ferroelectric materials. In general, the procedure consists of three steps:
- Converging the charge density and analyzing the electronic structure (from FLEUR)
- Preparing the files necessary for WANNIER90 (from FLEUR)
- Obtaining WFs (from WANNIER90)
Once you obtain WFs, we can proceed to post-processings. In this example, we will do the band interpolation.
Let's first define the current directory as EXERCISE_PbTiO3 as an environment variable.
EXERCISE_PbTiO3=$PWD
and check whether it is properly defined
echo $EXERCISE_PbTiO3
By the way, in this exercise, we provide you an already-converged charge density for PbTiO3. The charge density and other essential files are found in Wannier/data/PbTiO3, which we define as $(DATA_PbTiO3)
cd PbTiO3
DATA_PbTiO3=$PWD
and check again the address is right
echo $DATA_PbTiO3
Okay then let's go back to $EXERCISE_PbTiO3 and start the tutorial.
cd $EXERCISE_PbTiO3
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. You can find an inpgen file for PbTiO3 in DATA_PbTiO3. Let's have a look
cat $DATA_PbTiO3/inpgen.PbTiO3
You may copy and paste the above to generate inputs for the FLEUR calculation and converge the charge density. For me, it took about 25 iterations until the charge density fully converges for the default setting from the inpgen. Depending on where you run the calculation, it might take a while. So, in this tutorial, we provide you an an already-converged charge density assuming that you succeeded to converge.
1.2. Electronic structure analysis
Before proceeding further, it is important to analyze the band structure and their orbital characters because from this information we dicide which initial projections to use for the WFs. This is especially important when you wannierize a new class of materials for the first time. Once you know the recipe, you can apply same/similar projections to materials in the same class.
Again, to save the time, I have calculated the band structure and local density of state (LDOS), which can be found in $DATA_PbTiO3. We will use a python script plots_PbTiO3.ipynb. Please open this and check the band structure and local DOS.
If you have done that, you will see that the result can be summarized as the following:
- A band located at ~ -7 eV comes from Pb s-orbital
- A group of bands from -5 eV to 0 originates in O p-orbitals
- A group of bands from 1 eV to 8 eV mainly comes from Ti d-orbitals and O p-orbitals.
What we are going to do is construct WFs that describe the second group by using p-like WFs located near the O atoms. There are in total 9 bands, and we will use px, py, pz orbitals for three O atoms (3x3=9).
2. Preparing files necessary for Wannierization
To geneate files necessary for the wannierization from FLEUR, we need cdn.hdf, inp.xml, sym.xml, and kpts.xml. These can be found at $DATA_PbTiO3. Let's make another directory for the wannierization and copy these files.
mkdir wann
cd wann
cp $DATA_PbTiO3/{cdn.hdf,inp.xml,sym.xml,kpts.xml} .
2.1. kpoints
To prepare for Wannier functions, we use an equally-distant kpoint mesh that spans the full Brillouin zone. This does not necessarily need to be same as the kpts used in the self-consistency calculation. In most of the cases, 8x8x8 k-mesh is sufficient as long as the final WFs do not have significant overlap over 8 unit cell distance. An additional kPointList can be generated by using inpgen:
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 one we defined above.
2.2. Defining initial projections
In order to define initial projections, we need to provide proj file. For example, proj file looks like this
4 4
1 -3 1 0
0.000000 0.000000 0.000000 0.000000 1.00
1 -3 2 0
0.000000 0.000000 0.000000 0.000000 1.00
1 -3 3 0
0.000000 0.000000 0.000000 0.000000 1.00
1 -3 4 0
0.000000 0.000000 0.000000 0.000000 1.00
...
and you can make your own proj file if you like. Each column defines how the function behaves in their angular and radial parts. Detailed explanation of the file format can be found in here.
But you can simply use FLEUR to define, e.g. px, py, pz orbitals localized at O atoms, by providing projgen_inp. To define px, py, pz orbitals for O atoms, make a new file and name it as projgen_inp in the directory wann:
touch projgen_inp
Then open projgen_inp and copy and paste the following:
O 1 1 0
O 1 2 0
O 1 3 0
The first column specifies the atom, the second and third column are for l and m values (FLEUR follows the convention from WANNIER90. For the definition of l and m, please check Table 3.1 and 3.2 of the official documentation for WANNIER90). The last column specifies the radial dependence, where 0 chooses FLEUR radial function corresponding to the chosen angular momentum.
Next, we modify inp.xml. Open inp.xml and find a section
<output dos="F" band="F" slice="F">
...
</output>
and 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="14" maxSpinUp="22"/>
<jobList> projgen stopopt </jobList>
</wannier>
...
</output>
````
Here, ```bandSelection``` specifies the Bloch states that we are going to use to construct WFs. ```minSpinUp``` and ```maxSpinUp``` define indices of the Bloch states, which are labeled in the increasing order in energy.
```<bandSelection minSpinUp="14" maxSpinUp="22"/>``` means that we use from 14th Bloch state to 22nd Bloch state (in total 9 Bloch states). This is because
1. We defined local orbitals, 5d states for Pb and 3s3p states for Ti, 5+4=9 (Please check ``` <lo type="SCLO" ... > ``` in ```inp.xml```.). So, we have 9 semicore states far below the Fermi energy.
2. There are 3 O s-states (one s-state from each O atoms).
3. We do not include Pb s-state, which is located at ~ -9 eV.
In ```jobList```, as the name indicates we specify jobs related to wannierization and post-processings. Here, we use ```projgen``` in order to generate ```proj```. ```stopopt``` ensures that FLEUR finishes the job as soon as ```projgen``` is over and does not enter into self-consistency cycles. Then run FLEUR
```bash
fleur_MPI
It should be finished in a 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 necessary to construct the initial projection as well as the overlap matrix . bkpts contains vectors connecting neighboring kpoints, which are necessary for evaluating finite-difference expressions in . WF1.win is the input file for the WANNIER90 program.
These files can be geneated by specifying prepwan90 in inp.xml. Open inp.xml. Remove projgen and stopopt that we specified in the previous step, and insert prepwan90 switch:
<wannier>
<bandSelection minSpinUp="14" maxSpinUp="22"/>
<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
To run WANNIER90 we need to provide (initial projection) and (overlap of the Bloch states between neighboring kpoints). These can be generated by specifying matrixamn and matrixmmn in jobList.
From this point of the calculation including post-processings, we should keep using the same eigenstates to make the gauge convention cosistent over the evaluation of various physical quantities. All physical quantities need to be evaluated within the same guage choice. eig66 switch ensures that FLEUR does not calculate the eigenvector anymore and keep using eig.hdf which we generated before.
Thus, open inp.xml and modify the following:
<output dos="F" band="F" slice="F" wannier="T" eig66="T">
<wannier>
<bandSelection minSpinUp="14" maxSpinUp="22"/>
<jobList> matrixamn matrixmmn </jobList>
</wannier>
...
</output>
Note that we set eig66="T". Run FLEUR with -eig hdf switch so that it uses the existing eig.hdf file.
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. Before running the program, let's quickly check WF1.win, which specifies options and parameters for obtaining WFs. But for now, we do not need to change anything, except for num_iter. Let's keep it to 500 and run
wannier90.x WF1
If you open WF1.wout it 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, you can continue more iterations from the point where you stopped. 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.
4. Post-processing: Band interpolation
Post-processings in general consist of three separate steps. Let's assume that we want to calculate an observable . What we have to do is:
- Evaluation of in the Bloch basis.
- Obtaining in the Wannier basis by Fourier transformation together with .
- Interpolation in k-space by Fourier-transforming .
FLEUR provides several routines for the steps 1 and 2. WANNIER90 also has an user-friendly interface for this. Especially, the band interpolation can be done very easily within WANNIER90. To obtain the band structrue from the WFs, you need to open WF1.win and change a few keys. First, comment out restart=wannierize. Instead, uncomment restart=plot and bands_plot=true. We also need to provide high-symmetry lines at the end of WF1.win.
For example, WF1.win should look like this
length_unit=Bohr
num_wann= 9
num_iter= 100
num_bands= 9
search_shells= 200
!optional parameters for wannierization
!num_cg_steps=
!trial_step=
!fixed_step=
!restart=wannierise
! 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 X 0.5 0.0 0.0
X 0.5 0.0 0.0 M 0.5 0.5 0.0
M 0.5 0.5 0.0 G 0.0 0.0 0.0
G 0.0 0.0 0.0 Z 0.0 0.0 0.5
Z 0.0 0.0 0.5 R 0.5 0.0 0.5
R 0.5 0.0 0.5 A 0.5 0.5 0.5
A 0.5 0.5 0.5 Z 0.0 0.0 0.5
end kpoint_path
Then run
wannier90.x WF1
The band structure is written in WF1_band.dat.
head -20 WF1_band.dat
The first column is for the kpoint and the second column is for the energy eigenvalues. Let's check if the band structure obtained from the WFs agrees with the DFT band structure. Check section 3 of plots_PbTiO3.ipynb.
Does your band structure calculated from the Wannier basis agree with the FLEUR band structure? Then you're done. Congratulations!