Coverage for tbkit/kspace.py: 100%
1376 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 06:13 +0000
1from __future__ import annotations
3from concurrent.futures import ThreadPoolExecutor
4from itertools import product
5from typing import Literal, Sequence, overload
7import numpy as np
8from numpy.typing import ArrayLike, NDArray
9import matplotlib.pyplot as plt
10from matplotlib.figure import Figure
11from matplotlib.collections import LineCollection, PolyCollection
12from matplotlib.colors import LogNorm, Normalize
13from mpl_toolkits.mplot3d.art3d import Poly3DCollection
14import scipy.linalg as LA
15import scipy.sparse as sp
16from scipy.integrate import cumulative_trapezoid
17from scipy.special import expit, spence
18import tbkit.error_handling as error_handling
19import tbkit.dos as dos
20import tbkit.occupation as occupation
21import tbkit.neighbours as neighbours
22import tbkit.values as values
23from tbkit.lattice import Lattice
24from tbkit.transport import _green_surface_bulk
27PI = np.pi
29#: Pauli matrices (plus the identity, key ``'0'``), for building spinful
30#: hoppings/onsite terms (spin-orbit coupling, Zeeman splitting, ...) when
31#: :class:`KSpace` is constructed with ``spin=True``.
32PAULI = {
33 '0': np.eye(2, dtype='c16'),
34 'x': np.array([[0., 1.], [1., 0.]], dtype='c16'),
35 'y': np.array([[0., -1j], [1j, 0.]], dtype='c16'),
36 'z': np.array([[1., 0.], [0., -1.]], dtype='c16'),
37}
40#################################
41# CLASS KSPACE
42#################################
45def reciprocal_vectors(prim_vec: Sequence[tuple[float, ...]]) -> list[tuple[float, ...]]:
46 r'''
47 Get the reciprocal lattice vectors :math:`\mathbf{b}_i` such that
48 :math:`\mathbf{a}_i\cdot\mathbf{b}_j = 2\pi\delta_{ij}`, lying in the
49 span of the :math:`\mathbf{a}_i`.
51 :param prim_vec: List of one/two/three tuples. Primitive vectors (see class **lattice**).
53 :returns:
54 * **rec_vec** -- List of one/two/three tuples. Reciprocal vectors, with
55 as many components as the primitive vectors.
56 '''
57 if len(prim_vec[0]) == 3:
58 # b = 2 pi (A^+)^T: the dual basis within the span of the a_i
59 rec = 2 * PI * np.linalg.pinv(np.array(prim_vec, dtype='f8')).T
60 return [tuple(float(c) for c in b) for b in rec]
61 if len(prim_vec) == 1:
62 ax, ay = prim_vec[0]
63 norm2 = ax ** 2 + ay ** 2
64 return [(2*PI*ax/norm2, 2*PI*ay/norm2)]
65 (a1x, a1y), (a2x, a2y) = prim_vec
66 area = a1x * a2y - a1y * a2x
67 b1 = (2*PI*a2y/area, -2*PI*a2x/area)
68 b2 = (-2*PI*a1y/area, 2*PI*a1x/area)
69 return [b1, b2]
72class KSpace():
73 r'''
74 Build and solve the Tight-Binding Bloch Hamiltonian :math:`H(\mathbf{k})`
75 of a periodic lattice defined by the class **lattice**.
77 Hoppings are defined between orbitals of the unit cell, separated by a
78 lattice vector :math:`\mathbf{R} = n_1\mathbf{a}_1+n_2\mathbf{a}_2`:
80 .. math::
82 H_{ij}(\mathbf{k}) = \sum_{\mathbf{R}} t_{ij}(\mathbf{R})\,
83 e^{i\mathbf{k}\cdot\mathbf{R}}
85 :param lat: **lattice** class instance. Only *unit_cell* and *prim_vec*
86 are used (the instance need not call *get_lattice*).
87 :param spin: Boolean. Default value False. If True, every site of
88 *unit_cell* carries a spin-1/2 degree of freedom (*norb* doubles to
89 ``2*len(unit_cell)``, ordered site-major: orbitals ``2*i, 2*i+1``
90 are the up/down components of site *i*). *set_onsite* and
91 *set_hopping* then accept 2x2 (spin) matrices in addition to plain
92 numbers, to build spin-orbit coupling or Zeeman terms -- see
93 :data:`PAULI` for ready-made Pauli matrices.
95 Example usage::
97 # graphene, nearest-neighbor hopping t
98 DX, DY = 0.5 * 3 ** 0.5, 0.5
99 unit_cell = [{'tag': 'a', 'r0': (0., 0.)}, {'tag': 'b', 'r0': (DX, DY)}]
100 prim_vec = [(2*DX, 0.), (DX, 1.5)]
101 lat = Lattice(unit_cell=unit_cell, prim_vec=prim_vec)
102 gra = KSpace(lat)
103 gra.set_hopping([{'i': 0, 'j': 1, 'R': (0, 0), 't': 1.},
104 {'i': 0, 'j': 1, 'R': (-1, 0), 't': 1.},
105 {'i': 0, 'j': 1, 'R': (0, -1), 't': 1.}])
106 '''
108 workers = 1 # threads diagonalizing over k, see set_workers
110 def __init__(self, lat: Lattice, spin: bool = False) -> None:
111 error_handling.lat(lat)
112 error_handling.boolean(spin, 'spin')
113 error_handling.independent(lat.prim_vec)
114 self.lat = lat
115 self.dim = len(lat.prim_vec)
116 self.space_dim = len(lat.prim_vec[0])
117 self.spin = spin
118 self.n_sites = len(lat.unit_cell)
119 self.norb = 2*self.n_sites if spin else self.n_sites
120 self.tags = np.array([dic['tag'] for dic in lat.unit_cell])
121 self.onsite = np.zeros(self.norb, 'c16')
122 # spin-off-diagonal onsite terms (in-plane Zeeman, onsite Rashba):
123 # they have no place on the diagonal `onsite` array.
124 self._onsite_offdiag = np.zeros((self.norb, self.norb), 'c16')
125 self._hop_const: list = [] # list of (i, j, R_cartesian (np.ndarray), t)
126 # value functions (see tbkit.values): (i, j, R_cartesian, function,
127 # hermitian), with arrays of sites i, j and bond vectors R
128 self._hop_values: list = []
129 #: Parameters of the value functions (see *set_params*).
130 self.params: dict = {}
131 self._nonreciprocal = False # set_hopping(hermitian=False) was used
132 self._overlap_hop: list = [] # overlaps, as _hop: list of (i, j, R_cartesian, s)
133 self.rec_vec = reciprocal_vectors(lat.prim_vec)
134 # Orthonormal basis (columns) of the span of prim_vec, in which k is
135 # given: the identity when the lattice fills its space, a1/|a1| for a
136 # chain, the Gram-Schmidt basis of (a1, a2) for a 2D lattice in 3D.
137 self.k_basis = np.linalg.qr(np.array(lat.prim_vec, dtype='f8').T)[0]
138 if self.dim == self.space_dim:
139 self.k_basis = np.eye(self.space_dim)
140 elif self.dim == 1:
141 a1 = np.asarray(lat.prim_vec[0], dtype='f8')
142 self.k_basis = (a1 / np.linalg.norm(a1))[:, None]
143 else:
144 # keep the first axis along a1 (QR may flip signs)
145 signs = np.sign(np.diag(np.array(lat.prim_vec, dtype='f8') @ self.k_basis))
146 self.k_basis = self.k_basis * np.where(signs == 0, 1., signs)[None, :]
147 #: Reciprocal vectors in the k coordinates used by *get_ham* (rows).
148 self.rec_vec_k = np.array(self.rec_vec, dtype='f8') @ self.k_basis
149 self.ks = np.array([]) # k-points of the last band-structure calculation
150 self.ks_dist = np.array([]) # cumulative distance along the k-path
151 self.nodes = np.array([]) # positions, along ks_dist, of the k-path nodes
152 self.en = np.array([]) # bands, shape (len(ks), norb)
154 @property
155 def _hop(self) -> list:
156 '''
157 Private. The hoppings, as a list of (i, j, R_cartesian, t): those
158 given as numbers, then the value functions evaluated with *params*.
159 Assigning it replaces both.
160 '''
161 return self._hops(self.params)
163 @_hop.setter
164 def _hop(self, hops: list) -> None:
165 self._hop_const = list(hops)
166 self._hop_values = []
168 def _hops(self, params: dict) -> list:
169 '''
170 Private method. *_hop*, with the value functions evaluated with *params*.
171 '''
172 hops = list(self._hop_const)
173 if not self._hop_values:
174 return hops
175 tau = np.array([dic['r0'] for dic in self.lat.unit_cell], dtype='f8').reshape(self.n_sites, -1)
176 for i, j, R, func, hermitian in self._hop_values:
177 site_i = values.sites(tau[i], self.tags[i], i)
178 site_j = values.sites(tau[j] + R, self.tags[j], j)
179 t = values.evaluate(func, (site_i, site_j), params, len(i), self.spin)
180 if self.spin:
181 pairs = [(2*ii + a, 2*jj + b, RR, tt[a, b], np.conj(tt[a, b]))
182 for ii, jj, RR, tt in zip(i, j, R, t) for a in range(2) for b in range(2)]
183 else:
184 pairs = [(ii, jj, RR, tt, np.conj(tt)) for ii, jj, RR, tt in zip(i, j, R, t)]
185 for ii, jj, RR, tt, tc in pairs:
186 hops.append((int(ii), int(jj), RR, tt))
187 if hermitian:
188 hops.append((int(jj), int(ii), -RR, tc))
189 return hops
191 def set_params(self, **params) -> None:
192 '''
193 Set the parameters of the value functions of *set_hopping*: every
194 method (*get_bands*, *chern_number*, *finite_ham*, ...) then uses
195 these values. *get_ham* also takes them per call.
197 :param params: Parameter values, added to (or replacing those of) *params*.
199 Example usage::
201 ks.set_params(m=0.5)
202 ks.chern_number()
203 '''
204 self.params.update(params)
206 def set_onsite(self, dict_onsite: dict[str, complex | Sequence[complex]]) -> None:
207 '''
208 Set the onsite energies, by sublattice tag.
210 :param dict_onsite: Dictionary. key: tag, val: onsite energy
211 (a plain number), or, if ``spin=True``, a plain number (applied
212 equally to both spins), a pair ``(E_up, E_down)`` of numbers (a
213 spin splitting along z), or a 2x2 complex matrix (a general spin
214 structure, e.g. an in-plane Zeeman field built from :data:`PAULI`).
215 A complex onsite energy (gain/loss), or a non-Hermitian 2x2
216 block, makes :math:`H(\\mathbf{k})` non-Hermitian: *get_bands*
217 then returns complex energies, sorted by their real part.
219 Example usage::
221 kag.set_onsite({'a': 1., 'b': -1.})
222 # spinful: same onsite energy for both spins on 'a', a Zeeman
223 # splitting along z on 'b':
224 kag_spin.set_onsite({'a': 1., 'b': (1., -1.)})
225 # spinful: an in-plane Zeeman field on 'a':
226 kag_spin.set_onsite({'a': Bx*PAULI['x']})
227 '''
228 error_handling.set_onsite_kspace(dict_onsite, self.lat.tags, self.spin)
229 for tag, val in dict_onsite.items():
230 sites = np.where(self.tags == tag)[0]
231 if self.spin:
232 if np.ndim(val) == 0:
233 block = val * PAULI['0']
234 elif np.ndim(val) == 2:
235 block = np.asarray(val, 'c16')
236 else:
237 block = np.diag(np.asarray(val, 'c16'))
238 self.onsite[2*sites] = block[0, 0]
239 self.onsite[2*sites + 1] = block[1, 1]
240 for site in sites:
241 self._onsite_offdiag[2*site, 2*site + 1] = block[0, 1]
242 self._onsite_offdiag[2*site + 1, 2*site] = block[1, 0]
243 else:
244 self.onsite[sites] = val
246 def set_hopping(self, list_hop: list[dict], hermitian: bool = True) -> None:
247 r'''
248 Set the hoppings between orbitals of the unit cell.
250 Only one representative of each hopping needs to be given: its
251 Hermitian conjugate (:math:`j\to i`, :math:`\mathbf{R}\to-\mathbf{R}`)
252 is added automatically (unless ``hermitian=False``).
254 :param list_hop: List of dictionaries with keys ('i', 'j', 'R', 't'):
256 * 'i', 'j': Positive integers. Site indices within the unit cell
257 (following the order of *unit_cell*).
258 * 'R': Tuple of one/two integers :math:`(n_1, n_2)`. Lattice vector
259 :math:`\mathbf{R}=n_1\mathbf{a}_1+n_2\mathbf{a}_2` separating the
260 two sites.
261 * 't': Complex number, or, if ``spin=True``, either a complex
262 number (spin-independent hopping) or a 2x2 complex matrix (a
263 general, possibly spin-mixing, hopping -- e.g. built from
264 :data:`PAULI` for Rashba or intrinsic spin-orbit coupling).
265 Or a value function ``t(site_i, site_j, **params)``,
266 evaluated by ``get_ham(k, **params)`` (and by every other
267 method, with the parameters of *set_params*): *site_i* and
268 *site_j* are structured arrays of the bonds' end points,
269 *site_i* in the home cell and *site_j* in cell
270 :math:`\mathbf{R}` (fields 'x', 'y', ('z'), 'tag', 'index' --
271 the site index within *unit_cell*). It is called once with
272 every bond of the dictionaries that share it, and returns one
273 value, or one per bond (or, if ``spin=True``, a 2x2 matrix,
274 or one per bond).
276 *list_hop* may instead use the neighbour-order form of
277 *System.set_hopping*: dictionaries with keys ('n', 't') and
278 optionally 'ang' and/or 'tag' -- 'n' the neighbour order (1st,
279 2nd, ... shortest distance between orbitals of the infinite
280 lattice), 'ang' a bond angle in degrees, 'tag' a sublattice pair
281 such as ``'ab'``. The bonds are found from *unit_cell* and
282 *prim_vec* (in 1D, 2D and 3D) and turned into the explicit form
283 by *tbkit.neighbours.neighbour_hoppings* (see there for the bond
284 orientations; a negative angle addresses the reversed bonds).
285 A list uses one form or the other, not both.
286 :param hermitian: Boolean. Default value True. If False, the Hermitian
287 conjugates are *not* added: each dictionary is one matrix element
288 :math:`H_{ij}(\mathbf{R})` only, so that non-reciprocal
289 (non-Hermitian) hoppings can be built -- give the reverse hopping
290 explicitly, with its own amplitude.
292 Example usage::
294 # graphene, by neighbour order: nearest (t) and next-nearest
295 # (t2) neighbours, as in System.set_hopping
296 gra.set_hopping([{'n': 1, 't': t}, {'n': 2, 't': t2}])
297 # square lattice, anisotropic nearest neighbours
298 sq.set_hopping([{'n': 1, 'ang': 0., 't': tx}, {'n': 1, 'ang': 90., 't': ty}])
300 # 1D chain, nearest-neighbor hopping t between the only orbital
301 # and its right neighbor:
302 chain.set_hopping([{'i': 0, 'j': 0, 'R': (1,), 't': 1.}])
303 # spinful: spin-independent hopping t, plus a Rashba-like
304 # spin-flip term of strength alpha:
305 chain_spin.set_hopping([{'i': 0, 'j': 0, 'R': (1,),
306 't': t*PAULI['0'] + 1j*alpha*PAULI['y']}])
307 # value function: a staggered hopping, swept without rebuilding
308 chain.set_hopping([{'n': 1, 't': lambda si, sj, t, d: t + d * (-1) ** si['index']}])
309 ham = chain.get_ham((0.,), t=1., d=0.2)
310 '''
311 error_handling.boolean(hermitian, 'hermitian')
312 if error_handling.hopping_form(list_hop):
313 list_hop = neighbours.neighbour_hoppings(self.lat, list_hop, self.spin, hermitian)
314 error_handling.set_hopping_kspace(list_hop, self.n_sites, self.dim, self.spin, values=True)
315 if not hermitian:
316 self._nonreciprocal = True
317 funcs: dict = {} # the bonds of each value function: id -> (function, [(i, j, R_cart)])
318 for dic in list_hop:
319 R_cart = np.zeros(self.space_dim)
320 for n, vec in zip(dic['R'], self.lat.prim_vec):
321 R_cart += n * np.array(vec)
322 i, j, t = dic['i'], dic['j'], dic['t']
323 if callable(t):
324 funcs.setdefault(id(t), (t, []))[1].append((i, j, R_cart))
325 elif self.spin:
326 block = t*PAULI['0'] if np.ndim(t) == 0 else np.asarray(t, 'c16')
327 for a in range(2):
328 for b in range(2):
329 self._hop_const.append((2*i+a, 2*j+b, R_cart, block[a, b]))
330 if hermitian:
331 self._hop_const.append((2*j+b, 2*i+a, -R_cart, np.conj(block[a, b])))
332 else:
333 self._hop_const.append((i, j, R_cart, t))
334 if hermitian:
335 self._hop_const.append((j, i, -R_cart, np.conj(t)))
336 for func, bonds in funcs.values():
337 i, j, R = zip(*bonds)
338 self._hop_values.append((np.array(i), np.array(j), np.array(R), func, hermitian))
340 def set_overlap(self, list_hop: list[dict]) -> None:
341 r'''
342 Set the overlaps :math:`s_{ij}(\mathbf{R}) = \langle i,\mathbf{0}|j,\mathbf{R}\rangle`
343 of a non-orthogonal basis, in the format of *set_hopping* (key 't'
344 for the overlap; the Hermitian conjugates are added):
346 .. math::
348 S(\mathbf{k}) = \mathbb{1} + \sum_{\mathbf{R}} s(\mathbf{R})\,e^{i\mathbf{k}\cdot\mathbf{R}}\, ,
350 after which the bands solve :math:`H(\mathbf{k})v = E\,S(\mathbf{k})v`.
351 (With graphene's nearest-neighbour overlap :math:`s`, the
352 :math:`\pi` bands become :math:`E = (\epsilon \pm t|f|)/(1 \pm s|f|)`:
353 no longer symmetric about :math:`\epsilon`.) The topological tools use
354 the Lowdin-orthonormalized states :math:`S^{1/2}v`.
356 :param list_hop: List of dictionaries ('i', 'j', 'R', 't'), or by
357 neighbour order ('n', 't', optionally 'ang', 'tag'), see
358 *set_hopping* (numbers only, no value functions).
359 '''
360 if error_handling.hopping_form(list_hop):
361 list_hop = neighbours.neighbour_hoppings(self.lat, list_hop, self.spin)
362 error_handling.set_hopping_kspace(list_hop, self.n_sites, self.dim, self.spin)
363 saved = self._hop_const, self._hop_values, self._nonreciprocal
364 self._hop = []
365 self.set_hopping(list_hop)
366 self._overlap_hop += self._hop_const
367 self._hop_const, self._hop_values, self._nonreciprocal = saved
369 def get_overlap(self, k: ArrayLike) -> NDArray[np.complex128]:
370 r'''
371 Get the overlap matrix :math:`S(\mathbf{k})` (the identity for an
372 orthogonal basis), see *set_overlap*.
374 :param k: k point (see *get_ham*).
376 :returns:
377 * **S** -- Complex ndarray, shape (norb, norb).
378 '''
379 error_handling.k_vector(k, 'k', self.dim)
380 k_cart = self.k_basis @ np.asarray(k, dtype='f8')
381 return np.eye(self.norb) + self._bloch_sum(self._overlap_hop, k_cart[None])[0][0]
383 def clear_hopping(self) -> None:
384 '''
385 Clear the hoppings set by *set_hopping* (and the overlaps of *set_overlap*).
386 '''
387 self._hop = []
388 self._nonreciprocal = False
389 self._overlap_hop = []
391 def get_ham(self, k: ArrayLike, **params) -> NDArray[np.complex128]:
392 r'''
393 Get the dense Bloch Hamiltonian :math:`H(\mathbf{k})`.
395 :param k: Tuple/list/ndarray of *dim* real numbers. When the lattice
396 fills its space (2D in the plane, 3D in space), the
397 :math:`\mathbf{k}` point in the same Cartesian frame as *prim_vec*.
398 In 1D, the crystal momentum *along* the primitive vector (so the
399 Brillouin zone spans :math:`2\pi/|\mathbf{a}_1|`), which
400 coincides with :math:`k_x` for a chain aligned with :math:`x`.
401 For a 2D lattice in 3D space (a slab), the components along the
402 orthonormal basis *k_basis* of its plane (the first along
403 :math:`\mathbf{a}_1`). *rec_vec_k* holds the reciprocal vectors in
404 these coordinates.
405 :param params: Values of the parameters of the value functions of
406 *set_hopping*, for this call only: they override those of
407 *set_params*. Each function gets the ones its signature names
408 (all of them with ``**kwargs``); the others are ignored.
410 :returns:
411 * **ham** -- Complex ndarray, shape (norb, norb).
413 Example usage::
415 hams = [ks.get_ham((0., 0.), m=m) for m in np.linspace(-1., 1., 21)]
416 '''
417 error_handling.k_vector(k, 'k', self.dim)
418 return self._bloch_ham((self.k_basis @ np.asarray(k, dtype='f8'))[None],
419 params={**self.params, **params})[0]
421 def get_ham_peierls(self, k: ArrayLike, A: ArrayLike) -> NDArray[np.complex128]:
422 r'''
423 Get the Bloch Hamiltonian in a uniform vector potential
424 :math:`\mathbf{A}` (:math:`\hbar = e = 1`): every hopping picks up the
425 Peierls phase :math:`e^{i\mathbf{A}\cdot\mathbf{d}}` of its bond
426 vector :math:`\mathbf{d} = \mathbf{R} + \boldsymbol\tau_j - \boldsymbol\tau_i`
427 (see *System.set_peierls_phase*). A time-dependent :math:`\mathbf{A}(t)`
428 is a uniform electric field :math:`-\partial_t\mathbf{A}`, e.g. light
429 (see *tbkit.floquet*).
431 :param k: k point (see *get_ham*).
432 :param A: Tuple/ndarray of space_dim real numbers. Vector potential.
434 :returns:
435 * **ham** -- Complex ndarray, shape (norb, norb).
436 '''
437 error_handling.k_vector(k, 'k', self.dim)
438 A = np.asarray(A, dtype='f8')
439 error_handling.ndarray(A, 'A', self.space_dim)
440 return self._bloch_ham((self.k_basis @ np.asarray(k, dtype='f8'))[None], A)[0]
442 def _bloch_ham(
443 self, k_cart: NDArray[np.float64], A: NDArray[np.float64] | None = None, params: dict | None = None,
444 ) -> NDArray[np.complex128]:
445 '''
446 Private method. *get_ham* (or, with a vector potential *A*,
447 *get_ham_peierls*) at many k-points at once: *k_cart* has shape
448 (nk, space_dim) and the result (nk, norb, norb). The value functions
449 take *params* (default: *self.params*).
450 '''
451 hops = self._hops(self.params if params is None else params)
452 return self._bloch_sum(hops, k_cart, A=A)[0] + (np.diag(self.onsite) + self._onsite_offdiag)[None]
454 def get_bands(
455 self, ks: ArrayLike, eigenvec: bool = False,
456 ) -> NDArray[np.float64] | tuple[NDArray[np.float64], NDArray[np.complex128]]:
457 r'''
458 Diagonalize :math:`H(\mathbf{k})` over a set of k-points, and keep
459 the result for *plot_bands* (plotted against the cumulative distance
460 between consecutive k-points).
462 :param ks: ndarray, shape (nk, dim). k-points.
463 :param eigenvec: Boolean. Default value False. If True, also return
464 the eigenvectors.
466 :returns:
467 * **en** -- Real ndarray, shape (nk, norb). Band energies, sorted ascending.
468 Complex, and sorted by real part, if :math:`H(\mathbf{k})` is
469 non-Hermitian (complex onsite energies, see *set_onsite*).
470 * **vn** -- Complex ndarray, shape (nk, norb, norb), only if *eigenvec* is True.
471 vn[k, :, n] is the nth (right) eigenvector at ks[k].
472 '''
473 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
474 error_handling.ks(ks, self.dim)
475 out = self._diagonalize(ks, eigenvec)
476 self.ks = ks
477 self.en = out[0] if isinstance(out, tuple) else out
478 steps = np.linalg.norm(np.diff(ks, axis=0), axis=1)
479 self.ks_dist = np.concatenate([[0.], np.cumsum(steps)])[:len(ks)]
480 self.nodes = np.array([])
481 return out
483 @overload
484 def _diagonalize(
485 self, ks: NDArray[np.float64], eigenvec: Literal[False] = ...,
486 ) -> NDArray[np.float64]: ...
487 @overload
488 def _diagonalize(
489 self, ks: NDArray[np.float64], eigenvec: Literal[True],
490 ) -> tuple[NDArray[np.float64], NDArray[np.complex128]]: ...
491 @overload
492 def _diagonalize(
493 self, ks: NDArray[np.float64], eigenvec: bool,
494 ) -> NDArray[np.float64] | tuple[NDArray[np.float64], NDArray[np.complex128]]: ...
496 def _diagonalize(
497 self, ks: NDArray[np.float64], eigenvec: bool = False,
498 ) -> NDArray[np.float64] | tuple[NDArray[np.float64], NDArray[np.complex128]]:
499 r'''
500 Private method. Diagonalize :math:`H(\mathbf{k})` over the k-points
501 *ks* (shape (nk, dim)), without touching the stored band structure.
502 '''
503 return self._eigs(ks, eigenvec)
505 def set_workers(self, workers: int) -> None:
506 r'''
507 Set the number of threads that diagonalize :math:`H(\mathbf{k})`
508 over many k-points (*get_bands*, *mesh_bands*, *berry_curvature*,
509 the Wilson loops, ...). The k-points are split into chunks, each
510 diagonalized by one stacked LAPACK call, which releases the GIL, so
511 the chunks run in parallel. The results depend on *workers* only up
512 to rounding (and the phases of the eigenvectors).
514 :param workers: Positive integer. Number of threads (default 1, no threads).
516 Example usage::
518 ks.set_workers(4)
519 curv = ks.berry_curvature(0, nk=400)
520 '''
521 error_handling.positive_int(workers, 'workers')
522 self.workers = workers
524 def _hams(self, ks: NDArray[np.float64]) -> NDArray[np.complex128]:
525 '''
526 Private method. *get_ham* at the k-points *ks* (shape (nk, dim)),
527 shape (nk, norb, norb): one Bloch sum, unless a subclass overrides
528 *get_ham* (a Floquet model), which is then called per k-point.
529 '''
530 if type(self).get_ham is KSpace.get_ham:
531 return self._bloch_ham(ks @ self.k_basis.T)
532 return np.array([self.get_ham(k) for k in ks], dtype='c16')
534 def _overlaps(self, ks: NDArray[np.float64]) -> NDArray[np.complex128]:
535 '''
536 Private method. *get_overlap* at the k-points *ks*, shape (nk, norb, norb).
537 '''
538 return np.eye(self.norb) + self._bloch_sum(self._overlap_hop, ks @ self.k_basis.T)[0]
540 @overload
541 def _eigs(self, ks: NDArray[np.float64], eigenvec: Literal[False] = ...) -> NDArray: ...
542 @overload
543 def _eigs(
544 self, ks: NDArray[np.float64], eigenvec: Literal[True],
545 ) -> tuple[NDArray, NDArray[np.complex128]]: ...
546 @overload
547 def _eigs(
548 self, ks: NDArray[np.float64], eigenvec: bool,
549 ) -> NDArray | tuple[NDArray, NDArray[np.complex128]]: ...
551 def _eigs(
552 self, ks: NDArray[np.float64], eigenvec: bool = False,
553 ) -> NDArray | tuple[NDArray, NDArray[np.complex128]]:
554 '''
555 Private method. *_eig* at the k-points *ks* (shape (nk, dim)): the
556 eigenvalues, shape (nk, norb), and the eigenvectors, shape
557 (nk, norb, norb). The k-points go in chunks of bounded memory,
558 spread over *workers* threads (see *set_workers*).
559 '''
560 chunk = max(1, min(400000 // self.norb ** 2, -(-len(ks) // self.workers)))
561 parts = [ks[c0:c0 + chunk] for c0 in range(0, len(ks), chunk)] or [ks]
562 if self.workers > 1 and len(parts) > 1:
563 with ThreadPoolExecutor(self.workers) as pool:
564 out = list(pool.map(lambda part: self._eig_chunk(part, eigenvec), parts))
565 else:
566 out = [self._eig_chunk(part, eigenvec) for part in parts]
567 if eigenvec:
568 return np.concatenate([o[0] for o in out]), np.concatenate([o[1] for o in out])
569 return np.concatenate(out)
571 def _eig_chunk(
572 self, ks: NDArray[np.float64], eigenvec: bool,
573 ) -> NDArray | tuple[NDArray, NDArray[np.complex128]]:
574 r'''
575 Private method. One chunk of *_eigs*: the Hermitian solver when
576 *is_hermitian*, otherwise the general one, with the eigenvalues
577 sorted by real part, both stacked over the k-points. With an overlap
578 (see *set_overlap*), the generalized problem :math:`Hv = ES(\mathbf{k})v`:
579 the Hermitian one is reduced by the Cholesky factor
580 :math:`S = LL^\dagger` to :math:`L^{-1}HL^{-\dagger}` (eigenvectors
581 :math:`L^{-\dagger}v`, normalized to :math:`v^\dagger Sv = 1` as by
582 *scipy.linalg.eigh*); the non-Hermitian one has no stacked solver,
583 and is solved one k-point at a time.
584 '''
585 ham = self._hams(ks)
586 hermitian = self.is_hermitian()
587 en: NDArray
588 vn: NDArray
589 if self._overlap_hop and hermitian:
590 l_inv = np.linalg.inv(np.linalg.cholesky(self._overlaps(ks)))
591 l_inv_h = l_inv.conj().transpose(0, 2, 1)
592 ham = l_inv @ ham @ l_inv_h
593 if hermitian:
594 if not eigenvec:
595 return np.linalg.eigvalsh(ham)
596 en, vn = np.linalg.eigh(ham)
597 return en, (l_inv_h @ vn if self._overlap_hop else vn)
598 if self._overlap_hop:
599 solved = [LA.eig(h, s) for h, s in zip(ham, self._overlaps(ks))]
600 en = np.array([w for w, _ in solved], dtype='c16').reshape(len(ks), self.norb)
601 vn = np.array([v for _, v in solved], dtype='c16').reshape(ham.shape)
602 elif eigenvec:
603 en, vn = np.linalg.eig(ham)
604 else:
605 en = np.linalg.eigvals(ham)
606 ind = np.argsort(en.real, axis=1, kind='stable')
607 if not eigenvec:
608 return np.take_along_axis(en, ind, axis=1)
609 return np.take_along_axis(en, ind, axis=1), np.take_along_axis(vn, ind[:, None, :], axis=2)
611 def is_hermitian(self) -> bool:
612 '''
613 Check whether :math:`H(\\mathbf{k})` is Hermitian: it is, unless an
614 onsite term is not, or *set_hopping* was called with
615 ``hermitian=False`` (such a model is always treated as
616 non-Hermitian, and diagonalized with the general solver).
618 :returns:
619 * **hermitian** -- Boolean.
620 '''
621 return bool(not self._nonreciprocal and np.all(self.onsite.imag == 0.)
622 and np.array_equal(self._onsite_offdiag, self._onsite_offdiag.conj().T))
624 def k_path(
625 self, points: list[ArrayLike], nk: int,
626 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
627 r'''
628 Build a k-path through a list of high-symmetry points, and get the
629 associated bands.
631 :param points: List of at least two k-points (each a tuple/list of
632 *dim* real numbers), e.g. from *high_symmetry_path*.
633 :param nk: Positive integer. Number of k-points per path segment.
635 :returns:
636 * **ks_dist** -- Real ndarray. Cumulative distance along the path,
637 to be used as the x-axis of a band-structure plot.
638 * **en** -- Real ndarray, shape (len(ks_dist), norb). Band energies.
639 '''
640 error_handling.k_path_points(points, self.dim)
641 error_handling.positive_int(nk, 'nk')
642 ks, _, nodes = _path(points, nk)
643 self.get_bands(ks)
644 self.nodes = np.array(nodes)
645 return self.ks_dist, self.en
647 def mesh_grid(
648 self, nk: int | tuple[int, ...], monkhorst_pack: bool = False,
649 ) -> tuple[list[NDArray[np.float64]], NDArray[np.float64]]:
650 r'''
651 Build a uniform grid of fractional coordinates spanning the
652 Brillouin zone, and the corresponding Cartesian k-points. By
653 default the grid contains :math:`\Gamma`: :math:`r/q`,
654 :math:`r = 0, \dots, q-1` (each in [0, 1)). With *monkhorst_pack*,
655 it is the grid of Monkhorst and Pack (Phys. Rev. B 13, 5188 (1976)),
656 :math:`(2r - q - 1)/2q`, :math:`r = 1, \dots, q` (each in
657 (-1/2, 1/2)), symmetric about :math:`\Gamma`: the same points for
658 odd :math:`q`, shifted by half a step (off :math:`\Gamma` and the
659 zone boundary) for even :math:`q`. Both integrate exactly every
660 Fourier component :math:`e^{i\mathbf{k}\cdot\mathbf{R}}` with
661 :math:`|R_i| < q`.
663 :param nk: Positive integer, or tuple of *dim* positive integers.
664 Number of k-points along each reciprocal lattice vector.
665 :param monkhorst_pack: Boolean. Default value False. If True, the
666 Monkhorst-Pack grid.
668 :returns:
669 * **fracs** -- List of *dim* real ndarrays, shape (nk1, nk2) each
670 (or (nk1,) in 1D): fractional coordinates of the grid.
671 * **ks** -- Real ndarray, shape (nk1*nk2, dim) (or (nk1, dim) in 1D).
672 '''
673 error_handling.nk(nk, self.dim)
674 error_handling.boolean(monkhorst_pack, 'monkhorst_pack')
675 if isinstance(nk, int):
676 nk = (nk,) * self.dim
677 if monkhorst_pack:
678 axes = [(2 * np.arange(1, n + 1) - n - 1) / (2 * n) for n in nk]
679 else:
680 axes = [np.arange(n) / n for n in nk]
681 fracs = np.meshgrid(*axes, indexing='ij')
682 ks = sum(f.ravel()[:, None] * b[None, :] for f, b in zip(fracs, self.rec_vec_k))
683 if self.dim == 1:
684 return [fracs[0]], ks
685 return list(fracs), ks
687 def mesh_bands(self, nk: int | tuple[int, ...], monkhorst_pack: bool = False) -> NDArray[np.float64]:
688 '''
689 Diagonalize :math:`H(\\mathbf{k})` over a uniform mesh spanning the
690 Brillouin zone. Unlike *get_bands*, it leaves the band structure kept
691 for *plot_bands* untouched.
693 :param nk: Positive integer, or tuple of *dim* positive integers.
694 Number of k-points along each reciprocal lattice vector.
695 :param monkhorst_pack: Boolean. Default value False. If True, the
696 Monkhorst-Pack grid (see *mesh_grid*).
698 :returns:
699 * **en** -- Real ndarray, shape (nk1*nk2, norb). Band energies
700 over the mesh (flattened).
701 '''
702 _, ks = self.mesh_grid(nk, monkhorst_pack)
703 return self._diagonalize(ks)
705 def get_fermi_level(
706 self, n_electrons: float, nk: int | tuple[int, ...] = 30, temperature: float = 0.,
707 ) -> float:
708 '''
709 Get the Fermi level holding *n_electrons* electrons per unit cell
710 (one electron per band and k-point), from the bands over a uniform
711 Brillouin-zone mesh (see *tbkit.occupation.fermi_level*).
713 :param n_electrons: Positive real or zero, at most norb. Electrons per unit cell.
714 :param nk: Positive integer, or tuple. Default value 30. k-mesh.
715 :param temperature: Positive real or zero. Default value 0.
717 :returns:
718 * **e_fermi** -- Real number.
719 '''
720 en = self.mesh_bands(nk)
721 return occupation.fermi_level(en, n_electrons, temperature,
722 weights=np.full(en.shape, 1. / len(en)))
724 def berry_curvature(
725 self, bands: int | list[int], nk: int | tuple[int, int] = 30,
726 plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
727 ) -> NDArray[np.float64]:
728 r'''
729 Get the Berry curvature of a group of bands over a uniform
730 Brillouin-zone mesh, using the gauge-invariant lattice method of
731 Fukui, Hatsugai and Suzuki (J. Phys. Soc. Jpn. 74, 1674 (2005)):
732 the flux through each mesh plaquette is minus the phase of the
733 product of the (Slater-determinant) overlaps between the occupied
734 subspaces at its four corners,
736 .. math::
738 \Omega_{\square} = -\arg\prod_{\langle\mathbf{k}\mathbf{k}'\rangle\in\partial\square}
739 \det\langle u(\mathbf{k})|u(\mathbf{k}')\rangle\, ,
741 the corners being visited counterclockwise. This is the discrete
742 version of :math:`\Omega = \partial_{k_x}A_y - \partial_{k_y}A_x`,
743 :math:`\mathbf{A} = i\langle u|\nabla_{\mathbf{k}}u\rangle`, and does
744 not depend on the orientation (handedness) of *prim_vec*. The local
745 curvature is that of the periodic gauge of *get_ham* (orbitals on the
746 cell origin), which differs point by point from the one with the
747 orbital positions included (see *quantum_geometric_tensor*); their
748 integral, the Chern number, is the same. For a
749 non-Hermitian :math:`H(\mathbf{k})` (see *set_onsite*), the bands
750 are ordered by the real part of their energy and the overlaps are
751 taken between right eigenvectors (the "RR" curvature of
752 *biorthogonal_berry_curvature*). That ordering follows the bands
753 continuously only if a real line gap separates *bands* from the
754 others over the whole zone: without one (an imaginary line gap,
755 or exceptional points), a ValueError is raised instead of a
756 meaningless result -- see *biorthogonal_berry_curvature*.
758 :param bands: Positive integer, or list of positive integers. Band
759 index, or indices of a group of bands (e.g. all occupied bands
760 below a gap).
761 :param nk: Positive integer, or tuple of 2 positive integers.
762 Default value 30. Number of k-points along each of the two
763 reciprocal lattice vectors spanning the plane.
764 :param plane: Tuple of two distinct integers. Default value (0, 1).
765 3D lattices only: the reciprocal vectors
766 :math:`\mathbf{b}_i, \mathbf{b}_j` spanning the plane of the
767 Brillouin zone. The flux is oriented along the remaining
768 primitive vector :math:`\mathbf{a}_l` (along +z in 2D).
769 :param k_fixed: Real number. Default value 0. 3D lattices only: the
770 fractional coordinate of the plane along the remaining reciprocal
771 vector :math:`\mathbf{b}_l` (e.g. :math:`k_z` in units of
772 :math:`|\mathbf{b}_3|` for ``plane=(0, 1)``).
774 :returns:
775 * **curv** -- Real ndarray, shape (nk1, nk2). Berry curvature
776 (flux through each plaquette, in radians). Summing *curv* and
777 dividing by :math:`2\pi` gives the Chern number, see
778 *chern_number*.
779 '''
780 error_handling.dim_min(self.dim, 2)
781 if isinstance(bands, int):
782 bands = [bands]
783 error_handling.band_indices(bands, self.norb)
784 if isinstance(nk, int):
785 nk = (nk, nk)
786 error_handling.nk(nk, 2)
787 error_handling.plane(plane, self.dim)
788 error_handling.real_number(k_fixed, 'k_fixed')
789 n1, n2 = nk
790 ks = self._plane_mesh(nk, plane, k_fixed)
791 if not self.is_hermitian():
792 # bands labelled by Re E are continuous only across a real line gap
793 en = self._eigs(ks.reshape(-1, self.dim))
794 error_handling.line_gap(en.real, bands, 'real')
795 error_handling.band_continuity(en.reshape(n1, n2, self.norb), bands)
796 v1 = self._subspaces(ks.reshape(-1, self.dim), bands).reshape(n1, n2, self.norb, len(bands))
797 return self._flux(v1, plane)
799 def _flux(self, v1: NDArray[np.complex128], plane: tuple[int, int]) -> NDArray[np.float64]:
800 '''
801 Private method. Fukui-Hatsugai-Suzuki flux through each plaquette of
802 the mesh of *_plane_mesh*, from orthonormal bases *v1* of a subspace
803 at its nodes (shape (n1, n2, norb, n)), oriented as in *berry_curvature*.
804 '''
805 # the other corners of each plaquette: k + b_i/n1, k + b_i/n1 + b_j/n2, k + b_j/n2
806 v2 = np.roll(v1, -1, axis=0)
807 v3 = np.roll(v2, -1, axis=1)
808 v4 = np.roll(v1, -1, axis=1)
810 def overlap(a, b):
811 return np.linalg.det(a.conj().swapaxes(-1, -2) @ b)
812 link = overlap(v1, v2) * overlap(v2, v3) * overlap(v3, v4) * overlap(v4, v1)
813 curv = -np.angle(link)
814 # The corners k, k+b_i/n1, k+b_i/n1+b_j/n2, k+b_j/n2 run
815 # counterclockwise about the normal only if (b_i x b_j).normal > 0;
816 # otherwise every flux comes out with the opposite sign.
817 return self._plane_orientation(plane) * curv
819 def _plane_mesh(
820 self, nk: tuple[int, int], plane: tuple[int, int], k_fixed: float,
821 ) -> NDArray[np.float64]:
822 '''
823 Private method. k-points, shape (n1, n2, dim), of the uniform mesh of
824 the Brillouin-zone plane spanned by b_i, b_j (plane = (i, j)), at
825 fractional coordinate k_fixed along the remaining reciprocal vector.
826 '''
827 n1, n2 = nk
828 i, j = plane
829 rest = [l for l in range(self.dim) if l not in plane]
830 fracs = np.zeros((n1, n2, self.dim))
831 fracs[:, :, i] = (np.arange(n1) / n1)[:, None]
832 fracs[:, :, j] = (np.arange(n2) / n2)[None, :]
833 for l in rest:
834 fracs[:, :, l] = k_fixed
835 return fracs @ self.rec_vec_k
837 def _plane_orientation(self, plane: tuple[int, int]) -> float:
838 '''
839 Private method. +1 if (b_i, b_j) is right-handed about the plane's
840 normal (+z in 2D; the remaining primitive vector in 3D), else -1.
841 '''
842 b = [np.zeros(3) for _ in range(2)]
843 for bb, idx in zip(b, plane):
844 bb[:self.space_dim] = self.rec_vec[idx]
845 if self.dim == 3:
846 normal = np.array(self.lat.prim_vec[3 - sum(plane)], dtype='f8')
847 else:
848 a = [np.zeros(3), np.zeros(3)]
849 for aa, vec in zip(a, self.lat.prim_vec):
850 aa[:self.space_dim] = vec
851 normal = np.cross(a[0], a[1])
852 if abs(normal[2]) > 1e-12:
853 normal = np.array([0., 0., 1.])
854 return float(np.sign(np.dot(np.cross(b[0], b[1]), normal)))
856 def _subspace(self, k: NDArray[np.float64], bands: list[int]) -> NDArray[np.complex128]:
857 '''
858 Private method. Orthonormal basis (columns) of the eigenvectors of
859 *bands* at *k* (in the periodic gauge, orbital positions ignored).
860 With an overlap S(k), the Lowdin-orthonormalized S^(1/2) v.
861 '''
862 return self._subspaces(np.asarray(k, dtype='f8')[None], bands)[0]
864 def _subspaces(self, ks: NDArray[np.float64], bands: list[int]) -> NDArray[np.complex128]:
865 '''
866 Private method. *_subspace* at the k-points *ks* (shape (nk, dim)),
867 shape (nk, norb, len(bands)). S^(1/2) comes from the eigenvalues of
868 the (Hermitian, positive) overlap.
869 '''
870 _, vn = self._eigs(ks, eigenvec=True)
871 if self._overlap_hop:
872 sv, us = np.linalg.eigh(self._overlaps(ks))
873 vn = (us * np.sqrt(sv)[:, None, :]) @ (us.conj().transpose(0, 2, 1) @ vn)
874 return vn[:, :, bands]
876 def biorthogonal_berry_curvature(
877 self, bands: int | list[int], nk: int | tuple[int, int] = 30, kind: str = 'LR',
878 gap: str = 'real', plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
879 ) -> NDArray[np.float64]:
880 r'''
881 Get the Berry curvature of a group of bands of a (non-Hermitian)
882 model from its right and left eigenvectors, :math:`H|R_n\rangle =
883 E_n|R_n\rangle` and :math:`H^\dagger|L_n\rangle = E_n^*|L_n\rangle`,
884 normalized by :math:`\langle L_m|R_n\rangle = \delta_{mn}` (Shen, Zhen
885 and Fu, Phys. Rev. Lett. 120, 146402 (2018)). The four connections
887 .. math::
889 \mathbf{A}^{\alpha\beta} = i\langle u^\alpha|\nabla_{\mathbf{k}}u^\beta\rangle\, ,
890 \qquad \alpha, \beta \in \{L, R\}\, ,
892 have different curvatures point by point, but their integrals, the
893 Chern numbers :math:`C^{LR} = C^{RL} = C^{RR} = C^{LL}`, are the same
894 integer. The flux through each plaquette is computed as in
895 *berry_curvature* (Fukui-Hatsugai-Suzuki), with the links
896 :math:`\det\langle u^\alpha(\mathbf{k})|u^\beta(\mathbf{k}')\rangle` (for
897 RR and LL, the phase does not depend on the normalization). For a
898 Hermitian model all four are the curvature of *berry_curvature*.
900 The bands are ordered at every k-point by the real part of their
901 energy (*gap* = 'real') or by its imaginary part ('imaginary'), and
902 *bands* must be separated from the others by a line gap of that kind
903 over the whole zone -- a line :math:`\mathrm{Re}\,E = E_0` (or
904 :math:`\mathrm{Im}\,E = E_0`) that no band crosses. Otherwise a
905 ValueError is raised: the gap closes at exceptional points (see
906 *tbkit.exceptional.find_exceptional_points*), where the Chern number
907 can change.
909 :param bands: Positive integer, or list of positive integers. Band
910 indices in the order set by *gap*.
911 :param nk: Positive integer, or tuple of 2 positive integers.
912 Default value 30. See *berry_curvature*.
913 :param kind: 'LR', 'RL', 'RR' or 'LL'. Default value 'LR'.
914 :math:`\alpha\beta` above.
915 :param gap: 'real' or 'imaginary'. Default value 'real'. Kind of line
916 gap, and the order of the bands.
917 :param plane: Tuple of two integers. Default value (0, 1). 3D only, see
918 *berry_curvature*.
919 :param k_fixed: Real number. Default value 0. 3D only, see *berry_curvature*.
921 :returns:
922 * **curv** -- Real ndarray, shape (nk1, nk2). Flux through each
923 plaquette, in radians (sign convention of *berry_curvature*).
924 '''
925 error_handling.dim_min(self.dim, 2)
926 if isinstance(bands, int):
927 bands = [bands]
928 error_handling.band_indices(bands, self.norb)
929 if isinstance(nk, int):
930 nk = (nk, nk)
931 error_handling.nk(nk, 2)
932 error_handling.biorthogonal_kind(kind)
933 error_handling.gap_kind(gap)
934 error_handling.plane(plane, self.dim)
935 error_handling.real_number(k_fixed, 'k_fixed')
936 error_handling.no_overlap(self._overlap_hop)
937 n1, n2 = nk
938 ks = self._plane_mesh(nk, plane, k_fixed)
939 right = np.zeros((n1, n2, self.norb, len(bands)), 'c16')
940 left = np.zeros_like(right)
941 en = np.zeros((n1, n2, self.norb), 'c16')
942 for i1 in range(n1):
943 for i2 in range(n2):
944 w, vl, vr = LA.eig(self.get_ham(ks[i1, i2]), left=True, right=True)
945 order = np.argsort(w.real if gap == 'real' else w.imag, kind='stable')
946 en[i1, i2] = w[order]
947 r, l = vr[:, order][:, bands], vl[:, order][:, bands]
948 # biorthonormal: <L_m|R_n> = delta_mn
949 right[i1, i2] = r
950 left[i1, i2] = l @ np.linalg.inv(l.conj().T @ r).conj().T
951 error_handling.line_gap(en.reshape(-1, self.norb).real if gap == 'real'
952 else en.reshape(-1, self.norb).imag, bands, gap)
953 error_handling.band_continuity(en, bands)
954 bra, ket = {'LR': (left, right), 'RL': (right, left),
955 'RR': (right, right), 'LL': (left, left)}[kind]
956 # one link per oriented edge, U = det<u^alpha(k)|u^beta(k + dk)>, and
957 # its inverse when the edge is run backwards: <L_a|R_b> and <L_b|R_a>
958 # are not conjugate, so using both would spoil the exact
959 # cancellation of the shared edges (for RR, 1/U has the phase of U*).
960 link = np.zeros((2, n1, n2), 'c16')
961 for i1 in range(n1):
962 for i2 in range(n2):
963 link[0, i1, i2] = np.linalg.det(bra[i1, i2].conj().T @ ket[(i1 + 1) % n1, i2])
964 link[1, i1, i2] = np.linalg.det(bra[i1, i2].conj().T @ ket[i1, (i2 + 1) % n2])
965 flux = (link[0] * np.roll(link[1], -1, axis=0)
966 / (np.roll(link[0], -1, axis=1) * link[1]))
967 return self._plane_orientation(plane) * -np.angle(flux)
969 def biorthogonal_chern_number(
970 self, bands: int | list[int], nk: int | tuple[int, int] = 30, kind: str = 'LR',
971 gap: str = 'real', plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
972 ) -> float:
973 r'''
974 Get the Chern number of a line-gapped group of bands of a
975 (non-Hermitian) model, :math:`C^{\alpha\beta} = \frac{1}{2\pi}\int_{BZ}
976 \Omega^{\alpha\beta}\,d^2k`, the same integer for the four kinds
977 LR, RL, RR, LL (see *biorthogonal_berry_curvature*). It equals the
978 Hermitian Chern number of *chern_number* as long as the line gap stays
979 open.
981 :param bands: See *biorthogonal_berry_curvature*.
982 :param nk: See *biorthogonal_berry_curvature*.
983 :param kind: 'LR', 'RL', 'RR' or 'LL'. Default value 'LR'.
984 :param gap: 'real' or 'imaginary'. Default value 'real'.
985 :param plane: See *berry_curvature*.
986 :param k_fixed: See *berry_curvature*.
988 :returns:
989 * **chern** -- Real number, close to an integer.
990 '''
991 return self.biorthogonal_berry_curvature(bands, nk, kind, gap, plane, k_fixed).sum() / (2*PI)
993 def chern_number(
994 self, bands: int | list[int], nk: int | tuple[int, int] = 30,
995 plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
996 ) -> float:
997 r'''
998 Get the Chern number of a group of bands:
1000 .. math::
1002 C = \frac{1}{2\pi}\int_{BZ} \Omega(\mathbf{k})\, d^2k
1004 an integer (up to the numerical precision set by *nk*) for a group
1005 of bands that is isolated from the rest of the spectrum by a gap
1006 everywhere in the Brillouin zone. See *berry_curvature*. For a
1007 non-Hermitian model, *bands* (ordered by Re E) must be isolated by a
1008 real line gap, or a ValueError is raised; *biorthogonal_chern_number*
1009 also handles imaginary line gaps and the left eigenvectors.
1011 :param bands: Positive integer, or list of positive integers. Band
1012 index, or indices of a group of bands (e.g. all occupied bands
1013 below a gap).
1014 :param nk: Positive integer, or tuple of 2 positive integers.
1015 Default value 30. Number of k-points along each reciprocal
1016 lattice vector.
1017 :param plane: Tuple of two integers. Default value (0, 1). 3D lattices
1018 only, see *berry_curvature*.
1019 :param k_fixed: Real number. Default value 0. 3D lattices only, see
1020 *berry_curvature*: the Chern number of the plane at this
1021 fractional coordinate (it jumps across a Weyl point).
1023 :returns:
1024 * **chern** -- Real number, close to an integer.
1025 '''
1026 return self.berry_curvature(bands, nk, plane, k_fixed).sum() / (2*PI)
1028 def sector_chern_numbers(
1029 self, op, bands: int | list[int], nk: int | tuple[int, int] = 30,
1030 plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
1031 ) -> tuple[float, float]:
1032 r'''
1033 Get the Chern numbers of the two sectors into which a Hermitian
1034 operator :math:`O` splits a group of bands: at every k-point, the
1035 projected operator :math:`P(\mathbf{k})\,O\,P(\mathbf{k})`
1036 (:math:`P` the projector on *bands*) is diagonalized within the
1037 bands, and its eigenvectors with positive and negative eigenvalues
1038 span the two sectors, whose Chern numbers :math:`C_\pm` follow as in
1039 *berry_curvature*. They are well defined as long as the spectrum of
1040 :math:`POP` keeps a gap around zero over the whole plane, even if
1041 :math:`O` does not commute with :math:`H` (Prodan, Phys. Rev. B 80,
1042 125327 (2009)); otherwise a ValueError is raised. See
1043 *spin_chern_number* and *mirror_chern_number*.
1045 :param op: Complex ndarray, shape (norb, norb), or callable of k
1046 returning one: the Hermitian operator :math:`O`.
1047 :param bands: Band index, or list of band indices.
1048 :param nk: Positive integer, or tuple of 2 positive integers. Default value 30.
1049 :param plane: Tuple of two integers. Default value (0, 1). 3D only, see *berry_curvature*.
1050 :param k_fixed: Real number. Default value 0. 3D only, see *berry_curvature*.
1052 :returns:
1053 * **c_plus** -- Real number, close to an integer: Chern number
1054 of the sector with :math:`POP > 0`.
1055 * **c_minus** -- Real number: that of the sector with :math:`POP < 0`.
1056 :math:`C_+ + C_-` is the Chern number of *bands*.
1057 '''
1058 bands, nk, ks = self._check_sectors(bands, nk, plane, k_fixed)
1059 if callable(op):
1060 ops = np.array([self._operator(op, k) for k in ks])
1061 else:
1062 ops = self._operator(op, ks[0])[None]
1063 error_handling.hermitian_operator(ops)
1064 return self._sector_cherns(ops, bands, nk, ks, plane)
1066 def _check_sectors(
1067 self, bands, nk, plane, k_fixed,
1068 ) -> tuple[list[int], tuple[int, int], NDArray[np.float64]]:
1069 '''
1070 Private method. Validate the arguments of the sector Chern numbers;
1071 return them with the k-points of the plane mesh, shape (n1 n2, dim).
1072 '''
1073 error_handling.dim_min(self.dim, 2)
1074 if isinstance(bands, int):
1075 bands = [bands]
1076 error_handling.band_indices(bands, self.norb)
1077 if isinstance(nk, int):
1078 nk = (nk, nk)
1079 error_handling.nk(nk, 2)
1080 error_handling.plane(plane, self.dim)
1081 error_handling.real_number(k_fixed, 'k_fixed')
1082 error_handling.hermitian_kspace(self.is_hermitian(), self._overlap_hop)
1083 return bands, nk, self._plane_mesh(nk, plane, k_fixed).reshape(-1, self.dim)
1085 def _sector_cherns(
1086 self, ops: NDArray[np.complex128], bands: list[int], nk: tuple[int, int],
1087 ks: NDArray[np.float64], plane: tuple[int, int],
1088 ) -> tuple[float, float]:
1089 '''
1090 Private method. *sector_chern_numbers* with the Hermitian operator
1091 evaluated at the k-points *ks* (shape (n1 n2, norb, norb), or
1092 (1, norb, norb) if constant).
1093 '''
1094 n1, n2 = nk
1095 v = self._subspaces(ks, bands)
1096 w, u = np.linalg.eigh(v.conj().transpose(0, 2, 1) @ ops @ v)
1097 n_minus = error_handling.projected_gap(w)
1098 sectors = v @ u
1099 c_minus = self._flux(sectors[:, :, :n_minus].reshape(n1, n2, self.norb, -1), plane)
1100 c_plus = self._flux(sectors[:, :, n_minus:].reshape(n1, n2, self.norb, -1), plane)
1101 return float(c_plus.sum() / (2*PI)), float(c_minus.sum() / (2*PI))
1103 def spin_chern_number(
1104 self, bands: list[int], nk: int | tuple[int, int] = 30, s_z=None,
1105 plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
1106 ) -> float:
1107 r'''
1108 Get the spin Chern number of a group of bands (Prodan, Phys. Rev. B
1109 80, 125327 (2009); Sheng, Weng, Sheng and Haldane, Phys. Rev. Lett.
1110 97, 036808 (2006)),
1112 .. math::
1114 C_s = \frac{C_+ - C_-}{2}\, ,
1116 with :math:`C_\pm` the Chern numbers of the two sectors of the
1117 projected spin operator :math:`Ps_zP` (see *sector_chern_numbers*).
1118 When :math:`s_z` is conserved they are the Chern numbers of the two
1119 spins; Rashba coupling mixes the spins, but :math:`C_s` stays
1120 quantized, and unchanged, as long as the spectrum of :math:`Ps_zP`
1121 keeps a gap around zero. With time reversal, :math:`C_s` modulo 2 is
1122 the :math:`\mathbb{Z}_2` invariant of *z2_invariant*.
1124 :param bands: List of band indices (e.g. the occupied bands).
1125 :param nk: Positive integer, or tuple of 2 positive integers. Default value 30.
1126 :param s_z: Complex ndarray, shape (norb, norb), or callable of k.
1127 Default value None: :math:`\sigma_z` on every site of a model
1128 built with ``spin=True``. The spin operator.
1129 :param plane: Tuple of two integers. Default value (0, 1). 3D only, see *berry_curvature*.
1130 :param k_fixed: Real number. Default value 0. 3D only, see *berry_curvature*.
1132 :returns:
1133 * **c_s** -- Real number, close to an integer.
1134 '''
1135 if s_z is None:
1136 error_handling.spinful(self.spin)
1137 s_z = np.kron(np.eye(self.norb // 2), PAULI['z'])
1138 c_plus, c_minus = self.sector_chern_numbers(s_z, bands, nk, plane, k_fixed)
1139 return (c_plus - c_minus) / 2
1141 def mirror_chern_number(
1142 self, mirror, bands: list[int], nk: int | tuple[int, int] = 30,
1143 plane: tuple[int, int] = (0, 1), k_fixed: float = 0., tol: float = 1e-8,
1144 ) -> float:
1145 r'''
1146 Get the mirror Chern number of a group of bands on a mirror-invariant
1147 plane of the Brillouin zone (Teo, Fu and Kane, Phys. Rev. B 78,
1148 045426 (2008); Hsieh et al., Nat. Commun. 3, 982 (2012)),
1150 .. math::
1152 n_M = \frac{C_{+} - C_{-}}{2}\, ,
1154 with :math:`C_\pm` the Chern numbers of the bands of mirror
1155 eigenvalue :math:`+i` and :math:`-i` (:math:`M^2 = -1`, spinful
1156 electrons) or :math:`+1` and :math:`-1` (:math:`M^2 = +1`). The
1157 mirror must commute with :math:`H(\mathbf{k})` at every point of the
1158 plane: in 3D, the plane fixed by the reflection (e.g.
1159 :math:`k_z = 0` or :math:`\pi` for :math:`z\to-z`); in 2D, the whole
1160 zone, for a reflection of the third axis. A topological crystalline
1161 insulator such as SnTe has :math:`n_M = 2` and every
1162 :math:`\mathbb{Z}_2` index zero (:math:`\nu = n_M` modulo 2 on the plane).
1164 :param mirror: Complex ndarray, shape (norb, norb), or callable of k:
1165 the unitary mirror operator :math:`M`, with :math:`M^2 = \pm1`.
1166 :param bands: List of band indices (e.g. the occupied bands).
1167 :param nk: Positive integer, or tuple of 2 positive integers. Default value 30.
1168 :param plane: Tuple of two integers. Default value (0, 1). 3D only, see *berry_curvature*.
1169 :param k_fixed: Real number. Default value 0. 3D only: the
1170 mirror-invariant plane (0 or 0.5 for a reflection of the remaining axis).
1171 :param tol: Positive real. Default value 1e-8. Tolerance of the
1172 checks :math:`M^2 = \pm1` and :math:`[M, H(\mathbf{k})] = 0`.
1174 :returns:
1175 * **n_m** -- Real number, close to an integer.
1176 '''
1177 bands, nk, ks = self._check_sectors(bands, nk, plane, k_fixed)
1178 error_handling.positive_real(tol, 'tol')
1179 if callable(mirror):
1180 mats = np.array([self._operator(mirror, k) for k in ks])
1181 else:
1182 mats = self._operator(mirror, ks[0])[None]
1183 sign = error_handling.mirror_square(mats @ mats, tol)
1184 hams = self._hams(ks)
1185 error_handling.commutes(float(np.max(np.abs(mats @ hams - hams @ mats))), tol)
1186 # Hermitian, with eigenvalue +1 on the sector M = +i (M = +1 if M^2 = 1)
1187 herm = mats if sign > 0 else -1j * mats
1188 c_plus, c_minus = self._sector_cherns(herm, bands, nk, ks, plane)
1189 return (c_plus - c_minus) / 2
1191 # ------------------------------------------------------------------
1192 # Anomalous and spin Hall conductivities (Kubo formula)
1193 # ------------------------------------------------------------------
1195 def _check_static(self) -> None:
1196 '''
1197 Private method. Methods that differentiate the stored hoppings call
1198 it first (a Floquet model overrides it to refuse).
1199 '''
1201 def _bloch_sum(
1202 self, hops: list, k_cart: NDArray[np.float64], directions: Sequence[NDArray[np.float64]] = (),
1203 positions: bool = False, A: NDArray[np.float64] | None = None,
1204 ) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
1205 r'''
1206 Private method. The Bloch sum :math:`\sum t\,e^{i\mathbf{k}\cdot\mathbf{d}}`
1207 of a hopping list over many k-points at once (*k_cart*, shape
1208 (nk, space_dim)), and its derivatives
1209 :math:`\sum i(\mathbf{u}\cdot\mathbf{d})\,t\,e^{i\mathbf{k}\cdot\mathbf{d}}`
1210 along the Cartesian unit vectors *directions*, with
1211 :math:`\mathbf{d} = \mathbf{R}` (plus :math:`\boldsymbol\tau_j-\boldsymbol\tau_i`
1212 with *positions*). A vector potential *A* multiplies each term by the
1213 Peierls phase :math:`e^{i\mathbf{A}\cdot(\mathbf{R}+\boldsymbol\tau_j-\boldsymbol\tau_i)}`
1214 (see *get_ham_peierls*). Shapes (nk, norb, norb) and
1215 (len(directions), nk, norb, norb). Every Bloch matrix of the model
1216 (*get_ham*, *get_overlap*, *get_ham_peierls*, *_bloch_derivatives*)
1217 is built here.
1218 '''
1219 n, nk = self.norb, len(k_cart)
1220 i = np.array([h[0] for h in hops], dtype=int)
1221 j = np.array([h[1] for h in hops], dtype=int)
1222 d = np.array([h[2] for h in hops], dtype='f8').reshape(len(hops), self.space_dim)
1223 t = np.array([h[3] for h in hops], dtype='c16')
1224 if positions or A is not None:
1225 tau = self.orbital_positions()
1226 bond = d + tau[j] - tau[i]
1227 if positions:
1228 d = bond
1229 phase = k_cart @ d.T
1230 if A is not None:
1231 phase = phase + bond @ A
1232 w = np.exp(1j * phase) * t[None, :]
1233 # accumulate term m of k-point p in entry (p, i_m, j_m): repeated
1234 # (i, j) pairs add up, as they must
1235 index = (np.arange(nk)[:, None] * n * n + i * n + j).ravel()
1237 def scatter(v):
1238 size = nk * n * n
1239 return (np.bincount(index, v.real.ravel(), size)
1240 + 1j * np.bincount(index, v.imag.ravel(), size)).reshape(nk, n, n)
1241 mat = scatter(w)
1242 dmat = np.array([scatter(w * (1j * (d @ u))[None, :]) for u in directions],
1243 dtype='c16').reshape(len(directions), nk, n, n)
1244 return mat, dmat
1246 def _bloch_derivatives(
1247 self, ks: NDArray[np.float64], directions: list[NDArray[np.float64]], positions: bool,
1248 params: dict | None = None,
1249 ) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
1250 r'''
1251 Private method. :math:`H(\mathbf{k})` and its exact derivatives
1252 :math:`\partial H/\partial k_u` along the Cartesian unit vectors
1253 *directions*, at the k-points *ks* (shape (nk, dim), coordinates of
1254 *get_ham*). With an overlap, the Lowdin-orthogonalized
1255 :math:`H' = S^{-1/2}HS^{-1/2}` (the eigenvectors :math:`S^{1/2}v` of
1256 the topology tools) and its derivative, with
1257 :math:`\partial S^{-1/2}` from the Daleckii-Krein formula in the
1258 eigenbasis :math:`S = \sum_a s_a|a\rangle\langle a|`:
1259 :math:`\langle a|\partial S^{-1/2}|b\rangle = \langle a|\partial S|b\rangle
1260 (s_a^{-1/2}-s_b^{-1/2})/(s_a-s_b)`, :math:`-s_a^{-3/2}/2` when :math:`s_a = s_b`.
1261 The value functions take *params* (default: *self.params*).
1262 '''
1263 k_cart = ks @ self.k_basis.T
1264 hops = self._hop if params is None else self._hops(params)
1265 ham, dham = self._bloch_sum(hops, k_cart, directions, positions)
1266 ham = ham + (np.diag(self.onsite) + self._onsite_offdiag)[None]
1267 if not self._overlap_hop:
1268 return ham, dham
1269 s, ds = self._bloch_sum(self._overlap_hop, k_cart, directions, positions)
1270 s = s + np.eye(self.norb)[None]
1271 sv, vs = np.linalg.eigh(s)
1272 error_handling.positive_overlap(sv.min())
1273 r = sv ** -0.5
1274 diff = sv[:, :, None] - sv[:, None, :]
1275 close = np.abs(diff) < 1e-10
1276 lk = np.where(close, -0.5 * sv[:, :, None] ** -1.5,
1277 (r[:, :, None] - r[:, None, :]) / np.where(close, 1., diff))
1278 vh = vs.conj().transpose(0, 2, 1)
1279 s_inv = (vs * r[:, None, :]) @ vh
1280 ds_inv = vs @ ((vh @ ds @ vs) * lk) @ vh
1281 dham = ds_inv @ ham @ s_inv + s_inv @ dham @ s_inv + s_inv @ ham @ ds_inv
1282 return s_inv @ ham @ s_inv, dham
1284 def _kubo(
1285 self, ks: NDArray[np.float64], e_fermi: NDArray[np.float64], temperature: float,
1286 directions: list[NDArray[np.float64]], pairs: list[tuple[int, int]], positions: bool,
1287 spin_op: NDArray[np.complex128] | None, response: str = 'hall',
1288 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
1289 r'''
1290 Private method. The Kubo curvature
1291 :math:`F(\mathbf{k}) = \sum_{n\neq m}P_{nm}\,\mathrm{Im}[X_{nm}Y_{mn}]/(E_n-E_m)^2`
1292 at each k-point, Fermi energy and pair (X, Y) of velocities (X the
1293 spin current :math:`\{s, v\}/2` if *spin_op*), with the weights
1294 :math:`P_{nm}` of *response* (see *_response_weights*), e.g.
1295 :math:`P_{nm} = -(f_n-f_m)` for the Hall conductivity
1296 (:math:`F = \sum_n f_n\Omega_n`). Shape (nk, len(e_fermi), len(pairs)).
1297 Pairs of degenerate states (equal weights) are skipped, so :math:`F`
1298 stays finite at band crossings below or above the Fermi level. Also
1299 returns, per k-point, the occupation-independent size of the
1300 curvature, :math:`\sum_{n\neq m}|\mathrm{Im}[X_{nm}Y_{mn}]|/(E_n-E_m)^2`,
1301 which locates the hot spots of the adaptive refinement.
1302 '''
1303 out = np.zeros((len(ks), len(e_fermi), len(pairs)))
1304 hot = np.zeros(len(ks))
1305 chunk = max(1, 400000 // self.norb ** 2)
1306 for c0 in range(0, len(ks), chunk):
1307 ham, dham = self._bloch_derivatives(ks[c0:c0 + chunk], directions, positions)
1308 en, vec = np.linalg.eigh(ham)
1309 vh = vec.conj().transpose(0, 2, 1)
1310 vel = [vh @ dh @ vec for dh in dham]
1311 de = en[:, :, None] - en[:, None, :]
1312 degenerate = np.abs(de) < 1e-9
1313 inv2 = np.where(degenerate, 0., 1. / np.where(degenerate, 1., de) ** 2)
1314 if spin_op is not None:
1315 s_e = vh @ spin_op @ vec
1316 for p, (a, b) in enumerate(pairs):
1317 x = vel[a] if spin_op is None else 0.5 * (s_e @ vel[a] + vel[a] @ s_e)
1318 w = np.imag(x * vel[b].transpose(0, 2, 1)) * inv2
1319 hot[c0:c0 + chunk] += np.abs(w).sum(axis=(1, 2))
1320 for e, mu in enumerate(e_fermi):
1321 weights = _response_weights(en, float(mu), temperature, response)
1322 out[c0:c0 + chunk, e, p] = np.einsum('knm,knm->k', weights, w)
1323 return out, hot
1325 def _hall(
1326 self, e_fermi, temperature, nk, plane, k_fixed, positions, refine, refine_fraction,
1327 spin_axis, response='hall',
1328 ):
1329 '''
1330 Private method. Validate, set the geometry and mesh, integrate the
1331 Kubo curvature (with the adaptive refinement), for
1332 *hall_conductivity* (spin_axis None) and *spin_hall_conductivity*,
1333 and, with the weights of *response*, for *orbital_magnetization*,
1334 *anomalous_nernst_conductivity* and *thermal_hall_conductivity*.
1335 '''
1336 self._check_static()
1337 error_handling.dim_min(self.dim, 2)
1338 error_handling.hermitian_model(self.is_hermitian())
1339 error_handling.fermi_energies(e_fermi)
1340 if response in ('nernst', 'thermal'):
1341 error_handling.positive_real(temperature, 'temperature')
1342 else:
1343 error_handling.positive_real_zero(temperature, 'temperature')
1344 error_handling.plane(plane, self.dim)
1345 error_handling.k_fixed_hall(k_fixed, self.dim, spin_axis is None)
1346 error_handling.boolean(positions, 'positions')
1347 error_handling.positive_int(refine, 'refine')
1348 error_handling.refine_fraction(refine_fraction)
1349 full = self.dim == 3 and k_fixed is None
1350 axes = [0, 1, 2] if full else list(plane)
1351 if isinstance(nk, int):
1352 nk = (nk,) * len(axes)
1353 error_handling.nk(nk, len(axes))
1354 spin_op = None
1355 if spin_axis is not None:
1356 error_handling.spinful(self.spin)
1357 error_handling.spin_axis(spin_axis)
1358 spin_op = np.kron(np.eye(self.n_sites), PAULI[spin_axis] / 2.)
1359 # Directions (Cartesian unit vectors) of the velocities, the pairs
1360 # (u, w) whose curvature Omega_uw = Omega.(u x w) is integrated, and
1361 # the prefactor turning the mesh average into the conductivity.
1362 if full:
1363 directions = list(np.eye(3))
1364 pairs = [(1, 2), (2, 0), (0, 1)]
1365 volume = abs(np.linalg.det(np.array(self.lat.prim_vec, dtype='f8')))
1366 prefactor = 2 * PI / volume
1367 elif self.dim == 3:
1368 i, j = plane
1369 a_l = np.array(self.lat.prim_vec[3 - i - j], dtype='f8')
1370 b_i, b_j = (np.array(self.rec_vec[p], dtype='f8') for p in plane)
1371 u = b_i / np.linalg.norm(b_i)
1372 directions = [u, np.cross(a_l / np.linalg.norm(a_l), u)]
1373 pairs = [(0, 1)]
1374 prefactor = np.linalg.norm(np.cross(b_i, b_j)) / (2 * PI)
1375 else:
1376 directions = [self.k_basis[:, 0], self.k_basis[:, 1]]
1377 pairs = [(0, 1)]
1378 det = np.linalg.det(self.rec_vec_k)
1379 prefactor = self._plane_orientation((0, 1)) * np.sign(det) * abs(det) / (2 * PI)
1380 e_arr = np.atleast_1d(np.asarray(e_fermi, dtype='f8')).ravel()
1381 grids = np.meshgrid(*[np.arange(n) / n for n in nk], indexing='ij')
1382 fracs = np.zeros((grids[0].size, self.dim))
1383 for ax, g in zip(axes, grids):
1384 fracs[:, ax] = g.ravel()
1385 for ax in range(self.dim):
1386 if ax not in axes:
1387 fracs[:, ax] = k_fixed
1388 args = (e_arr, temperature, directions, pairs, positions, spin_op, response)
1389 coarse, size = self._kubo(fracs @ self.rec_vec_k, *args)
1390 total = coarse.sum(axis=0)
1391 if refine > 1:
1392 # Wang, Yates, Souza and Vanderbilt (2006): resample the cells
1393 # where the curvature peaks on a refine^d submesh (midpoints of
1394 # the sub-cells), and replace their coarse value by the average.
1395 # The cells are ranked by the size of the curvature whatever the
1396 # occupations, so that cells straddling the Fermi surface are
1397 # refined together with their neighbours.
1398 n_hot = int(np.ceil(refine_fraction * len(fracs)))
1399 hot = np.argsort(size)[::-1][:n_hot]
1400 sub = np.array(np.meshgrid(*[(np.arange(refine) + 0.5) / refine - 0.5] * len(axes),
1401 indexing='ij')).reshape(len(axes), -1).T
1402 offsets = np.zeros((len(sub), self.dim))
1403 for c, (ax, n) in enumerate(zip(axes, nk)):
1404 offsets[:, ax] = sub[:, c] / n
1405 fine = (fracs[hot][:, None, :] + offsets[None, :, :]).reshape(-1, self.dim)
1406 fine_val = self._kubo(fine @ self.rec_vec_k, *args)[0].reshape(n_hot, len(sub), len(e_arr), len(pairs))
1407 total = total + (fine_val.mean(axis=1) - coarse[hot]).sum(axis=0)
1408 sigma = prefactor * total / len(fracs)
1409 if not full:
1410 sigma = sigma[:, 0]
1411 if np.ndim(e_fermi) == 0:
1412 return sigma[0] if full else float(sigma[0])
1413 return sigma.reshape(np.shape(e_fermi) + sigma.shape[1:])
1415 def hall_conductivity(
1416 self, e_fermi: float | ArrayLike, temperature: float = 0., nk: int | tuple[int, ...] = 60,
1417 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1418 refine: int = 1, refine_fraction: float = 0.05,
1419 ) -> float | NDArray[np.float64]:
1420 r'''
1421 Get the intrinsic anomalous Hall conductivity at any Fermi level --
1422 in a gap, in a band, or at a band crossing -- from the Kubo formula
1423 over a uniform Brillouin-zone mesh,
1425 .. math::
1427 \sigma_{xy} = \frac{e^2}{\hbar}\int_{BZ}\frac{d^dk}{(2\pi)^d}\,
1428 F(\mathbf{k})\, ,\qquad
1429 F = -\sum_{n\neq m}(f_n - f_m)\,
1430 \frac{\mathrm{Im}\left[\langle n|\partial_x H|m\rangle\langle m|\partial_y H|n\rangle\right]}
1431 {(E_n - E_m)^2} = \sum_n f_n\,\Omega_n\, ,
1433 with :math:`f_n = f(E_n)` the Fermi-Dirac occupations and
1435 .. math::
1437 \Omega_n = -2\,\mathrm{Im}\sum_{m\neq n}
1438 \frac{\langle n|\partial_x H|m\rangle\langle m|\partial_y H|n\rangle}{(E_n-E_m)^2}
1439 = -2\,\mathrm{Im}\langle\partial_x u_n|\partial_y u_n\rangle
1441 the Berry curvature of band *n* (:math:`\partial_\mu H = \partial H/\partial k_\mu`,
1442 the convention :math:`\Omega = \partial_xA_y - \partial_yA_x`,
1443 :math:`\mathbf{A} = i\langle u|\nabla_{\mathbf{k}}u\rangle` of
1444 *berry_curvature*). The form with :math:`f_n - f_m` stays finite
1445 where bands touch: degenerate pairs have equal occupations and
1446 drop out, instead of cancelling between two diverging
1447 :math:`\Omega_n`. :math:`\partial H/\partial\mathbf{k}` is built
1448 analytically from the stored hoppings, :math:`\sum_{\mathbf{R}}
1449 i\mathbf{d}\,t\,e^{i\mathbf{k}\cdot\mathbf{d}}` with the Cartesian bond
1450 vector :math:`\mathbf{d}`.
1452 **Sign.** In a gap, :math:`\sigma_{xy} = C\,e^2/h` with :math:`C`
1453 the value *chern_number* returns for the occupied bands: this is the
1454 sign of TKNN (Phys. Rev. Lett. 49, 405 (1982)), whose
1455 :math:`\sigma_H = (e^2/h)\,\frac{i}{2\pi}\int d^2k\,(\langle\partial_1u|\partial_2u\rangle
1456 - \langle\partial_2u|\partial_1u\rangle)` has the same integrand,
1457 :math:`\Omega/2\pi`. The adiabatic (semiclassical) response of a
1458 charge :math:`q` of either sign, :math:`\dot{\mathbf{r}} = \partial_{\mathbf{k}}E
1459 - \dot{\mathbf{k}}\times\boldsymbol\Omega` with
1460 :math:`\hbar\dot{\mathbf{k}} = q\mathbf{E}`, gives the opposite sign
1461 for the tensor of Ohm's law :math:`j_x = \sigma^{\mathrm{Ohm}}_{xy}E_y`:
1462 :math:`\sigma^{\mathrm{Ohm}}_{xy} = -(q^2/\hbar)\int F\,d^dk/(2\pi)^d`,
1463 the convention of Wang, Yates, Souza and Vanderbilt (Phys. Rev. B
1464 74, 195118 (2006)) and of the Kubo-Bastin formula (*tbkit.kpm.hall_conductivity*
1465 returns, likewise, minus it). So the value returned is
1466 :math:`\sigma^{\mathrm{Ohm}}_{yx} = -\sigma^{\mathrm{Ohm}}_{xy}`; e.g. a
1467 massive Dirac cone :math:`v(k_x\sigma_x + k_y\sigma_y) + m\sigma_z`
1468 gives :math:`+\frac{e^2}{2h}\frac{m}{|E_F|}` outside its gap (the
1469 literature's :math:`\sigma_{xy} = -\frac{e^2}{2h}\frac{m}{|E_F|}`), and
1470 :math:`\frac{e^2}{2h}\mathrm{sgn}(m)` inside.
1472 **Orbital positions.** The velocity is :math:`\partial H/\partial\mathbf{k}`
1473 in the gauge of *positions*: with ``positions=True`` (the default)
1474 the bond vector is :math:`\mathbf{d} = \mathbf{R} + \boldsymbol\tau_j
1475 - \boldsymbol\tau_i`, the matrix of the physical velocity
1476 :math:`i[H, \mathbf{r}]` of a tight-binding model whose position
1477 operator is diagonal on its orbitals (the one *kpm.hall_conductivity*
1478 uses on a torus). With ``positions=False``, every orbital sits on
1479 its cell origin (the periodic gauge of *get_ham*). The two differ by
1480 :math:`i(E_n-E_m)\langle n|\boldsymbol\tau|m\rangle` in the
1481 velocity's matrix elements, so the band curvatures differ by
1482 :math:`\nabla\times\langle n|\boldsymbol\tau|n\rangle`, a total
1483 derivative: the conductivities agree whenever the bands below
1484 :math:`E_F` are all full (in a gap, at any temperature well below
1485 the gap) and differ in general at partial filling, by the circulation
1486 :math:`\frac{1}{2\pi}\oint_{FS}\langle n|\boldsymbol\tau|n\rangle\cdot d\mathbf{k}`
1487 around the Fermi surface. The default is the physical one.
1489 **Units and dimensions.** In 2D, :math:`\sigma` is in units of
1490 :math:`e^2/h` (:math:`\frac{1}{2\pi}\int d^2k\,F`). In 3D, with
1491 *k_fixed* given, the value of the plane (the conductance of that
1492 slice of the Brillouin zone, in :math:`e^2/h`, oriented like
1493 *chern_number* along :math:`\mathbf{a}_l`); with ``k_fixed=None``,
1494 the full 3D conductivity as the Hall (pseudo)vector
1495 :math:`(\sigma_{yz}, \sigma_{zx}, \sigma_{xy}) = \frac{1}{(2\pi)^2}\int d^3k\,\mathbf{F}`
1496 in units of :math:`e^2/(h\cdot\mathrm{length})` (length the unit of
1497 *prim_vec*): :math:`\frac{1}{2\pi}\sum_l\bar C_l\,\mathbf{b}_l`, e.g.
1498 :math:`C/c` along :math:`z` for a stack of Chern layers :math:`c` apart.
1500 **Convergence.** The integrand peaks where bands nearly touch close
1501 to :math:`E_F` (hot spots), and jumps across the Fermi surface. A
1502 *temperature* smooths the latter. *refine* resamples the fraction
1503 *refine_fraction* of mesh cells with the largest :math:`|F|` on a
1504 finer submesh (*refine* points per direction), in the spirit of the
1505 adaptive mesh refinement of Wang, Yates, Souza and Vanderbilt
1506 (2006); the cells are ranked by the size of the curvature whatever
1507 the occupations. It is meant for metals: in a gap the plain uniform
1508 mesh already converges exponentially, and refining part of it
1509 costs a little accuracy (about :math:`4\cdot10^{-4}` for the
1510 Haldane model on a 60 x 60 mesh). With an overlap (*set_overlap*), :math:`H` is first
1511 Lowdin-orthogonalized, :math:`S^{-1/2}HS^{-1/2}`, as by the topology
1512 tools (exact derivatives of :math:`S^{-1/2}` included).
1514 :param e_fermi: Real number, or array of real numbers. Fermi energies
1515 (all of them are computed from one diagonalization of the mesh).
1516 :param temperature: Positive real or zero. Default value 0. In
1517 energy units (:math:`k_B = 1`).
1518 :param nk: Positive integer, or tuple (two integers for a plane,
1519 three for the full 3D conductivity). Default value 60. k-mesh.
1520 :param plane: Tuple of two integers. Default value (0, 1). 3D only:
1521 the reciprocal vectors spanning the plane (see *berry_curvature*).
1522 :param k_fixed: Real number or None. Default value None. 3D only:
1523 the fractional coordinate of the plane along the remaining
1524 reciprocal vector; None for the full 3D conductivity. Ignored in 2D.
1525 :param positions: Boolean. Default value True. Include the orbital
1526 positions in the bond vectors (see above).
1527 :param refine: Positive integer. Default value 1 (no refinement).
1528 Submesh points per direction in the refined cells.
1529 :param refine_fraction: Real in (0, 1]. Default value 0.05. Fraction
1530 of the mesh cells that are refined.
1532 :returns:
1533 * **sigma** -- Real number (or ndarray shaped like *e_fermi*); in
1534 3D with ``k_fixed=None``, ndarray of shape (3,) (or
1535 ``e_fermi.shape + (3,)``).
1537 Example usage::
1539 e_f = np.linspace(-3., 3., 301)
1540 sigma = hal.hall_conductivity(e_f, nk=120) # plateau at C in the gap
1541 '''
1542 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1543 refine_fraction, None)
1545 def spin_hall_conductivity(
1546 self, e_fermi: float | ArrayLike, temperature: float = 0., nk: int | tuple[int, int] = 60,
1547 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1548 spin_axis: str = 'z', refine: int = 1, refine_fraction: float = 0.05,
1549 ) -> float | NDArray[np.float64]:
1550 r'''
1551 Get the intrinsic spin Hall conductivity of a spinful model
1552 (``spin=True``): the Kubo formula of *hall_conductivity* with the
1553 charge current along :math:`x` replaced by the spin current
1554 (Sinova et al., Phys. Rev. Lett. 92, 126603 (2004))
1556 .. math::
1558 j^{s}_x = \tfrac12\{s, v_x\}\, ,\qquad s = \tfrac{\hbar}{2}\sigma_{a}\, ,
1560 i.e. :math:`\sigma^{s}_{xy} = \frac{e}{2\pi}\cdot\frac{1}{2\pi}\int d^2k\,F^s`,
1561 :math:`F^s = -\sum_{n\neq m}(f_n-f_m)\,\mathrm{Im}[\langle n|j^s_x|m\rangle\langle m|\partial_yH|n\rangle]/(E_n-E_m)^2`
1562 (:math:`j^s_x` with :math:`s = \sigma_a/2`, in units of :math:`\hbar`).
1563 The unit is :math:`e/2\pi`: one :math:`e/\hbar` from the electric
1564 field, one :math:`\hbar/2` from the spin, and :math:`1/2\pi` from the
1565 Brillouin-zone integral, as for :math:`e^2/h = (e^2/\hbar)/2\pi`.
1566 When :math:`s_z` is conserved, the two spins decouple and
1567 :math:`j^s_x = \pm\frac{\hbar}{2}v_x` in each, so in a gap
1568 :math:`\sigma^s_{xy} = \frac{e}{2\pi}\cdot\frac{C_\uparrow - C_\downarrow}{2}`,
1569 :math:`C_\sigma` the Chern number of each spin sector (the sign of
1570 *hall_conductivity* and *chern_number*): :math:`\pm e/2\pi` for the
1571 Kane-Mele quantum spin Hall insulator (Kane and Mele 2005). Rashba
1572 coupling breaks the conservation of :math:`s_z`, and the value is
1573 then no longer quantized (the :math:`\mathbb{Z}_2` invariant, see
1574 *z2_invariant*, still is).
1576 :param e_fermi: See *hall_conductivity*.
1577 :param temperature: See *hall_conductivity*.
1578 :param nk: Positive integer, or tuple of 2 integers. Default value 60.
1579 :param plane: See *hall_conductivity*.
1580 :param k_fixed: Real number or None. Default value None. 3D only: the
1581 plane (required in 3D, where only planes are computed).
1582 :param positions: Boolean. Default value True. See *hall_conductivity*.
1583 :param spin_axis: 'x', 'y' or 'z'. Default value 'z'. Spin component
1584 :math:`a` carried by the current.
1585 :param refine: See *hall_conductivity*.
1586 :param refine_fraction: See *hall_conductivity*.
1588 :returns:
1589 * **sigma_s** -- Real number (or ndarray shaped like *e_fermi*), in
1590 units of :math:`e/2\pi`.
1591 '''
1592 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1593 refine_fraction, spin_axis)
1595 # ------------------------------------------------------------------
1596 # Berry-phase response beyond the Hall conductivity
1597 # ------------------------------------------------------------------
1599 def orbital_magnetization(
1600 self, e_fermi: float | ArrayLike, temperature: float = 0., nk: int | tuple[int, ...] = 60,
1601 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1602 refine: int = 1, refine_fraction: float = 0.05,
1603 ) -> float | NDArray[np.float64]:
1604 r'''
1605 Get the orbital magnetization in the modern theory (Thonhauser,
1606 Ceresoli, Vanderbilt and Resta, Phys. Rev. Lett. 95, 137205 (2005);
1607 Xiao, Shi and Niu, Phys. Rev. Lett. 95, 137204 (2005); Ceresoli et
1608 al., Phys. Rev. B 74, 024408 (2006)), at any Fermi level and
1609 temperature:
1611 .. math::
1613 M = \frac{e}{h}\,\frac{1}{2\pi}\int_{BZ} d^2k\sum_n\left[f_n\,m_n
1614 + g_n\,\Omega_n\right],\qquad
1615 m_n = \mathrm{Im}\langle\partial_xu_n|(H-E_n)|\partial_yu_n\rangle\, ,
1617 with :math:`g = k_BT\ln(1+e^{-(E-\mu)/k_BT})` (:math:`\mu - E` below
1618 :math:`\mu` and 0 above, at :math:`T = 0`). The first term is the
1619 self-rotation of the wavepackets, the second the circulation of
1620 the Berry-curvature currents at the edges. At :math:`T = 0` both
1621 combine into :math:`\sum_nf_n\,\mathrm{Im}\langle\partial_xu_n|(H+E_n-2\mu)|\partial_yu_n\rangle`,
1622 the formula of Thonhauser et al. It is evaluated, like
1623 *hall_conductivity*, from :math:`\langle n|\partial H|m\rangle`
1624 (no derivatives of the eigenstates), so it is finite where bands
1625 touch.
1627 **Sign and units.** The sign is that of *hall_conductivity*, so
1628 that the Streda formula reads :math:`\partial M/\partial\mu = \sigma`:
1629 in a gap (:math:`T = 0`), :math:`M` is linear in :math:`\mu` with
1630 slope the Chern number :math:`C`. In 2D, :math:`M` (a magnetic
1631 moment per area, i.e. a current) is in units of :math:`e/h` times
1632 the energy unit; with ``k_fixed=None`` in 3D, the vector
1633 :math:`(M_x, M_y, M_z)` per length. *positions* matters as for
1634 *hall_conductivity*, and the physical choice is the default.
1636 :param e_fermi: See *hall_conductivity*.
1637 :param temperature: See *hall_conductivity*.
1638 :param nk: See *hall_conductivity*.
1639 :param plane: See *hall_conductivity*.
1640 :param k_fixed: See *hall_conductivity*.
1641 :param positions: See *hall_conductivity*.
1642 :param refine: See *hall_conductivity*.
1643 :param refine_fraction: See *hall_conductivity*.
1645 :returns:
1646 * **M** -- Real number (or ndarray shaped like *e_fermi*); in 3D
1647 with ``k_fixed=None``, ndarray of shape (3,) (or ``e_fermi.shape + (3,)``).
1649 Example usage::
1651 mu = np.linspace(-0.5, 0.5, 11) # in the gap of a Chern insulator
1652 M = hal.orbital_magnetization(mu) # slope C = hall_conductivity(mu)
1653 '''
1654 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1655 refine_fraction, None, 'magnetization')
1657 def anomalous_nernst_conductivity(
1658 self, e_fermi: float | ArrayLike, temperature: float, nk: int | tuple[int, ...] = 60,
1659 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1660 refine: int = 1, refine_fraction: float = 0.05,
1661 ) -> float | NDArray[np.float64]:
1662 r'''
1663 Get the intrinsic anomalous Nernst (transverse thermoelectric)
1664 conductivity :math:`\alpha_{xy}`, :math:`j_x = \alpha_{xy}(-\partial_yT)`,
1665 from the Berry curvature weighted by the entropy of each state
1666 (Xiao, Yao, Fang and Niu, Phys. Rev. Lett. 97, 026603 (2006)):
1668 .. math::
1670 \alpha_{xy} = \frac{ek_B}{h}\,\frac{1}{2\pi}\int_{BZ}d^2k\sum_n s_n\,\Omega_n\, ,
1671 \qquad s = -f\ln f - (1-f)\ln(1-f)\, .
1673 Only the states within a few :math:`k_BT` of :math:`\mu` carry
1674 entropy, so it vanishes in a gap at low temperature. It obeys the
1675 Mott relation exactly,
1676 :math:`\alpha_{xy}(\mu, T) = \frac{1}{eT}\int dE\,(E-\mu)\left(-\frac{\partial f}{\partial E}\right)\sigma_{xy}(E)`,
1677 with :math:`\sigma_{xy}(E)` the :math:`T = 0` *hall_conductivity*
1678 at Fermi energy :math:`E`, which fixes the sign. The pair form of
1679 *hall_conductivity* keeps it finite where bands touch.
1681 **Units.** :math:`ek_B/h` in 2D (:math:`k_B = 1`, temperatures in
1682 energy units), per length for the 3D vector (``k_fixed=None``).
1684 :param e_fermi: See *hall_conductivity*.
1685 :param temperature: Positive real number. In energy units (:math:`k_B = 1`).
1686 :param nk: See *hall_conductivity*.
1687 :param plane: See *hall_conductivity*.
1688 :param k_fixed: See *hall_conductivity*.
1689 :param positions: See *hall_conductivity*.
1690 :param refine: See *hall_conductivity*.
1691 :param refine_fraction: See *hall_conductivity*.
1693 :returns:
1694 * **alpha** -- Real number (or ndarray, as *hall_conductivity*).
1696 Example usage::
1698 alpha = hal.anomalous_nernst_conductivity(np.linspace(-3, 3, 301), 0.05)
1699 '''
1700 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1701 refine_fraction, None, 'nernst')
1703 def thermal_hall_conductivity(
1704 self, e_fermi: float | ArrayLike, temperature: float, nk: int | tuple[int, ...] = 60,
1705 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1706 refine: int = 1, refine_fraction: float = 0.05,
1707 ) -> float | NDArray[np.float64]:
1708 r'''
1709 Get the intrinsic thermal Hall conductivity of the electrons
1710 :math:`\kappa_{xy}`, :math:`j^Q_x = \kappa_{xy}(-\partial_yT)`
1711 (Qin, Niu and Shi, Phys. Rev. Lett. 107, 236601 (2011)):
1713 .. math::
1715 \kappa_{xy} = \frac{k_B^2T}{h}\,\frac{1}{2\pi}\int_{BZ}d^2k
1716 \sum_n c_2(x_n)\,\Omega_n\, ,\qquad
1717 c_2(x) = \int_x^\infty y^2\left(-\frac{\partial f}{\partial y}\right)dy\, ,
1719 :math:`x_n = (E_n-\mu)/k_BT` (:math:`c_2 = \pi^2/3` deep below
1720 :math:`\mu`, 0 far above). Equivalently, exactly,
1721 :math:`\kappa_{xy}(\mu, T) = \frac{1}{e^2T}\int dE\,(E-\mu)^2\left(-\frac{\partial f}{\partial E}\right)\sigma_{xy}(E)`
1722 with the :math:`T = 0` *hall_conductivity*. In a gap at low
1723 temperature it is the Wiedemann-Franz value
1724 :math:`\kappa_{xy} = \frac{\pi^2k_B^2}{3e^2}T\sigma_{xy}`: for a Chern
1725 insulator the quantized :math:`C\,\pi^2k_B^2T/3h`.
1727 **Units.** :math:`k_B^2/h` times the temperature (in energy units,
1728 :math:`k_B = 1`) in 2D, so that :math:`\kappa/T = \pi^2C/3` on a
1729 Chern plateau; per length for the 3D vector (``k_fixed=None``).
1730 The sign is that of *hall_conductivity*.
1732 :param e_fermi: See *hall_conductivity*.
1733 :param temperature: Positive real number. In energy units (:math:`k_B = 1`).
1734 :param nk: See *hall_conductivity*.
1735 :param plane: See *hall_conductivity*.
1736 :param k_fixed: See *hall_conductivity*.
1737 :param positions: See *hall_conductivity*.
1738 :param refine: See *hall_conductivity*.
1739 :param refine_fraction: See *hall_conductivity*.
1741 :returns:
1742 * **kappa** -- Real number (or ndarray, as *hall_conductivity*).
1744 Example usage::
1746 kappa = hal.thermal_hall_conductivity(0., 0.01) # = pi^2/3 C T in the gap
1747 '''
1748 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1749 refine_fraction, None, 'thermal')
1751 def _axion_rate(
1752 self, ks: NDArray[np.float64], bands: list[int], params: dict, param: str,
1753 positions: bool,
1754 ) -> float:
1755 r'''
1756 Private method. The mesh average, over *ks*, of
1757 :math:`\mathrm{tr}[F_{\lambda x}F_{yz} + F_{\lambda y}F_{zx} + F_{\lambda z}F_{xy}]`
1758 for the bands *bands* at the parameters *params*, the derivative
1759 along *param* by central differences. The non-Abelian curvature is
1760 :math:`F_{\mu\nu} = i(X_\mu^\dagger X_\nu - X_\nu^\dagger X_\mu)`,
1761 :math:`(X_\mu)_{lm} = \langle l|\partial_\mu H|m\rangle/(E_m-E_l)`
1762 (:math:`m` in *bands*, :math:`l` not), i.e.
1763 :math:`Q|\partial_\mu u_m\rangle` in the eigenbasis: gauge covariant.
1764 '''
1765 value = params[param]
1766 step = 1e-5 * max(1., abs(value))
1767 others = [b for b in range(self.norb) if b not in bands]
1768 total, gap = 0., np.inf
1769 chunk = max(1, 400000 // self.norb ** 2)
1770 for c0 in range(0, len(ks), chunk):
1771 k = ks[c0:c0 + chunk]
1772 ham, dham = self._bloch_derivatives(k, list(np.eye(3)), positions, params)
1773 plus = self._bloch_derivatives(k, [], positions, dict(params, **{param: value + step}))[0]
1774 minus = self._bloch_derivatives(k, [], positions, dict(params, **{param: value - step}))[0]
1775 en, vec = np.linalg.eigh(ham)
1776 occ, emp = vec[:, :, bands], vec[:, :, others]
1777 de = en[:, others][:, :, None] - en[:, bands][:, None, :]
1778 gap = min(gap, np.abs(de).min())
1779 emp_h = emp.conj().transpose(0, 2, 1)
1780 x = [-(emp_h @ d @ occ) / de for d in [(plus - minus) / (2 * step), *dham]]
1781 f = {}
1782 for a, b in ((0, 1), (0, 2), (0, 3), (1, 2), (2, 3), (3, 1)):
1783 xa_h = x[a].conj().transpose(0, 2, 1)
1784 xb_h = x[b].conj().transpose(0, 2, 1)
1785 f[a, b] = 1j * (xa_h @ x[b] - xb_h @ x[a])
1786 trace = (np.einsum('kab,kba->k', f[0, 1], f[2, 3]) + np.einsum('kab,kba->k', f[0, 2], f[3, 1])
1787 + np.einsum('kab,kba->k', f[0, 3], f[1, 2]))
1788 total += trace.real.sum()
1789 error_handling.path_gap(gap, param, value)
1790 return total / len(ks)
1792 def axion_angle(
1793 self, bands: list[int], param: str, values: ArrayLike, nk: int | tuple[int, int, int] = 16,
1794 theta0: float = 0., positions: bool = True,
1795 ) -> NDArray[np.float64]:
1796 r'''
1797 Get the axion angle :math:`\theta` of a 3D insulator along a path of
1798 a parameter of its value functions (*set_hopping*), from the second
1799 Chern form (Qi, Hughes and Zhang, Phys. Rev. B 78, 195424 (2008);
1800 Essin, Moore and Vanderbilt, Phys. Rev. Lett. 102, 146805 (2009)):
1802 .. math::
1804 \frac{d\theta}{d\lambda} = \frac{1}{2\pi}\int_{BZ}d^3k\;
1805 \mathrm{tr}\left[F_{\lambda x}F_{yz} + F_{\lambda y}F_{zx} + F_{\lambda z}F_{xy}\right],
1807 :math:`F_{\mu\nu}` the non-Abelian Berry curvature of the bands
1808 *bands* in :math:`(\lambda, k_x, k_y, k_z)` space, so that
1809 :math:`\theta(\lambda) = \theta_0 + \int_{\lambda_0}^{\lambda}\frac{d\theta}{d\lambda'}d\lambda'`.
1810 This is the variation of the Chern-Simons axion coupling, computed
1811 from gauge-invariant quantities only (no smooth gauge is needed),
1812 and :math:`\theta_0` is its value at ``values[0]``: e.g. 0
1813 for a path starting from an atomic insulator. A closed path (a
1814 pumping cycle) changes :math:`\theta` by :math:`2\pi C_2`, with
1815 :math:`C_2` the second Chern number of the 4D family. At a point
1816 with time-reversal or inversion symmetry, :math:`\theta` is 0 or
1817 :math:`\pi` (mod :math:`2\pi`): :math:`\pi` for a strong topological
1818 insulator, cf. *z2_indices_3d*. The sign of :math:`\theta` follows
1819 the orientation :math:`(\lambda, k_x, k_y, k_z)` above.
1821 The path must stay gapped: the bands *bands* must be separated
1822 from the others at every k-point and every value (a ValueError
1823 names the value where they touch). The integrand is smooth, so the
1824 uniform k-mesh converges exponentially; the :math:`\lambda`
1825 integral is a trapezoid over *values* (error :math:`O(\Delta\lambda^2)`).
1826 :math:`\partial H/\partial\lambda` is a central difference, and
1827 :math:`\partial H/\partial\mathbf{k}` is exact (see
1828 *hall_conductivity*, whose *positions* convention applies).
1830 :param bands: List of integers. The occupied bands (e.g. the lower
1831 half of the spectrum).
1832 :param param: String. Name of the parameter of the value functions
1833 that the path varies (see *set_params*; the other parameters keep
1834 their values in *params*).
1835 :param values: 1D array of at least two real numbers, increasing
1836 or decreasing. The parameter values along the path.
1837 :param nk: Positive integer, or tuple of 3 positive integers.
1838 Default value 16. k-mesh.
1839 :param theta0: Real number. Default value 0. :math:`\theta` at ``values[0]``.
1840 :param positions: Boolean. Default value True. See *hall_conductivity*.
1842 :returns:
1843 * **theta** -- ndarray, shape ``(len(values),)``: :math:`\theta`
1844 at each value, not reduced modulo :math:`2\pi`.
1846 Example usage::
1848 phi = np.linspace(0., 2 * np.pi, 41)
1849 theta = ks.axion_angle([0, 1], 'phi', phi) # theta[-1] = 2 pi C2
1850 '''
1851 self._check_static()
1852 error_handling.dim_exact(self.dim, 3)
1853 error_handling.hermitian_model(self.is_hermitian())
1854 error_handling.band_indices(bands, self.norb)
1855 error_handling.path_parameter(param)
1856 error_handling.path_values(values)
1857 if isinstance(nk, int):
1858 nk = (nk,) * 3
1859 error_handling.nk(nk, 3)
1860 error_handling.real_number(theta0, 'theta0')
1861 error_handling.boolean(positions, 'positions')
1862 values = np.asarray(values, dtype='f8')
1863 fracs = np.array(np.meshgrid(*[np.arange(n) / n for n in nk], indexing='ij')).reshape(3, -1).T
1864 ks = fracs @ self.rec_vec_k
1865 volume = abs(np.linalg.det(self.rec_vec_k))
1866 rates = [self._axion_rate(ks, bands, dict(self.params, **{param: float(v)}), param, positions)
1867 for v in values]
1868 return theta0 + cumulative_trapezoid(np.array(rates) * volume / (2 * PI), values, initial=0.)
1870 # ------------------------------------------------------------------
1871 # Finite samples and non-Hermitian band theory
1872 # ------------------------------------------------------------------
1874 def _hop_cells(self) -> list[tuple[int, int, tuple[int, ...], complex]]:
1875 r'''
1876 Private method. The stored hoppings with their lattice vector as
1877 integers :math:`(n_1, ...)` (from :math:`\mathbf{R}\cdot\mathbf{b}_i/2\pi`).
1878 '''
1879 rec = np.array(self.rec_vec, dtype='f8')
1880 return [(i, j, tuple(int(n) for n in np.rint(rec @ R / (2 * PI))), t)
1881 for i, j, R, t in self._hop]
1883 def _finite_entries(
1884 self, n_cells: int | tuple[int, ...], periodic: bool,
1885 ) -> tuple[tuple[int, ...], int, list[tuple[NDArray, NDArray, int, int, tuple, NDArray, complex]]]:
1886 '''
1887 Private method. For each stored hopping, the row and column indices
1888 it fills in a finite sample of *n_cells* cells (see *finite_ham*),
1889 with its orbitals, lattice vector (integers and Cartesian) and amplitude.
1890 '''
1891 if isinstance(n_cells, int):
1892 n_cells = (n_cells,) * self.dim
1893 n_tot = int(np.prod(n_cells))
1894 cells = np.array(np.meshgrid(*[np.arange(n) for n in n_cells], indexing='ij'))
1895 cells = cells.reshape(self.dim, -1).T # (n_tot, dim)
1896 strides = np.cumprod((1,) + tuple(n_cells[:-1]))
1897 index = cells @ strides
1898 order = np.argsort(index)
1899 cells = cells[order]
1900 dims = np.array(n_cells)
1901 entries = []
1902 for (i, j, R, t), (_, _, R_cart, _) in zip(self._hop_cells(), self._hop):
1903 target = cells + np.array(R)[None, :]
1904 if periodic:
1905 target = target % dims
1906 keep = np.ones(n_tot, bool)
1907 else:
1908 keep = np.all((target >= 0) & (target < dims), axis=1)
1909 rows = np.arange(n_tot)[keep] * self.norb + i
1910 cols = (target[keep] @ strides) * self.norb + j
1911 entries.append((rows, cols, i, j, R, R_cart, t))
1912 return n_cells, n_tot, entries
1914 def finite_ham(
1915 self, n_cells: int | tuple[int, ...], periodic: bool = False, sparse: bool = False,
1916 ) -> NDArray[np.complex128] | sp.csr_matrix:
1917 r'''
1918 Get the real-space Hamiltonian of a finite sample of
1919 :math:`N_1\times N_2\times\dots` unit cells of the model -- open
1920 boundaries (a flake, a finite chain), or periodic ones (a torus).
1921 Orbital *o* of cell :math:`(n_1, n_2, n_3)` is row
1922 :math:`((n_3 N_2 + n_2) N_1 + n_1)\,\mathrm{norb} + o`.
1924 :param n_cells: Positive integer, or tuple of *dim* positive integers.
1925 :param periodic: Boolean. Default value False. Periodic boundary
1926 conditions (the spectrum is then that of *get_bands* on the
1927 :math:`N_1\times N_2\times\dots` k-mesh), else open ones.
1928 :param sparse: Boolean. Default value False. Return a sparse CSR
1929 matrix (for large samples, e.g. with *tbkit.kpm*) instead of a
1930 dense array.
1932 :returns:
1933 * **ham** -- Complex ndarray (or CSR matrix), shape (N*norb, N*norb), N the number of cells.
1934 '''
1935 error_handling.nk(n_cells, self.dim)
1936 error_handling.boolean(periodic, 'periodic')
1937 error_handling.boolean(sparse, 'sparse')
1938 _, n_tot, entries = self._finite_entries(n_cells, periodic)
1939 onsite = np.diag(self.onsite) + self._onsite_offdiag
1940 if sparse:
1941 onsite_rows, onsite_cols = np.nonzero(onsite)
1942 rows = [c * self.norb + onsite_rows for c in range(n_tot)]
1943 cols = [c * self.norb + onsite_cols for c in range(n_tot)]
1944 vals = [onsite[onsite_rows, onsite_cols]] * n_tot
1945 for r, c, _, _, _, _, t in entries:
1946 rows.append(r)
1947 cols.append(c)
1948 vals.append(np.full(len(r), t, dtype='c16'))
1949 n = n_tot * self.norb
1950 return sp.coo_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
1951 shape=(n, n)).tocsr()
1952 ham = np.zeros((n_tot * self.norb, n_tot * self.norb), 'c16')
1953 for cell in range(n_tot):
1954 ham[cell*self.norb:(cell+1)*self.norb, cell*self.norb:(cell+1)*self.norb] += onsite
1955 for r, c, _, _, _, _, t in entries:
1956 np.add.at(ham, (r, c), t)
1957 return ham
1959 def finite_velocity(
1960 self, n_cells: int | tuple[int, ...], periodic: bool = False, sparse: bool = False,
1961 ) -> list[NDArray[np.complex128] | sp.csr_matrix]:
1962 r'''
1963 Get the velocity operators :math:`\mathbf{v} = i[H, \mathbf{r}]`
1964 (:math:`\hbar = 1`) of the finite sample of *finite_ham*, from the
1965 Cartesian bond vectors of the hoppings rather than from the site
1966 coordinates:
1968 .. math::
1970 (v_\alpha)_{ab} = i\,H_{ab}\,d_\alpha\, ,\qquad
1971 \mathbf{d} = \mathbf{R} + \boldsymbol\tau_j - \boldsymbol\tau_i\, ,
1973 so that they are well defined on a torus (``periodic=True``), where
1974 the position operator is not: a bond that wraps around the sample
1975 keeps its short displacement. They are the real-space counterpart of
1976 :math:`\partial H/\partial\mathbf{k}` with ``positions=True``
1977 (*hall_conductivity*), and the input of *tbkit.kpm.hall_conductivity*.
1979 :param n_cells: See *finite_ham*.
1980 :param periodic: Boolean. Default value False. See *finite_ham*.
1981 :param sparse: Boolean. Default value False. See *finite_ham*.
1983 :returns:
1984 * **vel** -- List of *space_dim* complex ndarrays (or CSR
1985 matrices), shape (N*norb, N*norb): :math:`v_x, v_y, (v_z)`.
1986 '''
1987 error_handling.nk(n_cells, self.dim)
1988 error_handling.boolean(periodic, 'periodic')
1989 error_handling.boolean(sparse, 'sparse')
1990 _, n_tot, entries = self._finite_entries(n_cells, periodic)
1991 tau = self.orbital_positions()
1992 n = n_tot * self.norb
1993 vel = []
1994 for alpha in range(self.space_dim):
1995 rows, cols, vals = [np.zeros(0, int)], [np.zeros(0, int)], [np.zeros(0, 'c16')]
1996 for r, c, i, j, _, R_cart, t in entries:
1997 d = R_cart[alpha] + tau[j, alpha] - tau[i, alpha]
1998 rows.append(r)
1999 cols.append(c)
2000 vals.append(np.full(len(r), 1j * t * d, dtype='c16'))
2001 mat = sp.coo_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
2002 shape=(n, n)).tocsr()
2003 vel.append(mat if sparse else mat.toarray())
2004 return vel
2006 def get_ham_beta(self, beta: complex) -> NDArray[np.complex128]:
2007 r'''
2008 Get the non-Bloch Hamiltonian of a 1D model,
2009 :math:`H(\beta) = \sum_{n} t(n)\, \beta^{n}`: the Bloch Hamiltonian
2010 continued off the unit circle, :math:`\beta = e^{ik}` with complex
2011 :math:`k`. On the generalized Brillouin zone (see *gbz*), its
2012 spectrum is the open-chain spectrum of a non-Hermitian chain.
2014 :param beta: Complex number, nonzero.
2016 :returns:
2017 * **ham** -- Complex ndarray, shape (norb, norb).
2018 '''
2019 error_handling.dim_exact(self.dim, 1)
2020 error_handling.nonzero_number(beta, 'beta')
2021 ham = np.diag(self.onsite).astype('c16') + self._onsite_offdiag
2022 for i, j, (n,), t in self._hop_cells():
2023 ham[i, j] += t * complex(beta) ** n
2024 return ham
2026 def spectral_winding(self, e_ref: complex = 0., nk: int = 400) -> float:
2027 r'''
2028 Get the spectral winding number of a 1D (non-Hermitian) model about
2029 the reference energy :math:`E_r`,
2031 .. math::
2033 W(E_r) = \frac{1}{2\pi i}\oint_0^{2\pi} dk\,
2034 \partial_k \ln\det\left[H(k) - E_r\right]\, ,
2036 the number of times the periodic-chain spectrum winds around
2037 :math:`E_r` in the complex plane. :math:`W \neq 0` for some
2038 :math:`E_r` is the criterion for the non-Hermitian skin effect: the
2039 open-chain eigenstates then pile up at one end (Okuma, Kawabata,
2040 Shiozaki and Sato, Phys. Rev. Lett. 124, 086801 (2020)).
2042 :param e_ref: Complex number. Default value 0. Reference energy.
2043 :param nk: Positive integer. Default value 400. k-points.
2045 :returns:
2046 * **W** -- Real number, close to an integer (it is ill defined if
2047 :math:`E_r` lies on the periodic spectrum).
2048 '''
2049 error_handling.dim_exact(self.dim, 1)
2050 error_handling.number(e_ref, 'e_ref')
2051 error_handling.positive_int(nk, 'nk')
2052 ks = np.arange(nk + 1) / nk * self.rec_vec_k[0, 0]
2053 dets = np.array([np.linalg.det(self.get_ham([k]) - e_ref * np.eye(self.norb)) for k in ks])
2054 steps = np.angle(dets[1:] / dets[:-1])
2055 return float(np.sum(steps) / (2 * PI))
2057 def gbz(self, energies: ArrayLike) -> NDArray[np.complex128]:
2058 r'''
2059 Get the generalized Brillouin zone of a 1D non-Hermitian model
2060 (Yao and Wang, Phys. Rev. Lett. 121, 086803 (2018); Yokomizo and
2061 Murakami, Phys. Rev. Lett. 123, 066404 (2019)): for each energy
2062 :math:`E`, the roots of :math:`\det[H(\beta)-E] = 0` sorted by modulus,
2063 :math:`|\beta_1|\le\dots\le|\beta_{p+s}|` (:math:`p` the order of its
2064 pole at :math:`\beta = 0`), and the middle pair
2065 :math:`\beta_p, \beta_{p+1}`. The open-chain spectrum, in the limit of
2066 a long chain, is where :math:`|\beta_p| = |\beta_{p+1}|`; those
2067 :math:`\beta` trace the generalized Brillouin zone, on which the
2068 non-Bloch bands :math:`H(\beta)` (see *get_ham_beta*) replace the
2069 Bloch ones. For a Hermitian model it is the unit circle.
2071 :param energies: Complex array of energies (e.g. open-chain eigenvalues
2072 from *finite_ham*).
2074 :returns:
2075 * **beta** -- Complex ndarray, shape (len(energies), 2): :math:`\beta_p, \beta_{p+1}`.
2076 '''
2077 error_handling.dim_exact(self.dim, 1)
2078 energies = np.atleast_1d(np.asarray(energies, dtype='c16'))
2079 error_handling.ndarray_empty(energies, 'energies')
2080 powers = [n for _, _, (n,), _ in self._hop_cells()] + [0]
2081 low, high = self.norb * min(powers), self.norb * max(powers)
2082 n_samples = 2 * (high - low) + 8
2083 theta = 2 * PI * np.arange(n_samples) / n_samples
2084 hams = np.array([self.get_ham_beta(np.exp(1j * th)) for th in theta])
2085 out = np.zeros((len(energies), 2), 'c16')
2086 exps = np.arange(low, high + 1)
2087 for n, e in enumerate(energies):
2088 dets = np.linalg.det(hams - e * np.eye(self.norb)[None])
2089 # Laurent coefficients c_m of det[H(beta) - E] = sum_m c_m beta^m
2090 coef = np.array([np.mean(dets * np.exp(-1j * m * theta)) for m in exps])
2091 nonzero = np.abs(coef) > 1e-12 * np.max(np.abs(coef))
2092 first, last = np.argmax(nonzero), len(coef) - 1 - np.argmax(nonzero[::-1])
2093 p = -exps[first]
2094 roots = np.roots(coef[first:last + 1][::-1])
2095 roots = roots[np.argsort(np.abs(roots))]
2096 error_handling.gbz_roots(len(roots), p)
2097 out[n] = roots[p - 1], roots[p]
2098 return out
2100 # ------------------------------------------------------------------
2101 # Berry phases, Wannier centres, Z2, symmetries, quantum geometry
2102 # ------------------------------------------------------------------
2104 def orbital_positions(self) -> NDArray[np.float64]:
2105 r'''
2106 Get the positions :math:`\boldsymbol\tau` of the orbitals within the
2107 unit cell (each site twice, for its two spin components, if
2108 ``spin=True``).
2110 :returns:
2111 * **tau** -- Real ndarray, shape (norb, space_dim).
2112 '''
2113 tau = np.array([dic['r0'] for dic in self.lat.unit_cell], dtype='f8')
2114 return np.repeat(tau, 2, axis=0) if self.spin else tau
2116 def _loop(
2117 self, bands: list[int], nk: int, direction: int, fracs_perp: dict[int, float],
2118 positions: bool,
2119 ) -> NDArray[np.complex128]:
2120 '''
2121 Private method. Wilson loop matrix of *bands* along the closed path
2122 k0 -> k0 + b_direction (nk steps), at fixed fractional coordinates
2123 *fracs_perp* along the other reciprocal vectors.
2124 '''
2125 frac = np.zeros(self.dim)
2126 for d, f in fracs_perp.items():
2127 frac[d] = f
2128 b = self.rec_vec_k[direction]
2129 k0 = frac @ self.rec_vec_k
2130 # each link carries exp(-i dk.tau): the Bloch sums then include the
2131 # orbital positions, so that the phases locate the Wannier centres
2132 # in space rather than on the cell origin.
2133 dk_space = self.k_basis @ (b / nk)
2134 phase = np.exp(-1j * self.orbital_positions() @ dk_space) if positions \
2135 else np.ones(self.norb)
2136 vs = self._subspaces(k0 + np.arange(nk)[:, None] * b[None] / nk, bands)
2137 links = vs.conj().transpose(0, 2, 1) @ (phase[:, None] * np.roll(vs, -1, axis=0))
2138 # keep only the unitary part of each link (its singular values
2139 # tend to 1 as nk grows): the product then stays unitary
2140 u, _, vh = np.linalg.svd(links)
2141 wilson = np.eye(len(bands), dtype='c16')
2142 for unitary in u @ vh:
2143 wilson = wilson @ unitary
2144 return wilson
2146 def _perp(self, direction: int, k_perp: float | tuple | None) -> dict[int, float]:
2147 '''
2148 Private method. Fractional coordinates of a loop along the other
2149 reciprocal vectors, from *k_perp* (a number in 2D, a pair in 3D).
2150 '''
2151 others = [d for d in range(self.dim) if d != direction]
2152 if k_perp is None:
2153 k_perp = (0.,) * len(others)
2154 if not isinstance(k_perp, tuple):
2155 k_perp = (k_perp,)
2156 error_handling.k_perp(k_perp, len(others))
2157 return dict(zip(others, k_perp))
2159 def berry_phase(
2160 self, bands: int | list[int], nk: int = 100, direction: int = 0,
2161 k_perp: float | tuple[float, float] | None = None, positions: bool = True,
2162 ) -> float:
2163 r'''
2164 Get the Berry phase of a group of bands along a closed loop across the
2165 Brillouin zone, :math:`\mathbf{k}\to\mathbf{k}+\mathbf{b}_{direction}`,
2167 .. math::
2169 \gamma = \oint \mathbf{A}\cdot d\mathbf{k} = -\arg\det W\, ,\qquad
2170 W = \prod_{m} \langle u(\mathbf{k}_m)|u(\mathbf{k}_{m+1})\rangle\, ,
2172 with :math:`\mathbf{A} = i\langle u|\nabla_{\mathbf{k}}u\rangle` (the
2173 convention of *berry_curvature*). In 1D it is the Zak phase: its
2174 value :math:`\gamma/2\pi`, modulo 1, is the centre of the band's
2175 Wannier functions in units of :math:`\mathbf{a}_{direction}` (the
2176 polarization).
2178 :param bands: Band index, or list of band indices.
2179 :param nk: Positive integer. Default value 100. Number of k-points on the loop.
2180 :param direction: Integer. Default value 0. Reciprocal vector the loop follows.
2181 :param k_perp: Real number (2D) or pair of real numbers (3D). Default
2182 value None (zero). Fractional coordinates of the loop along the
2183 other reciprocal vectors.
2184 :param positions: Boolean. Default value True. Include the orbital
2185 positions in the Bloch functions (the physical choice: the phase
2186 then locates the Wannier centres in the crystal). If False, all
2187 orbitals are placed on the cell origin (the periodic gauge of
2188 *get_ham*); the two differ by :math:`\sum_n\mathbf{b}\cdot\boldsymbol\tau`
2189 weighted by the orbital content, and agree on quantized
2190 differences between phases.
2192 :returns:
2193 * **gamma** -- Real number in :math:`(-\pi, \pi]`.
2194 '''
2195 bands = self._check_loop(bands, nk, direction, positions)
2196 wilson = self._loop(bands, nk, direction, self._perp(direction, k_perp), positions)
2197 return float(-np.angle(np.linalg.det(wilson)))
2199 def wannier_centers(
2200 self, bands: int | list[int], nk: int = 100, direction: int = 0,
2201 k_perp: float | tuple[float, float] | None = None, positions: bool = True,
2202 ) -> NDArray[np.float64]:
2203 r'''
2204 Get the (hybrid) Wannier centres of a group of bands: minus the
2205 phases of the eigenvalues of the Wilson loop (see *berry_phase*),
2206 divided by :math:`2\pi`. They are the positions, in units of
2207 :math:`\mathbf{a}_{direction}` and modulo 1, of the maximally
2208 localized Wannier functions along that direction (hybrid: still
2209 Bloch waves along the others, at *k_perp*). Their sum is
2210 :math:`\gamma/2\pi`, modulo 1.
2212 :param bands: Band index, or list of band indices.
2213 :param nk: Positive integer. Default value 100. Number of k-points on the loop.
2214 :param direction: Integer. Default value 0.
2215 :param k_perp: See *berry_phase*.
2216 :param positions: Boolean. Default value True. See *berry_phase*.
2218 :returns:
2219 * **centers** -- Real ndarray, sorted, in :math:`[0, 1)`.
2220 '''
2221 bands = self._check_loop(bands, nk, direction, positions)
2222 wilson = self._loop(bands, nk, direction, self._perp(direction, k_perp), positions)
2223 return np.sort((-np.angle(np.linalg.eigvals(wilson)) / (2*PI)) % 1.)
2225 def wannier_flow(
2226 self, bands: int | list[int], nk: int = 100, nk_perp: int = 51, direction: int = 0,
2227 flow: int | None = None, k_fixed: float = 0., positions: bool = True,
2228 k_range: tuple[float, float] = (0., 1.),
2229 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
2230 r'''
2231 Get the flow of the hybrid Wannier centres (see *wannier_centers*)
2232 along *direction*, as the loop is moved across the Brillouin zone
2233 along the reciprocal vector *flow*. In a Chern insulator the centres
2234 wind by :math:`C` cells over a full period; in a
2235 :math:`\mathbb{Z}_2` topological insulator, the Kramers pairs swap
2236 partners over half a period (see *z2_invariant*).
2238 :param bands: Band index, or list of band indices.
2239 :param nk: Positive integer. Default value 100. k-points on each loop.
2240 :param nk_perp: Positive integer. Default value 51. Number of loops.
2241 :param direction: Integer. Default value 0. Direction of the loops.
2242 :param flow: Integer. Default value None (the first other direction).
2243 Reciprocal vector along which the loops are moved.
2244 :param k_fixed: Real number. Default value 0. 3D only: fractional
2245 coordinate along the remaining reciprocal vector.
2246 :param positions: Boolean. Default value True. See *berry_phase*.
2247 :param k_range: Pair of real numbers. Default value (0, 1). Range of
2248 the fractional coordinate along *flow* (both ends included).
2250 :returns:
2251 * **fracs** -- Real ndarray, shape (nk_perp,). Fractional coordinates along *flow*.
2252 * **centers** -- Real ndarray, shape (nk_perp, len(bands)). Sorted centres, in [0, 1).
2253 '''
2254 bands = self._check_loop(bands, nk, direction, positions)
2255 error_handling.dim_min(self.dim, 2)
2256 error_handling.positive_int(nk_perp, 'nk_perp')
2257 others = [d for d in range(self.dim) if d != direction]
2258 if flow is None:
2259 flow = others[0]
2260 error_handling.flow(flow, direction, self.dim)
2261 error_handling.real_number(k_fixed, 'k_fixed')
2262 error_handling.lims(k_range)
2263 fracs = np.linspace(k_range[0], k_range[1], nk_perp)
2264 centers = np.zeros((nk_perp, len(bands)))
2265 for n, f in enumerate(fracs):
2266 perp = {d: (f if d == flow else k_fixed) for d in others}
2267 wilson = self._loop(bands, nk, direction, perp, positions)
2268 centers[n] = np.sort((-np.angle(np.linalg.eigvals(wilson)) / (2*PI)) % 1.)
2269 return fracs, centers
2271 def z2_invariant(
2272 self, bands: list[int], nk: int = 100, nk_perp: int = 51, direction: int = 0,
2273 flow: int | None = None, k_fixed: float = 0.,
2274 ) -> int:
2275 r'''
2276 Get the :math:`\mathbb{Z}_2` invariant of a time-reversal-symmetric
2277 group of bands (an even number, e.g. all occupied bands of the
2278 Kane-Mele model) from the flow of its hybrid Wannier centres over
2279 half the Brillouin zone (Soluyanov and Vanderbilt, Phys. Rev. B 83,
2280 235401 (2011); Yu, Qi, Bernevig, Fang and Dai, Phys. Rev. B 84,
2281 075119 (2011)): follow the midpoint of the largest gap between the
2282 centres, and count, modulo 2, how many centres it jumps over. It is
2283 odd when the Kramers partners exchange along the flow -- a
2284 topological insulator -- and even otherwise.
2286 Unlike the parity criterion (*parity_z2*), it needs no inversion
2287 symmetry, so it also works with Rashba coupling.
2289 :param bands: List of band indices (an even number of bands).
2290 :param nk: Positive integer. Default value 100. k-points on each loop.
2291 :param nk_perp: Positive integer. Default value 51. Number of loops
2292 over the half zone (the flow must be resolved finely enough that
2293 no two centres cross between neighbouring loops unnoticed).
2294 :param direction: Integer. Default value 0. See *wannier_flow*.
2295 :param flow: Integer. Default value None. See *wannier_flow*.
2296 :param k_fixed: Real number. Default value 0. 3D only: a time-reversal
2297 invariant plane has *k_fixed* 0 or 0.5 (the weak and strong
2298 indices follow from the six such planes).
2300 :returns:
2301 * **nu** -- 0 (trivial) or 1 (topological).
2302 '''
2303 error_handling.even_bands(bands)
2304 _, centers = self.wannier_flow(bands, nk, nk_perp, direction, flow, k_fixed,
2305 positions=False, k_range=(0., 0.5))
2306 gap = [self._largest_gap(c) for c in centers]
2307 jumps = sum(self._between(gap[n], gap[n + 1], centers[n + 1])
2308 for n in range(len(gap) - 1))
2309 return int(jumps % 2)
2311 def z2_indices_3d(
2312 self, bands: list[int], nk: int = 100, nk_perp: int = 51,
2313 ) -> tuple[int, int, int, int]:
2314 r'''
2315 Get the four :math:`\mathbb{Z}_2` indices
2316 :math:`(\nu_0;\nu_1\nu_2\nu_3)` of a 3D time-reversal-invariant
2317 insulator (Fu, Kane and Mele, Phys. Rev. Lett. 98, 106803 (2007);
2318 Moore and Balents, Phys. Rev. B 75, 121306(R) (2007)) from the
2319 :math:`\mathbb{Z}_2` invariants of its six time-reversal-invariant
2320 planes :math:`k_i = 0` and :math:`k_i = 1/2` (fractional coordinates
2321 along :math:`\mathbf{b}_i`), each computed by *z2_invariant*:
2323 .. math::
2325 \nu_0 = \nu(k_i = 0) + \nu(k_i = 1/2) \bmod 2\, ,\qquad
2326 \nu_i = \nu(k_i = 1/2)\, .
2328 The strong index :math:`\nu_0` comes out the same for the three
2329 :math:`i`, which is checked (a ValueError asks for a finer flow
2330 otherwise). :math:`\nu_0 = 1` is a strong topological insulator; a
2331 weak one has :math:`\nu_0 = 0` and some :math:`\nu_i = 1`, i.e.
2332 stacked quantum spin Hall layers normal to
2333 :math:`\nu_1\mathbf{b}_1 + \nu_2\mathbf{b}_2 + \nu_3\mathbf{b}_3`.
2334 No inversion symmetry is needed (compare *parity_z2*).
2336 :param bands: List of band indices (an even number of bands).
2337 :param nk: Positive integer. Default value 100. See *z2_invariant*.
2338 :param nk_perp: Positive integer. Default value 51. See *z2_invariant*.
2340 :returns:
2341 * **indices** -- Tuple of four integers, 0 or 1: :math:`(\nu_0, \nu_1, \nu_2, \nu_3)`.
2342 '''
2343 error_handling.dim_exact(self.dim, 3)
2344 nu = np.zeros((3, 2), int)
2345 for i in range(3):
2346 direction, flow = [d for d in range(3) if d != i]
2347 for n, k_fixed in enumerate((0., 0.5)):
2348 nu[i, n] = self.z2_invariant(bands, nk, nk_perp, direction, flow, k_fixed)
2349 strong = nu.sum(axis=1) % 2
2350 error_handling.strong_index(strong)
2351 return int(strong[0]), int(nu[0, 1]), int(nu[1, 1]), int(nu[2, 1])
2353 def entanglement_spectrum(
2354 self, ks: ArrayLike, region: list[int], bands: int | list[int],
2355 ) -> NDArray[np.float64]:
2356 r'''
2357 Get the k-resolved entanglement spectrum of a group of bands
2358 (Peschel, J. Phys. A 36, L205 (2003); Li and Haldane, Phys. Rev. Lett.
2359 101, 010504 (2008); Fidkowski, Phys. Rev. Lett. 104, 130502 (2010)):
2360 the eigenvalues :math:`\xi_n(\mathbf{k})` of the correlation matrix
2361 :math:`P(\mathbf{k})` of the filled *bands*, restricted to the
2362 orbitals of *region*,
2364 .. math::
2366 C_A(\mathbf{k}) = \left[P(\mathbf{k})\right]_{A}\, ,\qquad
2367 P(\mathbf{k}) = \sum_{n\in\mathrm{bands}}|u_n(\mathbf{k})\rangle\langle u_n(\mathbf{k})|\, ,
2369 between 0 and 1. For a ribbon or a supercell
2370 (*tbkit.moire.supercell*) cut in two halves across its width, the
2371 states of the whole are either inside the region
2372 (:math:`\xi \approx 1`) or outside (:math:`\xi \approx 0`), except
2373 for those localized at the cuts: a topological phase leaves
2374 entanglement modes inside :math:`(0, 1)` -- a mode pinned at 1/2 in
2375 a chiral-symmetric chain, a branch that flows from 0 to 1 across
2376 the zone in a Chern insulator. The single-particle entanglement
2377 energies are :math:`\epsilon = \ln[(1-\xi)/\xi]`.
2379 :param ks: ndarray, shape (nk, dim). k-points (see *get_bands*).
2380 :param region: List of orbital indices (the region :math:`A`).
2381 :param bands: Band index, or list of band indices (the filled bands).
2383 :returns:
2384 * **xi** -- Real ndarray, shape (nk, len(region)), sorted ascending at every k.
2385 '''
2386 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
2387 error_handling.ks(ks, self.dim)
2388 error_handling.region(region, self.norb)
2389 if isinstance(bands, int):
2390 bands = [bands]
2391 error_handling.band_indices(bands, self.norb)
2392 error_handling.hermitian_kspace(self.is_hermitian(), self._overlap_hop)
2393 v = self._subspaces(ks, bands)[:, region, :]
2394 return np.linalg.eigvalsh(v @ v.conj().transpose(0, 2, 1))
2396 @staticmethod
2397 def _largest_gap(centers: NDArray[np.float64]) -> float:
2398 '''
2399 Private method. Midpoint of the largest gap between sorted centres on the circle [0, 1).
2400 '''
2401 ext = np.append(centers, centers[0] + 1.)
2402 n = np.argmax(np.diff(ext))
2403 return float((ext[n] + ext[n + 1]) / 2 % 1.)
2405 @staticmethod
2406 def _between(a: float, b: float, centers: NDArray[np.float64]) -> int:
2407 '''
2408 Private method. Number of centres on the shorter arc between a and b.
2409 '''
2410 d = (b - a) % 1.
2411 start, length = (a, d) if d <= 0.5 else (b, 1. - d)
2412 return int(np.sum((centers - start) % 1. < length))
2414 def _check_loop(self, bands, nk, direction, positions) -> list[int]:
2415 '''
2416 Private method. Validate the arguments shared by the Wilson-loop methods.
2417 '''
2418 if isinstance(bands, int):
2419 bands = [bands]
2420 error_handling.band_indices(bands, self.norb)
2421 error_handling.positive_int(nk, 'nk')
2422 error_handling.direction(direction, self.dim)
2423 error_handling.boolean(positions, 'positions')
2424 return bands
2426 def _trims(self) -> list[NDArray[np.float64]]:
2427 '''
2428 Private method. Time-reversal-invariant momenta: fractional
2429 coordinates in {0, 1/2}^dim, as k points.
2430 '''
2431 grid = np.array(np.meshgrid(*[[0., 0.5]] * self.dim, indexing='ij')).reshape(self.dim, -1).T
2432 return [f @ self.rec_vec_k for f in grid]
2434 def _operator(self, op, k: NDArray[np.float64]) -> NDArray[np.complex128]:
2435 '''
2436 Private method. Matrix of a symmetry operator at k (op is a matrix, or a callable of k).
2437 '''
2438 mat = np.asarray(op(k) if callable(op) else op, dtype='c16')
2439 error_handling.operator(mat, self.norb)
2440 return mat
2442 def parity_z2(self, inversion, bands: list[int]) -> int:
2443 r'''
2444 Get the :math:`\mathbb{Z}_2` invariant of an inversion- and
2445 time-reversal-symmetric model from the parity eigenvalues of its
2446 occupied bands at the time-reversal-invariant momenta
2447 :math:`\Gamma_i` (Fu and Kane, Phys. Rev. B 76, 045302 (2007)):
2449 .. math::
2451 (-1)^\nu = \prod_i \delta_i\, ,\qquad
2452 \delta_i = \prod_{m=1}^{N} \xi_{2m}(\Gamma_i)\, ,
2454 one parity :math:`\xi = \pm1` per Kramers pair. In 3D, :math:`\nu`
2455 is the strong index.
2457 :param inversion: Complex ndarray, shape (norb, norb), or callable of
2458 k returning one: the inversion operator :math:`P`, with
2459 :math:`P H(\mathbf{k}) P^{-1} = H(-\mathbf{k})` (e.g. swapping the
2460 two sublattices of graphene, and the identity on spin).
2461 :param bands: List of band indices (the occupied Kramers pairs).
2463 :returns:
2464 * **nu** -- 0 (trivial) or 1 (topological).
2465 '''
2466 error_handling.even_bands(bands)
2467 error_handling.band_indices(bands, self.norb)
2468 product = 1
2469 for k in self._trims():
2470 v = self._subspace(k, bands)
2471 xi = np.linalg.eigvals(v.conj().T @ self._operator(inversion, k) @ v)
2472 error_handling.parities(xi)
2473 n_minus = int(np.sum(xi.real < 0))
2474 product *= (-1) ** (n_minus // 2)
2475 return 0 if product == 1 else 1
2477 def symmetry_error(
2478 self, op, k_map: str = 'minus', antiunitary: bool = False, anti: bool = False,
2479 nk: int = 5, seed: int = 0,
2480 ) -> float:
2481 r'''
2482 Check a symmetry of the Bloch Hamiltonian,
2484 .. math::
2486 H(g\mathbf{k}) = \pm\, U\, \mathcal{K}H(\mathbf{k})\mathcal{K}^{-1} U^\dagger\, ,
2488 with :math:`g\mathbf{k} = -\mathbf{k}` or :math:`\mathbf{k}`,
2489 :math:`\mathcal{K}` complex conjugation if *antiunitary*, and the
2490 minus sign if *anti*: time reversal is antiunitary with
2491 :math:`\mathbf{k}\to-\mathbf{k}` and the plus sign, particle-hole
2492 antiunitary with :math:`\mathbf{k}\to-\mathbf{k}` and the minus sign,
2493 chiral unitary with :math:`\mathbf{k}\to\mathbf{k}` and the minus sign,
2494 inversion unitary with :math:`\mathbf{k}\to-\mathbf{k}` and the plus
2495 sign.
2497 :param op: Complex ndarray, shape (norb, norb), or callable of k: :math:`U`.
2498 :param k_map: 'minus' or 'identity'. Default value 'minus'.
2499 :param antiunitary: Boolean. Default value False.
2500 :param anti: Boolean. Default value False.
2501 :param nk: Positive integer. Default value 5. Points per direction of
2502 the (randomly shifted) test mesh.
2503 :param seed: Integer. Default value 0. Seed of the random shift.
2505 :returns:
2506 * **error** -- Real number. Largest deviation over the mesh (0 for a symmetry).
2507 '''
2508 error_handling.k_map(k_map)
2509 error_handling.boolean(antiunitary, 'antiunitary')
2510 error_handling.boolean(anti, 'anti')
2511 error_handling.positive_int(nk, 'nk')
2512 rng = np.random.default_rng(seed)
2513 fracs = np.array(np.meshgrid(*[np.arange(nk) / nk] * self.dim, indexing='ij')).reshape(self.dim, -1).T
2514 fracs = fracs + rng.uniform(0., 1. / nk, self.dim)
2515 sign = -1. if anti else 1.
2516 err = 0.
2517 for f in fracs:
2518 k = f @ self.rec_vec_k
2519 gk = -k if k_map == 'minus' else k
2520 ham = self.get_ham(k)
2521 u = self._operator(op, k)
2522 rhs = u @ (ham.conj() if antiunitary else ham) @ u.conj().T
2523 err = max(err, float(np.max(np.abs(self.get_ham(gk) - sign * rhs))))
2524 return err
2526 def tenfold_class(
2527 self, time_reversal=None, particle_hole=None, chiral=None, nk: int = 5, tol: float = 1e-8,
2528 ) -> str:
2529 r'''
2530 Get the Altland-Zirnbauer symmetry class (the "tenfold way") of the
2531 model from its time-reversal :math:`T = U_T\mathcal{K}`, particle-hole
2532 :math:`C = U_C\mathcal{K}` and chiral :math:`S = U_S` symmetries. Each
2533 given operator is first checked (see *symmetry_error*) to be a
2534 symmetry; the class then follows from :math:`T^2 = U_TU_T^* = \pm1`,
2535 :math:`C^2 = \pm1`, and the presence of :math:`S`. Any two of the
2536 three imply the third (:math:`S = TC`), which need not be given.
2538 :param time_reversal: Matrix :math:`U_T`, or None (absent).
2539 :param particle_hole: Matrix :math:`U_C`, or None (absent).
2540 :param chiral: Matrix :math:`U_S`, or None (absent).
2541 :param nk: Positive integer. Default value 5. See *symmetry_error*.
2542 :param tol: Positive real. Default value 1e-8. Tolerance of the checks.
2544 :returns:
2545 * **name** -- String: 'A', 'AIII', 'AI', 'BDI', 'D', 'DIII', 'AII', 'CII', 'C' or 'CI'.
2546 '''
2547 error_handling.positive_real(tol, 'tol')
2548 k0 = np.zeros(self.dim)
2549 for op, k_map, antiu, anti in ((time_reversal, 'minus', True, False),
2550 (particle_hole, 'minus', True, True),
2551 (chiral, 'identity', False, True)):
2552 if op is not None:
2553 error_handling.is_symmetry(self.symmetry_error(op, k_map, antiu, anti, nk), tol)
2554 u_t = None if time_reversal is None else self._operator(time_reversal, k0)
2555 u_c = None if particle_hole is None else self._operator(particle_hole, k0)
2556 # two of the three imply the third: C = S T, T = S C
2557 if chiral is not None and (u_t is None) != (u_c is None):
2558 u_s = self._operator(chiral, k0)
2559 if u_t is not None:
2560 u_c = u_s @ u_t
2561 elif u_c is not None:
2562 u_t = u_s @ u_c
2563 squares = []
2564 for u in (u_t, u_c):
2565 if u is None:
2566 squares.append(0)
2567 continue
2568 square = u @ u.conj()
2569 error_handling.square(square, tol)
2570 squares.append(int(np.sign(square[0, 0].real)))
2571 has_s = chiral is not None or (u_t is not None and u_c is not None)
2572 table = {(0, 0, False): 'A', (0, 0, True): 'AIII', (1, 0, False): 'AI',
2573 (1, 1, True): 'BDI', (0, 1, False): 'D', (-1, 1, True): 'DIII',
2574 (-1, 0, False): 'AII', (-1, -1, True): 'CII', (0, -1, False): 'C',
2575 (1, -1, True): 'CI'}
2576 return table[(squares[0], squares[1], has_s)]
2578 def quantum_geometric_tensor(
2579 self, bands: int | list[int], k: ArrayLike, dk: float = 1e-4, positions: bool = True,
2580 ) -> NDArray[np.complex128]:
2581 r'''
2582 Get the quantum geometric tensor of a group of bands at :math:`\mathbf{k}`,
2584 .. math::
2586 Q_{\mu\nu} = \mathrm{Tr}\left[P\,\partial_\mu P\,\partial_\nu P\right]
2587 = \sum_{n\in\mathrm{bands}}\langle\partial_\mu u_n|(1-P)|\partial_\nu u_n\rangle\, ,
2589 with :math:`P(\mathbf{k})` the projector on the bands (gauge
2590 invariant, so degeneracies within the group are harmless). Its real
2591 part is the quantum metric :math:`g_{\mu\nu}` (Provost and Vallee,
2592 1980), minus twice its imaginary part the Berry curvature
2593 :math:`\Omega_{\mu\nu} = -2\,\mathrm{Im}\,Q_{\mu\nu}` (the convention of
2594 *berry_curvature*, per unit k-area here). They obey
2595 :math:`\sqrt{\det g} \ge |\Omega_{xy}|/2` in 2D.
2597 :param bands: Band index, or list of band indices.
2598 :param k: k point (see *get_ham*).
2599 :param dk: Positive real. Default value 1e-4. Finite-difference step.
2600 :param positions: Boolean. Default value True. Include the orbital
2601 positions in the Bloch functions. The local metric and curvature
2602 depend on this choice (the Berry curvature of *berry_curvature* is
2603 that of ``positions=False``); integrals such as the Chern number
2604 do not.
2606 :returns:
2607 * **Q** -- Complex ndarray, shape (dim, dim), in the k coordinates of *get_ham*.
2608 '''
2609 if isinstance(bands, int):
2610 bands = [bands]
2611 error_handling.band_indices(bands, self.norb)
2612 error_handling.k_vector(k, 'k', self.dim)
2613 error_handling.positive_real(dk, 'dk')
2614 error_handling.boolean(positions, 'positions')
2615 k = np.asarray(k, dtype='f8')
2616 tau = self.orbital_positions()
2618 def proj(kk):
2619 v = self._subspace(kk, bands)
2620 if positions:
2621 v = np.exp(-1j * tau @ (self.k_basis @ kk))[:, None] * v
2622 return v @ v.conj().T
2624 p0 = proj(k)
2625 dp = []
2626 for mu in range(self.dim):
2627 e = np.zeros(self.dim)
2628 e[mu] = dk
2629 dp.append((proj(k + e) - proj(k - e)) / (2 * dk))
2630 return np.array([[np.trace(p0 @ dp[mu] @ dp[nu]) for nu in range(self.dim)]
2631 for mu in range(self.dim)])
2633 def surface_spectral_function(
2634 self, ks: ArrayLike, energies: ArrayLike, direction: int, side: int = 1,
2635 eta: float = 1e-2, bulk: bool = False, max_iter: int = 10000,
2636 ) -> NDArray[np.float64]:
2637 r'''
2638 Get the spectral function of the surface of a semi-infinite crystal,
2640 .. math::
2642 A_s(\mathbf{k}_\parallel, E) = -\frac{1}{\pi}\,
2643 \mathrm{Im}\,\mathrm{Tr}\, G_{00}(\mathbf{k}_\parallel, E + i\eta)\, ,
2645 with :math:`G_{00}` the Green's function of the outermost unit cell:
2646 the crystal fills the cells :math:`n\,\mathbf{a}_d`,
2647 :math:`n = 0, 1, 2, \dots` along the primitive vector
2648 :math:`\mathbf{a}_d` (``lat.prim_vec[direction]``), times *side*, and
2649 is periodic along the others. It comes from the Lopez Sancho-Rubio
2650 decimation (see *transport.surface_green*) of the chain of principal
2651 layers (as many cells as the longest hopping along
2652 :math:`\mathbf{a}_d` spans), so the crystal is truly semi-infinite: no
2653 finite slab, and no states of the opposite surface. It is the
2654 tight-binding picture of an ARPES map -- surface states (edge states,
2655 in 2D; the end state of a 1D chain) show up as sharp lines inside the
2656 bulk gaps. With *bulk*, the same quantity for a unit cell deep inside
2657 the crystal: the bulk bands projected on the surface Brillouin zone.
2659 :param ks: Real array, shape (nk, dim). k-points, in the coordinates
2660 of *get_ham* (e.g. from a *k_path* of the surface Brillouin zone).
2661 Only :math:`\mathbf{k}_\parallel` matters: the result does not
2662 depend on the component along the reciprocal vector
2663 :math:`\mathbf{b}_d`.
2664 :param energies: Real number or real array. Energies :math:`E`.
2665 :param direction: Integer, between 0 and dim-1. The primitive vector
2666 normal to the surface (the one along which the crystal is
2667 semi-infinite).
2668 :param side: +1 or -1. Default value 1. The crystal extends along
2669 :math:`+\mathbf{a}_d` (the surface faces :math:`-\mathbf{a}_d`) or
2670 along :math:`-\mathbf{a}_d` (the opposite surface).
2671 :param eta: Positive real. Default value 1e-2. Broadening: the
2672 width of the lines, and the decay length (about :math:`v/\eta`
2673 cells) over which the bulk is resolved.
2674 :param bulk: Boolean. Default value False. If True, the spectral
2675 function of a bulk unit cell instead of the surface one.
2676 :param max_iter: Positive integer. Default value 10000. Maximum
2677 number of decimation steps.
2679 :returns:
2680 * **A** -- Real ndarray, shape (nk, len(energies)).
2682 Example usage::
2684 # the (001) surface of a 3D model, along a path of the surface zone
2685 A = ks.surface_spectral_function(k_path, np.linspace(-1, 1, 201), direction=2)
2686 '''
2687 error_handling.no_overlap(self._overlap_hop)
2688 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
2689 error_handling.ks(ks, self.dim)
2690 error_handling.frequencies(energies, 'energies')
2691 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
2692 error_handling.direction(direction, self.dim)
2693 error_handling.surface_side(side)
2694 error_handling.positive_real(eta, 'eta')
2695 error_handling.boolean(bulk, 'bulk')
2696 error_handling.positive_int(max_iter, 'max_iter')
2697 # H_m(k): the hoppings from a cell to the cell m a_d away, with the
2698 # phase of their in-plane part only (independent of k along b_d)
2699 a_dir = np.array(self.lat.prim_vec[direction], dtype='f8')
2700 layers: dict[int, list] = {}
2701 for (i, j, cell, t), (_, _, R, _) in zip(self._hop_cells(), self._hop):
2702 layers.setdefault(cell[direction], []).append((i, j, R - cell[direction] * a_dir, t))
2703 k_cart = ks @ self.k_basis.T
2704 blocks: dict[int, NDArray] = {m: self._bloch_sum(hops, k_cart)[0] for m, hops in layers.items()}
2705 blocks[0] = blocks.get(0, 0.) + (np.diag(self.onsite) + self._onsite_offdiag)[None]
2706 width = max([abs(m) for m in layers] + [1])
2707 norb, nk = self.norb, len(ks)
2709 def layer(shift):
2710 # block (p, q) couples cell p of a principal layer to cell q of
2711 # the layer *shift* cells further into the crystal
2712 out = np.zeros((nk, width, norb, width, norb), 'c16')
2713 for p in range(width):
2714 for q in range(width):
2715 m = side * (q - p + shift)
2716 if m in blocks:
2717 out[:, p, :, q, :] = blocks[m]
2718 return out.reshape(nk, width * norb, width * norb)
2719 h0, alpha, beta = layer(0), layer(width), layer(-width)
2720 spec = np.empty((nk, len(energies)))
2721 for n in range(nk):
2722 gs, gb = _green_surface_bulk(energies + 1j * eta, h0[n], alpha[n], beta[n], 1e-12,
2723 max_iter, 'surface_spectral_function')
2724 g = gb if bulk else gs
2725 spec[n] = -np.trace(g[:, :norb, :norb], axis1=1, axis2=2).imag / PI
2726 return spec
2728 # ------------------------------------------------------------------
2729 # Projections, constant-energy contours and plot helpers
2730 # ------------------------------------------------------------------
2732 def spin_operator(self, axis: str) -> NDArray[np.complex128]:
2733 r'''
2734 Get the Pauli matrix :math:`\sigma_{axis}` acting on the spin of
2735 every site, :math:`\mathbb{1}_{sites}\otimes\sigma_{axis}`, for a model
2736 built with ``spin=True`` (for *band_weights* and *spin_texture*).
2738 :param axis: 'x', 'y' or 'z'.
2740 :returns:
2741 * **op** -- Complex ndarray, shape (norb, norb).
2742 '''
2743 error_handling.spinful(self.spin)
2744 error_handling.spin_axis(axis)
2745 return np.kron(np.eye(self.n_sites), PAULI[axis])
2747 def _projector(self, projector) -> NDArray[np.complex128]:
2748 '''
2749 Private method. The (norb, norb) operator of *band_weights* from an
2750 orbital index, a list of them, a sublattice tag or a matrix.
2751 '''
2752 error_handling.projector(projector, self.norb, self.tags)
2753 if isinstance(projector, np.ndarray):
2754 return projector.astype('c16')
2755 if isinstance(projector, str):
2756 sites = np.nonzero(self.tags == projector)[0]
2757 projector = np.concatenate([2 * sites, 2 * sites + 1]) if self.spin else sites
2758 diag = np.zeros(self.norb)
2759 diag[projector] = 1.
2760 return np.diag(diag).astype('c16')
2762 def _expectations(self, ks: NDArray[np.float64], ops: list) -> NDArray[np.float64]:
2763 r'''
2764 Private method. :math:`\mathrm{Re}\langle u_n(\mathbf{k})|O|u_n(\mathbf{k})\rangle
2765 /\langle u_n|u_n\rangle` for each operator of *ops*, band and k-point:
2766 shape (len(ops), nk, norb), from one diagonalization. Within a group
2767 of degenerate bands (energies within 1e-8 of each other, relative to
2768 the largest one), each band gets the average over the group: the
2769 trace of O over the degenerate subspace divided by its dimension,
2770 which does not depend on the eigenvectors the solver returns.
2771 '''
2772 en, vn = self._eigs(ks, eigenvec=True)
2773 norm = np.sum(np.abs(vn) ** 2, axis=1)
2774 values = np.array([np.einsum('kin,ij,kjn->kn', vn.conj(), op, vn).real / norm for op in ops])
2775 # label the degenerate groups (consecutive bands, sorted by energy)
2776 tol = 1e-8 * max(1., np.abs(en).max())
2777 new_group = np.concatenate([np.ones((len(en), 1), bool),
2778 np.abs(np.diff(en, axis=1)) > tol], axis=1)
2779 groups = np.cumsum(new_group.ravel()) - 1
2780 counts = np.bincount(groups)
2781 return np.array([(np.bincount(groups, v.ravel()) / counts)[groups].reshape(v.shape)
2782 for v in values])
2784 def band_weights(self, projector, ks: ArrayLike | None = None) -> NDArray[np.float64]:
2785 r'''
2786 Get the weight of each band on a set of orbitals, or the expectation
2787 value of an operator, for "fat band" plots (see *plot_bands*):
2789 .. math::
2791 w_n(\mathbf{k}) = \langle u_n(\mathbf{k})|P|u_n(\mathbf{k})\rangle\, ,
2793 with :math:`P` the projector on the chosen orbitals (a weight
2794 between 0 and 1; the weights on all the sublattices add up to 1), or
2795 any Hermitian operator, e.g. ``ks.spin_operator('z')`` for the spin
2796 polarization (between -1 and 1). Within a group of degenerate bands
2797 (e.g. Kramers pairs with inversion and time reversal), the weight of
2798 a single eigenvector depends on the choice of eigenvectors; each band
2799 of the group gets instead the average over the group, which does not
2800 (a Kramers pair has zero spin polarization).
2801 For a non-Hermitian model, the right eigenvectors are used; with an
2802 overlap (*set_overlap*), :math:`v^\dagger Pv/v^\dagger v`.
2804 :param projector: Orbital index, list of orbital indices, sublattice
2805 tag (all the orbitals, both spins, of the sites with that tag),
2806 or a (norb, norb) Hermitian matrix.
2807 :param ks: Real array, shape (nk, dim). Default value None: the
2808 k-points of the last *get_bands* or *k_path*.
2810 :returns:
2811 * **weights** -- Real ndarray, shape (nk, norb), ordered as the bands.
2813 Example usage::
2815 ks.k_path(points, nk=100)
2816 ks.plot_bands(weights=ks.band_weights('a'))
2817 '''
2818 op = self._projector(projector)
2819 if ks is None:
2820 error_handling.empty_ndarray(self.en, 'get_bands or k_path')
2821 ks = self.ks
2822 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
2823 error_handling.ks(ks, self.dim)
2824 return self._expectations(ks, [op])[0]
2826 def spin_texture(self, ks: ArrayLike, band: int) -> NDArray[np.float64]:
2827 r'''
2828 Get the spin expectation values
2829 :math:`\langle\boldsymbol\sigma\rangle_n(\mathbf{k}) =
2830 \langle u_n|\boldsymbol\sigma|u_n\rangle` of a band, for a model built
2831 with ``spin=True``. Their length is 1 if the band is fully polarized,
2832 less if spin-orbit coupling entangles the spin with the orbitals.
2833 Within a group of degenerate bands, each band gets the average spin
2834 of the group (see *band_weights*): zero for a Kramers pair, with both
2835 inversion and time reversal.
2837 :param ks: Real array, shape (nk, dim). k-points.
2838 :param band: Integer. Band index.
2840 :returns:
2841 * **spin** -- Real ndarray, shape (nk, 3):
2842 :math:`(\langle\sigma_x\rangle, \langle\sigma_y\rangle, \langle\sigma_z\rangle)`.
2843 '''
2844 error_handling.spinful(self.spin)
2845 error_handling.band_index(band, self.norb)
2846 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
2847 error_handling.ks(ks, self.dim)
2848 ops = [self.spin_operator(a) for a in 'xyz']
2849 return self._expectations(ks, ops)[:, :, band].T
2851 def _fold(self, elements: NDArray[np.float64]) -> NDArray[np.float64]:
2852 '''
2853 Private method. Bring segments or triangles (fractional coordinates,
2854 shape (n, nv, dim)) into the first Brillouin zone (the Wigner-Seitz
2855 cell of the reciprocal lattice, by their centres), in the k
2856 coordinates of *get_ham*.
2857 '''
2858 elements = elements - np.floor(elements.mean(axis=1) + 0.5)[:, None, :]
2859 cart = elements @ self.rec_vec_k
2860 shifts = np.array(list(product((-1, 0, 1), repeat=self.dim))) @ self.rec_vec_k
2861 dist = np.linalg.norm(cart.mean(axis=1)[:, None, :] - shifts[None], axis=2)
2862 return cart - shifts[np.argmin(dist, axis=1)][:, None, :]
2864 def fermi_surface(
2865 self, energy: float = 0., nk: int | tuple[int, ...] | None = None,
2866 bands: int | list[int] | None = None,
2867 ) -> list[NDArray[np.float64]]:
2868 r'''
2869 Get the constant-energy surfaces :math:`E_n(\mathbf{k}) = E` (the
2870 Fermi surface at :math:`E = E_F`): lines in a 2D Brillouin zone,
2871 surfaces in 3D. The bands are interpolated linearly in the triangles
2872 (tetrahedra in 3D) of a uniform mesh, the same decomposition as
2873 *dos.tetrahedron_dos*, and the level set of each simplex is a segment
2874 (a triangle or a quadrilateral, split in two, in 3D): marching
2875 triangles and tetrahedra, exact for the interpolated bands. Each
2876 piece is brought into the first Brillouin zone, the Wigner-Seitz
2877 cell around :math:`\Gamma`. For a non-Hermitian model, the real part
2878 of the energies is used.
2880 :param energy: Real number. Default value 0. The energy :math:`E`.
2881 :param nk: Positive integer, or tuple of *dim* positive integers.
2882 Default value None: 100 in 2D, 40 in 3D. Number of k-points along
2883 each reciprocal lattice vector.
2884 :param bands: Band index, or list of band indices. Default value
2885 None: all bands.
2887 :returns:
2888 * **surfaces** -- List, one entry per band of *bands*: a real
2889 ndarray of shape (n, 2, 2) in 2D (n segments, each two k-points)
2890 or (n, 3, 3) in 3D (n triangles, each three k-points), in the
2891 k coordinates of *get_ham*. Empty (n = 0) if the band does not
2892 cross *energy*.
2894 Example usage::
2896 # the Fermi surface of the half-filled square lattice
2897 segments = sq.fermi_surface(0.)[0]
2898 '''
2899 error_handling.dim_min(self.dim, 2)
2900 error_handling.real_number(energy, 'energy')
2901 if nk is None:
2902 nk = 100 if self.dim == 2 else 40
2903 error_handling.nk(nk, self.dim)
2904 if isinstance(nk, int):
2905 nk = (nk,) * self.dim
2906 if bands is None:
2907 bands = list(range(self.norb))
2908 elif isinstance(bands, int):
2909 bands = [bands]
2910 error_handling.band_indices(bands, self.norb)
2911 _, ks = self.mesh_grid(nk)
2912 en = self._eigs(ks).real.reshape(*nk, self.norb)
2913 return [self._fold(_level_set(en[..., n], energy)) for n in bands]
2915 def _k_labels(self) -> list[str]:
2916 '''
2917 Private method. Axis labels of the k coordinates of *get_ham*:
2918 '$k_x$' for a coordinate along the Cartesian x axis (the slab of a
2919 3D model along x and z gives '$k_x$', '$k_z$'), '$k_1$', ... otherwise.
2920 '''
2921 labels = []
2922 for n, col in enumerate(self.k_basis.T):
2923 axis = np.nonzero(np.isclose(np.abs(col), 1.))[0]
2924 labels.append('$k_{}$'.format('xyz'[axis[0]] if len(axis) else n + 1))
2925 return labels
2927 def _zone_axes(self, ax, fs: float) -> None:
2928 '''
2929 Private method. Draw the first Brillouin zone of a 2D model on *ax*,
2930 and label and scale its axes.
2931 '''
2932 zone = _zone_polygon(self.rec_vec_k)
2933 ax.plot(*np.vstack([zone, zone[:1]]).T, 'k', lw=1)
2934 ax.set_aspect('equal')
2935 pad = 0.05 * np.ptp(zone, axis=0).max()
2936 ax.set_xlim(zone[:, 0].min() - pad, zone[:, 0].max() + pad)
2937 ax.set_ylim(zone[:, 1].min() - pad, zone[:, 1].max() + pad)
2938 ax.set_xlabel(self._k_labels()[0], fontsize=fs)
2939 ax.set_ylabel(self._k_labels()[1], fontsize=fs)
2941 def plot_fermi_surface(
2942 self, energy: float = 0., nk: int | tuple[int, ...] | None = None,
2943 bands: int | list[int] | None = None, lw: float = 2., alpha: float = 0.6,
2944 fs: float = 20, figsize: tuple[float, float] | None = None,
2945 ) -> Figure:
2946 '''
2947 Plot the constant-energy contours of *fermi_surface* (one color per
2948 band, with a legend if several bands cross the energy): lines inside
2949 the first Brillouin zone (drawn in black) in 2D, surfaces in 3D.
2951 :param energy: Real number. Default value 0. The energy.
2952 :param nk: See *fermi_surface*.
2953 :param bands: See *fermi_surface*.
2954 :param lw: Positive number. Default value 2. Linewidth (2D).
2955 :param alpha: Real number in (0, 1]. Default value 0.6. Opacity of
2956 the surfaces (3D).
2957 :param fs: Positive number. Default value 20. Fontsize.
2958 :param figsize: Tuple. Default value None. Figure size.
2960 :returns:
2961 * **fig** -- Figure.
2962 '''
2963 error_handling.positive_real(lw, 'lw')
2964 error_handling.positive_real(alpha, 'alpha')
2965 error_handling.positive_real(fs, 'fs')
2966 error_handling.tuple_2elem(figsize, 'figsize')
2967 surfaces = self.fermi_surface(energy, nk, bands)
2968 bands = list(range(self.norb)) if bands is None else np.atleast_1d(bands).tolist()
2969 fig = plt.figure(figsize=figsize)
2970 if self.dim == 2:
2971 ax = fig.add_subplot()
2972 for n, (band, segs) in enumerate(zip(bands, surfaces)):
2973 if len(segs):
2974 ax.add_collection(LineCollection(list(segs), colors='C{}'.format(n % 10), lw=lw,
2975 label='band {}'.format(band)))
2976 self._zone_axes(ax, fs)
2977 else:
2978 ax = fig.add_subplot(projection='3d')
2979 for n, (band, tris) in enumerate(zip(bands, surfaces)):
2980 if len(tris):
2981 ax.add_collection3d(Poly3DCollection(tris, facecolor='C{}'.format(n % 10),
2982 alpha=min(alpha, 1.), edgecolor='none',
2983 label='band {}'.format(band)))
2984 extent = np.abs(np.concatenate([t.reshape(-1, 3) for t in surfaces] +
2985 [self.rec_vec_k / 2])).max()
2986 ax.set(xlim=(-extent, extent), ylim=(-extent, extent), zlim=(-extent, extent))
2987 ax.set_box_aspect((1, 1, 1))
2988 ax.set_xlabel(self._k_labels()[0], fontsize=fs)
2989 ax.set_ylabel(self._k_labels()[1], fontsize=fs)
2990 ax.set_zlabel(self._k_labels()[2], fontsize=fs)
2991 if sum(len(s) > 0 for s in surfaces) > 1:
2992 ax.legend(fontsize=0.6 * fs)
2993 ax.set_title('$E = {:g}$'.format(energy), fontsize=fs)
2994 plt.draw()
2995 return fig
2997 def plot_spin_texture(
2998 self, band: int, energy: float = 0., nk: int = 100, n_arrows: int = 40,
2999 cmap: str = 'RdBu_r', lw: float = 1., fs: float = 20,
3000 figsize: tuple[float, float] | None = None,
3001 ) -> Figure:
3002 r'''
3003 Plot the spin texture of a band of a 2D spinful model on its
3004 constant-energy contour (see *fermi_surface*), zoomed on the
3005 contour: arrows of
3006 :math:`(\langle\sigma_x\rangle, \langle\sigma_y\rangle)` along the
3007 contour, colored by :math:`\langle\sigma_z\rangle` (see *spin_texture*).
3008 Rashba coupling, for instance, locks the spin perpendicular to
3009 :math:`\mathbf{k}`, winding once around the contour.
3011 :param band: Integer. Band index.
3012 :param energy: Real number. Default value 0. Energy of the contour.
3013 :param nk: Positive integer. Default value 100. Mesh of the contour.
3014 :param n_arrows: Positive integer. Default value 40. Approximate
3015 number of arrows, evenly spaced along the contour. A fully
3016 polarized in-plane spin is drawn 80% as long as their spacing.
3017 :param cmap: Default value 'RdBu_r'. Colormap of :math:`\langle\sigma_z\rangle`.
3018 :param lw: Positive number. Default value 1. Linewidth of the contour.
3019 :param fs: Positive number. Default value 20. Fontsize.
3020 :param figsize: Tuple. Default value None. Figure size.
3022 :returns:
3023 * **fig** -- Figure.
3024 '''
3025 error_handling.dim_2(self.dim)
3026 error_handling.spinful(self.spin)
3027 error_handling.band_index(band, self.norb)
3028 error_handling.positive_int(nk, 'nk')
3029 error_handling.positive_int(n_arrows, 'n_arrows')
3030 error_handling.positive_real(lw, 'lw')
3031 error_handling.positive_real(fs, 'fs')
3032 error_handling.tuple_2elem(figsize, 'figsize')
3033 segs = self.fermi_surface(energy, nk, [band])[0]
3034 error_handling.contour_found(len(segs), energy)
3035 # arrows at segment midpoints, at least one contour length / n_arrows apart
3036 mids = segs.mean(axis=1)
3037 spacing = np.linalg.norm(segs[:, 1] - segs[:, 0], axis=1).sum() / n_arrows
3038 keep = [0]
3039 for i in range(1, len(mids)):
3040 if np.min(np.linalg.norm(mids[keep] - mids[i], axis=1)) >= spacing:
3041 keep.append(i)
3042 spin = self.spin_texture(mids[keep], band)
3043 fig, ax = plt.subplots(figsize=figsize)
3044 ax.add_collection(LineCollection(list(segs), colors='0.6', lw=lw))
3045 arrows = ax.quiver(mids[keep, 0], mids[keep, 1], spin[:, 0], spin[:, 1], spin[:, 2],
3046 cmap=cmap, norm=Normalize(-1., 1.), pivot='mid', edgecolor='k',
3047 linewidth=0.5, angles='xy', scale_units='xy', scale=1.25 / spacing)
3048 fig.colorbar(arrows, ax=ax).set_label(r'$\langle\sigma_z\rangle$', fontsize=fs)
3049 self._zone_axes(ax, fs)
3050 # zoom on the contour (within the zone)
3051 low, high = segs.reshape(-1, 2).min(axis=0), segs.reshape(-1, 2).max(axis=0)
3052 pad = 0.1 * (high - low).max() + spacing
3053 ax.set_xlim(max(ax.get_xlim()[0], low[0] - pad), min(ax.get_xlim()[1], high[0] + pad))
3054 ax.set_ylim(max(ax.get_ylim()[0], low[1] - pad), min(ax.get_ylim()[1], high[1] + pad))
3055 ax.set_title('Spin texture of band {}, $E = {:g}$'.format(band, energy), fontsize=fs)
3056 plt.draw()
3057 return fig
3059 def plot_berry_curvature(
3060 self, bands: int | list[int], nk: int | tuple[int, int] = 60,
3061 plane: tuple[int, int] = (0, 1), k_fixed: float = 0., cmap: str = 'RdBu_r',
3062 fs: float = 20, figsize: tuple[float, float] | None = None,
3063 ) -> Figure:
3064 r'''
3065 Plot the Berry curvature of *berry_curvature* as a density,
3066 :math:`\Omega(\mathbf{k})` = flux of each plaquette / its area. In
3067 2D, the plaquettes are brought into the first Brillouin zone (drawn
3068 in black), in the k coordinates of *get_ham*. In 3D, they fill the
3069 parallelogram spanned by the two reciprocal vectors of *plane*, in
3070 orthonormal coordinates of the plane (:math:`k_1` along
3071 :math:`\mathbf{b}_i`). The title gives the Chern number, the total
3072 flux over :math:`2\pi`.
3074 :param bands: See *berry_curvature*.
3075 :param nk: Default value 60. See *berry_curvature*.
3076 :param plane: See *berry_curvature*.
3077 :param k_fixed: See *berry_curvature*.
3078 :param cmap: Default value 'RdBu_r'. Colormap, centred on zero.
3079 :param fs: Positive number. Default value 20. Fontsize.
3080 :param figsize: Tuple. Default value None. Figure size.
3082 :returns:
3083 * **fig** -- Figure.
3084 '''
3085 error_handling.positive_real(fs, 'fs')
3086 error_handling.tuple_2elem(figsize, 'figsize')
3087 curv = self.berry_curvature(bands, nk, plane, k_fixed)
3088 n1, n2 = curv.shape
3089 # the corners of each plaquette, counterclockwise, in fractional coordinates
3090 f1, f2 = np.meshgrid(np.arange(n1) / n1, np.arange(n2) / n2, indexing='ij')
3091 corner = np.stack([f1.ravel(), f2.ravel()], axis=-1)[:, None, :]
3092 quads = corner + np.array([[0., 0.], [1. / n1, 0.], [1. / n1, 1. / n2], [0., 1. / n2]])
3093 b_i, b_j = self.rec_vec_k[plane[0]], self.rec_vec_k[plane[1]]
3094 if self.dim == 2:
3095 quads = self._fold(quads)
3096 else:
3097 # orthonormal coordinates of the plane, the first along b_i
3098 u = b_i / np.linalg.norm(b_i)
3099 v = b_j - (b_j @ u) * u
3100 b_i, b_j = np.array([b_i @ u, 0.]), np.array([b_j @ u, b_j @ v / np.linalg.norm(v)])
3101 quads = quads @ np.array([b_i, b_j])
3102 area = abs(b_i[0] * b_j[1] - b_i[1] * b_j[0]) / (n1 * n2)
3103 omega = curv.ravel() / area
3104 vmax = np.abs(omega).max() or 1.
3105 fig, ax = plt.subplots(figsize=figsize)
3106 patches = PolyCollection(quads, array=omega, cmap=cmap, norm=Normalize(-vmax, vmax),
3107 edgecolors='face', linewidths=0.2)
3108 ax.add_collection(patches)
3109 fig.colorbar(patches, ax=ax).set_label(r'$\Omega(\mathbf{k})$', fontsize=fs)
3110 if self.dim == 2:
3111 self._zone_axes(ax, fs)
3112 else:
3113 ax.autoscale_view()
3114 ax.set_aspect('equal')
3115 ax.set_xlabel('$k_1$', fontsize=fs)
3116 ax.set_ylabel('$k_2$', fontsize=fs)
3117 ax.set_title('Berry curvature, $C = {:.3f}$'.format(curv.sum() / (2 * PI)), fontsize=fs)
3118 fig.set_layout_engine('tight')
3119 plt.draw()
3120 return fig
3122 def plot_wannier_flow(
3123 self, bands: int | list[int], nk: int = 100, nk_perp: int = 51, direction: int = 0,
3124 flow: int | None = None, k_fixed: float = 0., positions: bool = True,
3125 k_range: tuple[float, float] = (0., 1.), ms: float = 4., c: str = 'b', fs: float = 20,
3126 figsize: tuple[float, float] | None = None,
3127 ) -> Figure:
3128 r'''
3129 Plot the flow of the hybrid Wannier centres (the Wilson-loop
3130 spectrum) of *wannier_flow*: the centres along *direction*, in
3131 units of :math:`\mathbf{a}_{direction}`, as the loop moves along
3132 the reciprocal vector *flow*. They wind :math:`C` times in a Chern
3133 insulator, and the Kramers pairs swap partners in a
3134 :math:`\mathbb{Z}_2` topological insulator.
3136 :param bands: See *wannier_flow*.
3137 :param nk: See *wannier_flow*.
3138 :param nk_perp: See *wannier_flow*.
3139 :param direction: See *wannier_flow*.
3140 :param flow: See *wannier_flow*.
3141 :param k_fixed: See *wannier_flow*.
3142 :param positions: See *wannier_flow*.
3143 :param k_range: See *wannier_flow*.
3144 :param ms: Positive number. Default value 4. Marker size.
3145 :param c: Default value 'b'. Marker color.
3146 :param fs: Positive number. Default value 20. Fontsize.
3147 :param figsize: Tuple. Default value None. Figure size.
3149 :returns:
3150 * **fig** -- Figure.
3151 '''
3152 error_handling.positive_real(ms, 'ms')
3153 error_handling.positive_real(fs, 'fs')
3154 error_handling.tuple_2elem(figsize, 'figsize')
3155 fracs, centers = self.wannier_flow(bands, nk, nk_perp, direction, flow, k_fixed,
3156 positions, k_range)
3157 if flow is None:
3158 flow = [d for d in range(self.dim) if d != direction][0]
3159 fig, ax = plt.subplots(figsize=figsize)
3160 ax.plot(fracs, centers, 'o', c=c, ms=ms)
3161 ax.set_xlim(fracs[0], fracs[-1])
3162 ax.set_ylim(0., 1.)
3163 ax.set_xlabel(r'$k_{0}/|\mathbf{{b}}_{0}|$'.format(flow + 1), fontsize=fs)
3164 ax.set_ylabel(r'Wannier centres ($\mathbf{{a}}_{}$)'.format(direction + 1), fontsize=fs)
3165 fig.set_layout_engine('tight')
3166 plt.draw()
3167 return fig
3169 def plot_surface_spectral_function(
3170 self, points: list[ArrayLike], energies: ArrayLike, direction: int, nk: int = 60,
3171 side: int = 1, eta: float = 1e-2, bulk: bool = False,
3172 node_labels: list[str] | None = None, log: bool = True, cmap: str = 'magma',
3173 fs: float = 20, figsize: tuple[float, float] | None = None,
3174 ) -> Figure:
3175 r'''
3176 Plot the surface spectral function of *surface_spectral_function*
3177 along a path through the surface Brillouin zone (built as in
3178 *k_path*), as an ARPES-like map of :math:`E` against
3179 :math:`\mathbf{k}_\parallel`. Surface states are the sharp lines in
3180 the bulk gaps.
3182 :param points: List of at least two k-points, the nodes of the path
3183 (see *k_path*).
3184 :param energies: Real array. Energies :math:`E`.
3185 :param direction: See *surface_spectral_function*.
3186 :param nk: Positive integer. Default value 60. Number of k-points per path segment.
3187 :param side: See *surface_spectral_function*.
3188 :param eta: See *surface_spectral_function*.
3189 :param bulk: See *surface_spectral_function*.
3190 :param node_labels: List of strings. Default value None. Labels of *points*.
3191 :param log: Boolean. Default value True. Logarithmic color scale, over
3192 four decades below the maximum (the surface states are much
3193 sharper than the bulk continuum).
3194 :param cmap: Default value 'magma'. Colormap.
3195 :param fs: Positive number. Default value 20. Fontsize.
3196 :param figsize: Tuple. Default value None. Figure size.
3198 :returns:
3199 * **fig** -- Figure.
3200 '''
3201 error_handling.k_path_points(points, self.dim)
3202 error_handling.positive_int(nk, 'nk')
3203 error_handling.boolean(log, 'log')
3204 error_handling.positive_real(fs, 'fs')
3205 error_handling.tuple_2elem(figsize, 'figsize')
3206 if node_labels is not None:
3207 error_handling.ndarray(np.array(node_labels), 'node_labels', len(points))
3208 ks, dist, nodes = _path(points, nk)
3209 spec = self.surface_spectral_function(ks, energies, direction, side, eta, bulk)
3210 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
3211 vmax = spec.max() if spec.max() > 0 else 1.
3212 norm = LogNorm(vmax * 1e-4, vmax) if log else Normalize(0., vmax)
3213 fig, ax = plt.subplots(figsize=figsize)
3214 mesh = ax.pcolormesh(dist, energies, np.maximum(spec, vmax * 1e-4).T, cmap=cmap,
3215 norm=norm, shading='nearest')
3216 fig.colorbar(mesh, ax=ax).set_label(r'$A(\mathbf{k}_\parallel, E)$', fontsize=fs)
3217 for node in nodes:
3218 ax.axvline(node, color='w', lw=0.5)
3219 if node_labels is not None:
3220 ax.set_xticks(nodes)
3221 ax.set_xticklabels(node_labels, fontsize=fs)
3222 ax.set_ylabel('$E$', fontsize=fs)
3223 ax.set_title('Bulk' if bulk else 'Surface', fontsize=fs)
3224 fig.set_layout_engine('tight')
3225 plt.draw()
3226 return fig
3228 def plot_dos(
3229 self,
3230 nk: int | tuple[int, ...] = 30,
3231 broadening: float = 0.05,
3232 kernel: str = 'gaussian',
3233 e_grid: ArrayLike | None = None,
3234 fs: float = 20,
3235 lw: float = 2.,
3236 figsize: tuple[float, float] | None = None,
3237 ) -> Figure:
3238 '''
3239 Plot the density of states, obtained by diagonalizing
3240 :math:`H(\\mathbf{k})` over a uniform Brillouin-zone mesh: broadened
3241 (see *tbkit.dos.density_of_states*), or by the linear tetrahedron
3242 method (``kernel='tetrahedron'``, see *tbkit.dos.tetrahedron_dos*),
3243 which needs no broadening and keeps band edges and van Hove
3244 singularities sharp.
3246 :param nk: Positive integer, or tuple of *dim* positive integers.
3247 Default value 30. Number of k-points along each reciprocal
3248 lattice vector.
3249 :param broadening: Positive real number. Default value 0.05. Kernel
3250 width (unused by the tetrahedron method).
3251 :param kernel: String. Default value 'gaussian'. 'gaussian',
3252 'lorentzian' or 'tetrahedron'.
3253 :param e_grid: Real ndarray. Default value None. Energies at which to
3254 evaluate the density of states.
3255 :param fs: Positive number. Default value 20. Fontsize.
3256 :param lw: Positive number. Default value 2. Linewidth.
3257 :param figsize: Tuple. Default value None. Figure size.
3259 :returns:
3260 * **fig** -- Figure.
3261 '''
3262 error_handling.positive_real(fs, 'fs')
3263 error_handling.positive_real(lw, 'lw')
3264 error_handling.tuple_2elem(figsize, 'figsize')
3265 error_handling.kspace_dos_kernel(kernel)
3266 en = self.mesh_bands(nk)
3267 if kernel != 'tetrahedron':
3268 return dos._plot_density_of_states(en, e_grid, broadening, kernel, fs, lw, figsize)
3269 nk = (nk,) * self.dim if isinstance(nk, int) else nk
3270 e_grid, rho = dos.tetrahedron_dos(en.reshape(*nk, self.norb), e_grid)
3271 return dos._plot_dos_curve(e_grid, rho, fs, lw, figsize)
3273 def plot_bands(
3274 self,
3275 node_labels: list[str] | None = None,
3276 fs: float = 20,
3277 lw: float = 2.,
3278 ms: float = 0.,
3279 c: str = 'b',
3280 lims: tuple[float, float] | None = None,
3281 figsize: tuple[float, float] | None = None,
3282 weights: ArrayLike | None = None,
3283 style: str = 'color',
3284 cmap: str | None = None,
3285 ) -> Figure:
3286 '''
3287 Plot the band structure computed by *k_path* or *get_bands*, against
3288 the cumulative distance along the k-points (vertical lines mark the
3289 nodes of a *k_path*). For a non-Hermitian model, the real part of
3290 the energies is plotted. With *weights* (e.g. from *band_weights*:
3291 an orbital, sublattice or spin projection), a "fat band" plot: the
3292 bands are colored by the weights (with a colorbar), or drawn with
3293 dots sized by their absolute value.
3295 :param node_labels: List of strings. Default value None. Labels of the
3296 high-symmetry points passed to *k_path*.
3297 :param fs: Positive number. Default value 20. Fontsize.
3298 :param lw: Positive number. Default value 2. Linewidth.
3299 :param ms: Positive number. Default value 0. Marker size.
3300 :param c: Default value 'b'. Line color.
3301 :param lims: List. Default value None. Energy plot limits.
3302 :param figsize: Tuple. Default value None. Figure size.
3303 :param weights: Real array, shape (nk, norb). Default value None.
3304 One weight per k-point and band, e.g. from *band_weights*.
3305 :param style: String. Default value 'color'. 'color': the bands
3306 colored by *weights*; 'size': dots of size proportional to
3307 :math:`|w|` (*ms*, or 8 if *ms* is 0, for the largest) on thin
3308 lines of color *c*.
3309 :param cmap: Default value None ('viridis', or 'RdBu_r', centred on
3310 zero, if some weights are negative). Colormap of *style* 'color'.
3312 :returns:
3313 * **fig** -- Figure.
3314 '''
3315 error_handling.empty_ndarray(self.en, 'get_bands or k_path')
3316 error_handling.positive_real(fs, 'fs')
3317 error_handling.positive_real(lw, 'lw')
3318 error_handling.lims(lims)
3319 error_handling.tuple_2elem(figsize, 'figsize')
3320 fig, ax = plt.subplots(figsize=figsize)
3321 if weights is None:
3322 for n in range(self.norb):
3323 ax.plot(self.ks_dist, self.en[:, n].real, c=c, lw=lw, marker='o', ms=ms)
3324 else:
3325 self._plot_fat_bands(fig, ax, weights, style, cmap, c, lw, ms)
3326 for node in self.nodes:
3327 ax.axvline(node, color='k', lw=0.5)
3328 ax.set_xlim(self.ks_dist[0], self.ks_dist[-1])
3329 if lims is not None:
3330 ax.set_ylim(lims)
3331 if node_labels is not None:
3332 error_handling.ndarray(np.array(node_labels), 'node_labels', len(self.nodes))
3333 ax.set_xticks(self.nodes)
3334 ax.set_xticklabels(node_labels, fontsize=fs)
3335 ax.set_ylabel('$E$', fontsize=fs)
3336 for label in ax.yaxis.get_majorticklabels():
3337 label.set_fontsize(fs)
3338 fig.set_layout_engine('tight')
3339 plt.draw()
3340 return fig
3342 def _plot_fat_bands(self, fig, ax, weights, style, cmap, c, lw, ms) -> None:
3343 '''
3344 Private method. The bands of *plot_bands* colored (or with dots
3345 sized) by *weights*.
3346 '''
3347 error_handling.band_weights(weights, self.en.shape)
3348 error_handling.weight_style(style)
3349 weights = np.asarray(weights, dtype='f8')
3350 signed = weights.min() < 0
3351 vmax = np.abs(weights).max() or 1.
3352 if style == 'size':
3353 size = (ms or 8.) ** 2
3354 for n in range(self.norb):
3355 ax.plot(self.ks_dist, self.en[:, n].real, c=c, lw=lw / 4)
3356 ax.scatter(self.ks_dist, self.en[:, n].real, s=size * np.abs(weights[:, n]) / vmax,
3357 c=c, lw=0)
3358 return
3359 norm = Normalize(-vmax if signed else 0., vmax)
3360 cmap = cmap or ('RdBu_r' if signed else 'viridis')
3361 for n in range(self.norb):
3362 pts = np.stack([self.ks_dist, self.en[:, n].real], axis=-1)
3363 lines = LineCollection(list(np.stack([pts[:-1], pts[1:]], axis=1)), cmap=cmap, norm=norm, lw=lw)
3364 lines.set_array((weights[:-1, n] + weights[1:, n]) / 2)
3365 ax.add_collection(lines)
3366 ax.autoscale_view()
3367 fig.colorbar(lines, ax=ax)
3369 def show(self) -> None:
3370 '''
3371 Emulate Matplotlib method plt.show().
3372 '''
3373 plt.show()
3376def ribbon(
3377 lat: Lattice,
3378 list_hop: list[dict],
3379 width: int,
3380 direction: int = 1,
3381 onsite: dict | None = None,
3382 spin: bool = False,
3383) -> KSpace:
3384 r'''
3385 Cut a ribbon out of a 2D periodic model: periodic along one primitive
3386 vector, finite (open boundary, *width* unit cells) along the other.
3387 This is the standard way to see edge states in a band structure (e.g.
3388 the zero-energy edge band of a zigzag graphene ribbon, or the helical
3389 edge states of a Kane-Mele ribbon). Cut out of a 3D model, it gives a
3390 slab: periodic along the two other primitive vectors, the way to see
3391 surface states (e.g. the Fermi arcs of a Weyl semimetal).
3393 :param lat: **Lattice** class instance (2D or 3D, i.e. two or three primitive
3394 vectors). Only *unit_cell* and *prim_vec* are used.
3395 :param list_hop: List of dictionaries, in the same format passed to
3396 *KSpace.set_hopping* -- the hoppings of the periodic (2D) model
3397 that the ribbon is cut from.
3398 :param width: Positive integer. Number of unit cells across the ribbon.
3399 :param direction: 0 or 1 (or 2, in 3D). Default value 1. Which primitive
3400 vector (``lat.prim_vec[direction]``) becomes finite; the others stay
3401 periodic.
3402 :param onsite: Dictionary. Default value None. Onsite energies, in the
3403 same format passed to *KSpace.set_onsite* -- applied identically
3404 on every row of the ribbon.
3405 :param spin: Boolean. Default value False. See *KSpace*.
3407 :returns:
3408 * **rib** -- **KSpace** instance, 1D-periodic (2D-periodic for a slab), with
3409 ``width * len(lat.unit_cell)`` sites (each site of *lat*,
3410 repeated once per row across the ribbon; row *w*'s copy of site
3411 *i* is orbital ``w*len(lat.unit_cell) + i``).
3413 Example usage::
3415 # zigzag graphene ribbon, 20 unit cells wide
3416 list_hop = [{'i': 0, 'j': 1, 'R': (0, 0), 't': 1.},
3417 {'i': 0, 'j': 1, 'R': (-1, 0), 't': 1.},
3418 {'i': 0, 'j': 1, 'R': (0, -1), 't': 1.}]
3419 rib = ribbon(lat, list_hop, width=20)
3420 '''
3421 error_handling.lat(lat)
3422 error_handling.dim_min(len(lat.prim_vec), 2)
3423 error_handling.positive_int(width, 'width')
3424 error_handling.direction(direction, len(lat.prim_vec))
3425 periodic = [d for d in range(len(lat.prim_vec)) if d != direction]
3426 n_sites = len(lat.unit_cell)
3427 a_dir = np.array(lat.prim_vec[direction])
3428 new_unit_cell = []
3429 for w in range(width):
3430 for dic in lat.unit_cell:
3431 r0 = np.array(dic['r0']) + w*a_dir
3432 new_unit_cell.append({'tag': dic['tag'], 'r0': tuple(float(c) for c in r0)})
3433 new_lat = Lattice(unit_cell=new_unit_cell, prim_vec=[lat.prim_vec[d] for d in periodic])
3434 rib = KSpace(new_lat, spin=spin)
3435 new_list_hop = []
3436 for dic in list_hop:
3437 w2_shift = dic['R'][direction]
3438 n_periodic = tuple(dic['R'][d] for d in periodic)
3439 for w in range(width):
3440 w2 = w + w2_shift
3441 if 0 <= w2 < width:
3442 new_list_hop.append({'i': w*n_sites + dic['i'],
3443 'j': w2*n_sites + dic['j'],
3444 'R': n_periodic,
3445 't': dic['t']})
3446 rib.set_hopping(new_list_hop)
3447 if onsite is not None:
3448 rib.set_onsite(onsite)
3449 return rib
3452def magnetic_supercell(
3453 lat: Lattice,
3454 list_hop: list[dict],
3455 p: int,
3456 q: int,
3457 onsite: dict | None = None,
3458 spin: bool = False,
3459) -> KSpace:
3460 r'''
3461 Thread a uniform magnetic field through a 2D periodic model, with a
3462 rational flux :math:`\alpha = p/q` flux quanta per unit cell, and return
3463 the Bloch Hamiltonian of its magnetic unit cell: *q* cells along
3464 :math:`\mathbf{a}_1`. Each hopping gets its Peierls phase (see
3465 *System.set_peierls_phase*)
3467 .. math::
3469 \phi_{ij} = \frac{2\pi}{\Phi_0}\int_{\mathbf{r}_i}^{\mathbf{r}_j}\mathbf{A}\cdot d\mathbf{l}
3470 = 2\pi\alpha\,\frac{s_{1,i}+s_{1,j}}{2}\,(s_{2,j}-s_{2,i})
3472 in the gauge :math:`\mathbf{A} = \alpha\Phi_0\, s_1\nabla s_2`, with
3473 :math:`(s_1, s_2)` the fractional coordinates of a point along
3474 :math:`(\mathbf{a}_1, \mathbf{a}_2)`, periodic along :math:`\mathbf{a}_2`.
3475 A translation by :math:`q\mathbf{a}_1` is a symmetry only up to the
3476 gauge transformation :math:`e^{-2\pi i p s_2}`, which the hoppings
3477 crossing the magnetic cell boundary carry. The bands of the result are
3478 the Hofstadter subbands, and their *chern_number* the TKNN integers.
3480 :param lat: **Lattice** class instance, 2D (two primitive vectors in the plane).
3481 :param list_hop: List of dictionaries, the hoppings of the field-free
3482 model, in the format of *KSpace.set_hopping*.
3483 :param p: Integer. Numerator of the flux per unit cell.
3484 :param q: Positive integer. Denominator of the flux per unit cell: the
3485 magnetic cell has *q* unit cells.
3486 :param onsite: Dictionary. Default value None. Onsite energies, in the
3487 format of *KSpace.set_onsite*, applied to every copy.
3488 :param spin: Boolean. Default value False. See *KSpace*.
3490 :returns:
3491 * **mag** -- **KSpace** instance, with ``q * len(lat.unit_cell)``
3492 sites: copy *m* of site *i* is site ``m*len(lat.unit_cell) + i``,
3493 displaced by :math:`m\mathbf{a}_1`.
3495 Example usage::
3497 # square lattice, flux 1/3 per plaquette: three Hofstadter bands
3498 mag = magnetic_supercell(lattices.square(), list_hop, p=1, q=3)
3499 chern = [mag.chern_number(bands=[n]) for n in range(3)]
3500 '''
3501 error_handling.lat(lat)
3502 error_handling.planar_cell(lat.prim_vec)
3503 error_handling.set_hopping_kspace(list_hop, len(lat.unit_cell), 2, spin)
3504 error_handling.integer(p, 'p')
3505 error_handling.positive_int(q, 'q')
3506 alpha = p / q
3507 n_sites = len(lat.unit_cell)
3508 a = np.array(lat.prim_vec, dtype='f8')
3509 frac = np.linalg.solve(a.T, np.array([dic['r0'] for dic in lat.unit_cell], dtype='f8').T).T
3510 new_unit_cell = []
3511 for m in range(q):
3512 for dic in lat.unit_cell:
3513 r0 = np.array(dic['r0'], dtype='f8') + m * a[0]
3514 new_unit_cell.append({'tag': dic['tag'], 'r0': (float(r0[0]), float(r0[1]))})
3515 new_lat = Lattice(unit_cell=new_unit_cell, prim_vec=[tuple(q * a[0]), lat.prim_vec[1]])
3516 mag = KSpace(new_lat, spin=spin)
3517 new_list_hop = []
3518 for dic in list_hop:
3519 n1, n2 = dic['R']
3520 for m in range(q):
3521 s1_i, s2_i = m + frac[dic['i'], 0], frac[dic['i'], 1]
3522 s1_j, s2_j = m + n1 + frac[dic['j'], 0], n2 + frac[dic['j'], 1]
3523 big, m_j = divmod(m + n1, q)
3524 phase = 2 * PI * alpha * 0.5 * (s1_i + s1_j) * (s2_j - s2_i) - 2 * PI * p * big * s2_j
3525 new_list_hop.append({'i': m * n_sites + dic['i'], 'j': m_j * n_sites + dic['j'],
3526 'R': (big, n2), 't': dic['t'] * np.exp(1j * phase)})
3527 mag.set_hopping(new_list_hop)
3528 if onsite is not None:
3529 mag.set_onsite(onsite)
3530 return mag
3533#################################
3534# K-PATHS AND CONSTANT-ENERGY CONTOURS
3535#################################
3538def _response_weights(
3539 en: NDArray[np.float64], mu: float, temperature: float, response: str,
3540) -> NDArray[np.float64]:
3541 r'''
3542 Private. The weights :math:`P_{nm}` of the Kubo pair sum of
3543 *KSpace._kubo*, shape (nk, n, n), from the bands *en* (nk, n). With
3544 :math:`x_n = (E_n-\mu)/T`, the conductivities are
3545 :math:`\sum_n w_n\Omega_n`, i.e. :math:`P_{nm} = -(w_n - w_m)`, with
3546 :math:`w = f` ('hall'), the entropy per state
3547 :math:`s = -f\ln f-(1-f)\ln(1-f)` ('nernst') and
3548 :math:`c_2/T = T\int_x^\infty y^2(-f'(y))\,dy` ('thermal'). The orbital
3549 magnetization ('magnetization'),
3550 :math:`\sum_n[f_n\,\mathrm{Im}\langle\partial_xu_n|(H-E_n)|\partial_yu_n\rangle + g_n\Omega_n]`
3551 with :math:`g = T\ln(1+e^{-x})` (:math:`\max(\mu-E, 0)` at :math:`T = 0`),
3552 has :math:`P_{nm} = \frac12(f_n+f_m)(E_m-E_n) - (g_n-g_m)`.
3553 '''
3554 if response == 'magnetization':
3555 f = occupation.fermi_dirac(en, mu, temperature)
3556 if temperature == 0:
3557 g = np.maximum(mu - en, 0.)
3558 else:
3559 g = temperature * np.logaddexp(0., -(en - mu) / temperature)
3560 return 0.5 * (f[:, :, None] + f[:, None, :]) * (en[:, None, :] - en[:, :, None]) \
3561 - (g[:, :, None] - g[:, None, :])
3562 if response == 'hall':
3563 w = occupation.fermi_dirac(en, mu, temperature)
3564 else:
3565 a = np.abs(en - mu) / temperature
3566 if response == 'nernst':
3567 w = np.log1p(np.exp(-a)) + a * expit(-a)
3568 else:
3569 # int_a^inf y^2 (-f') dy = a^2 f(a) + 2 [a ln(1 + e^-a) - Li2(-e^-a)]
3570 # for a >= 0, and pi^2/3 minus it below the Fermi level
3571 # (Li2(z) = spence(1 - z) in scipy's convention)
3572 c2 = a ** 2 * expit(-a) + 2 * (a * np.log1p(np.exp(-a)) - spence(1 + np.exp(-a)))
3573 w = temperature * np.where(en >= mu, c2, PI ** 2 / 3 - c2)
3574 return -(w[:, :, None] - w[:, None, :])
3577def _path(
3578 points: list[ArrayLike], nk: int,
3579) -> tuple[NDArray[np.float64], NDArray[np.float64], list[float]]:
3580 '''
3581 Private function. The k-points of a path through *points* (*nk* per
3582 segment, the last point included), their cumulative distance, and the
3583 positions of the nodes along it.
3584 '''
3585 pts = np.atleast_2d(np.asarray(points, dtype='f8'))
3586 ks = np.concatenate([np.linspace(pts[i], pts[i+1], nk, endpoint=False)
3587 for i in range(len(pts) - 1)] + [pts[-1:]])
3588 dist = np.concatenate([[0.], np.cumsum(np.linalg.norm(np.diff(ks, axis=0), axis=1))])
3589 return ks, dist, dist[::nk][:len(pts)-1].tolist() + [dist[-1]]
3592def _marching_table(dim: int) -> dict[int, list[list[tuple[int, int]]]]:
3593 '''
3594 Private function. For each pattern of the vertices of a triangle
3595 (*dim* = 2) or tetrahedron (*dim* = 3) above a level (bit m set if
3596 vertex m is above), the pieces of the level set, as lists of the edges
3597 (a, b) that they cross: one segment in 2D; one triangle, or two
3598 triangles forming a quadrilateral when two vertices are above, in 3D.
3599 '''
3600 table = {}
3601 for pattern in range(1, 2 ** (dim + 1) - 1):
3602 up = [m for m in range(dim + 1) if pattern >> m & 1]
3603 down = [m for m in range(dim + 1) if not pattern >> m & 1]
3604 if dim == 3 and len(up) == 2:
3605 (a, b), (c, d) = up, down
3606 table[pattern] = [[(a, c), (a, d), (b, d)], [(a, c), (b, d), (b, c)]]
3607 else:
3608 table[pattern] = [[(a, b) for a in up for b in down]]
3609 return table
3612def _level_set(values: NDArray[np.float64], level: float) -> NDArray[np.float64]:
3613 '''
3614 Private function. The level set ``values == level`` of a function on a
3615 periodic 2D or 3D mesh (*values* of shape (n1, n2) or (n1, n2, n3)),
3616 linear in each simplex of *dos._mesh_simplices*: segments (shape
3617 (n, 2, 2)) or triangles (shape (n, 3, 3)), in fractional coordinates
3618 (not folded: between 0 and 1 + 1/n_d).
3619 '''
3620 shape, dim = values.shape, values.ndim
3621 offsets = dos._mesh_simplices(dim)
3622 vals = dos._simplex_values(values - level, offsets).reshape(dim + 1, -1)
3623 pattern = sum((vals[m] > 0).astype(int) << m for m in range(dim + 1))
3624 pieces = [np.empty((0, dim, dim))]
3625 for pat, polygons in _marching_table(dim).items():
3626 sel = np.nonzero(pattern == pat)[0]
3627 if not len(sel):
3628 continue
3629 simplex, cell = np.divmod(sel, int(np.prod(shape)))
3630 corner = np.stack(np.unravel_index(cell, shape), axis=-1)
3631 pos = corner[None] + offsets[simplex].transpose(1, 0, 2)
3632 v = vals[:, sel]
3633 for polygon in polygons:
3634 # the crossing on edge (a, b): v[a] > 0 >= v[b] or the reverse
3635 pieces.append(np.stack([pos[a] + (v[a] / (v[a] - v[b]))[:, None] * (pos[b] - pos[a])
3636 for a, b in polygon], axis=1))
3637 return np.concatenate(pieces) / np.array(shape)
3640def _zone_polygon(rec: NDArray[np.float64]) -> NDArray[np.float64]:
3641 '''
3642 Private function. Vertices, counterclockwise, of the first Brillouin
3643 zone (Wigner-Seitz cell) of the 2D reciprocal lattice with basis *rec*
3644 (rows): the points at least as close to the origin as to any other
3645 reciprocal lattice vector, where two of the bisectors meet.
3646 '''
3647 gs = np.array([n @ rec for n in product(range(-2, 3), repeat=2) if any(n)])
3648 half = np.sum(gs ** 2, axis=1) / 2
3649 tol = 1e-9 * half.max()
3650 found = []
3651 for i in range(len(gs)):
3652 for j in range(i + 1, len(gs)):
3653 mat = gs[[i, j]]
3654 if abs(np.linalg.det(mat)) > tol:
3655 k = np.linalg.solve(mat, half[[i, j]])
3656 if np.all(gs @ k <= half + tol):
3657 found.append(k)
3658 verts = np.unique(np.round(np.array(found), 9), axis=0)
3659 return verts[np.argsort(np.arctan2(verts[:, 1], verts[:, 0]))]
3662def _gauss_reduce(u: NDArray[np.float64], v: NDArray[np.float64]):
3663 '''
3664 Private function. Gauss-reduced basis of a 2D lattice: the two
3665 shortest independent vectors, |u| <= |v|, at an obtuse (or right)
3666 angle, so that the edges of the Wigner-Seitz cell are normal to
3667 u, v and u + v.
3668 '''
3669 while True:
3670 if u @ u > v @ v:
3671 u, v = v, u
3672 m = np.round(u @ v / (u @ u))
3673 if m == 0:
3674 break
3675 v = v - m * u
3676 return (u, v - u) if u @ v > 0 else (u, v)
3679def _successive_minima(vecs: NDArray[np.float64]):
3680 '''
3681 Private function. The shortest lattice vectors of a 3D lattice with
3682 basis *vecs* (rows): the three successive minima, the shells of
3683 lattice vectors sorted by length (list of arrays), and their lengths.
3684 '''
3685 pts = np.array([n @ vecs for n in product(range(-3, 4), repeat=3) if any(n)])
3686 lengths = np.linalg.norm(pts, axis=1)
3687 # by length, then (among equal lengths) the most positive first: x, then y, then z
3688 order = np.lexsort((-pts[:, 2], -pts[:, 1], -pts[:, 0], np.round(lengths / lengths.min(), 6)))
3689 pts, lengths = pts[order], lengths[order]
3690 minima: list = []
3691 for p in pts:
3692 if np.linalg.matrix_rank(np.array(minima + [p]), tol=1e-8 * lengths[0]) == len(minima) + 1:
3693 minima.append(p)
3694 if len(minima) == 3:
3695 break
3696 bounds = np.nonzero(np.diff(lengths) > 1e-6 * lengths[0])[0] + 1
3697 return np.array(minima), np.split(pts, bounds), [s[0] for s in np.split(lengths, bounds)]
3700def _orthogonal_axes(shell: NDArray[np.float64]) -> NDArray[np.float64]:
3701 '''
3702 Private function. Three mutually orthogonal unit vectors among the
3703 vectors of *shell* (the conventional cubic axes).
3704 '''
3705 x = shell[0] / np.linalg.norm(shell[0])
3706 y = next(s for s in shell if abs(s @ x) < 1e-6 * np.linalg.norm(s))
3707 y = y / np.linalg.norm(y)
3708 return np.array([x, y, np.cross(x, y)])
3711def cubic_path(
3712 points: dict[str, ArrayLike], names: list[str], a: float = 1., axes: ArrayLike | None = None,
3713) -> tuple[list[NDArray[np.float64]], list[str]]:
3714 r'''
3715 Build a path through k-points given in the conventional coordinates of
3716 a cubic crystal: in units of :math:`2\pi/a` along the cubic axes, with
3717 :math:`a` the side of the conventional cube, as tables of high-symmetry
3718 points list them (e.g. :math:`X = (0, 1, 0)` and
3719 :math:`L = (\tfrac12, \tfrac12, \tfrac12)` for fcc). *high_symmetry_path*
3720 builds the standard simple cubic, fcc and bcc paths with it; use it
3721 for another path, or other points, of the same lattices.
3723 :param points: Dictionary. key: label, val: tuple of 3 real numbers,
3724 the point in units of :math:`2\pi/a`.
3725 :param names: List of at least two labels (keys of *points*), the
3726 nodes of the path in order.
3727 :param a: Positive real. Default value 1. Side of the conventional cube.
3728 :param axes: Real array, shape (3, 3). Default value None: the
3729 Cartesian axes. Rows: the three orthonormal cubic axes, in the
3730 coordinates of *KSpace.get_ham* (for a crystal whose cube is rotated).
3732 :returns:
3733 * **points** -- List of real ndarrays, the k-points, for *KSpace.k_path*.
3734 * **labels** -- List of strings, *names*, for *node_labels*.
3736 Example usage::
3738 # fcc, along Gamma - X - U and K - Gamma - L
3739 pts = {'G': (0, 0, 0), 'X': (0, 1, 0), 'U': (.25, 1, .25),
3740 'K': (.75, .75, 0), 'L': (.5, .5, .5)}
3741 points, labels = cubic_path(pts, ['G', 'X', 'U'], a=3.6)
3742 '''
3743 error_handling.cubic_path(points, names, axes)
3744 error_handling.positive_real(a, 'a')
3745 axes = np.eye(3) if axes is None else np.asarray(axes, dtype='f8')
3746 return [2 * PI / a * np.asarray(points[n], dtype='f8') @ axes for n in names], list(names)
3749def _bravais_3d(vecs: NDArray[np.float64]):
3750 '''
3751 Private function. The points and labels of the high-symmetry path of a
3752 3D lattice (basis *vecs*, rows), or (None, None) if it is not one of
3753 the simple cubic, fcc, bcc, simple tetragonal or hexagonal lattices.
3754 '''
3755 minima, shells, radii = _successive_minima(vecs)
3756 lengths = np.linalg.norm(minima, axis=1)
3757 pairs = [(0, 1), (0, 2), (1, 2)]
3758 cosine = {p: abs(minima[p[0]] @ minima[p[1]]) / (lengths[p[0]] * lengths[p[1]]) for p in pairs}
3759 equal = {p: abs(lengths[p[0]] - lengths[p[1]]) < 1e-6 * lengths[2] for p in pairs}
3760 right = {p: cosine[p] < 1e-6 for p in pairs}
3761 G = r'$\Gamma$'
3762 if all(right.values()):
3763 if all(equal.values()): # simple cubic
3764 return cubic_path({G: (0, 0, 0), 'X': (0, .5, 0), 'M': (.5, .5, 0), 'R': (.5, .5, .5)},
3765 [G, 'X', 'M', G, 'R', 'X'], lengths[0], minima / lengths[:, None])
3766 if sum(equal.values()) != 1: # orthorhombic
3767 return None, None
3768 # simple tetragonal: a, a, c
3769 (i, j), = [p for p in pairs if equal[p]]
3770 k = 3 - i - j
3771 a, c = lengths[i], lengths[k]
3772 axes = np.array([minima[i] / a, minima[j] / a, minima[k] / c])
3773 pts = {G: (0, 0, 0), 'X': (0, .5, 0), 'M': (.5, .5, 0), 'Z': (0, 0, .5),
3774 'R': (0, .5, .5), 'A': (.5, .5, .5)}
3775 names = [G, 'X', 'M', G, 'Z', 'R', 'A', 'Z']
3776 scale = 2 * PI * np.array([1 / a, 1 / a, 1 / c])
3777 return [(scale * pts[n]) @ axes for n in names], names
3778 for (i, j) in pairs:
3779 k = 3 - i - j
3780 if equal[(i, j)] and abs(cosine[(i, j)] - .5) < 1e-6 and \
3781 all(right[p] for p in pairs if p != (i, j)):
3782 # hexagonal: reciprocal vectors of (a1, a2, c), in-plane ones at 120 degrees
3783 rec = 2 * PI * np.linalg.inv(minima[[i, j, k]]).T
3784 b1, b2 = _gauss_reduce(rec[0], rec[1])
3785 M, K, A = b1 / 2, (2 * b1 + b2) / 3, rec[2] / 2
3786 hex_pts = {G: np.zeros(3), 'M': M, 'K': K, 'A': A, 'L': M + A, 'H': K + A}
3787 names = [G, 'M', 'K', G, 'A', 'L', 'H', 'A']
3788 return [hex_pts[n] for n in names], names
3789 if not all(equal.values()):
3790 return None, None
3791 if len(shells[0]) == 12: # fcc: 12 nearest neighbours at a/sqrt(2), then 6 at a
3792 return cubic_path({G: (0, 0, 0), 'X': (0, 1, 0), 'W': (.5, 1, 0), 'K': (.75, .75, 0),
3793 'L': (.5, .5, .5), 'U': (.25, 1, .25)},
3794 [G, 'X', 'W', 'K', G, 'L', 'U', 'W', 'L', 'K'],
3795 radii[1], _orthogonal_axes(shells[1]))
3796 if len(shells[0]) == 8 and all(abs(c - 1 / 3) < 1e-6 for c in cosine.values()):
3797 # bcc: 8 nearest neighbours at a sqrt(3)/2, then 6 at a
3798 return cubic_path({G: (0, 0, 0), 'H': (0, 0, 1), 'N': (.5, 0, .5), 'P': (.5, .5, .5)},
3799 [G, 'H', 'N', G, 'P', 'H'], radii[1], _orthogonal_axes(shells[1]))
3800 return None, None
3803def high_symmetry_path(lat: Lattice) -> tuple[list[NDArray[np.float64]], list[str]]:
3804 r'''
3805 Get a standard path through the high-symmetry points of the Brillouin
3806 zone of a lattice, for *KSpace.k_path* and the labels of
3807 *KSpace.plot_bands*. The Bravais lattice is recognized from
3808 *lat.prim_vec*, whatever the choice of primitive vectors and the
3809 orientation of the lattice:
3811 * 1D: :math:`-X, \Gamma, X` (the whole zone, :math:`X = \mathbf{b}/2`).
3812 * 2D (the five Bravais lattices): square :math:`\Gamma X M \Gamma`;
3813 rectangular :math:`\Gamma X S Y \Gamma`; hexagonal
3814 :math:`\Gamma M K \Gamma`; centred rectangular and oblique
3815 :math:`\Gamma X H C H_1 Y \Gamma`, along the edges of the
3816 zone: :math:`X, Y, C` are the midpoints of the edges normal to the
3817 reciprocal vectors :math:`\mathbf{b}_1, \mathbf{b}_2,
3818 \mathbf{b}_1+\mathbf{b}_2` of the reduced basis (with
3819 :math:`|\mathbf{b}_1| \le |\mathbf{b}_2|`), :math:`H` and
3820 :math:`H_1` the corners between them.
3821 * 3D: simple cubic :math:`\Gamma X M \Gamma R X`; fcc
3822 :math:`\Gamma X W K \Gamma L U W L K`; bcc :math:`\Gamma H N \Gamma P H`;
3823 simple tetragonal :math:`\Gamma X M \Gamma Z R A Z`; hexagonal
3824 :math:`\Gamma M K \Gamma A L H A`. These are the paths of Setyawan
3825 and Curtarolo (Comput. Mater. Sci. 49, 299 (2010)), up to their
3826 first discontinuity (a path here is one continuous line).
3828 The other 3D lattices (orthorhombic, monoclinic, body-centred
3829 tetragonal, rhombohedral, ...) raise a ValueError: give their points
3830 to *k_path* directly.
3832 :param lat: **lattice** class instance (only *prim_vec* is used).
3834 :returns:
3835 * **points** -- List of real ndarrays, the k-points, in the
3836 coordinates of *KSpace.get_ham*.
3837 * **labels** -- List of strings (LaTeX), for *node_labels*.
3839 Example usage::
3841 points, labels = high_symmetry_path(lat)
3842 ks.k_path(points, nk=60)
3843 ks.plot_bands(node_labels=labels)
3844 '''
3845 ks = KSpace(lat)
3846 vecs = np.array(lat.prim_vec, dtype='f8') @ ks.k_basis # in the k coordinates
3847 G = r'$\Gamma$'
3848 if ks.dim == 1:
3849 x = PI / vecs[0]
3850 return [-x, np.zeros(1), x], ['$-X$', G, 'X']
3851 if ks.dim == 3:
3852 points, labels = _bravais_3d(vecs)
3853 error_handling.bravais_lattice(labels, 3)
3854 return points, labels
3855 rec = 2 * PI * np.linalg.inv(vecs).T
3856 b1, b2 = _gauss_reduce(rec[0], rec[1])
3857 l1, l2, l3 = (np.linalg.norm(b) for b in (b1, b2, b1 + b2))
3858 close = lambda a, b: abs(a - b) < 1e-6 * max(abs(a), abs(b))
3859 cosine = b1 @ b2 / (l1 * l2)
3860 if abs(cosine) < 1e-6:
3861 if close(l1, l2):
3862 return [np.zeros(2), b1 / 2, (b1 + b2) / 2, np.zeros(2)], [G, 'X', 'M', G]
3863 return [np.zeros(2), b1 / 2, (b1 + b2) / 2, b2 / 2, np.zeros(2)], [G, 'X', 'S', 'Y', G]
3864 if close(l1, l2) and close(cosine, -.5):
3865 return [np.zeros(2), b1 / 2, (2 * b1 + b2) / 3, np.zeros(2)], [G, 'M', 'K', G]
3866 # centred rectangular and oblique: around half of the hexagonal zone
3867 b3 = b1 + b2
3868 corner = lambda u, v: np.linalg.solve(np.array([u, v]), np.array([u @ u, v @ v]) / 2)
3869 return ([np.zeros(2), b1 / 2, corner(b1, b3), b3 / 2, corner(b2, b3), b2 / 2, np.zeros(2)],
3870 [G, 'X', 'H', 'C', '$H_1$', 'Y', G])