Coverage for tbkit/system.py: 100%

1from __future__ import annotations 

2 

3from typing import Callable 

4 

5import numpy as np 

6from numpy.typing import NDArray 

7import scipy.sparse as sparse 

8import scipy.linalg as LA 

9import numpy.random as rand 

10import numpy.char as npc 

11from math import sin, cos 

12import tbkit.error_handling as error_handling 

13from tbkit.lattice import Lattice, COOR_DTYPE 

14 

15 

16PI = np.pi 

17ATOL = 1e-3 

18HOP_DTYPE = [('n', 'u2'), ('i', 'u4'), ('j', 'u4'), 

19 ('ang', 'f8'), ('tag', 'U2'), ('t', 'c16')] 

20 

21 

22class System(): 

23 ''' 

24 Solve the Tight-Binding eigenvalue problem of a lattice defined  

25 by the class **lattice**. 

26 

27 :param lat: **lattice** class instance. 

28 ''' 

29 

30 def __init__(self, lat: Lattice) -> None: 

31 error_handling.lat(lat) 

32 self.lat = lat 

33 self.sites = self.lat.sites # used to check if sites changes 

34 self.coor_hop = np.array([], dtype=COOR_DTYPE) 

35 self.vec_hop = np.array([], dtype=[('dis', 'f8'), ('ang', 'f8')]) # Hopping distances and angles 

36 self.dist_uni = np.array([], 'f8') # Different hopping distances 

37 self.store_hop = {} # Store the relevant hoppings (dynamic programming) 

38 self.hop = np.array([], dtype=HOP_DTYPE) # Hoppings to build-up the Hamiltonian 

39 self.onsite = np.array([], 'c16') # Onsite energies 

40 self.ham = sparse.csr_matrix(([], ([], [])), shape=(self.lat.sites, self.lat.sites)) # Hamiltonian 

41 self.en = np.array([], 'c16') # Eigenenergies 

42 self.rn = np.array([], 'c16') # Right eigenvectors: H |rn> = en |rn> 

43 self.ln = np.array([], 'c16') # Left eigenvectors: <ln| H = en <ln| 

44 self.intensity = np.array([], 'f8') # Intensities (|rn|**2) 

45 self.pola = np.array([], 'f8') # sublattices polarisation (|rn^{(S)}|**2) 

46 self.petermann = np.array([], 'f8') # Inverse Participation Ratio 

47 self.nmax = 0 # number of different hoppings 

48 

49 def clear_hopping(self) -> None: 

50 ''' 

51 Clear structured array *hop*. 

52 ''' 

53 self.hop = np.array([], dtype=HOP_DTYPE) 

54 

55 def get_distances(self) -> None: 

56 ''' 

57 Private method. 

58 Get distances and angles of the edges. 

59 ''' 

60 error_handling.sites(self.lat.sites) 

61 dif_x = self.lat.coor['x'] - self.lat.coor['x'].reshape(self.lat.sites, 1) 

62 dif_y = self.lat.coor['y'] - self.lat.coor['y'].reshape(self.lat.sites, 1) 

63 dist = np.sqrt(dif_x ** 2 + dif_y ** 2) 

64 ang = (180 / PI * np.arctan2(dif_y, dif_x)) 

65 self.vec_hop = np.zeros(dist.shape, dtype=[('dis', 'f8'), ('ang', 'f8')]) 

66 self.vec_hop['dis'] = dist 

67 self.vec_hop['ang'] = ang 

68 self.dist_uni = np.unique(self.vec_hop['dis'].round(4)) 

69 

70 def print_distances(self, n: int = 1) -> None: 

71 r''' 

72 Print distances and positive angles (in degrees) :math:`\phi_+\in[0, 180)` 

73 of the nth shortest edges. Negative angles are given by: 

74 :math:`\phi_-= \phi_+-180` and :math:`\phi_+\in[-180, 0)`. 

75 

76 :param n: Positive integer. Number of shortest edges. 

77 ''' 

78 error_handling.sites(self.lat.sites) 

79 self.get_distances() 

80 self.nmax = len(self.dist_uni) - 1 

