Coverage for tbkit/moire.py: 100%

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 :ivar prim: **KSpace** instance. The primitive model.

62 :ivar matrix: Integer ndarray, shape (dim, dim). Supercell matrix *M*.

63 :ivar cells: Real ndarray, shape (N, space_dim). Origins

64 :math:`\mathbf{r}_m` of the primitive cells in the supercell.

65 :ivar n_cells: Integer. :math:`N = |\det M|`.

66 '''

68 def __init__(self, lat: Lattice, spin: bool, prim: KSpace, matrix: NDArray, cells: NDArray) -> None:

69 super().__init__(lat, spin)

70 self.prim = prim

71 self.matrix = matrix

72 self.cells = cells

73 self.n_cells = len(cells)

76class MoireKSpace(KSpace):

77 r'''

78 A **KSpace** of a twisted bilayer (see *twisted_bilayer*): the orbitals

79 of layer 0 (tags 'a', 'b', rotated by :math:`-\theta/2`) then those of

80 layer 1 (tags 'c', 'd', rotated by :math:`+\theta/2`, a height *d*

81 above).

83 :ivar theta: Real number. Twist angle, in radians.

84 :ivar layer: Integer ndarray. Layer (0 or 1) of every orbital.

85 :ivar z: Real ndarray. Height of every orbital.

86 :ivar n_layer_cells: Integer. Primitive cells per layer in the supercell.

87 '''

89 def __init__(self, lat: Lattice, theta: float, layer: NDArray, z: NDArray, n_layer_cells: int) -> None:

90 super().__init__(lat)

91 self.theta = theta

92 self.layer = layer

93 self.z = z

94 self.n_layer_cells = n_layer_cells

97def _inner_cells(matrix: NDArray) -> NDArray:

98 '''

99 Private function. Integer lattice points n (rows) inside the supercell

100 spanned by the rows of *matrix*: n M^-1 in [0, 1)^dim, sorted.

101 '''

102 dim = len(matrix)

103 corners = np.array(np.meshgrid(*[[0, 1]] * dim, indexing='ij')).reshape(dim, -1).T @ matrix

104 ranges = [np.arange(corners[:, d].min(), corners[:, d].max() + 1) for d in range(dim)]

105 pts = np.array(np.meshgrid(*ranges, indexing='ij')).reshape(dim, -1).T

106 frac = pts @ np.linalg.inv(matrix)

107 eps = 1e-9

108 keep = np.all((frac > -eps) & (frac < 1 - eps), axis=1)

109 return pts[keep]

110

111

112def supercell(ks: KSpace, matrix: ArrayLike) -> SupercellKSpace:

113 r'''

114 Get the supercell of a periodic model spanned by

115 :math:`\mathbf{A}_i = \sum_j M_{ij}\mathbf{a}_j`, a **KSpace** whose unit

116 cell holds the :math:`N = |\det M|` primitive cells inside it. The bands

117 of the supercell are those of the primitive model folded into the

118 smaller Brillouin zone; *unfold* undoes the folding. Break the

119 translation symmetry afterwards by changing the supercell's

120 ``onsite`` array (one entry per orbital), e.g. with random energies.

121

122 :param ks: **KSpace** instance (any dimension; spinful or not).

123 :param matrix: Integer array, shape (dim, dim), nonzero determinant

124 (an integer *n* in 1D is also accepted).

125

126 :returns:

127 * **sc** -- :class:`SupercellKSpace` instance, ``N*ks.norb`` orbitals:

128 orbital *o* of cell *m* (origin ``sc.cells[m]``) is ``m*ks.norb + o``.

129

130 Example usage::

131

132 # a 3 x 3 supercell of graphene, 18 orbitals

133 sc = supercell(gra, [[3, 0], [0, 3]])

