.. This is a raw version of the Example documentation. This file
   was compiled into a pdf file and html pages.
   Please consult :file:`../FPLO../DOC/pyfplo/pyfplo.doc`
   or open :file:`../FPLO../DOC/pyfplo/html/index.html`

.. |refcalculate3dTIInvariants| replace:: `~pyfplo.slabify.Slabify.calculate3dTIInvariants`
.. |refkscale| replace:: `~pyfplo.slabify.Slabify.kscale`


3D topological insulator
========================

.. Contents::
   :local:
   :backlinks: top

.. note:: You need to use the newer **xfbp**\ /\ **xfplo** version,
   which comes with `pyfplo` in order for the ``cmd`` scripts to work
   properly.

Here we discuss the calculation of the topological indices for 3d
topological insulators. The algorithm is based on Wannier centers as
in :ref:`Sec2dTI`. The only difference is that we need to calculate
four :math:`Z_2` invariants for 4 planes in the 3d BZ. In principle
there are 6 different planes forming the faces of a parallelepiped
whose corners are the eight distinct time reversal invariant momenta
(TRIM) of any chosen smallest primitive reciprocal cell. However, if
the electronic spectrum is gapped above a certain band there are only
4 independent :math:`Z_2` invariants, which are usually grouped into
four topological indices
:math:`\nu_0;(\nu_1\nu_2\nu_3)`. :math:`\nu_0=1` indicates a strong
topological insulator while any :math:`\nu_{1,2,3}=1` indicate
conditions on the number of surface state pairs for certain surface
orientations. If :math:`\nu_0=0` but some :math:`\nu_{1,2,3}\neq 0` it
is a weak topological insulator. If all indices are zero it is a
trivial insulator. The mapping of :math:`Z_2` invariants to indices is
as follows: if any two parallel planes of the parallelepiped have
differing invariants :math:`\nu_0=1`. We take the two planes spanned
by the first two reciprocal lattice vectors (:math:`g_{1,2}`) called
:math:`z_0` if the plane goes through the origin and :math:`z_1` if it
goes through :math:`g_{3}`. The three weak indices are identical to
the :math:`Z_2` invariants of the three planes, which do not go
through the origin, which we call :math:`x_1` (spanned by
:math:`g_{2,3}` through :math:`g_1`), :math:`y_1` (spanned by
:math:`g_{3,1}` through :math:`g_2`) and :math:`z_1` (spanned by
:math:`g_{1,2}` through :math:`g_3`). The weak indices of course
depend on the chosen TRIM points.

The algorithm selects the TRIM points and performs four :math:`Z_2`
calculations.  It tries to determine the invariants via the automatic
procedure explained in :ref:`Sec2dTI`. The results including
reliability and electronic gap estimates are printed to the output and
convenience files :file:`Z2_3dTI_homo...cmd` for use in
**xfbp** are created.

This algorithm works for centro symmetric and non-centro symmetric
compounds. A version of this algorithm is linked directly into
**fplo**. The difference is that when used from **fplo** a all-orbital
Wannier basis is created instead of the reduces Wannier models, which
usually are the basis of `pyfplo.slabify`. This leads to differently
locking Wannier center curves, since there are more orbitals (semi
core and such). If surface state calculations are planned it is anyway
necessary to create a Wannier model first. The resulting
:math:`Z_2`-calculations are faster than in **fplo**.

For illustration we chose a system, which is centro symmetric, since
in this case we know the invariants from the parity algorithm, which
allows comparison. This example also shows some fo the complications.
We will use both the **fplo** version and the `pyfplo` version of the
Wannier center algorithm to understand the pecularities.

The tutorial files are in
:file:`FPLO.../DOC/pyfplo/Examples/slabify/3dTI/Bi2Se3`.

To start with change into this directory and have a look into the
topological insulator submenu in **fedit**. You will see that "Force
wannier center" is switched on, to force the use of this algorithm
despite the fact that this compound has an inversion center.

Now, have a look at :file:`makewandef.py`. It illustrates how to
create a :file:`=.wandef` quickly. There is a particularity in the
system, which is that the Bi6p and Se4p orbitals are separated from
lower and higher bands through two gaps. This makes creating a Wannier
model especially easy. We use the option *lbands* and *ubands* to
define an energy window which excludes all other bands.
Please execute ::

  python makewandef.py

You should have a =.wandef file now. Next run **fplo**. E.g.::

  fplo.... > out

to produce a bandstructure, :file:`+wancoeff` and the output from the
**fplo** Wannier center calculation. This will take some time, since the
**fplo** version of the WC algorithm is not the fastest, due to the
larger number of bands needed. Output related to the TI calculation
can be found via::

  grep TI: out