81 error_handling.positive_int_lim(n, 'n', self.nmax) 

82 print('\n{} different distances between sites:'.format(self.nmax)) 

83 print('\nDistances between sites:') 

84 for i, d in enumerate(self.dist_uni[1: n+1]): 

85 if i == 0: 

86 hop_name = 'st' 

87 elif i == 1: 

88 hop_name = 'nd' 

89 elif i == 2: 

90 hop_name = 'rd' 

91 else: 

92 hop_name = 'th' 

93 print('{}{} hopping, length: {:.3f}'.format(i+1, hop_name, d)) 

94 print('\twith positive angles:') 

95 positive_ang = self.vec_hop['ang'][np.isclose(d, self.vec_hop['dis'], atol=ATOL) & 

96 (self.vec_hop['ang'] >= 0.) & 

97 (self.vec_hop['ang'] < 180.)] 

98 print('\t', np.unique(positive_ang.round(4))) 

99 

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

101 ''' 

102 Set onsite energies. 

103 

104 :param on: Array. Sublattice onsite energies. 

105 

106 Example usage:: 

107 

108 # Line-Centered Square lattice 

109 sys.set_onsite({'a': -1j, 'b': -2j}) 

110 ''' 

111 error_handling.sites(self.lat.sites) 

112 error_handling.set_onsite(dict_onsite, self.lat.tags) 

113 self.onsite = np.zeros(self.lat.sites, 'c16') 

114 for tag, on in dict_onsite.items(): 

115 self.onsite[self.lat.coor['tag'] ==tag] = on 

116 

117 def fill_store_hop(self, n: int) -> None: 

118 ''' 

119 Private method. 

120 

121 Store in *store_hop* indices (with :math:`i < j`), positive angles, and tags 

122 of a given type of hopping. 

123 ''' 

124 ind = np.argwhere(np.isclose(self.dist_uni[n], self.vec_hop['dis'], atol=ATOL)) 

125 ind_up = ind[ind[:, 1] > ind[:, 0]] 

126 hop = np.zeros(len(ind_up), dtype=HOP_DTYPE[:-1]) 

127 hop['i'] = ind_up[:, 0] 

128 hop['j'] = ind_up[:, 1] 

129 hop['ang'] = self.vec_hop['ang'][ind_up[:, 0], ind_up[:, 1]] 

130 hop['tag'] = npc.add(self.lat.coor['tag'][ind_up[:, 0]], 

131 self.lat.coor['tag'][ind_up[:, 1]]) 

132 self.store_hop[n] = hop 

133 

134 def set_hopping(self, list_hop: list[dict], upper_part: bool = True) -> None: 