134 '''

135 error_handling.kspace(ks, KSpace)

136 if isinstance(matrix, int):

137 matrix = [[matrix]]

138 matrix = np.asarray(matrix)

139 error_handling.supercell_matrix(matrix, ks.dim)

140 matrix = matrix.astype(int)

141 a = np.array(ks.lat.prim_vec, dtype='f8')

142 cells = _inner_cells(matrix)

143 origins = cells @ a

144 cell = []

145 for r in origins:

146 for dic in ks.lat.unit_cell:

147 cell.append({'tag': dic['tag'], 'r0': tuple(float(c) for c in np.array(dic['r0']) + r)})

148 big = matrix @ a

149 lat = Lattice(unit_cell=cell, prim_vec=[tuple(float(c) for c in v) for v in big])

150 sc = SupercellKSpace(lat, ks.spin, ks, matrix, origins)

151 n_cells = len(cells)

152 lookup = {tuple(c): m for m, c in enumerate(cells)}

153 inv = np.linalg.inv(matrix)

154 rec = np.array(ks.rec_vec, dtype='f8')

155

156 def fold(hops):

157 out = []

158 for i, j, R_cart, t in hops:

159 R = np.rint(rec @ R_cart / (2 * PI)).astype(int)

160 target = cells + R[None, :]

161 shift = np.floor(target @ inv + 1e-9).astype(int)

162 inner = target - shift @ matrix

163 for m in range(n_cells):

164 out.append((m * ks.norb + i, lookup[tuple(inner[m])] * ks.norb + j,

165 shift[m] @ big, t))

166 return out

167

168 sc._hop = fold(ks._hop)

169 sc._overlap_hop = fold(ks._overlap_hop)

170 sc._nonreciprocal = ks._nonreciprocal

171 sc.onsite = np.tile(ks.onsite, n_cells)

172 sc._onsite_offdiag = np.kron(np.eye(n_cells), ks._onsite_offdiag)

173 return sc

174

175

176def unfold(sc: SupercellKSpace, ks: ArrayLike) -> tuple[NDArray[np.float64], NDArray[np.float64]]:

177 r'''

178 Unfold the bands of a supercell onto k-points of the primitive

179 Brillouin zone (see the module docstring): diagonalize the supercell

180 Hamiltonian at each primitive :math:`\mathbf{k}` (which folds onto

181 itself modulo the supercell reciprocal lattice), and weigh each

182 eigenstate by its overlap with the primitive Bloch states at

183 :math:`\mathbf{k}`. The weights at one :math:`\mathbf{k}` add up to the

184 number of primitive orbitals; for a pristine supercell they are 1 on

185 the primitive bands :math:`E_n(\mathbf{k})` and 0 on the others.

186

187 :param sc: :class:`SupercellKSpace` instance (Hermitian, without overlap).

188 :param ks: Real array, shape (nk, dim). Primitive k-points, in the

189 coordinates of ``sc.prim.get_ham``.

190

191 :returns:

192 * **en** -- Real ndarray, shape (nk, sc.norb). Supercell energies at each k.

193 * **weights** -- Real ndarray, shape (nk, sc.norb). Spectral weights :math:`W_J(\mathbf{k})`.

194 '''

195 error_handling.supercell_model(sc, SupercellKSpace)

196 error_handling.hermitian_kspace(sc.is_hermitian(), sc._overlap_hop)

197 ks = np.atleast_2d(np.asarray(ks, dtype='f8'))

198 error_handling.ks(ks, sc.dim)

199 k_cart = ks @ sc.prim.k_basis.T

200 en, vec = np.linalg.eigh(sc._bloch_ham(k_cart))

201 norb = sc.prim.norb

202 psi = vec.reshape(len(ks), sc.n_cells, norb, sc.norb)

203 phase = np.exp(-1j * k_cart @ sc.cells.T) # (nk, N)

204 amp = np.einsum('km,kmos->kos', phase, psi)

205 weights = (np.abs(amp) ** 2).sum(axis=1) / sc.n_cells

206 return en, weights

207

208

209def spectral_function(

210 sc: SupercellKSpace, ks: ArrayLike, energies: ArrayLike, broadening: float = 0.05,

211 kernel: str = 'lorentzian',

212) -> NDArray[np.float64]:

213 r'''

214 Get the unfolded spectral function of a supercell,

215

