# A basic tutorial

- The bulk band structure
- The bulk Fermi surface
- Fermi surface cuts
- Bulk projected bands
- Finite slab with 10 unit cells
- Finite slab with 10 unit cells (doubled in-plane cell)
- Finite slab with 10 unit cells (doubled in-plane cell), one atom removed
- Finite slab with 10 unit cells (doubled in-plane cell), 3 atoms removed
- Semi infinite slab
- Semi infinite slab, doubled planar cell

Note

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

This example tries to explain how `pyfplo.slabify` works in detail. The tutorial files are in `FPLO.../DOC/pyfplo/Examples/slabify/model` where `FPLO...` stands for your version’s FPLO directory, e.g. `FPLO22.00-62`.

We use a hand written Hamiltonian file (`+hamdata`) containing some model data. Usually this is created by the Wannier function module of fplo. The model defines a single orbital tight binding model on a cubic lattice.

## The bulk band structure

In a first step we plot the 3d bulk band structure. The python script `3d/slabify.py` is shown in the following

     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 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='3d'
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=[
            ['$~G',[0,0,0]],
            ['X',[0.5,0,0]],
            ['M',[0.5,0.5,0]],
            ['$~G',[0,0,0]],
            ['Z',[0,0,0.5]],
            ['R',[0.5,0,0.5]],
            ['L',[0.5,0.5,0.5]],
            ['Z',[0,0,0.5]],
            ]
        bp.calculateBandPlotMesh(s.dirname)

        s.calculateBandStructure(bp,suffix='_my_suffix');

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

    if __name__ == '__main__':

        work()

Line 8 can be uncommented and edited to make python search for pyfplo in a particular location (also see Setup). Line 14 shows which `pyfplo` version was loaded and from where. If you uncomment line 16, the script is showing an error message if the `pyfplo` version does not match a particular version.

Please run the script (from `3d/README.rst`):

    # run the script as in

    ./slabify.py | tee out
    xfbp bands.cmd

    # and have a look at the output and into slabifyres/

    #alternatively run

    python ./slabify.py | tee out
    xfbp bands.cmd

Now for the actual script. Line 25 defines a convenient variable pointing to the Hamiltonian file. Line 27 makes `s` an instance of `slabify.Slabify`. In line 30 we tell it that we want a 3d object and since we dont set any other options it will be the simple cubic cell as defined in `+hamdata`. Line 33 now reads the Hamiltonian data and sets up the structure. After this step there will be several files in the output directory `slabifyres/` called `=.in_...`:

> =.in_step_1_3d_enlarged:
>   This is the cell after the first structure manipulation step, which is the application of the `enlarge` matrix.
> =.in_final_PLlayer:
>   I the final result.

Of course in our case the two are identical, since we did not give any structure options except for `s.object=3d`. Furthremore, there are bandstructure files. Ignore `+sweights_sf_my_suffix` for now.

On with the code: Line 35 make `bp` an instance of `common.BandPlot` while lines 36.. set the high symmetry points. Note, that the units are `slabify.Slabify.kscale`. One can set `kscale` to something else after the call to `s.prepare(hamdata)`. Line 46 sets up `bp` and writes the file `+points` to the directory defined in `slabify.Slabify.dirname`. If this variable is set by the user, it must be before the call to `slabify.Slabify.prepare`. Finally, in line 48 the band structure and band weights are calculated. The `suffix` argument helps to give the files custom made names.

Finally, if you sucesfully installed xfbp and executed `xfbp bands.cmd` you saw the simple bulk band structure (That is what you should see.).

Fig. 1 The 3d band structure of the model.

It is strongly encouraged to also study the `.cmd` scripts to learn tricks. If you want to use other software for plotting you can use the `pyfplo.common.BandPlot.readBands` and `pyfplo.common.BandWeights.readBandWeights` methods to deal with the data as you please.

## The bulk Fermi surface