135 r''' 

136 Set lattice hoppings. 

137 

138 :param list_hop: List of Dictionaries. 

139 Dictionary with keys ('n', 'ang', 'tag', 't') where: 

140 

141 * 'n' Positive integer, type of hoppings: 

142 

143 * 'n': 1 for nearest neighbours. 

144 * 'n': 2 for next-nearest neighbours.  

145 * 'n': 3 for next-next-nearest neighbours.  

146 * etc... 

147 

148 * 'ang' value, float, angle, in deg, of the hoppings. (optional). 

149 

150 Hopping angles are given by the method *print_distances*. 

151 

152 * If :math:`ang \in[0, 180)`, fill the Hamiltonian upper part. 

153 * If :math:`ang \in[-180, 0)`, fill the Hamiltonian lower part. 

154 

155 * 'tag' string of length 2 (optional). 

156 

157 Hopping tags. 

158 

159 * 't' Complex number. 

160 

161 Hopping value. 

162 

163 :param upper_part: Boolean. Default value True. 

164 

165 * True get hoppings with (:math:`i<j`) *i.e.* fill the Hamiltonian upper part. 

166 * False get hoppings with (:math:`i>j`) *i.e.* fill the Hamiltonian lower part. 

167 

168 Example usage:: 

169 

170 # fill upper part: 

171 sys.set_hopping([{'n': 1, t: 1.}]) 

172 # fill lower part: 

173 sys.set_hopping([{'n': 1, t: 1.}], upper_part=False) 

174 # fill upper part: specifying the angles: 

175 sys.set_hopping([{'n': 1, 'ang': 0., t: 1.}, {'n': 1, 'ang': 90, t: 2.}]) 

176 # fill lower part: 

177 sys.set_hopping([{'n': 1, 'ang': -180., t: 1.}, {'n': 1, 'ang': -90, t: 2.}], upper_part=False) 

178 # fill upper part: specifying the tags: 

179 sys.set_hopping([{'n': 1, 'tag': 'ab', t: 1.}, {'n': 1, 'tag': 'ba', t: 2.}]) 

180 # fill lower part: 

181 sys.set_hopping([{'n': 1, 'tag': 'ab', t: 1.}, {'n': 1, 'tag': 'ba', t: 2.}], upper_part=False) 

182 # fill upper part: specifying the angles and tags: 

183 sys.set_hopping([{'n': 1, 'ang': 0., 'tag': 'ab', t: 1.},  

184 {'n': 1, 'ang': 0., 'tag': 'ba', t: 2.}, 

185 {'n': 1, 'ang': 90., 'tag': 'ab', t: 3.},  

186 {'n': 1, 'ang': 90., 'tag': 'ba', t: 4.}]) 

187 # fill lower part: 

188 sys.set_hopping([{'n': 1, 'ang': 0., 'tag': 'ab', t: 1.},  

189 {'n': 1, 'ang': 0., 'tag': 'ba', t: 2.}, 

190 {'n': 1, 'ang': 90., 'tag': 'ab', t: 3.},  

191 {'n': 1, 'ang': 90., 'tag': 'ba', t: 4.}]), upper_part=False) 

192 

193 .. note:: 

194 

195 A Hermitian hopping matrix can be build-up only using  

196 its upper part OR only using its lower part. The full matrix is then 

197 automatic built by Hermitian conjugaison.  

198  

199 If both upper AND lower parts are used to build up the hopping matrix. 

200 non Hermitian conjugaison is not performed *i.e.* non-Hermitian hopping matrix 

201 can be built. 

202 ''' 

203 error_handling.sites(self.lat.sites) 

204 error_handling.boolean(upper_part, 'upper_part') 

205 self.get_distances() 

206 self.nmax = len(self.dist_uni) - 1 

207 error_handling.set_hopping(list_hop, self.nmax) 

208 list_n = np.unique([dic['n'] for dic in list_hop]) 

209 # fill, if needed self.store_hop 

210 self.check_sites() 

211 for n in list_n: 

212 if n not in self.store_hop: 

213 self.fill_store_hop(n) 

214 # fill self.hop 

215 for dic in list_hop: 

216 if len(dic) == 2: 

217 size = len(self.store_hop[dic['n']]) 

218 if upper_part: 

219 mask = (self.hop['n'] == dic['n']) & (self.hop['i'] < self.hop['j']) 

220 else: 

221 mask = (self.hop['n'] == dic['n']) & (self.hop['i'] > self.hop['j']) 

222 if np.sum(mask): 

223 self.hop = self.hop[np.logical_not(mask)] 

224 ind = np.ones(size, bool) 

225 hop = self.set_given_hopping(dic['n'], size, dic, ind, upper_part=upper_part) 

226 elif len(dic) == 3 and 'ang' in dic: 

227 error_handling.angle(dic['ang'], np.unique(self.store_hop[dic['n']]['ang']), upper_part) 

228 if dic['ang'] >= 0: 

229 ang_store = dic['ang'] 

230 else: 

231 ang_store = dic['ang'] + 180. 

232 size = np.sum(np.isclose(ang_store, self.store_hop[dic['n']]['ang'], atol=ATOL)) 

233 mask = (self.hop['n'] == dic['n']) & np.isclose(self.hop['ang'], dic['ang'], atol=ATOL) 

234 if np.sum(mask): 

235 self.hop = self.hop[np.logical_not(mask)] 

236 ind = np.isclose(ang_store, self.store_hop[dic['n']]['ang'], atol=ATOL) 

237 error_handling.index(ind, dic) 

238 hop = self.set_given_hopping(dic['n'], size, dic, ind, upper_part=upper_part) 

