sisl.physics.GaussKronrodGrid

class sisl.physics.GaussKronrodGrid

Bases: AdaptiveBrillouinZone

Adaptive Gauss-Kronrod cubature of the Brillouin zone

The integrated Brillouin zone is split in boxes. Each box is integrated with a tensor product of the \(n\)-point Kronrod rule with the \(g\)-point Gauss rule embedded (\(g\) is gauss, \(n = 2g+1\)). A box uses \(n^d\) k-points for \(d\) integrated directions. The error of a box is the largest absolute element of the difference between the Kronrod and the Gauss results. The box with the largest error is halved along the direction with the largest error, until the summed error is below the tolerance.

This is efficient for integrands with few sharp features (e.g. small broadening), mainly in 1 or 2 dimensions. In 1 dimension, gauss=10 often needs fewer k-points than gauss=7. gauss=3 has fewer k-points in a box (\(7^3 = 343\) in 3 dimensions, \(15^3 = 3375\) for gauss=7), so it can halve boxes with a small max_k.

The integration domain is the reduced Brillouin zone \(]-1/2 ; 1/2]\) along the integrated directions. The k-point component along the other directions is 0.

With time-reversal symmetry (trs), only \([0 ; 1/2]\) along the first integrated direction is integrated, and the weights are doubled. As for MonkhorstPack, the result is only correct when \(f(\mathbf k)\) and \(f(-\mathbf k)\) have the same integral, e.g. a Green function symmetrized as \((\mathbf G + \mathbf G^T)/2\) after the integration.

Parameters:

parent (object or ndarray) – An object with associated parent.cell and parent.rcell or an array of floats which may be turned into a Lattice
nkpt – initial number of k-points along each lattice direction, per reciprocal lattice vector (as for MonkhorstPack). A direction with 0 k-points is not integrated. A scalar is the same value for the periodic directions of the parent lattice (pbc); the other directions are not integrated. If None, it is \(n\) (one box along each periodic direction). The value is rounded up to a multiple of \(n\) (one box has \(n\) k-points along a direction). With trs, the integrated half of the direction gets half the k-points. The initial k-points are always integrated, also when they exceed max_k.
max_k – a scalar is the maximum number of integrand evaluations for one integration. A sequence of 3 values is the maximum k-point density along each lattice direction, in k-points per reciprocal lattice vector: a box of width \(h\) (in reduced units) is not halved along a direction when \(n/(h/2)\) is larger than the value of that direction. The box is then halved along the direction with the next largest error. The total number is then not limited. Values for directions that are not integrated are not used.
trs – whether time-reversal symmetry is used. If None, it is true for (un)-polarized calculations, and false for non-colinear and higher spin configurations (as for MonkhorstPack).
atol – absolute tolerance of the integral (largest absolute element error)
rtol – relative tolerance of the integral, relative to the largest absolute element. For a sequence (oplist) result, the errors and the largest element are taken over all elements. Elements with very different magnitudes should then be scaled, or use atol.
gauss – number of Gauss nodes \(g\) of the Gauss-Kronrod rule, one of 3, 7 (QUADPACK qk15) or 10 (QUADPACK qk21). The Kronrod rule has \(n = 2g+1\) nodes.

Examples

The Green function of a 1D chain, with a small broadening, integrated along the first lattice direction only:

>>> H = Hamiltonian(geom.sc(1.0, Atom(1, R=1.1)))
>>> H.set_nsc([3, 1, 1])
>>> H.construct([[0.1, 1.1], [0.0, -1.0]])
>>> gk = GaussKronrodGrid(H, [15, 0, 0], atol=1e-8)
>>> z = 0.5 + 1e-3j
>>> def G(k):
...     return np.linalg.inv(z * H.Sk(k, format="array") - H.Hk(k, format="array"))
>>> Gavg = gk.integrate(G)
>>> print(f"{Gavg[0, 0]:.6f}")
0.000069-0.516398j
>>> len(gk.last_bz())  # the k-points that were used
255

Integrate the first two directions, and limit the k-point density to that of a 1000 x 200 grid:

>>> gk = GaussKronrodGrid(H, [15, 15, 0], max_k=[1000, 200, 1])

A 3D integration with the 7-point Kronrod rule, one box has 343 k-points:

>>> gk = GaussKronrodGrid(Lattice(1), gauss=3)
>>> print(f"{gk.integrate(lambda k: (k**2).sum()):.10f}")
0.2500000000
>>> len(gk.last_bz())
343

Plotting

`plot`	Plotting functions for the `BrillouinZone` class.
`plot.bands`([H, extra_vars, ...])	Creates a `BandsData` object and then plots a `BandsPlot` from it.
`plot.fatbands`([H, ...])	Creates a `BandsData` object and then plots a `FatbandsPlot` from it.
`plot.sites`([axes, sites, ...])	Builds a `SitesPlot` by setting the value of "sites_obj" to the current object.