Next we create a bulk Fermi surface of the model. We created an appropriate `=.in` with symmetry P1 and a simple cubic lattice with correct lattice constants (as in `+hamdata`). Next, we opened xfplo in fermi surface mode, defined a mesh, exported it, saved the settings in `=.xef`, stopped before the automatic fplo run and quit the program. Now, we got `=.kp` and `=.xef`.

Please change into `FS/` and run the script (from `FS/README.rst`):

    # Here we created an appropriate =.in with the correct
    # lattice. The we used xfplo -fs to setup and export a k-mesh to =.kp.
    # We use slabify to calculate the Fermisurface corresponding to the
    # model data in ../+hamdata.
    #
    # run

    ./slabify.py | tee out
    xfplo =.xef

    # to see the result.

You should see something like this.

Fig. 2 The 3d Fermi surface extracted via slabify.

Let’s explain the script `FS/slabify.py` now.

     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
    #! /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

    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()
        # set output directory to current directory
        s.dirname='.'

        s.object='3d'
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.readBandPlotMesh('./=.kp')

        s.calculateBandStructure(bp);

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

    if __name__ == '__main__':

        work()

In line 31 we explicitely set `slabify.Slabify.dirname` to the local directory. This will put everything into the same directory where `slabify.py` was executed. The setup continues as in The bulk band structure until line 42 where instead of setting a path through high symmetry points the k-point file from xfplo is used.

## Fermi surface cuts

There is an option to create Fermi surface cuts.

Please change into `cuts/` and run the script (from `cuts/README.rst`):

    #run

    ./slabify.py | tee out
    xfbp cuts.cmd

    #have a look at slabifyres/ and cuts.png

You should see something like this.

Fig. 3 Cuts through the 3d Fermi surface of the model.

Let’s explain the script `cuts/slabify.py` now.

     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
    #! /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

    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='3d'
        s.printStructureSettings()

        s.prepare(hamdata)

        fso=sla.FermiSurfaceOptions()
        fso.setMesh(100,[-0.5,0.5],100,[-0.5,0.5])

        n=10
        for ikz in range(0,n+1):
            kz=ikz*0.5/n
            fso.setPlane([1,0,0],[0,1,0],[0,0,kz])
            fso.fermienergy=0
            suffix='_kz={0:05.2f}'.format(kz)
            # we need to set forcerecalculation=True since
            # the k-plane changes with every loop
            s.calculateFermiSurfaceCuts(fso,suffix=suffix,forcerecalculation=True);

            fso.fermienergy=-3
            suffix='_kz={0:05.2f}_ef={1:05.2f}'.format(kz,fso.fermienergy)
            # Now we do NOT need to set forcerecalculation=True since
            # only the Fermi energy changed.
            # The actual example is so fast that you do not see the impact.
            s.calculateFermiSurfaceCuts(fso,suffix=suffix,forcerecalculation=False);

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

    if __name__ == '__main__':

        work()

In lines 35-36 we set up a `slabify.FermiSurfaceOptions` instance with default axes and define a 100x100 2d mesh which stretches from `-0.5` to `0.5` in units of `slabify.Slabify.kscale` which happens to be `2*pi/a` by default.

Next a loop is done over various `kz`-values. In line 41 the plane (axes and origin) is specified explicitely. The origin `[0,0,kz]` puts this plane parallel to `kx`, `ky` through the `kz`-value. Line 42 sets the Fermi energy and line 43 a filename suffix. In line 46 `Slabify.calculateFermiSurfaceCuts` is called with the flag `forcerecalculation=True`, which is explained under the link above. After this line the files

> `+cut_band_sf`,

> `+cut_bweights_sf` and

> `+cuts_kz=00.00.spin1_kz=...`

are created in `slabifyres/`.

Line 49 sets a different Fermi energy and line 54 recalculates the files `+cuts_kz=00.00.spin1_kz=...` (iso lines) but not the other files because of `forcerecalculation=False`. In a real example with lots of orbitals this avoids the expensive diagonalization step. There are also `+cutswithweights` files if the wds option is specified.