239 elif len(dic) == 3 and 'tag' in dic: 

240 if upper_part: 

241 tag_store = dic['tag'] 

242 else: 

243 tag_store = dic['tag'][::-1] 

244 size = np.sum(self.store_hop[dic['n']]['tag'] == tag_store) 

245 mask = (self.hop['n'] == dic['n']) & (self.hop['tag'] == dic['tag']) 

246 if upper_part: 

247 mask = self.hop['n'] == dic['n'] & (self.hop['tag'] == dic['tag']) & (self.hop['i'] < self.hop['j']) 

248 else: 

249 mask = self.hop['n'] == dic['n'] & (self.hop['tag'] == dic['tag']) & (self.hop['i'] > self.hop['j']) 

250 if np.sum(mask): 

251 self.hop = self.hop[np.logical_not(mask)] 

252 ind = self.store_hop[dic['n']]['tag'] == tag_store 

253 error_handling.index(ind, dic) 

254 hop = self.set_given_hopping(dic['n'], size, dic, ind, upper_part=upper_part) 

255 else: 

256 error_handling.angle(dic['ang'], np.unique(self.store_hop[dic['n']]['ang']), upper_part=upper_part) 

257 error_handling.tag(dic['tag'], np.unique(self.store_hop[dic['n']]['tag'])) 

258 if dic['ang'] >= 0: 

259 ang_store = dic['ang'] 

260 else: 

261 ang_store = dic['ang'] + 180. 

262 if upper_part: 

263 tag_store = dic['tag'] 

264 else: 

265 tag_store = dic['tag'][::-1] 

266 size = np.sum((self.store_hop[dic['n']]['tag'] == tag_store) & 

267 (np.isclose(ang_store, self.store_hop[dic['n']]['ang'], atol=ATOL))) 

268 bool1 = (self.hop['n'] == dic['n']) & (self.hop['tag'] == dic['tag']) 

269 bool2 = np.isclose(self.hop['ang'], dic['ang'], atol=ATOL) 

270 mask = bool1 & bool2 

271 if np.sum(mask): 

272 self.hop = self.hop[np.logical_not(mask)] 

273 ind = ((self.store_hop[dic['n']]['tag'] == tag_store) & 

274 (np.isclose(ang_store, self.store_hop[dic['n']]['ang'], atol=1))) 

275 error_handling.index(ind, dic) 

276 hop = self.set_given_hopping(dic['n'], size, dic, ind, upper_part=upper_part) 

277 self.hop = np.concatenate([self.hop, hop]) 

278 

279 def check_sites(self) -> None: 

280 ''' 

281 Private method. 

282 Check if the number of sites was changed after calling the  

283 method system.set_hopping(). 

284 ''' 

285 if self.sites != self.lat.sites: 

286 self.store_hop = {} 

287 self.sites = self.lat.sites 

288 

289 def set_given_hopping( 

290 self, n: int, size: int, dic: dict, mask: NDArray, upper_part: bool, 

291 ) -> NDArray: 

292 ''' 

293 Private method. 

294 Fill self.hop.  

295 

296 :param n: Integer. Hopping type. 

297 :param size: Integer. Number of hoppings. 

298 :param doc: Dictionary. Hopping dictionary. 

299 :param mask: np.ndarray. Mask. 

300 :param upper_part: Boolean. If True, self.hop['i'] < self.hop['j']. 

301 ''' 

302 hop = np.empty(size, dtype=HOP_DTYPE) 

303 hop['n'] = dic['n'] 

304 hop['t'] = dic['t'] 

305 if upper_part: 

306 hop['i'] = self.store_hop[n]['i'][mask] 

307 hop['j'] = self.store_hop[n]['j'][mask] 

308 hop['ang'] = self.store_hop[n]['ang'][mask] 

309 hop['tag'] = self.store_hop[n]['tag'][mask] 

310 else: 

311 hop['i'] = self.store_hop[n]['j'][mask] 

312 hop['j'] = self.store_hop[n]['i'][mask] 

313 hop['ang'] = self.store_hop[n]['ang'][mask] - 180 

