Coverage for tbkit/kspace.py: 100%

218 statements  

« prev     ^ index     » next       coverage.py v7.16.0, created at 2026-09-22 13:16 +0100

1from __future__ import annotations 

2 

3from typing import Sequence 

4 

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 

13 

14 

15PI = np.pi 

16 

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} 

26 

27 

28################################# 

29# CLASS KSPACE 

30################################# 

31 

32 

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}`. 

37 

38 :param prim_vec: List of one/two tuples. Primitive vectors (see class **lattice**). 

39 

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] 

52 

53 

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**. 

58 

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`: 

61 

62 .. math:: 

63 

64 H_{ij}(\mathbf{k}) = \sum_{\mathbf{R}} t_{ij}(\mathbf{R})\, 

65 e^{i\mathbf{k}\cdot\mathbf{R}} 

66 

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. 

76 

77 Example usage:: 

78 

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 ''' 

89 

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) 

106 

107 def set_onsite(self, dict_onsite: dict[str, complex | Sequence[complex]]) -> None: 

108 ''' 

109 Set the onsite energies, by sublattice tag. 

110 

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). 

115 

116 Example usage:: 

117 

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 

132 

133 def set_hopping(self, list_hop: list[dict]) -> None: 

134 r''' 

135 Set the hoppings between orbitals of the unit cell. 

136 

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. 

140 

141 :param list_hop: List of dictionaries with keys ('i', 'j', 'R', 't'): 

142 

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). 

152 

153 Example usage:: 

154 

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))) 

179 

180 def clear_hopping(self) -> None: 

181 ''' 

182 Clear the hoppings set by *set_hopping*. 

183 ''' 

184 self._hop = [] 

185 

186 def get_ham(self, k: ArrayLike) -> NDArray[np.complex128]: 

187 r''' 

188 Get the dense Bloch Hamiltonian :math:`H(\mathbf{k})`. 

189 

190 :param k: Tuple/list/ndarray of one/two real numbers. :math:`\mathbf{k}` point, 

191 in the same Cartesian frame as *prim_vec*. 

192 

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 

203 

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. 

209 

210 :param ks: ndarray, shape (nk, dim). k-points. 

211 :param eigenvec: Boolean. Default value False. If True, also return 

212 the eigenvectors. 

213 

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 

235 

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. 

242 

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. 

246 

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 

263 

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. 

271 

272 :param nk: Positive integer, or tuple of *dim* positive integers. 

273 Number of k-points along each reciprocal lattice vector. 

274 

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 

292 

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. 

297 

298 :param nk: Positive integer, or tuple of *dim* positive integers. 

299 Number of k-points along each reciprocal lattice vector. 

300 

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) 

307 

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. 

318 

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. 

325 

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 

360 

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: 

364 

365 .. math:: 

366 

367 C = \frac{1}{2\pi}\int_{BZ} \Omega(\mathbf{k})\, d^2k 

368 

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*. 

372 

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. 

379 

380 :returns: 

381 * **chern** -- Real number, close to an integer. 

382 ''' 

383 return self.berry_curvature(bands, nk).sum() / (2*PI) 

384 

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*. 

399 

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. 

410 

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 

435 

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*. 

448 

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. 

457 

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 

484 

485 def show(self) -> None: 

486 ''' 

487 Emulate Matplotlib method plt.show(). 

488 ''' 

489 plt.show() 

490 

491 

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). 

506 

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*. 

520 

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``). 

526 

527 Example usage:: 

528 

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