pyfplo.slabify

This is a collection of routines to map Wannier Hamiltonians onto larger structures, create Fermisurfaces, cuts, spectral densities and more. Some objects which are needed are defined in pyfplo.common but are also accessible from this module via aliases of the same name.

To get help on all of those objects use:

import pyfplo.common as com
help(com)

or if pydoc is installed: pydoc pyfplo.common

Please have a look at the Examples.

To better understand the structure manipulation performed consult Structure Manipulation Algorithm.

Slabify

class Slabify

The main object to access all procedures. There are low level procedures to extract the tight-binding Hamiltonian directly as well as high level routines, which perform various tasks.

Note

do not use =.in_... files from Slabify output in context of xfplo fermi surfaces.

Here is how Slabify works.

The structure is taken from the Wannier Hamiltonian file, which usually is called +hamdata. It also contains the Hamiltonin matrix elements and optionally spin operator matrix elements.

There are several options to change the structure into something new. After the structure setup the Hamiltonian matrix elements are mapped onto this new structure.

Three structure types are available (see object):

'3d' (bulk), 'slab' (finite slab) and 'semislab' (semi inifinite slab).

Depending in the structure type various operations can be performed to obtain the new structure. The resulting structure after each particular operation is saved in =.in_... files in the local result directory (stored in Slabify.dirname, default 'slabifyres'). The operations available for the different structure types are listed next:

  • 3d/slab/semislab

    enlargement

    Use the matrix in enlarge to construct a larger 3d lattice basis.

  • slab/semislab:

    layering

    Define zaxis to determine the direction perpendicular to the surface. The 3d cell will be transformed such that the a- and b-axis are perpendicular to zaxis. Note, that the c-axis might be inclined towards the zaxis after this step.

    anchoring

    Define anchor in relative z-coordinates of the layered cell. At this position the 3d-solid is cut (vaccum being inserted). Load the =.in... files into xfplo to orient yourself and to find a proper anchor.

  • slab

    cutting layers

    After anchoring the third (z) coordinate is in absolute units. Define cutlayersat as the lowest and highest z-coordinate of the slab. Everything outside will be removed. If the interval is inverted cutting will not be applied.

    cutting atoms

    Supply a list of atoms (sites) in cutatoms to be removed from the result of cutting layers.

(see Structure Manipulation Algorithm)

Note

All data members can be set by simple assignement. A returned data member is returned as copy:

s=sla.Slabify()
a=s.zaxis # a copy of s.zaxis
s.zaxis=[1,1,0] # set s.zaxis

Note

Do not use =.in_... files from the Slabify output in context of the xfplo fermi surfaces.

hamdataCell()
Returns:cell – The individual vectors are the columns of the returned matrix.
Return type:3x3 numpy.ndarray

Return the primitive unit cell vectors of the underlying data (see prepare).

hamdataCCell()
Returns:cell – The individual vectors are the columns of the returned matrix.
Return type:3x3 numpy.ndarray

Return the conventional unit cell vectors of the underlying data (see prepare). For rhombohedral lattices the conventional lattice is the trigonal (hexagonal) cell if the spacegroup setting is hexagonal or the same as the primitive rhombohedral cell if the spacegroup setting is rhombohedral.

hamdataRCell()
Returns:cell – The individual vectors are the columns of the returned matrix.
Return type:3x3 numpy.ndarray

Return the primitive reciprocal unit cell vectors of the underlying data (see prepare). The vectors are in units of kscale.

For better understanding the following code does the same thing:

# assume s is a Slabify instance, LA is numpy.linalg, np is numpy
A=s.hamdataCell()
# Next the reciprocal cell is the inverse-transpose of A times 2*Pi
# Use internal scale though
G=  LA.inv(A.T)  *(2*np.pi)  / s.kscale
# now G is equal to s.hamdataRCell()
layerCell()
Returns:cell – The individual vectors are the columns of the returned matrix.
Return type:3x3 numpy.ndarray

Return the primitive unit cell vectors of the primary layer.

layerSites()
Returns:sites
Return type:a list of Sites

Return a copy of the list of Sites of the primary layer.

prepare(hamdatafilename)
Parameters:hamdatafilename (str) – name of (WF) Hamiltonian data file (usually +hamdata)

Setup structure and map hopping data.

printStructureSettings()

Print a summary of all structure settings.

calculateBandStructure(bandplot, suffix='')
Parameters:

Calculate the bandstructure for 'slab' and '3d' along the path defined by bandplot and produce corresponding files with a file suffix appended.

calculateBulkProjectedEDC(bandplot, energymesh, zaxis=[0, 0, 1], nz=30, kzs=None, suffix='', query=False)
Parameters:
  • bandplot (BandPlot) – see help of pyfplo.common.BandPlot
  • energymesh (EnergyContour) – see help for EnergyContour
  • zaxis (3-vector of float) – a vector in cartesian coordinates, which describes the projection direction along which the spectral density is k-integrated. This direction must be parallel to a reciprocal lattice vector, since we need periodicity in the projection direction in order to define a proper integration interval.
  • nz (int) – the number of integration intervals
  • kzs (sequence of float) – specify this instead of nz to build your own non-uniform integration mesh. You could for instance have a fine mesh in some region and a course one in the rest of the integration interval.
  • suffix (str) – a suffix to alter the file name.
  • query (bool) – if this is True the function returns a float without calculating much, indicating the length of the integration interval in units of kscale. Use this to build your own non-uniform integration net. Call calculateBulkProjectedEDC with kzs set to a list or numpy ndarray of kz-values in the desired interval. You can do what you want. The queried interval length just indicates the periodicity in the zaxis direction in said units.
Returns:

kzlength – The length of the integration interval along the projection direction

Return type:

float

Calculate bulk projected band energy distribution curves. This is the 1d-integral of the spectral density over a single period in k-space along a direction, indicated by zaxis. The bandplot input defines a set of high-symmetry points, which should be in a 2d-projection plane perpendicular to the zaxis. (The latter restriction is not really needed.)

The integration interval is from 0 to some number, which is determined by inspecting the zaxis and the lattice. If a user defined mesh is given in kzs the integration goes from k+zaxis*kzs[0] to k+zaxis*kzs[len(kzs)-1], where k is a k-point along the path spanned by the high symmetry points in bandplot.

calculateBulkProjectedFS(fso, zaxis=[0, 0, 1], nz=30, kzs=None, suffix='', query=False)
Parameters:
  • fso (FermiSurfaceOptions) – see help for FermiSurfaceOptions
  • zaxis (3-vector of float) – a vector in cartesian coordinates, which describes the projection direction along which the spectral density is k-integrated. This direction must be parallel to a reciprocal lattice vector, since we need periodicity in the projection direction in order to define a proper integration interval.
  • nz (int) – the number of integration intervals
  • kzs (sequence of float) – specify this instead of nz to build your own non-uniform integration mesh. You could for instance have a fine mesh in some region and a course one in the rest of the integration interval.
  • suffix (str) – a suffix to alter the file name.
  • query (bool) – if this is True the function returns without calculating much with a float indicating the length of the integration interval in units of kscale. Use this to build your own non-uniform integration net. Call calculateBulkProjectedFS with kzs set to a list or numpy ndarray of kz-values in the desired interval. You can do what you want. The queried interval length just indicates the periodicity in the zaxis direction in said units.
Returns:

kzlength – The length of the integration interval along the projection direction

Return type:

float

Calculate bulk projected Fermi surface projections. This is the 1d-integral of the spectral density over a single period in k-space along a direction, indicated by zaxis. The fso input defines a 2d mesh in a plane perpendicular to the projection axis, a Fermi energy and an imagnary energy part. The integration interval is from 0 to some number, which is determined by inspecting the zaxis and the lattice. If a user defined mesh is given in kzs the integration goes from k+zaxis*kzs[0] to k+zaxis*kzs[len(kzs)-1], where k is a k-point in the 2d mesh (from fso).

calculateFermiSurfaceCuts(fso, wds=None, bandplot=None, suffix='', forcerecalculation=False)
Parameters:

Construct Fermi surface cuts. The cuts are calculated in two steps.

diagonalization
For each point of the grid defined by FermiSurfaceOptions the Hamiltonian is diagonalized. The results are written to files +cut_band_sf and +cut_bweights_sf. This step is costly.
iso line determination
In this step the files mentioned above are read back in and from it the iso-lines corresponding to the current Fermi energy (FermiSurfaceOptions.fermienergy) are determined. The results are written to +cuts_spin1 (and +cuts_spin2), suffix is appended to these files basename. One can change the fermi energy or the weights definitions and re-run the iso line step without diagonalization by just running calculateFermiSurfaceCuts again with forcerecalculation=False
calculateEDC(bandplot, energymesh, penetrationdepth=-1.0, greenoptions=None, suffix='')
Parameters:
  • bandplot (BandPlot) – see help for pyfplo.common.BandPlot
  • energymesh (EnergyContour) – see help for EnergyContour
  • penetrationdepth (float) –

    define to which depth the spectral density is collected. penetrationdepth is measured from the atom closest to the vacuum. When interpreting output messages: orbital indices are increasing with the z-coordinate of the atom:

    Positive values mean depth in +hamdata length units (usually Bohr radii).

    Negative values mean depth in number of blocks. A block is the cell in =.in_final_PLlayer

  • greenoptions (GreenOptions) – see help for GreenOptions
  • suffix (str) – the output file suffix

Calculate energy distribution curves (EDC) along the path defined by bandplot. These are k,energy resolved surface spectral densities.

calculateFermiSurfaceSpectralDensity(fso, penetrationdepth=-1.0, greenoptions=None, suffix='')
Parameters:
  • fso (FermiSurfaceOptions) – see help for FermiSurfaceOptions
  • penetrationdepth (float) –

    define to which depth the spectral density is collected. penetrationdepth is measured from the atom closest to the vacuum. When interpreting output messages: orbital indices are increasing with the z-coordinate of the atom.

    Positive values mean depth in +hamdata length units (usually Bohr radii).

    Negative values mean depth in number of blocks. A block is the cell in =.in_final_PLlayer

  • greenoptions (GreenOptions) – see help for GreenOptions
  • suffix (str) – the output file suffix

Calculate the surface Fermi surface (k,k resolved spectral density) for a particular Fermi energy.

orbitalIndicesByDepth(lowerdepthdatalimit=1e+30, upperdepthdatalimit=1e+30)
Returns:list of orbital indices
Return type:numpy.ndarray

Return a list of orbital (WF) indices for orbitals which belong to layers of maximum depth lowerdepthdatalimit and upperdepthdatalimit measured from the lower (upper) end of the slab.

orbitalNames(orbitalindices=None)
Parameters:orbitalindices (sequence of int) – a sequence of valid orbital indices as e.g. returned by orbitalIndicesByDepth.
Returns:list of orbital names
Return type:list

Return a list of orbitalnames. Optionally, give a list of orbital indices to narrow down the returned list.

Example

import pyfplo.slabify as sla

s=sla.Slabify()
...
s.prepare('+hamdata')

# This is a list of all orbitals upto a depth of 5 Bohr radii
# at the upper side of a slab.
upper=s.orbitalNames(s.orbitalIndicesByDepth(-1,5))
# or simply
upper=s.orbitalNamesByDepth(-1,5)

# This is a list of all orbitals upto a depth of 5 Bohr radii
#at the lower side of a slab.
lower=s.orbitalNames(s.orbitalIndicesByDepth(5,-1))
#
# Here comes the rest of the orbitals.
# (A nice python trick to get the rest list)
rest=list(set(s.orbitalNames())-set(upper)-set(lower))
# this can also be done with index lists
iupper=s.orbitalIndicesByDepth(-1,5)
ilower=s.orbitalIndicesByDepth(10,-1)
iall=s.orbitalIndicesByDepth()
irest=list(set(iall)-set(iupper)-set(ilower))
print 'the orbital indices lists:',ilower,irest,iupper
upper=s.orbitalNames(iupper)
lower=s.orbitalNames(ilower)
rest=s.orbitalNames(irest)
print 'the orbital lists:',lower,rest,upper
orbitalNamesByDepth(lowerdepthdatalimit=1e+30, upperdepthdatalimit=1e+30)
Returns:list of orbital names
Return type:list

Convenience function which is equivalent to:

orbitalNames(orbitalIndicesByDepth(lowerdepthdatalimit,lowerdepthdatalimit))
orbitalIndicesBySite(isite)
Parameters:isite (int) – the site number
Returns:list of orbital indices
Return type:numpy.ndarray

Return a list of orbital (WF) indices for orbitals which belong to site number isite.

wannierCenterMatrix()
Returns:a list of 3 matrices each containing the cartesian component of the Wannier centers (site vectors) in the diagonal.
Return type:list of 3 numpy.ndarray

If \(\vec{s}\) is the site in the unit cell on which the Wannier function \(w_{\vec{R}\vec{s}n}\) sits, wannierCenterMatrix returns \(\vec{M}_{\vec{s}^\prime{}n^\prime,\vec{s}n} =\delta_{\vec{s}^\prime,\vec{s}}\delta_{n^\prime,n}\vec{s}\) where \(n^\prime\), \(n\) are orbital indices.

calculateBerryCurvatureOnBox(box, homo, suffix='', fullF=False, toldegen=1e-10)
Parameters:
  • box (BoxMesh) – see help for BoxMesh
  • homo (int) – The berry curvature is calculated for a set of bands. This sets the band number of the highest band included in this set.
  • suffix (str) – the output file suffix
  • fullF (int) – if the position operator (basis connection) is contained in +hamdata the full Berry curvature is calculated if fullF is True
  • toldegen (float) – levels closer than this tolerance are considered to be a degenerate subspace in the non-Abelian correction to the Berry curvature

Output of the Berry curvature on a box mesh (optionally integrated over the third box direction). The non-Abelian expression is used for degenerate subspaces and a symmetry-restoring correction in periodic gauge is applied. The curvature corresponds to \(A^k=-i\left\langle u\mid \nabla u\right\rangle\).

Now, a correction for the non-Abelian Berry curvature is applied for degenerate subspaces. Which levels are considered degenerate is controlled by toldegen. If it is too small, spurious results can be obtained. If it is very large, some of the actual curvature will be mising. toldegen should be smaller than the smallest physical gap in the band structure.

See berryCurvature for details about the Berry curvature calculation.

calculateChernNumberInSphere(center, radius=0.1, nsubdiv=20, homo=1, suffix='', nradius=1, fullF=False, toldegen=1e-10)
Parameters:
  • center (3-vector of float) – the sphere center
  • radius (float) – the sphere radius
  • nsubdiv (int) – the mesh subdivision level
  • homo (int) – The berry curvature is calculated for a set of bands. This sets the band number of the highest band included in this set
  • suffix (str) – the output file suffix
  • nradius (int) – if given and >1 include a radial-mesh in the output of the Berry curvature. The result is volumetric instead of 2d-on-sphere data.
  • fullF (int) – if the position operator (basis connection) is contained in +hamdata the full Berry curvature is calculated if fullF is True
  • toldegen (float) – levels closer than this tolerance are considered to be a degenerate subspace in the non-Abelian correction to the Berry curvature

Calculate the Chern number within a sphere. If you suspect the presence of a Weyl point first it is good to find its position. Then a small sphere can be put around the Weyl point position. The integral of the Berry curvature over this sphere gives the associates Chern number (Chirality). The output shows the Chern numbers for the set of the lowest n bands if it is larger than some tolerance, where n runs over all bands. This means that it is not the band wise Chern number but the Chern number for various sets of occupied bands with corresponding homo. If the shown number is far from integer try a finer nsubdiv. The number will converge to an integer eventually. Usually these spheres have to be rather samll, especially if the Berry curvature is very structured. Note, that center and radius are in units of kscale.