## Bulk projected bands

On can calculate the bulk projected bands (BPB) as energy distribution curves (EDC) or Fermi surface projections.

Please change into `bpb/` and run the script (from `bpb/README.rst`):

    #run

    ./slabify.py | tee out
    xfbp bpbedc.cmd
    xfbp bpbfs.cmd

    #have a look at slabifyres/ and *.png

You should see something like this.

Fig. 4 The bulk projected bands of the model.

Let’s explain the script `bpb/slabify.py` now.

     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
    #! /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

    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='3d'
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=[
            ['$~G',[0,0,0]],
            ['X',[0.5,0,0]],
            ['M',[0.5,0.5,0]],
            ['$~G',[0,0,0]],
            ]
        bp.ndiv=100
        bp.calculateBandPlotMesh(s.dirname)

        Ne=200
        Nkz=100

        ec=sla.EnergyContour(Ne,-15,25)
        print( ec)
        s.calculateBulkProjectedEDC(bp,ec,[0,0,1],nz=Nkz)

        Nk=200
        fso=sla.FermiSurfaceOptions()
        fso.setMesh(Nk,[-0.5,0.5],Nk,[-0.5,0.5])
        fso.setPlane([1,0,0],[0,1,0],[0,0,0])
        fso.fermienergyim=2./Nk*10

        # now make projected fermi surfaces
        s.calculateBulkProjectedFS(fso,zaxis=[0,0,1],nz=Nkz)

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

    if __name__ == '__main__':

        work()

In line 34-42 we define a path through the 2d BZ and set the maximum number of subdivisions for the path segments to 150. In line 47 we define an `EnergyContour` for the EDC. In line 49 the bulk projected EDC is calculated (`Slabify.calculateBulkProjectedEDC`) by defining a projection axis zaxis and setting the number of integration intervals for the projetion. The code finds out which is the shortest period in the projection direction and integrates over this period via linear interpolation. A complex representation of the delta functions is used, which necessitates an imaginary part in `EnergyContour` (automatic in our case). The result is written to `slabifyres/+bpb_fs.spin1`.

In lines 52-59 we set up a 2d mesh in a plane perpendicular to the projection axis, set an imaginary energy part (needs some experimenting for good results) and do the calculation (`Slabify.calculateBulkProjectedFS`). The result is written to ‘slabifyres/+bpfs_sf.spin1’.

## Finite slab with 10 unit cells

We construct a finite slab with 10 unit cells. Please change into `slab10/` and run the script (from `slab10/README.rst`):

    # run the script as in

    ./slabify.py | tee out
    xfbp weigths.cmd

    # and have a look at the output and into slabifyres/

    #alternatively run

    python ./slabify.py | tee out
    xfbp weigths.cmd

You should see something like this.

Fig. 5 The 10-cell finite slab band structure.

Let’s explain the script `slab10/slabify.py` now.

     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
    #! /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

    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='slab'
        s.anchor=-0.001
        s.numberoflayers=10
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=[
            ['$~G',[0,0,0]],
            ['X',[0.5,0,0]],
            ['M',[0.5,0.5,0]],
            ['$~G',[0,0,0]],
            ]
        bp.calculateBandPlotMesh(s.dirname)

        s.calculateBandStructure(bp);

        # all orbitals up to 10 length units at the upper surface
        upper=s.orbitalNamesByDepth(-1,10)
        # all orbitals up to 10 length units at the lower surface
        lower=s.orbitalNamesByDepth(10,-1)
        # a nice python trick to get the rest list
        rest=list(set(s.orbitalNames())-set(upper)-set(lower))

        wds=sla.WeightDefinitions()
        wds.add('bulk').addLabels(rest)
        wds.add('lower').addLabels(lower)
        wds.add('upper').addLabels(upper)

        bw=sla.BandWeights(s.dirname+'/+bweights_sf')
        bw.addWeights(wds,s.dirname+'/+bwsum_sf')

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

    if __name__ == '__main__':

        work()

