Model Diagnostics¶
chinook is designed with the capacity to do fairly extensive characterization of the tight-binding model being used. These scripts contain useful tools for understanding the model in more detail
Density of States¶
Several different approaches can be taken to the execution of density of states. One can follow Blöchl , partitioning the Brillouin zone into a cubic mesh, with each cube further divided into a group of identical tetrahedra. This approach provides the ability to perform matrix diagonalization over a fairly small number of k-points, as one interpolates within the tetrahedra to construct the full density of states. This is executed using dos_tetra defined below.
Alternatively, one can perform diagonalization explicitly only at the nodes of a mesh defined over the Brillouin zone, and apply some broadened peak at the eigenvalues of the Hamiltonian at each point. Generally, for reasonably narrow Gaussian broadening this requires a fairly dense k-mesh. In practice, the large supercells where diagonalization becomes costly are also associated with much smaller Brillouin zones, allowing for a smaller k-space sampling. We find this method to perform better in most cases. This is executed using the dos_broad function defined below.
Related tools are also available to find the Fermi level, given a specified electronic occupation (dos.find_EF using gaussian, and dos.EF_find using tetrahedra).
-
dos.band_contribution(eigenvals, w_domain, volume)¶ Compute the contribution over a single tetrahedron, from a single band, to the density of states
args:
- eigenvals: numpy array of float, energy values at corners
- w_domain: numpy array of float, energy domain
- volume: int, number of tetrahedra in the total mesh
return:
- DOS: numpy array of float, same length as w_domain
-
dos.def_dE(Eband)¶ If energy broadening is not passed for density-of-states calculation, compute a reasonable value based on the energy between adjacent energies in the tight-binding calculation
args:
- Eband: numpy array of float, indicating band energies
return:
dE: float, energy broadening, as the smallest average energy spacing over all bands.
-
dos.dos_broad(TB, NK, NE=None, dE=None, origin=array([0., 0., 0.]))¶ Energy-broadened discrete density of states calculation. The Hamiltonian is diagonalized over the kmesh defined by NK and states are summed, as energy-broadened Gaussian peaks, rather than delta functions.
args:
- TB: tight-binding model object
- NK: int, or tuple of int, indicating number of k-points
kwargs:
- NE: int, number of energy bins for final output
- dE: float, energy broadening of peaks, eV
- origin: numpy array of 3 float, indicating the origin of the mesh to be used,
relevant for example in kz-specific contributions to density of states
return:
- DOS: numpy array of float, density-of-states in states/eV
- Elin: numpy array of float, energy domain in eV
-
dos.dos_func(energy, epars)¶ Piecewise function for calculation of density of states
args:
- energy: numpy array of float (energy domain)
- epars: tuple of parameters: e[0],e[1],e[2],e[3],V_T,V_G being the ranked band energies for the tetrahedron,
as well as the volume of both the tetrahedron and the Brillouin zone, all float
return:
- numpy array of float giving DOS contribution from this tetrahedron
-
dos.dos_tetra(TB, NE, NK)¶ Generate a tetrahedra mesh of k-points which span the BZ with even distribution Diagonalize over this mesh and then compute the resulting density of states as prescribed in the above paper. The result is plotted, and DOS returned
args:
- TB: tight-binding model object
- NE: int, number of energy points
- NK: int or list of 3 int – number of k-points in mesh
return:
- Elin: linear energy array of float, spanning the range of the eigenspectrum
- DOS: numpy array of float, same length as Elin, density of states
-
dos.error_function(x0, x, sigma)¶ Integral over the gaussian function, evaluated from -infinity to x, using the scipy implementation of the error function
args:
- x0: float, centre of Gaussian, in eV
- x: numpy array of float, energy domain eV
- sigma: float, width of Gaussian, in eV
return:
- analytical form of integral
-
dos.find_EF_broad_dos(TB, NK, occ, NE=None, dE=None, origin=array([0., 0., 0.]))¶ Find the Fermi level of a model Hamiltonian, for a designated electronic occupation. Note this is evaluated at T=0, so EF is well-defined.
args:
- TB: tight-binding model object
- NK: int, or tuple of int, indicating number of k-points
- occ: float, desired electronic occupation
kwargs:
- NE: int, number of energy bins for final output
- dE: float, energy spacing of bins, in eV
- origin: numpy array of 3 float, indicating the origin of the mesh to be used,
relevant for example in kz-specific contributions to density of states
return:
- EF: float, Fermi level in eV
-
dos.find_EF_tetra_dos(TB, occ, dE, NK)¶ Use the tetrahedron-integration method to establish the Fermi-level, for a given electron occupation.
args:
- TB: instance of tight-binding model object from TB_lib
- occ: float, desired electronic occupation
- dE: estimate of energy precision desired for evaluation of the
Fermi-level (in eV)
- NK: int or iterable of 3 int, number of k points in mesh.
return:
EF: float, Fermi Energy for the desired occupation, to within dE of actual value.
-
dos.gaussian(x0, x, sigma)¶ Evaluate a normalized Gaussian function.
args:
- x0: float, centre of peak, in eV
- x: numpy array of float, energy domain in eV
- sigma: float, width of Gaussian, in eV
return:
- numpy array of float, gaussian evaluated.
-
dos.n_func(energy, epars)¶ Piecewise function for evaluating contribution of tetrahedra to electronic occupation number
args:
- energy: numpy array of float, energy domain
- epars: tuple of parameters: e[0],e[1],e[2],e[3],V_T,V_G being the ranked band energies for the tetrahedron,
as well as the volume of both the tetrahedron and the Brillouin zone, all float
return:
- numpy array of float, same length as energy, providing contribution of
tetrahedra to the occupation function
-
dos.n_tetra(TB, dE, NK, plot=True)¶ This function, also from the algorithm of Blochl, gives the integrated DOS at every given energy (so from bottom of bandstructure up to its top. This makes for very convenient and precise evaluation of the Fermi level, given an electron number)
args:
- TB: tight-binding model object
- dE: float, energy spacing (meV)
- NK: int, iterable of 3 int. number of k-points in mesh
- plot: bool, to plot or not to plot the calculated array
return:
- Elin: linear energy array of float, spanning the range of the eigenspectrum
- n_elect: numpy array of float, same length as Elin, integrated DOS
at each energy, i.e. total number of electrons occupied at each energy
-
dos.ne_broad_analytical(TB, NK, NE=None, dE=None, origin=array([0., 0., 0.]), plot=True)¶ Analytical evaluation of the occupation function. Uses scipy’s errorfunction executable to evaluate the analytical form of a Gaussian-broadened state’s contribution to the total occupation, at each energy
args:
- TB: tight-binding model object
- NK: int, or tuple of int, indicating number of k-points
kwargs:
- NE: int, number of energy bins for final output
- dE: float, energy spacing of bins, in eV
- origin: numpy array of 3 float, indicating the origin of the mesh to be used,
relevant for example in kz-specific contributions to density of states
- plot: boolean, default to True, if false, suppress plot output
return:
- nE: numpy array of float, occupied states
- Elin: numpy array of float, energy domain in eV
-
dos.ne_broad_numerical(TB, NK, NE=None, dE=None, origin=array([0., 0., 0.]))¶ Occupation function, as a numerical integral over the density of states function.
args:
- TB: tight-binding model object
- NK: int, or tuple of int, indicating number of k-points
kwargs:
- NE: int, number of energy bins for final output
- dE: float, energy spacing of bins, in eV
- origin: numpy array of 3 float, indicating the origin of the mesh to be used,
relevant for example in kz-specific contributions to density of states
return:
- ne: numpy array of float, integrated density-of-states at each energy
- Elin: numpy array of float, energy domain in eV
-
dos.pdos_tetra(TB, NE, NK, proj)¶ Partial density of states calculation. Follows same tetrahedra method, weighting the contribution of a given tetrahedra by the average projection onto the indicated user-defined projection. The average here taken as the sum over projection at the 4 vertices of the tetrahedra.
args:
- TB: tight-binding model object
- NE: int, number of energy bins
- NK: int, or iterable of 3 int, indicating the number of k-points
along each of the axes of the Brillouin zone
- proj: numpy array of float, 1D or 2D, c.f. proj_mat.
return:
- Elin: numpy array of float, with length NE, spanning the
range of the tight-binding bandstructure
- pDOS: numpy array of float, len NE, projected density of states
- DOS: numpy array of float, len NE, full density of states
-
dos.proj_avg(eivecs, proj_matrix)¶ Calculate the expectation value of the projection operator, for each of the eigenvectors, at each of the vertices, and then sum over the vertices. We use numpy.einsum to perform matrix multiplication and contraction.
args:
- eivecs: numpy array of complex float, 4xNxM, with M number of eigenvectors,
N basis dimension
- proj_matrix: numpy array of complex float, NxN in size
return:
- numpy array of M float, indicating the average projection over the 4
corners of the tetrahedron
-
dos.proj_mat(proj, lenbasis)¶ Define projection matrix for fast evaluation of the partial density of states weighting. As the projector here is diagonal, and represents a Hermitian matrix, it is by definition a real matrix operator.
args:
- proj: numpy array, either 1D (indices of projection), or 2D (indices of
projection and weight of projection)
- lenbasis: int, size of the orbital basis
return:
- projector: numpy array of float, lenbasis x lenbasis
Fermi Surface¶
We use a modified version of the method of marching tetrahedra to find the Fermi surface within the reciprocal lattice. The standard definition of the reciprocal space mesh runs over the primitive parallelepiped defined by the recriprocal lattice vectors. Shifts to the lattice origin are provided as an option. The loci of the Fermi surface are found for each tetrahedra, and used to assemble a continuous set of triangular patches which ultimately construct the Fermi surface for each band which crosses the Fermi level.
-
FS_tetra.EF_tetra(TB, NK, EF, degen=False, origin=None)¶ Generate a tetrahedra mesh of k-points which span the BZ with even distribution Diagonalize over this mesh and then compute the resulting density of states as prescribed in the above paper.
args:
- TB: tight-binding model object
- NK: int or list,tuple of 3 int indicating number of k-points in mesh
- EF: float, Fermi energy, or energy of interest
kwargs:
- degen: bool, flag for whether the bands are two-fold degenerate, as for
Kramers degeneracy
- origin: numpy array of 3 float, corresponding to the desired centre of the
plotted Brillouin zone
return:
- surfaces: dictionary of dictionaries: Each key-value pair corresponds
to a different band index. For each case, the value is a dictionary with key-value pairs:
- ‘pts’: numpy array of Nx3 float, the N coordinates of EF crossing
- ‘tris’: numpy array of Nx3 int, the triangulation of the surface
-
FS_tetra.FS_generate(TB, Nk, EF, degen=False, origin=None, ax=None)¶ Wrapper function for computing Fermi surface triangulation, and then plotting the result.
args:
- TB: tight-binding model object
- Nk: int, or tuple/list of 3 int, number of k-points in Brillouin zone mesh
- EF: float, Fermi energy, or constant energy level of interest
kwargs:
- degen: bool, flag for whether the bands are two-fold degenerate, as for
Kramers degeneracy
- origin: numpy array of 3 float, corresponding to the desired centre of the
plotted Brillouin zone
- ax: matplotlib Axes, option for plotting onto existing Axes
return:
- surfaces: dictionary of dictionaries: Each key-value pair corresponds
to a different band index. For each case, the value is a dictionary with key-value pairs:
- ‘pts’: numpy array of Nx3 float, the N coordinates of EF crossing
- ‘tris’: numpy array of Nx3 int, the triangulation of the surface
- ax: matplotlib Axes, for further modification
-
FS_tetra.fermi_surface_2D(TB, npts=100, kfix=(2, 0), energy=0, shift=array([0, 0, 0]), do_plot=True)¶ Generate a 2D contour of the Fermi surface, projected into one of the 3 cardinal planes. User specifies which b-vector to be normal to, and its fixed value. The user also specifies the ‘Fermi’ energy, and can shift the centre of the plot away from the origin if desired.
- args:
- TB: tight-binding model object
- npts: int or tuple of 2-int, number of kpoints in grid
- kfix: tuple of two numeric. First is b-vector index (0-base), second is the fixed value (float)
- energy: float, Fermi energy
- shift: numpy array of 3 float, shift vector, in units or b-vectors
- do_plot: boolean, option to suppress plot and only return the FS contours
- returns:
- ax: if do_plot, then a figure is generated and the axes object returned
- FS: if not do_plot, then a dictionary of contours, with keys indicating the associated band index, and
- values being the arrays of K points is returned
-
FS_tetra.get_kpts(TB, kfix, npts=100, shift=array([0, 0, 0]))¶ Get k-grid for Brillouin zone sampling
- args:
- TB: tight-binding model object
- kfix: tuple of two numeric. First is b-vector index (0-base), second is the fixed value (float)
- npts: int or tuple of 2-int, number of kpoints in grid
- shift: numpy array of 3 float, shift vector, in units or b-vectors
-
FS_tetra.heron(vert)¶ Heron’s algorithm for calculation of triangle area, defined by only the vertices
args:
- vert: numpy array of 3x3 indicating vertices of triangle
return:
- float, area of triangle
-
FS_tetra.sim_tri(vert)¶ Take 4 vertices of a quadrilateral and split into two alternative triangulations of the corners. Return the vertices of the triangulation which has the more similar areas between the two triangles decomposed.
args:
- vert: (4 by 3 numpy array (or list) of float) in some coordinate frame
return:
- tris[0] , tris[1]: two numpy arrays of size 3 by 3 float containing
the coordinates of a triangulation
Operators¶
A number of tools are included for characterizing the expectation values of observables, as projected onto the eigenstates of the model Hamiltonian. The main function here is O_path, which will compute the expectation value of an indicated operator (represented by the user as a Hermitian numpy array of complex float with the same dimensions as the orbital basis). The resulting values are then displayed in the form of a colourmap applied to the bandstructure calculation along the desired path in momentum space. Several common operators are defined to facilitate these calculations, such as spin \(\hat S_i\) and \(\left<\vec{L}\cdot\vec{S}\right>\) . In addition, fatbs uses O_path to produce a “fat bands” type spaghetti plot, where the colourscale reflects the orbital projection onto the bandstructure.
Created on Mon Mar 23 20:08:46 2020
@author: ryanday
-
operator_library.FS(TB, ktuple, Ef, tol, ax=None)¶ A simplified form of Fermi surface extraction, for proper calculation of this, chinook.FS_tetra.py is preferred. This finds all points in kmesh within a tolerance of the constant energy level.
args:
- TB: tight-binding model object
- ktuple: tuple of k limits, len (3), First and second should be iterable,
define the limits and mesh of k for kx,ky, the third is constant, float for kz
- Ef: float, energy of interest, eV
- tol: float, energy tolerance, float
- ax: matplotlib Axes, option for plotting onto existing Axes
return:
- pts: numpy array of len(N) x 3 indicating x,y, band index
- TB.Eband: numpy array of float, energy spectrum
- TB.Evec: numpy array of complex float, eigenvectors
- ax: matplotlib Axes, for further user modification
-
operator_library.LSmat(TB, axis=None)¶ Generate an arbitary L.S type matrix for a given basis. Uses same Yproj as the HSO in the chinook.H_library, but is otherwise different, as it supports projection onto specific axes, in addition to the full vector dot product operator.
Otherwise, the LiSi matrix is computed with i the axis index. To do this, a linear combination of L+S+,L-S-,L+S-,L-S+,LzSz terms are used to compute.
In the factors dictionary, the weight of these terms is defined. The keys are tuples of (L+/-/z,S+/-/z) in a bit of a cryptic way. For L, range (0,1,2) ->(-1,0,1) and for S range (-1,0,1) = S1-S2 with S1/2 = +/- 1 here
L+,L-,Lz matrices are defined for each l shell in the basis, transformed into the basis of cubic harmonics. The nonzero terms will then just be used along with the spin and weighted by the factor value, and slotted into a len(basis)xlen(basis) matrix HSO
args:
- TB: tight-binding object, as defined in TB_lib.py
- axis: axis for calculation as either ‘x’,’y’,’z’,None,
or float (angle in the x-y plane)
return:
- HSO: (len(basis)xlen(basis)) numpy array of complex float
-
operator_library.LdotS(TB, axis=None, ax=None, colourbar=True)¶ Wrapper for O_path for computing L.S along a vector projection of interest, or none at all.
args:
- TB: tight-binding obect
kwargs:
- axis: numpy array of 3 float, indicating axis, or None for full L.S
- ax: matplotli.Axes object for plotting
- colourbar: bool, display colourbar on plot
return:
- O: numpy array of Nxlen(basis) float, expectation value of operator
on each band over the kpath of TB.Kobj.
-
operator_library.Lm(l)¶ L- operator in the l,m_l basis, organized with (0,0) = |l,l>, (2*l,2*l) = |l,-l> The nonzero elements are on the upper diagonal
arg:
- l: int orbital angular momentum
return:
- M: numpy array (2*l+1,2*l+1) of real float
-
operator_library.Lp(l)¶ L+ operator in the l,m_l basis, organized with (0,0) = |l,l>, (2*l,2*l) = |l,-l> The nonzero elements are on the upper diagonal
arg:
- l: int orbital angular momentum
return:
- M: numpy array (2*l+1,2*l+1) of real float
-
operator_library.Lz(l)¶ Lz operator in the l,m_l basis
arg:
- l: int orbital angular momentum
return:
- numpy array (2*l+1,2*l+1)
-
operator_library.O_path(Operator, TB, Kobj=None, vlims=None, Elims=None, degen=False, plot=True, ax=None, colourbar=True, colourmap=None)¶ Compute and plot the expectation value of an user-defined operator along a k-path Option of summing over degenerate bands (for e.g. fat bands) with degen boolean flag
args:
- Operator: matrix representation of the operator (numpy array len(basis), len(basis) of complex float)
- TB: Tight binding object from TB_lib
kwargs:
- Kobj: Momentum object, as defined in chinook.klib.py
- vlims: tuple of 2 float, limits of the colourscale for plotting,
if default value passed, will compute a reasonable range
- Elims: tuple of 2 float, limits of vertical scale for plotting
- plot: bool, default to True, plot, or not plot the result
- degen: bool, True if bands are degenerate, sum over adjacent bands
- ax: matplotlib Axes, option for plotting onto existing Axes
- colourbar: bool, plot colorbar on axes, default to True
- colourmap: matplotlib colourmap,i.e. LinearSegmentedColormap
return:
- O_vals: the numpy array of float, (len Kobj x len basis) expectation values
- ax: matplotlib Axes, allowing for user to further modify
-
operator_library.O_surf(O, TB, ktuple, Ef, tol, vlims=(0, 0), ax=None)¶ Compute and plot the expectation value of an user-defined operator over a surface of constant-binding energy
Option of summing over degenerate bands (for e.g. fat bands) with degen boolean flag
args:
- O: matrix representation of the operator (numpy array len(basis), len(basis) of complex float)
- TB: Tight binding object from chinook.TB_lib.py
- ktuple: momentum range for mesh:
- ktuple[0] = (x0,xn,n),ktuple[1]=(y0,yn,n),ktuple[2]=kz
kwargs:
- vlims: limits for the colourscale (optional argument), will choose
a reasonable set of limits if none passed by user
- ax: matplotlib Axes, option for plotting onto existing Axes
return:
- pts: the numpy array of expectation values, of shape Nx3, with first
two dimensions the kx,ky coordinates of the point, and the third the expectation value.
- ax: matplotlib Axes, allowing for further user modifications
-
operator_library.S_vec(LB, vec)¶ Spin operator along an arbitrary direction can be written as n.S = nx Sx + ny Sy + nz Sz
args:
- LB: int, length of basis
- vec: numpy array of 3 float, direction of spin projection
return:
- numpy array of complex float (LB by LB), spin operator matrix
-
operator_library.Sz(TB, ax=None, colourbar=True)¶ Wrapper for O_path for computing Sz along a vector projection of interest, or none at all.
args:
- TB: tight-binding obect
kwargs:
- ax: matplotlib.Axes plotting object
- colourbar: bool, display colourbar on plot
return:
- O: numpy array of Nxlen(basis) float, expectation value of operator
on each band over the kpath of TB.Kobj.
-
operator_library.colourmaps()¶ Plot utility, define a few colourmaps which scale to transparent at their zero values
-
operator_library.degen_Ovals(Oper_exp, Energy)¶ In the presence of degeneracy, we want to average over the evaluated orbital expectation values–numerically, the degenerate subspace can be arbitrarily diagonalized during numpy.linalg.eigh. All degeneracies are found, and the expectation values averaged.
args:
- Oper_exp: numpy array of float, operator expectations
- Energy: numpy array of float, energy eigenvalues.
-
operator_library.fatbs(proj, TB, Kobj=None, vlims=None, Elims=None, degen=False, ax=None, colourbar=True, plot=True)¶ Fat band projections. User denotes which orbital index projection is of interest Projection passed either as an Nx1 or Nx2 array of float. If Nx2, first column is the indices of the desired orbitals, the second column is the weight. If Nx1, then the weights are all taken to be eqaul
args:
- proj: iterable of projections, to be passed as either a 1-dimensional
with indices of projection only, OR, 2-dimensional, with the second column giving the amplitude of projection (for linear-combination projection)
- TB: tight-binding object
kwargs:
- Kobj: Momentum object, as defined in chinook.klib.py
- vlims: tuple of 2 float, limits of the colorscale for plotting, default to (0,1)
- Elims: tuple of 2 float, limits of vertical scale for plotting
- plot: bool, default to True, plot, or not plot the result
- degen: bool, True if bands are degenerate, sum over adjacent bands
- ax: matplotlib Axes, option for plotting onto existing Axes
- colorbar: bool, plot colorbar on axes, default to True
return:
- Ovals: numpy array of float, len(Kobj.kpts)*len(TB.basis)
-
operator_library.is_numeric(a)¶ Quick check if object is numeric
args:
- a: numeric, float/int
return:
- bool, if numeric True, else False
-
operator_library.operator_projected_fermi_surface(TB, matrix, npts=100, kfix=(2, 0), energy=0, shift=array([0, 0, 0]), degen=True, fig=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, scale=20)¶ Simple 2D-projected FS with operator expectation values plot over the FS contours.
- args:
- TB: tight-binding object
- matrix: numpy array of complex float, operator matrix
- npts: number of k-points along axes of BZ
- kfix: fixed index of BZ. First value is projected reciprocal lattice vector (0,1,2), second is value (inverse Angstrom)
- energy: float, fixed value of energy (EF = 0 )
- shift: numpy array of 3 float. Shift of centre of plot
- degen: boolean, average over degenerate bands
- fig: matplotlib figure to plot on top of
- cmap: colourmap
- scale: multiplier for scatterplot point sizes
-
operator_library.surface_proj(basis, length)¶ Operator for computing surface-projection of eigenstates. User passes the orbital basis and an extinction length (1/e) length for the ‘projection onto surface’. The operator is diagonal with exponenential suppression based on depth.
For use with SLAB geometry only
args:
- basis: list, orbital objects
- cutoff: float, cutoff length
return:
- M: numpy array of float, shape len(TB.basis) x len(TB.basis)
Orbital Visualization¶
-
orbital_plotting.col_phase(vals)¶ Define the phase of a complex number
args:
- vals: complex float, or numpy array of complex float
return:
- float, or numpy array of float of same shape as vals, from -pi to pi
-
orbital_plotting.make_angle_mesh(n)¶ Quick utility function for generating an angular mesh over spherical surface
args:
- n: int, number of divisions of the angular space
return:
- th: numpy array of 2n float from 0 to pi
- ph: numpy array of 4n float from 0 to 2pi
-
orbital_plotting.rephase_wavefunctions(vecs, index=-1)¶ The wavefunction at different k-points can choose an arbitrary phase, as can a subspace of degenerate eigenstates. As such, it is often advisable to choose a global phase definition when comparing several different vectors. The user here passes a set of vectors, and they are rephased. The user has the option of specifying which basis index they would like to set the phasing. It is essential however that the projection onto at least one basis element is non-zero over the entire set of vectors for this rephasing to work.
args:
- vecs: numpy array of complex float, ordered as rows:vector index, columns: basis index
kwargs:
- index: int, optional choice of basis phase selection
return:
- rephase: numpy array of complex float of same shape as vecs
-
class
orbital_plotting.wavefunction(basis, vector)¶ This class acts to reorganize basis and wavefunction information in a more suitable data structure than the native orbital class, or the sake of plotting orbital wavefunctions. The relevant eigenvector can be redefined, so long as it represents a projection onto the same orbital basis set as defined previously.
args:
- basis: list of orbital objects
- vector: numpy array of complex float, eigenvector projected onto the basis orbitals
-
calc_Ylm(th, ph)¶ Calculate all spherical harmonics needed for present calculation
return:
- numpy array of complex float, of shape (len(self.harmonics),len(th))
-
find_centres()¶ Create a Pointer array of basis indices and the centres of these basis orbitals.
return:
- all_centres: list of numpy array of length 3, indicating unique positions in the basis set
- centre_pointers: list of int, indicating the indices of position array, associated with the
location of the related orbital in real space.
-
find_harmonics()¶ Create a pointer array of basis indices and the associated spherical harmonics, as well as aa more convenient vector form of the projections themselves, as lists of complex float
return:
- all_lm: list of int, l,m pairs of all spherical harmonics relevant to calculation
- lm_pointers: list of int, pointer indices relating each basis orbital projection to the
lm pairs in all_lm
- projectors: list of arrays of complex float, providing the complex projection of basis
onto the related spherical harmonics
-
plot_wavefunction(vertices, triangulations, colours, plot_ax=None, cbar_ax=None)¶ Plotting function, for visualizing orbitals.
args:
- vertices: numpy array of float, shape (len(centres), len(th)*len(ph), 3) locations of vertices
- triangulations: numpy array of int, indicating the vertices connecting each surface patch
- colours: numpy array of float, of shape (len(centres),len(triangles)) encoding the orbital phase for each surface patch of the plotting
- plot_ax: matplotlib Axes, for plotting on existing axes
- cbar_ax: matplotlib Axes, for use in drawing colourbar
return:
- plots: list of plotted surfaces
- plot_ax: matplotlib Axes, for further modifications
-
redefine_vector(vector)¶ Update vector definition
args:
- vector: numpy array of complex float, same length as self.vector
-
triangulate_wavefunction(n, plotting=True, ax=None)¶ Plot the wavefunction stored in the class attributes as self.vector as a projection over the basis of spherical harmonics. The radial wavefunctions are not explicitly included, in the event of multiple basis atom sites, the length scale is set by the mean interatomic distance. The wavefunction phase is encoded in the colourscale of the mesh plot. The user sets the smoothness of the orbital projection by the integer argument n
args:
- n: int, number of angles in the mesh: Theta from 0 to pi is divided 2n times, and
Phi from 0 to 2pi is divided 4n times
kwargs:
- plotting: boolean, turn on/off to display plot
- ax: matplotlib Axes, for plotting on existing plot
return:
- vertices: numpy array of float, shape (len(centres), len(th)*len(ph), 3) locations of vertices
- triangulations: numpy array of int, indicating the vertices connecting each surface patch
- colours: numpy array of float, of shape (len(centres),len(triangles)) encoding the orbital phase for each surface patch of the plotting
- ax: matplotlib Axes, for further modifications
tetrahedra¶
Created on Tue Sep 4 16:27:40 2018
@author: rday
-
tetrahedra.corners()¶ Establish the shortest main diagonal of a cube of points, so as to establish the main diagonal for tetrahedral partitioning of the cube
return:
main: tuple of 2 integers indicating the cube coordinates
cube: numpy array of 8 corners (8x3) float
-
tetrahedra.gen_mesh(avec, N)¶ Generate a mesh of points in 3-dimensional momentum space over the first Brillouin zone. These are defined first in terms of recirocal lattice vectors,
i.e. from 0->1 along each, and then are multiplied by the rec. latt. vectors themselves. Note that this implicitly provides a mesh which is not centred at zero, but has an origin at the rec. latt. vector (0,0,0)
args:
- avec: numpy array of 3x3 float, lattice vectors
- N: int, or tuple of 3 int, indicating the number of points along
each of the reciprocal lattice vectors
-
tetrahedra.mesh_tetra(avec, N)¶ An equivalent definition of a spanning grid over the Brillouin zone is just one which spans the reciprocal cell unit cell. Translational symmetry imposes that this partitioning is equivalent to the conventional definition of the Brillouin zone, with the very big advantage that we can define a rectilinear grid which spans this volume in a way which can not be done for most Bravais lattices in R3.
args:
- avec: numpy array of 3x3 float, lattice vectors
- N: int, or iterable of 3 int which define the density of the mesh
over the Brillouin zone.
return:
- pts: numpy array of Mx3 float, indicating the points in momentum space
at the vertices of the mesh
- mesh_tet: numpy array of Lx4 int, indicating the L-tetrahedra
which partition the grid
-
tetrahedra.neighbours(point)¶ For an unit cube, we can define the set of 3 nearest neighbours by performing the requisite modular sum along one of the three Cartesian axes. In this way, for an input point, we can extract its neighbours easily.
args:
- point: numpy array of 3 int, all either 0 or 1
return:
- numpy array of 3x3 int, indicating the neighbours of point on the
unit cube.
-
tetrahedra.not_point(point)¶ Inverse of point, defined in an N-dimensional binary coordinate frame
args:
- point: int or numpy array of int between 0 and 1
return:
- numpy array of int, NOT gate applied to the binary vector point
-
tetrahedra.propagate(i, Nr, Nc)¶ Distribute the generic corner numbering convention defined for a cube at the origin to a cube starting at some arbitrary point in our grid. Excludes the edge points as starting points, so that all cubes are within the grid.
args:
- i: int, index of origin
- Nr: int, number of rows in grid
- Nc: int, number of columns in grid
return:
- **numpy array of int, len 8 corresponding to the re-numbering of the
corners of the cube.
-
tetrahedra.tet_inds()¶ Generate, for a single cube, the tetrahedral designations, for the following conventional numbering:
6 o —- o 7 / / |- 4 o —- o5 o 3
- | /
0 o —- o 1
with 2 hidden from view (below 6, and behind the line-segment connecting 4-5). Here drawn with x along horizontal, z into plane, y vertical Defining the real-index spacing between adjacent cubes in a larger array, we can apply this simple prescription to define the spanning tetrahedra over the larger k-mesh
return:
- tetra_inds: numpy array of integer (6x4), with each
row containing the index of the 4 tetrahedral vertices. Together, for of a set of neighbouring points on a grid, we divide into a set of covering tetrahedra which span the volume of the cube.