Now, a correction for the non-Abelian Berry curvature is applied for degenerate subspaces. Which levels are considered degenerate is controlled by toldegen. If it is too small, spurious results can be obtained. If it is very large, some of the actual curvature will be mising. toldegen should be smaller than the smallest physical gap in the band structure.

See berryCurvature for details about the Berry curvature calculation.

calculateZ2Invariant(gamma0, gamma1, gamma2, Nint, Nky, homos, efhomo=0, xmesh=None, ymesh=None, suffix='', gauge='periodic')
Parameters:
  • gamma0 (3-vector of float) – marks the center of the 2D plane in units of kscale
  • gamma1 (3-vector of float) – in units of kscale The integration direction is \(\Gamma_0 \to \Gamma_1\)
  • gamma2 (3-vector of float) – in units of kscale The ky-parameter of the Wannier centers runs along \(\Gamma_0 \to \Gamma_2\).
  • Nint (int) – gives the subdivision for the integration direction
  • Nky (int) – the subdivision of the ky-parameter direction.
  • homos (list or single int) – the band indices, which form the highest occuied bands.
  • efhomo (int) – One can specify the band below the Fermi energy as efhomo which is only used in output for orientation.
  • xmesh (sequence of float) – if given it must contain a grid in [-1,1], which will be used as integration grid instead of the default equidistant one. If given Nint is ignored.
  • ymesh (sequence of float) – if given it must contain a grid in [0,1], which will be used as ky-paramter grid instead of the default equidistant one. If given Nky is ignored. If the first point equals 0 and/or the last equals 1, the points are slightly shifted inwards to avoid particular problem cases.
  • suffix (str) – a convenience suffix to be appended to the created files.
  • gauge (str) – Can be 'periodic' (default) or 'relative'. The relative gauge seems to preserve the symmetry of the Wannier center curves, while the periodic does not. In pyfplo version <= 18.00 the periodic gauge was implemented. Formally, the relative gauge should be correct. This option was not tested a lot. The topological invariants should NOT depend on the gauge choice though, unless there is no gap. It seems that the Wannier center curves spread out more evenly in relative gauge, which makes the automatic index determination less reliable. So, visual checks should always be performed.

Calculate the Z2 invariant for a plane spanned by the TRIM points \(\Gamma_0\), \(\Gamma_1\) and \(\Gamma_2\) via Wannier centers.

All TRIM points must be actual points in the BZ (not directions). The plane spanned by the TRIMS shall only contain \(\Gamma_0\) in the center and the other TRIMS at the corners and mide-edges. No other TRIMS shall fall inside the planar cell. It is neccesary that the whole plane is gapped for the results to make sense. The TRIM points are always given with respect to the default 3d cell as defined by the data in +hamdata. In other words if you use enlarge the TRIM points must not be given with respect to the enlarged cell! See calculate3dTIInvariants for more explanations.

Example with xmesh/ymesh:

import numpy as np
...
s=sla.Slabify()
...
Gp=s.hamdataRCell()
G=[0,0,0]
Z=Gp.dot(np.array([0,0,1])/2.)
X=Gp.dot(np.array([1,0,0])/2.)
...
Nint=20
Nky=200
xmesh=map(lambda x: np.sign(x)*x**2,np.linspace(-1,1,Nint))
if False: # or
   xmesh=map(lambda x: x**3,np.linspace(-1,1,Nint))

ymesh=np.append(np.linspace(0,0.2,int(0.7*Nky),endpoint=False),
                np.linspace(0.2,1.,int(0.3*Nky)))
if False: # or
    ymesh=map(lambda x: np.sign(x)*x**2,np.linspace(0,1,Nky))

s.calculateZ2Invariant(G,X,Z,Nint,Nky,[14,16,18,20],xmesh,ymesh)
calculate3dTIInvariants(Nint, Nky, homos, efhomo=0, gauge='periodic')
Parameters:
  • Nint (int) – the number of intervals along the integration direction.
  • Nky (int) – the number of intervals along the ky-parameter direction.
  • homos (list or single int) – a list or a single homo can be given. A homo is the number of the highest band, assumed to be occupied. To find out the band numbers (homos) first make a 3d band calculation using calculateBandStructure and load the resulting +band... into xfbp: right click on a desired band and the set numbers and band numbers of the bands nearby are displayed. The band number is the important one, not the set number!
  • efhomo (int) – One can specify the band below the Fermi energy as efhomo which is only used in output for orientation.
  • gauge (str) – Can be 'periodic' (default) or 'relative'. The relative gauge seems to preserve the symmetry of the Wannier center curves, while the periodic does not. In pyfplo version <= 18.00 the periodic gauge was implemented. Formally, the relative gauge should be correct. This option was not tested a lot. The topological invariants should NOT depend on the gauge choice though, unless there is no gap. It seems that the Wannier center curves spread out more evenly in relative gauge, which makes the automatic index determination less reliable. So, visual checks should always be performed.

Use Wannier centers to calculate the Z2 invariants for a 3d TI.

For the whole thing to make sense the bulk band structure must be gapped above the band numbered by the homos throughout the whole BZ. The output will be a set of files, which have to be inspected by the user to decide how many Wannier centers cross a chosen horizontal reference line. There is also an automatic algorithm to determine the invariants (printed to output) This is, however, not always correct: a sufficiently large number of integration (Nint) and parameter points (Nky) are needed to obtain a valid result. Especially, if the wannier center curves vary fast in some parameter regions Nky must be sufficiently high for the automatic algorithm for the determination of the Z2 invariants to be correct. The reason are many small gaps and/or highly fluctuating Berry curvature.

The programm creates data files +Z2_homo..._... containing the Wannier centers where the last suffix indicates in which plane we are: _z0 is a (1/2 1/2 0) plane through the origin in primitive reciprocal basis, while _x1, _y1 and _z1 denote (0 1/2 1/2), (1/2 0 1/2) and (0 1/2 1/2) planes through (1/2 0 0), (0 1/2 0) and (0 0 1/2) as also printed to the output. If the invariants for z0 and z1 differ it is a strong TI. The x1, y1, z1 invariants give the three weak indices \(\nu_{1,2,3}\). Note, that the weak indices depend on the chosen planes.

Additionally files +zgap_homo_..._ containing a reference line which follows the largest gap [A. A. Soluyanov, PRB 83, 235401 (2011)] are created together with convenience files Z2_3dTI_homo*.cmd which load all data into xfbp:

xfbp Z2_3dTI_homo16.cmd

The reference line is printed with blue weights if the number of centers crossed so far is even and in red weights if the number is odd. If the last data point is odd the invariant is non-trivial. This algorithm only works for a sufficiently large grid although it is much better for smaller number of points when humans might not yet see how the centers are connected. Since this algo does not always work as wished we modified it such that some of the larger gaps are followed in separate curves (only one is shown in Z2_3dTI_homo...cmd). At the end we call the invariant odd if a majority of these gap-following curves indicate odd-ness. The reliability of the results are printed in the output table. Afterall, it is always better to check how the results converge with the grid spacing.