314 hop['tag'] = npc.add(self.lat.coor['tag'][hop['i']], 

315 self.lat.coor['tag'][hop['j']]) 

316 return hop 

317 

318 def set_hopping_manual(self, dict_hop: dict[tuple[int, int], complex], upper_part: bool = True) -> None: 

319 ''' 

320 Set hoppings manually. 

321 

322 :param dict_hop: Dictionary of hoppings. 

323 key: hopping indices, val: hopping values. 

324 

325 :parameter upper_part: Boolean.  

326 

327 * True, fill the Hamiltonian upper part. 

328 * False, fill the Hamiltonian lower part.  

329 ''' 

330 hop = np.zeros(len(dict_hop), dtype=HOP_DTYPE) 

331 i = [h[0] for h in dict_hop.keys()] 

332 j = [h[1] for h in dict_hop.keys()] 

333 t = [val for val in dict_hop.values()] 

334 hop['i'], hop['j']= i, j 

335 hop['t'] = t 

336 hop['tag'] = npc.add(self.lat.coor['tag'][i], 

337 self.lat.coor['tag'][j]) 

338 ang = 180 / PI * np.arctan2(self.lat.coor['y'][j]-self.lat.coor['y'][i], 

339 self.lat.coor['x'][j]-self.lat.coor['x'][i]) 

340 if upper_part: 

341 ang[ang < 0] += 180 

342 else: 

343 ang[ang >= 0] -= 180 

344 hop['ang'] = ang 

345 self.hop = np.concatenate([self.hop, hop]) 

346 

347 def set_hopping_dis(self, alpha: complex) -> None: 

348 ''' 

349 Set uniform hopping disorder.  

350 

351 :param alpha: Complex or Real number. Disorder stength. 

352 

353 Example usage:: 

354 

355 sys.set_hopping_dis(alpha=0.1) 

356 

357 ''' 

358 error_handling.empty_hop(self.hop) 

359 error_handling.number(alpha, 'alpha') 

360 self.hop['t'] *= 1. + alpha * rand.uniform(-1., 1., len(self.hop)) 

361 

362 def set_peierls_phase(self, phase: Callable[..., NDArray]) -> None: 

363 r''' 

364 Apply the Peierls substitution to the existing hoppings, to 

365 capture the effect of an orbital magnetic field: 

366 

367 .. math:: 

368 

369 t_{ij} \to t_{ij}\, e^{i\phi_{ij}}\, ,\quad 

370 \phi_{ij} = \frac{2\pi}{\Phi_0}\int_{\mathbf{r}_i}^{\mathbf{r}_j} 

371 \mathbf{A}\cdot d\mathbf{l} 

372 

373 where :math:`\mathbf{A}` is the vector potential, integrated along 

374 the straight bond from site :math:`i` to site :math:`j`. 

375 

376 Must be called after *set_hopping* / *set_hopping_manual* (it 

377 rescales the existing hoppings *in place*) and before *get_ham*. 

378 For a uniform perpendicular field, use the convenience method 

379 *set_magnetic_field* instead. 

380 

381 :param phase: Callable. ``phase(xi, yi, xj, yj)`` returns 

382 :math:`\phi_{ij}`, the (real-valued) Peierls phase for the bond 

383 from :math:`(x_i, y_i)` to :math:`(x_j, y_j)`. Called with 

384 Numpy arrays (one value per hopping in *sys.hop*). 

385 

386 .. note:: 

387 

388 *get_ham* automatically assigns the reversed bond its complex 

389 conjugate, so the Hamiltonian stays Hermitian as long as 

390 *phase* is antisymmetric under swapping :math:`i` and 

391 :math:`j` -- true for the line integral of any vector 

392 potential, since reversing the integration path negates it. 

393 

394 Example usage:: 

395 

396 # Peierls phase from a uniform field via the Landau gauge 

397 # A = (0, B x): captures the same physics as set_magnetic_field, 

398 # just in a different (equally valid) gauge. 

399 B = 0.05 

400 sys.set_peierls_phase(lambda xi, yi, xj, yj: B * (xj - xi) * (xi + xj) / 2) 

401 ''' 