216 .. math::

217

218 A(\mathbf{k}, \omega) = \sum_J W_J(\mathbf{k})\, g(\omega - E_J)\, ,

219

220 with the weights of *unfold* and a normalized kernel *g* (Lorentzian or

221 Gaussian of width *broadening*), so that :math:`\int A\,d\omega` is

222 the number of primitive orbitals. For a disordered supercell, average

223 it over disorder configurations.

224

225 :param sc: :class:`SupercellKSpace` instance.

226 :param ks: Real array, shape (nk, dim). Primitive k-points.

227 :param energies: Real array. Energies :math:`\omega`.

228 :param broadening: Positive real. Default value 0.05. Kernel width.

229 :param kernel: 'lorentzian' or 'gaussian'. Default value 'lorentzian'.

230

231 :returns:

232 * **A** -- Real ndarray, shape (nk, len(energies)).

233 '''

234 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))

235 error_handling.ndarray_empty(energies, 'energies')

236 error_handling.positive_real(broadening, 'broadening')

237 error_handling.dos_kernel(kernel)

238 en, weights = unfold(sc, ks)

239 diff = energies[None, :, None] - en[:, None, :]

240 if kernel == 'gaussian':

241 g = np.exp(-diff ** 2 / (2 * broadening ** 2)) / (broadening * np.sqrt(2 * PI))

242 else:

243 g = (broadening / PI) / (diff ** 2 + broadening ** 2)

244 return np.einsum('kej,kj->ke', g, weights)

245

246

247def pz_hopping(

248 d: ArrayLike, t_pi: float = -2.7, t_sigma: float = 0.48, a0: float = 0.142,

249 d0: float = 0.335, delta: float = 0.0453,

250) -> NDArray[np.float64]:

251 r'''

252 Get the hopping between two :math:`p_z` orbitals a vector

253 :math:`\mathbf{d}` apart, in the Slater-Koster form of Moon and Koshino

254 (Phys. Rev. B 85, 195458 (2012)), used for twisted bilayer graphene:

255

256 .. math::

257

258 t(\mathbf{d}) = V_{pp\pi}(d)\Big[1-\Big(\frac{d_z}{d}\Big)^2\Big]

259 + V_{pp\sigma}(d)\Big(\frac{d_z}{d}\Big)^2\, ,\quad

260 V_{pp\pi} = t_\pi e^{-(d-a_0)/\delta}\, ,\quad

261 V_{pp\sigma} = t_\sigma e^{-(d-d_0)/\delta}\, .

262

263 In-plane it is :math:`t_\pi` between nearest neighbours and decays with

264 distance; between the layers of AA-stacked graphene it is

265 :math:`t_\sigma`. Defaults: graphene, in eV and nm (:math:`\delta =

266 0.184\,a`, :math:`a = 0.246` nm).

267

268 :param d: Real array, shape (n, 3). Bond vectors.

269 :param t_pi: Real number. Default value -2.7.

270 :param t_sigma: Real number. Default value 0.48.

271 :param a0: Positive real. Default value 0.142. In-plane nearest-neighbour distance.

272 :param d0: Positive real. Default value 0.335. Interlayer distance.

273 :param delta: Positive real. Default value 0.0453. Decay length.

274

275 :returns:

276 * **t** -- Real ndarray, shape (n,).

277 '''

278 d = np.atleast_2d(np.asarray(d, dtype='f8'))

279 error_handling.bond_vectors(d)

280 error_handling.real_number(t_pi, 't_pi')

281 error_handling.real_number(t_sigma, 't_sigma')

282 for val, name in ((a0, 'a0'), (d0, 'd0'), (delta, 'delta')):

283 error_handling.positive_real(val, name)

284 r = np.linalg.norm(d, axis=1)

285 c2 = (d[:, 2] / r) ** 2

286 return t_pi * np.exp(-(r - a0) / delta) * (1 - c2) + t_sigma * np.exp(-(r - d0) / delta) * c2

287

288

289def _twist_matrices(m: int, lattice: str) -> tuple[NDArray, NDArray]:

290 '''

