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

.. |MoTe2| replace:: MoTe\ :sub:`2`

.. |BandPlot| replace:: `~pyfplo.common.BandPlot`

.. |findWeylPoints| replace:: `~pyfplo.slabify.Slabify.findWeylPoints`

.. |kscale| replace:: `~pyfplo.slabify.Slabify.kscale`

.. |BoxMesh| replace:: `~pyfplo.slabify.BoxMesh`

.. |BoxMesh.origin| replace:: `~pyfplo.slabify.BoxMesh.origin`

.. |WeylPoint| replace:: `~pyfplo.slabify.WeylPoint`

.. |WeylPoints| replace:: `WeylPoints <pyfplo.slabify.WeylPoint>`

.. [Wang15] Zhijun Wang et.a. arxiv:1511.07440 (2015)
               
.. [Fukui05] [Takahiro Fukui et.al. J. Phys. Soc. Jpn. 74, 1674 (2005)]
                         
.. _examplswsm:

Weyl semi metals
================

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


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

In this example we discuss methods for Weyl semi metals. We
demonstrate how to find Weyl points, how to prove that they are indeed
Weyl points and how to calculate surface spectra to analyse potential
topological surface states (Fermi arcs).

We start with |MoTe2|. This example is based on [Wang15]_.
The tutorial files are in
:file:`FPLO.../DOC/pyfplo/Examples/slabify/Weylpoints/MoTe2`.
They contain a converged full relativistic calculation.



Figure :ref:`FigMoTe2str` shows the bulk unit cell and band structure.
There is clearly an isolated band complex between -6 and +5 eV, which
makes the creation of a Wannier function model especially easy.




.. _FigMoTe2str:
.. figure:: pics/Weylpoints/MoTe2/strband.png
   :figwidth: 12cm
   :width: 12cm
   :align: center

   The bulk unit cell and band structure of |MoTe2|.


With a
right mouse click close to the highest and lowest band in **xfbp**
(when :file:`+band` is loaded) we determined that this band complex
has 88 bands. From the available valence orbitals we guess that these
must be the 4d-bands of the 4 Mo atoms and the 5p-bands of the 8 Te
atoms, which makes a total of 4*10+8*6=88 orbitals/bands. We could
also have used band weights to determine this. Fortunately, none of
these 88 bands has a section somewhere in the BZ where the sum of the
weights of the considered orbitals is exactly zero. In other words
there is no band section which needs other orbitals to get a non-zero
WF projector. This allows to build a very simple WF model. The example
directory contains the file :file:`makewandef.py` which has the
following content


.. literalinclude:: ../../DOC/pyfplo/Examples/slabify/Weylpoints/MoTe2/makewandef.py
   :linenos:

Line 23 defines the projectors for the 4d-orbitals of Mo sites
1--4. Note, that we use the python expression `range(1,4+1)` to obtain
the list `[1,2,3,4]`. Similarily, line 26 defines the projectors for
the 5p-orbitals of Te sites 5-12. The energy window encompasses the
isolated band complex and we use  a zero upper energy window tail to
keep the higher bands out of the desired Hilbert space.

Please execute this script::
  
  python makewandef.py

after which you should have the new file :file:`=.wandef`. We run
**fplo** now to obtain the necessary data for the WF calculation in
:file:`+wancoeff`::

   fplo.... > out

Convince yourself that the file :file:`+wancoeff` has been created and
re run **fplo** to calculate the Wannier functions::

   fplo.... > outwf

Now, the file :file:`+hamdata` was created, which contains the WF
model. In order to check the quality of the WF fit execute::

  xfbp wband.xpy

You will see something like :ref:`FigMoTe2WFfit`. The black curves are
the DFT results, the red curves are the exact WF transformed bands and
the green curves the bands resulting from the real space WF model. 

.. _FigMoTe2WFfit:
.. figure:: pics/Weylpoints/MoTe2/wbands.png
   :figwidth: 12cm
   :width: 12cm
   :align: center

   The WF fit.


If you use a zoom function you will notice that the green curves
deviate from the exact transform around the Fermi level by let's say
20meV. We could do better by using larger WF cutoffs and more SCF
k-points, but for our purpose this is good enough.