402 error_handling.empty_hop(self.hop) 

403 error_handling.is_callable(phase, 'phase') 

404 xi = self.lat.coor['x'][self.hop['i']] 

405 yi = self.lat.coor['y'][self.hop['i']] 

406 xj = self.lat.coor['x'][self.hop['j']] 

407 yj = self.lat.coor['y'][self.hop['j']] 

408 self.hop['t'] = self.hop['t'] * np.exp(1j * phase(xi, yi, xj, yj)) 

409 

410 def set_magnetic_field(self, alpha: float) -> None: 

411 r''' 

412 Set a uniform perpendicular magnetic field via the Peierls 

413 substitution (see *set_peierls_phase*), using the symmetric gauge 

414 :math:`\mathbf{A} = \frac{B}{2}(-y, x)`: 

415 

416 .. math:: 

417 

418 \phi_{ij} = \pi\alpha\,(x_iy_j - x_jy_i) 

419 

420 :param alpha: Real number. Flux density :math:`B/\Phi_0`, in flux 

421 quanta per unit area (in the lattice's length units) -- *i.e.* 

422 the flux through a region of area :math:`S` is 

423 :math:`\alpha S` flux quanta. 

424 

425 Example usage:: 

426 

427 # one flux quantum per 100 unit cells of a lattice with 

428 # lattice constant 1: 

429 sys.set_magnetic_field(alpha=0.01) 

430 ''' 

431 error_handling.real_number(alpha, 'alpha') 

432 self.set_peierls_phase(lambda xi, yi, xj, yj: PI * alpha * (xi*yj - xj*yi)) 

433 

434 def set_onsite_dis(self, alpha: complex) -> None: 

435 ''' 

436 Set uniform onsite disorder.  

437 

438 :param alpha: Complex or Real number. Disorder stength. 

439 

440 Example usage:: 

441 

442 sys.set_onsite_dis(alpha=0.1) 

443 

444 ''' 

445 error_handling.empty_onsite(self.onsite) 

446 error_handling.number(alpha, 'alpha') 

447 self.onsite += alpha * rand.uniform(-1., 1., self.lat.sites) 

448 

449 def set_onsite_def(self, onsite_def: dict[int, complex]) -> None: 

450 ''' 

451 Set specific onsite energies. 

452 

453 :param onsite_def: Dictionary.  

454 key: site indices, val: onsite values. 

455 

456 Example usage:: 

457 

458 set_onsite_def(0: 1., 1: -1j) 

459 ''' 

460 error_handling.empty_onsite(self.onsite) 

461 error_handling.set_onsite_def(onsite_def, self.lat.sites) 

462 for i, o in onsite_def.items(): 

463 self.onsite[i] = o 

464 

465 def set_hopping_def(self, hopping_def: dict[tuple[int, int], complex]) -> None: 

466 ''' 

467 Set specific hoppings.  

468 

469 :param hopping_def: Dictionary of hoppings.  

470 key: hopping indices, val: hopping values.  

471 

472 Example usage:: 

473 

474 sys.set_hopping_def({(0, 1): 1., (1, 2): -1j}) 

475 ''' 

476 error_handling.empty_hop(self.hop) 

477 error_handling.set_hopping_def(self.hop, hopping_def, self.lat.sites) 

478 for key, val in hopping_def.items(): 

479 cond = (self.hop['i'] == key[0]) & (self.hop['j'] == key[1]) 

480 self.hop['t'][cond] = val 

481 self.hop['ang'] = self.vec_hop['ang'][key[0], key[1]] 

482 self.hop['tag'] = npc.add(self.lat.coor['tag'][key[0]], 

483 self.lat.coor['tag'][key[1]]) 

484 

485 def set_new_hopping(self, list_hop: list[dict], ind: NDArray) -> None: 

486 ''' 

487 Private method. 

488 Set new hoppings. 

489 

490 :param list_hop: List of Dictionary (see set_hopping definition). 

491 :param ind: List. List of indices. 

492 ''' 

493 for dic in list_hop: 

494 if len(dic) == 2: 

495 self.hop['t'][ind] = dic['t'] 

