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.).
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.
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.
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.
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.
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
andb
axis perpendicular to thezaxis
. (thec
axis it not necessarily parallel to thez
axis) This creates the elementary building block for'slabs'
/'semislabs'
.- Then this block is repeated
numberoflayers
of times inc
-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 thatanchor
becomes thez=0
-plane. After this the cell is cut atz=0
andz=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 (largestz
) 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 inAl(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.
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.
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.
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.
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.
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..
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.
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.
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.