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