This time we make a slab structure (free standing x,y-periodic slab) in line 30. The `zaxis` is not set and hence has its default value `[0,0,1]`. Line 32 sets `numberoflayers`. We set the `anchor` to slightly below 0 in line 31. What happens now:

> - The default 3d unit cell is enlarged to produce the “enlarged 3d” unit cell in `=.in_step_1_3d_enlarged` (identical to the default here).
> - The enlarged 3d unit cell will be transformed into a 3d unit cell which has the `a` and `b` axis perpendicular to the `zaxis`. (the `c` axis it not necessarily parallel to the `z` axis) This creates the elementary building block for `'slabs'`/`'semislabs'`.
> - Then this block is repeated `numberoflayers` of times in `c`-direction to get the 3d “layered cell” (`=.in_step_2_3d_layered`).
> - Next, the “layered cell” will be anchored, which means that the `z`-periodic 3d “layered cell” will be shifted such that `anchor` becomes the `z=0`-plane. After this the cell is cut at `z=0` and `z=1`, which results in a free standing slab of the length of one “layered cell” unit cell. The result is found in `=.in_step_3_slab_anchored` For `'semislab'`s the resulting cell is used as surface block (layer) such that the upper boundary (largest `z`) represents the surface to the vaccum. Internally this block is repeated indefinitely at the lower boundary of the unit cell to form a semi-infinite slab.

The next option is only needed for `'slab'`s.

> - (Not used here) After this, `cutlayersat` is applied to cut away some layers at the top and bottom of the slab.

The next option is currently only used for `'slab'`s, which means that we cannot remove selected atoms for `'semislab'`s yet.

> - (Not used here) Finally, all atoms in the list `cutatoms` are removed from the slab. The number in this list must be taken from the site numbers in `=.in_step_4_slab_cutlayers`. Use xfplo and point the mouse at the atoms: the status bar will show the site number as in `Al(S3,W3,T3)...`. `S3` means site 3. The result after this step is `=.in_final_PLlayer`

We setup the high symmetry points an calculate bands and weights up to line 46.

What comes next is adding all orbitals at the two sides and in the middle to get band weights to color the band structure with. In the current slab this is not very interesting. But it becomes interesting later. Line 51 defines a list of orbitals which sit up to 10 length units (usually Bohr radii) deep at the upper side of the slab. Note the `-1` in the first argument. This is explained in `orbitalNamesByDepth` and `orbitalIndicesByDepth`. It just takes orbitals at the lower side which lie below the lower vaccum boundary, which is no orbitals. The same is done for the `lower` variable. The line 55 shows a neat trick to calculate the rest (which is not orbitals in the upper or lower set of orbitals).

Line 58-61 defines three added weights named `'bulk'`, `'lower'` and `'upper'` and line 63-64 read the bandweights file and write the added weights to `+bwsum_sf`.

Now, we can color the band weights according to their “surface character” which is shown in The 10-cell finite slab band structure.. Nothing interesting to see here. (But later!)

## Finite slab with 10 unit cells (doubled in-plane cell)

We construct a finite slab with 10 unit cells, doubled in the `a`, `b`-plane.

Please change into `slab10x2/` and run the script (from `slab10x2/README.rst`):

    # run the script as in

    ./slabify.py | tee out
    xfbp weigths.cmd

    # and have a look at the output and into slabifyres/

    #alternatively run

    python ./slabify.py | tee out
    xfbp weigths.cmd

You should see something like this.

Fig. 6 The 10-cell slab with doubled in-plane unit cell.

