PbTiO3
First we generate the input files using the input generator. We use the following inpgen.inp file for the input generator:
PbTiO3
&input film =f /
7.33438000 0.00000000 0.00000000
0.00000000 7.33438000 0.00000000
0.00000000 0.00000000 7.83312000
1.0
1.0 1.0 1.0
5
8 0.5 0.5 0.1040250
82 0.0 0.0 0.0
22 0.5 0.5 0.534678
8 0.5 0.0 0.612905
8 0.0 0.5 0.612905
Run inpgen2 -f inpgen.inp in order to obtain the inp.xml input file. Run fleur_MPI in order to converge the charge density. Check that the charge density is converged.
After converging the charge density we start with the construction of MLWFs. For this purpose, generate a new subdirectory Wann in your working directory. Copy inp.xml, cdn.hdf and sym.xml into the Wann subdirectory. Change into the Wann subdirectory. When you inspect the inp.xml file you realize that the input generator has defined the following local orbitals to treat the semicore states:
<species name="Pb-1" element="Pb" atomicNumber="82" flipSpinPhi=".00000000" flipSpinTheta=".00000000" flipSpinScale="F">
...
<lo type="SCLO" l="2" n="5" eDeriv="0"/>
for Pb and
<species name="Ti-1" element="Ti" atomicNumber="22" flipSpinPhi=".00000000" flipSpinTheta=".00000000" flipSpinScale="F">
...
<lo type="SCLO" l="0" n="3" eDeriv="0"/>
<lo type="SCLO" l="1" n="3" eDeriv="0"/>
for Ti. Thus, there are 10 (l=2) + 2 (l=0) + 6 (l=1)=18 local orbitals. Since this is a nonmagnetic calculation without SOC, there are 9 spin-up bands associated with these 18 local orbitals. Fleur considers only the spin-up bands, because the spin-down bands are identical to the spin-up bands in this case. Therefore, one might guess at first that the first 9 bands describe semicore states corresponding to these local orbitals. However, the three oxygen S-bands are lower in energy than the Pd-d semicore states.
When you plot the band structure you find that the occupied bands below the band gap may be split into well-separated groups of bands. These groups of bands are the following:
Bands 1-4: Titanium semicore states.
Bands 5-7: Oxygen s states.
Bands 8-12: Pb d semicore states.
Bands 13: Pb s state.
Bands 14-22: Oxygen p states.
We expect that the displacement of the MLWF centers of the oxygen p MLWFs contributes significantly to the ferroelectric polarization. Therefore, we compute the MLWFs for the oxygen p states first.
Therefore, we add the section "wannier" to inp.xml, where we specify minSpinUp="14" and maxSpinUp="22" in the bandSelection:
<output dos="F" band="F" vacdos="F" slice="F" mcd="F" wannier="T">
<wannier>
<bandSelection minSpinUp="14" maxSpinUp="22"/>
<jobList> projgen prepwan90 stopopt </jobList>
</wannier>
</output>
Note that we added wannier="T" to the line output and we added a jobList with the keywords projgen, prepwan90, and stopopt. In order to make sure that fleur computes the first 22 bands we specify numbands="22" in the line cutoffs in inp.xml:
<calculationSetup>
<cutoffs Kmax="3.50000000" Gmax="10.50000000" GmaxXC="10.50000000" numbands="22"/>
Create the file projgen_inp with the following contents:
O 1 0 0
Create a file kpts.xml with a uniform equidistant k-mesh. Running fleur_MPI should produce a proj file, which looks like this:
9 9
1 1 1 0
0.000000 0.000000 0.000000 0.000000 1.00
1 1 2 0
0.000000 0.000000 0.000000 0.000000 1.00
1 1 3 0
0.000000 0.000000 0.000000 0.000000 1.00
4 1 1 0
0.000000 0.000000 0.000000 0.000000 1.00
4 1 2 0
0.000000 0.000000 0.000000 0.000000 1.00
4 1 3 0
0.000000 0.000000 0.000000 0.000000 1.00
5 1 1 0
0.000000 0.000000 0.000000 0.000000 1.00
5 1 2 0
0.000000 0.000000 0.000000 0.000000 1.00
5 1 3 0
0.000000 0.000000 0.000000 0.000000 1.00
Modify the jobList in inp.xml as follows:
<wannier>
<bandSelection minSpinUp="14" maxSpinUp="22"/>
<jobList> matrixmmn matrixamn </jobList>
</wannier>
Run wannier90.x WF1 and converge the MLWFs. The final state should look similar to this:
Final State
WF centre and spread 1 ( 3.667190, 3.667190, 1.004554 ) 4.76227375
WF centre and spread 2 ( 3.667190, 3.667190, 1.012155 ) 7.50473611
WF centre and spread 3 ( 3.667190, 3.667190, 1.012303 ) 7.50631920
WF centre and spread 4 ( 3.667190, 0.000000, -2.874997 ) 7.88991260
WF centre and spread 5 ( 3.667190, -0.000000, -2.651064 ) 8.97821707
WF centre and spread 6 ( 3.667190, 0.000000, -3.085763 ) 4.66023727
WF centre and spread 7 ( -0.000000, 3.667190, -2.874901 ) 7.89444578
WF centre and spread 8 ( 0.000000, 3.667190, -2.651588 ) 8.96285758
WF centre and spread 9 ( 0.000000, 3.667190, -3.085915 ) 4.67065915
Sum of centres and spreads ( 22.003140, 22.003140,-14.195216 ) 62.82965850
Now, we may compute the contribution of the displacement of the MLWF centers
to the ferroelectric polarization (using the key dipole2) and the
ferroelectric polarization computed from nominal charges associated with
ionic sites (using the key dipole3).
Change the jobList in inp.xml as follows:
```xml
and set up the file [dipole](wannierdipole.md) like this:
```text
0
Additionally, set up the file IONS like this:
3
Ti 4.0
Pb 2.0
O -2.0
Note also the explanations given for the key dipole2 and for the key dipole3. Running fleur.x generates the file dipole_out, which lists the contributions due to the displacement of the Oxygen-p MLWFs with respect to their atomic host sites: The polarization due to the displacement of Oxygen-P orbitals is (0.0,0.0,-42.181) uC/cm2. Additionally, the file polarization_out has been generated, which lists the contributions computed by assuming nominal ionic charges: The polarization due to the displacements of ions is: (0.0,0.0,-55.41) uC/cm2.
Now, we have to calculate the contributions of the remaining orbitals. It is convenient to create new directories for these calculations:
mkdir TiP; cp * TiP mkdir OS; cp * OS mkdir PbD; cp * PbD mkdir PbS; cp * PbS
We proceed in the same way as for the MLWFs of the oxygen p states.
At the end we sum up all contributions:
Total polarization:
-55.41 (Ions)
-42.18 (Oxygen-P)
12.41 (Lead-S)
13.47 (Lead-D)
-16.58 (Oxygen-S)
0.96 (Ti-P)
---------------------------
-87.32 uC/cm2