496 elif len(dic) == 3 and 'ang' in dic: 

497 self.hop['t'][ind & (self.hop['ang'] == dic['ang'])] = dic['t'] 

498 elif len(dic) == 3 and 'tag' in dic: 

499 self.hop['t'][ind & (self.hop['tag'] == dic['tag'])] = dic['t'] 

500 else: 

501 self.hop['t'][ind & (self.hop['tag'] == dic['tag']) 

502 & (self.hop['ang'] == dic['ang'])] = dic['t'] 

503 

504 def find_square(self, xlims: tuple[float, float], ylims: tuple[float, float]) -> NDArray: 

505 ''' 

506 Private method. 

507 Find hoppings within the square. 

508 

509 :param xlims: List or Tuple. :math:`x` interval. 

510 :param ylims: List or Tuple. :math:`y` interval. 

511 ''' 

512 error_handling.lims(xlims) 

513 error_handling.lims(ylims) 

514 in1 = (self.lat.coor['x'][self.hop['i']] >= xlims[0]) & \ 

515 (self.lat.coor['y'][self.hop['i']] >= ylims[0]) & \ 

516 (self.lat.coor['x'][self.hop['j']] >= xlims[0]) & \ 

517 (self.lat.coor['y'][self.hop['j']] >= ylims[0]) 

518 in2 = (self.lat.coor['x'][self.hop['i']] <= xlims[1]) & \ 

519 (self.lat.coor['y'][self.hop['i']] <= ylims[1]) & \ 

520 (self.lat.coor['x'][self.hop['j']] <= xlims[1]) & \ 

521 (self.lat.coor['y'][self.hop['j']] <= ylims[1]) 

522 return in1 * in2 

523 

524 def find_ellipse(self, rx: float, ry: float, x0: float, y0: float) -> NDArray: 

525 ''' 

526 Private method. 

527 Find hoppings within the ellipse. 

528 

529 :param rx: Positive Float. Radius along :math:`x`.  

530 :param ry: Positive Float. Radius along :math:`y`. 

531 :param x0: Float. Defalut value 0. :math:`x` center.  

532 :param y0: Float. Defalut value 0. :math:`x` center. 

533 ''' 

534 in1 = (self.lat.coor['x'][self.hop['i']] - x0) ** 2 / rx ** 2 + \ 

535 (self.lat.coor['y'][self.hop['i']] - y0) ** 2 / ry ** 2 <= 1. 

536 in2 = (self.lat.coor['x'][self.hop['j']] - x0) ** 2 / rx ** 2 + \ 

537 (self.lat.coor['y'][self.hop['j']] - y0) ** 2 / ry ** 2 <= 1. 

538 return in1 * in2 

539 

540 def change_hopping_square( 

541 self, list_hop: list[dict], xlims: tuple[float, float], ylims: tuple[float, float] = [-1., 1.], 

542 ) -> None: 

543 ''' 

544 Change hopping values. 

545 

546 :param list_hop: List of Dictionary (see set_hopping definition). 

547 :param xlims: List or Tuple. :math:`x` interval. 

548 :param ylims: List or Tuple. :math:`y` interval. 

549 ''' 

550 error_handling.empty_hop(self.hop) 

551 error_handling.set_hopping(list_hop, self.nmax) 

552 ind = self.find_square(xlims, ylims) 

553 self.set_new_hopping(list_hop, ind) 

554 

555 def change_hopping_ellipse( 

556 self, list_hop: list[dict], rx: float, ry: float, x0: float = 0., y0: float = 0., 

557 ) -> None: 

558 ''' 

559 Change hopping values. 

560 

561 :param list_hop: List of Dictionary (see set_hopping definition). 

562 :param rx: Positive Float. Radius along :math:`x`.  

563 :param ry: Positive Float. Radius along :math:`y`. 

564 :param x0: Float. Default value 0. :math:`x` center.  

565 :param y0: Float. Default value 0. :math:`y` center. 

566 ''' 

567 error_handling.empty_hop(self.hop) 

568 error_handling.set_hopping(list_hop, self.nmax) 

