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 fromSlabify
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 inSlabify.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 tozaxis
. 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 ofkscale
.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
ofSites
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.- bandplot (BandPlot) – see help of
-
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 forEnergyContour
- 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 ofkscale
. Use this to build your own non-uniform integration net. CallcalculateBulkProjectedEDC
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]
tok+zaxis*kzs[len(kzs)-1]
, where k is a k-point along the path spanned by the high symmetry points in bandplot.- bandplot (BandPlot) – see help of
-
calculateBulkProjectedFS
(fso, zaxis=[0, 0, 1], nz=30, kzs=None, suffix='', query=False)¶ Parameters: - fso (
FermiSurfaceOptions
) – see help forFermiSurfaceOptions
- 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 ofkscale
. Use this to build your own non-uniform integration net. CallcalculateBulkProjectedFS
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]
tok+zaxis*kzs[len(kzs)-1]
, where k is a k-point in the 2d mesh (from fso).- fso (
-
calculateFermiSurfaceCuts
(fso, wds=None, bandplot=None, suffix='', forcerecalculation=False)¶ Parameters: - fso (
FermiSurfaceOptions
) – see help forFermiSurfaceOptions
- wds (WeightDefinitions) – see help for
pyfplo.common.WeightDefinitions
- bandplot (BandPlot) – see help for
pyfplo.common.BandPlot
Onlylowerdepthdatalimit
andupperdepthdatalimit
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 runningcalculateFermiSurfaceCuts
again withforcerecalculation=False
- fso (
-
calculateEDC
(bandplot, energymesh, penetrationdepth=-1.0, greenoptions=None, suffix='')¶ Parameters: - bandplot (BandPlot) – see help for
pyfplo.common.BandPlot
- energymesh (
EnergyContour
) – see help forEnergyContour
- 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 forGreenOptions
- 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.
- bandplot (BandPlot) – see help for
-
calculateFermiSurfaceSpectralDensity
(fso, penetrationdepth=-1.0, greenoptions=None, suffix='')¶ Parameters: - fso (
FermiSurfaceOptions
) – see help forFermiSurfaceOptions
- 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 forGreenOptions
- suffix (str) – the output file suffix
Calculate the surface Fermi surface (k,k resolved spectral density) for a particular Fermi energy.
- fso (
-
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 byorbitalIndicesByDepth
.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 3numpy.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 forBoxMesh
- 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 isTrue
- 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.- box (
-
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 isTrue
- 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 singleint
) – 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 useenlarge
the TRIM points must not be given with respect to the enlarged cell! Seecalculate3dTIInvariants
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)
- gamma0 (3-vector of float) – marks the center of the 2D plane in units of
-
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 singleint
) – 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 usingcalculateBandStructure
and load the resulting+band...
intoxfbp
: 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 filesZ2_3dTI_homo*.cmd
which load all data intoxfbp
: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 thegauge
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 returnedWFSymOp
. See also makewfsymops. These indices are unique identifiers independend of the structure and are printed whenprepare
-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 isTrue
. The actual type ofdhk
isBerryCurvatureData
. - 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, whereH
,S
andB
are returned byhamAtKPoint
with options makesigma and makexcfield. This is approximate since to obtainB
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 isTrue
. The first 3 components ofbasisconnection
are the basis Berry connection and the last 3 the basis Berry curvature. To convert the connection inrelative
gauge to the position operator add the matrices returned bywannierCenterMatrix
. See argumentgauge
above to remember whygauge=forcerelative
needs to be used whenrelative
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 ofWFSymOp
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 isTrue
.
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 areTrue
a tuple is returned containing H and all additionally requested operators in the order defined by the make…-arguments in theargument 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 atuple
of return valuesExample
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 alist`of 3 complex `numpy.ndarrays
which contain dH/dk as returned byhamAtKPoint
- makef (int) – if
True
return band wise Berry curvature - basisconnection (
list
of 6 complexnumpy.ndarrays
) – basisconnection and curvature as returned byhamAtKPoint
- 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 anumpy.ndarray
of eigenvectorsC[:,i]
is thei
-th eigenvector. If makef isTrue
the third tuple elements is anumpy.ndarray
containing the Berry curvature whereF[:,n]
is the band-n
Berry curvature (a 3-vector). The dimension of the Hamiltonian isnvdim
.Return type: tuple
of return valuesDiagonalize 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. SeeberryCurvature
.dhk must be obtained from
hamAtKPoint
. New in version 19: Note, that in theperiodic
gauge (used for h and dhk) a correction for the position operator is used when calculating the Berry curvature. Ifgauge=forcerelative
was used inhamAtKPoint
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 valuesFor 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 valuesGiven 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.
- E (sequence of
-
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 alist`of 3 complex `numpy.ndarrays
which containdH/dk as returned byhamAtKPoint
- subspace (sequence of
int
) – a list of zeros and ones, which describes the wanted subspace orNone
for the whole space. - basisconnection (
list
of 6 complexnumpy.ndarrays
) – basisconnection and curvature as returned byhamAtKPoint
- 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
) – ifTrue
details of the curvature contributions are returned
Returns: (F,divergenttermsdetected) – F is a
numpy.ndarray
ofshape=(3,nvdim)
and contains the band wise Berry curvature. divergenttermsdetected is now alwaysFalse
. It is still there to not break older code. Sorry. If basisconnection is given andreturndetails==True
the returned tuple is(F,Fdetails,divergenttermsdetected)
, where Fdetails is adict
of the various terms F is made of. The dict values are againnumpy.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 fromhamAtKPoint
) the full curvature is calculated (Wang2006) not only the D-D part. Ifreturndetails==True
the separate terms F is made off are returned in adict
.If subspace is
None
the Berry curvature is returned for all bands. If it is a list of sizelen(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.- E (sequence of
-
findWeylPoints
(box, homos, tol=0.001)¶ Parameters: - box (
BoxMesh
) – a definition of the box and its mesh. - homos (
list
or singleint
) – 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 usecalculateChernNumberInSphere
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)
- box (
-
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 isTrue
- 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
ornumpy.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. Ifcutlayersat[1]<cutlayersat[0]
it is not applied. (see Structure Manipulation Algorithm)Type: list
of 2float
-
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 aftercutlayersat
was applied. Use xfplo on the resulting=.in
file to find the site numbers. (see Structure Manipulation Algorithm)Type: list
ofint
-
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)
wherea
is obtained from+hamdata
or more precisely from the conventional cell. One can setkscale
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
andzinterval=[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
andzinterval
for docs.
-
mesh
(scale, xfirst=False)¶ Parameters: - scale (float) – all points are multiplied with scale
- xfirst (int) – if
True
invert the loop order (make thex
-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
, theny
, the outer loop isz
. If xfirst isTrue
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
andk
run over all posible interval values andind
runs over all mesh indices depending on the loop order.
-
relToAbs
(*args)¶ Parameters: args (3 float
orlist
ornumpy.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
orlist
ornumpy.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 thex
-loop the outer loop).
Returns: xyz – a list of 3 indices such that the mesh position
xyz[0]
,xyz[1]
,xyz[2]
has pointmesh(ind)
Return type: numpy.ndarray
of 3int
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 (seemesh
).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 arepyfplo.slabify.Slabify.kscale
. (See the scale parameter inmesh
)
-
xinterval
¶ numpy.ndarray
: Interval along the first (x
) axis of the box as defined byxaxis
. The units often arepyfplo.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 byyaxis
. The units often arepyfplo.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 byzaxis
. The units often arepyfplo.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 projectionsType: 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: 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 definee0
,e1
andne
first or simply use the constructor.Type: float
orstr
-
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 the010/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) insemislab
mode the k-plane axis must be orthogonal!-
setPlane
(xaxis=[1, 0, 0], yaxis=[0, 1, 0], origin=[0, 0, 0])¶ Parameters: Convenience function: Define the k-plane.
-
setMesh
(nx=100, xinterval=[-0.5, 0.5], ny=100, yinterval=[-0.5, 0.5])¶ Parameters: 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 isx
. 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 andind
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
) orNone
Returns: band file context
Return type: 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. TheDensPlotContext.close()
method can be called explicitly.The best way to use it is in a
with
-statement. Then it is closed automatically after thewith
-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 ofFermiSurfaceOptions
are not orthogonal one wants to give two axes fo the same angle as arguments toopenDensPlotFile
to orient the data in the plotting plane. The origin in the plotting plane can also be given. The length ofplotxaxis
andplotyaxis
determines the length scale when plotted withxfbp
.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
.- filename (str) – the name of the
-
__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 vectorxaxis
. The units arepyfplo.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 vectoryaxis
. The units arepyfplo.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 asfermienergy
). The smaller this value the higher the k-point density must be. Often a good value isinterval-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. Forsemislab
the axes must be orthogonalType: 3-vector of float
-
yaxis
¶ the first (
y
) axis spanning the k-plane. It will be automatically normalized. Forsemislab
the axes must be orthogonalType: 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
- nsigiter (int) – see
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 (seeSlabify.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 pythonmap
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
andtau
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 whenSlabify.prepare
-ing the structure.
-
isinlittlegroup
¶ This is
True
if the operation is in the little group of the k-point for whichSlabify.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
andtau
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 bySlabify.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 optionmakedhk==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 complexnumpy.ndarrays
in that it can be indexed. The returnednumpy.ndarray
s are copies of the underlying data. So you cannot modify theBerryCurvatureData
data.
-
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