To check the validity of the result we change into
:file:`slabify/3d`::

      cd slabify/3d

and execute the two following commands::

  python bands.py
  xfbp bands3d.xpy

with the result :ref:`FigMoTe2WFbands3d`.


.. _FigMoTe2WFbands3d:
.. figure:: pics/Weylpoints/MoTe2/bands3d.png
   :figwidth: 12cm
   :width: 12cm
   :align: center

   The 3d bulk band structure from the WF model.


Have a look at the script :file:`bands.py`
   
.. literalinclude:: ../../DOC/pyfplo/Examples/slabify/Weylpoints/MoTe2/slabify/3d/bands.py
   :linenos:


In line 27-32 we read the high symmetry points from the **fplo** input
file :file:`=.in` which we then assign to the |BandPlot| object in
line 42. Lien 43 prepares the path through the BZ and line 44 does the
actual calculation. Do not forget to have a look at the **xfbp**
script :file:`bands3d.xpy` e.g. via the Edit->Script/Transformations
editor of **xfbp**.

Next, in **xfbp** use a right mouse click at the highest band between
X and S below the Fermi level to read off the band number of this
band. You will get a popup menu with a number of bands (which might
not be sorted). The highest number in this list will be 56. This will
be the band for which we check the Weyl points.


Now, that we are confident that the model looks like the DFT result we
will try to find Weyl points. For this we change into
:file:`../wpsearch`::

      cd ../wpsearch

and have a look at the script :file:`wpsearch.py`

.. literalinclude:: ../../DOC/pyfplo/Examples/slabify/Weylpoints/MoTe2/slabify/wpsearch/wpsearch.py
   :linenos:

In line 25-27 the default bulk 3d unit cell is setup. In line 30 we
remember the number of the homo i.e. the band which forms the lower
two legs of a potential Weyl point crossing. Line 32 contains the
number of points of the initial search grid for the algorithm. Please
consult |findWeylPoints| for a description of the search algorithm.
Line 33 defines the smallest bisection cell size.

In lines 35-47 we set up the primitive reciprocal unit cell vectors
(`x`, `y`, `z`) in units of |kscale| (usually
:math:`\frac{2\pi}{a_0}`).  and their length (`dx`, `dy` and `dz`) as
well as a subdivision in each direction which results in the most
isotropic mesh cells. These values are used to define a |BoxMesh| in
lines 54-58.  Note line 51, in which we assign -1/2 of the body
diagonal of a mesh cell to `k0`. This is used to shift the origin of
the |BoxMesh| such that we do not miss WPs which lay on a mesh cell
boundary. Note, that |BoxMesh.origin| shifts the mesh and is not an
absolute origin.  Finally, line 61 executes the Weyl point
search. First, for each cell of the |BoxMesh| the Berry curvature is
integrated over the cell surface, which is done with compact formulas
as described for 2d in [Fukui05]_. If the cell has a nonzero Chern
number it is registered for future adaptive search.
The resulting number of found Weyl points is shown in the output at the
end of the progress message::

  CPU TIME WeylScanner::scan box:  9.17 sec done  60% ETA=15.29 sec: 4 WP candidates


After the whole |BoxMesh| was scaned a refined search takes place, in
which each resulting cell is bisected, which leads to a small Mesh
with 8 sub cells. For these sub cells a new search as described
above is performed. Any cells with nonzero Chern number are registered
for the next bisection step. The bisecting stops if the sub cell size
is smaller than *tol*. The alogrithm can search for several *homos*
at once (just give a list). During the bisections previously found
WPs can vanish. This has to do with the fine structure of the Berry
curvature and the limited accuracy of the 8-point integration of the
Berry curvature. An example of this could looks like this::

        box size=   0.0155955, 4 WPs ...
                ... finer scan of the currently 4 WPs:


        box size=   0.0077977, 3 WPs ...
                ... finer scan of the currently 3 WPs:

where two bisection steps are shown. After the first bisection one WP
is lost. This also means that a large-*tol* result is indicative of
the existence WPs (or false positives). If a WP get's lots the
algorithm will redo the bisection step with a small shift of the
cell. This seems to help in many cases. Yet, it does not guarantee
that no loss occurs.

