A basic tutorial

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

_images/bands3d.png

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.

_images/FS3d.jpg

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.

_images/cuts.png

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.

_images/bpb.png

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.

_images/weislab10.png

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.

_images/weislab10x2.png

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.

_images/weislab10x2mod1.png

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.

_images/sweislab10x2mod1.png

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.

_images/weislab10x2mod2.png

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.

_images/sweislab10x2mod2.png

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

_images/wcuts.png

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.

_images/edc2.png

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

_images/fssemi.png

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

_images/edc3.png

Fig. 14 The energy distribution curves.

_images/fssemi1.png

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.