hamAtKPoint(kpoint, ms, gauge='relative', opindices=None, makedhk=False, makesigma=False, makexcfield=False, makebasisconnection=False, makewfsymops=False)
Parameters:
  • kpoint (3-vector of float) – The k-point in absolute cartesian coordinates.
  • ms (int) – The spin component. Full relaticistic: ms=0 else: ms in [0,nspin-1]
  • gauge (str) – Can be 'relative' (default) or 'periodic' or 'forcerelative'. There are two possible phase choices for the underlying Bloch sums. The default is the relative-distance gauge \(H_{s^{\prime}s}^{k}=\sum_{R}\mathrm{e}^{ik\left(R+s-s^{\prime}\right)}\left\langle w_{0s^{\prime}}\mid H\mid w_{Rs}\right\rangle\) The second is the periodic gauge which is needed if derivatives with respect to k are required: \(H_{s^{\prime}s}^{k}=\sum_{R}\mathrm{e}^{ikR}\left\langle w_{0s^{\prime}}\mid H\mid w_{Rs}\right\rangle\) If the velocity matrix is requested (makedhk=True) the gauge will be set to 'periodic' automatically unless the gauge is set to 'forcerelative' (this is new and allows to overwrite the old behaviour).
  • opindices (sequence of int) – A list of operation indices can be given to restrict the list of returned WFSymOp. See also makewfsymops. These indices are unique identifiers independend of the structure and are printed when prepare-ing the structure.
  • makedhk (int) – Return dH/dk in the wannier basis additionally to H. The return value will be a tuple (H,dhk,...) if this option is True. The actual type of dhk is BerryCurvatureData.
  • makesigma (int) – Return the 3 matrices \(\langle \sigma_{x,y,z}\rangle\) in the Wannier basis additionally to H. The return value will be a tuple (H,…,sigma,…) if this option is True.
  • makexcfield (int) – Return the 3 matrices \(\mu_B\langle B^{\mathrm{xc}}_{x,y,z}\rangle\) in the Wannier basis in eV additionally to H. The return value will be a tuple (H,…,xcfield,…) if this option is True. (Only for full relativistic.) In a spin polarized calculation `` H- ( B[0].dot(S[0])+B[1].dot(S[1])+B[2].dot(S[2]))`` will approximately be the non-polarized Hamiltonian, where H, S and B are returned by hamAtKPoint with options makesigma and makexcfield. This is approximate since to obtain B we need to factor the xc-Zeeman term into the spin matrices and the field matrices. The FPLO basis is not complete and the Wannier basis even less. An illustrative example is given in ./Examples/slabify/Fe/SP/slabify/xcfield.
  • makebasisconnection (int) – Return, additionally to H, the 6 matrices \(\vec{A}^{w,k}=-i\left\langle u_w\mid \nabla u_w\right\rangle\) and \(\nabla\times{}\vec{A}^{w,k}\) in the Wannier basis, which are the Berry connection and Berry curvature of the WF basis itself: The relation to the position operator is given in periodic gauge by \(-\vec{A}^{b,k}_{mn}=\langle \vec{r}\rangle^k_{mn}\) and in relative gauge by \(-\vec{A}^k_{mn}=\langle \vec{r}\rangle^k- \delta_{mn}\delta_{ts}\vec{s}\) where \(s\) and \(t\) are Wannier centers (sites) and \(\langle \vec{r}\rangle^k =\sum_R \mathrm{e}^{ik\left(R+\lambda(s-t)\right)} \langle w_{0t} \mid r\mid w_{Rs} \rangle\) where \(\lambda=0\) in the periodic gauge and \(\lambda=1\) in the relative gauge. The return value will be a tuple (H,...,basisconnection,...) if this option is True. The first 3 components of basisconnection are the basis Berry connection and the last 3 the basis Berry curvature. To convert the connection in relative gauge to the position operator add the matrices returned by wannierCenterMatrix. See argument gauge above to remember why gauge=forcerelative needs to be used when relative is required. Examples for the use of this argument are found in ./Examples/slabify/Fe/SP/slabify/3dR and ./Examples/slabify/Fe/SP/slabify/AHC.
  • makewfsymops (int) – If True return a list of WFSymOp instances, which describe the symmetry transformations of the Hamiltonian and the Wannier function Bloch sums. The return value will be a tuple (H,…,wfsymops) if this option is True.
Returns:

(H,…) – If all make… options are False it simply returns the Hamiltonian H (no tuple) at kpoint for spin ms. If any of the make… options are True a tuple is returned containing H and all additionally requested operators in the order defined by the make…-arguments in the argument list above. The order of the keyword arguments in an actual function call does not matter for the order of return values. If the requested operator is a vector operator it is returned as a list of numpy.ndarrays. The derivative dH/dk is in units of eV*aB.

Return type:

numpy.ndarray or a tuple of return values

Example

For convenience you can use kscale as in:

s=sla.Slabify()
k=[1,0,0]
ms=0
H=hamAtKPoint(k*s.kscale,ms) # or to get dH/dk
(H,dHk)=hamAtKPoint(k*s.kscale,ms,makedhk=True)
diagonalize(h, dhk=None, makef=False, basisconnection=None, toldegen=1e-10)
Parameters:
  • h (square complex matrix) – the Hamiltonian
  • dhk (BerryCurvatureData) – behaves as a list`of 3 complex `numpy.ndarrays which contain dH/dk as returned by hamAtKPoint
  • makef (int) – if True return band wise Berry curvature
  • basisconnection (list of 6 complex numpy.ndarrays) – basisconnection and curvature as returned by hamAtKPoint
  • toldegen (float) – levels closer than this tolerance are considered to be a degenerate subspace in the non-Abelian correction to the Berry curvature
Returns:

(E,C,…) – The first tuple value is a numpy.ndarray of the energies, the second a numpy.ndarray of eigenvectors C[:,i] is the i-th eigenvector. If makef is True the third tuple elements is a numpy.ndarray containing the Berry curvature where F[:,n] is the band-n Berry curvature (a 3-vector). The dimension of the Hamiltonian is nvdim.

Return type:

tuple of return values

Diagonalize the Hamiltonian h and return the eigenvalues and eigenvectors. If makef is True, dhk must be given. Then additionally the Berry curvatures for all bands is returned. If basisconnection is also supplied the full Berry curvature is returned. See berryCurvature.

dhk must be obtained from hamAtKPoint. New in version 19: Note, that in the periodic gauge (used for h and dhk) a correction for the position operator is used when calculating the Berry curvature. If gauge=forcerelative was used in hamAtKPoint this correction is zero. The curvature corresponds to \(A^k=-i\left\langle u\mid \nabla u\right\rangle\).

Now, a correction for the non-Abelian Berry curvature is applied for degenerate subspaces. Which levels are considered degenerate is controlled by toldegen. If it is too small, spurious results can be obtained. If it is very large, some of the actual curvature will be mising. toldegen should be smaller than the smallest physical gap in the band structure.

Example:

k=np.array([0.5,0,0])
(Hk,dHk)=s.hamAtKPoint(k*s.kscale,ms,makedhk=True)

(E,CC,F)=s.diagonalize(Hk,dhk=dHk,makef=True,
                       toldegen=1.0e-9)
print 'Berry curvature of homo',homo,': ',F[:,homo]
diagonalizeUnitary(U, evdegentol=1e-08)
Parameters:
  • U (square complex matrix) – the unitary matrix
  • evdegentol (float) – eigenvalues with a distance less than this are considered degenerate.
Returns:

(E,Z)E are the complex eigenvalues and Z the eigenvectors.

Return type:

tuple of return values

For unitary matrices LAPACK does not return orthogonal eigenvectors for degenerate subspaces. This method calls LAPACK for diagonalization and corrects the eigenvectosr afterwards.

coDiagonalize(E, C, Dk, evtol=1e-08, check=False)
Parameters:
  • E (sequence of float) – the Hamiltonian eigenvalues
  • C (square complex matrix) – the Hamiltonian eigenvectors
  • Dk (square complex matrix) – the Bloch sum symmetry representation matrix
  • evtol (float) – tolerance to determine degenerate subspaces of E
  • check (int) – if True perform paranoia checks
Returns:

(EU,CZ)EU are the complex symmetry eigenvalues and CZ the transformed eigenvectors which diagonalize the Hamiltonian and the symmetry.

Return type:

tuple of return values

Given Hamiltonian eigenvalues E and eigenvectors C such that \(H C=C E\) and a representation matrix Dk for the Wannier function Bloch sums with symmetry properties \(H=D^{+} H D\) determine \(U = C^{+} D C\).

Then diagonalize \(U\): \(U=Z E_{U} Z^{+}\) in each degenerate subspace of E and form the transformed eigenvectors \(C^{\prime}=C Z\). Now \(C^{\prime}\) diagonalizes \(H\) and \(U\). The eigenvalues \(E_{U}\) of \(U\) and the transformed eigenvectors \(C^{\prime}\) are returned.

evtol is used to determine which Hamiltonian eigenvalues E[i] are degenerate.

The symmetry representation matrices are return by hamAtKPoint.

