Coverage for tbkit/kspace.py: 100%
218 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-22 13:16 +0100
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-22 13:16 +0100
1from __future__ import annotations
3from typing import Sequence
5import numpy as np
6from numpy.typing import ArrayLike, NDArray
7import matplotlib.pyplot as plt
8from matplotlib.figure import Figure
9import scipy.linalg as LA
10import tbkit.error_handling as error_handling
11import tbkit.dos as dos
12from tbkit.lattice import Lattice
15PI = np.pi
17#: Pauli matrices (plus the identity, key ``'0'``), for building spinful
18#: hoppings/onsite terms (spin-orbit coupling, Zeeman splitting, ...) when
19#: :class:`KSpace` is constructed with ``spin=True``.
20PAULI = {
21 '0': np.eye(2, dtype='c16'),
22 'x': np.array([[0., 1.], [1., 0.]], dtype='c16'),
23 'y': np.array([[0., -1j], [1j, 0.]], dtype='c16'),
24 'z': np.array([[1., 0.], [0., -1.]], dtype='c16'),
25}
28#################################
29# CLASS KSPACE
30#################################
33def reciprocal_vectors(prim_vec: list[tuple[float, float]]) -> list[tuple[float, float]]:
34 r'''
35 Get the reciprocal lattice vectors :math:`\mathbf{b}_i` such that
36 :math:`\mathbf{a}_i\cdot\mathbf{b}_j = 2\pi\delta_{ij}`.
38 :param prim_vec: List of one/two tuples. Primitive vectors (see class **lattice**).
40 :returns:
41 * **rec_vec** -- List of one/two tuples. Reciprocal vectors.
42 '''
43 if len(prim_vec) == 1:
44 ax, ay = prim_vec[0]
45 norm2 = ax ** 2 + ay ** 2
46 return [(2*PI*ax/norm2, 2*PI*ay/norm2)]
47 (a1x, a1y), (a2x, a2y) = prim_vec
48 area = a1x * a2y - a1y * a2x
49 b1 = (2*PI*a2y/area, -2*PI*a2x/area)
50 b2 = (-2*PI*a1y/area, 2*PI*a1x/area)
51 return [b1, b2]
54class KSpace():
55 r'''
56 Build and solve the Tight-Binding Bloch Hamiltonian :math:`H(\mathbf{k})`
57 of a periodic lattice defined by the class **lattice**.
59 Hoppings are defined between orbitals of the unit cell, separated by a
60 lattice vector :math:`\mathbf{R} = n_1\mathbf{a}_1+n_2\mathbf{a}_2`:
62 .. math::
64 H_{ij}(\mathbf{k}) = \sum_{\mathbf{R}} t_{ij}(\mathbf{R})\,
65 e^{i\mathbf{k}\cdot\mathbf{R}}
67 :param lat: **lattice** class instance. Only *unit_cell* and *prim_vec*
68 are used (the instance need not call *get_lattice*).
69 :param spin: Boolean. Default value False. If True, every site of
70 *unit_cell* carries a spin-1/2 degree of freedom (*norb* doubles to
71 ``2*len(unit_cell)``, ordered site-major: orbitals ``2*i, 2*i+1``
72 are the up/down components of site *i*). *set_onsite* and
73 *set_hopping* then accept 2x2 (spin) matrices in addition to plain
74 numbers, to build spin-orbit coupling or Zeeman terms -- see
75 :data:`PAULI` for ready-made Pauli matrices.
77 Example usage::
79 # graphene, nearest-neighbor hopping t
80 DX, DY = 0.5 * 3 ** 0.5, 0.5
81 unit_cell = [{'tag': 'a', 'r0': (0., 0.)}, {'tag': 'b', 'r0': (DX, DY)}]
82 prim_vec = [(2*DX, 0.), (DX, 1.5)]
83 lat = Lattice(unit_cell=unit_cell, prim_vec=prim_vec)
84 gra = KSpace(lat)
85 gra.set_hopping([{'i': 0, 'j': 1, 'R': (0, 0), 't': 1.},
86 {'i': 0, 'j': 1, 'R': (-1, 0), 't': 1.},
87 {'i': 0, 'j': 1, 'R': (0, -1), 't': 1.}])
88 '''
90 def __init__(self, lat: Lattice, spin: bool = False) -> None:
91 error_handling.lat(lat)
92 error_handling.boolean(spin, 'spin')
93 self.lat = lat
94 self.dim = len(lat.prim_vec)
95 self.spin = spin
96 self.n_sites = len(lat.unit_cell)
97 self.norb = 2*self.n_sites if spin else self.n_sites
98 self.tags = np.array([dic['tag'] for dic in lat.unit_cell])
99 self.onsite = np.zeros(self.norb, 'c16')
100 self._hop = [] # list of (i, j, R_cartesian (np.ndarray), t)
101 self.rec_vec = reciprocal_vectors(lat.prim_vec)
102 self.ks = np.array([]) # k-points of the last band-structure calculation
103 self.ks_dist = np.array([]) # cumulative distance along the k-path
104 self.nodes = np.array([]) # positions, along ks_dist, of the k-path nodes
105 self.en = np.array([]) # bands, shape (len(ks), norb)
107 def set_onsite(self, dict_onsite: dict[str, complex | Sequence[complex]]) -> None:
108 '''
109 Set the onsite energies, by sublattice tag.
111 :param dict_onsite: Dictionary. key: tag, val: onsite energy
112 (a plain number), or, if ``spin=True``, either a plain number
113 (applied equally to both spins) or a pair ``(E_up, E_down)`` of
114 numbers (a spin splitting, e.g. a Zeeman term along z).
116 Example usage::
118 kag.set_onsite({'a': 1., 'b': -1.})
119 # spinful: same onsite energy for both spins on 'a', a Zeeman
120 # splitting on 'b':
121 kag_spin.set_onsite({'a': 1., 'b': (1., -1.)})
122 '''
123 error_handling.set_onsite_kspace(dict_onsite, self.lat.tags, self.spin)
124 for tag, val in dict_onsite.items():
125 sites = np.where(self.tags == tag)[0]
126 if self.spin:
127 e_up, e_down = (val, val) if isinstance(val, (int, float, complex)) else val
128 self.onsite[2*sites] = e_up
129 self.onsite[2*sites + 1] = e_down
130 else:
131 self.onsite[sites] = val
133 def set_hopping(self, list_hop: list[dict]) -> None:
134 r'''
135 Set the hoppings between orbitals of the unit cell.
137 Only one representative of each hopping needs to be given: its
138 Hermitian conjugate (:math:`j\to i`, :math:`\mathbf{R}\to-\mathbf{R}`)
139 is added automatically.
141 :param list_hop: List of dictionaries with keys ('i', 'j', 'R', 't'):
143 * 'i', 'j': Positive integers. Site indices within the unit cell
144 (following the order of *unit_cell*).
145 * 'R': Tuple of one/two integers :math:`(n_1, n_2)`. Lattice vector
146 :math:`\mathbf{R}=n_1\mathbf{a}_1+n_2\mathbf{a}_2` separating the
147 two sites.
148 * 't': Complex number, or, if ``spin=True``, either a complex
149 number (spin-independent hopping) or a 2x2 complex matrix (a
150 general, possibly spin-mixing, hopping -- e.g. built from
151 :data:`PAULI` for Rashba or intrinsic spin-orbit coupling).
153 Example usage::
155 # 1D chain, nearest-neighbor hopping t between the only orbital
156 # and its right neighbor:
157 chain.set_hopping([{'i': 0, 'j': 0, 'R': (1,), 't': 1.}])
158 # spinful: spin-independent hopping t, plus a Rashba-like
159 # spin-flip term of strength alpha:
160 chain_spin.set_hopping([{'i': 0, 'j': 0, 'R': (1,),
161 't': t*PAULI['0'] + 1j*alpha*PAULI['y']}])
162 '''
163 error_handling.set_hopping_kspace(list_hop, self.n_sites, self.dim, self.spin)
164 for dic in list_hop:
165 R_cart = np.zeros(2)
166 for n, a in zip(dic['R'], self.lat.prim_vec):
167 R_cart += n * np.array(a)
168 i, j, t = dic['i'], dic['j'], dic['t']
169 if self.spin:
170 block = t*PAULI['0'] if isinstance(t, (int, float, complex)) else np.asarray(t, 'c16')
171 for a in range(2):
172 for b in range(2):
173 self._hop.append((2*i+a, 2*j+b, R_cart, block[a, b]))
174 self._hop.append((2*j+b, 2*i+a, -R_cart, np.conj(block[a, b])))
175 else:
176 self._hop.append((i, j, R_cart, t))
177 if not (i == j and not np.any(R_cart)):
178 self._hop.append((j, i, -R_cart, np.conj(t)))
180 def clear_hopping(self) -> None:
181 '''
182 Clear the hoppings set by *set_hopping*.
183 '''
184 self._hop = []
186 def get_ham(self, k: ArrayLike) -> NDArray[np.complex128]:
187 r'''
188 Get the dense Bloch Hamiltonian :math:`H(\mathbf{k})`.
190 :param k: Tuple/list/ndarray of one/two real numbers. :math:`\mathbf{k}` point,
191 in the same Cartesian frame as *prim_vec*.
193 :returns:
194 * **ham** -- Complex ndarray, shape (norb, norb).
195 '''
196 error_handling.k_vector(k, 'k', self.dim)
197 k_cart = np.zeros(2)
198 k_cart[:self.dim] = k
199 ham = np.diag(self.onsite).astype('c16')
200 for i, j, R_cart, t in self._hop:
201 ham[i, j] += t * np.exp(1j * np.dot(k_cart, R_cart))
202 return ham
204 def get_bands(
205 self, ks: ArrayLike, eigenvec: bool = False,
206 ) -> NDArray[np.float64] | tuple[NDArray[np.float64], NDArray[np.complex128]]:
207 r'''
208 Diagonalize :math:`H(\mathbf{k})` over a set of k-points.
210 :param ks: ndarray, shape (nk, dim). k-points.
211 :param eigenvec: Boolean. Default value False. If True, also return
212 the eigenvectors.
214 :returns:
215 * **en** -- Real ndarray, shape (nk, norb). Band energies, sorted ascending.
216 * **vn** -- Complex ndarray, shape (nk, norb, norb), only if *eigenvec* is True.
217 vn[k, :, n] is the nth eigenvector at ks[k].
218 '''
219 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
220 self.ks = ks
221 self.en = np.zeros((len(ks), self.norb))
222 if eigenvec:
223 vn = np.zeros((len(ks), self.norb, self.norb), 'c16')
224 for i, k in enumerate(ks):
225 ham = self.get_ham(k)
226 if eigenvec:
227 en, v = LA.eigh(ham)
228 vn[i] = v
229 else:
230 en = LA.eigvalsh(ham)
231 self.en[i] = en
232 if eigenvec:
233 return self.en, vn
234 return self.en
236 def k_path(
237 self, points: list[ArrayLike], nk: int,
238 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
239 r'''
240 Build a k-path through a list of high-symmetry points, and get the
241 associated bands.
243 :param points: List of at least two k-points (each a tuple/list of
244 one/two real numbers).
245 :param nk: Positive integer. Number of k-points per path segment.
247 :returns:
248 * **ks_dist** -- Real ndarray. Cumulative distance along the path,
249 to be used as the x-axis of a band-structure plot.
250 * **en** -- Real ndarray, shape (len(ks_dist), norb). Band energies.
251 '''
252 error_handling.k_path_points(points, self.dim)
253 error_handling.positive_int(nk, 'nk')
254 points = np.atleast_2d(np.asarray(points, dtype='f8'))
255 segments = [np.linspace(points[i], points[i+1], nk, endpoint=False)
256 for i in range(len(points) - 1)]
257 ks = np.concatenate(segments + [points[-1:]])
258 steps = np.linalg.norm(np.diff(ks, axis=0), axis=1)
259 self.ks_dist = np.concatenate([[0.], np.cumsum(steps)])
260 self.nodes = self.ks_dist[::nk][:len(points)-1].tolist() + [self.ks_dist[-1]]
261 en = self.get_bands(ks)
262 return self.ks_dist, en
264 def mesh_grid(
265 self, nk: int | tuple[int, int],
266 ) -> tuple[list[NDArray[np.float64]], NDArray[np.float64]]:
267 '''
268 Private method. Build a uniform grid of fractional coordinates
269 spanning the Brillouin zone (each in [0, 1)), and the corresponding
270 Cartesian k-points.
272 :param nk: Positive integer, or tuple of *dim* positive integers.
273 Number of k-points along each reciprocal lattice vector.
275 :returns:
276 * **fracs** -- List of *dim* real ndarrays, shape (nk1, nk2) each
277 (or (nk1,) in 1D): fractional coordinates of the grid.
278 * **ks** -- Real ndarray, shape (nk1*nk2, dim) (or (nk1, dim) in 1D).
279 '''
280 error_handling.nk(nk, self.dim)
281 if isinstance(nk, int):
282 nk = (nk,) * self.dim
283 rec_vec = [np.array(b) for b in self.rec_vec]
284 if self.dim == 1:
285 f1 = np.arange(nk[0]) / nk[0]
286 ks = f1[:, None] * rec_vec[0][None, :self.dim]
287 return [f1], ks
288 f1, f2 = np.meshgrid(np.arange(nk[0])/nk[0], np.arange(nk[1])/nk[1], indexing='ij')
289 ks = (f1.ravel()[:, None] * rec_vec[0][None, :]
290 + f2.ravel()[:, None] * rec_vec[1][None, :])
291 return [f1, f2], ks
293 def mesh_bands(self, nk: int | tuple[int, int]) -> NDArray[np.float64]:
294 '''
295 Diagonalize :math:`H(\\mathbf{k})` over a uniform mesh spanning the
296 Brillouin zone.
298 :param nk: Positive integer, or tuple of *dim* positive integers.
299 Number of k-points along each reciprocal lattice vector.
301 :returns:
302 * **en** -- Real ndarray, shape (nk1*nk2, norb). Band energies
303 over the mesh (flattened).
304 '''
305 _, ks = self.mesh_grid(nk)
306 return self.get_bands(ks)
308 def berry_curvature(
309 self, bands: int | list[int], nk: int | tuple[int, int] = 30,
310 ) -> NDArray[np.float64]:
311 r'''
312 Get the Berry curvature of a group of bands over a uniform
313 Brillouin-zone mesh, using the gauge-invariant lattice method of
314 Fukui, Hatsugai and Suzuki (J. Phys. Soc. Jpn. 74, 1674 (2005)):
315 the flux through each mesh plaquette is minus the phase of the
316 product of the (Slater-determinant) overlaps between the occupied
317 subspaces at its four corners.
319 :param bands: Positive integer, or list of positive integers. Band
320 index, or indices of a group of bands (e.g. all occupied bands
321 below a gap).
322 :param nk: Positive integer, or tuple of 2 positive integers.
323 Default value 30. Number of k-points along each reciprocal
324 lattice vector.
326 :returns:
327 * **curv** -- Real ndarray, shape (nk1, nk2). Berry curvature
328 (flux through each plaquette, in radians). Summing *curv* and
329 dividing by :math:`2\pi` gives the Chern number, see
330 *chern_number*.
331 '''
332 error_handling.dim_2(self.dim)
333 if isinstance(bands, int):
334 bands = [bands]
335 error_handling.band_indices(bands, self.norb)
336 if isinstance(nk, int):
337 nk = (nk, nk)
338 error_handling.nk(nk, 2)
339 n1, n2 = nk
340 _, ks = self.mesh_grid(nk)
341 ks = ks.reshape(n1, n2, 2)
342 v = np.zeros((n1, n2, self.norb, len(bands)), 'c16')
343 for i1 in range(n1):
344 for i2 in range(n2):
345 _, vn = LA.eigh(self.get_ham(ks[i1, i2]))
346 v[i1, i2] = vn[:, bands]
347 curv = np.zeros((n1, n2))
348 for i1 in range(n1):
349 for i2 in range(n2):
350 v1 = v[i1, i2]
351 v2 = v[(i1+1) % n1, i2]
352 v3 = v[(i1+1) % n1, (i2+1) % n2]
353 v4 = v[i1, (i2+1) % n2]
354 link = (np.linalg.det(v1.conj().T @ v2)
355 * np.linalg.det(v2.conj().T @ v3)
356 * np.linalg.det(v3.conj().T @ v4)
357 * np.linalg.det(v4.conj().T @ v1))
358 curv[i1, i2] = -np.angle(link)
359 return curv
361 def chern_number(self, bands: int | list[int], nk: int | tuple[int, int] = 30) -> float:
362 r'''
363 Get the Chern number of a group of bands:
365 .. math::
367 C = \frac{1}{2\pi}\int_{BZ} \Omega(\mathbf{k})\, d^2k
369 an integer (up to the numerical precision set by *nk*) for a group
370 of bands that is isolated from the rest of the spectrum by a gap
371 everywhere in the Brillouin zone. See *berry_curvature*.
373 :param bands: Positive integer, or list of positive integers. Band
374 index, or indices of a group of bands (e.g. all occupied bands
375 below a gap).
376 :param nk: Positive integer, or tuple of 2 positive integers.
377 Default value 30. Number of k-points along each reciprocal
378 lattice vector.
380 :returns:
381 * **chern** -- Real number, close to an integer.
382 '''
383 return self.berry_curvature(bands, nk).sum() / (2*PI)
385 def plot_dos(
386 self,
387 nk: int | tuple[int, int] = 30,
388 broadening: float = 0.05,
389 kernel: str = 'gaussian',
390 e_grid: ArrayLike | None = None,
391 fs: float = 20,
392 lw: float = 2.,
393 figsize: tuple[float, float] | None = None,
394 ) -> Figure:
395 '''
396 Plot the (broadened) density of states, obtained by diagonalizing
397 :math:`H(\\mathbf{k})` over a uniform Brillouin-zone mesh -- see
398 *tbkit.dos.density_of_states*.
400 :param nk: Positive integer, or tuple of *dim* positive integers.
401 Default value 30. Number of k-points along each reciprocal
402 lattice vector.
403 :param broadening: Positive real number. Default value 0.05. Kernel width.
404 :param kernel: String. Default value 'gaussian'. 'gaussian' or 'lorentzian'.
405 :param e_grid: Real ndarray. Default value None. Energies at which to
406 evaluate the density of states.
407 :param fs: Positive number. Default value 20. Fontsize.
408 :param lw: Positive number. Default value 2. Linewidth.
409 :param figsize: Tuple. Default value None. Figure size.
411 :returns:
412 * **fig** -- Figure.
413 '''
414 error_handling.positive_real(fs, 'fs')
415 error_handling.positive_real(lw, 'lw')
416 error_handling.tuple_2elem(figsize, 'figsize')
417 en = self.mesh_bands(nk)
418 e_grid, rho = dos.density_of_states(en, e_grid=e_grid,
419 broadening=broadening, kernel=kernel)
420 fig, ax = plt.subplots(figsize=figsize)
421 ax.plot(e_grid, rho, 'b', lw=lw)
422 ax.fill_between(e_grid, rho, color='b', alpha=0.2)
423 ax.set_xlim([e_grid[0], e_grid[-1]])
424 ax.set_ylim([0., None])
425 ax.set_title('Density of states', fontsize=fs)
426 ax.set_xlabel('$E$', fontsize=fs)
427 ax.set_ylabel(r'$\rho(E)$', fontsize=fs)
428 for label in ax.xaxis.get_majorticklabels():
429 label.set_fontsize(fs)
430 for label in ax.yaxis.get_majorticklabels():
431 label.set_fontsize(fs)
432 fig.set_layout_engine('tight')
433 plt.draw()
434 return fig
436 def plot_bands(
437 self,
438 node_labels: list[str] | None = None,
439 fs: float = 20,
440 lw: float = 2.,
441 ms: float = 0.,
442 c: str = 'b',
443 lims: tuple[float, float] | None = None,
444 figsize: tuple[float, float] | None = None,
445 ) -> Figure:
446 '''
447 Plot the band structure computed by *k_path* or *get_bands*.
449 :param node_labels: List of strings. Default value None. Labels of the
450 high-symmetry points passed to *k_path*.
451 :param fs: Positive number. Default value 20. Fontsize.
452 :param lw: Positive number. Default value 2. Linewidth.
453 :param ms: Positive number. Default value 0. Marker size.
454 :param c: Default value 'b'. Line color.
455 :param lims: List. Default value None. Energy plot limits.
456 :param figsize: Tuple. Default value None. Figure size.
458 :returns:
459 * **fig** -- Figure.
460 '''
461 error_handling.empty_ndarray(self.en, 'get_bands or k_path')
462 error_handling.positive_real(fs, 'fs')
463 error_handling.positive_real(lw, 'lw')
464 error_handling.lims(lims)
465 error_handling.tuple_2elem(figsize, 'figsize')
466 fig, ax = plt.subplots(figsize=figsize)
467 for n in range(self.norb):
468 ax.plot(self.ks_dist, self.en[:, n], c=c, lw=lw, marker='o', ms=ms)
469 for node in self.nodes:
470 ax.axvline(node, color='k', lw=0.5)
471 ax.set_xlim([self.ks_dist[0], self.ks_dist[-1]])
472 if lims is not None:
473 ax.set_ylim(lims)
474 if node_labels is not None:
475 error_handling.ndarray(np.array(node_labels), 'node_labels', len(self.nodes))
476 ax.set_xticks(self.nodes)
477 ax.set_xticklabels(node_labels, fontsize=fs)
478 ax.set_ylabel('$E$', fontsize=fs)
479 for label in ax.yaxis.get_majorticklabels():
480 label.set_fontsize(fs)
481 fig.set_layout_engine('tight')
482 plt.draw()
483 return fig
485 def show(self) -> None:
486 '''
487 Emulate Matplotlib method plt.show().
488 '''
489 plt.show()
492def ribbon(
493 lat: Lattice,
494 list_hop: list[dict],
495 width: int,
496 direction: int = 1,
497 onsite: dict | None = None,
498 spin: bool = False,
499) -> KSpace:
500 r'''
501 Cut a ribbon out of a 2D periodic model: periodic along one primitive
502 vector, finite (open boundary, *width* unit cells) along the other.
503 This is the standard way to see edge states in a band structure (e.g.
504 the zero-energy edge band of a zigzag graphene ribbon, or the helical
505 edge states of a Kane-Mele ribbon).
507 :param lat: **Lattice** class instance (2D, i.e. two primitive
508 vectors). Only *unit_cell* and *prim_vec* are used.
509 :param list_hop: List of dictionaries, in the same format passed to
510 *KSpace.set_hopping* -- the hoppings of the periodic (2D) model
511 that the ribbon is cut from.
512 :param width: Positive integer. Number of unit cells across the ribbon.
513 :param direction: 0 or 1. Default value 1. Which primitive vector
514 (``lat.prim_vec[direction]``) becomes finite; the other stays
515 periodic.
516 :param onsite: Dictionary. Default value None. Onsite energies, in the
517 same format passed to *KSpace.set_onsite* -- applied identically
518 on every row of the ribbon.
519 :param spin: Boolean. Default value False. See *KSpace*.
521 :returns:
522 * **rib** -- **KSpace** instance, 1D-periodic, with
523 ``width * len(lat.unit_cell)`` sites (each site of *lat*,
524 repeated once per row across the ribbon; row *w*'s copy of site
525 *i* is orbital ``w*len(lat.unit_cell) + i``).
527 Example usage::
529 # zigzag graphene ribbon, 20 unit cells wide
530 list_hop = [{'i': 0, 'j': 1, 'R': (0, 0), 't': 1.},
531 {'i': 0, 'j': 1, 'R': (-1, 0), 't': 1.},
532 {'i': 0, 'j': 1, 'R': (0, -1), 't': 1.}]
533 rib = ribbon(lat, list_hop, width=20)
534 '''
535 error_handling.lat(lat)
536 error_handling.dim_2(len(lat.prim_vec))
537 error_handling.positive_int(width, 'width')
538 error_handling.direction(direction)
539 periodic = 1 - direction
540 n_sites = len(lat.unit_cell)
541 a_dir = np.array(lat.prim_vec[direction])
542 new_unit_cell = []
543 for w in range(width):
544 for dic in lat.unit_cell:
545 r0 = np.array(dic['r0']) + w*a_dir
546 new_unit_cell.append({'tag': dic['tag'], 'r0': (float(r0[0]), float(r0[1]))})
547 new_lat = Lattice(unit_cell=new_unit_cell, prim_vec=[lat.prim_vec[periodic]])
548 rib = KSpace(new_lat, spin=spin)
549 new_list_hop = []
550 for dic in list_hop:
551 w2_shift = dic['R'][direction]
552 n_periodic = dic['R'][periodic]
553 for w in range(width):
554 w2 = w + w2_shift
555 if 0 <= w2 < width:
556 new_list_hop.append({'i': w*n_sites + dic['i'],
557 'j': w2*n_sites + dic['j'],
558 'R': (n_periodic,),
559 't': dic['t']})
560 rib.set_hopping(new_list_hop)
561 if onsite is not None:
562 rib.set_onsite(onsite)
563 return rib