[Wang15] | (1, 2, 3, 4, 5, 6) Zhijun Wang et.a. arxiv:1511.07440 (2015)

[Fukui05] | [Takahiro Fukui et.al. J. Phys. Soc. Jpn. 74, 1674 (2005)]

# Weyl semi metals

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 `FPLO.../DOC/pyfplo/Examples/slabify/Weylpoints/MoTe2`. They contain a converged full relativistic calculation.

Figure The bulk unit cell and band structure of MoTe2. 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.

Fig. 33 The bulk unit cell and band structure of MoTe2.

With a right mouse click close to the highest and lowest band in xfbp (when `+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 `makewandef.py` which has the following content

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    #! /usr/bin/env python
    from __future__ import print_function
    import pyfplo.wanniertools as wt

    # ===================================================================
    def main():

        wdc=wt.WanDefCreator(rcutoff=25,wftol=0.001,coeffformat='bin',
                             wfgriddirections=[[1,0,0],[0,1,0],[0,0,1]],
                             wfgridsubdiv=[1,1,1],savespininfo=False,
                             keeprunning=True,opendxinterface=False)

        emin=-6
        emax= 4
        delower=1
        deupper=0

        wdc.add(wt.MultipleOrbitalWandef('Mo',range(1,4+1),['4db'],
                                         emin=emin,emax=emax,
                                         delower=delower,deupper=deupper))
        wdc.add(wt.MultipleOrbitalWandef('Te',range(5,12+1),['5pb'],
                                         emin=emin,emax=emax,
                                         delower=delower,deupper=deupper))

        wdc.writeFile('=.wandef')

        return

    # ===================================================================
    # ===================================================================
    # ===================================================================
    # ===================================================================
    if __name__ == '__main__':
        main()

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 `=.wandef`. We run fplo now to obtain the necessary data for the WF calculation in `+wancoeff`:

    fplo.... > out

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

    fplo.... > outwf

Now, the 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 The WF fit.. 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.

Fig. 34 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 `slabify/3d`:

    cd slabify/3d

and execute the two following commands:

    python bands.py
    xfbp bands3d.xpy

with the result The 3d bulk band structure from the WF model..

Fig. 35 The 3d bulk band structure from the WF model.

Have a look at the script `bands.py`

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    #! /usr/bin/env python

    from __future__ import print_function
    import sys

    # If your pyfplo is not found you could also
    # explicitly specify the pyfplo version path:
    #sys.path.insert(0,"/home/magru/FPLO/FPLO22.00-62/PYTHON/doc");

    import numpy as np
    import pyfplo.slabify as sla
    import pyfplo.fploio as fploio

    print('\npyfplo version=: {0}\nfrom: {1}\n'.format(sla.version,sla.__file__))
    # protect against wrong version
    #if fedit.version!='22.00': raise RuntimeError('pyfplo version is incorrect.')

    # ===================================================================
    #
    # ===================================================================

    def work():

        p=fploio.INParser()
        p.parseFile('../../=.in')
        d=p()('special_sympoints')
        l=[]
        for i in range(d.size()):
            l.append([ d[i]('label').S,d[i]('kpoint').listD])

        hamdata='../../+hamdata'

        s=sla.Slabify()
        s.object='3d'
        s.printStructureSettings()
        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=l
        bp.calculateBandPlotMesh(s.dirname)
        s.calculateBandStructure(bp)

    # ===================================================================
    #
    # ===================================================================

    if __name__ == '__main__':

        work()

In line 27-32 we read the high symmetry points from the fplo input 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 `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 `../wpsearch`:

    cd ../wpsearch

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

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    #! /usr/bin/env python
    # ===================================================================
    # file:   auto.py
    # author: k.koepernik@ifw-dresden.de
    # date:   19 Jun 2017
    from __future__ import print_function
    import sys
    import numpy as np
    import numpy.linalg as LA

    import pyfplo.slabify as sla

    print( '\npyfplo version=: {0}\nfrom: {1}\n'.format(sla.version,sla.__file__))

    # ===================================================================
    #
    # ===================================================================
    def work():

        hamdata='../../+hamdata'
        s=sla.Slabify()
        s.object='3d'
        s.prepare(hamdata)

        homo=56

        nk=10
        tol=1e-4

        # setup primitive reciprocal cell BoxMesh
        G=s.hamdataRCell()
        print( 'G=\n',G)
        x=G[:,0]
        y=G[:,1]
        z=G[:,2]
        dx=LA.norm(x)
        dy=LA.norm(y)
        dz=LA.norm(z)
        # with most isotropic mesh subdivision
        nx=max(2,nk)
        ny=max(2,int(dy/dx*nx))
        nz=max(2,int(dz/dx*nx))

        # shift the origin by a small amount
        k0=-(x/(nx-1)+y/(ny-1)+z/(nz-1))*0.5

        # define the BoxMesh
        box=sla.BoxMesh()
        box.setBox(xaxis=x,yaxis=y,zaxis=z,origin=k0)
        box.setMesh(nx=nx,xinterval=[-dx/2,dx/2],
                    ny=ny,yinterval=[-dy/2,dy/2],
                    nz=nz,zinterval=[-dz/2,dz/2])

        # and try to find Weyl points
        s.findWeylPoints(box,[homo],tol)

        try:
            from weylpoints import wps
        except:
            raise RuntimeError('file weylpoints.py does not exist')

        # the conventional cell
        A=s.hamdataCCell()
        # its lattice parameters
        a=LA.norm(A[:,0])
        b=LA.norm(A[:,1])
        c=LA.norm(A[:,2])
        boa=b/a
        coa=c/a
        # print in units of 2p/a_i
        for i,w in enumerate(wps):
            print( 'WP{0:4d}: {1:10.6f} {2:10.6f} {3:10.6f} {4:5.1f} {5:10.6f} {6:10d} {7}'\
                .format(i,w.k[0],w.k[1]*boa,w.k[2]*coa,w.chirality,w.energy,
                        w.homo,w.spin))

        return

    # ===================================================================
    #
    # ===================================================================

    if __name__ == '__main__':
        work()

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 \(\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 `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 `origin` matters. As you can see it can be quite tricky.

After the search is over the 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 \(\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 `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 \(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 `cherninsphere.py`.

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    #! /usr/bin/env python
    # ===================================================================
    # file:   auto.py
    # author: k.koepernik@ifw-dresden.de
    # date:   19 Jun 2017
    from __future__ import print_function
    import sys
    import numpy as np
    import numpy.linalg as LA

    import pyfplo.slabify as sla

    print('\npyfplo version=: {0}\nfrom: {1}\n'.format(sla.version,sla.__file__))

    # ===================================================================
    #
    # ===================================================================
    def work(iwp):

        hamdata='../../+hamdata'
        s=sla.Slabify()
        s.object='3d'
        s.prepare(hamdata)

        try:
            from weylpoints import wps
        except:
            raise RuntimeError('file weylpoints.py does not exist')

        wp=wps[iwp]
        s.calculateChernNumberInSphere(center=wp.k,
                                       radius=wp.radius,homo=wp.homo,
                                       nsubdiv=10,nradius=1)

        return

    # ===================================================================
    #
    # ===================================================================

    if __name__ == '__main__':

        if len(sys.argv)>1 and sys.argv[1]=='-h':
            print("usage: {} [-h] -i wp-number".format(sys.argv[0]))
            sys.exit(0)
        try:
            i=int(sys.argv.index("-i"))
            iwp=int(sys.argv[i+1])
        except:
            raise RuntimeError('need option -i wp-number')

        work(iwp)

        sys.exit(0)

In lines 24-27 we have the usual 3d setup. In lines 30-33 we try to import `weylpoints.py` (result of `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 `weylpoints.py` as well as the radius of the sphere. The radius saved in `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 `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 The hedgehog (monopole) of WP 0.. (If your `opendx` works as mine does you can use `Ctrl+F` for reset and `Ctrl+R` for rotate.)

Fig. 36 The hedgehog (monopole) of WP 0.

Please note that all vectors point inside, which indicates negative chirality which fits the data stored in `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 `semi`:

    cd ../semi

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

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    #! /usr/bin/env python

    from __future__ import print_function
    import sys

    # If your pyfplo is not found you could also
    # explicitly specify the pyfplo version path:
    #sys.path.insert(0,"/home/magru/FPLO/FPLO22.00-62/PYTHON/doc");

    import numpy as np
    import numpy.linalg as LA
    import pyfplo.slabify as sla

    print( '\npyfplo version=: {0}\nfrom: {1}\n'.format(sla.version,sla.__file__))
    # protect against wrong version
    #if fedit.version!='22.00': raise RuntimeError('pyfplo version is incorrect.')

    # ===================================================================
    #
    # ===================================================================

    def work():

        hamdata='../../+hamdata'

        s=sla.Slabify()
        s.object='semislab'
        s.zaxis=[0,0,-1]
        s.numberoflayers=1
        s.anchor=0.001
        s.printStructureSettings()
        s.prepare(hamdata)

        # the conventional cell
        A=s.hamdataCCell()
        # its lattice parameters
        a=LA.norm(A[:,0])
        b=LA.norm(A[:,1])
        aob=a/b

        Nk=200
        fso=sla.FermiSurfaceOptions()
        fso.setMesh(Nk,[0,0.5],Nk,[0,0.5*aob])
        fso.setPlane([1,0,0],[0,1,0],[0,0,0])
        fso.fermienergy=-0.02
        fso.fermienergyim=1./Nk*1.
        print( fso)

        s.calculateFermiSurfaceSpectralDensity(fso,penetrationdepth=-1.)

    # ===================================================================
    #
    # ===================================================================

    if __name__ == '__main__':

        work()

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 `slabifyres` contains the 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 this Figure. 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.

Fig. 37 The surface Fermi surface spectral function for :\(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.