It is in general a good idea to
try several mesh sizes (*nk=*), especially also in small
steps. E.g. we could test for *nk=* 10, 11, 20, 21, 40, 41 and so on.
Currently, there is no better way. It is **not save** to assume
that a large *nk* will detect all Weyl points! It is also very likely
that the  |BoxMesh.origin| matters. As you can see it can be quite tricky.

After the search is over the file :file:`weylpoints.py` is created
which contains a `list` of |WeylPoints|
At the end of the script (lines 63-66) this file is imported to make
the WPs available and in lines 68-80 the WPs are printed in units of :math:`\frac{2\pi}{a_i}`.

One last word about the *tol* parameter. If is is to large the
accuracy of the resulting WP will be not so great. This is
disadvantageous for future tests oft these WPs.

We will now edit the script and put *nk=10*. Then we execute::

  python -u wpsearch.py | tee out

The results are written to `stdout` and into the file :file:`out`.
You will notice that no WPs where found.

Next, set *nk* to 20 and rerun... again no WPs.

Next, set *nk* to 24 and rerun. You will have found 4 WPs after the
initial scan and also after bisections. The output looks like::


  WP[000]: k=    -0.10121221      0.02978658      0.00001523 Chi=-1.00 E=    0.055099\
     homo=56 spin=1
  WP[001]: k=    -0.10119098     -0.02980890      0.00001523 Chi= 1.00 E=    0.055192\
     homo=56 spin=1
  WP[002]: k=     0.10116975      0.02976426      0.00001523 Chi= 1.00 E=    0.055209\
     homo=56 spin=1
  WP[003]: k=     0.10119098     -0.02980890      0.00001523 Chi=-1.00 E=    0.055192\
     homo=56 spin=1
  
  ------------------------------------------------------------------------
  Slabify::findWeylPoints: END
  ------------------------------------------------------------------------
  
  WP   0:  -0.101212   0.054301   0.000061  -1.0   0.055099         56 1
  WP   1:  -0.101191  -0.054342   0.000061   1.0   0.055192         56 1
  WP   2:   0.101170   0.054260   0.000061   1.0   0.055209         56 1
  WP   3:   0.101191  -0.054342   0.000061  -1.0   0.055192         56 1


You will notice that the point coordinates, chirality and energy match
the results of [Wang15]_. You also notice that time reversal related
WPs have the same chirality (WPs 1 and 2 or 0 and 3) while WPs
connected by mirror symemtry have opposite chirality (WPs 0 and 1 or 2
and 3 ...).  In principle we could have restricted the |BoxMesh| to
:math:`x,y,z\ge 0` in lines 56-58. We did not do this to make the
script most generic.