First, the invariants from parities are printed and then the Wannier
center algorithm is executed. Various files are created.
**Note, that the indices from parities can differ from the Wannier
center results; first, because the WC algorithm is tricky or second,
if the electronic gap for a certain homo is zero, in which case both
algorithms makes no sense.**

You will see the following output from the Wannier centers::

  TI: Invariants:
  TI:   homo  invariants  E(k=0)                estim.gaps [eV] for
  TI:                                  z0        z1        x1        y1
  TI:     114  ? 0;(000)  -1.09567   0.05585  0.01917  0.01917  0.01917 
  TI:     116  ? 1;(111)  -0.79835   0.01306  0.03873  0.03873  0.03873 
  TI:     118 *  1;(000)  -0.22789   0.44699  0.43513  0.43513  0.43513 
  TI:     120  ? 1;(111)   0.28647   0.15850  0.14845  0.14845  0.14845 
  TI:     122  ? 0;(111)   1.42635   0.02117  0.17378  0.17378  0.17378 
  TI:     124  ? 0;(000)   1.44752   0.21087  0.29692  0.29692  0.29692 
  TI: 
  TI:   homo  invariants  E(k=0)                  Z2 (reliability)
  TI:                                 z0        z1        x1        y1
  TI:     114  ? 0;(000)  -1.09567    0 ( 77%)  0 ( 77%)  0 ( 77%)  0 ( 77%)
  TI:     116  ? 1;(111)  -0.79835    0 ( 50%)  1 ( 73%)  1 ( 73%)  1 ( 73%)
  TI:     118 *  1;(000)  -0.22789    1 (100%)  0 (100%)  0 (100%)  0 (100%)
  TI:     120  ? 1;(111)   0.28647    0 ( 92%)  1 ( 65%)  1 ( 65%)  1 ( 65%)
  TI:     122  ? 0;(111)   1.42635    1 ( 69%)  1 (100%)  1 (100%)  1 (100%)
  TI:     124  ? 0;(000)   1.44752    0 ( 80%)  0 (100%)  0 (100%)  0 (100%)


The first table shows homos, indices the energy at the Gamma point and
estimates of the electronic gap (taken from the values on the finite
grid used for the WC algorithm), while the second shows the
:math:`Z_2` invariants and their estimated reliability for the four
planes of the parallelepiped mentioned above. If the reliability is
less than 100% a question mark is printed after the homo. The homo,
which is likely the highest occupied band is marked by a * for
orientation. This is of course only reliable if the system has a true
gap.

In our case we see that there are rather small estimated gaps for lots
of bands. Fortunately above the Fermi energy (homo 118) the gap is
sizable. Not, surprisingly the reliabilities are all very high
(actually 100%). We also find lots of question marks.
Let's have a look at the :ref:`FigBi2Se3FPLObands`.


.. _FigBi2Se3FPLObands:
.. figure:: pics/3dTI/Bi2Se3/fplo/wbands.png
   :figwidth: 17cm
   :width: 17cm
   :align: center

   fplo and Wannier model bands (red/green)

We see the two gaps, separating the p-bands as discussed above and
that nearly all bands are separated by gaps (although small). This
means that topological invariants have a meaning. Note, that it is not
needed that the gap is visible in the DOS. All one needs is that the
bands of homo and homo+1 never cross (warped gap). 

   
Since this table is not ment for straight forward consumption we have
to inspect the Wannier centers by hand. Run::

  xfbp Z2_3dTI_homo118.cmd

and so forth. You will get several pictures. Let's start with
:ref:`FigBi2Se3FPLOWChomo118`. As already pointed out in the table his
band has a nice electronic gap above itself, happens to be the last
occupied band below the Fermi level and has high reliability of the
automatically determined :math:`Z_2` invariants (reference line
winding number).

.. _FigBi2Se3FPLOWChomo118:
.. figure:: pics/3dTI/Bi2Se3/fplo/Z2_3dTI_homo118.png
   :figwidth: 10cm
   :width: 10cm
   :align: center

   Wannier centers and reference line for homo 118

It can be clearly seen that the :math:`x_1`, :math:`y_1` and
:math:`z_1` planes have zero :math:`Z_2` invariant, while plane
:math:`z_0` is non-trivial. To really verify this use the right mouse
click in **xfbp** close to the higest Wannier center curve in the
:math:`z_0` panel close to :math:`\theta=1` and :math:`k_y=0`. It will
show that only one center (called Set118 in **xfbp**) comes down from
it's degenerate value of +-1. A clean straight reference line can be
drawn for :math:`\theta=0.9` (not shown), which only crosses this
center and hence crosses an odd number of Wannier centers, which
results in :math:`Z_2=1`. This and the fact that we can visually
connect the Wannier center curves in a reasonable smooth way convinces
us that the topological indices are 1;(000). (Remember :math:`z_0\neq
z_1 \to \nu_0=1` and :math:`\nu_{1}=Z_2(x_1), \ldots`.)  In fact due
to symmetry the tree planes :math:`x_1`, :math:`y_1` and :math:`z_1`
give the same result in this compound.

