# pyfplo.slabify

- Slabify
- BoxMesh
- EnergyContour
- FermiSurfaceOptions
- DensPlotContext
- GreenOptions
- WeylPoint
- BfieldConfig
- WFSymOp
- BerryCurvatureData
- Site
- Data

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: | bandplot (BandPlot) – see help of `pyfplo.common.BandPlot` suffix (str) – a file suffix for convenience

    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: | fso (`FermiSurfaceOptions`) – see help for `FermiSurfaceOptions` wds (WeightDefinitions) – see help for `pyfplo.common.WeightDefinitions` bandplot (BandPlot) – see help for `pyfplo.common.BandPlot` Only `lowerdepthdatalimit` and `upperdepthdatalimit` and bandplot restrictions (`pyfplo.common.BandPlot.setOutputRestrictions`) are used suffix (str) – a string which is appended the the end of the output file `+cuts_spin1` (and `+cuts_spin2`) basename. forcerecalculation (`bool`) – If a re-calculation is wanted use this argument or follow the file deletion hint of the program output.

    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: | nsigiter (int) – see `nsigiter` sigitertol (float) – see `sigitertol` sigitermethod (str) – see `sigitermethod`

  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.ndarray`s 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`