Coverage for tbkit/moire.py: 100%
200 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +0000
1r"""
2Supercells, band unfolding and moiré bilayers.
4* :func:`supercell` repeats a periodic model (:class:`tbkit.kspace.KSpace`)
5 in a larger cell spanned by integer combinations of its primitive vectors,
6 :math:`\mathbf{A}_i = \sum_j M_{ij}\mathbf{a}_j`: the natural home of
7 disorder, defects and superstructures that the primitive cell cannot hold.
8* :func:`unfold` and :func:`spectral_function` project the supercell
9 eigenstates back onto the primitive Brillouin zone (Ku, Berlijn and Lee,
10 Phys. Rev. Lett. 104, 216401 (2010); Popescu and Zunger, Phys. Rev. B
11 85, 085201 (2012)): the weight of state :math:`|\mathbf{K}J\rangle` at the
12 primitive :math:`\mathbf{k}` (folding onto :math:`\mathbf{K}`) is
14 .. math::
16 W_J(\mathbf{k}) = \frac{1}{N}\sum_o\Big|\sum_{m=1}^{N}
17 e^{-i\mathbf{k}\cdot\mathbf{r}_m}\,\psi_J(m, o)\Big|^2\, ,
19 :math:`\mathbf{r}_m` the origins of the :math:`N` primitive cells in the
20 supercell and :math:`o` the primitive orbitals, so that
21 :math:`A(\mathbf{k},\omega) = \sum_J W_J(\mathbf{k})\,\delta(\omega - E_J)`.
22 A pristine supercell unfolds exactly onto the primitive bands
23 (:math:`W = 1` on them, 0 elsewhere); disorder smears them.
24* :func:`twisted_bilayer` builds the commensurate supercells of two
25 honeycomb (twisted bilayer graphene) or square layers twisted by
26 :math:`\theta`, with every hopping a function of the distance
27 (:func:`pz_hopping`, the Slater-Koster :math:`p_z` form of Moon and
28 Koshino, Phys. Rev. B 85, 195458 (2012), by default).
29 :func:`magic_angle_parameter` gives the Bistritzer-MacDonald
30 :math:`\alpha = w/\hbar v k_\theta` and Fermi velocity
31 :math:`v^*/v = (1-3\alpha^2)/(1+6\alpha^2)` of the same model.
33Every Bloch matrix comes from :meth:`KSpace._bloch_ham <tbkit.kspace.KSpace>`
34(the supercells are ordinary **KSpace** models: *k_path*, *mesh_bands*,
35*finite_ham*, the topology tools all apply).
36"""
37from __future__ import annotations
39from typing import Callable
41import numpy as np
42from numpy.typing import ArrayLike, NDArray
43from scipy.integrate import quad
44from scipy.spatial import cKDTree
45from scipy.special import j0
47import tbkit.error_handling as error_handling
48from tbkit.kspace import KSpace
49from tbkit.lattice import Lattice
52PI = np.pi
55class SupercellKSpace(KSpace):
56 r'''
57 A **KSpace** made of :math:`N` copies of a primitive model (see
58 *supercell*). Orbital *o* of primitive cell *m* is orbital
59 ``m*prim.norb + o``.
61 Attributes
62 ----------
63 prim : **KSpace** instance
64 The primitive model.
65 matrix : Integer ndarray, shape (dim, dim)
66 Supercell matrix *M*.
67 cells : Real ndarray, shape (N, space_dim)
68 Origins
69 :math:`\mathbf{r}_m` of the primitive cells in the supercell.
70 n_cells : Integer
71 :math:`N = |\det M|`.
72 '''
74 def __init__(self, lat: Lattice, spin: bool, prim: KSpace, matrix: NDArray, cells: NDArray) -> None:
75 super().__init__(lat, spin)
76 self.prim = prim
77 self.matrix = matrix
78 self.cells = cells
79 self.n_cells = len(cells)
82class MoireKSpace(KSpace):
83 r'''
84 A **KSpace** of a twisted bilayer (see *twisted_bilayer*): the orbitals
85 of layer 0 (tags 'a', 'b', rotated by :math:`-\theta/2`) then those of
86 layer 1 (tags 'c', 'd', rotated by :math:`+\theta/2`, a height *d*
87 above).
89 Attributes
90 ----------
91 theta : Real number
92 Twist angle, in radians.
93 layer : Integer ndarray
94 Layer (0 or 1) of every orbital.
95 z : Real ndarray
96 Height of every orbital.
97 n_layer_cells : Integer
98 Primitive cells per layer in the supercell.
99 '''
101 def __init__(self, lat: Lattice, theta: float, layer: NDArray, z: NDArray, n_layer_cells: int) -> None:
102 super().__init__(lat)
103 self.theta = theta
104 self.layer = layer
105 self.z = z
106 self.n_layer_cells = n_layer_cells
109def _inner_cells(matrix: NDArray) -> NDArray:
110 '''
111 Private function. Integer lattice points n (rows) inside the supercell
112 spanned by the rows of *matrix*: n M^-1 in [0, 1)^dim, sorted.
113 '''
114 dim = len(matrix)
115 corners = np.array(np.meshgrid(*[[0, 1]] * dim, indexing='ij')).reshape(dim, -1).T @ matrix
116 ranges = [np.arange(corners[:, d].min(), corners[:, d].max() + 1) for d in range(dim)]
117 pts = np.array(np.meshgrid(*ranges, indexing='ij')).reshape(dim, -1).T
118 frac = pts @ np.linalg.inv(matrix)
119 eps = 1e-9
120 keep = np.all((frac > -eps) & (frac < 1 - eps), axis=1)
121 return pts[keep]
124def supercell(ks: KSpace, matrix: ArrayLike) -> SupercellKSpace:
125 r'''
126 Get the supercell of a periodic model spanned by
127 :math:`\mathbf{A}_i = \sum_j M_{ij}\mathbf{a}_j`, a **KSpace** whose unit
128 cell holds the :math:`N = |\det M|` primitive cells inside it. The bands
129 of the supercell are those of the primitive model folded into the
130 smaller Brillouin zone; *unfold* undoes the folding. Break the
131 translation symmetry afterwards by changing the supercell's
132 ``onsite`` array (one entry per orbital), e.g. with random energies.
134 Parameters
135 ----------
136 ks : **KSpace** instance (any dimension; spinful or not)
137 matrix
138 Integer array, shape (dim, dim), nonzero determinant
139 (an integer *n* in 1D is also accepted).
141 Returns
142 -------
143 sc
144 :class:`SupercellKSpace` instance, ``N*ks.norb`` orbitals:
145 orbital *o* of cell *m* (origin ``sc.cells[m]``) is ``m*ks.norb + o``.
147 Examples
148 --------
149 Example usage::
151 # a 3 x 3 supercell of graphene, 18 orbitals
152 sc = supercell(gra, [[3, 0], [0, 3]])
153 '''
154 error_handling.kspace(ks, KSpace)
155 if isinstance(matrix, int):
156 matrix = [[matrix]]
157 matrix = np.asarray(matrix)
158 error_handling.supercell_matrix(matrix, ks.dim)
159 matrix = matrix.astype(int)
160 a = np.array(ks.lat.prim_vec, dtype='f8')
161 cells = _inner_cells(matrix)
162 origins = cells @ a
163 cell = []
164 for r in origins:
165 for dic in ks.lat.unit_cell:
166 cell.append({'tag': dic['tag'], 'r0': tuple(float(c) for c in np.array(dic['r0']) + r)})
167 big = matrix @ a
168 lat = Lattice(unit_cell=cell, prim_vec=[tuple(float(c) for c in v) for v in big])
169 sc = SupercellKSpace(lat, ks.spin, ks, matrix, origins)
170 n_cells = len(cells)
171 lookup = {tuple(c): m for m, c in enumerate(cells)}
172 inv = np.linalg.inv(matrix)
173 rec = np.array(ks.rec_vec, dtype='f8')
175 def fold(hops):
176 out = []
177 for i, j, R_cart, t in hops:
178 R = np.rint(rec @ R_cart / (2 * PI)).astype(int)
179 target = cells + R[None, :]
180 shift = np.floor(target @ inv + 1e-9).astype(int)
181 inner = target - shift @ matrix
182 for m in range(n_cells):
183 out.append((m * ks.norb + i, lookup[tuple(inner[m])] * ks.norb + j,
184 shift[m] @ big, t))
185 return out
187 sc._hop = fold(ks._hop)
188 sc._overlap_hop = fold(ks._overlap_hop)
189 sc._nonreciprocal = ks._nonreciprocal
190 sc.onsite = np.tile(ks.onsite, n_cells)
191 sc._onsite_offdiag = np.kron(np.eye(n_cells), ks._onsite_offdiag)
192 return sc
195def unfold(sc: SupercellKSpace, ks: ArrayLike) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
196 r'''
197 Unfold the bands of a supercell onto k-points of the primitive
198 Brillouin zone (see the module docstring): diagonalize the supercell
199 Hamiltonian at each primitive :math:`\mathbf{k}` (which folds onto
200 itself modulo the supercell reciprocal lattice), and weigh each
201 eigenstate by its overlap with the primitive Bloch states at
202 :math:`\mathbf{k}`. The weights at one :math:`\mathbf{k}` add up to the
203 number of primitive orbitals; for a pristine supercell they are 1 on
204 the primitive bands :math:`E_n(\mathbf{k})` and 0 on the others.
206 Parameters
207 ----------
208 sc : :class:`SupercellKSpace` instance (Hermitian, without overlap)
209 ks : Real array, shape (nk, dim)
210 Primitive k-points, in the
211 coordinates of ``sc.prim.get_ham``.
213 Returns
214 -------
215 en : Real ndarray, shape (nk, sc.norb)
216 Supercell energies at each k.
217 weights : Real ndarray, shape (nk, sc.norb)
218 Spectral weights :math:`W_J(\mathbf{k})`.
219 '''
220 error_handling.supercell_model(sc, SupercellKSpace)
221 error_handling.hermitian_kspace(sc.is_hermitian(), sc._overlap_hop)
222 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))
223 error_handling.ks(ks, sc.dim)
224 k_cart = ks @ sc.prim.k_basis.T
225 en, vec = np.linalg.eigh(sc._bloch_ham(k_cart))
226 norb = sc.prim.norb
227 psi = vec.reshape(len(ks), sc.n_cells, norb, sc.norb)
228 phase = np.exp(-1j * k_cart @ sc.cells.T) # (nk, N)
229 amp = np.einsum('km,kmos->kos', phase, psi)
230 weights = (np.abs(amp) ** 2).sum(axis=1) / sc.n_cells
231 return en, weights
234def spectral_function(
235 sc: SupercellKSpace, ks: ArrayLike, energies: ArrayLike, broadening: float = 0.05,
236 kernel: str = 'lorentzian',
237) -> NDArray[np.float64]:
238 r'''
239 Get the unfolded spectral function of a supercell,
241 .. math::
243 A(\mathbf{k}, \omega) = \sum_J W_J(\mathbf{k})\, g(\omega - E_J)\, ,
245 with the weights of *unfold* and a normalized kernel *g* (Lorentzian or
246 Gaussian of width *broadening*), so that :math:`\int A\,d\omega` is
247 the number of primitive orbitals. For a disordered supercell, average
248 it over disorder configurations.
250 Parameters
251 ----------
252 sc : :class:`SupercellKSpace` instance
253 ks : Real array, shape (nk, dim)
254 Primitive k-points.
255 energies : Real array
256 Energies :math:`\omega`.
257 broadening : Positive real, default 0.05
258 Kernel width.
259 kernel : 'lorentzian' or 'gaussian', default 'lorentzian'
261 Returns
262 -------
263 A : Real ndarray, shape (nk, len(energies))
264 '''
265 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
266 error_handling.ndarray_empty(energies, 'energies')
267 error_handling.positive_real(broadening, 'broadening')
268 error_handling.dos_kernel(kernel)
269 en, weights = unfold(sc, ks)
270 diff = energies[None, :, None] - en[:, None, :]
271 if kernel == 'gaussian':
272 g = np.exp(-diff ** 2 / (2 * broadening ** 2)) / (broadening * np.sqrt(2 * PI))
273 else:
274 g = (broadening / PI) / (diff ** 2 + broadening ** 2)
275 return np.einsum('kej,kj->ke', g, weights)
278def pz_hopping(
279 d: ArrayLike, t_pi: float = -2.7, t_sigma: float = 0.48, a0: float = 0.142,
280 d0: float = 0.335, delta: float = 0.0453,
281) -> NDArray[np.float64]:
282 r'''
283 Get the hopping between two :math:`p_z` orbitals a vector
284 :math:`\mathbf{d}` apart, in the Slater-Koster form of Moon and Koshino
285 (Phys. Rev. B 85, 195458 (2012)), used for twisted bilayer graphene:
287 .. math::
289 t(\mathbf{d}) = V_{pp\pi}(d)\Big[1-\Big(\frac{d_z}{d}\Big)^2\Big]
290 + V_{pp\sigma}(d)\Big(\frac{d_z}{d}\Big)^2\, ,\quad
291 V_{pp\pi} = t_\pi e^{-(d-a_0)/\delta}\, ,\quad
292 V_{pp\sigma} = t_\sigma e^{-(d-d_0)/\delta}\, .
294 In-plane it is :math:`t_\pi` between nearest neighbours and decays with
295 distance; between the layers of AA-stacked graphene it is
296 :math:`t_\sigma`. Defaults: graphene, in eV and nm (:math:`\delta =
297 0.184\,a`, :math:`a = 0.246` nm).
299 Parameters
300 ----------
301 d : Real array, shape (n, 3)
302 Bond vectors.
303 t_pi : Real number, default -2.7
304 t_sigma : Real number, default 0.48
305 a0 : Positive real, default 0.142
306 In-plane nearest-neighbour distance.
307 d0 : Positive real, default 0.335
308 Interlayer distance.
309 delta : Positive real, default 0.0453
310 Decay length.
312 Returns
313 -------
314 t : Real ndarray, shape (n,)
315 '''
316 d = np.atleast_2d(np.asarray(d, dtype='f8'))
317 error_handling.bond_vectors(d)
318 error_handling.real_number(t_pi, 't_pi')
319 error_handling.real_number(t_sigma, 't_sigma')
320 for val, name in ((a0, 'a0'), (d0, 'd0'), (delta, 'delta')):
321 error_handling.positive_real(val, name)
322 r = np.linalg.norm(d, axis=1)
323 c2 = (d[:, 2] / r) ** 2
324 return t_pi * np.exp(-(r - a0) / delta) * (1 - c2) + t_sigma * np.exp(-(r - d0) / delta) * c2
327def _twist_matrices(m: int, lattice: str) -> tuple[NDArray, NDArray]:
328 '''
329 Private function. Supercell matrices of the two layers: the rows are the
330 moire vectors in units of each layer's primitive vectors.
331 '''
332 if lattice == 'honeycomb': # R60 a1 = a2, R60 a2 = a2 - a1
333 return (np.array([[m, m + 1], [-(m + 1), 2 * m + 1]]),
334 np.array([[m + 1, m], [-m, 2 * m + 1]]))
335 # square: R90 a1 = a2, R90 a2 = -a1
336 return (np.array([[m, m + 1], [-(m + 1), m]]), np.array([[m + 1, m], [-m, m + 1]]))
339def _layer_lattice(lattice: str, a: float) -> tuple[NDArray, NDArray, list[str]]:
340 '''
341 Private function. Primitive vectors, orbital positions and tags of one layer.
342 '''
343 if lattice == 'honeycomb':
344 dx = 0.5 * np.sqrt(3) * a
345 return np.array([[2 * dx, 0.], [dx, 1.5 * a]]), np.array([[0., 0.], [dx, 0.5 * a]]), ['a', 'b']
346 return np.array([[a, 0.], [0., a]]), np.array([[0., 0.]]), ['a']
349def commensurate_angle(m: int, lattice: str = 'honeycomb') -> float:
350 r'''
351 Get the commensurate twist angle :math:`\theta_m` of two layers: the
352 rotation that maps the lattice vector :math:`(m+1)\mathbf{a}_1 + m\mathbf{a}_2`
353 onto :math:`m\mathbf{a}_1 + (m+1)\mathbf{a}_2`,
355 .. math::
357 \cos\theta_m = \frac{3m^2+3m+1/2}{3m^2+3m+1}\ \text{(honeycomb)},\qquad
358 \cos\theta_m = \frac{2m(m+1)}{2m^2+2m+1}\ \text{(square)}
360 (Lopes dos Santos, Peres and Castro Neto, Phys. Rev. B 86, 155449
361 (2012)). The moiré cell holds :math:`3m^2+3m+1` (honeycomb) or
362 :math:`2m^2+2m+1` (square) primitive cells per layer; :math:`m = 31`
363 gives graphene's first magic angle, :math:`\theta \approx 1.05°`.
365 Parameters
366 ----------
367 m : Positive integer
368 lattice : 'honeycomb' or 'square', default 'honeycomb'
370 Returns
371 -------
372 theta : Real number, in radians
373 '''
374 error_handling.positive_int(m, 'm')
375 error_handling.twist_lattice(lattice)
376 a1, _ = _twist_matrices(m, lattice)
377 prim, _, _ = _layer_lattice(lattice, 1.)
378 v, w = a1[0] @ prim, a1[0][::-1] @ prim
379 return float(abs(np.arctan2(v[0] * w[1] - v[1] * w[0], v @ w)))
382def twisted_bilayer(
383 m: int, lattice: str = 'honeycomb', a: float = 0.142, d: float = 0.335,
384 hopping: Callable[[NDArray], NDArray] | None = None, cutoff: float | None = None,
385 interlayer_scale: float = 1.,
386) -> MoireKSpace:
387 r'''
388 Get the commensurate twisted bilayer of two honeycomb or square layers,
389 twisted by :math:`\theta_m` (see *commensurate_angle*) about a common
390 site (AA stacking at the origin): a **KSpace** with the moiré cell as
391 unit cell, and every hopping, within a layer and between the layers,
392 given by the function *hopping* of the bond vector
393 :math:`\mathbf{d} = (d_x, d_y, d_z)`, for all bonds shorter than
394 *cutoff* (periodic images included).
396 Near the first magic angle (:math:`\theta \approx 1.05°` for graphene,
397 :math:`m = 31`, 11908 orbitals) the two bands at charge neutrality
398 become almost flat (Bistritzer and MacDonald, PNAS 108, 12233 (2011);
399 Cao et al., Nature 556, 43 (2018)). The physics depends on the ratio
400 :math:`\alpha = w/\hbar v k_\theta` only (see *magic_angle_parameter*),
401 so *interlayer_scale* (which multiplies every interlayer hopping)
402 moves the magic angle to larger angles and smaller cells.
404 Parameters
405 ----------
406 m : Positive integer
407 See *commensurate_angle*.
408 lattice : 'honeycomb' or 'square', default 'honeycomb'
409 a : Positive real
410 Default value 0.142 (graphene, nm). Nearest-neighbour distance.
411 d : Positive real
412 Default value 0.335 (nm). Interlayer distance.
413 hopping : Callable
414 Default value None (*pz_hopping*). Maps bond
415 vectors, shape (n, 3), to hoppings, shape (n,).
416 cutoff : Positive real
417 Default value None (3.8 *a*, between the
418 honeycomb neighbour shells at :math:`\sqrt{13}a` and :math:`4a`, so
419 that no shell is cut by rounding). Longest bond.
420 interlayer_scale : Positive real, default 1
421 Factor on the
422 interlayer hoppings.
424 Returns
425 -------
426 tbl : :class:`MoireKSpace` instance
427 '''
428 error_handling.positive_int(m, 'm')
429 error_handling.twist_lattice(lattice)
430 error_handling.positive_real(a, 'a')
431 error_handling.positive_real(d, 'd')
432 if hopping is None:
433 hopping = pz_hopping
434 error_handling.is_callable(hopping, 'hopping')
435 if cutoff is None:
436 cutoff = 3.8 * a
437 error_handling.positive_real(cutoff, 'cutoff')
438 error_handling.positive_real(interlayer_scale, 'interlayer_scale')
439 theta = commensurate_angle(m, lattice)
440 prim, tau, tags = _layer_lattice(lattice, a)
441 mats = _twist_matrices(m, lattice)
442 pos_parts: list = []
443 z_parts: list = []
444 layer_parts: list = []
445 cell_tags: list[str] = []
446 for l, (mat, sign) in enumerate(zip(mats, (-1., 1.))):
447 c, s = np.cos(sign * theta / 2), np.sin(sign * theta / 2)
448 rot = np.array([[c, -s], [s, c]])
449 origins = _inner_cells(mat) @ prim
450 p = (origins[:, None, :] + tau[None]).reshape(-1, 2) @ rot.T
451 pos_parts.append(p)
452 z_parts.append(np.full(len(p), l * d))
453 layer_parts.append(np.full(len(p), l))
454 cell_tags += [chr(ord(t) + 2 * l) for t in tags] * len(origins)
455 big = mat @ prim @ rot.T
456 if l == 0:
457 moire = big
458 pos, z, layer = np.concatenate(pos_parts), np.concatenate(z_parts), np.concatenate(layer_parts)
459 error_handling.commensurate(moire, big)
460 cell = [{'tag': t, 'r0': (float(x), float(y))} for t, (x, y) in zip(cell_tags, pos)]
461 tbl = MoireKSpace(Lattice(unit_cell=cell, prim_vec=[tuple(float(c) for c in v) for v in moire]),
462 theta, layer, z, len(pos) // (2 * len(tau)))
463 # bonds to the periodic images within the cutoff
464 area = abs(np.linalg.det(moire))
465 n_img = [int(np.ceil(cutoff / (area / np.linalg.norm(moire[1 - i])))) + 1 for i in range(2)]
466 tree = cKDTree(pos)
467 hops = []
468 for l1 in range(-n_img[0], n_img[0] + 1):
469 for l2 in range(-n_img[1], n_img[1] + 1):
470 if (l1, l2) < (0, 0):
471 continue # the conjugates of these are added below
472 R = l1 * moire[0] + l2 * moire[1]
473 pairs = tree.sparse_distance_matrix(cKDTree(pos + R), cutoff, output_type='ndarray')
474 i, j = pairs['i'], pairs['j']
475 keep = i < j if (l1, l2) == (0, 0) else np.ones(len(i), bool)
476 i, j = i[keep], j[keep]
477 vec = np.column_stack([pos[j] + R - pos[i], z[j] - z[i]])
478 short = np.linalg.norm(vec, axis=1) <= cutoff
479 i, j, vec = i[short], j[short], vec[short]
480 if len(vec) == 0:
481 continue
482 t = np.asarray(hopping(vec))
483 error_handling.hopping_values(t, len(vec))
484 t = np.where(layer[i] != layer[j], interlayer_scale * t, t)
485 for ii, jj, tt in zip(i, j, t):
486 hops.append((int(ii), int(jj), R, tt))
487 hops.append((int(jj), int(ii), -R, np.conj(tt)))
488 tbl._hop = hops
489 return tbl
492def magic_angle_parameter(
493 theta: float, a: float = 0.142, d: float = 0.335, t: float = -2.7,
494 interlayer: Callable[[NDArray], NDArray] | None = None, interlayer_scale: float = 1.,
495) -> tuple[float, float, float]:
496 r'''
497 Get the Bistritzer-MacDonald parameters of twisted bilayer graphene
498 for a tight-binding model with nearest-neighbour hopping *t* and the
499 interlayer hopping *interlayer* (a function of the bond vector): the
500 interlayer tunnelling :math:`w = \frac{1}{\Omega}\int d^2r\,
501 t_\perp(\mathbf{r}, d)\,e^{-i\mathbf{K}\cdot\mathbf{r}}`
502 (:math:`\Omega = \sqrt3a_{\mathrm{lat}}^2/2` the cell area,
503 :math:`|\mathbf{K}| = 4\pi/3a_{\mathrm{lat}}`, :math:`a_{\mathrm{lat}} = \sqrt3 a`),
504 the Dirac velocity :math:`\hbar v = \frac32|t|a`, the ratio
505 :math:`\alpha = w/\hbar vk_\theta` (:math:`k_\theta = 2|\mathbf{K}|\sin\frac\theta2`),
506 and the renormalized velocity at the Dirac points of the moiré bands
507 (Bistritzer and MacDonald, PNAS 108, 12233 (2011)),
509 .. math::
511 \frac{v^*}{v} = \frac{1-3\alpha^2}{1+6\alpha^2}\, ,
513 which vanishes at the first magic angle, :math:`\alpha = 1/\sqrt3`.
515 Parameters
516 ----------
517 theta : Positive real
518 Twist angle, in radians.
519 a : Positive real, default 0.142
520 Nearest-neighbour distance.
521 d : Positive real, default 0.335
522 Interlayer distance.
523 t : Real number, default -2.7
524 Nearest-neighbour hopping.
525 interlayer : Callable
526 Default value None (*pz_hopping*). Interlayer hopping.
527 interlayer_scale : Positive real, default 1
528 See *twisted_bilayer*.
530 Returns
531 -------
532 w : Real number
533 Interlayer tunnelling.
534 alpha : Real number
535 v_ratio : Real number
536 :math:`v^*/v`.
537 '''
538 error_handling.positive_real(theta, 'theta')
539 error_handling.positive_real(a, 'a')
540 error_handling.positive_real(d, 'd')
541 error_handling.real_number(t, 't')
542 error_handling.positive_real(interlayer_scale, 'interlayer_scale')
543 if interlayer is None:
544 interlayer = pz_hopping
545 error_handling.is_callable(interlayer, 'interlayer')
546 a_lat = np.sqrt(3) * a
547 k_d = 4 * PI / (3 * a_lat)
548 omega = np.sqrt(3) / 2 * a_lat ** 2
550 def integrand(r):
551 return r * float(np.real(interlayer(np.array([[r, 0., d]]))[0])) * j0(k_d * r)
552 w = interlayer_scale * 2 * PI * quad(integrand, 0., 20 * d, limit=200)[0] / omega
553 k_theta = 2 * k_d * np.sin(theta / 2)
554 alpha = abs(w) / (1.5 * abs(t) * a * k_theta)
555 return float(abs(w)), float(alpha), float((1 - 3 * alpha ** 2) / (1 + 6 * alpha ** 2))