\(k\)-point calculations

`apply`	Loop over all k-points by applying parent methods for all k.
`apply.array`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.average`(method)	Dispatch the method by an adaptive integration (`integrate`)
`apply.dataarray`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.grid`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.iter`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.list`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.ndarray`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.none`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.oplist`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.sum`(method)	Raise a TypeError, the dispatch needs fixed k-points
`apply.xarray`(method)	Raise a TypeError, the dispatch needs fixed k-points

Methods

`copy`([parent])	Create a copy of this object, optionally changing the parent
`in_primitive`(k)	Move the k-point into the primitive point(s) ]-0.5 ; 0.5]
`integrate`(func)	Integrate func over the Brillouin zone with adaptive Gauss-Kronrod boxes
`iter`([ret_weight])	An iterator for the k-points and (possibly) the weights
`last_bz`()	The k-points and weights of the last `integrate` call
`merge`(bzs[, weight_scale, parent])	Merge several BrillouinZone objects into one
`param_circle`(parent, N_or_dk, kR, normal, origin)	Create a parameterized k-point list where the k-points are generated on a circle around an origin
`parametrize`(parent, func, N, args, *kwargs)	Generate a new `BrillouinZone` object with k-points parameterized via the function func in N separations
`set_parent`(parent)	Update the parent associated to this object
`tocartesian`([k])	Transfer a k-point in reduced coordinates to the Cartesian coordinates
`toreduced`(k)	Transfer a k-point in Cartesian coordinates to the reduced coordinates
`volume`([ret_dim, axes])	Calculate the volume of the BrillouinZone, optionally only on some axes axes
`write`(sile, args, *kwargs)	Writes k-points to a `tableSile`.

Attributes

`cell`
`k`	Not defined for adaptive grids, raises TypeError
`rcell`
`trs`	Whether time-reversal symmetry is used
`weight`	Not defined for adaptive grids, raises TypeError

classmethod param_circle(parent, N_or_dk, kR, normal, origin, loop=False)

Create a parameterized k-point list where the k-points are generated on a circle around an origin

The generated circle is a perfect circle in the reciprocal space (Cartesian coordinates). To generate a perfect circle in units of the reciprocal lattice vectors one can generate the circle for a diagonal supercell with side-length \(2\pi\), see example below.

Parameters:

parent (Lattice, or LatticeChild) – the parent object
N_or_dk (int) – number of k-points generated using the parameterization (if an integer), otherwise it specifies the discretization length on the circle (in 1/Ang), If the latter case will use less than 2 points a warning will be raised and the number of points increased to 2.
kR (float) – radius of the k-point. In 1/Ang
normal (ndarray of float) – normal vector to determine the circle plane
origin (ndarray of float) – origin of the circle used to generate the circular parameterization
loop (bool) – whether the first and last point are equal

Examples

>>> lattice = Lattice([1, 1, 10, 90, 90, 60])
>>> bz = BrillouinZone.param_circle(lattice, 10, 0.05, [0, 0, 1], [1./3, 2./3, 0])

To generate a circular set of k-points in reduced coordinates (reciprocal

>>> lattice = Lattice([1, 1, 10, 90, 90, 60])
>>> bz = BrillouinZone.param_circle(lattice, 10, 0.05, [0, 0, 1], [1./3, 2./3, 0])
>>> bz_rec = BrillouinZone.param_circle(2*np.pi, 10, 0.05, [0, 0, 1], [1./3, 2./3, 0])
>>> bz.k[:, :] = bz_rec.k[:, :]

Returns:

with the parameterized k-points.

Return type:

BrillouinZone

Parameters:

N_or_dk (int | float)
kR (float)
normal (Sequence[float])
origin (Sequence[float])
loop (bool)

static in_primitive(k)

Move the k-point into the primitive point(s) ]-0.5 ; 0.5]

Parameters:: k (ndarray) – k-point(s) to move into the primitive cell
Returns:: all k-points moved into the primitive cell
Return type:: numpy.ndarray

static merge(bzs, weight_scale=1.0, parent=None)

Merge several BrillouinZone objects into one

The merging strategy only stores the new list of k-points and weights. Information retained in the merged objects will not be stored.

Parameters:

bzs (list-like of BrillouinZone objects) – each element is a BrillouinZone object with bzs[i].k and bzs[i].weight fields.
weight_scale (float | Sequence[float]) – these are matched item-wise with bzs and applied to. Internally itertools.zip_longest(fillvalue=weight_scale[-1]) will be used to extend for all bzs.
parent (object, optional) – Associated parent in the returned object, will default to bzs[0].parent