Note, that operations which are combined with time reversal do not really make sense in this context, so do not use this routine if the operation is combined time reversal.

berryCurvature(E, C, dhk, subspace=None, basisconnection=None, toldegen=1e-10, returndetails=False)
Parameters:
  • E (sequence of float) – the Hamiltonian eigenvalues
  • C (square complex matrix) – the Hamiltonian eigenvectors
  • dhk (BerryCurvatureData) – BerryCurvatureData, behaves as a list`of 3 complex `numpy.ndarrays which containdH/dk as returned by hamAtKPoint
  • subspace (sequence of int) – a list of zeros and ones, which describes the wanted subspace or None for the whole space.
  • basisconnection (list of 6 complex numpy.ndarrays) – basisconnection and curvature as returned by hamAtKPoint
  • toldegen (float) – levels closer than this tolerance are considered to be a degenerate subspace in the non-Abelian correction to the Berry curvature
  • returndetails (bool) – if True details of the curvature contributions are returned
Returns:

(F,divergenttermsdetected)F is a numpy.ndarray of shape=(3,nvdim) and contains the band wise Berry curvature. divergenttermsdetected is now always False. It is still there to not break older code. Sorry. If basisconnection is given and returndetails==True the returned tuple is (F,Fdetails,divergenttermsdetected), where Fdetails is a dict of the various terms F is made of. The dict values are again numpy.ndarray of the same shape as F.

Return type:

tuple

Calculate the Berry curvature from the Hamiltonian eigenvalues E and eigenvectors C and form dH/dk returned from hamAtKPoint. If basisconnection is given (as returned from hamAtKPoint) the full curvature is calculated (Wang2006) not only the D-D part. If returndetails==True the separate terms F is made off are returned in a dict.

If subspace is None the Berry curvature is returned for all bands. If it is a list of size len(E) the Berry curvature is calculated from the subspace for which subspace contains a nonzero value, preferably 1. This can be used to make plots of the k-resolved mirror Chern number.

If basisconnection is not provided, the resulting Berry curvature will contain a correction for the position operator in the periodic gauge. In the relative gauge (gauge=forcerelative) this correction is zero.

The curvature corresponds to \(A^k=-i\left\langle u\mid \nabla u\right\rangle\).

Now, a correction for the non-Abelian Berry curvature is applied for degenerate subspaces. Which levels are considered degenerate is controlled by toldegen. If it is too small, spurious results can be obtained. If it is very large, some of the actual curvature will be mising. toldegen should be smaller than the smallest physical gap in the band structure.

FPLO/...DOC/pyfplo/Examples/slabify/TCI/SnTe contains an example.

findWeylPoints(box, homos, tol=0.001)
Parameters:
  • box (BoxMesh) – a definition of the box and its mesh.
  • homos (list or single int) – the band indices, which form the highest occuied bands. If a list is given WPs for all these highest occuied bands are serached.
  • tol (float) – This sets the finest bi-section of the refinement part of the algorithm. What is a usefull value depends on the actual structuring of the Berry curvature.

This method tries to automatically find Weyl points using a modification of the algorithm described in [Takahiro Fukui et.al. J. Phys. Soc. Jpn. 74, 1674 (2005)] The user defines a box which spans a grid. On all grid points the Hamiltonian is diagonalized. For each micro cell of the grid the Berry curvature is integrated over the surface of the cell. This gives the chirality of the WP if the box contains a WP. For all cells, which have a non zero result subsequent bisections is performed until the box size falls below tol. Finally, all resulting boxes with non-zero chirality are written to the file weylpoints.py. This file can be loaded into a script for further processing (e.g. to use calculateChernNumberInSphere to make sure that is actually is a Weyl point.)

Note: that the alogrithm does not always find all Weyl points due to several reasons.

  • The box grid must be fine enough such that the Berry curvature variation is properly sampled by the micro cells.
  • The origin of the box matters if e.g. the WPs are pinned to symmetry planes The WPs ideally sit well within a micro cell. Since, one can not always know this in advance it is probably good to try two different origins. One at say (000) and another at \(-\Delta_k/2\), where \(\Delta_k\) is the body diagonal of a micro cell.
  • A micro cell might contain two WPs of opposite chirality, which would result in zero total chirality and hence these WPs are missed. The grid must be fine enough to avoid this.

Often symmetries can be used to complete the list of found WPs.

Warning: the algorithm can produce false positives. This is why it is always a good idea to check the validity of the results via calculateChernNumberInSphere.

(also see the Weyl semi metal example)

calculateMirrorChernNumbers(homo, nint, atneqk=False, showevs=False, evtol=1e-08, fullF=False, toldegen=1e-10)
Parameters:
  • homo (int) – the highest occupied band
  • nint (int) – the number of subdivision in the first k direction of the mirror plane
  • atneqk (int) – if True also calculate for plane at nonzero out-off plane position
  • showevs (int) – if True show mirror eigenvalues
  • evtol (float) – energies eigenvalues which are closer than this tolerance are considered degenerate when determining mirror subspaces
  • fullF (int) – if the position operator (basis connection) is contained in +hamdata the full Berry curvature is calculated if fullF is True
  • toldegen (float) – levels closer than this tolerance are considered to be a degenerate subspace in the non-Abelian correction to the Berry curvature

For each symmorphic mirror operation the mirror chern number for the mirror plane in k-space is calculated. The result will depend on the density of the inregration mesh which is controled by nint. The mirror eigenvalues are determined in the subspaces of degenerate energies. This increases the accuracy of the diagonalization of the unitary transformation. The Hamiltonian eigenvectors are transformed to belong to mirror eigenvalue subspaces from which the Berry curvature is determined. toldegen controls the application of the non-Abelian correction to the Berry curvature (see berryCurvature).

The mirror planes are determined automatically, which leads to a 2d reciprocal basis in the plane and an out off plane vector an integer multiple of which describes the position of all mirror planes. The mirror chern number is calculated for the plane through the origin and for the closest plane.

See berryCurvature for details about the Berry curvature calculation.

FPLO/...DOC/pyfplo/Examples/slabify/TCI/SnTe contains an example.

dirname

Local result directory. You can also set it

s=sla.Slabify()
s.dirname='.' # now output lands in the current directory
Type:str
nvdim

Return dimension of WF Hamiltonian.

Type:int
nspin

Return number of spins (always 1 for full-relativistic).

Type:int
object

The type of structure-object to be created

Values: '3d', 'slab' or 'semislab'

Type:str
enlarge

3x3 list or numpy.ndarray: A 3x3 integer matrix \(U\) for enlarging the cell. The rows of \(U\) represent the three new vectors. This matrix produces a new unit cell \(A^\prime\) out of the old \(A\) via \(A^\prime=U\cdot A\) where we assume that \(A\) is a column vector of the tree lattice vectors. \(A=\left(\begin{array}{c} \boldsymbol{a}_{1} \\ \boldsymbol{a}_{2} \\ \boldsymbol{a}_{3} \end{array}\right)\) (see Structure Manipulation Algorithm)

zaxis

The axis perpendicular to the surface (in cartesian coordinates) (see Structure Manipulation Algorithm)

Type:3-vector
numberoflayers

Number of layers to be constructed from the layered cell. For 'semislab's the smallest possible number should be chosen. This number is determined by the maximal inter-atom overlap in the Hamiltonian data. It is best to start with 1 and then increase it according to the program output (it will complain if needed). (see Structure Manipulation Algorithm)

Type:int
anchor

at which relative z-position in the layered unit cell do we cut the solid.

Type:float
Type:for 'slab' and 'semislab'
cutlayersat

Currently only for 'slab'! Lower and upper absolute z-coordinate at which to cut a finite slab. If cutlayersat[1]<cutlayersat[0] it is not applied. (see Structure Manipulation Algorithm)

Type:list of 2 float
cutatoms

Currently only for 'slab'! A list of atoms to be removed in the final step. Note that the site numbers refer to the cell after cutlayersat was applied. Use xfplo on the resulting =.in file to find the site numbers. (see Structure Manipulation Algorithm)