Next we have a look at all the other homos, just for learning purposes.

.. _FigBi2Se3FPLOWChomo114:
.. figure:: pics/3dTI/Bi2Se3/fplo/Z2_3dTI_homo114.png
   :figwidth: 10cm
   :width: 10cm
   :align: center

   Wannier centers and reference line for homo 114

The :ref:`FigBi2Se3FPLOWChomo114` are quite messy. However, you should
convince yourselfs that :math:`z_1` has a region (:math:`\theta\approx0.5`)
without any centers crossing :math:`\to` trivial, while for the
:math:`z_0`-plane an even number of curves cross any reference line. Hence we
get 0;(000).

.. _FigBi2Se3FPLOWChomo116:
.. figure:: pics/3dTI/Bi2Se3/fplo/Z2_3dTI_homo116.png
   :figwidth: 10cm
   :width: 10cm
   :align: center

   Wannier centers and reference line for homo 116

The :ref:`FigBi2Se3FPLOWChomo116` are also messy and reveal another
possible complication. First, note that both 114 and 116 yield curves
with fast varying (steep slopes) curves. These steep curves generate
only very few points on their steep section. Depending on the density
of the ky mesh there might even be zero points on such sections, which
would make it look as if certain curves end for some ky value to
reemerge for a later ky-value. Any algorithm including our own
judgment can fail in such cases. All we can do is to take a denser
ky-mesh. In this case the automatic reference line algorithm is
correct as can be visually verified. There is an even/odd number of
curves crossed by a suitably chosen reference curve for the
:math:`z_0`/:math:`z_1` planes
(:math:`x_1,y_1=z_1`). Hence we get 1;(111).  Note,
however, that the electronic gaps above homos 114 and 116 are rather
small, they could well be zero. If this homo were important we would
try to use much finer meshes to make absolutely sure that the gap is
finite.

Homo 122 also has a small electronic gap. If we assume that
it is finite we can proceed analysing the invariants. The
:ref:`FigBi2Se3FPLOWChomo120` and the :ref:`FigBi2Se3FPLOWChomo122`
both show that the indices printed in the output table are wrong (as
can be verified in this case due to the parity algorithm).

.. _FigBi2Se3FPLOWChomo120:
.. figure:: pics/3dTI/Bi2Se3/fplo/Z2_3dTI_homo120.png
   :figwidth: 10cm
   :width: 10cm
   :align: center

   Wannier centers and reference line for homo 120

For homo 120 the :math:`z_0` plane is clearly trivial since we can draw a
non crossing reference line at :math:`\theta=\pm1` (the automatic blue
reference line is just that). For the other three planes the algorithm
made a mistake around ky=0.33, where a single red dot appears. This
mistake is due to a steep section at one curve. One can find hand
drawn reference curves, which cross an even number of times and hence
these planes are trival. The result should be 0;(000). Please also
note that the reliability for homo 120 is high for the :math:`z_0` plane but
only 65% for the other three, which fits to the findings explained
above. This is of course only an indicator.

.. _FigBi2Se3FPLOWChomo122:
.. figure:: pics/3dTI/Bi2Se3/fplo/Z2_3dTI_homo122.png
   :figwidth: 10cm
   :width: 10cm
   :align: center

   Wannier centers and reference line for homo 122
   
For homo 122 the situation is even trickier. According to the
reliabilities the :math:`z_0` panel might have a problem. Indeed, the
:math:`z_1` plane (and :math:`x_1`, :math:`y_1`) is clearly
topological, but the :math:`z_0` plane was determined as non trivial
although it is not. Again due to small electronic gaps in the
:math:`z_0` plane we get very steep curves close to ky=0. In fact this
curve is so steep that we see only one point in the steep section. In
**xfbp** use the right mouse click at the dot in the upper left
corner. It will show that there are two curves. Given the rest of the
dots it is clear that one dot is connected with the declining curve
and the other must fall very steeply to connect to a curve close to
:math:`\theta=0.28`. Hence the reference line actually crosses two
curves and not one. Consequently, :math:`z_0\neq z_1 \to` 1;(111) in
contrast with the result printed in the output table.

Homo 124 is again a simple case and will not be discussed here.

In summary we saw that the **fplo** based version of the procedure is
rather slow, which is bad if a high grid density is needed to obtain
reliable results. Let's continue with the `pyfplo` version then. We
already set up the Wannier basis production in our first **fplo**
run. Now, we can proceed by re-running::

  fplo...