Let’s explain the script `slab10x2/slabify.py` now.

     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
    #! /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

    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='slab'
        s.enlarge=[ [1,1,0],[-1,1,0],[0,0,1]]
        s.anchor=-0.001
        s.numberoflayers=10
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=[
            ['$~G',[0,0,0]],
            ['X',[0.5,0,0]],
            ['M',[0.5,0.5,0]],
            ['$~G',[0,0,0]],
            ]
        bp.calculateBandPlotMesh(s.dirname)

        s.calculateBandStructure(bp);

        # all orbital indices up to 10 length units at the upper surface
        iupper=s.orbitalIndicesByDepth(-1,10)
        # all orbital indices up to 10 length units at the lower surface
        ilower=s.orbitalIndicesByDepth(10,-1)
        # get all orbital indices
        iall=s.orbitalIndicesByDepth()
        # a nice python trick to get the rest list
        irest=list(set(iall)-set(iupper)-set(ilower))

        print( 'the orbital indices lists:',ilower,irest,iupper)

        wds=sla.WeightDefinitions()
        wds.add('bulk').addLabels(s.orbitalNames(irest))
        wds.add('lower').addLabels(s.orbitalNames(ilower))
        wds.add('upper').addLabels(s.orbitalNames(iupper))

        bw=sla.BandWeights(s.dirname+'/+bweights_sf')
        bw.addWeights(wds,s.dirname+'/+bwsum_sf')

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

    if __name__ == '__main__':

        work()

We now used an `enlarge` matrix to double the 3d unit cell (Line 31). In contrast to the previous example we now use another route to get the neede orbital sets (Line 52-60). This also alters the line 66-68. Whatever method you use depends on what you want to do. The orbital-name based approach allows specific filtering, which was shown elsewhere.

This example was not very interesting either. Let’s go on.

## Finite slab with 10 unit cells (doubled in-plane cell), one atom removed

We construct a finite slab with 10 unit cells, doubled in the `a`, `b`-plane with one atom removed at the upper side.

Please change into `slab10x2mod1/` and run the script (from `slab10x2mod1/README.rst`):

    # run the script as in

    ./slabify.py | tee out
    xfbp weigths.cmd
    xfbp sweigths.cmd

    # and have a look at the output and into slabifyres/

    #alternatively run

    python ./slabify.py | tee out
    xfbp weigths.cmd
    xfbp sweigths.cmd

You should see something like this.

Fig. 7 10-cell slab, double in-plane cell, one atom removed.

Now, that is interesting. We see a surface state steming from the removed atom at the upper boundary.

Let’s explain the script `slab10x2mod1/slabify.py` now.

     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
    91
    92
    93
    94
    95
    #! /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

    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='slab'
        s.enlarge=[ [1,1,0],[-1,1,0],[0,0,1]]
        s.anchor=-0.001
        s.numberoflayers=10
        s.cutatoms=[20]
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=[
            ['$~G',[0,0,0]],
            ['X',[0.5,0,0]],
            ['M',[0.5,0.5,0]],
            ['$~G',[0,0,0]],
            ]
        bp.calculateBandPlotMesh(s.dirname)

        s.calculateBandStructure(bp);

        # all orbitals up to 10 length units at the upper surface
        upper=s.orbitalNamesByDepth(-1,10)
        # all orbitals up to 10 length units at the lower surface
        lower=s.orbitalNamesByDepth(10,-1)
        # a nice python trick to get the rest list
        rest=list(set(s.orbitalNames())-set(upper)-set(lower))

        wds=sla.WeightDefinitions()
        wds.add('bulk').addLabels(rest)
        wds.add('lower').addLabels(lower)
        wds.add('upper').addLabels(upper)

        bw=sla.BandWeights(s.dirname+'/+bweights_sf')
        bw.addWeights(wds,s.dirname+'/+bwsum_sf')

        # now we first create  BandWeights
        bw=sla.BandWeights(s.dirname+'/+sweights_sf')
        # such that we can read its header and get the labels
        labels=bw.header().labels
        l10=['band000000010']
        # and we get all but band number 10
        rest=list(filter(lambda x: x not in l10  ,labels))
        wds=sla.WeightDefinitions()
        wds.add('rest').addLabels(rest)
        wds.add('b10').addLabels(l10)

        print( wds)

        bw.addWeights(wds,s.dirname+'/+swsum_sf')

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

    if __name__ == '__main__':

        work()

