Coverage for tbkit/system.py: 100%
569 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +0000
1from __future__ import annotations
3from typing import Callable
5from numpy.typing import ArrayLike
7import numpy as np
8from numpy.typing import NDArray
9import scipy.sparse as sparse
10import scipy.linalg as LA
11import scipy.sparse.linalg as SLA
12from scipy.spatial import cKDTree
13import numpy.random as rand
14npc = np.char # `import numpy.char` needs NumPy >= 2.0
15from itertools import product
16from math import sin, cos
17import tbkit.error_handling as error_handling
18import tbkit.occupation as occupation
19import tbkit.values as values
20from tbkit.lattice import Lattice, COOR_DTYPE
23PI = np.pi
24ATOL = 1e-3
25HOP_DTYPE = [('n', 'u2'), ('i', 'u4'), ('j', 'u4'),
26 ('ang', 'f8'), ('tag', 'U2'), ('t', 'c16')]
29def _upper(ang: NDArray) -> NDArray:
30 '''
31 Private function.
32 Mask of the bond angles (in degrees, from np.arctan2) that point into the
33 half-plane :math:`[0, 180)`, within ATOL: a bond at -1e-14 degrees counts as
34 horizontal (0), one at 179.99999 as pointing the other way.
35 '''
36 return (ang >= -ATOL) & (ang < 180. - ATOL)
39class System():
40 r'''
41 Solve the Tight-Binding eigenvalue problem of a lattice defined
42 by the class **lattice**.
44 Parameters
45 ----------
46 lat : **lattice** class instance
47 periodic : Boolean, or tuple of booleans (one per primitive vector), default False
48 Periodic boundary conditions along the
49 primitive vectors :math:`\mathbf{a}_i`: the sample is then the
50 :math:`n_1\mathbf{a}_1\times n_2\mathbf{a}_2\,(\times n_3\mathbf{a}_3)`
51 box of *lat.get_lattice* wrapped into a torus (True), or a cylinder
52 (e.g. ``(True, False)``, periodic along :math:`\mathbf{a}_1` only).
53 *lat.get_lattice* must have been called; sites may be removed
54 afterwards (vacancies), but not moved. Every distance and bond
55 angle is then that of the shortest image, so the neighbour orders of
56 *set_hopping* include the bonds that wrap around, with their short
57 bond vector; a bond order reaching half-way around the torus is
58 refused (it would be ambiguous). A uniform magnetic field
59 (*set_magnetic_field*) and the local Chern marker need open
60 boundaries; *get_bott_index* is the Chern number of a torus.
62 Notes
63 -----
64 Up to *dense_max* sites (default 5000), the distances between all pairs
65 of sites are computed (and kept in *vec_hop*). Larger lattices only look
66 for the neighbours they need, with a k-d tree, up to the largest
67 neighbour order a call asks for: memory and time then grow linearly
68 with the number of sites (*vec_hop* stays empty).
69 '''
71 #: Largest lattice (number of sites) whose pairwise distances are all computed.
72 dense_max = 5000
74 def __init__(self, lat: Lattice, periodic: bool | tuple[bool, ...] = False) -> None:
75 error_handling.lat(lat)
76 self.periodic = error_handling.periodic(periodic, len(lat.prim_vec))
77 if any(self.periodic):
78 error_handling.torus_lattice(lat.n1)
79 self.lat = lat
80 self.sites = self.lat.sites # used to check if sites changes
81 self.coor_hop = np.array([], dtype=COOR_DTYPE)
82 self.vec_hop = np.array([], dtype=[('dis', 'f8'), ('ang', 'f8')]) # Hopping distances and angles
83 self.dist_uni = np.array([], 'f8') # Different hopping distances
84 self.store_hop: dict = {} # Store the relevant hoppings (dynamic programming)
85 self.hop = np.array([], dtype=HOP_DTYPE) # Hoppings to build-up the Hamiltonian
86 self.onsite = np.array([], 'c16') # Onsite energies
87 # Value functions (see tbkit.values): (directed bond keys, function)
88 # of the hoppings whose row of `hop` holds a factor 't', and (site
89 # indices, function) of the onsite terms added to `onsite`.
90 self._hop_values: list = []
91 self._onsite_values: list = []
92 self.ham = sparse.csr_matrix(([], ([], [])), shape=(self.lat.sites, self.lat.sites)) # Hamiltonian
93 self.en = np.array([], 'c16') # Eigenenergies
94 self.rn = np.array([], 'c16') # Right eigenvectors: H |rn> = en |rn>
95 self.ln = np.array([], 'c16') # Left eigenvectors: <ln| H = en <ln|
96 self.intensity = np.array([], 'f8') # Intensities (|rn|**2)
97 self.pola = np.array([], 'f8') # sublattices polarisation (|rn^{(S)}|**2)
98 self.ipr = np.array([], 'f8') # Inverse Participation Ratio
99 self.petermann = np.array([], 'f8') # Petermann factors
100 self.nmax = 0 # number of different hoppings
102 def clear_hopping(self) -> None:
103 '''
104 Clear structured array *hop*.
105 '''
106 self.hop = np.array([], dtype=HOP_DTYPE)
107 self._hop_values = []
109 def get_distances(self, n_orders: int = 1) -> None:
110 '''
111 Private method.
112 Get distances and angles of the edges. The angle is that of the
113 bond's projection on the (x, y) plane; in 3D space, a vertical bond
114 gets angle 0 pointing up (+z) and -180 pointing down.
116 Parameters
117 ----------
118 n_orders : Positive integer, default 1
119 Beyond *dense_max*
120 sites, the number of neighbour orders (distinct distances) needed.
121 '''
122 error_handling.sites(self.lat.sites)
123 if self.lat.sites > self.dense_max:
124 self._sparse_distances(n_orders)
125 return
126 self._pairs: NDArray | None = None
127 if any(self.periodic):
128 # shortest images, a block of rows at a time to bound the memory
129 coords = self._coords()
130 n = len(coords)
131 d = np.empty((n, n, coords.shape[1]))
132 step = max(1, 1000000 // n)
133 for r0 in range(0, n, step):
134 d[r0:r0 + step] = self._min_image(coords[None, :, :] - coords[r0:r0 + step, None, :])
135 dif = tuple(d[..., c] for c in range(d.shape[-1]))
136 else:
137 dif = self.lat.distances()
138 dist = np.sqrt(sum(d ** 2 for d in dif))
139 ang = (180 / PI * np.arctan2(dif[1], dif[0]))
140 if len(dif) == 3:
141 # 'ang' is the angle of the bond's projection on the (x, y) plane.
142 # A vertical bond has no such angle: orient it by dz instead, so
143 # that the sign of 'ang' still tells the two directions apart.
144 vertical = np.hypot(dif[0], dif[1]) < 1e-6 * np.maximum(dist, 1.)
145 ang[vertical & (dif[2] > 0)] = 0.
146 ang[vertical & (dif[2] <= 0)] = -180.
147 self.vec_hop = np.zeros(dist.shape, dtype=[('dis', 'f8'), ('ang', 'f8')])
148 self.vec_hop['dis'] = dist
149 self.vec_hop['ang'] = ang
150 self.dist_uni = np.unique(self.vec_hop['dis'].round(4))
152 def _coords(self) -> NDArray[np.float64]:
153 '''
154 Private method. Site coordinates, shape (sites, space_dim).
155 '''
156 return np.stack([self.lat.coor[f] for f in ('x', 'y', 'z')[:self.lat.space_dim]], axis=1)
158 def _box(self) -> NDArray[np.float64]:
159 '''
160 Private method. The vectors n_i a_i of the sample of *get_lattice*, shape (dim, space_dim).
161 '''
162 counts = np.array((self.lat.n1, self.lat.n2, self.lat.n3)[:len(self.lat.prim_vec)], dtype='f8')
163 return counts[:, None] * np.array(self.lat.prim_vec, dtype='f8')
165 def _torus(self) -> NDArray[np.float64]:
166 '''
167 Private method. The periodic vectors n_i a_i, shape (n_periodic, space_dim).
168 '''
169 return self._box()[list(self.periodic)]
171 def _min_image(self, d: NDArray[np.float64]) -> NDArray[np.float64]:
172 '''
173 Private method. Shortest images of displacements *d* (shape (..., space_dim))
174 under the periodic translations: the fractional coordinates along the
175 periodic vectors are rounded, then the neighbouring images compared
176 (rounding alone can miss the shortest one in a skewed cell). Only
177 displacements longer than half the shortest period need that
178 comparison: any other image of a shorter one is longer than it.
179 '''
180 if not any(self.periodic):
181 return d
182 torus = self._torus()
183 shape = d.shape
184 d = d.reshape(-1, shape[-1]) # 2D, for fast matrix products
185 frac = d @ np.linalg.pinv(self._box())[:, list(self.periodic)]
186 d = d - np.rint(frac) @ torus
187 sq = np.einsum('ij,ij->i', d, d)
188 far = np.flatnonzero(sq > (self._shortest_period() / 2) ** 2)
189 start = d[far]
190 best, best_sq = start.copy(), sq[far]
191 for m in product((-1, 0, 1), repeat=len(torus)):
192 if any(m):
193 cand = start + np.array(m, dtype='f8') @ torus
194 cand_sq = np.einsum('ij,ij->i', cand, cand)
195 shorter = cand_sq < best_sq - 1e-12
196 best[shorter], best_sq[shorter] = cand[shorter], cand_sq[shorter]
197 d[far] = best
198 return d.reshape(shape)
200 def _bond_vectors(self, i: NDArray, j: NDArray) -> NDArray[np.float64]:
201 '''
202 Private method. Bond vectors from sites *i* to sites *j*, shape
203 (len(i), space_dim): the shortest image on a torus.
204 '''
205 coords = self._coords()
206 return self._min_image(coords[j] - coords[i])
208 def _wrapped(self, i: NDArray, j: NDArray) -> NDArray:
209 '''
210 Private method. Mask of the bonds i -> j that wrap around the torus
211 (their shortest image is not the straight segment between the sites).
212 '''
213 coords = self._coords()
214 return np.linalg.norm(self._bond_vectors(i, j) - (coords[j] - coords[i]), axis=1) > ATOL
216 def _shortest_period(self) -> float:
217 '''
218 Private method. Length of the shortest nonzero periodic translation
219 (infinite with open boundaries).
220 '''
221 if not any(self.periodic):
222 return np.inf
223 torus = self._torus()
224 return min(np.linalg.norm(np.array(m, dtype='f8') @ torus)
225 for m in product(range(-2, 3), repeat=len(torus)) if any(m))
227 def _near_pairs(self, tree: cKDTree, coords: NDArray[np.float64], radius: float) -> NDArray:
228 '''
229 Private method. Pairs i < j of sites closer than *radius* (shortest
230 image on a torus, from the images of the sites under the periodic
231 translations of -1, 0, +1 periods), sorted.
232 '''
233 if not any(self.periodic):
234 return tree.query_pairs(radius, output_type='ndarray')
235 torus = self._torus()
236 shifts = np.array(list(product((-1, 0, 1), repeat=len(torus))), dtype='f8') @ torus
237 images = (coords[None, :, :] + shifts[:, None, :]).reshape(-1, coords.shape[1])
238 found = tree.sparse_distance_matrix(cKDTree(images), radius, output_type='ndarray')
239 i, j = found['i'], found['j'] % len(coords)
240 keep = i < j
241 return np.unique(np.column_stack([i[keep], j[keep]]), axis=0).reshape(-1, 2)
243 def _sparse_distances(self, n_orders: int) -> None:
244 '''
245 Private method. The pairs of sites up to the *n_orders*-th distinct
246 distance (or all of them), from a k-d tree: the search radius grows
247 until that many distinct distances are found, every pair within it
248 being found.
249 '''
250 coords = self._coords()
251 tree = cKDTree(coords)
252 near = tree.query(coords, k=2)[0][:, 1]
253 radius = max(np.min(near[near > 1e-9]) if np.any(near > 1e-9) else 1., 1e-6)
254 span = np.max(np.ptp(coords, axis=0)) * np.sqrt(coords.shape[1]) + radius
255 while True:
256 pairs = self._near_pairs(tree, coords, radius * (1. + 1e-6) + ATOL)
257 dis = np.linalg.norm(self._bond_vectors(pairs[:, 0], pairs[:, 1]), axis=1)
258 uni = np.unique(dis.round(4))
259 if len(uni) >= n_orders + 1 or radius > span:
260 break
261 radius *= 1.5
262 self.vec_hop = np.array([], dtype=[('dis', 'f8'), ('ang', 'f8')])
263 self._pairs, self._pair_dis = pairs, dis
264 # the (n+1)-th distance may be incomplete: keep the first n_orders
265 self.dist_uni = np.concatenate([[0.], uni[:n_orders] if radius <= span else uni])
267 def _angles(self, i: NDArray, j: NDArray) -> NDArray[np.float64]:
268 '''
269 Private method. Angles (degrees) of the bonds from sites i to sites j,
270 with the vertical-bond rule of *get_distances*.
271 '''
272 d = self._bond_vectors(i, j)
273 ang = 180 / PI * np.arctan2(d[:, 1], d[:, 0])
274 if d.shape[1] == 3:
275 vertical = np.hypot(d[:, 0], d[:, 1]) < 1e-6 * np.maximum(np.linalg.norm(d, axis=1), 1.)
276 ang[vertical & (d[:, 2] > 0)] = 0.
277 ang[vertical & (d[:, 2] <= 0)] = -180.
278 return ang
280 def print_distances(self, n: int = 1) -> None:
281 r'''
282 Print distances and positive angles (in degrees) :math:`\phi_+\in[0, 180)`
283 of the nth shortest edges. Negative angles are given by:
284 :math:`\phi_-= \phi_+-180` and :math:`\phi_-\in[-180, 0)`.
286 Parameters
287 ----------
288 n : Positive integer
289 Number of shortest edges.
290 '''
291 error_handling.sites(self.lat.sites)
292 error_handling.positive_int(n, 'n')
293 self.get_distances(n)
294 self.nmax = len(self.dist_uni) - 1
295 error_handling.positive_int_lim(n, 'n', self.nmax)
296 print('\n{} different distances between sites:'.format(self.nmax))
297 print('\nDistances between sites:')
298 for i, d in enumerate(self.dist_uni[1: n+1]):
299 if i == 0:
300 hop_name = 'st'
301 elif i == 1:
302 hop_name = 'nd'
303 elif i == 2:
304 hop_name = 'rd'
305 else:
306 hop_name = 'th'
307 print('{}{} hopping, length: {:.3f}'.format(i+1, hop_name, d))
308 print('\twith positive angles:')
309 print('\t', np.unique(self.get_bonds(i + 1)[2].round(4)))
311 def set_onsite(self, dict_onsite: dict[str, complex]) -> None:
312 '''
313 Set onsite energies.
315 Parameters
316 ----------
317 dict_onsite : Dictionary
318 key: sublattice tag, val: onsite
319 energy, or a value function ``onsite(site, **params)``, evaluated
320 by ``get_ham(**params)``. *site* is a structured array of the
321 sites of that tag (fields 'x', 'y', ('z'), 'tag', 'index'); the
322 function returns one energy, or one per site. *set_onsite_dis*
323 adds to it.
325 Examples
326 --------
327 Example usage::
329 # Line-Centered Square lattice
330 sys.set_onsite({'a': -1j, 'b': -2j})
331 # a gate voltage V across the sample, set by get_ham(V=...)
332 sys.set_onsite({'a': lambda site, V: V * site['x'] / 10, 'b': 0.})
333 sys.get_ham(V=0.5)
334 '''
335 error_handling.sites(self.lat.sites)
336 error_handling.set_onsite(dict_onsite, self.lat.tags)
337 self.onsite = np.zeros(self.lat.sites, 'c16')
338 self._onsite_values = []
339 for tag, on in dict_onsite.items():
340 if callable(on):
341 self._onsite_values.append((np.flatnonzero(self.lat.coor['tag'] == tag), on))
342 else:
343 self.onsite[self.lat.coor['tag'] ==tag] = on
345 def get_bonds(self, n: int) -> tuple[NDArray, NDArray, NDArray]:
346 '''
347 Private method.
348 Get every bond of the nth shortest length once, oriented from *i* to *j*
349 so that its angle lies in :math:`[0, 180)`.
351 The sign of *hop['ang']* is what marks a hopping as belonging to the
352 upper or the lower part of the Hamiltonian (see *get_ham*). Sites sorted
353 by (y, x), as *get_lattice* leaves them, give that orientation for
354 :math:`i < j`, but *rotation*, *clean_coor*, *+=*, ... reorder the
355 sites, so the orientation is taken from the angle itself.
357 Returns
358 -------
359 i, j : Integer ndarrays
360 Bond end points.
361 ang : Real ndarray
362 Bond angles, in degrees, in :math:`[0, 180)`.
363 '''
364 error_handling.torus_size(self.dist_uni[n], self._shortest_period(), n)
365 if self._pairs is not None:
366 pairs = self._pairs[np.isclose(self.dist_uni[n], self._pair_dis, atol=ATOL)]
367 pairs = pairs[np.lexsort((pairs[:, 1], pairs[:, 0]))]
368 i, j = pairs[:, 0], pairs[:, 1]
369 up = _upper(self._angles(i, j))
370 i, j = np.where(up, i, j), np.where(up, j, i)
371 order = np.lexsort((j, i))
372 i, j = i[order], j[order]
373 return i, j, np.maximum(self._angles(i, j), 0.)
374 ind = np.argwhere(np.isclose(self.dist_uni[n], self.vec_hop['dis'], atol=ATOL))
375 ind = ind[_upper(self.vec_hop['ang'][ind[:, 0], ind[:, 1]])]
376 ang = np.maximum(self.vec_hop['ang'][ind[:, 0], ind[:, 1]], 0.)
377 return ind[:, 0], ind[:, 1], ang
379 def fill_store_hop(self, n: int) -> None:
380 '''
381 Private method.
383 Store in *store_hop* indices, positive angles, and tags
384 of a given type of hopping (see *get_bonds*).
385 '''
386 i, j, ang = self.get_bonds(n)
387 hop = np.zeros(len(i), dtype=HOP_DTYPE[:-1])
388 hop['i'] = i
389 hop['j'] = j
390 hop['ang'] = ang
391 hop['tag'] = npc.add(self.lat.coor['tag'][i], self.lat.coor['tag'][j])
392 self.store_hop[n] = hop
394 def set_hopping(self, list_hop: list[dict], upper_part: bool = True) -> None:
395 r'''
396 Set lattice hoppings.
398 Parameters
399 ----------
400 list_hop : List of Dictionaries
401 Dictionary with keys ('n', 'ang', 'tag', 't') where:
403 * 'n' Positive integer, type of hoppings:
405 * 'n': 1 for nearest neighbours.
406 * 'n': 2 for next-nearest neighbours.
407 * 'n': 3 for next-next-nearest neighbours.
408 * etc...
410 * 'ang' value, float, angle, in deg, of the hoppings. (optional).
412 Hopping angles are given by the method *print_distances*.
414 * If :math:`ang \in[0, 180)`, fill the Hamiltonian upper part.
415 * If :math:`ang \in[-180, 0)`, fill the Hamiltonian lower part.
417 * 'tag' string of length 2 (optional).
419 Hopping tags.
421 * 't' Complex number, or value function.
423 Hopping value. A value function ``t(site_i, site_j, **params)``
424 is evaluated by ``get_ham(**params)``: *site_i* and *site_j*
425 are structured arrays of the bonds' end points (fields 'x',
426 'y', ('z'), 'tag', 'index'), and it returns one value, or
427 one per bond (see the example below). Modifiers applied
428 afterwards (*set_peierls_phase*, *set_hopping_dis*, ...)
429 multiply it.
430 upper_part : Boolean, default True
431 * True, get the hoppings with angles in :math:`[0, 180)`, *i.e.* fill
432 the Hamiltonian upper part.
433 * False, get the reversed hoppings, with angles in :math:`[-180, 0)`,
434 *i.e.* fill the Hamiltonian lower part.
436 Examples
437 --------
438 Example usage::
440 # fill upper part:
441 sys.set_hopping([{'n': 1, 't': 1.}])
442 # fill lower part:
443 sys.set_hopping([{'n': 1, 't': 1.}], upper_part=False)
444 # fill upper part: specifying the angles:
445 sys.set_hopping([{'n': 1, 'ang': 0., 't': 1.}, {'n': 1, 'ang': 90., 't': 2.}])
446 # fill lower part:
447 sys.set_hopping([{'n': 1, 'ang': -180., 't': 1.}, {'n': 1, 'ang': -90., 't': 2.}],
448 upper_part=False)
449 # fill upper part: specifying the tags:
450 sys.set_hopping([{'n': 1, 'tag': 'ab', 't': 1.}, {'n': 1, 'tag': 'ba', 't': 2.}])
451 # fill lower part:
452 sys.set_hopping([{'n': 1, 'tag': 'ab', 't': 1.}, {'n': 1, 'tag': 'ba', 't': 2.}],
453 upper_part=False)
454 # fill upper part: specifying the angles and tags:
455 sys.set_hopping([{'n': 1, 'ang': 0., 'tag': 'ab', 't': 1.},
456 {'n': 1, 'ang': 0., 'tag': 'ba', 't': 2.},
457 {'n': 1, 'ang': 90., 'tag': 'ab', 't': 3.},
458 {'n': 1, 'ang': 90., 'tag': 'ba', 't': 4.}])
459 # fill lower part:
460 sys.set_hopping([{'n': 1, 'ang': -180., 'tag': 'ab', 't': 1.},
461 {'n': 1, 'ang': -180., 'tag': 'ba', 't': 2.},
462 {'n': 1, 'ang': -90., 'tag': 'ab', 't': 3.},
463 {'n': 1, 'ang': -90., 'tag': 'ba', 't': 4.}],
464 upper_part=False)
465 # value function: a hopping that decays along x, swept without
466 # rebuilding the model
467 sys.set_hopping([{'n': 1, 't': lambda si, sj, t0, xi: t0 * np.exp(-si['x'] / xi)}])
468 for xi in (5., 10., 20.):
469 sys.get_ham(t0=1., xi=xi)
471 .. note::
473 A Hermitian hopping matrix can be built using only
474 its upper part OR only its lower part. The full matrix is then
475 built automatically by Hermitian conjugation.
477 If both the upper AND the lower parts are used to build the hopping
478 matrix, no Hermitian conjugation is performed, *i.e.* a non-Hermitian
479 hopping matrix can be built.
480 '''
481 error_handling.sites(self.lat.sites)
482 error_handling.boolean(upper_part, 'upper_part')
483 self.get_distances(error_handling.max_order(list_hop))
484 self.nmax = len(self.dist_uni) - 1
485 error_handling.set_hopping(list_hop, self.nmax)
486 list_n = np.unique([dic['n'] for dic in list_hop])
487 # fill, if needed self.store_hop
488 self.check_sites()
489 for n in list_n:
490 if n not in self.store_hop:
491 self.fill_store_hop(n)
492 # fill self.hop
493 for dic in list_hop:
494 if len(dic) == 2:
495 size = len(self.store_hop[dic['n']])
496 if upper_part:
497 mask = (self.hop['n'] == dic['n']) & (self.hop['ang'] >= 0)
498 else:
499 mask = (self.hop['n'] == dic['n']) & (self.hop['ang'] < 0)
500 if np.sum(mask):
501 self.hop = self.hop[np.logical_not(mask)]
502 ind = np.ones(size, bool)
503 hop = self.set_given_hopping(dic['n'], size, dic, ind, upper_part=upper_part)
504 elif len(dic) == 3 and 'ang' in dic:
505 error_handling.angle(dic['ang'], np.unique(self.store_hop[dic['n']]['ang']), upper_part)
506 if dic['ang'] >= 0:
507 ang_store = dic['ang']
508 else:
509 ang_store = dic['ang'] + 180.
510 size = np.sum(np.isclose(ang_store, self.store_hop[dic['n']]['ang'], atol=ATOL))
511 mask = (self.hop['n'] == dic['n']) & np.isclose(self.hop['ang'], dic['ang'], atol=ATOL)
512 if np.sum(mask):
513 self.hop = self.hop[np.logical_not(mask)]
514 ind = np.isclose(ang_store, self.store_hop[dic['n']]['ang'], atol=ATOL)
515 error_handling.index(ind, dic)
516 hop = self.set_given_hopping(dic['n'], size, dic, ind, upper_part=upper_part)
517 elif len(dic) == 3 and 'tag' in dic:
518 if upper_part:
519 tag_store = dic['tag']
520 else:
521 tag_store = dic['tag'][::-1]
522 size = np.sum(self.store_hop[dic['n']]['tag'] == tag_store)
523 mask = ((self.hop['n'] == dic['n']) & (self.hop['tag'] == dic['tag']))
524 if upper_part:
525 mask &= self.hop['ang'] >= 0
526 else:
527 mask &= self.hop['ang'] < 0
528 if np.sum(mask):
529 self.hop = self.hop[np.logical_not(mask)]
530 ind = self.store_hop[dic['n']]['tag'] == tag_store
531 error_handling.index(ind, dic)
532 hop = self.set_given_hopping(dic['n'], size, dic, ind, upper_part=upper_part)
533 else:
534 error_handling.angle(dic['ang'], np.unique(self.store_hop[dic['n']]['ang']), upper_part=upper_part)
535 error_handling.tag(dic['tag'], np.unique(self.store_hop[dic['n']]['tag']))
536 if dic['ang'] >= 0:
537 ang_store = dic['ang']
538 else:
539 ang_store = dic['ang'] + 180.
540 if upper_part:
541 tag_store = dic['tag']
542 else:
543 tag_store = dic['tag'][::-1]
544 size = np.sum((self.store_hop[dic['n']]['tag'] == tag_store) &
545 (np.isclose(ang_store, self.store_hop[dic['n']]['ang'], atol=ATOL)))
546 bool1 = (self.hop['n'] == dic['n']) & (self.hop['tag'] == dic['tag'])
547 bool2 = np.isclose(self.hop['ang'], dic['ang'], atol=ATOL)
548 mask = bool1 & bool2
549 if np.sum(mask):
550 self.hop = self.hop[np.logical_not(mask)]
551 ind = ((self.store_hop[dic['n']]['tag'] == tag_store) &
552 (np.isclose(ang_store, self.store_hop[dic['n']]['ang'], atol=ATOL)))
553 error_handling.index(ind, dic)
554 hop = self.set_given_hopping(dic['n'], size, dic, ind, upper_part=upper_part)
555 self.hop = np.concatenate([self.hop, hop])
556 self._assign_hopping(np.arange(len(self.hop) - len(hop), len(self.hop)), dic['t'])
558 def check_sites(self) -> None:
559 '''
560 Private method.
561 Check if the number of sites was changed after calling the
562 method system.set_hopping().
563 '''
564 if self.sites != self.lat.sites:
565 self.store_hop = {}
566 self.sites = self.lat.sites
568 def set_given_hopping(
569 self, n: int, size: int, dic: dict, mask: NDArray, upper_part: bool,
570 ) -> NDArray:
571 '''
572 Private method.
573 Fill self.hop.
575 Parameters
576 ----------
577 n : Integer
578 Hopping type.
579 size : Integer
580 Number of hoppings.
581 dic : Dictionary
582 Hopping dictionary.
583 mask : np.ndarray
584 Mask.
585 upper_part : Boolean
586 If True, keep the orientation of *store_hop*
587 (angles in [0, 180)), else reverse it.
588 '''
589 hop = np.empty(size, dtype=HOP_DTYPE)
590 hop['n'] = dic['n']
591 hop['t'] = 1. if callable(dic['t']) else dic['t']
592 if upper_part:
593 hop['i'] = self.store_hop[n]['i'][mask]
594 hop['j'] = self.store_hop[n]['j'][mask]
595 hop['ang'] = self.store_hop[n]['ang'][mask]
596 hop['tag'] = self.store_hop[n]['tag'][mask]
597 else:
598 hop['i'] = self.store_hop[n]['j'][mask]
599 hop['j'] = self.store_hop[n]['i'][mask]
600 hop['ang'] = self.store_hop[n]['ang'][mask] - 180
601 hop['tag'] = npc.add(self.lat.coor['tag'][hop['i']],
602 self.lat.coor['tag'][hop['j']])
603 return hop
605 def set_hopping_manual(self, dict_hop: dict[tuple[int, int], complex], upper_part: bool = True) -> None:
606 '''
607 Set hoppings manually.
609 Parameters
610 ----------
611 dict_hop : Dictionary of hoppings
612 key: hopping indices, val: hopping values.
613 upper_part : Boolean, default True
614 * True, fill the Hamiltonian upper part.
615 * False, fill the Hamiltonian lower part.
617 Examples
618 --------
619 Example usage::
621 sys.set_hopping_manual({(0, 1): 1., (1, 2): 2.})
622 '''
623 error_handling.sites(self.lat.sites)
624 error_handling.set_hopping_def(None, dict_hop, self.lat.sites, 'dict_hop', same_site=True)
625 hop = np.zeros(len(dict_hop), dtype=HOP_DTYPE)
626 i = [h[0] for h in dict_hop.keys()]
627 j = [h[1] for h in dict_hop.keys()]
628 t = [val for val in dict_hop.values()]
629 hop['i'], hop['j']= i, j
630 hop['t'] = t
631 hop['tag'] = npc.add(self.lat.coor['tag'][i],
632 self.lat.coor['tag'][j])
633 d = self._bond_vectors(np.array(i, int), np.array(j, int))
634 ang = 180 / PI * np.arctan2(d[:, 1], d[:, 0])
635 if upper_part:
636 ang[ang < 0] += 180
637 else:
638 ang[ang >= 0] -= 180
639 hop['ang'] = ang
640 self.hop = np.concatenate([self.hop, hop])
641 self._assign_hopping(np.arange(len(self.hop) - len(hop), len(self.hop)), hop['t'])
643 def _assign_hopping(self, rows: NDArray, t) -> None:
644 '''
645 Private method. Set the values of the rows *rows* (indices or mask)
646 of *hop* to *t*: number(s), or a value function (the rows then hold
647 the factor 1). It replaces any value function set on the same bonds.
648 '''
649 keys = self._bond_keys(rows)
650 self._hop_values = [(k[~np.isin(k, keys)], f) for k, f in self._hop_values]
651 self._hop_values = [(k, f) for k, f in self._hop_values if k.size]
652 if callable(t):
653 self.hop['t'][rows] = 1.
654 self._hop_values.append((keys, t))
655 else:
656 self.hop['t'][rows] = t
658 def _bond_keys(self, rows: NDArray | slice) -> NDArray[np.int64]:
659 '''
660 Private method. One integer per directed bond (i, j) of the rows *rows* of *hop*.
661 '''
662 return (self.hop['i'][rows].astype('i8') << 32) | self.hop['j'][rows].astype('i8')
664 def _sites(self) -> NDArray:
665 '''
666 Private method. All the sites, as passed to value functions (see *tbkit.values.sites*).
667 '''
668 return values.sites(self._coords(), self.lat.coor['tag'], np.arange(self.lat.sites))
670 def _hop_t(self, params: dict) -> NDArray[np.complex128]:
671 '''
672 Private method. The hopping values of *hop*, value functions evaluated with *params*.
673 '''
674 t = self.hop['t'].copy()
675 if not self._hop_values:
676 return t
677 keys, sites = self._bond_keys(slice(None)), self._sites()
678 for k, func in self._hop_values:
679 rows = np.flatnonzero(np.isin(keys, k))
680 t[rows] *= values.evaluate(func, (sites[self.hop['i'][rows]], sites[self.hop['j'][rows]]),
681 params, len(rows))
682 return t
684 def _onsite_e(self, params: dict) -> NDArray[np.complex128]:
685 '''
686 Private method. The onsite energies, value functions evaluated with *params*.
687 '''
688 onsite = self.onsite.copy()
689 if self._onsite_values:
690 sites = self._sites()
691 for idx, func in self._onsite_values:
692 onsite[idx] += values.evaluate(func, (sites[idx],), params, len(idx))
693 return onsite
695 def set_hopping_dis(self, alpha: complex) -> None:
696 '''
697 Set uniform hopping disorder.
699 Parameters
700 ----------
701 alpha : Complex or Real number
702 Disorder strength.
704 Examples
705 --------
706 Example usage::
708 sys.set_hopping_dis(alpha=0.1)
709 '''
710 error_handling.empty_hop(self.hop)
711 error_handling.number(alpha, 'alpha')
712 self.hop['t'] *= 1. + alpha * rand.uniform(-1., 1., len(self.hop))
714 def set_peierls_phase(self, phase: Callable[..., NDArray]) -> None:
715 r'''
716 Apply the Peierls substitution to the existing hoppings, to
717 capture the effect of an orbital magnetic field:
719 .. math::
721 t_{ij} \to t_{ij}\, e^{i\phi_{ij}}\, ,\quad
722 \phi_{ij} = \frac{2\pi}{\Phi_0}\int_{\mathbf{r}_i}^{\mathbf{r}_j}
723 \mathbf{A}\cdot d\mathbf{l}
725 where :math:`\mathbf{A}` is the vector potential, integrated along
726 the straight bond from site :math:`i` to site :math:`j`.
728 Must be called after *set_hopping* / *set_hopping_manual* (it
729 rescales the existing hoppings *in place*) and before *get_ham*.
730 For a uniform perpendicular field, use the convenience method
731 *set_magnetic_field* instead.
733 Parameters
734 ----------
735 phase : Callable
736 ``phase(xi, yi, xj, yj)`` returns
737 :math:`\phi_{ij}`, the (real-valued) Peierls phase for the bond
738 from :math:`(x_i, y_i)` to :math:`(x_j, y_j)`. Called with
739 Numpy arrays (one value per hopping in *sys.hop*). With periodic
740 boundaries, :math:`(x_j, y_j)` is the end of the short bond, i.e.
741 the image of site :math:`j` next to site :math:`i`, which lies
742 outside the sample for a bond that wraps around. The vector
743 potential must then be consistent with the torus: periodic (a
744 field of zero net flux), or uniform (a flux threading the holes
745 of the torus, i.e. twisted boundary conditions).
747 Notes
748 -----
749 .. note::
751 *get_ham* automatically assigns the reversed bond its complex
752 conjugate, so the Hamiltonian stays Hermitian as long as
753 *phase* is antisymmetric under swapping :math:`i` and
754 :math:`j` -- true for the line integral of any vector
755 potential, since reversing the integration path negates it.
757 Examples
758 --------
759 Example usage::
761 # Peierls phase from a uniform field via the Landau gauge
762 # A = (0, B x): the same field as set_magnetic_field(alpha),
763 # just in a different (equally valid) gauge.
764 alpha = 0.01
765 sys.set_peierls_phase(lambda xi, yi, xj, yj:
766 2 * np.pi * alpha * (yj - yi) * (xi + xj) / 2)
767 '''
768 error_handling.empty_hop(self.hop)
769 error_handling.is_callable(phase, 'phase')
770 xi = self.lat.coor['x'][self.hop['i']]
771 yi = self.lat.coor['y'][self.hop['i']]
772 d = self._bond_vectors(self.hop['i'].astype(int), self.hop['j'].astype(int))
773 xj, yj = xi + d[:, 0], yi + d[:, 1]
774 self.hop['t'] = self.hop['t'] * np.exp(1j * phase(xi, yi, xj, yj))
776 def set_magnetic_field(self, alpha: float) -> None:
777 r'''
778 Set a uniform perpendicular magnetic field via the Peierls
779 substitution (see *set_peierls_phase*), using the symmetric gauge
780 :math:`\mathbf{A} = \frac{B}{2}(-y, x)`:
782 .. math::
784 \phi_{ij} = \pi\alpha\,(x_iy_j - x_jy_i)
786 Parameters
787 ----------
788 alpha : Real number
789 Flux density :math:`B/\Phi_0`, in flux
790 quanta per unit area (in the lattice's length units) -- *i.e.*
791 the flux through a region of area :math:`S` is
792 :math:`\alpha S` flux quanta.
794 Examples
795 --------
796 Example usage::
798 # one flux quantum per 100 unit cells of a lattice with
799 # lattice constant 1:
800 sys.set_magnetic_field(alpha=0.01)
802 Open boundaries only: on a torus the symmetric gauge is not
803 periodic, and a uniform field needs a quantized total flux and
804 magnetic boundary conditions.
805 '''
806 error_handling.not_periodic(self.periodic, 'set_magnetic_field')
807 error_handling.real_number(alpha, 'alpha')
808 self.set_peierls_phase(lambda xi, yi, xj, yj: PI * alpha * (xi*yj - xj*yi))
810 def set_onsite_dis(self, alpha: complex) -> None:
811 '''
812 Set uniform onsite disorder.
814 Parameters
815 ----------
816 alpha : Complex or Real number
817 Disorder strength.
819 Examples
820 --------
821 Example usage::
823 sys.set_onsite_dis(alpha=0.1)
824 '''
825 error_handling.empty_onsite(self.onsite)
826 error_handling.number(alpha, 'alpha')
827 self.onsite += alpha * rand.uniform(-1., 1., self.lat.sites)
829 def set_onsite_def(self, onsite_def: dict[int, complex]) -> None:
830 '''
831 Set specific onsite energies.
833 Parameters
834 ----------
835 onsite_def : Dictionary
836 key: site indices, val: onsite values.
838 Examples
839 --------
840 Example usage::
842 sys.set_onsite_def({0: 1., 1: -1j})
843 '''
844 error_handling.empty_onsite(self.onsite)
845 error_handling.set_onsite_def(onsite_def, self.lat.sites)
846 for i, o in onsite_def.items():
847 self.onsite[i] = o
848 self._onsite_values = [(idx[~np.isin(idx, list(onsite_def))], f) for idx, f in self._onsite_values]
850 def set_hopping_def(self, hopping_def: dict[tuple[int, int], complex]) -> None:
851 '''
852 Set specific hoppings.
854 Parameters
855 ----------
856 hopping_def : Dictionary of hoppings
857 key: hopping indices :math:`(i, j)`, val: hopping value :math:`t_{ij}`.
858 If *sys.hop* stores the bond the other way round, as :math:`(j, i)`
859 (only one direction is stored for a Hermitian Hamiltonian), that
860 entry is set to :math:`t_{ij}^*`, so that :math:`H_{ij} = t_{ij}`.
861 A pair that is not a hopping of *sys.hop* is ignored.
863 Examples
864 --------
865 Example usage::
867 sys.set_hopping_def({(0, 1): 1., (1, 2): -1j})
868 '''
869 error_handling.empty_hop(self.hop)
870 error_handling.set_hopping_def(self.hop, hopping_def, self.lat.sites)
871 for key, val in hopping_def.items():
872 cond = (self.hop['i'] == key[0]) & (self.hop['j'] == key[1])
873 if cond.any():
874 self._assign_hopping(cond, val)
875 continue
876 cond = (self.hop['i'] == key[1]) & (self.hop['j'] == key[0])
877 self._assign_hopping(cond, np.conj(val))
879 def set_new_hopping(self, list_hop: list[dict], ind: NDArray) -> None:
880 '''
881 Private method.
882 Set new hoppings.
884 Parameters
885 ----------
886 list_hop : List of Dictionary (see set_hopping definition)
887 ind : List
888 List of indices.
889 '''
890 for dic in list_hop:
891 if len(dic) == 2:
892 self._assign_hopping(ind, dic['t'])
893 elif len(dic) == 3 and 'ang' in dic:
894 self._assign_hopping(ind & np.isclose(self.hop['ang'], dic['ang'], atol=ATOL), dic['t'])
895 elif len(dic) == 3 and 'tag' in dic:
896 self._assign_hopping(ind & (self.hop['tag'] == dic['tag']), dic['t'])
897 else:
898 self._assign_hopping(ind & (self.hop['tag'] == dic['tag'])
899 & np.isclose(self.hop['ang'], dic['ang'], atol=ATOL), dic['t'])
901 def find_square(self, xlims: tuple[float, float], ylims: tuple[float, float]) -> NDArray:
902 '''
903 Private method.
904 Find hoppings within the square.
906 Parameters
907 ----------
908 xlims : List or Tuple
909 :math:`x` interval.
910 ylims : List or Tuple
911 :math:`y` interval.
912 '''
913 error_handling.lims(xlims)
914 error_handling.lims(ylims)
915 in1 = (self.lat.coor['x'][self.hop['i']] >= xlims[0]) & \
916 (self.lat.coor['y'][self.hop['i']] >= ylims[0]) & \
917 (self.lat.coor['x'][self.hop['j']] >= xlims[0]) & \
918 (self.lat.coor['y'][self.hop['j']] >= ylims[0])
919 in2 = (self.lat.coor['x'][self.hop['i']] <= xlims[1]) & \
920 (self.lat.coor['y'][self.hop['i']] <= ylims[1]) & \
921 (self.lat.coor['x'][self.hop['j']] <= xlims[1]) & \
922 (self.lat.coor['y'][self.hop['j']] <= ylims[1])
923 return in1 * in2
925 def find_ellipse(self, rx: float, ry: float, x0: float, y0: float) -> NDArray:
926 '''
927 Private method.
928 Find hoppings within the ellipse.
930 Parameters
931 ----------
932 rx : Positive Float
933 Radius along :math:`x`.
934 ry : Positive Float
935 Radius along :math:`y`.
936 x0 : Float, default 0
937 :math:`x` center.
938 y0 : Float, default 0
939 :math:`y` center.
940 '''
941 in1 = (self.lat.coor['x'][self.hop['i']] - x0) ** 2 / rx ** 2 + \
942 (self.lat.coor['y'][self.hop['i']] - y0) ** 2 / ry ** 2 <= 1.
943 in2 = (self.lat.coor['x'][self.hop['j']] - x0) ** 2 / rx ** 2 + \
944 (self.lat.coor['y'][self.hop['j']] - y0) ** 2 / ry ** 2 <= 1.
945 return in1 * in2
947 def change_hopping_square(
948 self, list_hop: list[dict], xlims: tuple[float, float], ylims: tuple[float, float] = (-1., 1.),
949 ) -> None:
950 '''
951 Change the values of the hoppings whose two sites both lie in
952 the rectangle *xlims* x *ylims*.
954 Parameters
955 ----------
956 list_hop : List of Dictionary (see set_hopping definition)
957 xlims : List or Tuple
958 :math:`x` interval.
959 ylims : List or Tuple, default [-1, 1]
960 :math:`y` interval.
961 '''
962 error_handling.empty_hop(self.hop)
963 error_handling.set_hopping(list_hop, self.nmax)
964 ind = self.find_square(xlims, ylims)
965 self.set_new_hopping(list_hop, ind)
967 def change_hopping_ellipse(
968 self, list_hop: list[dict], rx: float, ry: float, x0: float = 0., y0: float = 0.,
969 ) -> None:
970 '''
971 Change the values of the hoppings whose two sites both lie in
972 the ellipse of radii *rx*, *ry* centered at (*x0*, *y0*).
974 Parameters
975 ----------
976 list_hop : List of Dictionary (see set_hopping definition)
977 rx : Positive Float
978 Radius along :math:`x`.
979 ry : Positive Float
980 Radius along :math:`y`.
981 x0 : Float, default 0
982 :math:`x` center.
983 y0 : Float, default 0
984 :math:`y` center.
985 '''
986 error_handling.empty_hop(self.hop)
987 error_handling.set_hopping(list_hop, self.nmax)
988 error_handling.positive_real(rx, 'rx')
989 error_handling.positive_real(ry, 'ry')
990 error_handling.real_number(x0, 'x0')
991 error_handling.real_number(y0, 'y0')
992 ind = self.find_ellipse(rx, ry, x0, y0)
993 self.set_new_hopping(list_hop, ind)
995 def get_coor_hop(self) -> None:
996 '''
997 Get the site coordinates in hopping space
998 only considering the nearest neighbours hoppings.
999 '''
1000 error_handling.empty_hop(self.hop)
1001 visited = np.zeros(self.lat.sites, 'u2')
1002 self.coor_hop = np.zeros(self.lat.sites, dtype=COOR_DTYPE)
1003 self.coor_hop['tag'] = self.lat.coor['tag']
1004 hop = self.hop[self.hop['n'] == 1]
1005 hop_down = np.copy(hop)
1006 hop_down['i'] = hop['j']
1007 hop_down['j'] = hop[ 'i']
1008 hop_down['ang'] = -180 + hop['ang']
1009 hop = np.concatenate([hop, hop_down])
1010 i_visit = np.min(hop['i'])
1011 while True:
1012 hs = hop[hop['i'] == i_visit]
1013 for h in hs:
1014 if visited[h['j']] == 2:
1015 continue
1016 self.coor_hop['x'][h['j']] = self.coor_hop['x'][i_visit] + \
1017 h['t'].real*cos(PI / 180 * h['ang'])
1018 self.coor_hop['y'][h['j']] = self.coor_hop['y'][i_visit] + \
1019 h['t'].real*sin(PI / 180 * h['ang'])
1020 visited[h['j']] = 1
1021 visited[i_visit] = 2
1022 explored = np.argwhere(visited == 1)
1023 if not explored.any():
1024 break
1025 i_visit = explored[0, 0]
1027 def get_ham(self, **params) -> None:
1028 '''
1029 Get the Tight-Binding Hamiltonian using sys.hop.
1031 Parameters
1032 ----------
1033 params
1034 Values of the parameters of the value functions given
1035 to *set_hopping* / *set_onsite*: each function gets the ones its
1036 signature names (all of them with ``**kwargs``); the others are
1037 ignored. Only the values are recomputed: sweeping a parameter
1038 costs one call per point, not a new model.
1040 Examples
1041 --------
1042 Example usage::
1044 for B in np.linspace(0., 0.1, 50):
1045 sys.get_ham(B=B)
1046 sys.get_eig()
1047 '''
1048 error_handling.empty_hop(self.hop)
1049 error_handling.hop_sites(self.hop, self.lat.sites)
1050 t = self._hop_t(params)
1051 if np.all(self.hop['ang'] >= 0) or np.all(self.hop['ang'] < 0):
1052 self.ham = sparse.csr_matrix((t, (self.hop['i'], self.hop['j'])),
1053 shape=(self.lat.sites, self.lat.sites)) \
1054 + sparse.csr_matrix((t.conj(), (self.hop['j'], self.hop['i'])),
1055 shape=(self.lat.sites, self.lat.sites))
1056 else:
1057 self.ham = sparse.csr_matrix((t, (self.hop['i'], self.hop['j'])),
1058 shape=(self.lat.sites, self.lat.sites))
1059 if self.onsite.size == self.lat.sites:
1060 self.ham += sparse.diags(self._onsite_e(params), 0)
1062 def get_eig(self, eigenvec: bool = False, left: bool = False) -> None:
1063 '''
1064 Get the eigenergies, eigenvectors and polarisation.
1066 Parameters
1067 ----------
1068 eigenvec : Boolean, default False
1069 If True, get the eigenvectors.
1070 left : Boolean, default False
1071 If True, get the left eigenvectors too.
1072 Relevant for non-Hermitian matrices.
1073 '''
1074 error_handling.empty_ham(self.ham)
1075 error_handling.boolean(eigenvec, 'eigenvec')
1076 error_handling.boolean(left, 'left')
1077 # Drop any left eigenvectors from a previous call: they belong to the
1078 # previous diagonalization and must not be re-sorted with this one's
1079 # ordering (get_petermann would then pair the wrong vectors).
1080 self.ln = np.array([], 'c16')
1081 if eigenvec:
1082 if (self.ham.conj().T != self.ham).nnz:
1083 if not left:
1084 self.en, self.rn = LA.eig(self.ham.toarray())
1085 else:
1086 # scipy.linalg.eig returns (w, vl, vr) -- the LEFT
1087 # eigenvectors come back before the right ones.
1088 self.en, self.ln, self.rn = LA.eig(self.ham.toarray(), left=left)
1089 ind = np.argsort(self.en.real)
1090 self.en = self.en[ind]
1091 self.rn = self.rn[:, ind]
1092 if self.ln.size:
1093 self.ln = self.ln[:, ind]
1094 else:
1095 self.en, self.rn = LA.eigh(self.ham.toarray())
1096 self._set_states()
1097 else:
1098 if (self.ham.conj().T != self.ham).nnz:
1099 self.en = LA.eigvals(self.ham.toarray())
1100 ind = np.argsort(self.en.real)
1101 self.en = self.en[ind]
1102 else:
1103 ham = self.ham.toarray()
1104 # A real symmetric matrix (no fluxes, no spin-orbit) has the same
1105 # eigenvalues from the real solver, about twice as fast.
1106 if not ham.imag.any():
1107 ham = ham.real
1108 self.en = LA.eigvalsh(ham)
1110 def _set_states(self) -> None:
1111 '''
1112 Private method. Intensities and sublattice polarizations of the
1113 eigenvectors *rn*.
1114 '''
1115 self.intensity = np.abs(self.rn) ** 2
1116 self.pola = np.zeros((self.rn.shape[1], len(self.lat.tags)))
1117 for i, tag in enumerate(self.lat.tags):
1118 self.pola[:, i] = np.sum(self.intensity[self.lat.coor['tag'] == tag, :], axis=0)
1120 def get_eig_sparse(self, n_eig: int = 10, sigma: complex = 0., eigenvec: bool = False) -> None:
1121 r'''
1122 Get the *n_eig* eigenenergies closest to *sigma*, and optionally
1123 their eigenvectors, with the sparse shift-invert Lanczos (Hermitian)
1124 or Arnoldi (non-Hermitian) solver of ARPACK. Unlike *get_eig*, it
1125 never builds the dense matrix, so it reaches lattices of
1126 :math:`10^5` sites and more -- e.g. the states near a band centre, or
1127 the zero modes of a large flake. The results replace *en* (sorted by
1128 real part) and, with *eigenvec*, *rn*, *intensity* and *pola*, as for
1129 *get_eig*.
1131 Parameters
1132 ----------
1133 n_eig : Positive integer, default 10
1134 Number of eigenpairs
1135 (at most sites - 2).
1136 sigma : Real or complex number, default 0
1137 Target energy.
1138 The factorization uses *sigma* shifted by about 1e-9 (relative to
1139 the largest matrix element), so that it may be exactly an
1140 eigenvalue, e.g. zero modes at the default *sigma* = 0; the
1141 eigenvalues themselves are not shifted.
1142 eigenvec : Boolean, default False
1143 '''
1144 error_handling.empty_ham(self.ham)
1145 error_handling.n_eig(n_eig, self.lat.sites)
1146 error_handling.number(sigma, 'sigma')
1147 error_handling.boolean(eigenvec, 'eigenvec')
1148 self.ln = np.array([], 'c16')
1149 ham = self.ham.tocsc()
1150 hermitian = not (self.ham.conj().T != self.ham).nnz
1151 # Never factorize H - sigma at sigma itself: when sigma is exactly an
1152 # eigenvalue, SuperLU usually raises, but with a two-fold one it
1153 # passes a negative size to BLAS, and Apple's Accelerate then aborts
1154 # the whole process (no exception to catch).
1155 target = sigma + 1e-9 * (1. + abs(sigma) + abs(ham).max())
1156 if hermitian:
1157 out = SLA.eigsh(ham, k=n_eig, sigma=float(np.real(target)), return_eigenvectors=eigenvec)
1158 else:
1159 out = SLA.eigs(ham, k=n_eig, sigma=target, return_eigenvectors=eigenvec)
1160 en, vec = out if eigenvec else (out, None)
1161 ind = np.argsort(en.real, kind='stable')
1162 self.en = en[ind]
1163 if vec is not None:
1164 self.rn = vec[:, ind]
1165 self._set_states()
1167 def get_green(self, energy: complex, eta: float = 1e-6) -> NDArray[np.complex128]:
1168 r'''
1169 Get the retarded Green's function (resolvent)
1171 .. math::
1173 G(E) = \left[(E + i\eta) - H\right]^{-1}\, .
1175 Parameters
1176 ----------
1177 energy : Real (or complex) number
1178 Energy :math:`E`.
1179 eta : Positive real or zero, default 1e-6
1180 Broadening :math:`\eta`.
1182 Returns
1183 -------
1184 G : Complex ndarray, shape (sites, sites)
1185 '''
1186 error_handling.empty_ham(self.ham)
1187 error_handling.number(energy, 'energy')
1188 error_handling.positive_real_zero(eta, 'eta')
1189 z = (energy + 1j * eta) * np.eye(self.ham.shape[0])
1190 return LA.inv(z - self.ham.toarray())
1192 def get_ldos(self, energies: ArrayLike, eta: float = 0.05) -> NDArray[np.float64]:
1193 r'''
1194 Get the local density of states
1196 .. math::
1198 \rho_i(E) = -\frac{1}{\pi}\,\mathrm{Im}\,G_{ii}(E)
1199 = \sum_n |\psi_n(i)|^2\,\frac{\eta/\pi}{(E - E_n)^2 + \eta^2}
1201 (the second form, from the eigenvectors, for a Hermitian
1202 Hamiltonian). Summed over sites it is the Lorentzian-broadened
1203 density of states.
1205 Parameters
1206 ----------
1207 energies : Real array
1208 Energies.
1209 eta : Positive real, default 0.05
1210 Broadening.
1212 Returns
1213 -------
1214 ldos : Real ndarray, shape (len(energies), sites)
1215 '''
1216 error_handling.empty_ham(self.ham)
1217 error_handling.positive_real(eta, 'eta')
1218 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1219 error_handling.ndarray_empty(energies, 'energies')
1220 ham = self.ham.toarray()
1221 if np.allclose(ham, ham.conj().T, atol=1e-14):
1222 en, vec = LA.eigh(ham)
1223 lor = (eta / PI) / ((energies[:, None] - en[None, :]) ** 2 + eta ** 2)
1224 return lor @ (np.abs(vec) ** 2).T
1225 eye = np.eye(len(ham))
1226 return np.array([-np.diag(LA.inv((e + 1j * eta) * eye - ham)).imag / PI
1227 for e in energies])
1229 def get_fermi_level(self, n_electrons: float, temperature: float = 0.) -> float:
1230 '''
1231 Get the Fermi level holding *n_electrons* electrons in the spectrum
1232 *en* (one electron per state; see *tbkit.occupation.fermi_level*).
1234 Parameters
1235 ----------
1236 n_electrons : Positive real or zero
1237 Number of electrons.
1238 temperature : Positive real or zero, default 0
1240 Returns
1241 -------
1242 e_fermi : Real number
1243 '''
1244 error_handling.empty_ndarray(self.en, 'sys.get_eig')
1245 return occupation.fermi_level(self.en, n_electrons, temperature)
1247 def get_occupations(self, e_fermi: float, temperature: float = 0.) -> NDArray[np.float64]:
1248 '''
1249 Get the Fermi-Dirac occupation of every eigenstate of *en*.
1251 Parameters
1252 ----------
1253 e_fermi : Real number
1254 Fermi level (see *get_fermi_level*).
1255 temperature : Positive real or zero, default 0
1257 Returns
1258 -------
1259 f : Real ndarray, same length as *en*
1260 '''
1261 error_handling.empty_ndarray(self.en, 'sys.get_eig')
1262 return occupation.fermi_dirac(self.en, e_fermi, temperature)
1264 def get_charge_density(self, e_fermi: float, temperature: float = 0.) -> NDArray[np.float64]:
1265 r'''
1266 Get the electron density on every site,
1267 :math:`n_i = \sum_n f(E_n)\,|\psi_n(i)|^2`.
1269 Parameters
1270 ----------
1271 e_fermi : Real number
1272 Fermi level.
1273 temperature : Positive real or zero, default 0
1275 Returns
1276 -------
1277 density : Real ndarray, shape (sites,)
1278 '''
1279 error_handling.empty_ndarray(self.rn, 'sys.get_eig(eigenvec=True)')
1280 return self.intensity @ self.get_occupations(e_fermi, temperature)
1282 def get_total_energy(self, e_fermi: float, temperature: float = 0.) -> float:
1283 r'''
1284 Get the band energy :math:`E = \sum_n f(E_n)\,E_n`.
1286 Parameters
1287 ----------
1288 e_fermi : Real number
1289 Fermi level.
1290 temperature : Positive real or zero, default 0
1292 Returns
1293 -------
1294 energy : Real number
1295 '''
1296 return float(np.sum(self.get_occupations(e_fermi, temperature) * self.en.real))
1298 def get_ipr(self) -> None:
1299 r'''
1300 Get the Inverse Participation Ratio:
1302 .. math::
1304 IPR_n = \sum_i|\psi_i^{n}|^4\, .
1305 '''
1306 error_handling.empty_ndarray(self.rn, 'sys.get_eig(eigenvec=True)')
1307 self.ipr = np.sum(self.intensity ** 2, axis=0)
1309 def get_petermann(self) -> None:
1310 r'''
1311 Get the Petermann factor:
1313 .. math::
1315 K_n = \frac{\langle\psi_L^{n}|\psi_L^{n}\rangle\langle\psi_R^{n}|\psi_R^{n}\rangle}{|\langle\psi_L^{n}|\psi_R^{n}\rangle|^2}\, .
1317 .. note::
1319 LA.eig fixes the norm such that :math:`\langle\psi_L^{n}|\psi_L^{n}\rangle = 1` and :math:`\langle\psi_R^{n}|\psi_R^{n}\rangle = 1`.
1320 '''
1321 if not (self.ham.conj().T != self.ham).nnz:
1322 self.petermann = np.ones(self.lat.sites)
1323 return
1324 error_handling.empty_ndarray(self.ln, 'sys.get_eig(eigenvec=True, left=True)')
1325 # LA.eig fixes each eigenvector's norm but not the relative phase
1326 # between the left and right ones, so the overlap carries an
1327 # arbitrary phase: take its modulus, not its real part.
1328 left_right = np.abs(np.sum(self.ln * np.conjugate(self.rn), axis=0))
1329 self.petermann = 1. / left_right ** 2
1331 def get_intensity_pola_max(self, tag_pola: str) -> NDArray[np.float64]:
1332 '''
1333 Get the state with largest polarization on one sublattice.
1335 Parameters
1336 ----------
1337 tag_pola : One-character string
1338 Sublattice tag.
1340 Returns
1341 -------
1342 intensity
1343 Intensity of max polarized state on *tag*.
1344 '''
1345 error_handling.empty_ndarray(self.rn, 'sys.get_eig(eigenvec=True)')
1346 error_handling.tag(tag_pola, self.lat.tags)
1347 i_tag = self.lat.tags == tag_pola
1348 ind = np.argmax(self.pola[:, i_tag])
1349 print('State with polarization: {:.5f}'.format(self.pola[ind, i_tag].item()))
1350 return self.intensity[:, ind]
1352 def get_intensity_pola_min(self, tag_pola: str) -> NDArray[np.float64]:
1353 '''
1354 Get the state with smallest polarization on one sublattice.
1356 Parameters
1357 ----------
1358 tag_pola : One-character string
1359 Sublattice tag.
1361 Returns
1362 -------
1363 intensity
1364 Intensity of min polarized state on *tag*.
1365 '''
1366 error_handling.empty_ndarray(self.rn, 'sys.get_eig(eigenvec=True)')
1367 error_handling.tag(tag_pola, self.lat.tags)
1368 i_tag = self.lat.tags == tag_pola
1369 ind = np.argmin(self.pola[:, i_tag])
1370 print('State with polarization: {:.5f}'.format(self.pola[ind, i_tag].item()))
1371 return self.intensity[:, ind]
1373 def get_intensity_en(self, lims: tuple[float, float]) -> NDArray[np.float64]:
1374 '''
1375 Get, if any, the intensity of the sum of the states whose
1376 energies (real part) lie between *lims[0]* and *lims[1]*.
1378 Parameters
1379 ----------
1380 lims : List
1381 lims[0] energy min, lims[1] energy max.
1383 Returns
1384 -------
1385 intensity
1386 Sum of the intensities between (lims[0], lims[1]).
1387 '''
1388 error_handling.empty_ndarray(self.rn, 'sys.get_eig(eigenvec=True)')
1389 error_handling.lims(lims)
1390 ind = np.ravel(np.where((self.en.real > lims[0]) & (self.en.real < lims[1])))
1391 print('{} states between {} and {}'.format(len(ind), lims[0], lims[1]))
1392 return np.sum(self.intensity[:, ind], axis=1)
1395 def get_local_chern_marker(self, e_fermi: float = 0., area: float | None = None) -> NDArray[np.float64]:
1396 r'''
1397 Get the local Chern marker of the states below *e_fermi* (Bianco and
1398 Resta, Phys. Rev. B 84, 241106 (2011)),
1400 .. math::
1402 C(\mathbf{r}_i) = -\frac{4\pi}{a}\,\mathrm{Im}\,
1403 \langle i|\,P\,x\,Q\,y\,P\,|i\rangle\, ,
1405 with :math:`P` the projector on the occupied states, :math:`Q = 1-P`,
1406 and :math:`a` the area per site. Averaged over the bulk of a finite
1407 flake -- even a disordered or amorphous one, where no Brillouin zone
1408 exists -- it gives the Chern number, in the sign convention of
1409 *KSpace.chern_number*; it sums to zero over the whole flake, the
1410 edges carrying the opposite value.
1412 Parameters
1413 ----------
1414 e_fermi : Real number, default 0
1415 Fermi energy: the states
1416 with energies below it are occupied.
1417 area : Positive real number
1418 Default value None: the area of
1419 the unit cell of *lat* divided by its number of sites.
1421 Returns
1422 -------
1423 marker : Real ndarray, shape (sites,)
1424 '''
1425 error_handling.not_periodic(self.periodic, 'get_local_chern_marker')
1426 error_handling.empty_ham(self.ham)
1427 error_handling.real_number(e_fermi, 'e_fermi')
1428 error_handling.hermitian(self.ham)
1429 if area is None:
1430 error_handling.planar_cell(self.lat.prim_vec)
1431 (a1x, a1y), (a2x, a2y) = self.lat.prim_vec
1432 area = abs(a1x * a2y - a1y * a2x) / len(self.lat.unit_cell)
1433 error_handling.positive_real(area, 'area')
1434 en, vec = LA.eigh(self.ham.toarray())
1435 occ = vec[:, en < e_fermi]
1436 proj = occ @ occ.conj().T
1437 rows = self._row_sites()
1438 x = self.lat.coor['x'][rows]
1439 y = self.lat.coor['y'][rows]
1440 # P x Q y P = P x y P - P x P y P, diagonal only
1441 pxq = proj * x[None, :] - (proj * x[None, :]) @ proj
1442 qyp = y[:, None] * proj - proj @ (y[:, None] * proj)
1443 per_row = -4 * PI / area * np.imag(np.sum(proj * (pxq @ qyp).T, axis=1))
1444 marker = np.zeros(self.lat.sites)
1445 np.add.at(marker, rows, per_row)
1446 return marker
1448 def get_bott_index(self, e_fermi: float = 0.) -> float:
1449 r'''
1450 Get the Bott index of the states below *e_fermi* of a 2D lattice
1451 with periodic boundaries along both primitive vectors (Loring and
1452 Hastings, EPL 92, 67004 (2010); see *tbkit.topology.bott_index*):
1453 the Chern number of a torus, in the sign convention of
1454 *KSpace.chern_number*, and an integer for any disordered sample
1455 whose occupied states are gapped or localized.
1457 Parameters
1458 ----------
1459 e_fermi : Real number, default 0
1460 Fermi energy.
1462 Returns
1463 -------
1464 bott : Real number, close to an integer
1466 Examples
1467 --------
1468 Example usage::
1470 lat = lattices.honeycomb()
1471 lat.get_lattice(12, 12)
1472 sys = System(lat, periodic=True)
1473 sys.set_hopping([{'n': 1, 't': 1.}])
1474 ...
1475 sys.get_ham()
1476 sys.get_bott_index()
1477 '''
1478 from tbkit.topology import bott_index
1479 error_handling.empty_ham(self.ham)
1480 error_handling.torus_2d(self.periodic)
1481 return bott_index(self.ham, self._coords()[self._row_sites()], self._box(), e_fermi)
1483 def _row_sites(self) -> NDArray:
1484 '''
1485 Private method. Site of every row of the Hamiltonian (one row per site here).
1486 '''
1487 return np.arange(self.lat.sites)
1490# Backward-compatible lowercase alias (pre-0.2 API).
1491system = System