569 error_handling.positive_real(rx, 'rx') 

570 error_handling.positive_real(ry, 'rx') 

571 error_handling.real_number(x0, 'x0') 

572 error_handling.real_number(y0, 'y0') 

573 ind = self.find_ellipse(rx, ry, x0, y0) 

574 self.set_new_hopping(list_hop, ind) 

575 

576 def get_coor_hop(self) -> None: 

577 ''' 

578 Get the site coordinates in hopping space 

579 only considering the nearest neighbours hoppings. 

580 ''' 

581 error_handling.empty_hop(self.hop) 

582 visited = np.zeros(self.lat.sites, 'u2') 

583 self.coor_hop = np.zeros(self.lat.sites, dtype=COOR_DTYPE) 

584 self.coor_hop['tag'] = self.lat.coor['tag'] 

585 hop = self.hop[self.hop['n'] == 1] 

586 hop_down = np.copy(hop) 

587 hop_down['i'] = hop['j'] 

588 hop_down['j'] = hop[ 'i'] 

589 hop_down['ang'] = -180 + hop['ang'] 

590 hop = np.concatenate([hop, hop_down]) 

591 i_visit = np.min(hop['i']) 

592 while True: 

593 hs = hop[hop['i'] == i_visit] 

594 for h in hs: 

595 if visited[h['j']] == 2: 

596 continue 

597 self.coor_hop['x'][h['j']] = self.coor_hop['x'][i_visit] + \ 

598 h['t'].real*cos(PI / 180 * h['ang']) 

599 self.coor_hop['y'][h['j']] = self.coor_hop['y'][i_visit] + \ 

600 h['t'].real*sin(PI / 180 * h['ang']) 

601 visited[h['j']] = 1 

602 visited[i_visit] = 2 

603 explored = np.argwhere(visited == 1) 

604 if not explored.any(): 

605 break 

606 i_visit = explored[0, 0] 

607 

608 def get_ham(self) -> None: 

609 ''' 

610 Get the Tight-Binding Hamiltonian using sys.hop. 

611 ''' 

612 error_handling.empty_hop(self.hop) 

613 error_handling.hop_sites(self.hop, self.lat.sites) 

614 if np.all(self.hop['ang'] >= 0) or np.all(self.hop['ang'] < 0): 

615 self.ham = sparse.csr_matrix((self.hop['t'], (self.hop['i'], self.hop['j'])), 

616 shape=(self.lat.sites, self.lat.sites)) \ 

617 + sparse.csr_matrix((self.hop['t'].conj(), (self.hop['j'], self.hop['i'])), 

618 shape=(self.lat.sites, self.lat.sites)) 

619 else: 

620 self.ham = sparse.csr_matrix((self.hop['t'], (self.hop['i'], self.hop['j'])), 

621 shape=(self.lat.sites, self.lat.sites)) 

622 if self.onsite.size == self.lat.sites: 

623 self.ham += sparse.diags(self.onsite, 0) 

624 

625 def get_eig(self, eigenvec: bool = False, left: bool = False) -> None: 

626 ''' 

627 Get the eigenergies, eigenvectors and polarisation. 

628 

629 :param eigenvec: Boolean. Default value False.  

630 If True, get the eigenvectors. 

631 :param left: Boolean. Default value False.  

632 If True, get the left eigenvectors too.  

633 Relevant for non-Hermitian matrices. 

634 ''' 

635 error_handling.empty_ham(self.ham) 

636 error_handling.boolean(eigenvec, 'eigenvec') 

637 error_handling.boolean(left, 'left') 

638 if eigenvec: 

639 if (self.ham.conj().T != self.ham).nnz: 

640 if not left: 

641 self.en, self.rn = LA.eig(self.ham.toarray()) 

642 else: 

643 self.en, self.rn, self.ln = LA.eig(self.ham.toarray(), left=left) 

644 ind = np.argsort(self.en.real) 

645 self.en = self.en[ind] 

646 self.rn = self.rn[:, ind] 

647 if self.ln.size: 

648 self.ln = self.ln[:, ind] 

649 else: 

650 self.en, self.rn = LA.eigh(self.ham.toarray())