Source code for sisl.physics.brillouinzone

# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at https://mozilla.org/MPL/2.0/.
"""Brillouin zone classes
=========================

The Brillouin zone objects are all special classes enabling easy manipulation
of an underlying physical quantity.

Quite often a physical quantity will be required to be averaged, or calculated individually
over a number of k-points. In this regard the Brillouin zone objects can help.

The BrillouinZone object allows direct looping of contained k-points while invoking
particular methods from the contained object.
This is best shown with an example:

>>> H = Hamiltonian(...)
>>> bz = BrillouinZone(H)
>>> bz.apply.array.eigh()

This will calculate eigenvalues for all k-points associated with the `BrillouinZone` and
return everything as an array. The `~sisl.physics.BrillouinZone.dispatch` property of
the `BrillouinZone` object has several use cases (here ``array`` is shown).

This may be extremely convenient when calculating band-structures:

>>> H = Hamiltonian(...)
>>> bs = BandStructure(H, [[0, 0, 0], [0.5, 0, 0]], 100)
>>> bs_eig = bs.apply.array.eigh()
>>> plt.plot(bs.lineark(), bs_eig)

and then you have all eigenvalues for all the k-points along the path.

Sometimes one may want to post-process the data for each k-point.
As an example lets post-process the DOS on a per k-point basis while
calculating the average:
 
>>> H = Hamiltonian(...)
>>> mp = MonkhorstPack(H, [10, 10, 10])
>>> E = np.linspace(-2, 2, 100)
>>> def wrap_DOS(eigenstate):
...    # Calculate the DOS for the eigenstates
...    DOS = eigenstate.DOS(E)
...    # Calculate the velocity for the eigenstates
...    v = eigenstate.velocity()
...    V = (v ** 2).sum(1)
...    return DOS.reshape(-1, 1) * v ** 2 / V.reshape(-1, 1)
>>> DOS = mp.apply.average.eigenstate(wrap=wrap_DOS, eta=True)

This will, calculate the Monkhorst pack k-averaged DOS split into 3 Cartesian
directions based on the eigenstates velocity direction. This method of manipulating
the result can be extremely powerful to calculate many quantities while running an
efficient `BrillouinZone` average. The `eta` flag will print, to stdout, a progress-bar.
The usage of the ``wrap`` method are also passed optional arguments, ``parent`` which is
``H`` in the above example. ``k`` and ``weight`` are the current k-point and weight of the
corresponding k-point. An example could be to manipulate the DOS depending on the k-point and
weight:

>>> H = Hamiltonian(...)
>>> mp = MonkhorstPack(H, [10, 10, 10])
>>> E = np.linspace(-2, 2, 100)
>>> def wrap_DOS(eigenstate, k, weight):
...    # Calculate the DOS for the eigenstates and weight by k_x and weight
...    return eigenstate.DOS(E) * k[0] * weight
>>> DOS = mp.apply.sum.eigenstate(wrap=wrap_DOS, eta=True)

When using wrap to calculate more than one quantity per eigenstate it may be advantageous
to use `~sisl.oplist` to handle cases of `BrillouinZone.apply.average` and `BrillouinZone.apply.sum`.

>>> H = Hamiltonian(...)
>>> mp = MonkhorstPack(H, [10, 10, 10])
>>> E = np.linspace(-2, 2, 100)
>>> def wrap_multiple(eigenstate):
...    # Calculate DOS/PDOS for eigenstates
...    DOS = eigenstate.DOS(E)
...    PDOS = eigenstate.PDOS(E)
...    # Calculate velocity for the eigenstates
...    v = eigenstate.velocity()
...    return oplist([DOS, PDOS, v])
>>> DOS, PDOS, v = mp.apply.average.eigenstate(wrap=wrap_multiple, eta=True)

Which does mathematical operations (averaging/summing) using `~sisl.oplist`.


Parallel calculations
---------------------

The ``apply`` method looping k-points may be explicitly parallelized.
To run parallel do:

>>> H = Hamiltonian(...)
>>> mp = MonkhorstPack(H, [10, 10, 10])
>>> with mp.apply.renew(pool=True) as par:
...     par.eigh()

This requires you also have the package ``pathos`` available.
The above will run in parallel using a default number of processors
in priority:

1. Environment variable ``SISL_NUM_PROCS``
2. Return value of ``os.cpu_count()``.

Note that this may interfere with BLAS implementation which defaults
to use all CPU's for threading. The total processors/threads that will
be created is ``SISL_NUM_PROCS * OMP_NUM_THREADS``. Try and ensure this is below
or equal to the actual core-count of your machine (or the number of requested
cores in a HPC environment).


Alternatively one can control the number of processors locally by doing:

>>> H = Hamiltonian(...)
>>> mp = MonkhorstPack(H, [10, 10, 10])
>>> with mp.apply.renew(pool=2) as par:
...     par.eigh()

which will request 2 processors (regardless of core-count).
As a last resort you can pass your own ``Pool`` of workers that
will be used for the parallel processing.

>>> from multiprocessing import Pool
>>> pool = Pool(4)
>>> H = Hamiltonian(...)
>>> mp = MonkhorstPack(H, [10, 10, 10])
>>> with mp.apply.renew(pool=pool) as par:
...     par.eigh()

The ``Pool`` should implement some standard methods that are
existing in the ``pathos`` enviroment such as ``Pool.restart`` and ``Pool.terminate``
and ``imap`` and ``uimap`` methods. See the ``pathos`` documentation for detalis.


   BrillouinZone
   MonkhorstPack
   BandStructure

"""

import itertools
from functools import reduce
from numbers import Integral, Real

import numpy as np
from numpy import argsort, dot, pi, sum