Now, we remove atom 20 in line 34. Have a look at `=.in_final_PLlayer`, you will see that an atom at the top is missing.

In lines 72-84 we do something new. We define added band weights for the file `+sweights_sf`, which contains fatband data of layer-coordinate versus k-path. The band weights are refering to the bands not orbitals. This file tells us in which layer a particular bands has what weight. By defining added weights for band 10 and the rest we can see in this Figure that the surface state in 10-cell slab, double in-plane cell, one atom removed. really lives in the first layer and the rest of the bands does not.

Fig. 8 Layer resolved weights of the bands.

## Finite slab with 10 unit cells (doubled in-plane cell), 3 atoms removed

We construct a finite slab with 10 unit cells, doubled in the `a`, `b`-plane with one atom removed at the upper side and two at the lower side.

Please change into `slab10x2mod2/` and run the script (from `slab10x2mod2/README.rst`):

    # run the scripts as in

    ./slabify.py | tee out
    ./cuts.py | tee out
    xfbp weigths.cmd
    xfbp sweigths.cmd
    xfbp wcuts.cmd

    # and have a look at the output and into slabifyres/

    # alternatively run

    python ./slabify.py | tee out
    python ./cuts.py | tee out
    xfbp weigths.cmd
    xfbp sweigths.cmd
    xfbp wcuts.cmd

You should see something like this.

Fig. 9 10-cell slab, doubled in-plane, 3 atoms removed.

Now, that is even more interesting. We see a surface state steming from the removed atom at the upper surface and two surface states from the atom removal at the lower boundary.

Let’s explain the script `slab10x2mod2/slabify.py` now.

     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
    91
    92
    93
    94
    95
    96
    97
    98
    99
    #! /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

    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='slab'
        s.enlarge=[ [1,1,0],[-1,1,0],[0,0,1]]
        s.anchor=-0.001
        s.numberoflayers=10
        s.cutatoms=[2,4,20]
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=[
            ['$~G',[0,0,0]],
            ['X',[0.5,0,0]],
            ['M',[0.5,0.5,0]],
            ['$~G',[0,0,0]],
            ]
        bp.calculateBandPlotMesh(s.dirname)

        s.calculateBandStructure(bp);

        # all orbitals up to 10 length units at the upper surface
        upper=s.orbitalNamesByDepth(-1,10)
        # all orbitals up to 10 length units at the lower surface
        lower=s.orbitalNamesByDepth(10,-1)
        # a nice python trick to get the rest list
        rest=list(set(s.orbitalNames())-set(upper)-set(lower))

        wds=sla.WeightDefinitions()
        wds.add('bulk').addLabels(rest)
        wds.add('lower').addLabels(lower)
        wds.add('upper').addLabels(upper)

        bw=sla.BandWeights(s.dirname+'/+bweights_sf')
        bw.addWeights(wds,s.dirname+'/+bwsum_sf')

        # now we first create  BandWeights
        bw=sla.BandWeights(s.dirname+'/+sweights_sf')
        # such that we can read its header and get the labels
        labels=bw.header().labels

        # and we get addweights for bands 8, 9, 10 and the rest

        ibands=[8,9,10]

        bands=['band{0:09d}'.format(x) for x in ibands]
        rest=list(filter(lambda x: x not in bands,labels))
        wds=sla.WeightDefinitions()
        wds.add('rest').addLabels(rest)
        for ib,i in enumerate(ibands):
            wds.add('b{0:02}'.format(i)).addLabels([bands[ib]])

        bw.addWeights(wds,s.dirname+'/+swsum_sf',vlevel=sla.Vlevel.All)

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

    if __name__ == '__main__':

        work()

Line 34 now removes sites 2, 4 and 20. In lines 80-88 we use a different technique to extract the layer-wise fatbands for 3 bands and the rest resulting in this Figure.