291 Private function. Supercell matrices of the two layers: the rows are the

292 moire vectors in units of each layer's primitive vectors.

293 '''

294 if lattice == 'honeycomb': # R60 a1 = a2, R60 a2 = a2 - a1

295 return (np.array([[m, m + 1], [-(m + 1), 2 * m + 1]]),

296 np.array([[m + 1, m], [-m, 2 * m + 1]]))

297 # square: R90 a1 = a2, R90 a2 = -a1

298 return (np.array([[m, m + 1], [-(m + 1), m]]), np.array([[m + 1, m], [-m, m + 1]]))

299

300

301def _layer_lattice(lattice: str, a: float) -> tuple[NDArray, NDArray, list[str]]:

302 '''

303 Private function. Primitive vectors, orbital positions and tags of one layer.

304 '''

305 if lattice == 'honeycomb':

306 dx = 0.5 * np.sqrt(3) * a

307 return np.array([[2 * dx, 0.], [dx, 1.5 * a]]), np.array([[0., 0.], [dx, 0.5 * a]]), ['a', 'b']

308 return np.array([[a, 0.], [0., a]]), np.array([[0., 0.]]), ['a']

309

310

311def commensurate_angle(m: int, lattice: str = 'honeycomb') -> float:

312 r'''

313 Get the commensurate twist angle :math:`\theta_m` of two layers: the

314 rotation that maps the lattice vector :math:`(m+1)\mathbf{a}_1 + m\mathbf{a}_2`

315 onto :math:`m\mathbf{a}_1 + (m+1)\mathbf{a}_2`,

316

317 .. math::

318

319 \cos\theta_m = \frac{3m^2+3m+1/2}{3m^2+3m+1}\ \text{(honeycomb)},\qquad

320 \cos\theta_m = \frac{2m(m+1)}{2m^2+2m+1}\ \text{(square)}

321

322 (Lopes dos Santos, Peres and Castro Neto, Phys. Rev. B 86, 155449

323 (2012)). The moiré cell holds :math:`3m^2+3m+1` (honeycomb) or

324 :math:`2m^2+2m+1` (square) primitive cells per layer; :math:`m = 31`

325 gives graphene's first magic angle, :math:`\theta \approx 1.05°`.

326

327 :param m: Positive integer.

328 :param lattice: 'honeycomb' or 'square'. Default value 'honeycomb'.

329

330 :returns:

331 * **theta** -- Real number, in radians.

332 '''

333 error_handling.positive_int(m, 'm')

334 error_handling.twist_lattice(lattice)

335 a1, _ = _twist_matrices(m, lattice)

336 prim, _, _ = _layer_lattice(lattice, 1.)

337 v, w = a1[0] @ prim, a1[0][::-1] @ prim

338 return float(abs(np.arctan2(v[0] * w[1] - v[1] * w[0], v @ w)))

339

340

341def twisted_bilayer(

342 m: int, lattice: str = 'honeycomb', a: float = 0.142, d: float = 0.335,

343 hopping: Callable[[NDArray], NDArray] | None = None, cutoff: float | None = None,

344 interlayer_scale: float = 1.,

345) -> MoireKSpace:

346 r'''

347 Get the commensurate twisted bilayer of two honeycomb or square layers,

348 twisted by :math:`\theta_m` (see *commensurate_angle*) about a common

349 site (AA stacking at the origin): a **KSpace** with the moiré cell as

350 unit cell, and every hopping, within a layer and between the layers,

351 given by the function *hopping* of the bond vector

352 :math:`\mathbf{d} = (d_x, d_y, d_z)`, for all bonds shorter than

353 *cutoff* (periodic images included).

354

355 Near the first magic angle (:math:`\theta \approx 1.05°` for graphene,

356 :math:`m = 31`, 11908 orbitals) the two bands at charge neutrality

357 become almost flat (Bistritzer and MacDonald, PNAS 108, 12233 (2011);

358 Cao et al., Nature 556, 43 (2018)). The physics depends on the ratio