import sisl._array as _a
from sisl._dispatcher import ClassDispatcher
from sisl._internal import set_module
from sisl.grid import Grid
from sisl.lattice import Lattice
from sisl.messages import SislError, deprecate_argument, info, progressbar, warn
from sisl.oplist import oplist
from sisl.quaternion import Quaternion
from sisl.unit import units
from sisl.utils import batched_indices
from sisl.utils.mathematics import cart2spher, fnorm

__all__ = ["BrillouinZone", "MonkhorstPack", "BandStructure", "linspace_bz"]


class BrillouinZoneDispatcher(ClassDispatcher):
    r"""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 :ref:`physics.brillouinzone` for further examples.
    """

    pass


@set_module("sisl.physics")
def linspace_bz(bz, stop=None, jumps=None, jump_dk=0.05):
    r"""Convert points from a BZ object into a linear spacing of maximum value `stop`

    Parameters
    ----------
    bz : BrillouinZone, or ndarray
       the object containing the k-points
    stop : int or None, optional
       maximum value in the linear space, or if None, will return the cumulative
       distance of the k-points in the Brillouin zone
    jumps: array_like, optional
       whether there are any jumps for the k-points that should not be taken into account
    jump_dk: float or array_like, optional
       how much total distance the jump points will take

    Returns
    -------

    """
    if isinstance(bz, BrillouinZone):
        cart = bz.tocartesian(bz.k)
    else:
        cart = bz
    # calculate vectors between each neighbouring points
    dcart = np.diff(cart, axis=0, prepend=cart[0].reshape(1, -1))
    # calculate distances
    dist = (dcart**2).sum(1) ** 0.5

    if jumps is not None:
        # calculate the total distance
        total_dist = dist.sum()

        # Zero out the jumps
        dist[jumps] = 0.0
        total_dist = dist.sum()
        # correct jumps
        dist[jumps] = total_dist * np.asarray(jump_dk)

    # convert to linear scale
    if stop is None:
        return np.cumsum(dist)

    total_dist = dist.sum() / stop
    # Scale to total length of `stop`
    return np.cumsum(dist) / total_dist


@set_module("sisl.physics")
class BrillouinZone:
    """A class to construct Brillouin zone related quantities

    It takes any object (which has access to cell-vectors) as an argument
    and can then return the k-points in non-reduced units from reduced units.

    The object associated with the BrillouinZone object *has* to implement
    at least two different properties:

    1. `cell` which is the lattice vector
    2. `rcell` which is the reciprocal lattice vectors.

    The object may also be an array of floats in which case an internal
    `Lattice` object will be created from the cell vectors (see `Lattice` for
    details).

    Parameters
    ----------
    parent : object or array_like
       An object with associated ``parent.cell`` and ``parent.rcell`` or
       an array of floats which may be turned into a `Lattice`
    k : array_like, optional
       k-points that this Brillouin zone represents
    weight : scalar or array_like, optional
       weights for the k-points.
    """