Let's now put *nk=40* and rerun... you got 4 WPs after the initial
scan but magically found 5 more (of which #3 and #4 are identical) in
the second bisection step::

    
  WP[000]: k=    -0.10151555     -0.00971942      0.00001142 Chi=-1.00 E=    0.016153\
     homo=56 spin=1
  WP[001]: k=    -0.10151555      0.00975960      0.00001142 Chi= 1.00 E=    0.016206\
     homo=56 spin=1
  WP[002]: k=    -0.10121507     -0.02978769     -0.00001142 Chi= 1.00 E=    0.055082\
     homo=56 spin=1
  WP[003]: k=    -0.10119629      0.02978769      0.00001142 Chi=-1.00 E=    0.055135\
     homo=56 spin=1
  WP[004]: k=    -0.10117751      0.02976761      0.00001142 Chi=-1.00 E=    0.055175\
     homo=56 spin=1
  WP[005]: k=     0.10117751     -0.02976761      0.00001142 Chi=-1.00 E=    0.055175\
     homo=56 spin=1
  WP[006]: k=     0.10119629      0.02980778      0.00001142 Chi= 1.00 E=    0.055162\
     homo=56 spin=1
  WP[007]: k=     0.10155311     -0.00973951      0.00001142 Chi= 1.00 E=    0.016018\
     homo=56 spin=1
  WP[008]: k=     0.10155311      0.00973951     -0.00001142 Chi=-1.00 E=    0.016018\
     homo=56 spin=1
  
  ------------------------------------------------------------------------
  Slabify::findWeylPoints: END
  ------------------------------------------------------------------------
  
  WP   0:  -0.101516  -0.017719   0.000046  -1.0   0.016153         56 1
  WP   1:  -0.101516   0.017792   0.000046   1.0   0.016206         56 1
  WP   2:  -0.101215  -0.054303  -0.000046   1.0   0.055082         56 1
  WP   3:  -0.101196   0.054303   0.000046  -1.0   0.055135         56 1
  WP   4:  -0.101178   0.054266   0.000046  -1.0   0.055175         56 1
  WP   5:   0.101178  -0.054266   0.000046  -1.0   0.055175         56 1
  WP   6:   0.101196   0.054340   0.000046   1.0   0.055162         56 1
  WP   7:   0.101553  -0.017755   0.000046   1.0   0.016018         56 1
  WP   8:   0.101553   0.017755  -0.000046  -1.0   0.016018         56 1


What must have happend is that the two WPs at low resolution in the
initial step must have looked like a cell with chirality larger than
one. This can easily happen because of the crude 8-point formula for
the chirality of a mesh cell.


The 9 WPs are basically two different ones if we do not count symmetry
equivalent WPs (look at the coordinates and energies).  We cannot tell
why this second WP was not been reported in [Wang15]_. Either they
missed it or it is a computational differences.

The occureance of duplicates is due to the adaptive search and
especially its internal cell shift to avoid loss of WPs.

 

After finding the WPs we need to check if they are false positives.
Such things actually exist, e.g. if the berry curvature is locally
cylindrical with outwards pointing vectors, in which case tiny
roundoff errors will give a nonzero Chern number although it is not a 3d
monopole. We do this with another script, which calculates the Chern
number in a small sphere around the WP and which also produces a
picture of the monopole (hedgehog). Have a look at the script
:file:`cherninsphere.py`.

.. literalinclude:: ../../DOC/pyfplo/Examples/slabify/Weylpoints/MoTe2/slabify/wpsearch/cherninsphere.py
   :linenos:

In lines 24-27 we have the usual 3d setup. In lines 30-33 we try to
import :file:`weylpoints.py` (result of :file:`wpsearch.py`) and in
lines 35-37 we pick a particular WP and calculate the Chern number in
a sphere around the WP center. The center and homo are taken from
:file:`weylpoints.py` as well as the radius of the sphere. The radius
saved in :file:`weylpoints.py` is the size of the smallest bisection
cell in which the WP was found. This is usually a sensible radius. In
some cases it might be better to take a larger or smaller radius. If a
smaller radius is chosen the WP might actually fall outside the
sphere, which is not good. You most likely can judge this from the
graphical representation (see below). The script takes a command line
option to specify which WP to look at. Note, that WPs are numbered
starting with 0.

Let's start with WP No. 0::

  python -u cherninsphere.py -i 0

In the resulting output you find lines::

  List of chern numbers with abs value larger than 0.01

  chern up to band 56:             -0.387852580115874
  
In our case only one chern number is written (for homo 56, which is
what we where looking at).However, the printed Chern number is far
from integer. So, let's look at the picture first. The program writes
files called :file:`berrycurvsphere.net` and so forth. You need to
have the program `opendx` installed for this to work. Please execute::

  dx -execute_on_change -image berrycurvsphere.net

From the menu `Options` in the main window chose `View Control`. In
the dialog click the `Reset` button at the bottom right. Then in the
`mode` combobox select `Rotate` and use the left mouse click and drag
to rotate the picture. You will see something like
:ref:`FigHedgehog0`. (If your `opendx` works as mine does you can use
`Ctrl+F` for reset and `Ctrl+R` for rotate.)

.. _FigHedgehog0:
.. figure:: pics/Weylpoints/MoTe2/hedgehog0.jpg
   :figwidth: 8cm
   :width: 8cm
   :align: center

   The hedgehog (monopole) of WP 0.

Please note that all vectors point inside, which indicates negative
chirality which fits the data stored in :file:`weylpoints.py` and also
the sign of the Chern number written to the output.

Now, in the control panel (menu `windows` -> `Open All Control
Panels`) use the arrows on the `AutoGlyph upper limit` to
decrease/increase the vector length limit.  If you click the arrows
for a longer time the changes will speed up.  If you go to the lowest
value (left arrow) the arrow will be of the same length since each
vector of larger length will be limited in length. If you however go
to the highest limit (right arrow) essentially no limit is applied and
you see the actual length of the Berry curvature field. You will
notice that the vectors are rather long around the poles and very
short around the equator. This is the reason for the non-integer Chern
number. The curvature field is highly anisotropic. This is also an
issue for the Weyl point finding algorithm (as you noticed it is not
straight forward to find all 8 WPs). In such a situation we can
increase the number of points on the sphere *nsubdiv* in line 37.

Please increase it to *nsubdiv=20* and rerun the script... the Chern
number is -0.632.

Please increase it to *nsubdiv=40* and rerun the script... the Chern
number is -0.911.

Please increase it to *nsubdiv=100* and rerun the script... the Chern
number is -0.993.

So, we finally converged to the result, which we already knew. This is
a somewhat extreme case. The more important point is the graphical
representation. If it shows a monopole (especially for limited vector
length) the Weyl point is valid.
Btw. WP No. 3 is even more extreme.


Please revert the script to *nsubdiv=10* and check the other Weyl
points: have a look at the graphical representation. We will have
confirmed that all WPs are actually Weyl points by confirming the monopole.
Indeed, it would have sufficed to check the two WPs with different
energies since the others are obtained by symmetry.

   

With this the WP search is complete. We can proceed e.g. with surface
spectra. Please change into :file:`semi`::

  cd ../semi

Have a look at the script :file:`fs.py` which calculates the surface
spectral function for :math:`E=0` in the surface Brillouin zone (the
Fermi surface equivalent).

.. literalinclude:: ../../DOC/pyfplo/Examples/slabify/Weylpoints/MoTe2/slabify/semi/fs.py
   :linenos:

In lines 29-35 we setup a semi infinite slab with a (00-1)
surface. We anchor the surface such that a |MoTe2| layer is the last
layer. After executing the script (after the structure setup was
executed, hence before the actual calculation starts) the directory
:file:`slabifyres` contains the file :file:`=.in_final_PLlayer`. If
you load it into **xfplo** you will see the cell whose semi-infinite
repetition will make up the semi-slab. Note, that the c-axis points
downwards! Hence the surface is the lower face of the cell while the
infinte side grows on the upper cell face (not shown). If you had
chosen the (001) surface the c-axis would have pointed upwards.
Note, that *numberoflayers* is set to 1. Our 3d unit cell is already
long enough to fulfill the condition that no hopping reaches further
then between two adjacent primary layers.

In lines 37-42 we setup the axis scales of the surface BZ and in lines
44-50 we set up the surface BZ and a mesh which represents the upper
right quadrant of the surface BZ. We set the Fermi surface to -20meV
as in [Wang15]_. The imaginary part is chosen in a resonable way such
that the Lorentzian width is comparable to the mesh distance. If we
would take a much smaller imaginary part, the picture would tend to
look spotty. A larger value will smear it out.

In line 52 the actual calculation takes place. After the calculation
finishes please execute::

  xfbp fs.xpy

to produce :ref:`this Figure <Figfs>`. Compare this to Fig. 3 of
[Wang15]_ and also read the discussion of the surface state
therein. Please note, that we show the whole quadrant, while [Wang15]_
shows a smaller part of the whole surface BZ.


   
.. _Figfs:
.. figure:: pics/Weylpoints/MoTe2/fs.png
   :figwidth: 8cm
   :width: 8cm
   :align: center

   The surface Fermi surface spectral function for ::math:`E=-20` meV.
   The two crosses denote the two WPs of this quadrant. 

FInally, we could also calculate energy distribution curves. This was,
however already discussed in other examples. With this the basic Weyl
point example shall come to it's end.