359 :math:`\alpha = w/\hbar v k_\theta` only (see *magic_angle_parameter*),

360 so *interlayer_scale* (which multiplies every interlayer hopping)

361 moves the magic angle to larger angles and smaller cells.

362

363 :param m: Positive integer. See *commensurate_angle*.

364 :param lattice: 'honeycomb' or 'square'. Default value 'honeycomb'.

365 :param a: Positive real. Default value 0.142 (graphene, nm). Nearest-neighbour distance.

366 :param d: Positive real. Default value 0.335 (nm). Interlayer distance.

367 :param hopping: Callable. Default value None (*pz_hopping*). Maps bond

368 vectors, shape (n, 3), to hoppings, shape (n,).

369 :param cutoff: Positive real. Default value None (3.8 *a*, between the

370 honeycomb neighbour shells at :math:`\sqrt{13}a` and :math:`4a`, so

371 that no shell is cut by rounding). Longest bond.

372 :param interlayer_scale: Positive real. Default value 1. Factor on the

373 interlayer hoppings.

374

375 :returns:

376 * **tbl** -- :class:`MoireKSpace` instance.

377 '''

378 error_handling.positive_int(m, 'm')

379 error_handling.twist_lattice(lattice)

380 error_handling.positive_real(a, 'a')

381 error_handling.positive_real(d, 'd')

382 if hopping is None:

383 hopping = pz_hopping

384 error_handling.is_callable(hopping, 'hopping')

385 if cutoff is None:

386 cutoff = 3.8 * a

387 error_handling.positive_real(cutoff, 'cutoff')

388 error_handling.positive_real(interlayer_scale, 'interlayer_scale')

389 theta = commensurate_angle(m, lattice)

390 prim, tau, tags = _layer_lattice(lattice, a)

391 mats = _twist_matrices(m, lattice)

392 pos_parts: list = []

393 z_parts: list = []

394 layer_parts: list = []

395 cell_tags: list[str] = []

396 for l, (mat, sign) in enumerate(zip(mats, (-1., 1.))):

397 c, s = np.cos(sign * theta / 2), np.sin(sign * theta / 2)

398 rot = np.array([[c, -s], [s, c]])

399 origins = _inner_cells(mat) @ prim

400 p = (origins[:, None, :] + tau[None]).reshape(-1, 2) @ rot.T

401 pos_parts.append(p)

402 z_parts.append(np.full(len(p), l * d))

403 layer_parts.append(np.full(len(p), l))

404 cell_tags += [chr(ord(t) + 2 * l) for t in tags] * len(origins)

405 big = mat @ prim @ rot.T

406 if l == 0:

407 moire = big

408 pos, z, layer = np.concatenate(pos_parts), np.concatenate(z_parts), np.concatenate(layer_parts)

409 error_handling.commensurate(moire, big)

410 cell = [{'tag': t, 'r0': (float(x), float(y))} for t, (x, y) in zip(cell_tags, pos)]

411 tbl = MoireKSpace(Lattice(unit_cell=cell, prim_vec=[tuple(float(c) for c in v) for v in moire]),

412 theta, layer, z, len(pos) // (2 * len(tau)))

413 # bonds to the periodic images within the cutoff

414 area = abs(np.linalg.det(moire))

415 n_img = [int(np.ceil(cutoff / (area / np.linalg.norm(moire[1 - i])))) + 1 for i in range(2)]

416 tree = cKDTree(pos)

417 hops = []

418 for l1 in range(-n_img[0], n_img[0] + 1):

419 for l2 in range(-n_img[1], n_img[1] + 1):

420 if (l1, l2) < (0, 0):

421 continue # the conjugates of these are added below

422 R = l1 * moire[0] + l2 * moire[1]

423 pairs = tree.sparse_distance_matrix(cKDTree(pos + R), cutoff, output_type='ndarray')

424 i, j = pairs['i'], pairs['j']

425 keep = i < j if (l1, l2) == (0, 0) else np.ones(len(i), bool)

426 i, j = i[keep], j[keep]

427 vec = np.column_stack([pos[j] + R - pos[i], z[j] - z[i]])

428 short = np.linalg.norm(vec, axis=1) <= cutoff

429 i, j, vec = i[short], j[short], vec[short]

430 if len(vec) == 0:

431 continue

432 t = np.asarray(hopping(vec))

433 error_handling.hopping_values(t, len(vec))

434 t = np.where(layer[i] != layer[j], interlayer_scale * t, t)

435 for ii, jj, tt in zip(i, j, t):

436 hops.append((int(ii), int(jj), R, tt))

437 hops.append((int(jj), int(ii), -R, np.conj(tt)))

438 tbl._hop = hops

439 return tbl

440

441

442def magic_angle_parameter(

443 theta: float, a: float = 0.142, d: float = 0.335, t: float = -2.7,

444 interlayer: Callable[[NDArray], NDArray] | None = None, interlayer_scale: float = 1.,

445) -> tuple[float, float, float]:

446 r'''