Type:list of int
options

A set of options (for debugging), see pyfplo.common.OptionSet

Example

import pyfplo.slabify as sla
s=sla.Slabify()
op=s.options
for n in op.names:
   print n,'is set to ',op[n]
Type:OptionSet
kscale

The common factor used for the k-points. Usually k-points are given as 2*pi/a *(kx,ky,kz) where a is obtained from +hamdata or more precisely from the conventional cell. One can set kscale to any value which is convenient.

Note: For rhombohedral lattices the conventional lattice is the trigonal (hexagonal) cell if the spacegroup setting is hexagonal or the same as the primitive rhombohedral cell if the spacegroup setting is rhombohedral. Consequently kscale depends on the setting for rhombohedral lattices.

Type:float
bfield

the optional B-field configuration.

Type:BfieldConfig
hasxcfield

if True the Bxc-data are available in +hamdata

hassigma

if True \(\langle \sigma_{x,y,z}\rangle\) is available in +hamdata

hasbasisconnection

if True basis connection and curvature are available in +hamdata

BoxMesh

class BoxMesh

Define a box mesh in 3 dimensions. For a planar mesh set nz=1 and zinterval=[0,0].

setBox(xaxis=[1, 0, 0], yaxis=[0, 1, 0], zaxis=[0, 0, 1], origin=[0, 0, 0])
Returns:self – for call chaining.
Return type:BoxMesh

Convenience function: Define the box. See origin, xaxis, yaxis, zaxis for docs.

setMesh(nx=100, ny=100, nz=100, xinterval=[-0.5, 0.5], yinterval=[-0.5, 0.5], zinterval=[-0.5, 0.5])
Returns:self – for call chaining.
Return type:BoxMesh

Convenience function: Define the mesh subdivisions of the box. See nx, ny, nz, xinterval, yinterval and zinterval for docs.

xmesh()
Returns:xmesh – the mesh coordinates along the xaxis
Return type:numpy.ndarray
ymesh()
Returns:ymesh – the mesh coordinates along the yaxis
Return type:numpy.ndarray
zmesh()
Returns:zmesh – the mesh coordinates along the zaxis
Return type:numpy.ndarray
mesh(scale, xfirst=False)
Parameters:
  • scale (float) – all points are multiplied with scale
  • xfirst (int) – if True invert the loop order (make the x-loop the outer loop).
Returns:

mesh – The flat list of all vectors of the mesh.

Return type:

numpy.ndarray

Return the box mesh in absolute coordinates as a single flattend continous numpy.ndarray of vectors. This routine actually first calculates the mesh. Hence it is strongly advised to use a local copy as in:

kpnts=box.mesh(scale,False)
someFunctionCall(...,kpnts,...)

When creating the mesh the innermost loop is x, then y, the outer loop is z. If xfirst is True the loop order is inverted. All points are scaled by scale.

The mesh is calculated as:

mesh[ind]=(xaxis*xinterval[i]+
           yaxis*yinterval[j]+
           zaxis*zinterval[k]+origin)*scale`

where i, j and k run over all posible interval values and ind runs over all mesh indices depending on the loop order.

relToAbs(*args)
Parameters:args (3 float or list or numpy.ndarray) – can be a 3-vector or three coordinates.
Returns:point – the mesh point in absolute coordinates.
Return type:numpy.ndarray

Take a point (x, y, z) in relative mesh coordinates and return the absolute coordinates. Relative means relative with respect to the normalized axes not the intervals:

abs=origin+x*xaxis+y*yaxis+z*zaxis

absToRel(*args)
Parameters:args (3 float or list or numpy.ndarray) – can be a 3-vector or three coordinates.
Returns:point – the mesh point in relative coordinates.
Return type:numpy.ndarray

Take a point in absolute coordinates and return the relative mesh coordinates (x, y, z). Relative means relative with respect to the normalized axes not the intervals:

abs=origin+x*xaxis+y*yaxis+z*zaxis

xyzFromIndex(ind, xfirst=False)
Parameters:
  • ind (int) – an index into the continous flattend mesh array of vectors.
  • xfirst (int) – if True invert the loop order (make the x-loop the outer loop).
Returns:

xyz – a list of 3 indices such that the mesh position xyz[0], xyz[1], xyz[2] has point mesh(ind)

Return type:

numpy.ndarray of 3 int

Given an index into the flattened mesh array return the individual indices of the three box axis: ix, iy, iz. The order of loops matter (see mesh).

xfirst=False: ind=ix+nx*(iy+ny*iz)

xfirst=True : ind=iz+nz*(iy+ny*ix)

__str__()

return printable representation. You do not need to call this explicitly. An object obj with this function provides usefull info when printed:

print(obj)
origin

numpy.ndarray: The origin(shift) of the box mesh given with respect to the global coordinate system, in some units, which usually are pyfplo.slabify.Slabify.kscale. (See the scale parameter in mesh)

xinterval

numpy.ndarray: Interval along the first (x) axis of the box as defined by xaxis. The units often are pyfplo.slabify.Slabify.kscale, which makes it a uniform scale no matter the orientation of the box.

yinterval

numpy.ndarray: Interval along the first (y) axis of the box as defined by yaxis. The units often are pyfplo.slabify.Slabify.kscale, which makes it a uniform scale no matter the orientation of the box.

zinterval

numpy.ndarray: Interval along the first (z) axis of the box as defined by zaxis. The units often are pyfplo.slabify.Slabify.kscale, which makes it a uniform scale no matter the orientation of the box.

nx

number of points in the first (x) direction.

Type:int
ny

number of points in the second (y) direction.

Type:int
nz

number of points into the thid (z) direction, which also is integrated over for planar projections

Type:int
xaxis

the first (x) axis spanning the box. It will be automatically normalized.

Type:numpy.ndarray
yaxis

the second (y) axis spanning the box. It will be automatically normalized.

Type:numpy.ndarray
zaxis

the third (z) axis spanning the box. It will be automatically normalized. This is integrated over.

Type:numpy.ndarray

EnergyContour

class EnergyContour(ne=100, e0=-1, e1=1, ime='auto')

Create a complex energy contour for energy distribution curves (EDC) Currently, only equidistant and parallel to real axis. Usage:

import pyfplo.slabify as sla
ec=sla.EnergyContour(100,-1.5,2.5,0.001)

or set individual properties:

import pyfplo.slabify as sla
ec=sla.EnergyContour()
ec.ne=100
...

or set individual mesh:

import pyfplo.slabify as sla
ec=sla.EnergyContour()
ec.setMesh([1,2,3,4+.1j])
Parameters:
  • ne (int) – see ne
  • e0 (float) – see e0
  • e1 (float) – see e1
  • ime ('auto' or float) – see ime

Set an equidistant energy contour paralell to the real axis. See the individual properties for further docs.

mesh()
Returns:mesh – the energy mesh
Return type:numpy.ndarray
setMesh(m)
Parameters:m (list or numpy.ndarray of complex) – a list of energy points in the complex plane

Set a hand crafted mesh. DO NOT set any other variable after this.

__str__()

return printable representation. You do not need to call this explicitly. An object obj with this function provides usefull info when printed:

print(obj)
ne

number of energy points

Type:int
e0

real part of starting point

Type:float
e1

real part of end point

Type:float
ime

imaginary part of the contour. ime can be a real number or the string 'auto' in which case ime is set to (e1-e0)/ne Of course you have to define e0, e1 and ne first or simply use the constructor.

Type:float or str

FermiSurfaceOptions

class FermiSurfaceOptions

Fermisurface controls.

Define the 2d k-mesh on which to calculate.

For '3d'/'slab' this mesh is used as a basis for interpolation to obtain the iso-lines. For 'semislab' it defines the pixel mesh on which to calculate the spectral density.

Note, that the slabify structure manipulation routines transform the original cell in such a way that the orientation of the resulting cell with respect to the original global cartesian system is unaltered. The k-mesh, is given in this global cartesian system. So a slab with zaxis=[1,0,0] has a surface BZ in the 010/001-plane and the axis must be set accordingly. Always try with a small mesh first to get your bearings. Also for spectral densities (pixelized pictures) in semislab mode the k-plane axis must be orthogonal!

setPlane(xaxis=[1, 0, 0], yaxis=[0, 1, 0], origin=[0, 0, 0])
Parameters:
  • xaxis (3-vector of float) – see xaxis
  • yaxis (3-vector of float) – see yaxis
  • origin (3-vector of float) – see origin

Convenience function: Define the k-plane.

setMesh(nx=100, xinterval=[-0.5, 0.5], ny=100, yinterval=[-0.5, 0.5])
Parameters:
  • nx (int) – see nx
  • xinterval (2-vector of float) – see xinterval
  • ny (int) – see ny
  • yinterval (2-vector of float) – see yinterval

Convenience function: Define the mesh subdivisions of the k-plane.

on()

Activate fermi-surface related routines.

off()

Deactivate fermi-surface related routines.

xmesh()
Returns:xmesh – the mesh coordinates along the xaxis
Return type:numpy.ndarray
ymesh()
Returns:ymesh – the mesh coordinates along the yaxis
Return type:numpy.ndarray
mesh(scale)
Parameters:scale (float) – all points are multiplied with scale
Returns:mesh – The flat list of all vectors of the mesh. x is running first.
Return type:numpy.ndarray

Return the box mesh in absolute coordinates as a single flattend continous numpy.ndarray of vectors. This routine actually first calculates the mesh. Hence it is strongly advised to use a local copy as in:

kpnts=fso.mesh(scale)
someFunctionCall(...,kpnts,...)

When creating the mesh the innermost loop is y, the outer loop is x. All points are scaled by scale.

The mesh is calculated as:

mesh[ind]=(xaxis*xmesh()[i]+yaxis*ymesh()[j]+origin)*scale`