Note that we do not use an output file (we don't want to overright the
TI results). This will produce among others the file :file:`+hamdata`,
which we will use for the `pyfplo` based version of the TI procedure.
If the fplo run is over we confirm the Wannier basis by checking::

  xfbp wband.cmd

which should show green bands (finite cut-off tight-binding Wannier
bands) on top of red bands (exact Wannier bands) on top of black bands
(**fplo** bands). Please, change into :file:`slabify/Z2` and execute::

  python -u Z2.py | tee out

The flag ``-u`` un-buffers the output, such that the progress is updated
in real time. The file's content is
  
.. literalinclude:: ../../DOC/pyfplo/Examples/slabify/3dTI/Bi2Se3/slabify/Z2/Z2.py
   :linenos:

In lines 34-47 we calculate the band structure of the Wannier model
along the same high symmetry points as were used in the fplo
calculation. With the help of::

  xfbp slabifyres/+band_sf


and a right mouse click we find out that the Fermi level is above homo
18 (just a coincidence that 118-18=100). This is used in line 56 to
define a set of *homos* and *efhomo*. The grid is chosen finer for
both the "integration" and the ky direction. See |refcalculate3dTIInvariants|.

The output table::



  Invariants:
    homo  invariants  E(k=0)                estim.gaps [eV] for
                                   z0        z1        x1        y1
       14  ? 0;(000)  -1.09645   0.02338  0.01563  0.01563  0.01563 
       16  ? 1;(111)  -0.79820   0.00531  0.02098  0.02098  0.02098 <-small
       18 *  1;(000)  -0.22063   0.42925  0.43611  0.43611  0.43611 
       20    0;(000)   0.27854   0.14357  0.06937  0.06937  0.06937 
       22    1;(111)   1.42038   0.00911  0.17176  0.17176  0.17176 <-small
       24    0;(000)   1.42948   0.20549  0.24304  0.24304  0.24304 
  
    homo  invariants  E(k=0)                  Z2 (reliability)
                                  z0        z1        x1        y1
       14  ? 0;(000)  -1.09645    0 (100%)  0 ( 93%)  0 ( 93%)  0 ( 93%)
       16  ? 1;(111)  -0.79820    0 (100%)  1 ( 86%)  1 ( 93%)  1 ( 93%)
       18 *  1;(000)  -0.22063    1 (100%)  0 (100%)  0 (100%)  0 (100%)
       20    0;(000)   0.27854    0 (100%)  0 (100%)  0 (100%)  0 (100%)
       22    1;(111)   1.42038    0 (100%)  1 (100%)  1 (100%)  1 (100%)
       24    0;(000)   1.42948    0 (100%)  0 (100%)  0 (100%)  0 (100%)

  
  
looks a bit different from the previous one. First of all
the gap values are different since there are always tiny shifts in a
reduced Wannier model compared to a full band structure and second we
used a finer grid. Furthermore, the reliabilities increased
considerably, again mostly due to the finer grid. Finally, all
topological indices are correct.

We can now load all the :file:`Z2_3dTI_homo...cmd` files into **xfbp**
and analyse the validity by hand.  We will look at three Wannier
center graphs. The :ref:`FigBi2Se3FPLOWChomo16` show a  quick change
in the :math:`z_0`-panel around ky=0.16, which is most likely due to the small
electronic gap. Furthermore, the steep sections in the :math:`z_1`-panel are
now represented by more points, which makes the automatic reference
line algorithm more reliable.


.. _FigBi2Se3FPLOWChomo16:
.. figure:: pics/3dTI/Bi2Se3/slabify/Z2_3dTI_homo16.png
   :figwidth: 10cm
   :width: 10cm
   :align: center

   Wannier centers and reference line for homo 16

The :ref:`FigBi2Se3FPLOWChomo20` has the correct reference line due to
the increased mesh size.

.. _FigBi2Se3FPLOWChomo20:
.. figure:: pics/3dTI/Bi2Se3/slabify/Z2_3dTI_homo20.png
   :figwidth: 10cm
   :width: 10cm
   :align: center

   Wannier centers and reference line for homo 20


The :ref:`FigBi2Se3FPLOWChomo22` now clearly shows that two WC curves
meet close to the upper left corner of the :math:`z_0`-panel, which makes this
plane trivial. Consequently, the invariants are correctly determined
as 1;(111). 
           
.. _FigBi2Se3FPLOWChomo22:
.. figure:: pics/3dTI/Bi2Se3/slabify/Z2_3dTI_homo22.png
   :figwidth: 10cm
   :width: 10cm
   :align: center

   Wannier centers and reference line for homo 22


In summary, the use of a reduce Wannier function model leads to a much
faster TI procedure, which allows to use finer grids, which in turn
increase the accuracy of the results. Additionally, one could
calculate surface spectra using the same model.