447 Get the Bistritzer-MacDonald parameters of twisted bilayer graphene

448 for a tight-binding model with nearest-neighbour hopping *t* and the

449 interlayer hopping *interlayer* (a function of the bond vector): the

450 interlayer tunnelling :math:`w = \frac{1}{\Omega}\int d^2r\,

451 t_\perp(\mathbf{r}, d)\,e^{-i\mathbf{K}\cdot\mathbf{r}}`

452 (:math:`\Omega = \sqrt3a_{\mathrm{lat}}^2/2` the cell area,

453 :math:`|\mathbf{K}| = 4\pi/3a_{\mathrm{lat}}`, :math:`a_{\mathrm{lat}} = \sqrt3 a`),

454 the Dirac velocity :math:`\hbar v = \frac32|t|a`, the ratio

455 :math:`\alpha = w/\hbar vk_\theta` (:math:`k_\theta = 2|\mathbf{K}|\sin\frac\theta2`),

456 and the renormalized velocity at the Dirac points of the moiré bands

457 (Bistritzer and MacDonald, PNAS 108, 12233 (2011)),

458

459 .. math::

460

461 \frac{v^*}{v} = \frac{1-3\alpha^2}{1+6\alpha^2}\, ,

462

463 which vanishes at the first magic angle, :math:`\alpha = 1/\sqrt3`.

464

465 :param theta: Positive real. Twist angle, in radians.

466 :param a: Positive real. Default value 0.142. Nearest-neighbour distance.

467 :param d: Positive real. Default value 0.335. Interlayer distance.

468 :param t: Real number. Default value -2.7. Nearest-neighbour hopping.

469 :param interlayer: Callable. Default value None (*pz_hopping*). Interlayer hopping.

470 :param interlayer_scale: Positive real. Default value 1. See *twisted_bilayer*.

471

472 :returns:

473 * **w** -- Real number. Interlayer tunnelling.

474 * **alpha** -- Real number.

475 * **v_ratio** -- Real number. :math:`v^*/v`.

476 '''

477 error_handling.positive_real(theta, 'theta')

478 error_handling.positive_real(a, 'a')

479 error_handling.positive_real(d, 'd')

480 error_handling.real_number(t, 't')

481 error_handling.positive_real(interlayer_scale, 'interlayer_scale')

482 if interlayer is None:

483 interlayer = pz_hopping

484 error_handling.is_callable(interlayer, 'interlayer')

485 a_lat = np.sqrt(3) * a

486 k_d = 4 * PI / (3 * a_lat)

487 omega = np.sqrt(3) / 2 * a_lat ** 2

488

489 def integrand(r):

490 return r * float(np.real(interlayer(np.array([[r, 0., d]]))[0])) * j0(k_d * r)

491 w = interlayer_scale * 2 * PI * quad(integrand, 0., 20 * d, limit=200)[0] / omega

492 k_theta = 2 * k_d * np.sin(theta / 2)

493 alpha = abs(w) / (1.5 * abs(t) * a * k_theta)

494 return float(abs(w)), float(alpha), float((1 - 3 * alpha ** 2) / (1 + 6 * alpha ** 2))