Coverage for tbkit/kspace.py: 100%
1465 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +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}
39# The symmetry kinds of KSpace.non_hermitian_class, in its argument order:
40# (complex conjugation, transposition, minus sign) of U op(H(k)) U^dagger =
41# +-H(+-k). TRS, PHS, CS, TRS^dagger, PHS^dagger, SLS, pseudo-Hermiticity.
42_NH_KINDS = ((1, 0, 0), (0, 1, 1), (1, 1, 1), (0, 1, 0), (1, 0, 1), (0, 0, 1), (1, 1, 0))
44# Names of KSpace.non_hermitian_class equivalent under H -> iH (Kawabata et
45# al. 2019, Tables XII and XIV), mapped to the representative of the 38 classes.
46_NH_CANONICAL = {
47 'D†': 'AI', 'C†': 'AII', 'AII + S-': 'AI + S-', 'DIII + S-+': 'BDI + S-+',
48 'DIII + S--': 'BDI + S--', 'CI + S-+': 'CII + S-+', 'CI + S--': 'CII + S--',
49 'A + η': 'AIII', 'AI + η+': 'BDI†', 'AI + η-': 'DIII†', 'AII + η+': 'CII†',
50 'AII + η-': 'CI†', 'D + η+': 'BDI', 'D + η-': 'DIII', 'C + η+': 'CII', 'C + η-': 'CI',
51}
54#################################
55# CLASS KSPACE
56#################################
59def reciprocal_vectors(prim_vec: Sequence[tuple[float, ...]]) -> list[tuple[float, ...]]:
60 r'''
61 Get the reciprocal lattice vectors :math:`\mathbf{b}_i` such that
62 :math:`\mathbf{a}_i\cdot\mathbf{b}_j = 2\pi\delta_{ij}`, lying in the
63 span of the :math:`\mathbf{a}_i`.
65 Parameters
66 ----------
67 prim_vec : List of one/two/three tuples
68 Primitive vectors (see class **lattice**).
70 Returns
71 -------
72 rec_vec : List of one/two/three tuples
73 Reciprocal vectors, with
74 as many components as the primitive vectors.
75 '''
76 if len(prim_vec[0]) == 3:
77 # b = 2 pi (A^+)^T: the dual basis within the span of the a_i
78 rec = 2 * PI * np.linalg.pinv(np.array(prim_vec, dtype='f8')).T
79 return [tuple(float(c) for c in b) for b in rec]
80 if len(prim_vec) == 1:
81 ax, ay = prim_vec[0]
82 norm2 = ax ** 2 + ay ** 2
83 return [(2*PI*ax/norm2, 2*PI*ay/norm2)]
84 (a1x, a1y), (a2x, a2y) = prim_vec
85 area = a1x * a2y - a1y * a2x
86 b1 = (2*PI*a2y/area, -2*PI*a2x/area)
87 b2 = (-2*PI*a1y/area, 2*PI*a1x/area)
88 return [b1, b2]
91class KSpace():
92 r'''
93 Build and solve the Tight-Binding Bloch Hamiltonian :math:`H(\mathbf{k})`
94 of a periodic lattice defined by the class **lattice**.
96 Hoppings are defined between orbitals of the unit cell, separated by a
97 lattice vector :math:`\mathbf{R} = n_1\mathbf{a}_1+n_2\mathbf{a}_2`:
99 .. math::
101 H_{ij}(\mathbf{k}) = \sum_{\mathbf{R}} t_{ij}(\mathbf{R})\,
102 e^{i\mathbf{k}\cdot\mathbf{R}}
104 Parameters
105 ----------
106 lat : **lattice** class instance
107 Only *unit_cell* and *prim_vec*
108 are used (the instance need not call *get_lattice*).
109 spin : Boolean, default False
110 If True, every site of
111 *unit_cell* carries a spin-1/2 degree of freedom (*norb* doubles to
112 ``2*len(unit_cell)``, ordered site-major: orbitals ``2*i, 2*i+1``
113 are the up/down components of site *i*). *set_onsite* and
114 *set_hopping* then accept 2x2 (spin) matrices in addition to plain
115 numbers, to build spin-orbit coupling or Zeeman terms -- see
116 :data:`PAULI` for ready-made Pauli matrices.
118 Examples
119 --------
120 Example usage::
122 # graphene, nearest-neighbor hopping t
123 DX, DY = 0.5 * 3 ** 0.5, 0.5
124 unit_cell = [{'tag': 'a', 'r0': (0., 0.)}, {'tag': 'b', 'r0': (DX, DY)}]
125 prim_vec = [(2*DX, 0.), (DX, 1.5)]
126 lat = Lattice(unit_cell=unit_cell, prim_vec=prim_vec)
127 gra = KSpace(lat)
128 gra.set_hopping([{'i': 0, 'j': 1, 'R': (0, 0), 't': 1.},
129 {'i': 0, 'j': 1, 'R': (-1, 0), 't': 1.},
130 {'i': 0, 'j': 1, 'R': (0, -1), 't': 1.}])
131 '''
133 workers = 1 # threads diagonalizing over k, see set_workers
134 # _bloch_sum groups the hoppings by lattice vector when they fill at
135 # least this fraction of the dense (n_R, norb, norb) blocks; sparser
136 # models (large supercells) are summed hopping by hopping
137 group_fill = 0.01
139 def __init__(self, lat: Lattice, spin: bool = False) -> None:
140 error_handling.lat(lat)
141 error_handling.boolean(spin, 'spin')
142 error_handling.independent(lat.prim_vec)
143 self.lat = lat
144 self.dim = len(lat.prim_vec)
145 self.space_dim = len(lat.prim_vec[0])
146 self.spin = spin
147 self.n_sites = len(lat.unit_cell)
148 self.norb = 2*self.n_sites if spin else self.n_sites
149 self.tags = np.array([dic['tag'] for dic in lat.unit_cell])
150 self.onsite = np.zeros(self.norb, 'c16')
151 # spin-off-diagonal onsite terms (in-plane Zeeman, onsite Rashba):
152 # they have no place on the diagonal `onsite` array.
153 self._onsite_offdiag = np.zeros((self.norb, self.norb), 'c16')
154 self._hop_const: list = [] # list of (i, j, R_cartesian (np.ndarray), t)
155 # value functions (see tbkit.values): (i, j, R_cartesian, function,
156 # hermitian), with arrays of sites i, j and bond vectors R
157 self._hop_values: list = []
158 #: Parameters of the value functions (see *set_params*).
159 self.params: dict = {}
160 self._nonreciprocal = False # set_hopping(hermitian=False) was used
161 self._overlap_hop: list = [] # overlaps, as _hop: list of (i, j, R_cartesian, s)
162 self.rec_vec = reciprocal_vectors(lat.prim_vec)
163 # Orthonormal basis (columns) of the span of prim_vec, in which k is
164 # given: the identity when the lattice fills its space, a1/|a1| for a
165 # chain, the Gram-Schmidt basis of (a1, a2) for a 2D lattice in 3D.
166 self.k_basis = np.linalg.qr(np.array(lat.prim_vec, dtype='f8').T)[0]
167 if self.dim == self.space_dim:
168 self.k_basis = np.eye(self.space_dim)
169 elif self.dim == 1:
170 a1 = np.asarray(lat.prim_vec[0], dtype='f8')
171 self.k_basis = (a1 / np.linalg.norm(a1))[:, None]
172 else:
173 # keep the first axis along a1 (QR may flip signs)
174 signs = np.sign(np.diag(np.array(lat.prim_vec, dtype='f8') @ self.k_basis))
175 self.k_basis = self.k_basis * np.where(signs == 0, 1., signs)[None, :]
176 #: Reciprocal vectors in the k coordinates used by *get_ham* (rows).
177 self.rec_vec_k = np.array(self.rec_vec, dtype='f8') @ self.k_basis
178 self.ks = np.array([]) # k-points of the last band-structure calculation
179 self.ks_dist = np.array([]) # cumulative distance along the k-path
180 self.nodes = np.array([]) # positions, along ks_dist, of the k-path nodes
181 self.en = np.array([]) # bands, shape (len(ks), norb)
183 @property
184 def _hop(self) -> list:
185 '''
186 Private. The hoppings, as a list of (i, j, R_cartesian, t): those
187 given as numbers, then the value functions evaluated with *params*.
188 Assigning it replaces both.
189 '''
190 return self._hops(self.params)
192 @_hop.setter
193 def _hop(self, hops: list) -> None:
194 self._hop_const = list(hops)
195 self._hop_values = []
197 def _hops(self, params: dict) -> list:
198 '''
199 Private method. *_hop*, with the value functions evaluated with *params*.
200 '''
201 hops = list(self._hop_const)
202 if not self._hop_values:
203 return hops
204 tau = np.array([dic['r0'] for dic in self.lat.unit_cell], dtype='f8').reshape(self.n_sites, -1)
205 for i, j, R, func, hermitian in self._hop_values:
206 site_i = values.sites(tau[i], self.tags[i], i)
207 site_j = values.sites(tau[j] + R, self.tags[j], j)
208 t = values.evaluate(func, (site_i, site_j), params, len(i), self.spin)
209 if self.spin:
210 pairs = [(2*ii + a, 2*jj + b, RR, tt[a, b], np.conj(tt[a, b]))
211 for ii, jj, RR, tt in zip(i, j, R, t) for a in range(2) for b in range(2)]
212 else:
213 pairs = [(ii, jj, RR, tt, np.conj(tt)) for ii, jj, RR, tt in zip(i, j, R, t)]
214 for ii, jj, RR, tt, tc in pairs:
215 hops.append((int(ii), int(jj), RR, tt))
216 if hermitian:
217 hops.append((int(jj), int(ii), -RR, tc))
218 return hops
220 def set_params(self, **params) -> None:
221 '''
222 Set the parameters of the value functions of *set_hopping*: every
223 method (*get_bands*, *chern_number*, *finite_ham*, ...) then uses
224 these values. *get_ham* also takes them per call.
226 Parameters
227 ----------
228 params
229 Parameter values, added to (or replacing those of) *params*.
231 Examples
232 --------
233 Example usage::
235 ks.set_params(m=0.5)
236 ks.chern_number()
237 '''
238 self.params.update(params)
240 def set_onsite(self, dict_onsite: dict[str, complex | Sequence[complex]]) -> None:
241 '''
242 Set the onsite energies, by sublattice tag.
244 Parameters
245 ----------
246 dict_onsite : Dictionary
247 key: tag, val: onsite energy
248 (a plain number), or, if ``spin=True``, a plain number (applied
249 equally to both spins), a pair ``(E_up, E_down)`` of numbers (a
250 spin splitting along z), or a 2x2 complex matrix (a general spin
251 structure, e.g. an in-plane Zeeman field built from :data:`PAULI`).
252 A complex onsite energy (gain/loss), or a non-Hermitian 2x2
253 block, makes :math:`H(\\mathbf{k})` non-Hermitian: *get_bands*
254 then returns complex energies, sorted by their real part.
256 Examples
257 --------
258 Example usage::
260 kag.set_onsite({'a': 1., 'b': -1.})
261 # spinful: same onsite energy for both spins on 'a', a Zeeman
262 # splitting along z on 'b':
263 kag_spin.set_onsite({'a': 1., 'b': (1., -1.)})
264 # spinful: an in-plane Zeeman field on 'a':
265 kag_spin.set_onsite({'a': Bx*PAULI['x']})
266 '''
267 error_handling.set_onsite_kspace(dict_onsite, self.lat.tags, self.spin)
268 for tag, val in dict_onsite.items():
269 sites = np.where(self.tags == tag)[0]
270 if self.spin:
271 if np.ndim(val) == 0:
272 block = val * PAULI['0']
273 elif np.ndim(val) == 2:
274 block = np.asarray(val, 'c16')
275 else:
276 block = np.diag(np.asarray(val, 'c16'))
277 self.onsite[2*sites] = block[0, 0]
278 self.onsite[2*sites + 1] = block[1, 1]
279 for site in sites:
280 self._onsite_offdiag[2*site, 2*site + 1] = block[0, 1]
281 self._onsite_offdiag[2*site + 1, 2*site] = block[1, 0]
282 else:
283 self.onsite[sites] = val
285 def set_hopping(self, list_hop: list[dict], hermitian: bool = True) -> None:
286 r'''
287 Set the hoppings between orbitals of the unit cell.
289 Only one representative of each hopping needs to be given: its
290 Hermitian conjugate (:math:`j\to i`, :math:`\mathbf{R}\to-\mathbf{R}`)
291 is added automatically (unless ``hermitian=False``).
293 Parameters
294 ----------
295 list_hop
296 List of dictionaries with keys ('i', 'j', 'R', 't'):
298 * 'i', 'j': Positive integers. Site indices within the unit cell
299 (following the order of *unit_cell*).
300 * 'R': Tuple of one/two integers :math:`(n_1, n_2)`. Lattice vector
301 :math:`\mathbf{R}=n_1\mathbf{a}_1+n_2\mathbf{a}_2` separating the
302 two sites.
303 * 't': Complex number, or, if ``spin=True``, either a complex
304 number (spin-independent hopping) or a 2x2 complex matrix (a
305 general, possibly spin-mixing, hopping -- e.g. built from
306 :data:`PAULI` for Rashba or intrinsic spin-orbit coupling).
307 Or a value function ``t(site_i, site_j, **params)``,
308 evaluated by ``get_ham(k, **params)`` (and by every other
309 method, with the parameters of *set_params*): *site_i* and
310 *site_j* are structured arrays of the bonds' end points,
311 *site_i* in the home cell and *site_j* in cell
312 :math:`\mathbf{R}` (fields 'x', 'y', ('z'), 'tag', 'index' --
313 the site index within *unit_cell*). It is called once with
314 every bond of the dictionaries that share it, and returns one
315 value, or one per bond (or, if ``spin=True``, a 2x2 matrix,
316 or one per bond).
318 *list_hop* may instead use the neighbour-order form of
319 *System.set_hopping*: dictionaries with keys ('n', 't') and
320 optionally 'ang' and/or 'tag' -- 'n' the neighbour order (1st,
321 2nd, ... shortest distance between orbitals of the infinite
322 lattice), 'ang' a bond angle in degrees, 'tag' a sublattice pair
323 such as ``'ab'``. The bonds are found from *unit_cell* and
324 *prim_vec* (in 1D, 2D and 3D) and turned into the explicit form
325 by *tbkit.neighbours.neighbour_hoppings* (see there for the bond
326 orientations; a negative angle addresses the reversed bonds).
327 A list uses one form or the other, not both.
328 hermitian : Boolean, default True
329 If False, the Hermitian
330 conjugates are *not* added: each dictionary is one matrix element
331 :math:`H_{ij}(\mathbf{R})` only, so that non-reciprocal
332 (non-Hermitian) hoppings can be built -- give the reverse hopping
333 explicitly, with its own amplitude.
335 Examples
336 --------
337 Example usage::
339 # graphene, by neighbour order: nearest (t) and next-nearest
340 # (t2) neighbours, as in System.set_hopping
341 gra.set_hopping([{'n': 1, 't': t}, {'n': 2, 't': t2}])
342 # square lattice, anisotropic nearest neighbours
343 sq.set_hopping([{'n': 1, 'ang': 0., 't': tx}, {'n': 1, 'ang': 90., 't': ty}])
345 # 1D chain, nearest-neighbor hopping t between the only orbital
346 # and its right neighbor:
347 chain.set_hopping([{'i': 0, 'j': 0, 'R': (1,), 't': 1.}])
348 # spinful: spin-independent hopping t, plus a Rashba-like
349 # spin-flip term of strength alpha:
350 chain_spin.set_hopping([{'i': 0, 'j': 0, 'R': (1,),
351 't': t*PAULI['0'] + 1j*alpha*PAULI['y']}])
352 # value function: a staggered hopping, swept without rebuilding
353 chain.set_hopping([{'n': 1, 't': lambda si, sj, t, d: t + d * (-1) ** si['index']}])
354 ham = chain.get_ham((0.,), t=1., d=0.2)
355 '''
356 error_handling.boolean(hermitian, 'hermitian')
357 if error_handling.hopping_form(list_hop):
358 list_hop = neighbours.neighbour_hoppings(self.lat, list_hop, self.spin, hermitian)
359 error_handling.set_hopping_kspace(list_hop, self.n_sites, self.dim, self.spin, values=True)
360 if not hermitian:
361 self._nonreciprocal = True
362 funcs: dict = {} # the bonds of each value function: id -> (function, [(i, j, R_cart)])
363 for dic in list_hop:
364 R_cart = np.zeros(self.space_dim)
365 for n, vec in zip(dic['R'], self.lat.prim_vec):
366 R_cart += n * np.array(vec)
367 i, j, t = dic['i'], dic['j'], dic['t']
368 if callable(t):
369 funcs.setdefault(id(t), (t, []))[1].append((i, j, R_cart))
370 elif self.spin:
371 block = t*PAULI['0'] if np.ndim(t) == 0 else np.asarray(t, 'c16')
372 for a in range(2):
373 for b in range(2):
374 self._hop_const.append((2*i+a, 2*j+b, R_cart, block[a, b]))
375 if hermitian:
376 self._hop_const.append((2*j+b, 2*i+a, -R_cart, np.conj(block[a, b])))
377 else:
378 self._hop_const.append((i, j, R_cart, t))
379 if hermitian:
380 self._hop_const.append((j, i, -R_cart, np.conj(t)))
381 for func, bonds in funcs.values():
382 i, j, R = zip(*bonds)
383 self._hop_values.append((np.array(i), np.array(j), np.array(R), func, hermitian))
385 def set_overlap(self, list_hop: list[dict]) -> None:
386 r'''
387 Set the overlaps :math:`s_{ij}(\mathbf{R}) = \langle i,\mathbf{0}|j,\mathbf{R}\rangle`
388 of a non-orthogonal basis, in the format of *set_hopping* (key 't'
389 for the overlap; the Hermitian conjugates are added):
391 .. math::
393 S(\mathbf{k}) = \mathbb{1} + \sum_{\mathbf{R}} s(\mathbf{R})\,e^{i\mathbf{k}\cdot\mathbf{R}}\, ,
395 after which the bands solve :math:`H(\mathbf{k})v = E\,S(\mathbf{k})v`.
396 (With graphene's nearest-neighbour overlap :math:`s`, the
397 :math:`\pi` bands become :math:`E = (\epsilon \pm t|f|)/(1 \pm s|f|)`:
398 no longer symmetric about :math:`\epsilon`.) The topological tools use
399 the Lowdin-orthonormalized states :math:`S^{1/2}v`.
401 Parameters
402 ----------
403 list_hop
404 List of dictionaries ('i', 'j', 'R', 't'), or by
405 neighbour order ('n', 't', optionally 'ang', 'tag'), see
406 *set_hopping* (numbers only, no value functions).
407 '''
408 if error_handling.hopping_form(list_hop):
409 list_hop = neighbours.neighbour_hoppings(self.lat, list_hop, self.spin)
410 error_handling.set_hopping_kspace(list_hop, self.n_sites, self.dim, self.spin)
411 saved = self._hop_const, self._hop_values, self._nonreciprocal
412 self._hop = []
413 self.set_hopping(list_hop)
414 self._overlap_hop += self._hop_const
415 self._hop_const, self._hop_values, self._nonreciprocal = saved
417 def get_overlap(self, k: ArrayLike) -> NDArray[np.complex128]:
418 r'''
419 Get the overlap matrix :math:`S(\mathbf{k})` (the identity for an
420 orthogonal basis), see *set_overlap*.
422 Parameters
423 ----------
424 k : k point (see *get_ham*)
426 Returns
427 -------
428 S : Complex ndarray, shape (norb, norb)
429 '''
430 error_handling.k_vector(k, 'k', self.dim)
431 k_cart = self.k_basis @ np.asarray(k, dtype='f8')
432 return np.eye(self.norb) + self._bloch_sum(self._overlap_hop, k_cart[None])[0][0]
434 def clear_hopping(self) -> None:
435 '''
436 Clear the hoppings set by *set_hopping* (and the overlaps of *set_overlap*).
437 '''
438 self._hop = []
439 self._nonreciprocal = False
440 self._overlap_hop = []
442 def get_ham(self, k: ArrayLike, **params) -> NDArray[np.complex128]:
443 r'''
444 Get the dense Bloch Hamiltonian :math:`H(\mathbf{k})`.
446 Parameters
447 ----------
448 k : Tuple/list/ndarray of *dim* real numbers
449 When the lattice
450 fills its space (2D in the plane, 3D in space), the
451 :math:`\mathbf{k}` point in the same Cartesian frame as *prim_vec*.
452 In 1D, the crystal momentum *along* the primitive vector (so the
453 Brillouin zone spans :math:`2\pi/|\mathbf{a}_1|`), which
454 coincides with :math:`k_x` for a chain aligned with :math:`x`.
455 For a 2D lattice in 3D space (a slab), the components along the
456 orthonormal basis *k_basis* of its plane (the first along
457 :math:`\mathbf{a}_1`). *rec_vec_k* holds the reciprocal vectors in
458 these coordinates.
459 params
460 Values of the parameters of the value functions of
461 *set_hopping*, for this call only: they override those of
462 *set_params*. Each function gets the ones its signature names
463 (all of them with ``**kwargs``); the others are ignored.
465 Returns
466 -------
467 ham : Complex ndarray, shape (norb, norb)
469 Examples
470 --------
471 Example usage::
473 hams = [ks.get_ham((0., 0.), m=m) for m in np.linspace(-1., 1., 21)]
474 '''
475 error_handling.k_vector(k, 'k', self.dim)
476 return self._bloch_ham((self.k_basis @ np.asarray(k, dtype='f8'))[None],
477 params={**self.params, **params})[0]
479 def get_ham_peierls(self, k: ArrayLike, A: ArrayLike) -> NDArray[np.complex128]:
480 r'''
481 Get the Bloch Hamiltonian in a uniform vector potential
482 :math:`\mathbf{A}` (:math:`\hbar = e = 1`): every hopping picks up the
483 Peierls phase :math:`e^{i\mathbf{A}\cdot\mathbf{d}}` of its bond
484 vector :math:`\mathbf{d} = \mathbf{R} + \boldsymbol\tau_j - \boldsymbol\tau_i`
485 (see *System.set_peierls_phase*). A time-dependent :math:`\mathbf{A}(t)`
486 is a uniform electric field :math:`-\partial_t\mathbf{A}`, e.g. light
487 (see *tbkit.floquet*).
489 Parameters
490 ----------
491 k : k point (see *get_ham*)
492 A : Tuple/ndarray of space_dim real numbers
493 Vector potential.
495 Returns
496 -------
497 ham : Complex ndarray, shape (norb, norb)
498 '''
499 error_handling.k_vector(k, 'k', self.dim)
500 A = np.asarray(A, dtype='f8')
501 error_handling.ndarray(A, 'A', self.space_dim)
502 return self._bloch_ham((self.k_basis @ np.asarray(k, dtype='f8'))[None], A)[0]
504 def _bloch_ham(
505 self, k_cart: NDArray[np.float64], A: NDArray[np.float64] | None = None, params: dict | None = None,
506 ) -> NDArray[np.complex128]:
507 '''
508 Private method. *get_ham* (or, with a vector potential *A*,
509 *get_ham_peierls*) at many k-points at once: *k_cart* has shape
510 (nk, space_dim) and the result (nk, norb, norb). The value functions
511 take *params* (default: *self.params*).
512 '''
513 hops = self._hops(self.params if params is None else params)
514 return self._bloch_sum(hops, k_cart, A=A)[0] + (np.diag(self.onsite) + self._onsite_offdiag)[None]
516 def get_bands(
517 self, ks: ArrayLike, eigenvec: bool = False,
518 ) -> NDArray[np.float64] | tuple[NDArray[np.float64], NDArray[np.complex128]]:
519 r'''
520 Diagonalize :math:`H(\mathbf{k})` over a set of k-points, and keep
521 the result for *plot_bands* (plotted against the cumulative distance
522 between consecutive k-points).
524 Parameters
525 ----------
526 ks : ndarray, shape (nk, dim)
527 k-points.
528 eigenvec : Boolean, default False
529 If True, also return
530 the eigenvectors.
532 Returns
533 -------
534 en : Real ndarray, shape (nk, norb)
535 Band energies, sorted ascending.
536 Complex, and sorted by real part, if :math:`H(\mathbf{k})` is
537 non-Hermitian (complex onsite energies, see *set_onsite*).
538 vn : Complex ndarray, shape (nk, norb, norb), only if *eigenvec* is True
539 vn[k, :, n] is the nth (right) eigenvector at ks[k].
540 '''
541 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
542 error_handling.ks(ks, self.dim)
543 out = self._diagonalize(ks, eigenvec)
544 self.ks = ks
545 self.en = out[0] if isinstance(out, tuple) else out
546 steps = np.linalg.norm(np.diff(ks, axis=0), axis=1)
547 self.ks_dist = np.concatenate([[0.], np.cumsum(steps)])[:len(ks)]
548 self.nodes = np.array([])
549 return out
551 @overload
552 def _diagonalize(
553 self, ks: NDArray[np.float64], eigenvec: Literal[False] = ...,
554 ) -> NDArray[np.float64]: ...
555 @overload
556 def _diagonalize(
557 self, ks: NDArray[np.float64], eigenvec: Literal[True],
558 ) -> tuple[NDArray[np.float64], NDArray[np.complex128]]: ...
559 @overload
560 def _diagonalize(
561 self, ks: NDArray[np.float64], eigenvec: bool,
562 ) -> NDArray[np.float64] | tuple[NDArray[np.float64], NDArray[np.complex128]]: ...
564 def _diagonalize(
565 self, ks: NDArray[np.float64], eigenvec: bool = False,
566 ) -> NDArray[np.float64] | tuple[NDArray[np.float64], NDArray[np.complex128]]:
567 r'''
568 Private method. Diagonalize :math:`H(\mathbf{k})` over the k-points
569 *ks* (shape (nk, dim)), without touching the stored band structure.
570 '''
571 return self._eigs(ks, eigenvec)
573 def set_workers(self, workers: int) -> None:
574 r'''
575 Set the number of threads that diagonalize :math:`H(\mathbf{k})`
576 over many k-points (*get_bands*, *mesh_bands*, *berry_curvature*,
577 the Wilson loops, ...). The k-points are split into chunks, each
578 diagonalized by one stacked LAPACK call, which releases the GIL, so
579 the chunks run in parallel. *surface_spectral_function* spreads its
580 k-points (and, with fewer k-points than threads, chunks of its
581 energies) over the threads in the same way. The results depend on
582 *workers* only up to rounding (and the phases of the eigenvectors).
584 Parameters
585 ----------
586 workers : Positive integer
587 Number of threads (default 1, no threads).
589 Examples
590 --------
591 Example usage::
593 ks.set_workers(4)
594 curv = ks.berry_curvature(0, nk=400)
595 '''
596 error_handling.positive_int(workers, 'workers')
597 self.workers = workers
599 def _hams(self, ks: NDArray[np.float64]) -> NDArray[np.complex128]:
600 '''
601 Private method. *get_ham* at the k-points *ks* (shape (nk, dim)),
602 shape (nk, norb, norb): one Bloch sum, unless a subclass overrides
603 *get_ham*, which is then called per k-point (a Floquet model
604 overrides *_hams* too, see *floquet.DrivenKSpace._hams*).
605 '''
606 if type(self).get_ham is KSpace.get_ham:
607 return self._bloch_ham(ks @ self.k_basis.T)
608 return np.array([self.get_ham(k) for k in ks], dtype='c16')
610 def _overlaps(self, ks: NDArray[np.float64]) -> NDArray[np.complex128]:
611 '''
612 Private method. *get_overlap* at the k-points *ks*, shape (nk, norb, norb).
613 '''
614 return np.eye(self.norb) + self._bloch_sum(self._overlap_hop, ks @ self.k_basis.T)[0]
616 @overload
617 def _eigs(self, ks: NDArray[np.float64], eigenvec: Literal[False] = ...) -> NDArray: ...
618 @overload
619 def _eigs(
620 self, ks: NDArray[np.float64], eigenvec: Literal[True],
621 ) -> tuple[NDArray, NDArray[np.complex128]]: ...
622 @overload
623 def _eigs(
624 self, ks: NDArray[np.float64], eigenvec: bool,
625 ) -> NDArray | tuple[NDArray, NDArray[np.complex128]]: ...
627 def _eigs(
628 self, ks: NDArray[np.float64], eigenvec: bool = False,
629 ) -> NDArray | tuple[NDArray, NDArray[np.complex128]]:
630 '''
631 Private method. *_eig* at the k-points *ks* (shape (nk, dim)): the
632 eigenvalues, shape (nk, norb), and the eigenvectors, shape
633 (nk, norb, norb). The k-points go in chunks of bounded memory,
634 spread over *workers* threads (see *set_workers*).
635 '''
636 chunk = max(1, min(400000 // self.norb ** 2, -(-len(ks) // self.workers)))
637 parts = [ks[c0:c0 + chunk] for c0 in range(0, len(ks), chunk)] or [ks]
638 if self.workers > 1 and len(parts) > 1:
639 with ThreadPoolExecutor(self.workers) as pool:
640 out = list(pool.map(lambda part: self._eig_chunk(part, eigenvec), parts))
641 else:
642 out = [self._eig_chunk(part, eigenvec) for part in parts]
643 if eigenvec:
644 return np.concatenate([o[0] for o in out]), np.concatenate([o[1] for o in out])
645 return np.concatenate(out)
647 def _eig_chunk(
648 self, ks: NDArray[np.float64], eigenvec: bool,
649 ) -> NDArray | tuple[NDArray, NDArray[np.complex128]]:
650 r'''
651 Private method. One chunk of *_eigs*: the Hermitian solver when
652 *is_hermitian*, otherwise the general one, with the eigenvalues
653 sorted by real part, both stacked over the k-points. With an overlap
654 (see *set_overlap*), the generalized problem :math:`Hv = ES(\mathbf{k})v`:
655 the Hermitian one is reduced by the Cholesky factor
656 :math:`S = LL^\dagger` to :math:`L^{-1}HL^{-\dagger}` (eigenvectors
657 :math:`L^{-\dagger}v`, normalized to :math:`v^\dagger Sv = 1` as by
658 *scipy.linalg.eigh*); the non-Hermitian one has no stacked solver,
659 and is solved one k-point at a time.
660 '''
661 ham = self._hams(ks)
662 hermitian = self.is_hermitian()
663 en: NDArray
664 vn: NDArray
665 if self._overlap_hop and hermitian:
666 l_inv = np.linalg.inv(np.linalg.cholesky(self._overlaps(ks)))
667 l_inv_h = l_inv.conj().transpose(0, 2, 1)
668 ham = l_inv @ ham @ l_inv_h
669 if hermitian:
670 if not eigenvec:
671 return np.linalg.eigvalsh(ham)
672 en, vn = np.linalg.eigh(ham)
673 return en, (l_inv_h @ vn if self._overlap_hop else vn)
674 if self._overlap_hop:
675 solved = [LA.eig(h, s) for h, s in zip(ham, self._overlaps(ks))]
676 en = np.array([w for w, _ in solved], dtype='c16').reshape(len(ks), self.norb)
677 vn = np.array([v for _, v in solved], dtype='c16').reshape(ham.shape)
678 elif eigenvec:
679 en, vn = np.linalg.eig(ham)
680 else:
681 en = np.linalg.eigvals(ham)
682 ind = np.argsort(en.real, axis=1, kind='stable')
683 if not eigenvec:
684 return np.take_along_axis(en, ind, axis=1)
685 return np.take_along_axis(en, ind, axis=1), np.take_along_axis(vn, ind[:, None, :], axis=2)
687 def is_hermitian(self) -> bool:
688 '''
689 Check whether :math:`H(\\mathbf{k})` is Hermitian: it is, unless an
690 onsite term is not, or *set_hopping* was called with
691 ``hermitian=False`` (such a model is always treated as
692 non-Hermitian, and diagonalized with the general solver).
694 Returns
695 -------
696 hermitian : Boolean
697 '''
698 return bool(not self._nonreciprocal and np.all(self.onsite.imag == 0.)
699 and np.array_equal(self._onsite_offdiag, self._onsite_offdiag.conj().T))
701 def k_path(
702 self, points: list[ArrayLike], nk: int,
703 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
704 r'''
705 Build a k-path through a list of high-symmetry points, and get the
706 associated bands.
708 Parameters
709 ----------
710 points
711 List of at least two k-points (each a tuple/list of
712 *dim* real numbers), e.g. from *high_symmetry_path*.
713 nk : Positive integer
714 Number of k-points per path segment.
716 Returns
717 -------
718 ks_dist : Real ndarray
719 Cumulative distance along the path,
720 to be used as the x-axis of a band-structure plot.
721 en : Real ndarray, shape (len(ks_dist), norb)
722 Band energies.
723 '''
724 error_handling.k_path_points(points, self.dim)
725 error_handling.positive_int(nk, 'nk')
726 ks, _, nodes = _path(points, nk)
727 self.get_bands(ks)
728 self.nodes = np.array(nodes)
729 return self.ks_dist, self.en
731 def mesh_grid(
732 self, nk: int | tuple[int, ...], monkhorst_pack: bool = False,
733 ) -> tuple[list[NDArray[np.float64]], NDArray[np.float64]]:
734 r'''
735 Build a uniform grid of fractional coordinates spanning the
736 Brillouin zone, and the corresponding Cartesian k-points. By
737 default the grid contains :math:`\Gamma`: :math:`r/q`,
738 :math:`r = 0, \dots, q-1` (each in [0, 1)). With *monkhorst_pack*,
739 it is the grid of Monkhorst and Pack (Phys. Rev. B 13, 5188 (1976)),
740 :math:`(2r - q - 1)/2q`, :math:`r = 1, \dots, q` (each in
741 (-1/2, 1/2)), symmetric about :math:`\Gamma`: the same points for
742 odd :math:`q`, shifted by half a step (off :math:`\Gamma` and the
743 zone boundary) for even :math:`q`. Both integrate exactly every
744 Fourier component :math:`e^{i\mathbf{k}\cdot\mathbf{R}}` with
745 :math:`|R_i| < q`.
747 Parameters
748 ----------
749 nk : Positive integer, or tuple of *dim* positive integers
750 Number of k-points along each reciprocal lattice vector.
751 monkhorst_pack : Boolean, default False
752 If True, the
753 Monkhorst-Pack grid.
755 Returns
756 -------
757 fracs
758 List of *dim* real ndarrays, shape (nk1, nk2) each
759 (or (nk1,) in 1D): fractional coordinates of the grid.
760 ks : Real ndarray, shape (nk1*nk2, dim) (or (nk1, dim) in 1D)
761 '''
762 error_handling.nk(nk, self.dim)
763 error_handling.boolean(monkhorst_pack, 'monkhorst_pack')
764 if isinstance(nk, int):
765 nk = (nk,) * self.dim
766 if monkhorst_pack:
767 axes = [(2 * np.arange(1, n + 1) - n - 1) / (2 * n) for n in nk]
768 else:
769 axes = [np.arange(n) / n for n in nk]
770 fracs = np.meshgrid(*axes, indexing='ij')
771 ks = sum(f.ravel()[:, None] * b[None, :] for f, b in zip(fracs, self.rec_vec_k))
772 if self.dim == 1:
773 return [fracs[0]], ks
774 return list(fracs), ks
776 def mesh_bands(self, nk: int | tuple[int, ...], monkhorst_pack: bool = False) -> NDArray[np.float64]:
777 '''
778 Diagonalize :math:`H(\\mathbf{k})` over a uniform mesh spanning the
779 Brillouin zone. Unlike *get_bands*, it leaves the band structure kept
780 for *plot_bands* untouched.
782 Parameters
783 ----------
784 nk : Positive integer, or tuple of *dim* positive integers
785 Number of k-points along each reciprocal lattice vector.
786 monkhorst_pack : Boolean, default False
787 If True, the
788 Monkhorst-Pack grid (see *mesh_grid*).
790 Returns
791 -------
792 en : Real ndarray, shape (nk1*nk2, norb)
793 Band energies
794 over the mesh (flattened).
795 '''
796 _, ks = self.mesh_grid(nk, monkhorst_pack)
797 return self._diagonalize(ks)
799 def get_fermi_level(
800 self, n_electrons: float, nk: int | tuple[int, ...] = 30, temperature: float = 0.,
801 ) -> float:
802 '''
803 Get the Fermi level holding *n_electrons* electrons per unit cell
804 (one electron per band and k-point), from the bands over a uniform
805 Brillouin-zone mesh (see *tbkit.occupation.fermi_level*).
807 Parameters
808 ----------
809 n_electrons : Positive real or zero, at most norb
810 Electrons per unit cell.
811 nk : Positive integer, or tuple, default 30
812 k-mesh.
813 temperature : Positive real or zero, default 0
815 Returns
816 -------
817 e_fermi : Real number
818 '''
819 en = self.mesh_bands(nk)
820 return occupation.fermi_level(en, n_electrons, temperature,
821 weights=np.full(en.shape, 1. / len(en)))
823 def berry_curvature(
824 self, bands: int | list[int], nk: int | tuple[int, int] = 30,
825 plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
826 ) -> NDArray[np.float64]:
827 r'''
828 Get the Berry curvature of a group of bands over a uniform
829 Brillouin-zone mesh, using the gauge-invariant lattice method of
830 Fukui, Hatsugai and Suzuki (J. Phys. Soc. Jpn. 74, 1674 (2005)):
831 the flux through each mesh plaquette is minus the phase of the
832 product of the (Slater-determinant) overlaps between the occupied
833 subspaces at its four corners,
835 .. math::
837 \Omega_{\square} = -\arg\prod_{\langle\mathbf{k}\mathbf{k}'\rangle\in\partial\square}
838 \det\langle u(\mathbf{k})|u(\mathbf{k}')\rangle\, ,
840 the corners being visited counterclockwise. This is the discrete
841 version of :math:`\Omega = \partial_{k_x}A_y - \partial_{k_y}A_x`,
842 :math:`\mathbf{A} = i\langle u|\nabla_{\mathbf{k}}u\rangle`, and does
843 not depend on the orientation (handedness) of *prim_vec*. The local
844 curvature is that of the periodic gauge of *get_ham* (orbitals on the
845 cell origin), which differs point by point from the one with the
846 orbital positions included (see *quantum_geometric_tensor*); their
847 integral, the Chern number, is the same. For a
848 non-Hermitian :math:`H(\mathbf{k})` (see *set_onsite*), the bands
849 are ordered by the real part of their energy and the overlaps are
850 taken between right eigenvectors (the "RR" curvature of
851 *biorthogonal_berry_curvature*). That ordering follows the bands
852 continuously only if a real line gap separates *bands* from the
853 others over the whole zone: without one (an imaginary line gap,
854 or exceptional points), a ValueError is raised instead of a
855 meaningless result -- see *biorthogonal_berry_curvature*.
857 Parameters
858 ----------
859 bands : Positive integer, or list of positive integers
860 Band
861 index, or indices of a group of bands (e.g. all occupied bands
862 below a gap).
863 nk : Positive integer, or tuple of 2 positive integers, default 30
864 Number of k-points along each of the two
865 reciprocal lattice vectors spanning the plane.
866 plane : Tuple of two distinct integers
867 Default value (0, 1).
868 3D lattices only: the reciprocal vectors
869 :math:`\mathbf{b}_i, \mathbf{b}_j` spanning the plane of the
870 Brillouin zone. The flux is oriented along the remaining
871 primitive vector :math:`\mathbf{a}_l` (along +z in 2D).
872 k_fixed : Real number, default 0
873 3D lattices only: the
874 fractional coordinate of the plane along the remaining reciprocal
875 vector :math:`\mathbf{b}_l` (e.g. :math:`k_z` in units of
876 :math:`|\mathbf{b}_3|` for ``plane=(0, 1)``).
878 Returns
879 -------
880 curv : Real ndarray, shape (nk1, nk2)
881 Berry curvature
882 (flux through each plaquette, in radians). Summing *curv* and
883 dividing by :math:`2\pi` gives the Chern number, see
884 *chern_number*.
885 '''
886 error_handling.dim_min(self.dim, 2)
887 if isinstance(bands, int):
888 bands = [bands]
889 error_handling.band_indices(bands, self.norb)
890 if isinstance(nk, int):
891 nk = (nk, nk)
892 error_handling.nk(nk, 2)
893 error_handling.plane(plane, self.dim)
894 error_handling.real_number(k_fixed, 'k_fixed')
895 n1, n2 = nk
896 ks = self._plane_mesh(nk, plane, k_fixed)
897 if not self.is_hermitian():
898 # bands labelled by Re E are continuous only across a real line gap
899 en = self._eigs(ks.reshape(-1, self.dim))
900 error_handling.line_gap(en.real, bands, 'real')
901 error_handling.band_continuity(en.reshape(n1, n2, self.norb), bands)
902 v1 = self._subspaces(ks.reshape(-1, self.dim), bands).reshape(n1, n2, self.norb, len(bands))
903 return self._flux(v1, plane)
905 def _flux(self, v1: NDArray[np.complex128], plane: tuple[int, int]) -> NDArray[np.float64]:
906 '''
907 Private method. Fukui-Hatsugai-Suzuki flux through each plaquette of
908 the mesh of *_plane_mesh*, from orthonormal bases *v1* of a subspace
909 at its nodes (shape (n1, n2, norb, n)), oriented as in *berry_curvature*.
910 '''
911 # the other corners of each plaquette: k + b_i/n1, k + b_i/n1 + b_j/n2, k + b_j/n2
912 v2 = np.roll(v1, -1, axis=0)
913 v3 = np.roll(v2, -1, axis=1)
914 v4 = np.roll(v1, -1, axis=1)
916 def overlap(a, b):
917 return np.linalg.det(a.conj().swapaxes(-1, -2) @ b)
918 link = overlap(v1, v2) * overlap(v2, v3) * overlap(v3, v4) * overlap(v4, v1)
919 curv = -np.angle(link)
920 # The corners k, k+b_i/n1, k+b_i/n1+b_j/n2, k+b_j/n2 run
921 # counterclockwise about the normal only if (b_i x b_j).normal > 0;
922 # otherwise every flux comes out with the opposite sign.
923 return self._plane_orientation(plane) * curv
925 def _plane_mesh(
926 self, nk: tuple[int, int], plane: tuple[int, int], k_fixed: float,
927 ) -> NDArray[np.float64]:
928 '''
929 Private method. k-points, shape (n1, n2, dim), of the uniform mesh of
930 the Brillouin-zone plane spanned by b_i, b_j (plane = (i, j)), at
931 fractional coordinate k_fixed along the remaining reciprocal vector.
932 '''
933 n1, n2 = nk
934 i, j = plane
935 rest = [l for l in range(self.dim) if l not in plane]
936 fracs = np.zeros((n1, n2, self.dim))
937 fracs[:, :, i] = (np.arange(n1) / n1)[:, None]
938 fracs[:, :, j] = (np.arange(n2) / n2)[None, :]
939 for l in rest:
940 fracs[:, :, l] = k_fixed
941 return fracs @ self.rec_vec_k
943 def _plane_orientation(self, plane: tuple[int, int]) -> float:
944 '''
945 Private method. +1 if (b_i, b_j) is right-handed about the plane's
946 normal (+z in 2D; the remaining primitive vector in 3D), else -1.
947 '''
948 b = [np.zeros(3) for _ in range(2)]
949 for bb, idx in zip(b, plane):
950 bb[:self.space_dim] = self.rec_vec[idx]
951 if self.dim == 3:
952 normal = np.array(self.lat.prim_vec[3 - sum(plane)], dtype='f8')
953 else:
954 a = [np.zeros(3), np.zeros(3)]
955 for aa, vec in zip(a, self.lat.prim_vec):
956 aa[:self.space_dim] = vec
957 normal = np.cross(a[0], a[1])
958 if abs(normal[2]) > 1e-12:
959 normal = np.array([0., 0., 1.])
960 return float(np.sign(np.dot(np.cross(b[0], b[1]), normal)))
962 def _subspace(self, k: NDArray[np.float64], bands: list[int]) -> NDArray[np.complex128]:
963 '''
964 Private method. Orthonormal basis (columns) of the eigenvectors of
965 *bands* at *k* (in the periodic gauge, orbital positions ignored).
966 With an overlap S(k), the Lowdin-orthonormalized S^(1/2) v.
967 '''
968 return self._subspaces(np.asarray(k, dtype='f8')[None], bands)[0]
970 def _subspaces(self, ks: NDArray[np.float64], bands: list[int]) -> NDArray[np.complex128]:
971 '''
972 Private method. *_subspace* at the k-points *ks* (shape (nk, dim)),
973 shape (nk, norb, len(bands)). S^(1/2) comes from the eigenvalues of
974 the (Hermitian, positive) overlap.
975 '''
976 _, vn = self._eigs(ks, eigenvec=True)
977 if self._overlap_hop:
978 sv, us = np.linalg.eigh(self._overlaps(ks))
979 vn = (us * np.sqrt(sv)[:, None, :]) @ (us.conj().transpose(0, 2, 1) @ vn)
980 return vn[:, :, bands]
982 def biorthogonal_berry_curvature(
983 self, bands: int | list[int], nk: int | tuple[int, int] = 30, kind: str = 'LR',
984 gap: str = 'real', plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
985 ) -> NDArray[np.float64]:
986 r'''
987 Get the Berry curvature of a group of bands of a (non-Hermitian)
988 model from its right and left eigenvectors, :math:`H|R_n\rangle =
989 E_n|R_n\rangle` and :math:`H^\dagger|L_n\rangle = E_n^*|L_n\rangle`,
990 normalized by :math:`\langle L_m|R_n\rangle = \delta_{mn}` (Shen, Zhen
991 and Fu, Phys. Rev. Lett. 120, 146402 (2018)). The four connections
993 .. math::
995 \mathbf{A}^{\alpha\beta} = i\langle u^\alpha|\nabla_{\mathbf{k}}u^\beta\rangle\, ,
996 \qquad \alpha, \beta \in \{L, R\}\, ,
998 have different curvatures point by point, but their integrals, the
999 Chern numbers :math:`C^{LR} = C^{RL} = C^{RR} = C^{LL}`, are the same
1000 integer. The flux through each plaquette is computed as in
1001 *berry_curvature* (Fukui-Hatsugai-Suzuki), with the links
1002 :math:`\det\langle u^\alpha(\mathbf{k})|u^\beta(\mathbf{k}')\rangle` (for
1003 RR and LL, the phase does not depend on the normalization). For a
1004 Hermitian model all four are the curvature of *berry_curvature*.
1006 The bands are ordered at every k-point by the real part of their
1007 energy (*gap* = 'real') or by its imaginary part ('imaginary'), and
1008 *bands* must be separated from the others by a line gap of that kind
1009 over the whole zone -- a line :math:`\mathrm{Re}\,E = E_0` (or
1010 :math:`\mathrm{Im}\,E = E_0`) that no band crosses. Otherwise a
1011 ValueError is raised: the gap closes at exceptional points (see
1012 *tbkit.exceptional.find_exceptional_points*), where the Chern number
1013 can change.
1015 Parameters
1016 ----------
1017 bands : Positive integer, or list of positive integers
1018 Band
1019 indices in the order set by *gap*.
1020 nk : Positive integer, or tuple of 2 positive integers, default 30
1021 See *berry_curvature*.
1022 kind : 'LR', 'RL', 'RR' or 'LL', default 'LR'
1023 :math:`\alpha\beta` above.
1024 gap : 'real' or 'imaginary', default 'real'
1025 Kind of line
1026 gap, and the order of the bands.
1027 plane : Tuple of two integers
1028 Default value (0, 1). 3D only, see
1029 *berry_curvature*.
1030 k_fixed : Real number, default 0
1031 3D only, see *berry_curvature*.
1033 Returns
1034 -------
1035 curv : Real ndarray, shape (nk1, nk2)
1036 Flux through each
1037 plaquette, in radians (sign convention of *berry_curvature*).
1038 '''
1039 error_handling.dim_min(self.dim, 2)
1040 if isinstance(bands, int):
1041 bands = [bands]
1042 error_handling.band_indices(bands, self.norb)
1043 if isinstance(nk, int):
1044 nk = (nk, nk)
1045 error_handling.nk(nk, 2)
1046 error_handling.biorthogonal_kind(kind)
1047 error_handling.gap_kind(gap)
1048 error_handling.plane(plane, self.dim)
1049 error_handling.real_number(k_fixed, 'k_fixed')
1050 error_handling.no_overlap(self._overlap_hop)
1051 n1, n2 = nk
1052 ks = self._plane_mesh(nk, plane, k_fixed)
1053 w, vr = np.linalg.eig(self._hams(ks.reshape(n1 * n2, -1)))
1054 order = np.argsort(w.real if gap == 'real' else w.imag, axis=1, kind='stable')
1055 en = np.take_along_axis(w, order, axis=1).reshape(n1, n2, self.norb)
1056 vr = np.take_along_axis(vr, order[:, None, :], axis=2)
1057 # biorthonormal: <L_m|R_n> = delta_mn, the rows of R^-1
1058 right = vr[:, :, bands].reshape(n1, n2, self.norb, len(bands))
1059 left = np.linalg.inv(vr).conj().transpose(0, 2, 1)[:, :, bands].reshape(right.shape)
1060 error_handling.line_gap(en.reshape(-1, self.norb).real if gap == 'real'
1061 else en.reshape(-1, self.norb).imag, bands, gap)
1062 error_handling.band_continuity(en, bands)
1063 bra, ket = {'LR': (left, right), 'RL': (right, left),
1064 'RR': (right, right), 'LL': (left, left)}[kind]
1065 # one link per oriented edge, U = det<u^alpha(k)|u^beta(k + dk)>, and
1066 # its inverse when the edge is run backwards: <L_a|R_b> and <L_b|R_a>
1067 # are not conjugate, so using both would spoil the exact
1068 # cancellation of the shared edges (for RR, 1/U has the phase of U*).
1069 bra_h = bra.conj().swapaxes(-1, -2)
1070 link = np.stack([np.linalg.det(bra_h @ np.roll(ket, -1, axis=0)),
1071 np.linalg.det(bra_h @ np.roll(ket, -1, axis=1))])
1072 flux = (link[0] * np.roll(link[1], -1, axis=0)
1073 / (np.roll(link[0], -1, axis=1) * link[1]))
1074 return self._plane_orientation(plane) * -np.angle(flux)
1076 def biorthogonal_chern_number(
1077 self, bands: int | list[int], nk: int | tuple[int, int] = 30, kind: str = 'LR',
1078 gap: str = 'real', plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
1079 ) -> float:
1080 r'''
1081 Get the Chern number of a line-gapped group of bands of a
1082 (non-Hermitian) model, :math:`C^{\alpha\beta} = \frac{1}{2\pi}\int_{BZ}
1083 \Omega^{\alpha\beta}\,d^2k`, the same integer for the four kinds
1084 LR, RL, RR, LL (see *biorthogonal_berry_curvature*). It equals the
1085 Hermitian Chern number of *chern_number* as long as the line gap stays
1086 open.
1088 Parameters
1089 ----------
1090 bands
1091 See *biorthogonal_berry_curvature*.
1092 nk
1093 See *biorthogonal_berry_curvature*.
1094 kind : 'LR', 'RL', 'RR' or 'LL', default 'LR'
1095 gap : 'real' or 'imaginary', default 'real'
1096 plane
1097 See *berry_curvature*.
1098 k_fixed
1099 See *berry_curvature*.
1101 Returns
1102 -------
1103 chern : Real number, close to an integer
1104 '''
1105 return self.biorthogonal_berry_curvature(bands, nk, kind, gap, plane, k_fixed).sum() / (2*PI)
1107 def chern_number(
1108 self, bands: int | list[int], nk: int | tuple[int, int] = 30,
1109 plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
1110 ) -> float:
1111 r'''
1112 Get the Chern number of a group of bands:
1114 .. math::
1116 C = \frac{1}{2\pi}\int_{BZ} \Omega(\mathbf{k})\, d^2k
1118 an integer (up to the numerical precision set by *nk*) for a group
1119 of bands that is isolated from the rest of the spectrum by a gap
1120 everywhere in the Brillouin zone. See *berry_curvature*. For a
1121 non-Hermitian model, *bands* (ordered by Re E) must be isolated by a
1122 real line gap, or a ValueError is raised; *biorthogonal_chern_number*
1123 also handles imaginary line gaps and the left eigenvectors.
1125 Parameters
1126 ----------
1127 bands : Positive integer, or list of positive integers
1128 Band
1129 index, or indices of a group of bands (e.g. all occupied bands
1130 below a gap).
1131 nk : Positive integer, or tuple of 2 positive integers, default 30
1132 Number of k-points along each reciprocal
1133 lattice vector.
1134 plane : Tuple of two integers
1135 Default value (0, 1). 3D lattices
1136 only, see *berry_curvature*.
1137 k_fixed : Real number, default 0
1138 3D lattices only, see
1139 *berry_curvature*: the Chern number of the plane at this
1140 fractional coordinate (it jumps across a Weyl point).
1142 Returns
1143 -------
1144 chern : Real number, close to an integer
1145 '''
1146 return self.berry_curvature(bands, nk, plane, k_fixed).sum() / (2*PI)
1148 def sector_chern_numbers(
1149 self, op, bands: int | list[int], nk: int | tuple[int, int] = 30,
1150 plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
1151 ) -> tuple[float, float]:
1152 r'''
1153 Get the Chern numbers of the two sectors into which a Hermitian
1154 operator :math:`O` splits a group of bands: at every k-point, the
1155 projected operator :math:`P(\mathbf{k})\,O\,P(\mathbf{k})`
1156 (:math:`P` the projector on *bands*) is diagonalized within the
1157 bands, and its eigenvectors with positive and negative eigenvalues
1158 span the two sectors, whose Chern numbers :math:`C_\pm` follow as in
1159 *berry_curvature*. They are well defined as long as the spectrum of
1160 :math:`POP` keeps a gap around zero over the whole plane, even if
1161 :math:`O` does not commute with :math:`H` (Prodan, Phys. Rev. B 80,
1162 125327 (2009)); otherwise a ValueError is raised. See
1163 *spin_chern_number* and *mirror_chern_number*.
1165 Parameters
1166 ----------
1167 op
1168 Complex ndarray, shape (norb, norb), or callable of k
1169 returning one: the Hermitian operator :math:`O`.
1170 bands : Band index, or list of band indices
1171 nk : Positive integer, or tuple of 2 positive integers, default 30
1172 plane : Tuple of two integers
1173 Default value (0, 1). 3D only, see *berry_curvature*.
1174 k_fixed : Real number, default 0
1175 3D only, see *berry_curvature*.
1177 Returns
1178 -------
1179 c_plus
1180 Real number, close to an integer: Chern number
1181 of the sector with :math:`POP > 0`.
1182 c_minus : Real number: that of the sector with :math:`POP < 0`
1183 :math:`C_+ + C_-` is the Chern number of *bands*.
1184 '''
1185 bands, nk, ks = self._check_sectors(bands, nk, plane, k_fixed)
1186 if callable(op):
1187 ops = np.array([self._operator(op, k) for k in ks])
1188 else:
1189 ops = self._operator(op, ks[0])[None]
1190 error_handling.hermitian_operator(ops)
1191 return self._sector_cherns(ops, bands, nk, ks, plane)
1193 def _check_sectors(
1194 self, bands, nk, plane, k_fixed,
1195 ) -> tuple[list[int], tuple[int, int], NDArray[np.float64]]:
1196 '''
1197 Private method. Validate the arguments of the sector Chern numbers;
1198 return them with the k-points of the plane mesh, shape (n1 n2, dim).
1199 '''
1200 error_handling.dim_min(self.dim, 2)
1201 if isinstance(bands, int):
1202 bands = [bands]
1203 error_handling.band_indices(bands, self.norb)
1204 if isinstance(nk, int):
1205 nk = (nk, nk)
1206 error_handling.nk(nk, 2)
1207 error_handling.plane(plane, self.dim)
1208 error_handling.real_number(k_fixed, 'k_fixed')
1209 error_handling.hermitian_kspace(self.is_hermitian(), self._overlap_hop)
1210 return bands, nk, self._plane_mesh(nk, plane, k_fixed).reshape(-1, self.dim)
1212 def _sector_cherns(
1213 self, ops: NDArray[np.complex128], bands: list[int], nk: tuple[int, int],
1214 ks: NDArray[np.float64], plane: tuple[int, int],
1215 ) -> tuple[float, float]:
1216 '''
1217 Private method. *sector_chern_numbers* with the Hermitian operator
1218 evaluated at the k-points *ks* (shape (n1 n2, norb, norb), or
1219 (1, norb, norb) if constant).
1220 '''
1221 n1, n2 = nk
1222 v = self._subspaces(ks, bands)
1223 w, u = np.linalg.eigh(v.conj().transpose(0, 2, 1) @ ops @ v)
1224 n_minus = error_handling.projected_gap(w)
1225 sectors = v @ u
1226 c_minus = self._flux(sectors[:, :, :n_minus].reshape(n1, n2, self.norb, -1), plane)
1227 c_plus = self._flux(sectors[:, :, n_minus:].reshape(n1, n2, self.norb, -1), plane)
1228 return float(c_plus.sum() / (2*PI)), float(c_minus.sum() / (2*PI))
1230 def spin_chern_number(
1231 self, bands: list[int], nk: int | tuple[int, int] = 30, s_z=None,
1232 plane: tuple[int, int] = (0, 1), k_fixed: float = 0.,
1233 ) -> float:
1234 r'''
1235 Get the spin Chern number of a group of bands (Prodan, Phys. Rev. B
1236 80, 125327 (2009); Sheng, Weng, Sheng and Haldane, Phys. Rev. Lett.
1237 97, 036808 (2006)),
1239 .. math::
1241 C_s = \frac{C_+ - C_-}{2}\, ,
1243 with :math:`C_\pm` the Chern numbers of the two sectors of the
1244 projected spin operator :math:`Ps_zP` (see *sector_chern_numbers*).
1245 When :math:`s_z` is conserved they are the Chern numbers of the two
1246 spins; Rashba coupling mixes the spins, but :math:`C_s` stays
1247 quantized, and unchanged, as long as the spectrum of :math:`Ps_zP`
1248 keeps a gap around zero. With time reversal, :math:`C_s` modulo 2 is
1249 the :math:`\mathbb{Z}_2` invariant of *z2_invariant*.
1251 Parameters
1252 ----------
1253 bands : List of band indices (e.g. the occupied bands)
1254 nk : Positive integer, or tuple of 2 positive integers, default 30
1255 s_z : Complex ndarray, shape (norb, norb), or callable of k
1256 Default value None: :math:`\sigma_z` on every site of a model
1257 built with ``spin=True``. The spin operator.
1258 plane : Tuple of two integers
1259 Default value (0, 1). 3D only, see *berry_curvature*.
1260 k_fixed : Real number, default 0
1261 3D only, see *berry_curvature*.
1263 Returns
1264 -------
1265 c_s : Real number, close to an integer
1266 '''
1267 if s_z is None:
1268 error_handling.spinful(self.spin)
1269 s_z = np.kron(np.eye(self.norb // 2), PAULI['z'])
1270 c_plus, c_minus = self.sector_chern_numbers(s_z, bands, nk, plane, k_fixed)
1271 return (c_plus - c_minus) / 2
1273 def mirror_chern_number(
1274 self, mirror, bands: list[int], nk: int | tuple[int, int] = 30,
1275 plane: tuple[int, int] = (0, 1), k_fixed: float = 0., tol: float = 1e-8,
1276 ) -> float:
1277 r'''
1278 Get the mirror Chern number of a group of bands on a mirror-invariant
1279 plane of the Brillouin zone (Teo, Fu and Kane, Phys. Rev. B 78,
1280 045426 (2008); Hsieh et al., Nat. Commun. 3, 982 (2012)),
1282 .. math::
1284 n_M = \frac{C_{+} - C_{-}}{2}\, ,
1286 with :math:`C_\pm` the Chern numbers of the bands of mirror
1287 eigenvalue :math:`+i` and :math:`-i` (:math:`M^2 = -1`, spinful
1288 electrons) or :math:`+1` and :math:`-1` (:math:`M^2 = +1`). The
1289 mirror must commute with :math:`H(\mathbf{k})` at every point of the
1290 plane: in 3D, the plane fixed by the reflection (e.g.
1291 :math:`k_z = 0` or :math:`\pi` for :math:`z\to-z`); in 2D, the whole
1292 zone, for a reflection of the third axis. A topological crystalline
1293 insulator such as SnTe has :math:`n_M = 2` and every
1294 :math:`\mathbb{Z}_2` index zero (:math:`\nu = n_M` modulo 2 on the plane).
1296 Parameters
1297 ----------
1298 mirror
1299 Complex ndarray, shape (norb, norb), or callable of k:
1300 the unitary mirror operator :math:`M`, with :math:`M^2 = \pm1`.
1301 bands : List of band indices (e.g. the occupied bands)
1302 nk : Positive integer, or tuple of 2 positive integers, default 30
1303 plane : Tuple of two integers
1304 Default value (0, 1). 3D only, see *berry_curvature*.
1305 k_fixed : Real number, default 0
1306 3D only: the
1307 mirror-invariant plane (0 or 0.5 for a reflection of the remaining axis).
1308 tol : Positive real, default 1e-8
1309 Tolerance of the
1310 checks :math:`M^2 = \pm1` and :math:`[M, H(\mathbf{k})] = 0`.
1312 Returns
1313 -------
1314 n_m : Real number, close to an integer
1315 '''
1316 bands, nk, ks = self._check_sectors(bands, nk, plane, k_fixed)
1317 error_handling.positive_real(tol, 'tol')
1318 if callable(mirror):
1319 mats = np.array([self._operator(mirror, k) for k in ks])
1320 else:
1321 mats = self._operator(mirror, ks[0])[None]
1322 sign = error_handling.mirror_square(mats @ mats, tol)
1323 hams = self._hams(ks)
1324 error_handling.commutes(float(np.max(np.abs(mats @ hams - hams @ mats))), tol)
1325 # Hermitian, with eigenvalue +1 on the sector M = +i (M = +1 if M^2 = 1)
1326 herm = mats if sign > 0 else -1j * mats
1327 c_plus, c_minus = self._sector_cherns(herm, bands, nk, ks, plane)
1328 return (c_plus - c_minus) / 2
1330 # ------------------------------------------------------------------
1331 # Anomalous and spin Hall conductivities (Kubo formula)
1332 # ------------------------------------------------------------------
1334 def _check_static(self) -> None:
1335 '''
1336 Private method. Methods that differentiate the stored hoppings call
1337 it first (a Floquet model overrides it to refuse).
1338 '''
1340 def _bloch_sum(
1341 self, hops: list, k_cart: NDArray[np.float64], directions: Sequence[NDArray[np.float64]] = (),
1342 positions: bool = False, A: NDArray[np.float64] | None = None,
1343 ) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
1344 r'''
1345 Private method. The Bloch sum :math:`\sum t\,e^{i\mathbf{k}\cdot\mathbf{d}}`
1346 of a hopping list over many k-points at once (*k_cart*, shape
1347 (nk, space_dim)), and its derivatives
1348 :math:`\sum i(\mathbf{u}\cdot\mathbf{d})\,t\,e^{i\mathbf{k}\cdot\mathbf{d}}`
1349 along the Cartesian unit vectors *directions*, with
1350 :math:`\mathbf{d} = \mathbf{R}` (plus :math:`\boldsymbol\tau_j-\boldsymbol\tau_i`
1351 with *positions*). A vector potential *A* multiplies each term by the
1352 Peierls phase :math:`e^{i\mathbf{A}\cdot(\mathbf{R}+\boldsymbol\tau_j-\boldsymbol\tau_i)}`
1353 (see *get_ham_peierls*). Shapes (nk, norb, norb) and
1354 (len(directions), nk, norb, norb). Every Bloch matrix of the model
1355 (*get_ham*, *get_overlap*, *get_ham_peierls*, *_bloch_derivatives*)
1356 is built here.
1358 When the hoppings fill at least *group_fill* of the blocks
1359 :math:`T_{\mathbf{R}}` (one norb x norb matrix per distinct
1360 :math:`\mathbf{R}`), the sum is one matrix product
1361 :math:`\sum_{\mathbf{R}} e^{i\mathbf{k}\cdot\mathbf{R}}\,T_{\mathbf{R}}`,
1362 with the orbital positions as the phases
1363 :math:`e^{-i\mathbf{k}\cdot\boldsymbol\tau_i}` and
1364 :math:`e^{i\mathbf{k}\cdot\boldsymbol\tau_j}` of its rows and
1365 columns. Otherwise, each hopping is summed on its own.
1366 '''
1367 n, nk = self.norb, len(k_cart)
1368 i = np.array([h[0] for h in hops], dtype=int)
1369 j = np.array([h[1] for h in hops], dtype=int)
1370 d = np.array([h[2] for h in hops], dtype='f8').reshape(len(hops), self.space_dim)
1371 t = np.array([h[3] for h in hops], dtype='c16')
1372 if positions or A is not None:
1373 tau = self.orbital_positions()
1374 bond = d + tau[j] - tau[i]
1375 if A is not None:
1376 # the Peierls phase of a bond does not depend on k
1377 t = t * np.exp(1j * (bond @ A))
1378 if len(hops):
1379 cells, cell = np.unique(d, axis=0, return_inverse=True)
1380 if len(hops) >= self.group_fill * len(cells) * n * n:
1381 return self._bloch_grouped(i, j, t, cells, cell.ravel(), k_cart, directions,
1382 tau if positions else None)
1383 if positions:
1384 d = bond
1385 w = np.exp(1j * (k_cart @ d.T)) * t[None, :]
1386 # accumulate term m of k-point p in entry (p, i_m, j_m): repeated
1387 # (i, j) pairs add up, as they must
1388 index = (np.arange(nk)[:, None] * n * n + i * n + j).ravel()
1390 def scatter(v):
1391 size = nk * n * n
1392 return (np.bincount(index, v.real.ravel(), size)
1393 + 1j * np.bincount(index, v.imag.ravel(), size)).reshape(nk, n, n)
1394 mat = scatter(w)
1395 dmat = np.array([scatter(w * (1j * (d @ u))[None, :]) for u in directions],
1396 dtype='c16').reshape(len(directions), nk, n, n)
1397 return mat, dmat
1399 def _bloch_grouped(
1400 self, i: NDArray[np.int_], j: NDArray[np.int_], t: NDArray[np.complex128],
1401 cells: NDArray[np.float64], cell: NDArray[np.int_], k_cart: NDArray[np.float64],
1402 directions: Sequence[NDArray[np.float64]], tau: NDArray[np.float64] | None,
1403 ) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
1404 r'''
1405 Private method. *_bloch_sum* as one matrix product over the distinct
1406 lattice vectors *cells* (hopping m connects to cells[cell[m]]), with
1407 the orbital positions *tau* in the phases (None: periodic gauge).
1408 The derivative along :math:`\mathbf{u}` weights each block by
1409 :math:`i\mathbf{u}\cdot\mathbf{R}`, plus
1410 :math:`i\mathbf{u}\cdot(\boldsymbol\tau_j-\boldsymbol\tau_i)` with *tau*.
1411 '''
1412 n, nk = self.norb, len(k_cart)
1413 blocks = np.zeros((len(cells), n * n), dtype='c16')
1414 np.add.at(blocks, (cell, i * n + j), t)
1415 phases = np.exp(1j * (k_cart @ cells.T))
1416 mat = (phases @ blocks).reshape(nk, n, n)
1417 dmat = np.empty((len(directions), nk, n, n), dtype='c16')
1418 for a, u in enumerate(directions):
1419 dmat[a] = ((phases * (1j * (cells @ u))[None, :]) @ blocks).reshape(nk, n, n)
1420 if tau is not None:
1421 x = tau @ u
1422 dmat[a] += 1j * (x[None, :] - x[:, None])[None] * mat
1423 if tau is not None:
1424 orb = np.exp(1j * (k_cart @ tau.T))
1425 row, col = orb.conj()[:, :, None], orb[:, None, :]
1426 mat = row * mat * col
1427 dmat = row[None] * dmat * col[None]
1428 return mat, dmat
1430 def _bloch_derivatives(
1431 self, ks: NDArray[np.float64], directions: list[NDArray[np.float64]], positions: bool,
1432 params: dict | None = None,
1433 ) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
1434 r'''
1435 Private method. :math:`H(\mathbf{k})` and its exact derivatives
1436 :math:`\partial H/\partial k_u` along the Cartesian unit vectors
1437 *directions*, at the k-points *ks* (shape (nk, dim), coordinates of
1438 *get_ham*). With an overlap, the Lowdin-orthogonalized
1439 :math:`H' = S^{-1/2}HS^{-1/2}` (the eigenvectors :math:`S^{1/2}v` of
1440 the topology tools) and its derivative, with
1441 :math:`\partial S^{-1/2}` from the Daleckii-Krein formula in the
1442 eigenbasis :math:`S = \sum_a s_a|a\rangle\langle a|`:
1443 :math:`\langle a|\partial S^{-1/2}|b\rangle = \langle a|\partial S|b\rangle
1444 (s_a^{-1/2}-s_b^{-1/2})/(s_a-s_b)`, :math:`-s_a^{-3/2}/2` when :math:`s_a = s_b`.
1445 The value functions take *params* (default: *self.params*).
1446 '''
1447 k_cart = ks @ self.k_basis.T
1448 hops = self._hop if params is None else self._hops(params)
1449 ham, dham = self._bloch_sum(hops, k_cart, directions, positions)
1450 ham = ham + (np.diag(self.onsite) + self._onsite_offdiag)[None]
1451 if not self._overlap_hop:
1452 return ham, dham
1453 s, ds = self._bloch_sum(self._overlap_hop, k_cart, directions, positions)
1454 s = s + np.eye(self.norb)[None]
1455 sv, vs = np.linalg.eigh(s)
1456 error_handling.positive_overlap(sv.min())
1457 r = sv ** -0.5
1458 diff = sv[:, :, None] - sv[:, None, :]
1459 close = np.abs(diff) < 1e-10
1460 lk = np.where(close, -0.5 * sv[:, :, None] ** -1.5,
1461 (r[:, :, None] - r[:, None, :]) / np.where(close, 1., diff))
1462 vh = vs.conj().transpose(0, 2, 1)
1463 s_inv = (vs * r[:, None, :]) @ vh
1464 ds_inv = vs @ ((vh @ ds @ vs) * lk) @ vh
1465 dham = ds_inv @ ham @ s_inv + s_inv @ dham @ s_inv + s_inv @ ham @ ds_inv
1466 return s_inv @ ham @ s_inv, dham
1468 def _kubo(
1469 self, ks: NDArray[np.float64], e_fermi: NDArray[np.float64], temperature: float,
1470 directions: list[NDArray[np.float64]], pairs: list[tuple[int, int]], positions: bool,
1471 spin_op: NDArray[np.complex128] | None, response: str = 'hall',
1472 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
1473 r'''
1474 Private method. The Kubo curvature
1475 :math:`F(\mathbf{k}) = \sum_{n\neq m}P_{nm}\,\mathrm{Im}[X_{nm}Y_{mn}]/(E_n-E_m)^2`
1476 at each k-point, Fermi energy and pair (X, Y) of velocities (X the
1477 spin current :math:`\{s, v\}/2` if *spin_op*), with the weights
1478 :math:`P_{nm}` of *response* (see *_response_weights*), e.g.
1479 :math:`P_{nm} = -(f_n-f_m)` for the Hall conductivity
1480 (:math:`F = \sum_n f_n\Omega_n`). Shape (nk, len(e_fermi), len(pairs)).
1481 Pairs of degenerate states (equal weights) are skipped, so :math:`F`
1482 stays finite at band crossings below or above the Fermi level. Also
1483 returns, per k-point, the occupation-independent size of the
1484 curvature, :math:`\sum_{n\neq m}|\mathrm{Im}[X_{nm}Y_{mn}]|/(E_n-E_m)^2`,
1485 which locates the hot spots of the adaptive refinement.
1486 '''
1487 out = np.zeros((len(ks), len(e_fermi), len(pairs)))
1488 hot = np.zeros(len(ks))
1489 chunk = max(1, 400000 // self.norb ** 2)
1490 for c0 in range(0, len(ks), chunk):
1491 ham, dham = self._bloch_derivatives(ks[c0:c0 + chunk], directions, positions)
1492 en, vec = np.linalg.eigh(ham)
1493 vh = vec.conj().transpose(0, 2, 1)
1494 vel = [vh @ dh @ vec for dh in dham]
1495 de = en[:, :, None] - en[:, None, :]
1496 degenerate = np.abs(de) < 1e-9
1497 inv2 = np.where(degenerate, 0., 1. / np.where(degenerate, 1., de) ** 2)
1498 if spin_op is not None:
1499 s_e = vh @ spin_op @ vec
1500 for p, (a, b) in enumerate(pairs):
1501 x = vel[a] if spin_op is None else 0.5 * (s_e @ vel[a] + vel[a] @ s_e)
1502 w = np.imag(x * vel[b].transpose(0, 2, 1)) * inv2
1503 hot[c0:c0 + chunk] += np.abs(w).sum(axis=(1, 2))
1504 for e, mu in enumerate(e_fermi):
1505 weights = _response_weights(en, float(mu), temperature, response)
1506 out[c0:c0 + chunk, e, p] = np.einsum('knm,knm->k', weights, w)
1507 return out, hot
1509 def _hall(
1510 self, e_fermi, temperature, nk, plane, k_fixed, positions, refine, refine_fraction,
1511 spin_axis, response='hall',
1512 ):
1513 '''
1514 Private method. Validate, set the geometry and mesh, integrate the
1515 Kubo curvature (with the adaptive refinement), for
1516 *hall_conductivity* (spin_axis None) and *spin_hall_conductivity*,
1517 and, with the weights of *response*, for *orbital_magnetization*,
1518 *anomalous_nernst_conductivity* and *thermal_hall_conductivity*.
1519 '''
1520 self._check_static()
1521 error_handling.dim_min(self.dim, 2)
1522 error_handling.hermitian_model(self.is_hermitian())
1523 error_handling.fermi_energies(e_fermi)
1524 if response in ('nernst', 'thermal'):
1525 error_handling.positive_real(temperature, 'temperature')
1526 else:
1527 error_handling.positive_real_zero(temperature, 'temperature')
1528 error_handling.plane(plane, self.dim)
1529 error_handling.k_fixed_hall(k_fixed, self.dim, spin_axis is None)
1530 error_handling.boolean(positions, 'positions')
1531 error_handling.positive_int(refine, 'refine')
1532 error_handling.refine_fraction(refine_fraction)
1533 full = self.dim == 3 and k_fixed is None
1534 axes = [0, 1, 2] if full else list(plane)
1535 if isinstance(nk, int):
1536 nk = (nk,) * len(axes)
1537 error_handling.nk(nk, len(axes))
1538 spin_op = None
1539 if spin_axis is not None:
1540 error_handling.spinful(self.spin)
1541 error_handling.spin_axis(spin_axis)
1542 spin_op = np.kron(np.eye(self.n_sites), PAULI[spin_axis] / 2.)
1543 # Directions (Cartesian unit vectors) of the velocities, the pairs
1544 # (u, w) whose curvature Omega_uw = Omega.(u x w) is integrated, and
1545 # the prefactor turning the mesh average into the conductivity.
1546 if full:
1547 directions = list(np.eye(3))
1548 pairs = [(1, 2), (2, 0), (0, 1)]
1549 volume = abs(np.linalg.det(np.array(self.lat.prim_vec, dtype='f8')))
1550 prefactor = 2 * PI / volume
1551 elif self.dim == 3:
1552 i, j = plane
1553 a_l = np.array(self.lat.prim_vec[3 - i - j], dtype='f8')
1554 b_i, b_j = (np.array(self.rec_vec[p], dtype='f8') for p in plane)
1555 u = b_i / np.linalg.norm(b_i)
1556 directions = [u, np.cross(a_l / np.linalg.norm(a_l), u)]
1557 pairs = [(0, 1)]
1558 prefactor = np.linalg.norm(np.cross(b_i, b_j)) / (2 * PI)
1559 else:
1560 directions = [self.k_basis[:, 0], self.k_basis[:, 1]]
1561 pairs = [(0, 1)]
1562 det = np.linalg.det(self.rec_vec_k)
1563 prefactor = self._plane_orientation((0, 1)) * np.sign(det) * abs(det) / (2 * PI)
1564 e_arr = np.atleast_1d(np.asarray(e_fermi, dtype='f8')).ravel()
1565 grids = np.meshgrid(*[np.arange(n) / n for n in nk], indexing='ij')
1566 fracs = np.zeros((grids[0].size, self.dim))
1567 for ax, g in zip(axes, grids):
1568 fracs[:, ax] = g.ravel()
1569 for ax in range(self.dim):
1570 if ax not in axes:
1571 fracs[:, ax] = k_fixed
1572 args = (e_arr, temperature, directions, pairs, positions, spin_op, response)
1573 coarse, size = self._kubo(fracs @ self.rec_vec_k, *args)
1574 total = coarse.sum(axis=0)
1575 if refine > 1:
1576 # Wang, Yates, Souza and Vanderbilt (2006): resample the cells
1577 # where the curvature peaks on a refine^d submesh (midpoints of
1578 # the sub-cells), and replace their coarse value by the average.
1579 # The cells are ranked by the size of the curvature whatever the
1580 # occupations, so that cells straddling the Fermi surface are
1581 # refined together with their neighbours.
1582 n_hot = int(np.ceil(refine_fraction * len(fracs)))
1583 hot = np.argsort(size)[::-1][:n_hot]
1584 sub = np.array(np.meshgrid(*[(np.arange(refine) + 0.5) / refine - 0.5] * len(axes),
1585 indexing='ij')).reshape(len(axes), -1).T
1586 offsets = np.zeros((len(sub), self.dim))
1587 for c, (ax, n) in enumerate(zip(axes, nk)):
1588 offsets[:, ax] = sub[:, c] / n
1589 fine = (fracs[hot][:, None, :] + offsets[None, :, :]).reshape(-1, self.dim)
1590 fine_val = self._kubo(fine @ self.rec_vec_k, *args)[0].reshape(n_hot, len(sub), len(e_arr), len(pairs))
1591 total = total + (fine_val.mean(axis=1) - coarse[hot]).sum(axis=0)
1592 sigma = prefactor * total / len(fracs)
1593 if not full:
1594 sigma = sigma[:, 0]
1595 if np.ndim(e_fermi) == 0:
1596 return sigma[0] if full else float(sigma[0])
1597 return sigma.reshape(np.shape(e_fermi) + sigma.shape[1:])
1599 def hall_conductivity(
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 intrinsic anomalous Hall conductivity at any Fermi level --
1606 in a gap, in a band, or at a band crossing -- from the Kubo formula
1607 over a uniform Brillouin-zone mesh,
1609 .. math::
1611 \sigma_{xy} = \frac{e^2}{\hbar}\int_{BZ}\frac{d^dk}{(2\pi)^d}\,
1612 F(\mathbf{k})\, ,\qquad
1613 F = -\sum_{n\neq m}(f_n - f_m)\,
1614 \frac{\mathrm{Im}\left[\langle n|\partial_x H|m\rangle\langle m|\partial_y H|n\rangle\right]}
1615 {(E_n - E_m)^2} = \sum_n f_n\,\Omega_n\, ,
1617 with :math:`f_n = f(E_n)` the Fermi-Dirac occupations and
1619 .. math::
1621 \Omega_n = -2\,\mathrm{Im}\sum_{m\neq n}
1622 \frac{\langle n|\partial_x H|m\rangle\langle m|\partial_y H|n\rangle}{(E_n-E_m)^2}
1623 = -2\,\mathrm{Im}\langle\partial_x u_n|\partial_y u_n\rangle
1625 the Berry curvature of band *n* (:math:`\partial_\mu H = \partial H/\partial k_\mu`,
1626 the convention :math:`\Omega = \partial_xA_y - \partial_yA_x`,
1627 :math:`\mathbf{A} = i\langle u|\nabla_{\mathbf{k}}u\rangle` of
1628 *berry_curvature*). The form with :math:`f_n - f_m` stays finite
1629 where bands touch: degenerate pairs have equal occupations and
1630 drop out, instead of cancelling between two diverging
1631 :math:`\Omega_n`. :math:`\partial H/\partial\mathbf{k}` is built
1632 analytically from the stored hoppings, :math:`\sum_{\mathbf{R}}
1633 i\mathbf{d}\,t\,e^{i\mathbf{k}\cdot\mathbf{d}}` with the Cartesian bond
1634 vector :math:`\mathbf{d}`.
1636 **Sign.** In a gap, :math:`\sigma_{xy} = C\,e^2/h` with :math:`C`
1637 the value *chern_number* returns for the occupied bands: this is the
1638 sign of TKNN (Phys. Rev. Lett. 49, 405 (1982)), whose
1639 :math:`\sigma_H = (e^2/h)\,\frac{i}{2\pi}\int d^2k\,(\langle\partial_1u|\partial_2u\rangle
1640 - \langle\partial_2u|\partial_1u\rangle)` has the same integrand,
1641 :math:`\Omega/2\pi`. The adiabatic (semiclassical) response of a
1642 charge :math:`q` of either sign, :math:`\dot{\mathbf{r}} = \partial_{\mathbf{k}}E
1643 - \dot{\mathbf{k}}\times\boldsymbol\Omega` with
1644 :math:`\hbar\dot{\mathbf{k}} = q\mathbf{E}`, gives the opposite sign
1645 for the tensor of Ohm's law :math:`j_x = \sigma^{\mathrm{Ohm}}_{xy}E_y`:
1646 :math:`\sigma^{\mathrm{Ohm}}_{xy} = -(q^2/\hbar)\int F\,d^dk/(2\pi)^d`,
1647 the convention of Wang, Yates, Souza and Vanderbilt (Phys. Rev. B
1648 74, 195118 (2006)) and of the Kubo-Bastin formula (*tbkit.kpm.hall_conductivity*
1649 returns, likewise, minus it). So the value returned is
1650 :math:`\sigma^{\mathrm{Ohm}}_{yx} = -\sigma^{\mathrm{Ohm}}_{xy}`; e.g. a
1651 massive Dirac cone :math:`v(k_x\sigma_x + k_y\sigma_y) + m\sigma_z`
1652 gives :math:`+\frac{e^2}{2h}\frac{m}{|E_F|}` outside its gap (the
1653 literature's :math:`\sigma_{xy} = -\frac{e^2}{2h}\frac{m}{|E_F|}`), and
1654 :math:`\frac{e^2}{2h}\mathrm{sgn}(m)` inside.
1656 **Orbital positions.** The velocity is :math:`\partial H/\partial\mathbf{k}`
1657 in the gauge of *positions*: with ``positions=True`` (the default)
1658 the bond vector is :math:`\mathbf{d} = \mathbf{R} + \boldsymbol\tau_j
1659 - \boldsymbol\tau_i`, the matrix of the physical velocity
1660 :math:`i[H, \mathbf{r}]` of a tight-binding model whose position
1661 operator is diagonal on its orbitals (the one *kpm.hall_conductivity*
1662 uses on a torus). With ``positions=False``, every orbital sits on
1663 its cell origin (the periodic gauge of *get_ham*). The two differ by
1664 :math:`i(E_n-E_m)\langle n|\boldsymbol\tau|m\rangle` in the
1665 velocity's matrix elements, so the band curvatures differ by
1666 :math:`\nabla\times\langle n|\boldsymbol\tau|n\rangle`, a total
1667 derivative: the conductivities agree whenever the bands below
1668 :math:`E_F` are all full (in a gap, at any temperature well below
1669 the gap) and differ in general at partial filling, by the circulation
1670 :math:`\frac{1}{2\pi}\oint_{FS}\langle n|\boldsymbol\tau|n\rangle\cdot d\mathbf{k}`
1671 around the Fermi surface. The default is the physical one.
1673 **Units and dimensions.** In 2D, :math:`\sigma` is in units of
1674 :math:`e^2/h` (:math:`\frac{1}{2\pi}\int d^2k\,F`). In 3D, with
1675 *k_fixed* given, the value of the plane (the conductance of that
1676 slice of the Brillouin zone, in :math:`e^2/h`, oriented like
1677 *chern_number* along :math:`\mathbf{a}_l`); with ``k_fixed=None``,
1678 the full 3D conductivity as the Hall (pseudo)vector
1679 :math:`(\sigma_{yz}, \sigma_{zx}, \sigma_{xy}) = \frac{1}{(2\pi)^2}\int d^3k\,\mathbf{F}`
1680 in units of :math:`e^2/(h\cdot\mathrm{length})` (length the unit of
1681 *prim_vec*): :math:`\frac{1}{2\pi}\sum_l\bar C_l\,\mathbf{b}_l`, e.g.
1682 :math:`C/c` along :math:`z` for a stack of Chern layers :math:`c` apart.
1684 **Convergence.** The integrand peaks where bands nearly touch close
1685 to :math:`E_F` (hot spots), and jumps across the Fermi surface. A
1686 *temperature* smooths the latter. *refine* resamples the fraction
1687 *refine_fraction* of mesh cells with the largest :math:`|F|` on a
1688 finer submesh (*refine* points per direction), in the spirit of the
1689 adaptive mesh refinement of Wang, Yates, Souza and Vanderbilt
1690 (2006); the cells are ranked by the size of the curvature whatever
1691 the occupations. It is meant for metals: in a gap the plain uniform
1692 mesh already converges exponentially, and refining part of it
1693 costs a little accuracy (about :math:`4\cdot10^{-4}` for the
1694 Haldane model on a 60 x 60 mesh). With an overlap (*set_overlap*), :math:`H` is first
1695 Lowdin-orthogonalized, :math:`S^{-1/2}HS^{-1/2}`, as by the topology
1696 tools (exact derivatives of :math:`S^{-1/2}` included).
1698 Parameters
1699 ----------
1700 e_fermi : Real number, or array of real numbers
1701 Fermi energies
1702 (all of them are computed from one diagonalization of the mesh).
1703 temperature : Positive real or zero, default 0
1704 In
1705 energy units (:math:`k_B = 1`).
1706 nk
1707 Positive integer, or tuple (two integers for a plane,
1708 three for the full 3D conductivity). Default value 60. k-mesh.
1709 plane : Tuple of two integers
1710 Default value (0, 1). 3D only:
1711 the reciprocal vectors spanning the plane (see *berry_curvature*).
1712 k_fixed : Real number or None, default None
1713 3D only:
1714 the fractional coordinate of the plane along the remaining
1715 reciprocal vector; None for the full 3D conductivity. Ignored in 2D.
1716 positions : Boolean, default True
1717 Include the orbital
1718 positions in the bond vectors (see above).
1719 refine : Positive integer
1720 Default value 1 (no refinement).
1721 Submesh points per direction in the refined cells.
1722 refine_fraction : Real in (0, 1], default 0.05
1723 Fraction
1724 of the mesh cells that are refined.
1726 Returns
1727 -------
1728 sigma
1729 Real number (or ndarray shaped like *e_fermi*); in
1730 3D with ``k_fixed=None``, ndarray of shape (3,) (or
1731 ``e_fermi.shape + (3,)``).
1733 Examples
1734 --------
1735 Example usage::
1737 e_f = np.linspace(-3., 3., 301)
1738 sigma = hal.hall_conductivity(e_f, nk=120) # plateau at C in the gap
1739 '''
1740 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1741 refine_fraction, None)
1743 def spin_hall_conductivity(
1744 self, e_fermi: float | ArrayLike, temperature: float = 0., nk: int | tuple[int, int] = 60,
1745 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1746 spin_axis: str = 'z', refine: int = 1, refine_fraction: float = 0.05,
1747 ) -> float | NDArray[np.float64]:
1748 r'''
1749 Get the intrinsic spin Hall conductivity of a spinful model
1750 (``spin=True``): the Kubo formula of *hall_conductivity* with the
1751 charge current along :math:`x` replaced by the spin current
1752 (Sinova et al., Phys. Rev. Lett. 92, 126603 (2004))
1754 .. math::
1756 j^{s}_x = \tfrac12\{s, v_x\}\, ,\qquad s = \tfrac{\hbar}{2}\sigma_{a}\, ,
1758 i.e. :math:`\sigma^{s}_{xy} = \frac{e}{2\pi}\cdot\frac{1}{2\pi}\int d^2k\,F^s`,
1759 :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`
1760 (:math:`j^s_x` with :math:`s = \sigma_a/2`, in units of :math:`\hbar`).
1761 The unit is :math:`e/2\pi`: one :math:`e/\hbar` from the electric
1762 field, one :math:`\hbar/2` from the spin, and :math:`1/2\pi` from the
1763 Brillouin-zone integral, as for :math:`e^2/h = (e^2/\hbar)/2\pi`.
1764 When :math:`s_z` is conserved, the two spins decouple and
1765 :math:`j^s_x = \pm\frac{\hbar}{2}v_x` in each, so in a gap
1766 :math:`\sigma^s_{xy} = \frac{e}{2\pi}\cdot\frac{C_\uparrow - C_\downarrow}{2}`,
1767 :math:`C_\sigma` the Chern number of each spin sector (the sign of
1768 *hall_conductivity* and *chern_number*): :math:`\pm e/2\pi` for the
1769 Kane-Mele quantum spin Hall insulator (Kane and Mele 2005). Rashba
1770 coupling breaks the conservation of :math:`s_z`, and the value is
1771 then no longer quantized (the :math:`\mathbb{Z}_2` invariant, see
1772 *z2_invariant*, still is).
1774 Parameters
1775 ----------
1776 e_fermi
1777 See *hall_conductivity*.
1778 temperature
1779 See *hall_conductivity*.
1780 nk : Positive integer, or tuple of 2 integers, default 60
1781 plane
1782 See *hall_conductivity*.
1783 k_fixed : Real number or None, default None
1784 3D only: the
1785 plane (required in 3D, where only planes are computed).
1786 positions : Boolean, default True
1787 See *hall_conductivity*.
1788 spin_axis : 'x', 'y' or 'z', default 'z'
1789 Spin component
1790 :math:`a` carried by the current.
1791 refine
1792 See *hall_conductivity*.
1793 refine_fraction
1794 See *hall_conductivity*.
1796 Returns
1797 -------
1798 sigma_s
1799 Real number (or ndarray shaped like *e_fermi*), in
1800 units of :math:`e/2\pi`.
1801 '''
1802 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1803 refine_fraction, spin_axis)
1805 # ------------------------------------------------------------------
1806 # Berry-phase response beyond the Hall conductivity
1807 # ------------------------------------------------------------------
1809 def orbital_magnetization(
1810 self, e_fermi: float | ArrayLike, temperature: float = 0., nk: int | tuple[int, ...] = 60,
1811 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1812 refine: int = 1, refine_fraction: float = 0.05,
1813 ) -> float | NDArray[np.float64]:
1814 r'''
1815 Get the orbital magnetization in the modern theory (Thonhauser,
1816 Ceresoli, Vanderbilt and Resta, Phys. Rev. Lett. 95, 137205 (2005);
1817 Xiao, Shi and Niu, Phys. Rev. Lett. 95, 137204 (2005); Ceresoli et
1818 al., Phys. Rev. B 74, 024408 (2006)), at any Fermi level and
1819 temperature:
1821 .. math::
1823 M = \frac{e}{h}\,\frac{1}{2\pi}\int_{BZ} d^2k\sum_n\left[f_n\,m_n
1824 + g_n\,\Omega_n\right],\qquad
1825 m_n = \mathrm{Im}\langle\partial_xu_n|(H-E_n)|\partial_yu_n\rangle\, ,
1827 with :math:`g = k_BT\ln(1+e^{-(E-\mu)/k_BT})` (:math:`\mu - E` below
1828 :math:`\mu` and 0 above, at :math:`T = 0`). The first term is the
1829 self-rotation of the wavepackets, the second the circulation of
1830 the Berry-curvature currents at the edges. At :math:`T = 0` both
1831 combine into :math:`\sum_nf_n\,\mathrm{Im}\langle\partial_xu_n|(H+E_n-2\mu)|\partial_yu_n\rangle`,
1832 the formula of Thonhauser et al. It is evaluated, like
1833 *hall_conductivity*, from :math:`\langle n|\partial H|m\rangle`
1834 (no derivatives of the eigenstates), so it is finite where bands
1835 touch.
1837 **Sign and units.** The sign is that of *hall_conductivity*, so
1838 that the Streda formula reads :math:`\partial M/\partial\mu = \sigma`:
1839 in a gap (:math:`T = 0`), :math:`M` is linear in :math:`\mu` with
1840 slope the Chern number :math:`C`. In 2D, :math:`M` (a magnetic
1841 moment per area, i.e. a current) is in units of :math:`e/h` times
1842 the energy unit; with ``k_fixed=None`` in 3D, the vector
1843 :math:`(M_x, M_y, M_z)` per length. *positions* matters as for
1844 *hall_conductivity*, and the physical choice is the default.
1846 Parameters
1847 ----------
1848 e_fermi
1849 See *hall_conductivity*.
1850 temperature
1851 See *hall_conductivity*.
1852 nk
1853 See *hall_conductivity*.
1854 plane
1855 See *hall_conductivity*.
1856 k_fixed
1857 See *hall_conductivity*.
1858 positions
1859 See *hall_conductivity*.
1860 refine
1861 See *hall_conductivity*.
1862 refine_fraction
1863 See *hall_conductivity*.
1865 Returns
1866 -------
1867 M
1868 Real number (or ndarray shaped like *e_fermi*); in 3D
1869 with ``k_fixed=None``, ndarray of shape (3,) (or ``e_fermi.shape + (3,)``).
1871 Examples
1872 --------
1873 Example usage::
1875 mu = np.linspace(-0.5, 0.5, 11) # in the gap of a Chern insulator
1876 M = hal.orbital_magnetization(mu) # slope C = hall_conductivity(mu)
1877 '''
1878 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1879 refine_fraction, None, 'magnetization')
1881 def anomalous_nernst_conductivity(
1882 self, e_fermi: float | ArrayLike, temperature: float, nk: int | tuple[int, ...] = 60,
1883 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1884 refine: int = 1, refine_fraction: float = 0.05,
1885 ) -> float | NDArray[np.float64]:
1886 r'''
1887 Get the intrinsic anomalous Nernst (transverse thermoelectric)
1888 conductivity :math:`\alpha_{xy}`, :math:`j_x = \alpha_{xy}(-\partial_yT)`,
1889 from the Berry curvature weighted by the entropy of each state
1890 (Xiao, Yao, Fang and Niu, Phys. Rev. Lett. 97, 026603 (2006)):
1892 .. math::
1894 \alpha_{xy} = \frac{ek_B}{h}\,\frac{1}{2\pi}\int_{BZ}d^2k\sum_n s_n\,\Omega_n\, ,
1895 \qquad s = -f\ln f - (1-f)\ln(1-f)\, .
1897 Only the states within a few :math:`k_BT` of :math:`\mu` carry
1898 entropy, so it vanishes in a gap at low temperature. It obeys the
1899 Mott relation exactly,
1900 :math:`\alpha_{xy}(\mu, T) = \frac{1}{eT}\int dE\,(E-\mu)\left(-\frac{\partial f}{\partial E}\right)\sigma_{xy}(E)`,
1901 with :math:`\sigma_{xy}(E)` the :math:`T = 0` *hall_conductivity*
1902 at Fermi energy :math:`E`, which fixes the sign. The pair form of
1903 *hall_conductivity* keeps it finite where bands touch.
1905 **Units.** :math:`ek_B/h` in 2D (:math:`k_B = 1`, temperatures in
1906 energy units), per length for the 3D vector (``k_fixed=None``).
1908 Parameters
1909 ----------
1910 e_fermi
1911 See *hall_conductivity*.
1912 temperature : Positive real number
1913 In energy units (:math:`k_B = 1`).
1914 nk
1915 See *hall_conductivity*.
1916 plane
1917 See *hall_conductivity*.
1918 k_fixed
1919 See *hall_conductivity*.
1920 positions
1921 See *hall_conductivity*.
1922 refine
1923 See *hall_conductivity*.
1924 refine_fraction
1925 See *hall_conductivity*.
1927 Returns
1928 -------
1929 alpha : Real number (or ndarray, as *hall_conductivity*)
1931 Examples
1932 --------
1933 Example usage::
1935 alpha = hal.anomalous_nernst_conductivity(np.linspace(-3, 3, 301), 0.05)
1936 '''
1937 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1938 refine_fraction, None, 'nernst')
1940 def thermal_hall_conductivity(
1941 self, e_fermi: float | ArrayLike, temperature: float, nk: int | tuple[int, ...] = 60,
1942 plane: tuple[int, int] = (0, 1), k_fixed: float | None = None, positions: bool = True,
1943 refine: int = 1, refine_fraction: float = 0.05,
1944 ) -> float | NDArray[np.float64]:
1945 r'''
1946 Get the intrinsic thermal Hall conductivity of the electrons
1947 :math:`\kappa_{xy}`, :math:`j^Q_x = \kappa_{xy}(-\partial_yT)`
1948 (Qin, Niu and Shi, Phys. Rev. Lett. 107, 236601 (2011)):
1950 .. math::
1952 \kappa_{xy} = \frac{k_B^2T}{h}\,\frac{1}{2\pi}\int_{BZ}d^2k
1953 \sum_n c_2(x_n)\,\Omega_n\, ,\qquad
1954 c_2(x) = \int_x^\infty y^2\left(-\frac{\partial f}{\partial y}\right)dy\, ,
1956 :math:`x_n = (E_n-\mu)/k_BT` (:math:`c_2 = \pi^2/3` deep below
1957 :math:`\mu`, 0 far above). Equivalently, exactly,
1958 :math:`\kappa_{xy}(\mu, T) = \frac{1}{e^2T}\int dE\,(E-\mu)^2\left(-\frac{\partial f}{\partial E}\right)\sigma_{xy}(E)`
1959 with the :math:`T = 0` *hall_conductivity*. In a gap at low
1960 temperature it is the Wiedemann-Franz value
1961 :math:`\kappa_{xy} = \frac{\pi^2k_B^2}{3e^2}T\sigma_{xy}`: for a Chern
1962 insulator the quantized :math:`C\,\pi^2k_B^2T/3h`.
1964 **Units.** :math:`k_B^2/h` times the temperature (in energy units,
1965 :math:`k_B = 1`) in 2D, so that :math:`\kappa/T = \pi^2C/3` on a
1966 Chern plateau; per length for the 3D vector (``k_fixed=None``).
1967 The sign is that of *hall_conductivity*.
1969 Parameters
1970 ----------
1971 e_fermi
1972 See *hall_conductivity*.
1973 temperature : Positive real number
1974 In energy units (:math:`k_B = 1`).
1975 nk
1976 See *hall_conductivity*.
1977 plane
1978 See *hall_conductivity*.
1979 k_fixed
1980 See *hall_conductivity*.
1981 positions
1982 See *hall_conductivity*.
1983 refine
1984 See *hall_conductivity*.
1985 refine_fraction
1986 See *hall_conductivity*.
1988 Returns
1989 -------
1990 kappa : Real number (or ndarray, as *hall_conductivity*)
1992 Examples
1993 --------
1994 Example usage::
1996 kappa = hal.thermal_hall_conductivity(0., 0.01) # = pi^2/3 C T in the gap
1997 '''
1998 return self._hall(e_fermi, temperature, nk, plane, k_fixed, positions, refine,
1999 refine_fraction, None, 'thermal')
2001 def _axion_rate(
2002 self, ks: NDArray[np.float64], bands: list[int], params: dict, param: str,
2003 positions: bool,
2004 ) -> float:
2005 r'''
2006 Private method. The mesh average, over *ks*, of
2007 :math:`\mathrm{tr}[F_{\lambda x}F_{yz} + F_{\lambda y}F_{zx} + F_{\lambda z}F_{xy}]`
2008 for the bands *bands* at the parameters *params*, the derivative
2009 along *param* by central differences. The non-Abelian curvature is
2010 :math:`F_{\mu\nu} = i(X_\mu^\dagger X_\nu - X_\nu^\dagger X_\mu)`,
2011 :math:`(X_\mu)_{lm} = \langle l|\partial_\mu H|m\rangle/(E_m-E_l)`
2012 (:math:`m` in *bands*, :math:`l` not), i.e.
2013 :math:`Q|\partial_\mu u_m\rangle` in the eigenbasis: gauge covariant.
2014 '''
2015 value = params[param]
2016 step = 1e-5 * max(1., abs(value))
2017 others = [b for b in range(self.norb) if b not in bands]
2018 total, gap = 0., np.inf
2019 chunk = max(1, 400000 // self.norb ** 2)
2020 for c0 in range(0, len(ks), chunk):
2021 k = ks[c0:c0 + chunk]
2022 ham, dham = self._bloch_derivatives(k, list(np.eye(3)), positions, params)
2023 plus = self._bloch_derivatives(k, [], positions, dict(params, **{param: value + step}))[0]
2024 minus = self._bloch_derivatives(k, [], positions, dict(params, **{param: value - step}))[0]
2025 en, vec = np.linalg.eigh(ham)
2026 occ, emp = vec[:, :, bands], vec[:, :, others]
2027 de = en[:, others][:, :, None] - en[:, bands][:, None, :]
2028 gap = min(gap, np.abs(de).min())
2029 emp_h = emp.conj().transpose(0, 2, 1)
2030 x = [-(emp_h @ d @ occ) / de for d in [(plus - minus) / (2 * step), *dham]]
2031 f = {}
2032 for a, b in ((0, 1), (0, 2), (0, 3), (1, 2), (2, 3), (3, 1)):
2033 xa_h = x[a].conj().transpose(0, 2, 1)
2034 xb_h = x[b].conj().transpose(0, 2, 1)
2035 f[a, b] = 1j * (xa_h @ x[b] - xb_h @ x[a])
2036 trace = (np.einsum('kab,kba->k', f[0, 1], f[2, 3]) + np.einsum('kab,kba->k', f[0, 2], f[3, 1])
2037 + np.einsum('kab,kba->k', f[0, 3], f[1, 2]))
2038 total += trace.real.sum()
2039 error_handling.path_gap(gap, param, value)
2040 return total / len(ks)
2042 def axion_angle(
2043 self, bands: list[int], param: str, values: ArrayLike, nk: int | tuple[int, int, int] = 16,
2044 theta0: float = 0., positions: bool = True,
2045 ) -> NDArray[np.float64]:
2046 r'''
2047 Get the axion angle :math:`\theta` of a 3D insulator along a path of
2048 a parameter of its value functions (*set_hopping*), from the second
2049 Chern form (Qi, Hughes and Zhang, Phys. Rev. B 78, 195424 (2008);
2050 Essin, Moore and Vanderbilt, Phys. Rev. Lett. 102, 146805 (2009)):
2052 .. math::
2054 \frac{d\theta}{d\lambda} = \frac{1}{2\pi}\int_{BZ}d^3k\;
2055 \mathrm{tr}\left[F_{\lambda x}F_{yz} + F_{\lambda y}F_{zx} + F_{\lambda z}F_{xy}\right],
2057 :math:`F_{\mu\nu}` the non-Abelian Berry curvature of the bands
2058 *bands* in :math:`(\lambda, k_x, k_y, k_z)` space, so that
2059 :math:`\theta(\lambda) = \theta_0 + \int_{\lambda_0}^{\lambda}\frac{d\theta}{d\lambda'}d\lambda'`.
2060 This is the variation of the Chern-Simons axion coupling, computed
2061 from gauge-invariant quantities only (no smooth gauge is needed),
2062 and :math:`\theta_0` is its value at ``values[0]``: e.g. 0
2063 for a path starting from an atomic insulator. A closed path (a
2064 pumping cycle) changes :math:`\theta` by :math:`2\pi C_2`, with
2065 :math:`C_2` the second Chern number of the 4D family. At a point
2066 with time-reversal or inversion symmetry, :math:`\theta` is 0 or
2067 :math:`\pi` (mod :math:`2\pi`): :math:`\pi` for a strong topological
2068 insulator, cf. *z2_indices_3d*. The sign of :math:`\theta` follows
2069 the orientation :math:`(\lambda, k_x, k_y, k_z)` above.
2071 The path must stay gapped: the bands *bands* must be separated
2072 from the others at every k-point and every value (a ValueError
2073 names the value where they touch). The integrand is smooth, so the
2074 uniform k-mesh converges exponentially; the :math:`\lambda`
2075 integral is a trapezoid over *values* (error :math:`O(\Delta\lambda^2)`).
2076 :math:`\partial H/\partial\lambda` is a central difference, and
2077 :math:`\partial H/\partial\mathbf{k}` is exact (see
2078 *hall_conductivity*, whose *positions* convention applies).
2080 Parameters
2081 ----------
2082 bands : List of integers
2083 The occupied bands (e.g. the lower
2084 half of the spectrum).
2085 param : String
2086 Name of the parameter of the value functions
2087 that the path varies (see *set_params*; the other parameters keep
2088 their values in *params*).
2089 values : 1D array of at least two real numbers, increasing or decreasing
2090 The parameter values along the path.
2091 nk : Positive integer, or tuple of 3 positive integers, default 16
2092 k-mesh.
2093 theta0 : Real number, default 0
2094 :math:`\theta` at ``values[0]``.
2095 positions : Boolean, default True
2096 See *hall_conductivity*.
2098 Returns
2099 -------
2100 theta
2101 ndarray, shape ``(len(values),)``: :math:`\theta`
2102 at each value, not reduced modulo :math:`2\pi`.
2104 Examples
2105 --------
2106 Example usage::
2108 phi = np.linspace(0., 2 * np.pi, 41)
2109 theta = ks.axion_angle([0, 1], 'phi', phi) # theta[-1] = 2 pi C2
2110 '''
2111 self._check_static()
2112 error_handling.dim_exact(self.dim, 3)
2113 error_handling.hermitian_model(self.is_hermitian())
2114 error_handling.band_indices(bands, self.norb)
2115 error_handling.path_parameter(param)
2116 error_handling.path_values(values)
2117 if isinstance(nk, int):
2118 nk = (nk,) * 3
2119 error_handling.nk(nk, 3)
2120 error_handling.real_number(theta0, 'theta0')
2121 error_handling.boolean(positions, 'positions')
2122 values = np.asarray(values, dtype='f8')
2123 fracs = np.array(np.meshgrid(*[np.arange(n) / n for n in nk], indexing='ij')).reshape(3, -1).T
2124 ks = fracs @ self.rec_vec_k
2125 volume = abs(np.linalg.det(self.rec_vec_k))
2126 rates = [self._axion_rate(ks, bands, dict(self.params, **{param: float(v)}), param, positions)
2127 for v in values]
2128 return theta0 + cumulative_trapezoid(np.array(rates) * volume / (2 * PI), values, initial=0.)
2130 # ------------------------------------------------------------------
2131 # Finite samples and non-Hermitian band theory
2132 # ------------------------------------------------------------------
2134 def _hop_cells(self) -> list[tuple[int, int, tuple[int, ...], complex]]:
2135 r'''
2136 Private method. The stored hoppings with their lattice vector as
2137 integers :math:`(n_1, ...)` (from :math:`\mathbf{R}\cdot\mathbf{b}_i/2\pi`).
2138 '''
2139 rec = np.array(self.rec_vec, dtype='f8')
2140 return [(i, j, tuple(int(n) for n in np.rint(rec @ R / (2 * PI))), t)
2141 for i, j, R, t in self._hop]
2143 def _finite_entries(
2144 self, n_cells: int | tuple[int, ...], periodic: bool,
2145 ) -> tuple[tuple[int, ...], int, list[tuple[NDArray, NDArray, int, int, tuple, NDArray, complex]]]:
2146 '''
2147 Private method. For each stored hopping, the row and column indices
2148 it fills in a finite sample of *n_cells* cells (see *finite_ham*),
2149 with its orbitals, lattice vector (integers and Cartesian) and amplitude.
2150 '''
2151 if isinstance(n_cells, int):
2152 n_cells = (n_cells,) * self.dim
2153 n_tot = int(np.prod(n_cells))
2154 cells = np.array(np.meshgrid(*[np.arange(n) for n in n_cells], indexing='ij'))
2155 cells = cells.reshape(self.dim, -1).T # (n_tot, dim)
2156 strides = np.cumprod((1,) + tuple(n_cells[:-1]))
2157 index = cells @ strides
2158 order = np.argsort(index)
2159 cells = cells[order]
2160 dims = np.array(n_cells)
2161 entries = []
2162 for (i, j, R, t), (_, _, R_cart, _) in zip(self._hop_cells(), self._hop):
2163 target = cells + np.array(R)[None, :]
2164 if periodic:
2165 target = target % dims
2166 keep = np.ones(n_tot, bool)
2167 else:
2168 keep = np.all((target >= 0) & (target < dims), axis=1)
2169 rows = np.arange(n_tot)[keep] * self.norb + i
2170 cols = (target[keep] @ strides) * self.norb + j
2171 entries.append((rows, cols, i, j, R, R_cart, t))
2172 return n_cells, n_tot, entries
2174 def finite_ham(
2175 self, n_cells: int | tuple[int, ...], periodic: bool = False, sparse: bool = False,
2176 ) -> NDArray[np.complex128] | sp.csr_matrix:
2177 r'''
2178 Get the real-space Hamiltonian of a finite sample of
2179 :math:`N_1\times N_2\times\dots` unit cells of the model -- open
2180 boundaries (a flake, a finite chain), or periodic ones (a torus).
2181 Orbital *o* of cell :math:`(n_1, n_2, n_3)` is row
2182 :math:`((n_3 N_2 + n_2) N_1 + n_1)\,\mathrm{norb} + o`.
2184 Parameters
2185 ----------
2186 n_cells : Positive integer, or tuple of *dim* positive integers
2187 periodic : Boolean, default False
2188 Periodic boundary
2189 conditions (the spectrum is then that of *get_bands* on the
2190 :math:`N_1\times N_2\times\dots` k-mesh), else open ones.
2191 sparse : Boolean, default False
2192 Return a sparse CSR
2193 matrix (for large samples, e.g. with *tbkit.kpm*) instead of a
2194 dense array.
2196 Returns
2197 -------
2198 ham
2199 Complex ndarray (or CSR matrix), shape (N*norb, N*norb), N the number of cells.
2200 '''
2201 error_handling.nk(n_cells, self.dim)
2202 error_handling.boolean(periodic, 'periodic')
2203 error_handling.boolean(sparse, 'sparse')
2204 _, n_tot, entries = self._finite_entries(n_cells, periodic)
2205 onsite = np.diag(self.onsite) + self._onsite_offdiag
2206 if sparse:
2207 onsite_rows, onsite_cols = np.nonzero(onsite)
2208 rows = [c * self.norb + onsite_rows for c in range(n_tot)]
2209 cols = [c * self.norb + onsite_cols for c in range(n_tot)]
2210 vals = [onsite[onsite_rows, onsite_cols]] * n_tot
2211 for r, c, _, _, _, _, t in entries:
2212 rows.append(r)
2213 cols.append(c)
2214 vals.append(np.full(len(r), t, dtype='c16'))
2215 n = n_tot * self.norb
2216 return sp.coo_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
2217 shape=(n, n)).tocsr()
2218 ham = np.zeros((n_tot * self.norb, n_tot * self.norb), 'c16')
2219 for cell in range(n_tot):
2220 ham[cell*self.norb:(cell+1)*self.norb, cell*self.norb:(cell+1)*self.norb] += onsite
2221 for r, c, _, _, _, _, t in entries:
2222 np.add.at(ham, (r, c), t)
2223 return ham
2225 def finite_velocity(
2226 self, n_cells: int | tuple[int, ...], periodic: bool = False, sparse: bool = False,
2227 ) -> list[NDArray[np.complex128] | sp.csr_matrix]:
2228 r'''
2229 Get the velocity operators :math:`\mathbf{v} = i[H, \mathbf{r}]`
2230 (:math:`\hbar = 1`) of the finite sample of *finite_ham*, from the
2231 Cartesian bond vectors of the hoppings rather than from the site
2232 coordinates:
2234 .. math::
2236 (v_\alpha)_{ab} = i\,H_{ab}\,d_\alpha\, ,\qquad
2237 \mathbf{d} = \mathbf{R} + \boldsymbol\tau_j - \boldsymbol\tau_i\, ,
2239 so that they are well defined on a torus (``periodic=True``), where
2240 the position operator is not: a bond that wraps around the sample
2241 keeps its short displacement. They are the real-space counterpart of
2242 :math:`\partial H/\partial\mathbf{k}` with ``positions=True``
2243 (*hall_conductivity*), and the input of *tbkit.kpm.hall_conductivity*.
2245 Parameters
2246 ----------
2247 n_cells
2248 See *finite_ham*.
2249 periodic : Boolean, default False
2250 See *finite_ham*.
2251 sparse : Boolean, default False
2252 See *finite_ham*.
2254 Returns
2255 -------
2256 vel
2257 List of *space_dim* complex ndarrays (or CSR
2258 matrices), shape (N*norb, N*norb): :math:`v_x, v_y, (v_z)`.
2259 '''
2260 error_handling.nk(n_cells, self.dim)
2261 error_handling.boolean(periodic, 'periodic')
2262 error_handling.boolean(sparse, 'sparse')
2263 _, n_tot, entries = self._finite_entries(n_cells, periodic)
2264 tau = self.orbital_positions()
2265 n = n_tot * self.norb
2266 vel = []
2267 for alpha in range(self.space_dim):
2268 rows, cols, vals = [np.zeros(0, int)], [np.zeros(0, int)], [np.zeros(0, 'c16')]
2269 for r, c, i, j, _, R_cart, t in entries:
2270 d = R_cart[alpha] + tau[j, alpha] - tau[i, alpha]
2271 rows.append(r)
2272 cols.append(c)
2273 vals.append(np.full(len(r), 1j * t * d, dtype='c16'))
2274 mat = sp.coo_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
2275 shape=(n, n)).tocsr()
2276 vel.append(mat if sparse else mat.toarray())
2277 return vel
2279 def get_ham_beta(self, beta: complex) -> NDArray[np.complex128]:
2280 r'''
2281 Get the non-Bloch Hamiltonian of a 1D model,
2282 :math:`H(\beta) = \sum_{n} t(n)\, \beta^{n}`: the Bloch Hamiltonian
2283 continued off the unit circle, :math:`\beta = e^{ik}` with complex
2284 :math:`k`. On the generalized Brillouin zone (see *gbz*), its
2285 spectrum is the open-chain spectrum of a non-Hermitian chain.
2287 Parameters
2288 ----------
2289 beta : Complex number, nonzero
2291 Returns
2292 -------
2293 ham : Complex ndarray, shape (norb, norb)
2294 '''
2295 error_handling.dim_exact(self.dim, 1)
2296 error_handling.nonzero_number(beta, 'beta')
2297 ham = np.diag(self.onsite).astype('c16') + self._onsite_offdiag
2298 for i, j, (n,), t in self._hop_cells():
2299 ham[i, j] += t * complex(beta) ** n
2300 return ham
2302 def spectral_winding(self, e_ref: complex = 0., nk: int = 400) -> float:
2303 r'''
2304 Get the spectral winding number of a 1D (non-Hermitian) model about
2305 the reference energy :math:`E_r`,
2307 .. math::
2309 W(E_r) = \frac{1}{2\pi i}\oint_0^{2\pi} dk\,
2310 \partial_k \ln\det\left[H(k) - E_r\right]\, ,
2312 the number of times the periodic-chain spectrum winds around
2313 :math:`E_r` in the complex plane. :math:`W \neq 0` for some
2314 :math:`E_r` is the criterion for the non-Hermitian skin effect: the
2315 open-chain eigenstates then pile up at one end (Okuma, Kawabata,
2316 Shiozaki and Sato, Phys. Rev. Lett. 124, 086801 (2020)).
2318 Parameters
2319 ----------
2320 e_ref : Complex number, default 0
2321 Reference energy.
2322 nk : Positive integer, default 400
2323 k-points.
2325 Returns
2326 -------
2327 W
2328 Real number, close to an integer (it is ill defined if
2329 :math:`E_r` lies on the periodic spectrum).
2330 '''
2331 error_handling.dim_exact(self.dim, 1)
2332 error_handling.number(e_ref, 'e_ref')
2333 error_handling.positive_int(nk, 'nk')
2334 ks = np.arange(nk + 1) / nk * self.rec_vec_k[0, 0]
2335 dets = np.array([np.linalg.det(self.get_ham([k]) - e_ref * np.eye(self.norb)) for k in ks])
2336 steps = np.angle(dets[1:] / dets[:-1])
2337 return float(np.sum(steps) / (2 * PI))
2339 def gbz(self, energies: ArrayLike) -> NDArray[np.complex128]:
2340 r'''
2341 Get the generalized Brillouin zone of a 1D non-Hermitian model
2342 (Yao and Wang, Phys. Rev. Lett. 121, 086803 (2018); Yokomizo and
2343 Murakami, Phys. Rev. Lett. 123, 066404 (2019)): for each energy
2344 :math:`E`, the roots of :math:`\det[H(\beta)-E] = 0` sorted by modulus,
2345 :math:`|\beta_1|\le\dots\le|\beta_{p+s}|` (:math:`p` the order of its
2346 pole at :math:`\beta = 0`), and the middle pair
2347 :math:`\beta_p, \beta_{p+1}`. The open-chain spectrum, in the limit of
2348 a long chain, is where :math:`|\beta_p| = |\beta_{p+1}|`; those
2349 :math:`\beta` trace the generalized Brillouin zone, on which the
2350 non-Bloch bands :math:`H(\beta)` (see *get_ham_beta*) replace the
2351 Bloch ones. For a Hermitian model it is the unit circle.
2353 Parameters
2354 ----------
2355 energies
2356 Complex array of energies (e.g. open-chain eigenvalues
2357 from *finite_ham*).
2359 Returns
2360 -------
2361 beta : Complex ndarray, shape (len(energies), 2): :math:`\beta_p, \beta_{p+1}`
2362 '''
2363 error_handling.dim_exact(self.dim, 1)
2364 energies = np.atleast_1d(np.asarray(energies, dtype='c16'))
2365 error_handling.ndarray_empty(energies, 'energies')
2366 powers = [n for _, _, (n,), _ in self._hop_cells()] + [0]
2367 low, high = self.norb * min(powers), self.norb * max(powers)
2368 n_samples = 2 * (high - low) + 8
2369 theta = 2 * PI * np.arange(n_samples) / n_samples
2370 hams = np.array([self.get_ham_beta(np.exp(1j * th)) for th in theta])
2371 out = np.zeros((len(energies), 2), 'c16')
2372 exps = np.arange(low, high + 1)
2373 for n, e in enumerate(energies):
2374 dets = np.linalg.det(hams - e * np.eye(self.norb)[None])
2375 # Laurent coefficients c_m of det[H(beta) - E] = sum_m c_m beta^m
2376 coef = np.array([np.mean(dets * np.exp(-1j * m * theta)) for m in exps])
2377 nonzero = np.abs(coef) > 1e-12 * np.max(np.abs(coef))
2378 first, last = np.argmax(nonzero), len(coef) - 1 - np.argmax(nonzero[::-1])
2379 p = -exps[first]
2380 roots = np.roots(coef[first:last + 1][::-1])
2381 roots = roots[np.argsort(np.abs(roots))]
2382 error_handling.gbz_roots(len(roots), p)
2383 out[n] = roots[p - 1], roots[p]
2384 return out
2386 # ------------------------------------------------------------------
2387 # Berry phases, Wannier centres, Z2, symmetries, quantum geometry
2388 # ------------------------------------------------------------------
2390 def orbital_positions(self) -> NDArray[np.float64]:
2391 r'''
2392 Get the positions :math:`\boldsymbol\tau` of the orbitals within the
2393 unit cell (each site twice, for its two spin components, if
2394 ``spin=True``).
2396 Returns
2397 -------
2398 tau : Real ndarray, shape (norb, space_dim)
2399 '''
2400 tau = np.array([dic['r0'] for dic in self.lat.unit_cell], dtype='f8')
2401 return np.repeat(tau, 2, axis=0) if self.spin else tau
2403 def _loop(
2404 self, bands: list[int], nk: int, direction: int, fracs_perp: dict[int, float],
2405 positions: bool,
2406 ) -> NDArray[np.complex128]:
2407 '''
2408 Private method. Wilson loop matrix of *bands* along the closed path
2409 k0 -> k0 + b_direction (nk steps), at fixed fractional coordinates
2410 *fracs_perp* along the other reciprocal vectors.
2411 '''
2412 frac = np.zeros(self.dim)
2413 for d, f in fracs_perp.items():
2414 frac[d] = f
2415 b = self.rec_vec_k[direction]
2416 k0 = frac @ self.rec_vec_k
2417 # each link carries exp(-i dk.tau): the Bloch sums then include the
2418 # orbital positions, so that the phases locate the Wannier centres
2419 # in space rather than on the cell origin.
2420 dk_space = self.k_basis @ (b / nk)
2421 phase = np.exp(-1j * self.orbital_positions() @ dk_space) if positions \
2422 else np.ones(self.norb)
2423 vs = self._subspaces(k0 + np.arange(nk)[:, None] * b[None] / nk, bands)
2424 links = vs.conj().transpose(0, 2, 1) @ (phase[:, None] * np.roll(vs, -1, axis=0))
2425 # keep only the unitary part of each link (its singular values
2426 # tend to 1 as nk grows): the product then stays unitary
2427 u, _, vh = np.linalg.svd(links)
2428 wilson = np.eye(len(bands), dtype='c16')
2429 for unitary in u @ vh:
2430 wilson = wilson @ unitary
2431 return wilson
2433 def _perp(self, direction: int, k_perp: float | tuple | None) -> dict[int, float]:
2434 '''
2435 Private method. Fractional coordinates of a loop along the other
2436 reciprocal vectors, from *k_perp* (a number in 2D, a pair in 3D).
2437 '''
2438 others = [d for d in range(self.dim) if d != direction]
2439 if k_perp is None:
2440 k_perp = (0.,) * len(others)
2441 if not isinstance(k_perp, tuple):
2442 k_perp = (k_perp,)
2443 error_handling.k_perp(k_perp, len(others))
2444 return dict(zip(others, k_perp))
2446 def berry_phase(
2447 self, bands: int | list[int], nk: int = 100, direction: int = 0,
2448 k_perp: float | tuple[float, float] | None = None, positions: bool = True,
2449 ) -> float:
2450 r'''
2451 Get the Berry phase of a group of bands along a closed loop across the
2452 Brillouin zone, :math:`\mathbf{k}\to\mathbf{k}+\mathbf{b}_{direction}`,
2454 .. math::
2456 \gamma = \oint \mathbf{A}\cdot d\mathbf{k} = -\arg\det W\, ,\qquad
2457 W = \prod_{m} \langle u(\mathbf{k}_m)|u(\mathbf{k}_{m+1})\rangle\, ,
2459 with :math:`\mathbf{A} = i\langle u|\nabla_{\mathbf{k}}u\rangle` (the
2460 convention of *berry_curvature*). In 1D it is the Zak phase: its
2461 value :math:`\gamma/2\pi`, modulo 1, is the centre of the band's
2462 Wannier functions in units of :math:`\mathbf{a}_{direction}` (the
2463 polarization).
2465 Parameters
2466 ----------
2467 bands : Band index, or list of band indices
2468 nk : Positive integer, default 100
2469 Number of k-points on the loop.
2470 direction : Integer, default 0
2471 Reciprocal vector the loop follows.
2472 k_perp : Real number (2D) or pair of real numbers (3D)
2473 Default
2474 value None (zero). Fractional coordinates of the loop along the
2475 other reciprocal vectors.
2476 positions : Boolean, default True
2477 Include the orbital
2478 positions in the Bloch functions (the physical choice: the phase
2479 then locates the Wannier centres in the crystal). If False, all
2480 orbitals are placed on the cell origin (the periodic gauge of
2481 *get_ham*); the two differ by :math:`\sum_n\mathbf{b}\cdot\boldsymbol\tau`
2482 weighted by the orbital content, and agree on quantized
2483 differences between phases.
2485 Returns
2486 -------
2487 gamma : Real number in :math:`(-\pi, \pi]`
2488 '''
2489 bands = self._check_loop(bands, nk, direction, positions)
2490 wilson = self._loop(bands, nk, direction, self._perp(direction, k_perp), positions)
2491 return float(-np.angle(np.linalg.det(wilson)))
2493 def wannier_centers(
2494 self, bands: int | list[int], nk: int = 100, direction: int = 0,
2495 k_perp: float | tuple[float, float] | None = None, positions: bool = True,
2496 ) -> NDArray[np.float64]:
2497 r'''
2498 Get the (hybrid) Wannier centres of a group of bands: minus the
2499 phases of the eigenvalues of the Wilson loop (see *berry_phase*),
2500 divided by :math:`2\pi`. They are the positions, in units of
2501 :math:`\mathbf{a}_{direction}` and modulo 1, of the maximally
2502 localized Wannier functions along that direction (hybrid: still
2503 Bloch waves along the others, at *k_perp*). Their sum is
2504 :math:`\gamma/2\pi`, modulo 1.
2506 Parameters
2507 ----------
2508 bands : Band index, or list of band indices
2509 nk : Positive integer, default 100
2510 Number of k-points on the loop.
2511 direction : Integer, default 0
2512 k_perp
2513 See *berry_phase*.
2514 positions : Boolean, default True
2515 See *berry_phase*.
2517 Returns
2518 -------
2519 centers : Real ndarray, sorted, in :math:`[0, 1)`
2520 '''
2521 bands = self._check_loop(bands, nk, direction, positions)
2522 wilson = self._loop(bands, nk, direction, self._perp(direction, k_perp), positions)
2523 return np.sort((-np.angle(np.linalg.eigvals(wilson)) / (2*PI)) % 1.)
2525 def wannier_flow(
2526 self, bands: int | list[int], nk: int = 100, nk_perp: int = 51, direction: int = 0,
2527 flow: int | None = None, k_fixed: float = 0., positions: bool = True,
2528 k_range: tuple[float, float] = (0., 1.),
2529 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
2530 r'''
2531 Get the flow of the hybrid Wannier centres (see *wannier_centers*)
2532 along *direction*, as the loop is moved across the Brillouin zone
2533 along the reciprocal vector *flow*. In a Chern insulator the centres
2534 wind by :math:`C` cells over a full period; in a
2535 :math:`\mathbb{Z}_2` topological insulator, the Kramers pairs swap
2536 partners over half a period (see *z2_invariant*).
2538 Parameters
2539 ----------
2540 bands : Band index, or list of band indices
2541 nk : Positive integer, default 100
2542 k-points on each loop.
2543 nk_perp : Positive integer, default 51
2544 Number of loops.
2545 direction : Integer, default 0
2546 Direction of the loops.
2547 flow : Integer
2548 Default value None (the first other direction).
2549 Reciprocal vector along which the loops are moved.
2550 k_fixed : Real number, default 0
2551 3D only: fractional
2552 coordinate along the remaining reciprocal vector.
2553 positions : Boolean, default True
2554 See *berry_phase*.
2555 k_range : Pair of real numbers
2556 Default value (0, 1). Range of
2557 the fractional coordinate along *flow* (both ends included).
2559 Returns
2560 -------
2561 fracs : Real ndarray, shape (nk_perp,)
2562 Fractional coordinates along *flow*.
2563 centers : Real ndarray, shape (nk_perp, len(bands))
2564 Sorted centres, in [0, 1).
2565 '''
2566 bands = self._check_loop(bands, nk, direction, positions)
2567 error_handling.dim_min(self.dim, 2)
2568 error_handling.positive_int(nk_perp, 'nk_perp')
2569 others = [d for d in range(self.dim) if d != direction]
2570 if flow is None:
2571 flow = others[0]
2572 error_handling.flow(flow, direction, self.dim)
2573 error_handling.real_number(k_fixed, 'k_fixed')
2574 error_handling.lims(k_range)
2575 fracs = np.linspace(k_range[0], k_range[1], nk_perp)
2576 centers = np.zeros((nk_perp, len(bands)))
2577 for n, f in enumerate(fracs):
2578 perp = {d: (f if d == flow else k_fixed) for d in others}
2579 wilson = self._loop(bands, nk, direction, perp, positions)
2580 centers[n] = np.sort((-np.angle(np.linalg.eigvals(wilson)) / (2*PI)) % 1.)
2581 return fracs, centers
2583 def z2_invariant(
2584 self, bands: list[int], nk: int = 100, nk_perp: int = 51, direction: int = 0,
2585 flow: int | None = None, k_fixed: float = 0.,
2586 ) -> int:
2587 r'''
2588 Get the :math:`\mathbb{Z}_2` invariant of a time-reversal-symmetric
2589 group of bands (an even number, e.g. all occupied bands of the
2590 Kane-Mele model) from the flow of its hybrid Wannier centres over
2591 half the Brillouin zone (Soluyanov and Vanderbilt, Phys. Rev. B 83,
2592 235401 (2011); Yu, Qi, Bernevig, Fang and Dai, Phys. Rev. B 84,
2593 075119 (2011)): follow the midpoint of the largest gap between the
2594 centres, and count, modulo 2, how many centres it jumps over. It is
2595 odd when the Kramers partners exchange along the flow -- a
2596 topological insulator -- and even otherwise.
2598 Unlike the parity criterion (*parity_z2*), it needs no inversion
2599 symmetry, so it also works with Rashba coupling.
2601 Parameters
2602 ----------
2603 bands : List of band indices (an even number of bands)
2604 nk : Positive integer, default 100
2605 k-points on each loop.
2606 nk_perp : Positive integer, default 51
2607 Number of loops
2608 over the half zone (the flow must be resolved finely enough that
2609 no two centres cross between neighbouring loops unnoticed).
2610 direction : Integer, default 0
2611 See *wannier_flow*.
2612 flow : Integer, default None
2613 See *wannier_flow*.
2614 k_fixed : Real number, default 0
2615 3D only: a time-reversal
2616 invariant plane has *k_fixed* 0 or 0.5 (the weak and strong
2617 indices follow from the six such planes).
2619 Returns
2620 -------
2621 nu : 0 (trivial) or 1 (topological)
2622 '''
2623 error_handling.even_bands(bands)
2624 _, centers = self.wannier_flow(bands, nk, nk_perp, direction, flow, k_fixed,
2625 positions=False, k_range=(0., 0.5))
2626 gap = [self._largest_gap(c) for c in centers]
2627 jumps = sum(self._between(gap[n], gap[n + 1], centers[n + 1])
2628 for n in range(len(gap) - 1))
2629 return int(jumps % 2)
2631 def z2_indices_3d(
2632 self, bands: list[int], nk: int = 100, nk_perp: int = 51,
2633 ) -> tuple[int, int, int, int]:
2634 r'''
2635 Get the four :math:`\mathbb{Z}_2` indices
2636 :math:`(\nu_0;\nu_1\nu_2\nu_3)` of a 3D time-reversal-invariant
2637 insulator (Fu, Kane and Mele, Phys. Rev. Lett. 98, 106803 (2007);
2638 Moore and Balents, Phys. Rev. B 75, 121306(R) (2007)) from the
2639 :math:`\mathbb{Z}_2` invariants of its six time-reversal-invariant
2640 planes :math:`k_i = 0` and :math:`k_i = 1/2` (fractional coordinates
2641 along :math:`\mathbf{b}_i`), each computed by *z2_invariant*:
2643 .. math::
2645 \nu_0 = \nu(k_i = 0) + \nu(k_i = 1/2) \bmod 2\, ,\qquad
2646 \nu_i = \nu(k_i = 1/2)\, .
2648 The strong index :math:`\nu_0` comes out the same for the three
2649 :math:`i`, which is checked (a ValueError asks for a finer flow
2650 otherwise). :math:`\nu_0 = 1` is a strong topological insulator; a
2651 weak one has :math:`\nu_0 = 0` and some :math:`\nu_i = 1`, i.e.
2652 stacked quantum spin Hall layers normal to
2653 :math:`\nu_1\mathbf{b}_1 + \nu_2\mathbf{b}_2 + \nu_3\mathbf{b}_3`.
2654 No inversion symmetry is needed (compare *parity_z2*).
2656 Parameters
2657 ----------
2658 bands : List of band indices (an even number of bands)
2659 nk : Positive integer, default 100
2660 See *z2_invariant*.
2661 nk_perp : Positive integer, default 51
2662 See *z2_invariant*.
2664 Returns
2665 -------
2666 indices : Tuple of four integers, 0 or 1: :math:`(\nu_0, \nu_1, \nu_2, \nu_3)`
2667 '''
2668 error_handling.dim_exact(self.dim, 3)
2669 nu = np.zeros((3, 2), int)
2670 for i in range(3):
2671 direction, flow = [d for d in range(3) if d != i]
2672 for n, k_fixed in enumerate((0., 0.5)):
2673 nu[i, n] = self.z2_invariant(bands, nk, nk_perp, direction, flow, k_fixed)
2674 strong = nu.sum(axis=1) % 2
2675 error_handling.strong_index(strong)
2676 return int(strong[0]), int(nu[0, 1]), int(nu[1, 1]), int(nu[2, 1])
2678 def entanglement_spectrum(
2679 self, ks: ArrayLike, region: list[int], bands: int | list[int],
2680 ) -> NDArray[np.float64]:
2681 r'''
2682 Get the k-resolved entanglement spectrum of a group of bands
2683 (Peschel, J. Phys. A 36, L205 (2003); Li and Haldane, Phys. Rev. Lett.
2684 101, 010504 (2008); Fidkowski, Phys. Rev. Lett. 104, 130502 (2010)):
2685 the eigenvalues :math:`\xi_n(\mathbf{k})` of the correlation matrix
2686 :math:`P(\mathbf{k})` of the filled *bands*, restricted to the
2687 orbitals of *region*,
2689 .. math::
2691 C_A(\mathbf{k}) = \left[P(\mathbf{k})\right]_{A}\, ,\qquad
2692 P(\mathbf{k}) = \sum_{n\in\mathrm{bands}}|u_n(\mathbf{k})\rangle\langle u_n(\mathbf{k})|\, ,
2694 between 0 and 1. For a ribbon or a supercell
2695 (*tbkit.moire.supercell*) cut in two halves across its width, the
2696 states of the whole are either inside the region
2697 (:math:`\xi \approx 1`) or outside (:math:`\xi \approx 0`), except
2698 for those localized at the cuts: a topological phase leaves
2699 entanglement modes inside :math:`(0, 1)` -- a mode pinned at 1/2 in
2700 a chiral-symmetric chain, a branch that flows from 0 to 1 across
2701 the zone in a Chern insulator. The single-particle entanglement
2702 energies are :math:`\epsilon = \ln[(1-\xi)/\xi]`.
2704 Parameters
2705 ----------
2706 ks : ndarray, shape (nk, dim)
2707 k-points (see *get_bands*).
2708 region : List of orbital indices (the region :math:`A`)
2709 bands : Band index, or list of band indices (the filled bands)
2711 Returns
2712 -------
2713 xi : Real ndarray, shape (nk, len(region)), sorted ascending at every k
2714 '''
2715 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
2716 error_handling.ks(ks, self.dim)
2717 error_handling.region(region, self.norb)
2718 if isinstance(bands, int):
2719 bands = [bands]
2720 error_handling.band_indices(bands, self.norb)
2721 error_handling.hermitian_kspace(self.is_hermitian(), self._overlap_hop)
2722 v = self._subspaces(ks, bands)[:, region, :]
2723 return np.linalg.eigvalsh(v @ v.conj().transpose(0, 2, 1))
2725 @staticmethod
2726 def _largest_gap(centers: NDArray[np.float64]) -> float:
2727 '''
2728 Private method. Midpoint of the largest gap between sorted centres on the circle [0, 1).
2729 '''
2730 ext = np.append(centers, centers[0] + 1.)
2731 n = np.argmax(np.diff(ext))
2732 return float((ext[n] + ext[n + 1]) / 2 % 1.)
2734 @staticmethod
2735 def _between(a: float, b: float, centers: NDArray[np.float64]) -> int:
2736 '''
2737 Private method. Number of centres on the shorter arc between a and b.
2738 '''
2739 d = (b - a) % 1.
2740 start, length = (a, d) if d <= 0.5 else (b, 1. - d)
2741 return int(np.sum((centers - start) % 1. < length))
2743 def _check_loop(self, bands, nk, direction, positions) -> list[int]:
2744 '''
2745 Private method. Validate the arguments shared by the Wilson-loop methods.
2746 '''
2747 if isinstance(bands, int):
2748 bands = [bands]
2749 error_handling.band_indices(bands, self.norb)
2750 error_handling.positive_int(nk, 'nk')
2751 error_handling.direction(direction, self.dim)
2752 error_handling.boolean(positions, 'positions')
2753 return bands
2755 def _trims(self) -> list[NDArray[np.float64]]:
2756 '''
2757 Private method. Time-reversal-invariant momenta: fractional
2758 coordinates in {0, 1/2}^dim, as k points.
2759 '''
2760 grid = np.array(np.meshgrid(*[[0., 0.5]] * self.dim, indexing='ij')).reshape(self.dim, -1).T
2761 return [f @ self.rec_vec_k for f in grid]
2763 def _operator(self, op, k: NDArray[np.float64]) -> NDArray[np.complex128]:
2764 '''
2765 Private method. Matrix of a symmetry operator at k (op is a matrix, or a callable of k).
2766 '''
2767 mat = np.asarray(op(k) if callable(op) else op, dtype='c16')
2768 error_handling.operator(mat, self.norb)
2769 return mat
2771 def parity_z2(self, inversion, bands: list[int]) -> int:
2772 r'''
2773 Get the :math:`\mathbb{Z}_2` invariant of an inversion- and
2774 time-reversal-symmetric model from the parity eigenvalues of its
2775 occupied bands at the time-reversal-invariant momenta
2776 :math:`\Gamma_i` (Fu and Kane, Phys. Rev. B 76, 045302 (2007)):
2778 .. math::
2780 (-1)^\nu = \prod_i \delta_i\, ,\qquad
2781 \delta_i = \prod_{m=1}^{N} \xi_{2m}(\Gamma_i)\, ,
2783 one parity :math:`\xi = \pm1` per Kramers pair. In 3D, :math:`\nu`
2784 is the strong index.
2786 Parameters
2787 ----------
2788 inversion
2789 Complex ndarray, shape (norb, norb), or callable of
2790 k returning one: the inversion operator :math:`P`, with
2791 :math:`P H(\mathbf{k}) P^{-1} = H(-\mathbf{k})` (e.g. swapping the
2792 two sublattices of graphene, and the identity on spin).
2793 bands : List of band indices (the occupied Kramers pairs)
2795 Returns
2796 -------
2797 nu : 0 (trivial) or 1 (topological)
2798 '''
2799 error_handling.even_bands(bands)
2800 error_handling.band_indices(bands, self.norb)
2801 product = 1
2802 for k in self._trims():
2803 v = self._subspace(k, bands)
2804 xi = np.linalg.eigvals(v.conj().T @ self._operator(inversion, k) @ v)
2805 error_handling.parities(xi)
2806 n_minus = int(np.sum(xi.real < 0))
2807 product *= (-1) ** (n_minus // 2)
2808 return 0 if product == 1 else 1
2810 def symmetry_error(
2811 self, op, k_map: str = 'minus', antiunitary: bool = False, anti: bool = False,
2812 nk: int = 5, seed: int = 0, transpose: bool = False,
2813 ) -> float:
2814 r'''
2815 Check a symmetry of the Bloch Hamiltonian,
2817 .. math::
2819 H(g\mathbf{k}) = \pm\, U\, \mathcal{K}H(\mathbf{k})\mathcal{K}^{-1} U^\dagger\, ,
2821 with :math:`g\mathbf{k} = -\mathbf{k}` or :math:`\mathbf{k}`,
2822 :math:`\mathcal{K}` complex conjugation if *antiunitary*, and the
2823 minus sign if *anti*: time reversal is antiunitary with
2824 :math:`\mathbf{k}\to-\mathbf{k}` and the plus sign, particle-hole
2825 antiunitary with :math:`\mathbf{k}\to-\mathbf{k}` and the minus sign,
2826 chiral unitary with :math:`\mathbf{k}\to\mathbf{k}` and the minus sign,
2827 inversion unitary with :math:`\mathbf{k}\to-\mathbf{k}` and the plus
2828 sign. With *transpose*, :math:`H^T` replaces :math:`\mathcal{K}H\mathcal{K}^{-1} = H^*`
2829 (and :math:`H^\dagger` with *antiunitary* too): the symmetries
2830 of non-Hermitian Hamiltonians, where :math:`H^T \neq H^*` (see
2831 *non_hermitian_class*).
2833 Parameters
2834 ----------
2835 op : Complex ndarray, shape (norb, norb), or callable of k: :math:`U`
2836 k_map : 'minus' or 'identity', default 'minus'
2837 antiunitary : Boolean, default False
2838 anti : Boolean, default False
2839 nk : Positive integer, default 5
2840 Points per direction of
2841 the (randomly shifted) test mesh.
2842 seed : Integer, default 0
2843 Seed of the random shift.
2844 transpose : Boolean, default False
2846 Returns
2847 -------
2848 error : Real number
2849 Largest deviation over the mesh (0 for a symmetry).
2850 '''
2851 error_handling.k_map(k_map)
2852 error_handling.boolean(antiunitary, 'antiunitary')
2853 error_handling.boolean(anti, 'anti')
2854 error_handling.positive_int(nk, 'nk')
2855 error_handling.boolean(transpose, 'transpose')
2856 rng = np.random.default_rng(seed)
2857 fracs = np.array(np.meshgrid(*[np.arange(nk) / nk] * self.dim, indexing='ij')).reshape(self.dim, -1).T
2858 fracs = fracs + rng.uniform(0., 1. / nk, self.dim)
2859 sign = -1. if anti else 1.
2860 err = 0.
2861 for f in fracs:
2862 k = f @ self.rec_vec_k
2863 gk = -k if k_map == 'minus' else k
2864 ham = self.get_ham(k)
2865 u = self._operator(op, k)
2866 ham = ham.conj() if antiunitary else ham
2867 rhs = u @ (ham.T if transpose else ham) @ u.conj().T
2868 err = max(err, float(np.max(np.abs(self.get_ham(gk) - sign * rhs))))
2869 return err
2871 def tenfold_class(
2872 self, time_reversal=None, particle_hole=None, chiral=None, nk: int = 5, tol: float = 1e-8,
2873 ) -> str:
2874 r'''
2875 Get the Altland-Zirnbauer symmetry class (the "tenfold way") of the
2876 model from its time-reversal :math:`T = U_T\mathcal{K}`, particle-hole
2877 :math:`C = U_C\mathcal{K}` and chiral :math:`S = U_S` symmetries. Each
2878 given operator is first checked (see *symmetry_error*) to be a
2879 symmetry; the class then follows from :math:`T^2 = U_TU_T^* = \pm1`,
2880 :math:`C^2 = \pm1`, and the presence of :math:`S`. Any two of the
2881 three imply the third (:math:`S = TC`), which need not be given.
2883 Parameters
2884 ----------
2885 time_reversal : Matrix :math:`U_T`, or None (absent)
2886 particle_hole : Matrix :math:`U_C`, or None (absent)
2887 chiral : Matrix :math:`U_S`, or None (absent)
2888 nk : Positive integer, default 5
2889 See *symmetry_error*.
2890 tol : Positive real, default 1e-8
2891 Tolerance of the checks.
2893 Returns
2894 -------
2895 name : String: 'A', 'AIII', 'AI', 'BDI', 'D', 'DIII', 'AII', 'CII', 'C' or 'CI'
2896 '''
2897 error_handling.positive_real(tol, 'tol')
2898 k0 = np.zeros(self.dim)
2899 for op, k_map, antiu, anti in ((time_reversal, 'minus', True, False),
2900 (particle_hole, 'minus', True, True),
2901 (chiral, 'identity', False, True)):
2902 if op is not None:
2903 error_handling.is_symmetry(self.symmetry_error(op, k_map, antiu, anti, nk), tol)
2904 u_t = None if time_reversal is None else self._operator(time_reversal, k0)
2905 u_c = None if particle_hole is None else self._operator(particle_hole, k0)
2906 # two of the three imply the third: C = S T, T = S C
2907 if chiral is not None and (u_t is None) != (u_c is None):
2908 u_s = self._operator(chiral, k0)
2909 if u_t is not None:
2910 u_c = u_s @ u_t
2911 elif u_c is not None:
2912 u_t = u_s @ u_c
2913 squares = []
2914 for u in (u_t, u_c):
2915 if u is None:
2916 squares.append(0)
2917 continue
2918 square = u @ u.conj()
2919 error_handling.square(square, tol)
2920 squares.append(int(np.sign(square[0, 0].real)))
2921 has_s = chiral is not None or (u_t is not None and u_c is not None)
2922 table = {(0, 0, False): 'A', (0, 0, True): 'AIII', (1, 0, False): 'AI',
2923 (1, 1, True): 'BDI', (0, 1, False): 'D', (-1, 1, True): 'DIII',
2924 (-1, 0, False): 'AII', (-1, -1, True): 'CII', (0, -1, False): 'C',
2925 (1, -1, True): 'CI'}
2926 return table[(squares[0], squares[1], has_s)]
2928 def non_hermitian_class(
2929 self, time_reversal=None, particle_hole=None, chiral=None, time_reversal_dagger=None,
2930 particle_hole_dagger=None, sublattice=None, pseudo_hermiticity=None, canonical: bool = False,
2931 nk: int = 5, tol: float = 1e-8,
2932 ) -> str:
2933 r'''
2934 Get the non-Hermitian symmetry class (the "38-fold way", Kawabata,
2935 Shiozaki, Ueda and Sato, Phys. Rev. X 9, 041015 (2019); Zhou and Lee,
2936 Phys. Rev. B 99, 235112 (2019)) of the model. A non-Hermitian
2937 Hamiltonian has :math:`H^* \neq H^T`, so each symmetry of the tenfold
2938 way (see *tenfold_class*) splits in two:
2940 .. math::
2942 \text{TRS: } U_T H^*(\mathbf{k}) U_T^\dagger = H(-\mathbf{k}),\qquad
2943 \text{TRS}^\dagger\text{: } U_{C_+} H^T(\mathbf{k}) U_{C_+}^\dagger = H(-\mathbf{k}),
2945 \text{PHS: } U_C H^T(\mathbf{k}) U_C^\dagger = -H(-\mathbf{k}),\qquad
2946 \text{PHS}^\dagger\text{: } U_{T_-} H^*(\mathbf{k}) U_{T_-}^\dagger = -H(-\mathbf{k}),
2948 \text{CS: } U_\Gamma H^\dagger(\mathbf{k}) U_\Gamma^\dagger = -H(\mathbf{k}),\qquad
2949 \text{SLS: } U_S H(\mathbf{k}) U_S^\dagger = -H(\mathbf{k}),\qquad
2950 \text{pH: } U_\eta H^\dagger(\mathbf{k}) U_\eta^\dagger = H(\mathbf{k}).
2952 Each given operator is first checked (see *symmetry_error*) to be a
2953 symmetry. The operators are then closed under products (two
2954 symmetries imply their product, which need not be given, as
2955 :math:`S = TC` in the tenfold way); a product of two operators of
2956 the same kind must agree with the given one up to a phase (no
2957 unitary symmetry left: block-diagonalize it first). The class
2958 follows from the kinds present, the squares :math:`UU^* = \pm1` of
2959 the four antiunitary kinds, and the commutation of
2960 :math:`S` (or :math:`\eta`) with them: :math:`U_SU = \epsilon\,UU_S^*`
2961 gives the subscript :math:`\epsilon` (for classes with TRS and PHS,
2962 the first refers to TRS, the second to PHS; in class AIII it is
2963 :math:`U_\Gamma U_S = \epsilon\,U_SU_\Gamma`).
2965 The name is that of the tables of Kawabata et al.: an AZ class
2966 ('AI', ...), an AZ\ :math:`^\dagger` class ('AI†', ...), an AZ class
2967 with SLS ('A + S', 'BDI + S-+', ...) when :math:`S` is present, or
2968 with pseudo-Hermiticity ('A + η', 'D + η+', ...) when :math:`\eta`
2969 is. These 54 names describe 38 classes: :math:`H \to iH` turns TRS
2970 into PHS\ :math:`^\dagger` and CS into pH, keeping the topology but
2971 exchanging real and imaginary line gaps, so that D† is AI, C† is AII,
2972 AII + S- is AI + S-, DIII + S-+ is BDI + S-+, A + η is AIII, and so
2973 on (their Tables XII and XIV). Since the tables differ for such
2974 partners, the name keeps the symmetries as given unless *canonical*.
2976 Parameters
2977 ----------
2978 time_reversal : Matrix :math:`U_T` (TRS), or None (absent)
2979 particle_hole : Matrix :math:`U_C` (PHS), or None (absent)
2980 chiral : Matrix :math:`U_\Gamma` (CS), or None (absent)
2981 time_reversal_dagger : Matrix :math:`U_{C_+}` (TRS†), or None (absent)
2982 particle_hole_dagger : Matrix :math:`U_{T_-}` (PHS†), or None (absent)
2983 sublattice : Matrix :math:`U_S` (SLS), or None (absent)
2984 pseudo_hermiticity : Matrix :math:`U_\eta` (pH), or None (absent)
2985 canonical : Boolean, default False
2986 Return the representative of
2987 the 38 classes instead (AZ names before AZ† ones, η classes as
2988 their equivalents without η).
2989 nk : Positive integer, default 5
2990 See *symmetry_error*.
2991 tol : Positive real, default 1e-8
2992 Tolerance of the checks.
2994 Returns
2995 -------
2996 name : String, for instance 'A', 'AI†', 'AIII + S-', 'BDI + S+-' or 'C + η+'
2998 Notes
2999 -----
3000 For a Hermitian model :math:`H^* = H^T` and :math:`H^\dagger = H`,
3001 so TRS is also TRS†, PHS also PHS†, CS also SLS, and :math:`\eta = 1`
3002 is a pseudo-Hermiticity: a Hermitian class BDI is 'BDI + S++' once
3003 all of them are given.
3004 '''
3005 error_handling.positive_real(tol, 'tol')
3006 error_handling.boolean(canonical, 'canonical')
3007 k0 = np.zeros(self.dim)
3008 given = {}
3009 ops = (time_reversal, particle_hole, chiral, time_reversal_dagger, particle_hole_dagger,
3010 sublattice, pseudo_hermiticity)
3011 for op, kind in zip(ops, _NH_KINDS):
3012 if op is not None:
3013 conj, transpose, minus = kind
3014 k_map = 'minus' if conj != transpose else 'identity'
3015 err = self.symmetry_error(op, k_map, bool(conj), bool(minus), nk, transpose=bool(transpose))
3016 error_handling.is_symmetry(err, tol)
3017 given[kind] = self._operator(op, k0)
3018 group = self._symmetry_group(given, tol)
3019 squares = {}
3020 for kind, u in group.items():
3021 if kind[0] != kind[1]:
3022 square = u @ u.conj()
3023 error_handling.square(square, tol)
3024 squares[kind] = int(np.sign(square[0, 0].real))
3025 t_p, c_m, c_p, t_m = (squares.get(kind, 0) for kind in ((1, 0, 0), (0, 1, 1), (0, 1, 0), (1, 0, 1)))
3026 az = {(1, 1): 'BDI', (-1, 1): 'DIII', (-1, -1): 'CII', (1, -1): 'CI', (1, 0): 'AI',
3027 (-1, 0): 'AII', (0, 1): 'D', (0, -1): 'C'}
3029 def sign(eps):
3030 return '+' if eps > 0 else '-'
3032 if (0, 0, 1) in group: # SLS
3033 if not squares:
3034 if (1, 1, 1) not in group:
3035 name = 'A + S'
3036 else:
3037 u_g, u_s = group[(1, 1, 1)], group[(0, 0, 1)]
3038 eps = np.trace((u_s @ u_g).conj().T @ u_g @ u_s).real
3039 name = 'AIII + S' + sign(eps)
3040 elif not (c_p or c_m):
3041 name = az[(t_p, 0)] + ' + S' + sign(t_p * t_m)
3042 elif not (t_p or t_m):
3043 name = az[(0, c_m)] + ' + S' + sign(c_p * c_m)
3044 else:
3045 name = az[(t_p, c_m)] + ' + S' + sign(t_p * t_m) + sign(c_p * c_m)
3046 elif (1, 1, 0) in group: # pseudo-Hermiticity
3047 if not squares:
3048 name = 'A + η'
3049 elif t_p:
3050 name = az[(t_p, 0)] + ' + η' + sign(t_p * c_p)
3051 else:
3052 name = az[(0, c_m)] + ' + η' + sign(t_m * c_m)
3053 elif t_p or c_m:
3054 name = az[(t_p, c_m)]
3055 elif c_p or t_m:
3056 name = az[(c_p, t_m)] + '†'
3057 else:
3058 name = 'AIII' if (1, 1, 1) in group else 'A'
3059 return _NH_CANONICAL.get(name, name) if canonical else name
3061 def _symmetry_group(self, given: dict, tol: float) -> dict:
3062 '''
3063 Private method. Close the symmetry operators *given* ({kind: matrix},
3064 kind = (conjugation, transposition, minus sign), see *_NH_KINDS*)
3065 under products, checking that two products of the same kind agree up
3066 to a phase. Unitary kinds come back normalized to :math:`U^2 = 1`.
3067 '''
3068 group = {(0, 0, 0): np.eye(self.norb, dtype='c16')}
3069 todo = list(given.items())
3070 while todo:
3071 kind, u = todo.pop()
3072 if kind in group:
3073 error_handling.symmetry_algebra(group[kind], u, tol)
3074 continue
3075 group[kind] = u
3076 for other, v in list(group.items()):
3077 product = tuple(a ^ b for a, b in zip(kind, other))
3078 # first one, then the other: U_2 f(U_1), f = conjugation for the antiunitary kinds
3079 todo.append((product, v @ (u.conj() if other[0] != other[1] else u)))
3080 todo.append((product, u @ (v.conj() if kind[0] != kind[1] else v)))
3081 for kind, u in group.items():
3082 if kind[0] == kind[1]:
3083 group[kind] = u / np.sqrt((u @ u)[0, 0])
3084 return group
3086 def quantum_geometric_tensor(
3087 self, bands: int | list[int], k: ArrayLike, dk: float = 1e-4, positions: bool = True,
3088 ) -> NDArray[np.complex128]:
3089 r'''
3090 Get the quantum geometric tensor of a group of bands at :math:`\mathbf{k}`,
3092 .. math::
3094 Q_{\mu\nu} = \mathrm{Tr}\left[P\,\partial_\mu P\,\partial_\nu P\right]
3095 = \sum_{n\in\mathrm{bands}}\langle\partial_\mu u_n|(1-P)|\partial_\nu u_n\rangle\, ,
3097 with :math:`P(\mathbf{k})` the projector on the bands (gauge
3098 invariant, so degeneracies within the group are harmless), from
3099 central differences. Many k-points are evaluated at once, with one
3100 batched diagonalization (see *set_workers*). Its real
3101 part is the quantum metric :math:`g_{\mu\nu}` (Provost and Vallee,
3102 1980), minus twice its imaginary part the Berry curvature
3103 :math:`\Omega_{\mu\nu} = -2\,\mathrm{Im}\,Q_{\mu\nu}` (the convention of
3104 *berry_curvature*, per unit k-area here). They obey
3105 :math:`\sqrt{\det g} \ge |\Omega_{xy}|/2` in 2D.
3107 Parameters
3108 ----------
3109 bands : Band index, or list of band indices
3110 k : k point (see *get_ham*), or real array of shape (nk, dim)
3111 One or many k-points.
3112 dk : Positive real, default 1e-4
3113 Finite-difference step.
3114 positions : Boolean, default True
3115 Include the orbital
3116 positions in the Bloch functions. The local metric and curvature
3117 depend on this choice (the Berry curvature of *berry_curvature* is
3118 that of ``positions=False``); integrals such as the Chern number
3119 do not.
3121 Returns
3122 -------
3123 Q : Complex ndarray, shape (dim, dim), in the k coordinates of *get_ham*
3124 Shape (nk, dim, dim) for many k-points.
3125 '''
3126 if isinstance(bands, int):
3127 bands = [bands]
3128 error_handling.band_indices(bands, self.norb)
3129 error_handling.k_points(k, 'k', self.dim)
3130 error_handling.positive_real(dk, 'dk')
3131 error_handling.boolean(positions, 'positions')
3132 k = np.asarray(k, dtype='f8')
3133 ks = np.atleast_2d(k)
3134 # k, then k + dk e_mu and k - dk e_mu for every mu
3135 shifts = np.concatenate([np.zeros((1, self.dim)), dk * np.eye(self.dim), -dk * np.eye(self.dim)])
3136 kk = (ks[:, None, :] + shifts[None]).reshape(-1, self.dim)
3137 v = self._subspaces(kk, bands)
3138 if positions:
3139 v = np.exp(-1j * (kk @ self.k_basis.T) @ self.orbital_positions().T)[:, :, None] * v
3140 proj = (v @ v.conj().transpose(0, 2, 1)).reshape(len(ks), 1 + 2 * self.dim, self.norb, self.norb)
3141 dp = (proj[:, 1:1 + self.dim] - proj[:, 1 + self.dim:]) / (2 * dk)
3142 q = np.einsum('kab,kmbc,knca->kmn', proj[:, 0], dp, dp)
3143 return q if k.ndim == 2 else q[0]
3145 def surface_spectral_function(
3146 self, ks: ArrayLike, energies: ArrayLike, direction: int, side: int = 1,
3147 eta: float = 1e-2, bulk: bool = False, max_iter: int = 10000,
3148 ) -> NDArray[np.float64]:
3149 r'''
3150 Get the spectral function of the surface of a semi-infinite crystal,
3152 .. math::
3154 A_s(\mathbf{k}_\parallel, E) = -\frac{1}{\pi}\,
3155 \mathrm{Im}\,\mathrm{Tr}\, G_{00}(\mathbf{k}_\parallel, E + i\eta)\, ,
3157 with :math:`G_{00}` the Green's function of the outermost unit cell:
3158 the crystal fills the cells :math:`n\,\mathbf{a}_d`,
3159 :math:`n = 0, 1, 2, \dots` along the primitive vector
3160 :math:`\mathbf{a}_d` (``lat.prim_vec[direction]``), times *side*, and
3161 is periodic along the others. It comes from the Lopez Sancho-Rubio
3162 decimation (see *transport.surface_green*) of the chain of principal
3163 layers (as many cells as the longest hopping along
3164 :math:`\mathbf{a}_d` spans), so the crystal is truly semi-infinite: no
3165 finite slab, and no states of the opposite surface. It is the
3166 tight-binding picture of an ARPES map -- surface states (edge states,
3167 in 2D; the end state of a 1D chain) show up as sharp lines inside the
3168 bulk gaps. With *bulk*, the same quantity for a unit cell deep inside
3169 the crystal: the bulk bands projected on the surface Brillouin zone.
3171 Parameters
3172 ----------
3173 ks : Real array, shape (nk, dim)
3174 k-points, in the coordinates
3175 of *get_ham* (e.g. from a *k_path* of the surface Brillouin zone).
3176 Only :math:`\mathbf{k}_\parallel` matters: the result does not
3177 depend on the component along the reciprocal vector
3178 :math:`\mathbf{b}_d`.
3179 energies : Real number or real array
3180 Energies :math:`E`.
3181 direction : Integer, between 0 and dim-1
3182 The primitive vector
3183 normal to the surface (the one along which the crystal is
3184 semi-infinite).
3185 side : +1 or -1, default 1
3186 The crystal extends along
3187 :math:`+\mathbf{a}_d` (the surface faces :math:`-\mathbf{a}_d`) or
3188 along :math:`-\mathbf{a}_d` (the opposite surface).
3189 eta : Positive real, default 1e-2
3190 Broadening: the
3191 width of the lines, and the decay length (about :math:`v/\eta`
3192 cells) over which the bulk is resolved.
3193 bulk : Boolean, default False
3194 If True, the spectral
3195 function of a bulk unit cell instead of the surface one.
3196 max_iter : Positive integer, default 10000
3197 Maximum
3198 number of decimation steps.
3200 Returns
3201 -------
3202 A : Real ndarray, shape (nk, len(energies))
3204 Examples
3205 --------
3206 Example usage::
3208 # the (001) surface of a 3D model, along a path of the surface zone
3209 A = ks.surface_spectral_function(k_path, np.linspace(-1, 1, 201), direction=2)
3210 '''
3211 error_handling.no_overlap(self._overlap_hop)
3212 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
3213 error_handling.ks(ks, self.dim)
3214 error_handling.frequencies(energies, 'energies')
3215 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
3216 error_handling.direction(direction, self.dim)
3217 error_handling.surface_side(side)
3218 error_handling.positive_real(eta, 'eta')
3219 error_handling.boolean(bulk, 'bulk')
3220 error_handling.positive_int(max_iter, 'max_iter')
3221 # H_m(k): the hoppings from a cell to the cell m a_d away, with the
3222 # phase of their in-plane part only (independent of k along b_d)
3223 a_dir = np.array(self.lat.prim_vec[direction], dtype='f8')
3224 layers: dict[int, list] = {}
3225 for (i, j, cell, t), (_, _, R, _) in zip(self._hop_cells(), self._hop):
3226 layers.setdefault(cell[direction], []).append((i, j, R - cell[direction] * a_dir, t))
3227 k_cart = ks @ self.k_basis.T
3228 blocks: dict[int, NDArray] = {m: self._bloch_sum(hops, k_cart)[0] for m, hops in layers.items()}
3229 blocks[0] = blocks.get(0, 0.) + (np.diag(self.onsite) + self._onsite_offdiag)[None]
3230 width = max([abs(m) for m in layers] + [1])
3231 norb, nk = self.norb, len(ks)
3233 def layer(shift):
3234 # block (p, q) couples cell p of a principal layer to cell q of
3235 # the layer *shift* cells further into the crystal
3236 out = np.zeros((nk, width, norb, width, norb), 'c16')
3237 for p in range(width):
3238 for q in range(width):
3239 m = side * (q - p + shift)
3240 if m in blocks:
3241 out[:, p, :, q, :] = blocks[m]
3242 return out.reshape(nk, width * norb, width * norb)
3243 h0, alpha, beta = layer(0), layer(width), layer(-width)
3244 # one task per k-point and chunk of energies, enough tasks for every thread
3245 n_e = max(1, min(len(energies), -(-self.workers // nk)))
3246 tasks = [(n, part) for n in range(nk) for part in np.array_split(np.arange(len(energies)), n_e)]
3248 def task(n_part):
3249 n, part = n_part
3250 gs, gb = _green_surface_bulk(energies[part] + 1j * eta, h0[n], alpha[n], beta[n], 1e-12,
3251 max_iter, 'surface_spectral_function')
3252 g = gb if bulk else gs
3253 return -np.trace(g[:, :norb, :norb], axis1=1, axis2=2).imag / PI
3254 if self.workers > 1 and len(tasks) > 1:
3255 with ThreadPoolExecutor(min(self.workers, len(tasks))) as pool:
3256 out = list(pool.map(task, tasks))
3257 else:
3258 out = [task(t) for t in tasks]
3259 spec = np.empty((nk, len(energies)))
3260 for (n, part), a in zip(tasks, out):
3261 spec[n, part] = a
3262 return spec
3264 # ------------------------------------------------------------------
3265 # Projections, constant-energy contours and plot helpers
3266 # ------------------------------------------------------------------
3268 def spin_operator(self, axis: str) -> NDArray[np.complex128]:
3269 r'''
3270 Get the Pauli matrix :math:`\sigma_{axis}` acting on the spin of
3271 every site, :math:`\mathbb{1}_{sites}\otimes\sigma_{axis}`, for a model
3272 built with ``spin=True`` (for *band_weights* and *spin_texture*).
3274 Parameters
3275 ----------
3276 axis : 'x', 'y' or 'z'
3278 Returns
3279 -------
3280 op : Complex ndarray, shape (norb, norb)
3281 '''
3282 error_handling.spinful(self.spin)
3283 error_handling.spin_axis(axis)
3284 return np.kron(np.eye(self.n_sites), PAULI[axis])
3286 def _projector(self, projector) -> NDArray[np.complex128]:
3287 '''
3288 Private method. The (norb, norb) operator of *band_weights* from an
3289 orbital index, a list of them, a sublattice tag or a matrix.
3290 '''
3291 error_handling.projector(projector, self.norb, self.tags)
3292 if isinstance(projector, np.ndarray):
3293 return projector.astype('c16')
3294 if isinstance(projector, str):
3295 sites = np.nonzero(self.tags == projector)[0]
3296 projector = np.concatenate([2 * sites, 2 * sites + 1]) if self.spin else sites
3297 diag = np.zeros(self.norb)
3298 diag[projector] = 1.
3299 return np.diag(diag).astype('c16')
3301 def _expectations(self, ks: NDArray[np.float64], ops: list) -> NDArray[np.float64]:
3302 r'''
3303 Private method. :math:`\mathrm{Re}\langle u_n(\mathbf{k})|O|u_n(\mathbf{k})\rangle
3304 /\langle u_n|u_n\rangle` for each operator of *ops*, band and k-point:
3305 shape (len(ops), nk, norb), from one diagonalization. Within a group
3306 of degenerate bands (energies within 1e-8 of each other, relative to
3307 the largest one), each band gets the average over the group: the
3308 trace of O over the degenerate subspace divided by its dimension,
3309 which does not depend on the eigenvectors the solver returns.
3310 '''
3311 en, vn = self._eigs(ks, eigenvec=True)
3312 norm = np.sum(np.abs(vn) ** 2, axis=1)
3313 values = np.array([np.einsum('kin,ij,kjn->kn', vn.conj(), op, vn).real / norm for op in ops])
3314 # label the degenerate groups (consecutive bands, sorted by energy)
3315 tol = 1e-8 * max(1., np.abs(en).max())
3316 new_group = np.concatenate([np.ones((len(en), 1), bool),
3317 np.abs(np.diff(en, axis=1)) > tol], axis=1)
3318 groups = np.cumsum(new_group.ravel()) - 1
3319 counts = np.bincount(groups)
3320 return np.array([(np.bincount(groups, v.ravel()) / counts)[groups].reshape(v.shape)
3321 for v in values])
3323 def band_weights(self, projector, ks: ArrayLike | None = None) -> NDArray[np.float64]:
3324 r'''
3325 Get the weight of each band on a set of orbitals, or the expectation
3326 value of an operator, for "fat band" plots (see *plot_bands*):
3328 .. math::
3330 w_n(\mathbf{k}) = \langle u_n(\mathbf{k})|P|u_n(\mathbf{k})\rangle\, ,
3332 with :math:`P` the projector on the chosen orbitals (a weight
3333 between 0 and 1; the weights on all the sublattices add up to 1), or
3334 any Hermitian operator, e.g. ``ks.spin_operator('z')`` for the spin
3335 polarization (between -1 and 1). Within a group of degenerate bands
3336 (e.g. Kramers pairs with inversion and time reversal), the weight of
3337 a single eigenvector depends on the choice of eigenvectors; each band
3338 of the group gets instead the average over the group, which does not
3339 (a Kramers pair has zero spin polarization).
3340 For a non-Hermitian model, the right eigenvectors are used; with an
3341 overlap (*set_overlap*), :math:`v^\dagger Pv/v^\dagger v`.
3343 Parameters
3344 ----------
3345 projector
3346 Orbital index, list of orbital indices, sublattice
3347 tag (all the orbitals, both spins, of the sites with that tag),
3348 or a (norb, norb) Hermitian matrix.
3349 ks : Real array, shape (nk, dim)
3350 Default value None: the
3351 k-points of the last *get_bands* or *k_path*.
3353 Returns
3354 -------
3355 weights : Real ndarray, shape (nk, norb), ordered as the bands
3357 Examples
3358 --------
3359 Example usage::
3361 ks.k_path(points, nk=100)
3362 ks.plot_bands(weights=ks.band_weights('a'))
3363 '''
3364 op = self._projector(projector)
3365 if ks is None:
3366 error_handling.empty_ndarray(self.en, 'get_bands or k_path')
3367 ks = self.ks
3368 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
3369 error_handling.ks(ks, self.dim)
3370 return self._expectations(ks, [op])[0]
3372 def spin_texture(self, ks: ArrayLike, band: int) -> NDArray[np.float64]:
3373 r'''
3374 Get the spin expectation values
3375 :math:`\langle\boldsymbol\sigma\rangle_n(\mathbf{k}) =
3376 \langle u_n|\boldsymbol\sigma|u_n\rangle` of a band, for a model built
3377 with ``spin=True``. Their length is 1 if the band is fully polarized,
3378 less if spin-orbit coupling entangles the spin with the orbitals.
3379 Within a group of degenerate bands, each band gets the average spin
3380 of the group (see *band_weights*): zero for a Kramers pair, with both
3381 inversion and time reversal.
3383 Parameters
3384 ----------
3385 ks : Real array, shape (nk, dim)
3386 k-points.
3387 band : Integer
3388 Band index.
3390 Returns
3391 -------
3392 spin
3393 Real ndarray, shape (nk, 3):
3394 :math:`(\langle\sigma_x\rangle, \langle\sigma_y\rangle, \langle\sigma_z\rangle)`.
3395 '''
3396 error_handling.spinful(self.spin)
3397 error_handling.band_index(band, self.norb)
3398 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
3399 error_handling.ks(ks, self.dim)
3400 ops = [self.spin_operator(a) for a in 'xyz']
3401 return self._expectations(ks, ops)[:, :, band].T
3403 def _fold(self, elements: NDArray[np.float64]) -> NDArray[np.float64]:
3404 '''
3405 Private method. Bring segments or triangles (fractional coordinates,
3406 shape (n, nv, dim)) into the first Brillouin zone (the Wigner-Seitz
3407 cell of the reciprocal lattice, by their centres), in the k
3408 coordinates of *get_ham*.
3409 '''
3410 elements = elements - np.floor(elements.mean(axis=1) + 0.5)[:, None, :]
3411 cart = elements @ self.rec_vec_k
3412 shifts = np.array(list(product((-1, 0, 1), repeat=self.dim))) @ self.rec_vec_k
3413 dist = np.linalg.norm(cart.mean(axis=1)[:, None, :] - shifts[None], axis=2)
3414 return cart - shifts[np.argmin(dist, axis=1)][:, None, :]
3416 def fermi_surface(
3417 self, energy: float = 0., nk: int | tuple[int, ...] | None = None,
3418 bands: int | list[int] | None = None,
3419 ) -> list[NDArray[np.float64]]:
3420 r'''
3421 Get the constant-energy surfaces :math:`E_n(\mathbf{k}) = E` (the
3422 Fermi surface at :math:`E = E_F`): lines in a 2D Brillouin zone,
3423 surfaces in 3D. The bands are interpolated linearly in the triangles
3424 (tetrahedra in 3D) of a uniform mesh, the same decomposition as
3425 *dos.tetrahedron_dos*, and the level set of each simplex is a segment
3426 (a triangle or a quadrilateral, split in two, in 3D): marching
3427 triangles and tetrahedra, exact for the interpolated bands. Each
3428 piece is brought into the first Brillouin zone, the Wigner-Seitz
3429 cell around :math:`\Gamma`. For a non-Hermitian model, the real part
3430 of the energies is used.
3432 Parameters
3433 ----------
3434 energy : Real number, default 0
3435 The energy :math:`E`.
3436 nk : Positive integer, or tuple of *dim* positive integers, default None: 100 in 2D, 40 in 3D
3437 Number of k-points along
3438 each reciprocal lattice vector.
3439 bands : Band index, or list of band indices
3440 Default value
3441 None: all bands.
3443 Returns
3444 -------
3445 surfaces
3446 List, one entry per band of *bands*: a real
3447 ndarray of shape (n, 2, 2) in 2D (n segments, each two k-points)
3448 or (n, 3, 3) in 3D (n triangles, each three k-points), in the
3449 k coordinates of *get_ham*. Empty (n = 0) if the band does not
3450 cross *energy*.
3452 Examples
3453 --------
3454 Example usage::
3456 # the Fermi surface of the half-filled square lattice
3457 segments = sq.fermi_surface(0.)[0]
3458 '''
3459 error_handling.dim_min(self.dim, 2)
3460 error_handling.real_number(energy, 'energy')
3461 if nk is None:
3462 nk = 100 if self.dim == 2 else 40
3463 error_handling.nk(nk, self.dim)
3464 if isinstance(nk, int):
3465 nk = (nk,) * self.dim
3466 if bands is None:
3467 bands = list(range(self.norb))
3468 elif isinstance(bands, int):
3469 bands = [bands]
3470 error_handling.band_indices(bands, self.norb)
3471 _, ks = self.mesh_grid(nk)
3472 en = self._eigs(ks).real.reshape(*nk, self.norb)
3473 return [self._fold(_level_set(en[..., n], energy)) for n in bands]
3475 def _k_labels(self) -> list[str]:
3476 '''
3477 Private method. Axis labels of the k coordinates of *get_ham*:
3478 '$k_x$' for a coordinate along the Cartesian x axis (the slab of a
3479 3D model along x and z gives '$k_x$', '$k_z$'), '$k_1$', ... otherwise.
3480 '''
3481 labels = []
3482 for n, col in enumerate(self.k_basis.T):
3483 axis = np.nonzero(np.isclose(np.abs(col), 1.))[0]
3484 labels.append('$k_{}$'.format('xyz'[axis[0]] if len(axis) else n + 1))
3485 return labels
3487 def _zone_axes(self, ax, fs: float) -> None:
3488 '''
3489 Private method. Draw the first Brillouin zone of a 2D model on *ax*,
3490 and label and scale its axes.
3491 '''
3492 zone = _zone_polygon(self.rec_vec_k)
3493 ax.plot(*np.vstack([zone, zone[:1]]).T, 'k', lw=1)
3494 ax.set_aspect('equal')
3495 pad = 0.05 * np.ptp(zone, axis=0).max()
3496 ax.set_xlim(zone[:, 0].min() - pad, zone[:, 0].max() + pad)
3497 ax.set_ylim(zone[:, 1].min() - pad, zone[:, 1].max() + pad)
3498 ax.set_xlabel(self._k_labels()[0], fontsize=fs)
3499 ax.set_ylabel(self._k_labels()[1], fontsize=fs)
3501 def plot_fermi_surface(
3502 self, energy: float = 0., nk: int | tuple[int, ...] | None = None,
3503 bands: int | list[int] | None = None, lw: float = 2., alpha: float = 0.6,
3504 fs: float = 20, figsize: tuple[float, float] | None = None,
3505 ) -> Figure:
3506 '''
3507 Plot the constant-energy contours of *fermi_surface* (one color per
3508 band, with a legend if several bands cross the energy): lines inside
3509 the first Brillouin zone (drawn in black) in 2D, surfaces in 3D.
3511 Parameters
3512 ----------
3513 energy : Real number, default 0
3514 The energy.
3515 nk
3516 See *fermi_surface*.
3517 bands
3518 See *fermi_surface*.
3519 lw : Positive number, default 2
3520 Linewidth (2D).
3521 alpha : Real number in (0, 1], default 0.6
3522 Opacity of
3523 the surfaces (3D).
3524 fs : Positive number, default 20
3525 Fontsize.
3526 figsize : Tuple, default None
3527 Figure size.
3529 Returns
3530 -------
3531 fig : Figure
3532 '''
3533 error_handling.positive_real(lw, 'lw')
3534 error_handling.positive_real(alpha, 'alpha')
3535 error_handling.positive_real(fs, 'fs')
3536 error_handling.tuple_2elem(figsize, 'figsize')
3537 surfaces = self.fermi_surface(energy, nk, bands)
3538 bands = list(range(self.norb)) if bands is None else np.atleast_1d(bands).tolist()
3539 fig = plt.figure(figsize=figsize)
3540 if self.dim == 2:
3541 ax = fig.add_subplot()
3542 for n, (band, segs) in enumerate(zip(bands, surfaces)):
3543 if len(segs):
3544 ax.add_collection(LineCollection(list(segs), colors='C{}'.format(n % 10), lw=lw,
3545 label='band {}'.format(band)))
3546 self._zone_axes(ax, fs)
3547 else:
3548 ax = fig.add_subplot(projection='3d')
3549 for n, (band, tris) in enumerate(zip(bands, surfaces)):
3550 if len(tris):
3551 ax.add_collection3d(Poly3DCollection(tris, facecolor='C{}'.format(n % 10),
3552 alpha=min(alpha, 1.), edgecolor='none',
3553 label='band {}'.format(band)))
3554 extent = np.abs(np.concatenate([t.reshape(-1, 3) for t in surfaces] +
3555 [self.rec_vec_k / 2])).max()
3556 ax.set(xlim=(-extent, extent), ylim=(-extent, extent), zlim=(-extent, extent))
3557 ax.set_box_aspect((1, 1, 1))
3558 ax.set_xlabel(self._k_labels()[0], fontsize=fs)
3559 ax.set_ylabel(self._k_labels()[1], fontsize=fs)
3560 ax.set_zlabel(self._k_labels()[2], fontsize=fs)
3561 if sum(len(s) > 0 for s in surfaces) > 1:
3562 ax.legend(fontsize=0.6 * fs)
3563 ax.set_title('$E = {:g}$'.format(energy), fontsize=fs)
3564 plt.draw()
3565 return fig
3567 def plot_spin_texture(
3568 self, band: int, energy: float = 0., nk: int = 100, n_arrows: int = 40,
3569 cmap: str = 'RdBu_r', lw: float = 1., fs: float = 20,
3570 figsize: tuple[float, float] | None = None,
3571 ) -> Figure:
3572 r'''
3573 Plot the spin texture of a band of a 2D spinful model on its
3574 constant-energy contour (see *fermi_surface*), zoomed on the
3575 contour: arrows of
3576 :math:`(\langle\sigma_x\rangle, \langle\sigma_y\rangle)` along the
3577 contour, colored by :math:`\langle\sigma_z\rangle` (see *spin_texture*).
3578 Rashba coupling, for instance, locks the spin perpendicular to
3579 :math:`\mathbf{k}`, winding once around the contour.
3581 Parameters
3582 ----------
3583 band : Integer
3584 Band index.
3585 energy : Real number, default 0
3586 Energy of the contour.
3587 nk : Positive integer, default 100
3588 Mesh of the contour.
3589 n_arrows : Positive integer, default 40
3590 Approximate
3591 number of arrows, evenly spaced along the contour. A fully
3592 polarized in-plane spin is drawn 80% as long as their spacing.
3593 cmap
3594 Default value 'RdBu_r'. Colormap of :math:`\langle\sigma_z\rangle`.
3595 lw : Positive number, default 1
3596 Linewidth of the contour.
3597 fs : Positive number, default 20
3598 Fontsize.
3599 figsize : Tuple, default None
3600 Figure size.
3602 Returns
3603 -------
3604 fig : Figure
3605 '''
3606 error_handling.dim_2(self.dim)
3607 error_handling.spinful(self.spin)
3608 error_handling.band_index(band, self.norb)
3609 error_handling.positive_int(nk, 'nk')
3610 error_handling.positive_int(n_arrows, 'n_arrows')
3611 error_handling.positive_real(lw, 'lw')
3612 error_handling.positive_real(fs, 'fs')
3613 error_handling.tuple_2elem(figsize, 'figsize')
3614 segs = self.fermi_surface(energy, nk, [band])[0]
3615 error_handling.contour_found(len(segs), energy)
3616 # arrows at segment midpoints, at least one contour length / n_arrows apart
3617 mids = segs.mean(axis=1)
3618 spacing = np.linalg.norm(segs[:, 1] - segs[:, 0], axis=1).sum() / n_arrows
3619 keep = [0]
3620 for i in range(1, len(mids)):
3621 if np.min(np.linalg.norm(mids[keep] - mids[i], axis=1)) >= spacing:
3622 keep.append(i)
3623 spin = self.spin_texture(mids[keep], band)
3624 fig, ax = plt.subplots(figsize=figsize)
3625 ax.add_collection(LineCollection(list(segs), colors='0.6', lw=lw))
3626 arrows = ax.quiver(mids[keep, 0], mids[keep, 1], spin[:, 0], spin[:, 1], spin[:, 2],
3627 cmap=cmap, norm=Normalize(-1., 1.), pivot='mid', edgecolor='k',
3628 linewidth=0.5, angles='xy', scale_units='xy', scale=1.25 / spacing)
3629 fig.colorbar(arrows, ax=ax).set_label(r'$\langle\sigma_z\rangle$', fontsize=fs)
3630 self._zone_axes(ax, fs)
3631 # zoom on the contour (within the zone)
3632 low, high = segs.reshape(-1, 2).min(axis=0), segs.reshape(-1, 2).max(axis=0)
3633 pad = 0.1 * (high - low).max() + spacing
3634 ax.set_xlim(max(ax.get_xlim()[0], low[0] - pad), min(ax.get_xlim()[1], high[0] + pad))
3635 ax.set_ylim(max(ax.get_ylim()[0], low[1] - pad), min(ax.get_ylim()[1], high[1] + pad))
3636 ax.set_title('Spin texture of band {}, $E = {:g}$'.format(band, energy), fontsize=fs)
3637 plt.draw()
3638 return fig
3640 def plot_berry_curvature(
3641 self, bands: int | list[int], nk: int | tuple[int, int] = 60,
3642 plane: tuple[int, int] = (0, 1), k_fixed: float = 0., cmap: str = 'RdBu_r',
3643 fs: float = 20, figsize: tuple[float, float] | None = None,
3644 ) -> Figure:
3645 r'''
3646 Plot the Berry curvature of *berry_curvature* as a density,
3647 :math:`\Omega(\mathbf{k})` = flux of each plaquette / its area. In
3648 2D, the plaquettes are brought into the first Brillouin zone (drawn
3649 in black), in the k coordinates of *get_ham*. In 3D, they fill the
3650 parallelogram spanned by the two reciprocal vectors of *plane*, in
3651 orthonormal coordinates of the plane (:math:`k_1` along
3652 :math:`\mathbf{b}_i`). The title gives the Chern number, the total
3653 flux over :math:`2\pi`.
3655 Parameters
3656 ----------
3657 bands
3658 See *berry_curvature*.
3659 nk
3660 Default value 60. See *berry_curvature*.
3661 plane
3662 See *berry_curvature*.
3663 k_fixed
3664 See *berry_curvature*.
3665 cmap
3666 Default value 'RdBu_r'. Colormap, centred on zero.
3667 fs : Positive number, default 20
3668 Fontsize.
3669 figsize : Tuple, default None
3670 Figure size.
3672 Returns
3673 -------
3674 fig : Figure
3675 '''
3676 error_handling.positive_real(fs, 'fs')
3677 error_handling.tuple_2elem(figsize, 'figsize')
3678 curv = self.berry_curvature(bands, nk, plane, k_fixed)
3679 n1, n2 = curv.shape
3680 # the corners of each plaquette, counterclockwise, in fractional coordinates
3681 f1, f2 = np.meshgrid(np.arange(n1) / n1, np.arange(n2) / n2, indexing='ij')
3682 corner = np.stack([f1.ravel(), f2.ravel()], axis=-1)[:, None, :]
3683 quads = corner + np.array([[0., 0.], [1. / n1, 0.], [1. / n1, 1. / n2], [0., 1. / n2]])
3684 b_i, b_j = self.rec_vec_k[plane[0]], self.rec_vec_k[plane[1]]
3685 if self.dim == 2:
3686 quads = self._fold(quads)
3687 else:
3688 # orthonormal coordinates of the plane, the first along b_i
3689 u = b_i / np.linalg.norm(b_i)
3690 v = b_j - (b_j @ u) * u
3691 b_i, b_j = np.array([b_i @ u, 0.]), np.array([b_j @ u, b_j @ v / np.linalg.norm(v)])
3692 quads = quads @ np.array([b_i, b_j])
3693 area = abs(b_i[0] * b_j[1] - b_i[1] * b_j[0]) / (n1 * n2)
3694 omega = curv.ravel() / area
3695 vmax = np.abs(omega).max() or 1.
3696 fig, ax = plt.subplots(figsize=figsize)
3697 patches = PolyCollection(quads, array=omega, cmap=cmap, norm=Normalize(-vmax, vmax),
3698 edgecolors='face', linewidths=0.2)
3699 ax.add_collection(patches)
3700 fig.colorbar(patches, ax=ax).set_label(r'$\Omega(\mathbf{k})$', fontsize=fs)
3701 if self.dim == 2:
3702 self._zone_axes(ax, fs)
3703 else:
3704 ax.autoscale_view()
3705 ax.set_aspect('equal')
3706 ax.set_xlabel('$k_1$', fontsize=fs)
3707 ax.set_ylabel('$k_2$', fontsize=fs)
3708 ax.set_title('Berry curvature, $C = {:.3f}$'.format(curv.sum() / (2 * PI)), fontsize=fs)
3709 fig.set_layout_engine('tight')
3710 plt.draw()
3711 return fig
3713 def plot_wannier_flow(
3714 self, bands: int | list[int], nk: int = 100, nk_perp: int = 51, direction: int = 0,
3715 flow: int | None = None, k_fixed: float = 0., positions: bool = True,
3716 k_range: tuple[float, float] = (0., 1.), ms: float = 4., c: str = 'b', fs: float = 20,
3717 figsize: tuple[float, float] | None = None,
3718 ) -> Figure:
3719 r'''
3720 Plot the flow of the hybrid Wannier centres (the Wilson-loop
3721 spectrum) of *wannier_flow*: the centres along *direction*, in
3722 units of :math:`\mathbf{a}_{direction}`, as the loop moves along
3723 the reciprocal vector *flow*. They wind :math:`C` times in a Chern
3724 insulator, and the Kramers pairs swap partners in a
3725 :math:`\mathbb{Z}_2` topological insulator.
3727 Parameters
3728 ----------
3729 bands
3730 See *wannier_flow*.
3731 nk
3732 See *wannier_flow*.
3733 nk_perp
3734 See *wannier_flow*.
3735 direction
3736 See *wannier_flow*.
3737 flow
3738 See *wannier_flow*.
3739 k_fixed
3740 See *wannier_flow*.
3741 positions
3742 See *wannier_flow*.
3743 k_range
3744 See *wannier_flow*.
3745 ms : Positive number, default 4
3746 Marker size.
3747 c
3748 Default value 'b'. Marker color.
3749 fs : Positive number, default 20
3750 Fontsize.
3751 figsize : Tuple, default None
3752 Figure size.
3754 Returns
3755 -------
3756 fig : Figure
3757 '''
3758 error_handling.positive_real(ms, 'ms')
3759 error_handling.positive_real(fs, 'fs')
3760 error_handling.tuple_2elem(figsize, 'figsize')
3761 fracs, centers = self.wannier_flow(bands, nk, nk_perp, direction, flow, k_fixed,
3762 positions, k_range)
3763 if flow is None:
3764 flow = [d for d in range(self.dim) if d != direction][0]
3765 fig, ax = plt.subplots(figsize=figsize)
3766 ax.plot(fracs, centers, 'o', c=c, ms=ms)
3767 ax.set_xlim(fracs[0], fracs[-1])
3768 ax.set_ylim(0., 1.)
3769 ax.set_xlabel(r'$k_{0}/|\mathbf{{b}}_{0}|$'.format(flow + 1), fontsize=fs)
3770 ax.set_ylabel(r'Wannier centres ($\mathbf{{a}}_{}$)'.format(direction + 1), fontsize=fs)
3771 fig.set_layout_engine('tight')
3772 plt.draw()
3773 return fig
3775 def plot_surface_spectral_function(
3776 self, points: list[ArrayLike], energies: ArrayLike, direction: int, nk: int = 60,
3777 side: int = 1, eta: float = 1e-2, bulk: bool = False,
3778 node_labels: list[str] | None = None, log: bool = True, cmap: str = 'magma',
3779 fs: float = 20, figsize: tuple[float, float] | None = None,
3780 ) -> Figure:
3781 r'''
3782 Plot the surface spectral function of *surface_spectral_function*
3783 along a path through the surface Brillouin zone (built as in
3784 *k_path*), as an ARPES-like map of :math:`E` against
3785 :math:`\mathbf{k}_\parallel`. Surface states are the sharp lines in
3786 the bulk gaps.
3788 Parameters
3789 ----------
3790 points : List of at least two k-points, the nodes of the path (see *k_path*)
3791 energies : Real array
3792 Energies :math:`E`.
3793 direction
3794 See *surface_spectral_function*.
3795 nk : Positive integer, default 60
3796 Number of k-points per path segment.
3797 side
3798 See *surface_spectral_function*.
3799 eta
3800 See *surface_spectral_function*.
3801 bulk
3802 See *surface_spectral_function*.
3803 node_labels : List of strings, default None
3804 Labels of *points*.
3805 log : Boolean, default True
3806 Logarithmic color scale, over
3807 four decades below the maximum (the surface states are much
3808 sharper than the bulk continuum).
3809 cmap
3810 Default value 'magma'. Colormap.
3811 fs : Positive number, default 20
3812 Fontsize.
3813 figsize : Tuple, default None
3814 Figure size.
3816 Returns
3817 -------
3818 fig : Figure
3819 '''
3820 error_handling.k_path_points(points, self.dim)
3821 error_handling.positive_int(nk, 'nk')
3822 error_handling.boolean(log, 'log')
3823 error_handling.positive_real(fs, 'fs')
3824 error_handling.tuple_2elem(figsize, 'figsize')
3825 if node_labels is not None:
3826 error_handling.ndarray(np.array(node_labels), 'node_labels', len(points))
3827 ks, dist, nodes = _path(points, nk)
3828 spec = self.surface_spectral_function(ks, energies, direction, side, eta, bulk)
3829 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
3830 vmax = spec.max() if spec.max() > 0 else 1.
3831 norm = LogNorm(vmax * 1e-4, vmax) if log else Normalize(0., vmax)
3832 fig, ax = plt.subplots(figsize=figsize)
3833 mesh = ax.pcolormesh(dist, energies, np.maximum(spec, vmax * 1e-4).T, cmap=cmap,
3834 norm=norm, shading='nearest')
3835 fig.colorbar(mesh, ax=ax).set_label(r'$A(\mathbf{k}_\parallel, E)$', fontsize=fs)
3836 for node in nodes:
3837 ax.axvline(node, color='w', lw=0.5)
3838 if node_labels is not None:
3839 ax.set_xticks(nodes)
3840 ax.set_xticklabels(node_labels, fontsize=fs)
3841 ax.set_ylabel('$E$', fontsize=fs)
3842 ax.set_title('Bulk' if bulk else 'Surface', fontsize=fs)
3843 fig.set_layout_engine('tight')
3844 plt.draw()
3845 return fig
3847 def plot_dos(
3848 self,
3849 nk: int | tuple[int, ...] = 30,
3850 broadening: float = 0.05,
3851 kernel: str = 'gaussian',
3852 e_grid: ArrayLike | None = None,
3853 fs: float = 20,
3854 lw: float = 2.,
3855 figsize: tuple[float, float] | None = None,
3856 ) -> Figure:
3857 '''
3858 Plot the density of states, obtained by diagonalizing
3859 :math:`H(\\mathbf{k})` over a uniform Brillouin-zone mesh: broadened
3860 (see *tbkit.dos.density_of_states*), or by the linear tetrahedron
3861 method (``kernel='tetrahedron'``, see *tbkit.dos.tetrahedron_dos*),
3862 which needs no broadening and keeps band edges and van Hove
3863 singularities sharp.
3865 Parameters
3866 ----------
3867 nk : Positive integer, or tuple of *dim* positive integers, default 30
3868 Number of k-points along each reciprocal
3869 lattice vector.
3870 broadening : Positive real number, default 0.05
3871 Kernel
3872 width (unused by the tetrahedron method).
3873 kernel : String, default 'gaussian'
3874 'gaussian',
3875 'lorentzian' or 'tetrahedron'.
3876 e_grid : Real ndarray, default None
3877 Energies at which to
3878 evaluate the density of states.
3879 fs : Positive number, default 20
3880 Fontsize.
3881 lw : Positive number, default 2
3882 Linewidth.
3883 figsize : Tuple, default None
3884 Figure size.
3886 Returns
3887 -------
3888 fig : Figure
3889 '''
3890 error_handling.positive_real(fs, 'fs')
3891 error_handling.positive_real(lw, 'lw')
3892 error_handling.tuple_2elem(figsize, 'figsize')
3893 error_handling.kspace_dos_kernel(kernel)
3894 en = self.mesh_bands(nk)
3895 if kernel != 'tetrahedron':
3896 return dos._plot_density_of_states(en, e_grid, broadening, kernel, fs, lw, figsize)
3897 nk = (nk,) * self.dim if isinstance(nk, int) else nk
3898 e_grid, rho = dos.tetrahedron_dos(en.reshape(*nk, self.norb), e_grid)
3899 return dos._plot_dos_curve(e_grid, rho, fs, lw, figsize)
3901 def plot_bands(
3902 self,
3903 node_labels: list[str] | None = None,
3904 fs: float = 20,
3905 lw: float = 2.,
3906 ms: float = 0.,
3907 c: str = 'b',
3908 lims: tuple[float, float] | None = None,
3909 figsize: tuple[float, float] | None = None,
3910 weights: ArrayLike | None = None,
3911 style: str = 'color',
3912 cmap: str | None = None,
3913 ) -> Figure:
3914 '''
3915 Plot the band structure computed by *k_path* or *get_bands*, against
3916 the cumulative distance along the k-points (vertical lines mark the
3917 nodes of a *k_path*). For a non-Hermitian model, the real part of
3918 the energies is plotted. With *weights* (e.g. from *band_weights*:
3919 an orbital, sublattice or spin projection), a "fat band" plot: the
3920 bands are colored by the weights (with a colorbar), or drawn with
3921 dots sized by their absolute value.
3923 Parameters
3924 ----------
3925 node_labels : List of strings, default None
3926 Labels of the
3927 high-symmetry points passed to *k_path*.
3928 fs : Positive number, default 20
3929 Fontsize.
3930 lw : Positive number, default 2
3931 Linewidth.
3932 ms : Positive number, default 0
3933 Marker size.
3934 c
3935 Default value 'b'. Line color.
3936 lims : List, default None
3937 Energy plot limits.
3938 figsize : Tuple, default None
3939 Figure size.
3940 weights : Real array, shape (nk, norb), default None
3941 One weight per k-point and band, e.g. from *band_weights*.
3942 style : String, default 'color'
3943 'color': the bands
3944 colored by *weights*; 'size': dots of size proportional to
3945 :math:`|w|` (*ms*, or 8 if *ms* is 0, for the largest) on thin
3946 lines of color *c*.
3947 cmap
3948 Default value None ('viridis', or 'RdBu_r', centred on
3949 zero, if some weights are negative). Colormap of *style* 'color'.
3951 Returns
3952 -------
3953 fig : Figure
3954 '''
3955 error_handling.empty_ndarray(self.en, 'get_bands or k_path')
3956 error_handling.positive_real(fs, 'fs')
3957 error_handling.positive_real(lw, 'lw')
3958 error_handling.lims(lims)
3959 error_handling.tuple_2elem(figsize, 'figsize')
3960 fig, ax = plt.subplots(figsize=figsize)
3961 if weights is None:
3962 for n in range(self.norb):
3963 ax.plot(self.ks_dist, self.en[:, n].real, c=c, lw=lw, marker='o', ms=ms)
3964 else:
3965 self._plot_fat_bands(fig, ax, weights, style, cmap, c, lw, ms)
3966 for node in self.nodes:
3967 ax.axvline(node, color='k', lw=0.5)
3968 ax.set_xlim(self.ks_dist[0], self.ks_dist[-1])
3969 if lims is not None:
3970 ax.set_ylim(lims)
3971 if node_labels is not None:
3972 error_handling.ndarray(np.array(node_labels), 'node_labels', len(self.nodes))
3973 ax.set_xticks(self.nodes)
3974 ax.set_xticklabels(node_labels, fontsize=fs)
3975 ax.set_ylabel('$E$', fontsize=fs)
3976 for label in ax.yaxis.get_majorticklabels():
3977 label.set_fontsize(fs)
3978 fig.set_layout_engine('tight')
3979 plt.draw()
3980 return fig
3982 def _plot_fat_bands(self, fig, ax, weights, style, cmap, c, lw, ms) -> None:
3983 '''
3984 Private method. The bands of *plot_bands* colored (or with dots
3985 sized) by *weights*.
3986 '''
3987 error_handling.band_weights(weights, self.en.shape)
3988 error_handling.weight_style(style)
3989 weights = np.asarray(weights, dtype='f8')
3990 signed = weights.min() < 0
3991 vmax = np.abs(weights).max() or 1.
3992 if style == 'size':
3993 size = (ms or 8.) ** 2
3994 for n in range(self.norb):
3995 ax.plot(self.ks_dist, self.en[:, n].real, c=c, lw=lw / 4)
3996 ax.scatter(self.ks_dist, self.en[:, n].real, s=size * np.abs(weights[:, n]) / vmax,
3997 c=c, lw=0)
3998 return
3999 norm = Normalize(-vmax if signed else 0., vmax)
4000 cmap = cmap or ('RdBu_r' if signed else 'viridis')
4001 for n in range(self.norb):
4002 pts = np.stack([self.ks_dist, self.en[:, n].real], axis=-1)
4003 lines = LineCollection(list(np.stack([pts[:-1], pts[1:]], axis=1)), cmap=cmap, norm=norm, lw=lw)
4004 lines.set_array((weights[:-1, n] + weights[1:, n]) / 2)
4005 ax.add_collection(lines)
4006 ax.autoscale_view()
4007 fig.colorbar(lines, ax=ax)
4009 def show(self) -> None:
4010 '''
4011 Emulate Matplotlib method plt.show().
4012 '''
4013 plt.show()
4016def ribbon(
4017 lat: Lattice,
4018 list_hop: list[dict],
4019 width: int,
4020 direction: int = 1,
4021 onsite: dict | None = None,
4022 spin: bool = False,
4023) -> KSpace:
4024 r'''
4025 Cut a ribbon out of a 2D periodic model: periodic along one primitive
4026 vector, finite (open boundary, *width* unit cells) along the other.
4027 This is the standard way to see edge states in a band structure (e.g.
4028 the zero-energy edge band of a zigzag graphene ribbon, or the helical
4029 edge states of a Kane-Mele ribbon). Cut out of a 3D model, it gives a
4030 slab: periodic along the two other primitive vectors, the way to see
4031 surface states (e.g. the Fermi arcs of a Weyl semimetal).
4033 Parameters
4034 ----------
4035 lat
4036 **Lattice** class instance (2D or 3D, i.e. two or three primitive
4037 vectors). Only *unit_cell* and *prim_vec* are used.
4038 list_hop
4039 List of dictionaries, in the same format passed to
4040 *KSpace.set_hopping* -- the hoppings of the periodic (2D) model
4041 that the ribbon is cut from.
4042 width : Positive integer
4043 Number of unit cells across the ribbon.
4044 direction : 0 or 1 (or 2, in 3D), default 1
4045 Which primitive
4046 vector (``lat.prim_vec[direction]``) becomes finite; the others stay
4047 periodic.
4048 onsite : Dictionary, default None
4049 Onsite energies, in the
4050 same format passed to *KSpace.set_onsite* -- applied identically
4051 on every row of the ribbon.
4052 spin : Boolean, default False
4053 See *KSpace*.
4055 Returns
4056 -------
4057 rib
4058 **KSpace** instance, 1D-periodic (2D-periodic for a slab), with
4059 ``width * len(lat.unit_cell)`` sites (each site of *lat*,
4060 repeated once per row across the ribbon; row *w*'s copy of site
4061 *i* is orbital ``w*len(lat.unit_cell) + i``).
4063 Examples
4064 --------
4065 Example usage::
4067 # zigzag graphene ribbon, 20 unit cells wide
4068 list_hop = [{'i': 0, 'j': 1, 'R': (0, 0), 't': 1.},
4069 {'i': 0, 'j': 1, 'R': (-1, 0), 't': 1.},
4070 {'i': 0, 'j': 1, 'R': (0, -1), 't': 1.}]
4071 rib = ribbon(lat, list_hop, width=20)
4072 '''
4073 error_handling.lat(lat)
4074 error_handling.dim_min(len(lat.prim_vec), 2)
4075 error_handling.positive_int(width, 'width')
4076 error_handling.direction(direction, len(lat.prim_vec))
4077 periodic = [d for d in range(len(lat.prim_vec)) if d != direction]
4078 n_sites = len(lat.unit_cell)
4079 a_dir = np.array(lat.prim_vec[direction])
4080 new_unit_cell = []
4081 for w in range(width):
4082 for dic in lat.unit_cell:
4083 r0 = np.array(dic['r0']) + w*a_dir
4084 new_unit_cell.append({'tag': dic['tag'], 'r0': tuple(float(c) for c in r0)})
4085 new_lat = Lattice(unit_cell=new_unit_cell, prim_vec=[lat.prim_vec[d] for d in periodic])
4086 rib = KSpace(new_lat, spin=spin)
4087 new_list_hop = []
4088 for dic in list_hop:
4089 w2_shift = dic['R'][direction]
4090 n_periodic = tuple(dic['R'][d] for d in periodic)
4091 for w in range(width):
4092 w2 = w + w2_shift
4093 if 0 <= w2 < width:
4094 new_list_hop.append({'i': w*n_sites + dic['i'],
4095 'j': w2*n_sites + dic['j'],
4096 'R': n_periodic,
4097 't': dic['t']})
4098 rib.set_hopping(new_list_hop)
4099 if onsite is not None:
4100 rib.set_onsite(onsite)
4101 return rib
4104def magnetic_supercell(
4105 lat: Lattice,
4106 list_hop: list[dict],
4107 p: int,
4108 q: int,
4109 onsite: dict | None = None,
4110 spin: bool = False,
4111) -> KSpace:
4112 r'''
4113 Thread a uniform magnetic field through a 2D periodic model, with a
4114 rational flux :math:`\alpha = p/q` flux quanta per unit cell, and return
4115 the Bloch Hamiltonian of its magnetic unit cell: *q* cells along
4116 :math:`\mathbf{a}_1`. Each hopping gets its Peierls phase (see
4117 *System.set_peierls_phase*)
4119 .. math::
4121 \phi_{ij} = \frac{2\pi}{\Phi_0}\int_{\mathbf{r}_i}^{\mathbf{r}_j}\mathbf{A}\cdot d\mathbf{l}
4122 = 2\pi\alpha\,\frac{s_{1,i}+s_{1,j}}{2}\,(s_{2,j}-s_{2,i})
4124 in the gauge :math:`\mathbf{A} = \alpha\Phi_0\, s_1\nabla s_2`, with
4125 :math:`(s_1, s_2)` the fractional coordinates of a point along
4126 :math:`(\mathbf{a}_1, \mathbf{a}_2)`, periodic along :math:`\mathbf{a}_2`.
4127 A translation by :math:`q\mathbf{a}_1` is a symmetry only up to the
4128 gauge transformation :math:`e^{-2\pi i p s_2}`, which the hoppings
4129 crossing the magnetic cell boundary carry. The bands of the result are
4130 the Hofstadter subbands, and their *chern_number* the TKNN integers.
4132 Parameters
4133 ----------
4134 lat : **Lattice** class instance, 2D (two primitive vectors in the plane)
4135 list_hop
4136 List of dictionaries, the hoppings of the field-free
4137 model, in the format of *KSpace.set_hopping*.
4138 p : Integer
4139 Numerator of the flux per unit cell.
4140 q : Positive integer
4141 Denominator of the flux per unit cell: the
4142 magnetic cell has *q* unit cells.
4143 onsite : Dictionary, default None
4144 Onsite energies, in the
4145 format of *KSpace.set_onsite*, applied to every copy.
4146 spin : Boolean, default False
4147 See *KSpace*.
4149 Returns
4150 -------
4151 mag
4152 **KSpace** instance, with ``q * len(lat.unit_cell)``
4153 sites: copy *m* of site *i* is site ``m*len(lat.unit_cell) + i``,
4154 displaced by :math:`m\mathbf{a}_1`.
4156 Examples
4157 --------
4158 Example usage::
4160 # square lattice, flux 1/3 per plaquette: three Hofstadter bands
4161 mag = magnetic_supercell(lattices.square(), list_hop, p=1, q=3)
4162 chern = [mag.chern_number(bands=[n]) for n in range(3)]
4163 '''
4164 error_handling.lat(lat)
4165 error_handling.planar_cell(lat.prim_vec)
4166 error_handling.set_hopping_kspace(list_hop, len(lat.unit_cell), 2, spin)
4167 error_handling.integer(p, 'p')
4168 error_handling.positive_int(q, 'q')
4169 alpha = p / q
4170 n_sites = len(lat.unit_cell)
4171 a = np.array(lat.prim_vec, dtype='f8')
4172 frac = np.linalg.solve(a.T, np.array([dic['r0'] for dic in lat.unit_cell], dtype='f8').T).T
4173 new_unit_cell = []
4174 for m in range(q):
4175 for dic in lat.unit_cell:
4176 r0 = np.array(dic['r0'], dtype='f8') + m * a[0]
4177 new_unit_cell.append({'tag': dic['tag'], 'r0': (float(r0[0]), float(r0[1]))})
4178 new_lat = Lattice(unit_cell=new_unit_cell, prim_vec=[tuple(q * a[0]), lat.prim_vec[1]])
4179 mag = KSpace(new_lat, spin=spin)
4180 new_list_hop = []
4181 for dic in list_hop:
4182 n1, n2 = dic['R']
4183 for m in range(q):
4184 s1_i, s2_i = m + frac[dic['i'], 0], frac[dic['i'], 1]
4185 s1_j, s2_j = m + n1 + frac[dic['j'], 0], n2 + frac[dic['j'], 1]
4186 big, m_j = divmod(m + n1, q)
4187 phase = 2 * PI * alpha * 0.5 * (s1_i + s1_j) * (s2_j - s2_i) - 2 * PI * p * big * s2_j
4188 new_list_hop.append({'i': m * n_sites + dic['i'], 'j': m_j * n_sites + dic['j'],
4189 'R': (big, n2), 't': dic['t'] * np.exp(1j * phase)})
4190 mag.set_hopping(new_list_hop)
4191 if onsite is not None:
4192 mag.set_onsite(onsite)
4193 return mag
4196#################################
4197# K-PATHS AND CONSTANT-ENERGY CONTOURS
4198#################################
4201def _response_weights(
4202 en: NDArray[np.float64], mu: float, temperature: float, response: str,
4203) -> NDArray[np.float64]:
4204 r'''
4205 Private. The weights :math:`P_{nm}` of the Kubo pair sum of
4206 *KSpace._kubo*, shape (nk, n, n), from the bands *en* (nk, n). With
4207 :math:`x_n = (E_n-\mu)/T`, the conductivities are
4208 :math:`\sum_n w_n\Omega_n`, i.e. :math:`P_{nm} = -(w_n - w_m)`, with
4209 :math:`w = f` ('hall'), the entropy per state
4210 :math:`s = -f\ln f-(1-f)\ln(1-f)` ('nernst') and
4211 :math:`c_2/T = T\int_x^\infty y^2(-f'(y))\,dy` ('thermal'). The orbital
4212 magnetization ('magnetization'),
4213 :math:`\sum_n[f_n\,\mathrm{Im}\langle\partial_xu_n|(H-E_n)|\partial_yu_n\rangle + g_n\Omega_n]`
4214 with :math:`g = T\ln(1+e^{-x})` (:math:`\max(\mu-E, 0)` at :math:`T = 0`),
4215 has :math:`P_{nm} = \frac12(f_n+f_m)(E_m-E_n) - (g_n-g_m)`.
4216 '''
4217 if response == 'magnetization':
4218 f = occupation.fermi_dirac(en, mu, temperature)
4219 if temperature == 0:
4220 g = np.maximum(mu - en, 0.)
4221 else:
4222 g = temperature * np.logaddexp(0., -(en - mu) / temperature)
4223 return 0.5 * (f[:, :, None] + f[:, None, :]) * (en[:, None, :] - en[:, :, None]) \
4224 - (g[:, :, None] - g[:, None, :])
4225 if response == 'hall':
4226 w = occupation.fermi_dirac(en, mu, temperature)
4227 else:
4228 a = np.abs(en - mu) / temperature
4229 if response == 'nernst':
4230 w = np.log1p(np.exp(-a)) + a * expit(-a)
4231 else:
4232 # int_a^inf y^2 (-f') dy = a^2 f(a) + 2 [a ln(1 + e^-a) - Li2(-e^-a)]
4233 # for a >= 0, and pi^2/3 minus it below the Fermi level
4234 # (Li2(z) = spence(1 - z) in scipy's convention)
4235 c2 = a ** 2 * expit(-a) + 2 * (a * np.log1p(np.exp(-a)) - spence(1 + np.exp(-a)))
4236 w = temperature * np.where(en >= mu, c2, PI ** 2 / 3 - c2)
4237 return -(w[:, :, None] - w[:, None, :])
4240def _path(
4241 points: list[ArrayLike], nk: int,
4242) -> tuple[NDArray[np.float64], NDArray[np.float64], list[float]]:
4243 '''
4244 Private function. The k-points of a path through *points* (*nk* per
4245 segment, the last point included), their cumulative distance, and the
4246 positions of the nodes along it.
4247 '''
4248 pts = np.atleast_2d(np.asarray(points, dtype='f8'))
4249 ks = np.concatenate([np.linspace(pts[i], pts[i+1], nk, endpoint=False)
4250 for i in range(len(pts) - 1)] + [pts[-1:]])
4251 dist = np.concatenate([[0.], np.cumsum(np.linalg.norm(np.diff(ks, axis=0), axis=1))])
4252 return ks, dist, dist[::nk][:len(pts)-1].tolist() + [dist[-1]]
4255def _marching_table(dim: int) -> dict[int, list[list[tuple[int, int]]]]:
4256 '''
4257 Private function. For each pattern of the vertices of a triangle
4258 (*dim* = 2) or tetrahedron (*dim* = 3) above a level (bit m set if
4259 vertex m is above), the pieces of the level set, as lists of the edges
4260 (a, b) that they cross: one segment in 2D; one triangle, or two
4261 triangles forming a quadrilateral when two vertices are above, in 3D.
4262 '''
4263 table = {}
4264 for pattern in range(1, 2 ** (dim + 1) - 1):
4265 up = [m for m in range(dim + 1) if pattern >> m & 1]
4266 down = [m for m in range(dim + 1) if not pattern >> m & 1]
4267 if dim == 3 and len(up) == 2:
4268 (a, b), (c, d) = up, down
4269 table[pattern] = [[(a, c), (a, d), (b, d)], [(a, c), (b, d), (b, c)]]
4270 else:
4271 table[pattern] = [[(a, b) for a in up for b in down]]
4272 return table
4275def _level_set(values: NDArray[np.float64], level: float) -> NDArray[np.float64]:
4276 '''
4277 Private function. The level set ``values == level`` of a function on a
4278 periodic 2D or 3D mesh (*values* of shape (n1, n2) or (n1, n2, n3)),
4279 linear in each simplex of *dos._mesh_simplices*: segments (shape
4280 (n, 2, 2)) or triangles (shape (n, 3, 3)), in fractional coordinates
4281 (not folded: between 0 and 1 + 1/n_d).
4282 '''
4283 shape, dim = values.shape, values.ndim
4284 offsets = dos._mesh_simplices(dim)
4285 vals = dos._simplex_values(values - level, offsets).reshape(dim + 1, -1)
4286 pattern = sum((vals[m] > 0).astype(int) << m for m in range(dim + 1))
4287 pieces = [np.empty((0, dim, dim))]
4288 for pat, polygons in _marching_table(dim).items():
4289 sel = np.nonzero(pattern == pat)[0]
4290 if not len(sel):
4291 continue
4292 simplex, cell = np.divmod(sel, int(np.prod(shape)))
4293 corner = np.stack(np.unravel_index(cell, shape), axis=-1)
4294 pos = corner[None] + offsets[simplex].transpose(1, 0, 2)
4295 v = vals[:, sel]
4296 for polygon in polygons:
4297 # the crossing on edge (a, b): v[a] > 0 >= v[b] or the reverse
4298 pieces.append(np.stack([pos[a] + (v[a] / (v[a] - v[b]))[:, None] * (pos[b] - pos[a])
4299 for a, b in polygon], axis=1))
4300 return np.concatenate(pieces) / np.array(shape)
4303def _zone_polygon(rec: NDArray[np.float64]) -> NDArray[np.float64]:
4304 '''
4305 Private function. Vertices, counterclockwise, of the first Brillouin
4306 zone (Wigner-Seitz cell) of the 2D reciprocal lattice with basis *rec*
4307 (rows): the points at least as close to the origin as to any other
4308 reciprocal lattice vector, where two of the bisectors meet.
4309 '''
4310 gs = np.array([n @ rec for n in product(range(-2, 3), repeat=2) if any(n)])
4311 half = np.sum(gs ** 2, axis=1) / 2
4312 tol = 1e-9 * half.max()
4313 found = []
4314 for i in range(len(gs)):
4315 for j in range(i + 1, len(gs)):
4316 mat = gs[[i, j]]
4317 if abs(np.linalg.det(mat)) > tol:
4318 k = np.linalg.solve(mat, half[[i, j]])
4319 if np.all(gs @ k <= half + tol):
4320 found.append(k)
4321 verts = np.unique(np.round(np.array(found), 9), axis=0)
4322 return verts[np.argsort(np.arctan2(verts[:, 1], verts[:, 0]))]
4325def _gauss_reduce(u: NDArray[np.float64], v: NDArray[np.float64]):
4326 '''
4327 Private function. Gauss-reduced basis of a 2D lattice: the two
4328 shortest independent vectors, |u| <= |v|, at an obtuse (or right)
4329 angle, so that the edges of the Wigner-Seitz cell are normal to
4330 u, v and u + v.
4331 '''
4332 while True:
4333 if u @ u > v @ v:
4334 u, v = v, u
4335 m = np.round(u @ v / (u @ u))
4336 if m == 0:
4337 break
4338 v = v - m * u
4339 return (u, v - u) if u @ v > 0 else (u, v)
4342def _successive_minima(vecs: NDArray[np.float64]):
4343 '''
4344 Private function. The shortest lattice vectors of a 3D lattice with
4345 basis *vecs* (rows): the three successive minima, the shells of
4346 lattice vectors sorted by length (list of arrays), and their lengths.
4347 '''
4348 pts = np.array([n @ vecs for n in product(range(-3, 4), repeat=3) if any(n)])
4349 lengths = np.linalg.norm(pts, axis=1)
4350 # by length, then (among equal lengths) the most positive first: x, then y, then z
4351 order = np.lexsort((-pts[:, 2], -pts[:, 1], -pts[:, 0], np.round(lengths / lengths.min(), 6)))
4352 pts, lengths = pts[order], lengths[order]
4353 minima: list = []
4354 for p in pts:
4355 if np.linalg.matrix_rank(np.array(minima + [p]), tol=1e-8 * lengths[0]) == len(minima) + 1:
4356 minima.append(p)
4357 if len(minima) == 3:
4358 break
4359 bounds = np.nonzero(np.diff(lengths) > 1e-6 * lengths[0])[0] + 1
4360 return np.array(minima), np.split(pts, bounds), [s[0] for s in np.split(lengths, bounds)]
4363def _orthogonal_axes(shell: NDArray[np.float64]) -> NDArray[np.float64]:
4364 '''
4365 Private function. Three mutually orthogonal unit vectors among the
4366 vectors of *shell* (the conventional cubic axes).
4367 '''
4368 x = shell[0] / np.linalg.norm(shell[0])
4369 y = next(s for s in shell if abs(s @ x) < 1e-6 * np.linalg.norm(s))
4370 y = y / np.linalg.norm(y)
4371 return np.array([x, y, np.cross(x, y)])
4374def cubic_path(
4375 points: dict[str, ArrayLike], names: list[str], a: float = 1., axes: ArrayLike | None = None,
4376) -> tuple[list[NDArray[np.float64]], list[str]]:
4377 r'''
4378 Build a path through k-points given in the conventional coordinates of
4379 a cubic crystal: in units of :math:`2\pi/a` along the cubic axes, with
4380 :math:`a` the side of the conventional cube, as tables of high-symmetry
4381 points list them (e.g. :math:`X = (0, 1, 0)` and
4382 :math:`L = (\tfrac12, \tfrac12, \tfrac12)` for fcc). *high_symmetry_path*
4383 builds the standard simple cubic, fcc and bcc paths with it; use it
4384 for another path, or other points, of the same lattices.
4386 Parameters
4387 ----------
4388 points : Dictionary
4389 key: label, val: tuple of 3 real numbers,
4390 the point in units of :math:`2\pi/a`.
4391 names
4392 List of at least two labels (keys of *points*), the
4393 nodes of the path in order.
4394 a : Positive real, default 1
4395 Side of the conventional cube.
4396 axes : Real array, shape (3, 3)
4397 Default value None: the
4398 Cartesian axes. Rows: the three orthonormal cubic axes, in the
4399 coordinates of *KSpace.get_ham* (for a crystal whose cube is rotated).
4401 Returns
4402 -------
4403 points : List of real ndarrays, the k-points, for *KSpace.k_path*
4404 labels : List of strings, *names*, for *node_labels*
4406 Examples
4407 --------
4408 Example usage::
4410 # fcc, along Gamma - X - U and K - Gamma - L
4411 pts = {'G': (0, 0, 0), 'X': (0, 1, 0), 'U': (.25, 1, .25),
4412 'K': (.75, .75, 0), 'L': (.5, .5, .5)}
4413 points, labels = cubic_path(pts, ['G', 'X', 'U'], a=3.6)
4414 '''
4415 error_handling.cubic_path(points, names, axes)
4416 error_handling.positive_real(a, 'a')
4417 axes = np.eye(3) if axes is None else np.asarray(axes, dtype='f8')
4418 return [2 * PI / a * np.asarray(points[n], dtype='f8') @ axes for n in names], list(names)
4421def _bravais_3d(vecs: NDArray[np.float64]):
4422 '''
4423 Private function. The points and labels of the high-symmetry path of a
4424 3D lattice (basis *vecs*, rows), or (None, None) if it is not one of
4425 the simple cubic, fcc, bcc, simple tetragonal or hexagonal lattices.
4426 '''
4427 minima, shells, radii = _successive_minima(vecs)
4428 lengths = np.linalg.norm(minima, axis=1)
4429 pairs = [(0, 1), (0, 2), (1, 2)]
4430 cosine = {p: abs(minima[p[0]] @ minima[p[1]]) / (lengths[p[0]] * lengths[p[1]]) for p in pairs}
4431 equal = {p: abs(lengths[p[0]] - lengths[p[1]]) < 1e-6 * lengths[2] for p in pairs}
4432 right = {p: cosine[p] < 1e-6 for p in pairs}
4433 G = r'$\Gamma$'
4434 if all(right.values()):
4435 if all(equal.values()): # simple cubic
4436 return cubic_path({G: (0, 0, 0), 'X': (0, .5, 0), 'M': (.5, .5, 0), 'R': (.5, .5, .5)},
4437 [G, 'X', 'M', G, 'R', 'X'], lengths[0], minima / lengths[:, None])
4438 if sum(equal.values()) != 1: # orthorhombic
4439 return None, None
4440 # simple tetragonal: a, a, c
4441 (i, j), = [p for p in pairs if equal[p]]
4442 k = 3 - i - j
4443 a, c = lengths[i], lengths[k]
4444 axes = np.array([minima[i] / a, minima[j] / a, minima[k] / c])
4445 pts = {G: (0, 0, 0), 'X': (0, .5, 0), 'M': (.5, .5, 0), 'Z': (0, 0, .5),
4446 'R': (0, .5, .5), 'A': (.5, .5, .5)}
4447 names = [G, 'X', 'M', G, 'Z', 'R', 'A', 'Z']
4448 scale = 2 * PI * np.array([1 / a, 1 / a, 1 / c])
4449 return [(scale * pts[n]) @ axes for n in names], names
4450 for (i, j) in pairs:
4451 k = 3 - i - j
4452 if equal[(i, j)] and abs(cosine[(i, j)] - .5) < 1e-6 and \
4453 all(right[p] for p in pairs if p != (i, j)):
4454 # hexagonal: reciprocal vectors of (a1, a2, c), in-plane ones at 120 degrees
4455 rec = 2 * PI * np.linalg.inv(minima[[i, j, k]]).T
4456 b1, b2 = _gauss_reduce(rec[0], rec[1])
4457 M, K, A = b1 / 2, (2 * b1 + b2) / 3, rec[2] / 2
4458 hex_pts = {G: np.zeros(3), 'M': M, 'K': K, 'A': A, 'L': M + A, 'H': K + A}
4459 names = [G, 'M', 'K', G, 'A', 'L', 'H', 'A']
4460 return [hex_pts[n] for n in names], names
4461 if not all(equal.values()):
4462 return None, None
4463 if len(shells[0]) == 12: # fcc: 12 nearest neighbours at a/sqrt(2), then 6 at a
4464 return cubic_path({G: (0, 0, 0), 'X': (0, 1, 0), 'W': (.5, 1, 0), 'K': (.75, .75, 0),
4465 'L': (.5, .5, .5), 'U': (.25, 1, .25)},
4466 [G, 'X', 'W', 'K', G, 'L', 'U', 'W', 'L', 'K'],
4467 radii[1], _orthogonal_axes(shells[1]))
4468 if len(shells[0]) == 8 and all(abs(c - 1 / 3) < 1e-6 for c in cosine.values()):
4469 # bcc: 8 nearest neighbours at a sqrt(3)/2, then 6 at a
4470 return cubic_path({G: (0, 0, 0), 'H': (0, 0, 1), 'N': (.5, 0, .5), 'P': (.5, .5, .5)},
4471 [G, 'H', 'N', G, 'P', 'H'], radii[1], _orthogonal_axes(shells[1]))
4472 return None, None
4475def high_symmetry_path(lat: Lattice) -> tuple[list[NDArray[np.float64]], list[str]]:
4476 r'''
4477 Get a standard path through the high-symmetry points of the Brillouin
4478 zone of a lattice, for *KSpace.k_path* and the labels of
4479 *KSpace.plot_bands*. The Bravais lattice is recognized from
4480 *lat.prim_vec*, whatever the choice of primitive vectors and the
4481 orientation of the lattice:
4483 * 1D: :math:`-X, \Gamma, X` (the whole zone, :math:`X = \mathbf{b}/2`).
4484 * 2D (the five Bravais lattices): square :math:`\Gamma X M \Gamma`;
4485 rectangular :math:`\Gamma X S Y \Gamma`; hexagonal
4486 :math:`\Gamma M K \Gamma`; centred rectangular and oblique
4487 :math:`\Gamma X H C H_1 Y \Gamma`, along the edges of the
4488 zone: :math:`X, Y, C` are the midpoints of the edges normal to the
4489 reciprocal vectors :math:`\mathbf{b}_1, \mathbf{b}_2,
4490 \mathbf{b}_1+\mathbf{b}_2` of the reduced basis (with
4491 :math:`|\mathbf{b}_1| \le |\mathbf{b}_2|`), :math:`H` and
4492 :math:`H_1` the corners between them.
4493 * 3D: simple cubic :math:`\Gamma X M \Gamma R X`; fcc
4494 :math:`\Gamma X W K \Gamma L U W L K`; bcc :math:`\Gamma H N \Gamma P H`;
4495 simple tetragonal :math:`\Gamma X M \Gamma Z R A Z`; hexagonal
4496 :math:`\Gamma M K \Gamma A L H A`. These are the paths of Setyawan
4497 and Curtarolo (Comput. Mater. Sci. 49, 299 (2010)), up to their
4498 first discontinuity (a path here is one continuous line).
4500 The other 3D lattices (orthorhombic, monoclinic, body-centred
4501 tetragonal, rhombohedral, ...) raise a ValueError: give their points
4502 to *k_path* directly.
4504 Parameters
4505 ----------
4506 lat : **lattice** class instance (only *prim_vec* is used)
4508 Returns
4509 -------
4510 points
4511 List of real ndarrays, the k-points, in the
4512 coordinates of *KSpace.get_ham*.
4513 labels : List of strings (LaTeX), for *node_labels*
4515 Examples
4516 --------
4517 Example usage::
4519 points, labels = high_symmetry_path(lat)
4520 ks.k_path(points, nk=60)
4521 ks.plot_bands(node_labels=labels)
4522 '''
4523 ks = KSpace(lat)
4524 vecs = np.array(lat.prim_vec, dtype='f8') @ ks.k_basis # in the k coordinates
4525 G = r'$\Gamma$'
4526 if ks.dim == 1:
4527 x = PI / vecs[0]
4528 return [-x, np.zeros(1), x], ['$-X$', G, 'X']
4529 if ks.dim == 3:
4530 points, labels = _bravais_3d(vecs)
4531 error_handling.bravais_lattice(labels, 3)
4532 return points, labels
4533 rec = 2 * PI * np.linalg.inv(vecs).T
4534 b1, b2 = _gauss_reduce(rec[0], rec[1])
4535 l1, l2, l3 = (np.linalg.norm(b) for b in (b1, b2, b1 + b2))
4536 close = lambda a, b: abs(a - b) < 1e-6 * max(abs(a), abs(b))
4537 cosine = b1 @ b2 / (l1 * l2)
4538 if abs(cosine) < 1e-6:
4539 if close(l1, l2):
4540 return [np.zeros(2), b1 / 2, (b1 + b2) / 2, np.zeros(2)], [G, 'X', 'M', G]
4541 return [np.zeros(2), b1 / 2, (b1 + b2) / 2, b2 / 2, np.zeros(2)], [G, 'X', 'S', 'Y', G]
4542 if close(l1, l2) and close(cosine, -.5):
4543 return [np.zeros(2), b1 / 2, (2 * b1 + b2) / 3, np.zeros(2)], [G, 'M', 'K', G]
4544 # centred rectangular and oblique: around half of the hexagonal zone
4545 b3 = b1 + b2
4546 corner = lambda u, v: np.linalg.solve(np.array([u, v]), np.array([u @ u, v @ v]) / 2)
4547 return ([np.zeros(2), b1 / 2, corner(b1, b3), b3 / 2, corner(b2, b3), b2 / 2, np.zeros(2)],
4548 [G, 'X', 'H', 'C', '$H_1$', 'Y', G])