where i, j interval values and ind runs over all mesh indices.

openDensPlotFile(filename, ms, plotorigin=[0, 0], plotxaxis=[1, 0], plotyaxis=[0, 1], progress=None)
Parameters:
  • filename (str) – the name of the xynz-type file
  • ms (int) – The spin component. Full relaticistic: ms=0 else: ms in [0,nspin-1]
  • plotorigin (2-vector of float) – a vector describing the origin in the plotting plane
  • plotxaxis (2-vector of float) – a vector describing the x-axis in the plotting plane
  • plotyaxis (2-vector of float) – a vector describing the y-axis in the plotting plane
  • progress – a progress message (str) or None
Returns:

band file context

Return type:

DensPlotContext

Low level routine. Return an object of type DensPlotContext for creation of FPLO xynz density plot files.

The returned object will open the file and organizes the proper file format. Its DensPlotContext.write() method can be used to write the actual data. If the object gets deleted (automatic if the scope is left) the file gets closed. The DensPlotContext.close() method can be called explicitly.

The best way to use it is in a with-statement. Then it is closed automatically after the with-block is exited:

with fso.openDensPlotFile(...) as f:
  for ...:
    f.write(...)
pass # now the file is closed.

If multiple files are written at the same time one can do the following:

f1=fso.openDensPlotFile(filename1,...)
f2=fso.openDensPlotFile(filename2,...)
for ...:
  f1.write(data1,...)
  f2.write(data2,...)
f1.close()
f2.close()
pass # now the files are closed.

Usually FermiSurfaceOptions defines a plane in k-space on which pixeled data are calculated. To plot them, we need to map this onto the plotting x,y-plane. If the x,y-axes of FermiSurfaceOptions are not orthogonal one wants to give two axes fo the same angle as arguments to openDensPlotFile to orient the data in the plotting plane. The origin in the plotting plane can also be given. The length of plotxaxis and plotyaxis determines the length scale when plotted with xfbp.

If progress is set to a string a progress message is written in subsequent calls to DensPlotContext.write().

see help of DensPlotContext.

An example can be found in ./Examples/slabify/densplot and ./Examples/slabify/Fe/SP/slabify/AHC/curvature.py.

__str__()

return printable representation. You do not need to call this explicitly. An object obj with this function provides usefull info when printed:

print(obj)
active

if True calculate the 3d/slab Fermi surface features.

Type:bool
nx

number of points into the first (x) direction.

Type:int
ny

number of points into the first (y) direction.

Type:int
xinterval

interval along the first (x) axis of the BZ-plane as defined by the in-plane vector xaxis. The units are pyfplo.slabify.Slabify.kscale, which makes it a uniform scale no matter the orientation of the k-plane.

Type:numpy.ndarray
yinterval

interval along the first (y) axis of the BZ-plane as defined by the in-plane vector yaxis. The units are pyfplo.slabify.Slabify.kscale, which makes it a uniform scale no matter the orientation of the k-plane.

Type:numpy.ndarray
fermienergy

at which Fermi energy? Units like in +hamdata, which by default is eV.

Type:float
fermienergyim

a finite imaginary part for the energy is needed for. semislab spectral functions (same units as fermienergy). The smaller this value the higher the k-point density must be. Often a good value is interval-length/max(Nx,Ny)

Type:float
origin

the origin of the k-plane in units of pyfplo.slabify.Slabify.kscale

Type:numpy.ndarray
xaxis

the first (x) axis spanning the k-plane. It will be automatically normalized. For semislab the axes must be orthogonal

Type:3-vector of float
yaxis

the first (y) axis spanning the k-plane. It will be automatically normalized. For semislab the axes must be orthogonal

Type:3-vector of float

DensPlotContext

class DensPlotContext

This class wraps data to easily manage the creation of density plot files. This class cannot be instantiated directly. It only is produced and returned via a call to FermiSurfaceOptions.openDensPlotFile() An example can be found in ./Examples/slabify/densplot and ./Examples/slabify/Fe/SP/slabify/AHC/curvature.py.

close()

Explicitely close the file. Usefull, if multiple files are used in the same loop, in which case the with-statement is not usefull. The underlying file gets closed when this object gets garbage collected (after its scope is exited). For cleanliness it is a good measure to always close files.

Method1:

with fso.openDensPlotFile(...) as f:
   doseomthing with f
# here f is closed

Method2:

f=fso.openDensPlotFile(...):
do something with f
f.close()
# here f is closed
write(components)
Parameters:components (sequence (list,tuple,..)) – sequence of values for all density data components at this k-point. There can be more than one value per k-point, e.g. three components of a vector field.

Write the density data components for the current k-vector to the file.

GreenOptions

class GreenOptions(nsigiter=30, sigitertol=0.001, sigitermethod='accel')

Settings which control the self energy (sigma) calculation.

Parameters:

Create a GreenOptions object with the settings given in the argument.

__str__()

return printable representation. You do not need to call this explicitly. An object obj with this function provides usefull info when printed:

print(obj)
nsigiter

maximum number of sigma (self energy) iteration loops.

Type:int
sigitermethod

sigma iteration method, 'accel' (default) or ''

Type:str
sigitertol

sigma iteration tolerance.

Type:float

WeylPoint

class WeylPoint(k, axis1=[1.0, 0.0, 0.0], axis2=[0.0, 1.0, 0.0], axis3=[0.0, 0.0, 1.0], chirality=1.0, radius=0.1, energy=0.0, homo=1, spin=1)

Collect information about a Weyl point. This is for easier book keeping.

for documentation consult the individual properties.

__str__()

return printable representation. You do not need to call this explicitly. An object obj with this function provides usefull info when printed:

print(obj)
k

k is the WP position. Usually it is given in units Slabify.kscale

Type:3-vector
axis1

the first axis to span a volume around the WP or for bandplot purposes. On being set it will be normalized.

Type:3-vector
axis2

the second axis to span a volume around the WP or for bandplot purposes. On being set it will be normalized.

Type:3-vector
axis3

the third axis to span a volume around the WP or for bandplot purposes. On being set it will be normalized.

Type:3-vector
chirality

Store the chirality

Type:float
radius

a radius within which a monopole was found.

Type:float
energy

float at which energy does the WP sit.