Fig. 10 Layer resolved weights of the bands.

The second script in this example is `slab10x2mod2/cuts.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
    #! /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

    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='slab'
        s.enlarge=[ [1,1,0],[-1,1,0],[0,0,1]]
        s.anchor=-0.001
        s.numberoflayers=10
        s.cutatoms=[2,4,20]
        s.printStructureSettings()

        s.prepare(hamdata)

        # all orbitals up to 10 length units at the upper surface
        upper=s.orbitalNamesByDepth(-1,10)
        # all orbitals up to 10 length units at the lower surface
        lower=s.orbitalNamesByDepth(10,-1)
        # a nice python trick to get the rest list
        rest=list(set(s.orbitalNames())-set(upper)-set(lower))

        wds=sla.WeightDefinitions()
        wds.add('bulk').addLabels(rest)
        wds.add('lower').addLabels(lower)
        wds.add('upper').addLabels(upper)

        fso=sla.FermiSurfaceOptions()
        fso.setMesh(200,[-0.5,0.5],200,[-0.5,0.5])
        fso.setPlane([1,0,0],[0,1,0],[0,0,1])
        fso.fermienergy=0
        s.calculateFermiSurfaceCuts(fso,wds=wds,suffix='E=0',
                                    forcerecalculation=True);

        fso.fermienergy=-1.
        s.calculateFermiSurfaceCuts(fso,wds=wds,suffix='E=-1.',
                                    forcerecalculation=False);

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

    if __name__ == '__main__':

        work()

Here we first define the addweights in lines 40-51. Then we give wds as argument to `Slabify.calculateFermiSurfaceCuts` which triggers the creation of the files `+cutswithweights...`. Note, that the flag `forcerecalculation=False` is used in the second call to `Slabify.calculateFermiSurfaceCuts`. The resulting files are plottet via `xfbp wcuts.cmd` and produce this Figure..

Fig. 11 Fermi surface cuts colored by bulk, lower and upper weights.

You can see that there are a lots of tiny wiggles in the bulk band region. This is so because we backfolded the planar unit cell. The only changes to the Hamiltonian are the cut-out atoms, which do not influence the bulk much. Together with the finit number of layers this is what you get. For infinite layers these bulk regions would fill up completely. What you also see however are the blue upper-surface surface-band for \(E=0\) and additionally one green lower-surface surface-band for \(E=-1\). Compare to the band weights.

## Semi infinite slab

Finally we create spectral density plots for a semi infinite slab. In this case not much interesing is seen, since the surface states where a consequence of cutting out and atom, which currently is not yet implemented for `'semislab'`s yet.

Please change into `semi/` and run the script (from `semi/README.rst`):

    # Here we set up a semi infinte slab with a normal in-plane unit cell.
    # run
    ./slabify.py | tee out
    xfbp edc.cmd
    xfbp fssemi.cmd

    # have a look at slabifyres/=.in_step... and slabifyres/=.in_final_PLlayer.
    # Try to understand how the structure settings work.

You should see something like this and something like this. Compare this to the bulk projected bands.

Fig. 12 Energy distribution curves for the semi-infinite slab.

Fig. 13 Surface Fermi surface spectral function.

Let’s explain the script `semi/slabify.py` now.

     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
    #! /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

    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.anchor=-0.001
        s.numberoflayers=4
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=[
            ['$~G',[0,0,0]],
            ['X',[0.5,0,0]],
            ['M',[0.5,0.5,0]],
            ['$~G',[0,0,0]],
            ]
        bp.ndiv=60
        bp.calculateBandPlotMesh(s.dirname)

        ec=sla.EnergyContour(200,-15,25)
        print( ec)

        bp.on()

        # penetration one block (primary layer)
        s.calculateEDC(bp,ec,penetrationdepth=-1,suffix='_1block')

        # penetration 10 blocks (primary layer)
        s.calculateEDC(bp,ec,penetrationdepth=-10,suffix='_10block')

        fso=sla.FermiSurfaceOptions()
        fso.setMesh(200,[-0.5,0.5],200,[-0.5,0.5])
        s.calculateFermiSurfaceSpectralDensity(fso,penetrationdepth=-1,
                                               suffix='_1block')
        s.calculateFermiSurfaceSpectralDensity(fso,penetrationdepth=-10,
                                               suffix='_10block')

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

    if __name__ == '__main__':

        work()

In line 32 the structure type is set. We u se `numberoflayers=4` in this case, although it is not really needed. See documentation of `numberoflayers`. In lines 39-47 a band structure path is setup for the energy distribution curve (EDC). Line 49 defines an energy contour (see `EnergyContour`) with an automatic imaginary part. Line 55 and 58 actually execute the EDC calculation (`Slabify.calculateEDC`) with two different `penetrationdepth`s. The resulting files are called `+akbl_sf.spin1` (or `spin2`), which we modified via the suffix argument. The name means: `akbl`=\(A_\mathrm{Bl}(\boldsymbol{k})\) = Bloch-spectral-function and `sf` = slabify.

`penetrationdepth` indicates to which depth the spectral density is collected. The larger the `penetrationdepth` the more bulk signal will be sampled.

Lines 60-65 setup and calculate the surface Fermi surface, which is basically the k,k-resolved spectral density for the 2d surface BZ (`Slabify.calculateFermiSurfaceSpectralDensity`). The resulting files are called `+fs_sf.spin1` (or `spin2`), which we modified via the suffix argument. The name means `fs` = Fermi-Surface-Spectral-Function and `sf` = slabify.

## Semi infinite slab, doubled planar cell

At very last, we double the planar unit cell for illustrative purpose.

Please change into `semix2/` and run the script (from `semix2/README.rst`):

    # Here we set up a semi infinte slab with a normal in-plane unit cell.
    # run
    ./slabify.py | tee out
    xfbp edc.cmd
    xfbp fssemi.cmd

    # have a look at slabifyres/=.in_step... and slabifyres/=.in_final_PLlayer.
    # Try to understand how the structure settings work.

You should see something like this and somehting like this. Compare this to the left panel of the finite slab result with removed atoms. It demonstrates the backfolding in the bulk projected bands region.

Fig. 14 The energy distribution curves.

Fig. 15 Surface Fermi-surface spectral function.

Let’s explain the script `semix2/slabify.py` now.

     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
    #! /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

    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'
        # make a larger 3d cell out of the simple cell
        s.enlarge=[ [1,1,0],[-1,1,0],[0,0,1]]
        s.anchor=-0.001
        s.numberoflayers=4
        s.printStructureSettings()

        s.prepare(hamdata)

        bp=sla.BandPlot()
        bp.points=[
            ['$~G',[0,0,0]],
            ['X',[0.5,0,0]],
            ['M',[0.5,0.5,0]],
            ['$~G',[0,0,0]],
            ]
        bp.ndiv=60
        bp.calculateBandPlotMesh(s.dirname)

        ec=sla.EnergyContour(200,-15,25)
        print( ec)

        bp.on()

        # penetration one block (primary layer)
        s.calculateEDC(bp,ec,penetrationdepth=-1,suffix='_1block')

        # penetration 10 blocks (primary layer)
        s.calculateEDC(bp,ec,penetrationdepth=-10,suffix='_10block')

        fso=sla.FermiSurfaceOptions()
        fso.setMesh(200,[-0.5,0.5],200,[-0.5,0.5])
        s.calculateFermiSurfaceSpectralDensity(fso,penetrationdepth=-1,
                                               suffix='_1block')
        s.calculateFermiSurfaceSpectralDensity(fso,penetrationdepth=-10,
                                               suffix='_10block')

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

    if __name__ == '__main__':

        work()

The only thing we changed was line 34.

With this the introductory tutorial shall end.