[docs] def __init__(self, parent, k=None, weight=None): self.set_parent(parent) # define a bz_attr as though it has not been set self._bz_attr = ("", None) # Gamma point if k is None: self._k = _a.zerosd([1, 3]) self._w = _a.onesd(1) else: self._k = _a.arrayd(k).reshape(-1, 3) self._w = _a.emptyd(len(k)) if weight is None: weight = 1.0 / len(self._k) self._w[:] = weight
apply = BrillouinZoneDispatcher( "apply", # Do not allow class dispatching type_dispatcher=None, obj_getattr=lambda obj, key: getattr(obj.parent, key), )
[docs] def set_parent(self, parent): """Update the parent associated to this object Parameters ---------- parent : object or array_like an object containing cell vectors """ try: # It probably has the supercell attached parent.cell parent.rcell self.parent = parent except Exception: self.parent = Lattice(parent)
def __str__(self): """String representation of the BrillouinZone""" parent = str(self._parent_lattice()).replace("\n", "\n ") return f"{self.__class__.__name__}{{nk: {len(self)},\n {parent}\n}}" def _parent_lattice(self): """Ensures to return the lattice of the parent, however it may be nested""" parent = self.parent if isinstance(parent, Lattice): return parent if isinstance(parent.lattice, Lattice): return parent.lattice raise ValueError( f"{self.__class__.__name__} could not extract the lattice object, it must be located elsewhere." ) def __getstate__(self): """Return dictionary with the current state""" return { "parent_class": self.parent.__class__, "parent": self.parent.__getstate__(), "k": self._k.copy(), "weight": self._w.copy(), } def __setstate__(self, state): """Reset state of the object""" self._k = state["k"] self._w = state["weight"] parent = state["parent_class"].__new__(state["parent_class"]) parent.__setstate__(state["parent"]) self.set_parent(parent)
[docs] @staticmethod def 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 : list-like or 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 ------- BrillouinZone: even if all objects are not BrillouinZone objects the returned object will be. """ if isinstance(weight_scale, Real): weight_scale = [weight_scale] # check for lengths (scales cannot be longer!) if len(bzs) < len(weight_scale): raise ValueError( "BrillouinZone.merge requires length of weight_scale to be smaller or equal to " "the objects." ) if parent is None: parent = bzs[0].parent k = [] w = [] for bz, scale in itertools.zip_longest( bzs, weight_scale, fillvalue=weight_scale[-1] ): k.append(bz.k) w.append(bz.weight * scale) return BrillouinZone(parent, np.concatenate(k), np.concatenate(w))
[docs] def volume(self, ret_dim=False, periodic=None): """Calculate the volume of the full Brillouin zone of the parent This will return the volume 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, optional also return the dimensionality of the system periodic : array_like of int, optional estimate the volume using only the directions indexed by this array. The default value is `(self.parent.nsc > 1).nonzero()[0]`. Returns ------- vol : the volume of the Brillouin zone. Units are Ang^D with D being the dimensionality. For 0D it will return 0. dimensionality : int the dimensionality of the volume """ # default periodic array if periodic is None: periodic = (self.parent.nsc > 1).nonzero()[0] dim = len(periodic) vol = 0.0 if dim == 3: vol = self.parent.volume elif dim == 2: vol = self.parent.area(*periodic) elif dim == 1: vol = self.parent.length[periodic[0]] if ret_dim: return vol, dim return vol
[docs] @staticmethod def 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 or list of 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` 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 """ if isinstance(N, Integral): k = np.empty([N, 3], np.float64) for i in range(N): k[i] = func(parent, N, i, *args, **kwargs) else: # N must be some-kind of list like thingy Nk = np.prod(N) k = np.empty([Nk, 3], np.float64) for i, indices in enumerate(itertools.product(*map(range, N))): k[i] = func(parent, N, indices, *args, **kwargs) return BrillouinZone(parent, k)
[docs] @staticmethod def param_circle(parent, N_or_dk, kR, normal, origin, loop=False): r"""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 :math:`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 4 points a warning will be raised and the number of points increased to 4. kR : float radius of the k-point. In 1/Ang normal : array_like of float normal vector to determine the circle plane origin : array_like of float origin of the circle used to generate the circular parameterization loop : bool, optional 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 ------- BrillouinZone with the parameterized k-points. """ if isinstance(N_or_dk, Integral): N = N_or_dk else: # Calculate the required number of points N = int(kR**2 * pi / N_or_dk + 0.5) if N < 2: N = 2 info( "BrillouinZone.param_circle increased the number of circle points to 2." ) # Conversion object bz = BrillouinZone(parent) normal = _a.asarrayd(normal) origin = _a.asarrayd(origin) k_n = bz.tocartesian(normal) k_o = bz.tocartesian(origin) # Generate a preset list of k-points on the unit-circle if loop: radians = _a.aranged(N) / (N - 1) * 2 * np.pi else: radians = _a.aranged(N) / N * 2 * np.pi k = _a.emptyd([N, 3]) k[:, 0] = np.cos(radians) k[:, 1] = np.sin(radians) k[:, 2] = 0.0 # Now generate the rotation _, theta, phi = cart2spher(k_n) if theta != 0: pv = _a.arrayd([k_n[0], k_n[1], 0]) pv /= fnorm(pv) q = Quaternion(phi, pv, rad=True) * Quaternion(theta, [0, 0, 1], rad=True) else: q = Quaternion(0.0, [0, 0, k_n[2] / abs(k_n[2])], rad=True) # Calculate k-points k = q.rotate(k) k *= kR / fnorm(k).reshape(-1, 1) k = bz.toreduced(k + k_o) # The sum of weights is equal to the BZ area W = np.pi * kR**2 w = np.repeat([W / N], N) return BrillouinZone(parent, k, w)
[docs] def copy(self, parent=None): """Create a copy of this object, optionally changing the parent Parameters ---------- parent : optional change the parent """ if parent is None: parent = self.parent bz = self.__class__(parent, self._k, self.weight) bz._k = self._k.copy() bz._w = self._w.copy() return bz
@property def k(self): """A list of all k-points (if available)""" return self._k @property def weight(self): """Weight of the k-points in the `BrillouinZone` object""" return self._w @property def cell(self): return self.parent.cell @property def rcell(self): return self.parent.rcell
[docs] def tocartesian(self, k): """Transfer a k-point in reduced coordinates to the Cartesian coordinates Parameters ---------- k : list of float k-point in reduced coordinates Returns ------- numpy.ndarray in units of 1/Ang """ return dot(k, self.rcell)
[docs] def toreduced(self, k): """Transfer a k-point in Cartesian coordinates to the reduced coordinates Parameters ---------- k : list of float k-point in Cartesian coordinates Returns ------- numpy.ndarray in units of reciprocal lattice vectors ]-0.5 ; 0.5] (if k is in the primitive cell) """ return dot(k, self.cell.T / (2 * pi))
[docs] @staticmethod def in_primitive(k): """Move the k-point into the primitive point(s) ]-0.5 ; 0.5] Parameters ---------- k : array_like k-point(s) to move into the primitive cell Returns ------- numpy.ndarray all k-points moved into the primitive cell """ k = _a.arrayd(k) % 1.0 # Ensure that we are in the interval ]-0.5; 0.5] k[k > 0.5] -= 1 return k
[docs] def iter(self, ret_weight=False): """An iterator for the k-points and (possibly) the weights Parameters ---------- ret_weight : bool, optional if true, also yield the weight for the respective k-point Yields ------ kpt : k-point weight : weight of k-point, only if `ret_weight` is true. """ if ret_weight: for i in range(len(self)): yield self.k[i], self.weight[i] else: yield from self.k
__iter__ = iter def __len__(self): return len(self._k)
[docs] def write(self, sile, *args, **kwargs): """Writes k-points to a `~sisl.io.tableSile`. This allows one to pass a `tableSile` or a file-name. """ from sisl.io import tableSile kw = np.concatenate((self.k, self.weight.reshape(-1, 1)), axis=1) if isinstance(sile, tableSile): sile.write_data(kw.T, *args, **kwargs) else: with tableSile(sile, "w") as fh: fh.write_data(kw.T, *args, **kwargs)
@set_module("sisl.physics") class MonkhorstPack(BrillouinZone): r"""Create a Monkhorst-Pack grid for the Brillouin zone Parameters ---------- parent : object or array_like An object with associated `parent.cell` and `parent.rcell` or an array of floats which may be turned into a `Lattice` nkpt : array_like of ints a list of number of k-points along each cell direction displacement : float or array_like of float, optional the displacement of the evenly spaced grid, a single floating number is the displacement for the 3 directions, else they are the individual displacements size : float or array_like of float, optional the size of the Brillouin zone sampled. This reduces the boundaries of the Brillouin zone around the displacement to the fraction specified. I.e. `size` must be of values :math:`]0 ; 1]`. Defaults to the entire BZ. Note that this will also reduce the weights such that the weights are normalized to the entire BZ. centered : bool, optional whether the k-points are :math:`\Gamma`-centered (for zero displacement) trs : bool, optional whether time-reversal symmetry exists in the Brillouin zone. Examples -------- >>> lattice = Lattice(3.) >>> MonkhorstPack(lattice, 10) # 10 x 10 x 10 (with TRS) >>> MonkhorstPack(lattice, [10, 5, 5]) # 10 x 5 x 5 (with TRS) >>> MonkhorstPack(lattice, [10, 5, 5], trs=False) # 10 x 5 x 5 (without TRS) """
[docs] def __init__( self, parent, nkpt, displacement=None, size=None, centered=True, trs=True ): super().__init__(parent) if isinstance(nkpt, Integral): nkpt = np.diag([nkpt] * 3) elif isinstance(nkpt[0], Integral): nkpt = np.diag(nkpt) # Now we have a matrix of k-points if np.any(nkpt - np.diag(np.diag(nkpt)) != 0): raise NotImplementedError( f"{self.__class__.__name__} with off-diagonal components is not implemented yet" ) if displacement is None: displacement = np.zeros(3, np.float64) elif isinstance(displacement, Real): displacement = _a.fulld(3, displacement) else: # transfer the displacement to the primitive cell displacement = _a.asarrayd(displacement) displacement = self.in_primitive(displacement) if size is None: size = _a.onesd(3) elif isinstance(size, Real): size = _a.fulld(3, size) else: size = _a.asarrayd(size) # Retrieve the diagonal number of values Dn = np.diag(nkpt).astype(np.int32) if np.any(Dn) == 0: raise ValueError( f"{self.__class__.__name__} *must* be initialized with " "diagonal elements different from 0." ) i_trs = -1 if trs: # Figure out which direction to TRS nmax = 0 for i in [0, 1, 2]: if displacement[i] in [0.0, 0.5] and Dn[i] > nmax: nmax = Dn[i] i_trs = i if nmax == 1: i_trs = -1 if i_trs == -1: # If we still haven't decided (say for weird displacements) # simply take the one with the maximum number of k-points. i_trs = np.argmax(Dn) # Calculate k-points and weights along all directions kw = [ self.grid(Dn[i], displacement[i], size[i], centered, i == i_trs) for i in (0, 1, 2) ] # Now figure out if we have a 0 point along the TRS direction if trs: # Figure out if the first value is zero if abs(kw[i_trs][0][0]) < 1e-10: # Find indices we want to delete ik1, ik2 = (i_trs + 1) % 3, (i_trs + 2) % 3 k1, k2 = kw[ik1][0], kw[ik2][0] k_dup = _a.emptyd([k1.size, k2.size, 2]) k_dup[:, :, 0] = k1.reshape(-1, 1) k_dup[:, :, 1] = k2.reshape(1, -1) # Figure out the duplicate values # To do this we calculate the norm matrix # Note for a 100 x 100 k-point sampling this will produce # a 100 ^ 4 matrix ~ 93 MB # For larger k-point samplings this is probably not so good (300x300 -> 7.5 GB) k_dup = k_dup.reshape(k1.size, k2.size, 1, 1, 2) + k_dup.reshape( 1, 1, k1.size, k2.size, 2 ) k_dup = ( (k_dup[..., 0] ** 2 + k_dup[..., 1] ** 2) ** 0.5 < 1e-10 ).nonzero() # At this point we have found all duplicate points, to only take one # half of the points we only take the lower half # Also, the Gamma point is *always* zero, so we shouldn't do <=! # Now check the case where one of the directions is (only) the Gamma-point if kw[ik1][0].size == 1 and kw[ik1][0][0] == 0.0: # We keep all indices for the ik1 direction (since it is the Gamma-point! rel = (k_dup[1] > k_dup[3]).nonzero()[0] elif kw[ik2][0].size == 1 and kw[ik2][0][0] == 0.0: # We keep all indices for the ik2 direction (since it is the Gamma-point! rel = (k_dup[0] > k_dup[2]).nonzero()[0] else: rel = np.logical_and(k_dup[0] > k_dup[2], k_dup[1] > k_dup[3]) k_dup = (k_dup[0][rel], k_dup[1][rel], k_dup[2][rel], k_dup[3][rel]) del rel, k1, k2 else: # To signal we can't do this k_dup = None self._k = _a.emptyd((kw[0][0].size, kw[1][0].size, kw[2][0].size, 3)) self._w = _a.onesd(self._k.shape[:-1]) for i in (0, 1, 2): k = kw[i][0].reshape(-1, 1, 1) w = kw[i][1].reshape(-1, 1, 1) self._k[..., i] = np.rollaxis(k, 0, i + 1) self._w[...] *= np.rollaxis(w, 0, i + 1) del kw # Now clean up a few of the points if trs and k_dup is not None: # Create the correct indices in the ravelled indices k = [0] * 3 k[ik1] = k_dup[2] k[ik2] = k_dup[3] k_del = np.ravel_multi_index(tuple(k), self._k.shape[:-1]) k[ik1] = k_dup[0] k[ik2] = k_dup[1] k_dup = np.ravel_multi_index(tuple(k), self._k.shape[:-1]) del k self._k.shape = (-1, 3) self._w.shape = (-1,) if trs and k_dup is not None: self._k = np.delete(self._k, k_del, 0) self._w[k_dup] += self._w[k_del] self._w = np.delete(self._w, k_del) del k_dup, k_del # Store information regarding size and diagonal elements # This information is basically only necessary when # we want to replace special k-points self._diag = Dn # vector self._displ = displacement # vector self._size = size # vector self._centered = centered self._trs = i_trs
@property def displacement(self): """Displacement for this Monkhorst-Pack grid""" return self._displ def __str__(self): """String representation of `MonkhorstPack`""" p = self._parent_lattice() return ( "{cls}{{nk: {nk:d}, size: [{size[0]:.5f} {size[1]:.5f} {size[0]:.5f}], trs: {trs}," "\n diagonal: [{diag[0]:d} {diag[1]:d} {diag[2]:d}], displacement: [{disp[0]:.5f} {disp[1]:.5f} {disp[2]:.5f}]," "\n {lattice}\n}}" ).format( cls=self.__class__.__name__, nk=len(self), size=self._size, trs={0: "A", 1: "B", 2: "C"}.get(self._trs, "no"), diag=self._diag, disp=self._displ, lattice=str(p).replace("\n", "\n "), ) def __getstate__(self): """Return dictionary with the current state""" state = super().__getstate__() state["diag"] = self._diag state["displ"] = self._displ state["size"] = self._size state["centered"] = self._centered state["trs"] = self._trs return state def __setstate__(self, state): """Reset state of the object""" super().__setstate__(state) self._diag = state["diag"] self._displ = state["displ"] self._size = state["size"] self._centered = state["centered"] self._trs = state["trs"]
[docs] def copy(self, parent=None): """Create a copy of this object, optionally changing the parent Parameters ---------- parent : optional change the parent """ if parent is None: parent = self.parent bz = self.__class__( parent, self._diag, self._displ, self._size, self._centered, self._trs >= 0 ) # this is required due to replace calls bz._k = self._k.copy() bz._w = self._w.copy() return bz
[docs] @classmethod def grid(cls, n, displ=0.0, size=1.0, centered=True, trs=False): r"""Create a grid of `n` points with an offset of `displ` and sampling `size` around `displ` The :math:`k`-points are :math:`\Gamma` centered. Parameters ---------- n : int number of points in the grid. If `trs` is ``True`` this may be smaller than `n` displ : float, optional the displacement of the grid size : float, optional the total size of the Brillouin zone to sample centered : bool, optional if the points are centered trs : bool, optional whether time-reversal-symmetry is applied Returns ------- k : numpy.ndarray the list of k-points in the Brillouin zone to be sampled w : numpy.ndarray weights for the k-points """ # First ensure that displ is in the Brillouin displ = displ % 1.0 if displ > 0.5: displ -= 1.0 if displ < -0.5: displ += 1.0 # Centered _only_ has effect IFF # displ == 0. and size == 1 # Otherwise we resort to other schemes if displ != 0.0 or size != 1.0: centered = False # size *per k-point* dsize = size / n # We create the full grid, then afterwards we figure out TRS n_half = n // 2 if n % 2 == 1: k = _a.aranged(-n_half, n_half + 1) * dsize + displ else: k = _a.aranged(-n_half, n_half) * dsize + displ if not centered: # Shift everything by halve the size each occupies k += dsize / 2 # Move k to the primitive cell and generate weights k = cls.in_primitive(k) w = _a.fulld(n, dsize) # Check for TRS points if trs and np.any(k < 0.0): # Make all positive to remove the double conting terms k_pos = np.fabs(k) # Sort k-points and weights idx = argsort(k_pos) # Re-arange according to k value k_pos = k_pos[idx] w = w[idx] # Find indices of all equivalent k-points (tolerance of 1e-10 in reciprocal units) # Use the dsize to estimate the difference in positions idx_same = (np.diff(k_pos) < dsize * 1e-3).nonzero()[0] # The above algorithm should never create more than two duplicates. # Hence we can simply remove all idx_same and double the weight for all # idx_same + 1. w[idx_same + 1] *= 2 # Delete the duplicated k-points (they are already sorted) k = np.delete(k_pos, idx_same, axis=0) w = np.delete(w, idx_same) else: # Sort them, because it makes more visual sense idx = argsort(k) k = k[idx] w = w[idx] # Return values return k, w
[docs] def replace(self, k, mp, displacement=False, as_index=False, check_vol=True): r"""Replace a k-point with a new set of k-points from a Monkhorst-Pack grid This method tries to replace an area corresponding to `mp.size` around the k-point `k` such that the k-points are replaced. This enables one to zoom in on specific points in the Brillouin zone for detailed analysis. Parameters ---------- k : array_like k-point in this object to replace, if `as_index` is true, it will be regarded as integer positions of the k-points to replace, otherwise the indices of the k-points will be located individually (in chunks of 200 MB). mp : MonkhorstPack object containing the replacement k-points. displacement : array_like or bool, optional the displacment of the `mp` k-points. Needed for doing *lots* of replacements due to efficiency. Defaults to not displace anything. The inserted k-points will be `mp.k + displacement`. If True, it will use `k` as the displacement vector. For multiple k-point replacements each k-point will be replaced my `mp` with k as the displacement. as_index : bool, optional whether `k` is input as reciprocal k-points, or as indices of k-points in this object. check_vol : bool, optional whether to check the volume of the replaced k-point(s); by default the volume of each k-point is determined by the original ``size`` and ``nkpt`` values. However, when doing replacements of k-points these values are not kept for the individual k-points that were replaced, so subsequent replacements of these points will cause errors that effectively are not valid. Examples -------- This example creates a zoomed-in view of the :math:`\Gamma`-point by replacing it with a 3x3x3 Monkhorst-Pack grid. >>> lattice = Lattice(1.) >>> mp = MonkhorstPack(lattice, [3, 3, 3]) >>> mp.replace([0, 0, 0], MonkhorstPack(lattice, [3, 3, 3], size=1./3)) This example creates a zoomed-in view of the :math:`\Gamma`-point by replacing it with a 4x4x4 Monkhorst-Pack grid. >>> lattice = Lattice(1.) >>> mp = MonkhorstPack(lattice, [3, 3, 3]) >>> mp.replace([0, 0, 0], MonkhorstPack(lattice, [4, 4, 4], size=1./3)) This example creates a zoomed-in view of the :math:`\Gamma`-point by replacing it with a 4x4x1 Monkhorst-Pack grid. >>> lattice = Lattice(1.) >>> mp = MonkhorstPack(lattice, [3, 3, 3]) >>> mp.replace([0, 0, 0], MonkhorstPack(lattice, [4, 4, 1], size=1./3)) Raises ------ SislError if the size of the replacement `MonkhorstPack` grid is not compatible with the k-point spacing in this object. """ # First we find all k-points within k +- mp.size # Those are the points we wish to remove. # Secondly we need to ensure that the k-points we remove are occupying *exactly* # the Brillouin zone we wish to replace. if not isinstance(mp, MonkhorstPack): raise ValueError("Object 'mp' is not a MonkhorstPack object") if check_vol: # We can easily figure out the BZ that each k-point is averaging k_vol = self._size / self._diag # Compare against the size of this one # Since we can remove more than one k-point, we require that the # size of the replacement MP is an integer multiple of the # k-point volumes. k_int = mp._size / k_vol if not np.allclose(np.rint(k_int), k_int): raise SislError( f"{self.__class__.__name__}.reduce could not replace k-point, BZ " "volume replaced is not equivalent to the inherent k-point volume." ) # the size of the k-points that will be added s_size2 = self._size / 2 mp_size2 = mp._size / 2 dk = np.where(mp_size2 < s_size2, mp_size2, s_size2) dk.shape = (1, 3) # determine indices of k-point inputs k = np.asarray(k) if as_index: idx = k.ravel() k = self.k[idx] else: # find k-points in batches of 200 MB k = self.in_primitive(k).reshape(-1, 3) idx = batched_indices( self.k, k, atol=dk, batch_size=200, diff_func=self.in_primitive )[0] if ( self._trs >= 0 ): # TRS along a given axis, we can search the mirrored values idx2 = batched_indices( self.k, -k, atol=dk, batch_size=200, diff_func=self.in_primitive )[0] idx = np.concatenate((idx, idx2)) # we may find 2 indices for gamm-point in this case... not useful idx = np.unique(idx) if len(idx) == 0: raise SislError( f"{self.__class__.__name__}.reduce found no k-points to replace. " f"Searched with precision: {dk.ravel()}" ) # Idea of fast replacements is attributed @ahkole in #454, but the resulting code needed some # changes since that code was not stable againts *wrong input*, i.e. k=[0, 0, 0] # replacements. # determine the displacement vector if isinstance(displacement, bool): if displacement: displacement = k else: displacement = None elif displacement is not None: # convert to array displacement = _a.asarray(displacement).reshape(-1, 3) if displacement is None: displ_nk = 1 else: # Ensure we are in the central k-grid displacement = self.in_primitive(displacement) displ_nk = len(displacement) # Now we have the k-points we need to remove # Figure out if the total weight is consistent total_weight = self.weight[idx].sum() replace_weight = mp.weight.sum() * displ_nk atol = min(total_weight, replace_weight) * 1e-4 if abs(total_weight - replace_weight) < atol: weight_factor = 1.0 elif abs(total_weight - replace_weight * 2) < atol: weight_factor = 2.0 if self._trs < 0: info( f"{self.__class__.__name__}.reduce assumes that the replaced k-point has double weights." ) else: # print("k-point to replace: ", k.ravel()) # print("delta-k: ", dk.ravel()) # print("Found k-indices that will be replaced:") # print(idx) # print("k-points replaced:") # print(self.k[idx, :]) # print("weights replaced:") # print(self.weight[idx]) # print(self.weight.min(), self.weight.max()) # print(mp.weight.min(), mp.weight.max()) # print("Summed weights vs. replaced summed weights: ") # print(total_weight, replace_weight) # print(mp) raise SislError( f"{self.__class__.__name__}.reduce found inconsistent replacement weights " f"self={total_weight} vs. mp={replace_weight}. " f"Replacement indices: {idx}." ) # delete and append new k-points and weights if displacement is None: self._k = np.concatenate((np.delete(self._k, idx, axis=0), mp._k), axis=0) else: self._k = np.concatenate( ( np.delete(self._k, idx, axis=0), self.in_primitive(mp.k + displacement.reshape(-1, 1, 3)).reshape( -1, 3 ), ), axis=0, ) self._w = np.concatenate( (np.delete(self._w, idx), np.tile(mp._w * weight_factor, displ_nk)) )
@set_module("sisl.physics") class BandStructure(BrillouinZone): """Create a path in the Brillouin zone for plotting band-structures etc. Parameters ---------- parent : object or array_like An object with associated `parent.cell` and `parent.rcell` or an array of floats which may be turned into a `Lattice` points : array_like of float a list of points that are the *corners* of the path. Define a discontinuity in the points by adding a `None` in the list. divisions : int or array_like of int number of divisions in each segment. If a single integer is passed it is the total number of points on the path (equally separated). If it is an array_like input it must have length one less than `point`, in this case the total number of points will be ``sum(divisions) + 1`` due to the end-point constraint. names : array_like of str the associated names of the points on the Brillouin Zone path jump_dk: float or array_like, optional Percentage of ``self.lineark()[-1]`` that is used as separation between discontinued jumps in the band-structure. For band-structures with disconnected jumps the `lineark` and `lineartick` methods returns a separation between the disconnected points according to this percentage. Default value is 5% of the total distance. Alternatively an array equal to the number of discontinuity jumps may be passed for individual percentages. Keyword only, argument. Examples -------- >>> lattice = Lattice(10) >>> bs = BandStructure(lattice, [[0] * 3, [0.5] * 3], 200) >>> bs = BandStructure(lattice, [[0] * 3, [0.5] * 3, [1.] * 3], 200) >>> bs = BandStructure(lattice, [[0] * 3, [0.5] * 3, [1.] * 3], 200, ['Gamma', 'M', 'Gamma']) A disconnected band structure may be created by having None as the element. Note that the number of names does not contain the empty points (they are simply removed). Such a band-structure may be useful when one is interested in a discontinuous band structure. >>> bs = BandStructure(lattice, [[0, 0, 0], [0, 0.5, 0], None, [0.5, 0, 0], [0.5, 0.5, 0]], 200) """
[docs] @deprecate_argument( "name", "names", "argument 'name' has been deprecated in favor of 'names', please update your code.", "0.15.0", ) def __init__(self, parent, *args, **kwargs): # points, divisions, names=None): super().__init__(parent) lattice = self._parent_lattice() points = kwargs.pop("points", None) if points is None: if len(args) > 0: points, *args = args else: raise ValueError(f"{self.__class__.__name__} 'points' argument missing") # Extract the jump indices # In that case it is a disconnected path def is_empty(p): try: return len(p) == 0 except Exception: return p is None jump_idx = [] _points = [] for i, p in enumerate(points): if is_empty(p): if i > 0: # we can't have a jump at the first index jump_idx.append(i) else: _points.append(lattice._fill(p, dtype=np.float64)) # convert to exact array and correct for removed indices jump_idx = _a.arrayi(jump_idx) - _a.arangei(len(jump_idx)) # fill with correct points self.points = _a.arrayd(_points) del _points # clean-up for clarity divisions = kwargs.pop("divisions", None) if divisions is None: if len(args) > 0: divisions, *args = args else: raise ValueError( f"{self.__class__.__name__} 'divisions' argument missing" ) names = kwargs.pop("names", None) if names is None: if len(args) > 0: names, *args = args if len(args) > 0: raise ValueError( f"{self.__class__.__name__} unknown arguments after parsing 'points', 'divisions' and 'names': {args}" ) # Store empty split size self._jump_dk = np.asarray(kwargs.pop("jump_dk", 0.05)) if len(kwargs) > 0: raise ValueError( f"{self.__class__.__name__} unknown keyword arguments after parsing [points, divisions, names, jump_dk]: {list(kwargs.keys())}" ) # remove erroneous jumps if len(self.points) in jump_idx: jump_idx = jump_idx[:-1] if self._jump_dk.size > 1 and jump_idx.size != self._jump_dk.size: raise ValueError( f"{self.__class__.__name__} got inconsistent argument lengths (jump_dk does not match jumps in points)" ) # The jump-idx is equal to using np.split(self.points, jump_idx) # which then returns continuous sections self._jump_idx = jump_idx # If the array has fewer points we try and determine if self.points.shape[1] < 3: if self.points.shape[1] != np.sum(self.parent.nsc > 1): raise ValueError("Could not determine the non-periodic direction") # fix the points where there are no periodicity for i in (0, 1, 2): if self.parent.nsc[i] == 1: self.points = np.insert(self.points, i, 0.0, axis=1) # Ensure the shape is correct self.points.shape = (-1, 3) # Now figure out what to do with the divisions if isinstance(divisions, Integral): if divisions < len(self.points): raise ValueError( f"Can not evenly split {len(self.points)} points into {divisions} divisions, ensure division>=len(points)" ) # Get length between different k-points with a total length # of division dists = np.diff( linspace_bz(self.tocartesian(self.points), jumps=jump_idx, jump_dk=0.0) ) # Get floating point divisions divs_r = dists * divisions / dists.sum() # Convert to integers divs = np.rint(divs_r).astype(np.int32) # ensure at least 1 point along each division # 1 division means only the starting point divs[divs == 0] = 1 divs[jump_idx - 1] = 1 divs_sum = divs.sum() while divs_sum != divisions - 1: # only check indices where divs > 1 idx = (divs > 1).nonzero()[0] dk = dists[idx] / divs[idx] if divs_sum >= divisions: divs[idx[np.argmin(dk)]] -= 1 else: divs[idx[np.argmax(dk)]] += 1 divs_sum = divs.sum() divisions = divs[:] elif len(divisions) + 1 != len(self.points): raise ValueError( f"inconsistent number of elements in 'points' and 'divisions' argument. One less 'divisions' elements." ) self.divisions = _a.arrayi(divisions).ravel() if names is None: self.names = "ABCDEFGHIJKLMNOPQRSTUVXYZ"[: len(self.points)] else: self.names = names if len(self.names) != len(self.points): raise ValueError( f"inconsistent number of elements in 'points' and 'names' argument" ) # Calculate points dpoint = np.diff(self.points, axis=0) k = _a.emptyd([self.divisions.sum() + 1, 3]) i = 0 for ik, (divs, dk) in enumerate(zip(self.divisions, dpoint)): k[i : i + divs, :] = ( self.points[ik] + dk * _a.aranged(divs).reshape(-1, 1) / divs ) i += divs k[-1] = self.points[-1] # sanity check that should always be obeyed assert i + 1 == len(k) self._k = k self._w = _a.fulld(len(self.k), 1 / len(self.k))
[docs] def copy(self, parent=None): """Create a copy of this object, optionally changing the parent Parameters ---------- parent : optional change the parent """ if parent is None: parent = self.parent bz = self.__class__( parent, self.points, self.divisions, self.names, jump_dk=self._jump_dk ) return bz
def __getstate__(self): """Return dictionary with the current state""" state = super().__getstate__() state["points"] = self.points.copy() state["divisions"] = self.divisions.copy() state["jump_idx"] = self._jump_idx.copy() state["names"] = list(self.names) state["jump_dk"] = self._jump_dk return state def __setstate__(self, state): """Reset state of the object""" super().__setstate__(state) self.points = state["points"] self.divisions = state["divisions"] self.names = state["names"] self._jump_dk = state["jump_dk"] self._jump_idx = state["jump_idx"]
[docs] def insert_jump(self, *arrays, value=np.nan): """Return a copy of `arrays` filled with `value` at indices of discontinuity jumps Arrays with `value` in jumps is easier to plot since those lines will be naturally discontinued. For band structures without discontinuity jumps in the Brillouin zone the `arrays` will be return as is. It will insert `value` along the first dimension matching the length of `self`. For each discontinuity jump an element will be inserted. This may be useful for plotting since `np.nan` gets interpreted as a discontinuity in the graph thus removing connections between the segments. Parameters ---------- *arrays : array_like arrays will get `value` inserted where there are jumps in the band structure value : optional the value to be inserted at the jump points in the data array Examples -------- Create a bandrstructure with a discontinuity. >>> gr = geom.graphene() >>> bs = BandStructure(gr, [[0, 0, 0], [0.5, 0, 0], None, [0, 0, 0], [0, 0.5, 0]], 4) >>> data = np.zeros([len(bs), 10]) >>> data_with_jump = bs.insert_jump(data) >>> assert data_with_jump.shape == (len(bs)+1, 10) >>> np.all(data_with_jump[2] == np.nan) True """ # quick return if nothing needs changed if len(self._jump_idx) == 0: if len(arrays) == 1: return arrays[0] return arrays nk = len(self) full_jumps = np.cumsum(self.divisions)[self._jump_idx - 1] def _insert(array): array = np.asarray(array) # ensure dtype is equivalent as input array nans = np.empty(len(full_jumps), dtype=array.dtype) nans.fill(value) axis = array.shape.index(nk) shape = list(1 for _ in array.shape) shape[axis] = -1 return np.insert(array, full_jumps, nans.reshape(shape), axis=axis) # convert all arrays = tuple(_insert(array) for array in arrays) if len(arrays) == 1: return arrays[0] return arrays
[docs] def lineartick(self): """The tick-marks corresponding to the linear-k values Returns ------- numpy.ndarray the positions in reciprocal space determined by the distance between points See Also -------- lineark : Routine used to calculate the tick-marks. """ return self.lineark(True)[1:3]
[docs] def tolinear(self, k, ret_index=False, tol=1e-4): """Convert a k-point into the equivalent linear k-point via the distance Finds the index of the k-point in `self.k` that is closests to `k`. The returned value is then the equivalent index in `lineark`. This is very useful for extracting certain points along the band structure. Parameters ---------- k : array_like the k-point(s) to locate in the linear values ret_index : bool, optional whether the indices are also returned tol : float, optional when the found k-point has a distance (in Cartesian coordinates) is differing by more than `tol` a warning will be issued. The tolerance is in units 1/Ang. """ # Faster than to do sqrt all the time tol = tol**2 # first convert to the cartesian coordinates (for proper distances) ks = self.tocartesian(np.atleast_2d(k)) kk = self.tocartesian(self.k) # find closest values def find(k): dist = ((kk - k) ** 2).sum(-1) idx = np.argmin(dist) if dist[idx] > tol: warn( f"{self.__class__.__name__}.tolinear could not find a k-point within given tolerance ({self.toreduced(k)})" ) return idx idxs = [find(k) for k in ks] if ret_index: return self.lineark()[idxs], idxs return self.lineark()[idxs]
[docs] def lineark(self, ticks=False): """A 1D array which corresponds to the delta-k values of the path This is mainly meant for plotting but may be useful for finding out distances in the reciprocal lattice. Examples -------- >>> p = BandStructure(...) >>> eigs = Hamiltonian.eigh(p) >>> for i in range(len(Hamiltonian)): ... plt.plot(p.lineark(), eigs[:, i]) >>> p = BandStructure(...) >>> eigs = Hamiltonian.eigh(p) >>> lk, kt, kl = p.lineark(True) >>> plt.xticks(kt, kl) >>> for i in range(len(Hamiltonian)): ... plt.plot(lk, eigs[:, i]) Parameters ---------- ticks : bool, optional if `True` the ticks for the points are also returned See Also -------- linspace_bz : converts k-points into a linear distance parameterization Returns ------- linear_k : numpy.ndarray the positions in reciprocal space determined by the distance between points ticks : numpy.ndarray linear k-positions of the points, only returned if `ticks` is ``True`` ticklabels : list of str labels at `ticks`, only returned if `ticks` is ``True`` """ cum_divs = np.cumsum(self.divisions) # Calculate points # First we also need to calculate the jumps dK = linspace_bz( self, jumps=cum_divs[self._jump_idx - 1], jump_dk=self._jump_dk ) # Get label tick, in case self.names is a single string 'ABCD' if ticks: # Get number of points xtick = np.zeros(len(self.points), dtype=int) xtick[1:] = cum_divs # Ensure the returned label_tick is a copy return dK, dK[xtick], [a for a in self.names] return dK