homo

The highest occupid band.

Type:int
spin

always 1).

Type:int
Type:the spin (1=up, 2=down, relativistic

BfieldConfig

class BfieldConfig

The class allows to add model spin-only magnetic fields. to the Hamiltonian. This class cannot be instantiated directly. It only is returned from objects, which have an BfieldConfig member variable (see Slabify.bfield). Example usage:

s=sla.Slabify()
bf=s.bfield
bf.setGlobalField([0,0,1.3])
#or
s=sla.Slabify()
s.bfield.setGlobalField([0,0,1.3])

For an example see ./Examples/slabify/Fe/NSP/slabify/addBfield.

setGlobalField(B)
Parameters:B (3-vector of float) – the magnetic field

Set constant global spin-only magnetic field. The last field set in the script (global or local) will win.

setLocalFields(listofcomponents)
Parameters:listofcomponents (list of pairs of a projector (list) and field 3-vector) –

list of pairs of a projector and field vector:

[
 [ [0,0,0,1,1,1,1,0,0,0],[0, 1.2,0] ],
 [ [1,1,1,0,0,0,0,0,0,0],[0,-1.2,0] ],
...
]

where the first internal list is a projector P and the second a field B. The projector must have the dimension of the Hamiltonian. It usually is a list of zeros and ones. The Hamiltonian is modified according to:

\(H=H_0+ \sum_i P_i \left(\vec{\sigma}\cdot\vec{B}_i\right) P_i\)

which means that the blocks where P is nonzero get a Zeeman term added. There can be several P, B pairs as indicated by the subscript i. You can e.g. use Slabify.orbitalNames() and python map to create the projectors.

Set constant local spin-only magnetic fields. The last field set in the script (global or local) will win.

__str__()

return printable representation. You do not need to call this explicitly. An object obj with this function provides usefull info when printed:

print(obj)

WFSymOp

class WFSymOp

A list of instances of this class is returned by hamAtKPoint if options are set accodingly. It provides information on the symmetry transfrormations of the Hamiltonian and Wannier function Bloch sums.

There is a script in FPLO/...DOC/pyfplo/Examples/slabify/symmetryops which after adjusting the path should run through without errors. It demonstrates the symmetry operations and also prints the eigenvalues of operations. Use this as a starting point for your own work.

__str__()

return printable representation. You do not need to call this explicitly. An object obj with this function provides usefull info when printed:

print(obj)
symbol

A symbolic representation of the symmetry operation. If in doubt refer to alpha and tau for the exact meaning.

index

This is a unique index, which is used to identify particular operations. A list of such indices can be used in a call to Slabify.hamAtKPoint to request a subset of symmetries. The indices are unique identifiers independend of the structure and are printed when Slabify.prepare-ing the structure.

isinlittlegroup

This is True if the operation is in the little group of the k-point for which Slabify.hamAtKPoint was called.

alpha

This is the (improper) rotational part (3x3-matrix ) of the symmetry operation. It acts on a real space vector (or k-point) as \(\vec{r}^{\prime} = \alpha \vec{r}\). Together with tau it forms the seitz symbol of the space group operation which acts like \(\{\alpha\mid\tau\} \vec{r} = \alpha \vec{r}+\tau\). alpha and tau are given in absolute cartesian coordinates.

tau

This is the fractional translational part of the space group operation. It describes the location of the operation in space and possible glide and screw components. If it is zero the operation is definitely symmorphic but if it is non-zero it can be symmorphic (operation does not sit at (000) ) or non-symmorpic (is a glide/screw). See alpha.

timerev

For full relativistic magnetic calculations some space group operations of the input space group are invalid. FPLO reduces the symmetry to a subset which leaves the magnetic field unchanged. These are all operations which leave the field axis invariant. Among these there can be operations which invert the axis. For these an additional time reversal puts the field back into its original direction. Such combined operations have timerev==True. Note, that these operations do not act like normal space group representations. If a normal space group operation acts like \(H^{\vec{k}}=D^{\vec{k}+} H^{\alpha \vec{k}} D^{\vec{k}}\), time reversed ops act like \(H^{\vec{k}}=\left(D^{\vec{k}+} H^{-\alpha \vec{k}} D^{\vec{k}}\right)^{*}\) where \(D^{\vec{k}}\) is a product of the time reversal representation matrix and the space group representation matrix in the Wannier orbital space.

Note: that for non-magnetic full relativistic calculations time reversal is also a symmetry. The corresponding representation matrix is to found among the list of WFSymOp instances returned by Slabify.hamAtKPoint. For this operation the complex conjugation also has to be applied.

isimproper

If True this operation is improper; a rotation times inversion.

Dk

This is the representation matrix of \(\hat{g}=\{\alpha\mid\tau\}\) (or pure time reversal \(\Theta\) if present or a combined operation \(\Theta\hat{g}\)) in the space of the Wannier orbitals at a given k-point. If the operation is in the little group of the k-point the Hamiltonian fulfills \(H^{\vec{k}}=D^{\vec{k}+} H^{\vec{k}} D^{\vec{k}}\) otherwise \(H^{\vec{k}}=D^{\vec{k}+} H^{\alpha \vec{k}} D^{\vec{k}}\). If the operation is combined with time reversal an additional complex conjugation must be applied. If time reversal is in the little group we get \(H^{\vec{k}}=\left(D^{\vec{k}+} H^{\vec{k}} D^{\vec{k}}\right)^{*}\) otherwise \(H^{\vec{k}}=\left(D^{\vec{k}+} H^{-\alpha \vec{k}} D^{\vec{k}}\right)^{*}\). Also see timerev.

The Wannier orbital Bloch sums transform as \(\hat{g} w^{\vec{k}} =w^{\alpha \vec{k}} D^{\vec{k}}\) where \(w\) is the row vector of all WFs and the multiplication with \(D^{\vec{k}}\) a matrix multiplication. For time reversed operations this reads \(\Theta\hat{g} w^{\vec{k}} =w^{-\alpha \vec{k}} D^{\vec{k}} K_0\) where \(K_0\) indicates complex conjugation of everything which stands right of this operator as in \(\Theta\hat{g} w^{\vec{k}} C^{\vec{k}} =w^{-\alpha \vec{k}} D^{\vec{k}} C^{\vec{k}*} K_0\)

The matrix \(D^{\vec{k}}\) contains the fractional translational part tau. For additional translations by lattice vectors additional phase factors \(\exp(-i\alpha\vec{k} \cdot \vec{R})\) must be added. However, such a situation usually does not occur. For time reversed operations the minus sign becomes a plus sign.

Note: that the chosen gauge of the Hamiltonian affects the periodicity of the Hamiltonian and WFs in k-space. In particular these objects are only k-periodic for the periodic gauge. In the relative gauge additional phase factors would occure.

equivalentSites

The site number gisa=equivalentSites[isa] is related to isa by this operation.

BerryCurvatureData

class BerryCurvatureData

This helper class is returned by Slabify.hamAtKPoint if option makedhk==True. It contains the k-gradient of the Hamiltonian matrix and some additional data which are needed to calculate the Berry curvature with proper symmetry for periodic gauge.

__getitem__()

For compatibility with older pyfplo versions this object behaves a bit like a list of 3 complex numpy.ndarrays in that it can be indexed. The returned numpy.ndarrays are copies of the underlying data. So you cannot modify the BerryCurvatureData data.

Site

Site

Data

For convenience:

version

copy of pyfplo.common.version

Version

copy of pyfplo.common.Version

c_elements

copy of pyfplo.common.c_elements

BandPlot

copy of pyfplo.common.BandPlot

BandWeights

copy of pyfplo.common.BandWeights

WeightDefinition

copy of pyfplo.common.WeightDefinition

WeightDefinitions

copy of pyfplo.common.WeightDefinitions

BandFileContext

copy of pyfplo.common.BandFileContext

Vlevel

copy of pyfplo.common.Vlevel

OptionSet

copy of pyfplo.common.OptionSet

BandHeader

copy of pyfplo.common.BandHeader