Returns:

even if all objects are not BrillouinZone objects the returned object will be.

Return type:

BrillouinZone

static parametrize(parent, func, N, *args, **kwargs)

Generate a new BrillouinZone object with k-points parameterized via the function func in N separations

Generator of a parameterized Brillouin zone object that contains a parameterized k-point list.

Parameters:

parent (Lattice, or LatticeChild) – the object that the returned object will contain as parent
func (callable) –
method that parameterizes the k-points, must at least accept three arguments, 1. parent: object 2. N: total number of k-points 3. i: current index of the k-point (starting from 0)

the function must return a k-point in 3 dimensions.
N (int | Sequence[int]) – number of k-points generated using the parameterization, or a list of integers that will be looped over. In this case arguments N and i in func will be lists accordingly.
*args – additional arguments passed directly to func
**kwargs – additional keyword arguments passed directly to func

Return type:

BrillouinZone

Examples

Simple linear k-points

>>> def func(sc, N, i):
...    return [i/N, 0, 0]
>>> bz = BrillouinZone.parametrize(1, func, 10)
>>> assert len(bz) == 10
>>> assert np.allclose(bz.k[-1, :], [9./10, 0, 0])

For double looping, say to create your own grid

>>> def func(sc, N, i):
...    return [i[0]/N[0], i[1]/N[1], 0]
>>> bz = BrillouinZone.parametrize(1, func, [10, 5])
>>> assert len(bz) == 50

apply.array(method): Raise a TypeError, the dispatch needs fixed k-points

apply.average(method): Dispatch the method by an adaptive integration (integrate)

apply.dataarray(method): Raise a TypeError, the dispatch needs fixed k-points

apply.grid(method): Raise a TypeError, the dispatch needs fixed k-points

apply.iter(method): Raise a TypeError, the dispatch needs fixed k-points

apply.list(method): Raise a TypeError, the dispatch needs fixed k-points

apply.ndarray(method): Raise a TypeError, the dispatch needs fixed k-points

apply.none(method): Raise a TypeError, the dispatch needs fixed k-points

apply.oplist(method): Raise a TypeError, the dispatch needs fixed k-points

apply.sum(method): Raise a TypeError, the dispatch needs fixed k-points

apply.xarray(method): Raise a TypeError, the dispatch needs fixed k-points

copy(parent=None)

Create a copy of this object, optionally changing the parent

Parameters:

parent (optional) – change the parent
bz (BrillouinZone)

Return type:

BrillouinZone

integrate(func)[source]

Integrate func over the Brillouin zone with adaptive Gauss-Kronrod boxes

The integration starts with the boxes from nkpt. Then it halves the box with the largest error, until the summed error of all boxes is below max(atol, rtol * max|I|), or until max_k stops it.

Parameters:: func (Callable) – the integrand, called as func(k) with k a reduced k-point with shape (3,). It must return a scalar, an array with a fixed shape, or a sequence (list, tuple or oplist, not numpy.ndarray) of these. A sequence is integrated element by element, and its elements can have different shapes. The errors include all elements. Each returned value is used at once and not stored, so func can return the same array (changed in place) at each call.
Returns:: the Brillouin zone average of func. It is an oplist when func returns a sequence.
Return type:: numpy.ndarray or oplist
Raises:: ValueError – if func returns sequences with different lengths, or arrays with different shapes

Notes

The integration keeps one integrated value for each box, and a few values for the box that it evaluates (6 values for \(d=1\), \(6+d\) for \(d>1\) integrated directions). The memory thus scales with the number of boxes times the size of a value of func.

Warns:: SislWarning – when max_k stops the integration before the error is below the tolerance. The returned value is then the last value. Also when func returns a value that is not finite (NaN or inf); the integration then stops at once.
Parameters:: func (Callable)
Return type:: TypeAliasForwardRef(’~numpy.ndarray’) | oplist

See also

FatbandsPlot: The plot class used to generate the plot.
BandsData: The class to which data is converted.

plot.sites(axes=['x', 'y', 'z'], sites=None, sites_style=[], sites_scale=1.0, sites_name='Sites', sites_colorscale=None, drawing_mode=None, show_cell=False, cell_style={}, nsc=(1, 1, 1), sites_ndim_scale=(1, 1, 1), dataaxis_1d=None, arrows=(), backend='plotly')

Builds a SitesPlot by setting the value of “sites_obj” to the current object.

Plots sites from an object that can be parsed into a geometry.

The only differences between this plot and a geometry plot is the naming of the inputs and the fact that there are no options to plot bonds.

Parameters:

axes (Axes) – The axes to project the sites to.
sites (AtomsIndex) – The sites to plot. If None, all sites are plotted.
sites_style (Sequence[AtomsStyleSpec]) – List of style specifications for the sites. See the showcase notebooks for examples.
sites_scale (float) – Scaling factor for the size of all sites.
sites_name (str) – Name to give to the trace that draws the sites.
sites_colorscale (Optional[Colorscale]) – Colorscale to use for the sites in case the color attribute is an array of values. If None, the default colorscale is used for each backend.
drawing_mode (Literal['scatter', 'balls', 'line', None]) – The method used to draw the sites.
show_cell (Literal['box', 'axes', False]) – Mode to display the reciprocal cell. If False, the cell is not displayed.
cell_style (StyleSpec) – Style specification for the reciprocal cell. See the showcase notebooks for examples.
nsc (tuple[int, int, int]) – Number of unit cells to display in each direction.
sites_ndim_scale (tuple[float, float, float]) – Scaling factor for the size of the sites for different dimensionalities (1D, 2D, 3D).
dataaxis_1d (Optional[Union[np.ndarray, Callable]]) – Only meaningful for 1D plots. The data to plot on the Y axis.
arrows (Sequence[AtomArrowSpec]) – List of arrow specifications to display. See the showcase notebooks for examples.
backend – The backend to use to generate the figure.

Return type:

SitesPlot

set_parent(parent)

Update the parent associated to this object

Parameters:: parent (object or ndarray) – an object containing cell vectors
Return type:: None

tocartesian(k=None)

Transfer a k-point in reduced coordinates to the Cartesian coordinates

Parameters:: k (TypeAliasForwardRef('~numpy.typing.ArrayLike') | None) – k-point in reduced coordinates, defaults to this objects k-points.
Returns:: in units of 1/Ang
Return type:: numpy.ndarray

toreduced(k)

Transfer a k-point in Cartesian coordinates to the reduced coordinates

Parameters:: k (list of float) – k-point in Cartesian coordinates
Returns:: in units of reciprocal lattice vectors ]-0.5 ; 0.5] (if k is in the primitive cell)
Return type:: numpy.ndarray

volume(ret_dim=False, axes=None)

Calculate the volume of the BrillouinZone, optionally only on some axes axes

This will return the volume of the Brillouin zone, depending on the dimensions of the system. Here the dimensions of the system is determined by how many dimensions have auxilliary supercells that can contribute to Brillouin zone integrals. Therefore the returned value will have differing units depending on dimensionality.

Parameters:

ret_dim (bool) – also return the dimensionality of the system
axes (Literal[0, 1, 2] | ~typing.Literal['a', 'b', 'c'] | ~collections.abc.Sequence[~typing.Literal[0, 1, 2, 'a', 'b', 'c']] | None) – estimate the volume using only the directions indexed by this array. The default axes are only the periodic ones (self.parent.pbc.nonzero()[0]). Hence the units might not necessarily be 1/Ang^3.

Returns:

vol – the volume of the Brillouin zone. Units are 1/Ang^D with D being the dimensionality. For 0D it will return 0.
dimensionality (int) – the dimensionality of the volume

Return type:

float | tuple[float, int]

write(sile, *args, **kwargs)

Writes k-points to a tableSile.

This allows one to pass a tableSile or a file-name.

Parameters:

bz (BrillouinZone)
sile (SileLike)

Return type:

None

apply

Loop over all k-points by applying parent methods for all k.

This allows potential for running and collecting various computationally heavy methods from a single point on all k-points.

The apply method will dispatch the parent methods through all k-points and passing k as arguments to the parent methods in a straight-forward manner.

For instance to iterate over all eigenvalues of a Hamiltonian

>>> H = Hamiltonian(...)
>>> bz = BrillouinZone(H)
>>> for ik, eigh in enumerate(bz.apply.eigh()):
...    # do something with eigh which corresponds to bz.k[ik]

By default the apply method exposes a set of dispatch methods:

apply.iter, the default iterator module
apply.average reduced result by averaging (using BrillouinZone.weight as the weight per k-point.
apply.sum reduced result without weighing
apply.array return a single array with all values; has len equal to number of k-points
apply.none, specialized method that is mainly useful when wrapping methods
apply.list same as apply.array but using Python list as return value
apply.oplist using sisl.oplist allows greater flexibility for mathematical operations element wise
apply.datarray if xarray is available one can retrieve an xarray.DataArray instance

Please see Brillouin zone for further examples.

property cell: ndarray

property k: Not defined for adaptive grids, raises TypeError

plot: Plotting functions for the BrillouinZone class.

property rcell: ndarray

property trs: bool: Whether time-reversal symmetry is used

property weight: Not defined for adaptive grids, raises TypeError