Coverage for tbkit/topology.py: 100%
120 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +0000
1r"""
2Topological diagnostics beyond the Brillouin-zone integrals of
3:class:`tbkit.kspace.KSpace`:
5* **Bott index** (*bott_index*): the Chern number of a finite sample on a
6 torus, from its projector on the occupied states and the exponentiated
7 position operators (Loring and Hastings, EPL 92, 67004 (2010)). It needs
8 no Brillouin zone, so it applies to disordered samples and quasicrystals,
9 and complements the local Chern marker of
10 :meth:`tbkit.system.System.get_local_chern_marker`.
11* **Entanglement spectrum** (*entanglement_spectrum*): the eigenvalues of
12 the correlation matrix of the occupied states restricted to a region of
13 a finite sample (Peschel 2003; Li and Haldane 2008). Its k-resolved
14 version for ribbons and supercells is
15 :meth:`tbkit.kspace.KSpace.entanglement_spectrum`.
16* **Weyl points** (*find_weyl_points*): the band touchings of a 3D model
17 and their chiralities, the Berry flux through a small sphere around each
18 one (Wan, Turner, Vishwanath and Savrasov, Phys. Rev. B 83, 205101
19 (2011)) -- the Hermitian counterpart of
20 :func:`tbkit.exceptional.find_exceptional_points`.
21"""
22from __future__ import annotations
24import numpy as np
25from numpy.typing import ArrayLike, NDArray
26import scipy.linalg as LA
27import scipy.sparse as sparse
29import tbkit.error_handling as error_handling
30from tbkit.kspace import KSpace
33PI = np.pi
34_MESH_SHIFT = (0.3183098861837907, 0.2718281828459045, 0.1414213562373095) # off symmetry planes
37def _occupied(ham, e_fermi: float) -> NDArray[np.complex128]:
38 '''
39 Private function. Validate a Hermitian Hamiltonian (dense or sparse)
40 and return its eigenvectors below *e_fermi* (columns).
41 '''
42 h = sparse.csr_matrix(ham, dtype='c16')
43 error_handling.square_matrix(h, 'ham')
44 error_handling.hermitian(h)
45 error_handling.real_number(e_fermi, 'e_fermi')
46 en, vec = LA.eigh(h.toarray())
47 return vec[:, en < e_fermi]
50def bott_index(ham, positions: ArrayLike, cell: ArrayLike, e_fermi: float = 0.) -> float:
51 r'''
52 Get the Bott index of the states below *e_fermi* of a finite sample with
53 periodic boundary conditions (a torus spanned by :math:`\mathbf{L}_1,
54 \mathbf{L}_2`; Loring and Hastings, EPL 92, 67004 (2010)):
56 .. math::
58 B = \frac{1}{2\pi}\,\mathrm{Im}\,\mathrm{Tr}\log\left(
59 \tilde V\tilde U\tilde V^\dagger\tilde U^\dagger\right)\, ,\qquad
60 \tilde U = P e^{2\pi i X_1} P\, ,\quad \tilde V = P e^{2\pi i X_2} P\, ,
62 with :math:`P` the projector on the occupied states and
63 :math:`X_{1,2}` the fractional coordinates of the orbitals along
64 :math:`\mathbf{L}_{1,2}`. The exponentials are well defined on the torus,
65 where the position operator is not. :math:`B` is an integer whenever the
66 occupied states are localized or gapped, equals the Chern number of a
67 clean sample (in the sign convention of *KSpace.chern_number*), and
68 keeps working with disorder until the mobility gap closes. The torus
69 comes from *KSpace.finite_ham* with ``periodic=True`` and the positions
70 from *tbkit.higher_order.flake_positions*.
72 Parameters
73 ----------
74 ham : Square Hermitian matrix (dense or sparse), shape (N, N)
75 positions : Real array, shape (N, d)
76 Position of the orbital of each row.
77 cell : Real array, shape (2, d)
78 The two vectors
79 :math:`\mathbf{L}_1, \mathbf{L}_2` spanning the torus (e.g.
80 :math:`N_1\mathbf{a}_1` and :math:`N_2\mathbf{a}_2`).
81 e_fermi : Real number, default 0
82 Fermi energy.
84 Returns
85 -------
86 bott : Real number, close to an integer
88 Examples
89 --------
90 Example usage::
92 n = 12
93 ham = ks.finite_ham((n, n), periodic=True)
94 pos = flake_positions(ks, (n, n))
95 cell = [n * np.array(a) for a in ks.lat.prim_vec]
96 bott = bott_index(ham, pos, cell)
97 '''
98 occ = _occupied(ham, e_fermi)
99 positions = np.asarray(positions, dtype='f8')
100 cell = np.asarray(cell, dtype='f8')
101 error_handling.site_positions(positions, len(occ))
102 error_handling.torus_cell(cell, positions.shape[1])
103 fracs = positions @ np.linalg.pinv(cell)
104 u = occ.conj().T @ (np.exp(2j * PI * fracs[:, 0])[:, None] * occ)
105 v = occ.conj().T @ (np.exp(2j * PI * fracs[:, 1])[:, None] * occ)
106 w = v @ u @ v.conj().T @ u.conj().T
107 return float(np.sum(np.angle(np.linalg.eigvals(w))) / (2 * PI))
110def entanglement_spectrum(ham, region: list[int], e_fermi: float = 0.) -> NDArray[np.float64]:
111 r'''
112 Get the entanglement spectrum of the states below *e_fermi* of a finite
113 sample: the eigenvalues :math:`\xi_n \in [0, 1]` of the correlation
114 matrix :math:`C_{ij} = \langle c_i^\dagger c_j\rangle` restricted to the
115 orbitals :math:`i, j` of *region* (Peschel, J. Phys. A 36, L205 (2003)).
116 The entanglement Hamiltonian of the region is quadratic, with
117 single-particle energies :math:`\epsilon_n = \ln[(1-\xi_n)/\xi_n]`, and
118 the entanglement entropy is
119 :math:`S = -\sum_n[\xi_n\ln\xi_n + (1-\xi_n)\ln(1-\xi_n)]`.
120 Its low-lying part mirrors the edge spectrum at the cut (Li and Haldane,
121 Phys. Rev. Lett. 101, 010504 (2008)): a topological phase leaves
122 entanglement modes inside (0, 1), e.g. a mode at exactly 1/2 at each cut
123 of a chiral-symmetric chain -- without any physical edge, so it is a
124 property of the ground state alone. See
125 *KSpace.entanglement_spectrum* for the k-resolved version.
127 Parameters
128 ----------
129 ham : Square Hermitian matrix (dense or sparse), shape (N, N)
130 region : List of row indices (the region :math:`A`)
131 e_fermi : Real number, default 0
132 Fermi energy.
134 Returns
135 -------
136 xi : Real ndarray, shape (len(region),), sorted ascending
137 '''
138 occ = _occupied(ham, e_fermi)
139 error_handling.region(region, len(occ))
140 occ = occ[region]
141 return np.linalg.eigvalsh(occ @ occ.conj().T)
144class WeylPoints():
145 r'''
146 Weyl points of a 3D model, see *find_weyl_points*.
148 Attributes
149 ----------
150 k : Real ndarray, shape (n, 3)
151 Positions, in the k coordinates of
152 *KSpace.get_ham*, folded to fractional coordinates in [-1/2, 1/2).
153 chirality : Integer ndarray, shape (n,)
154 Berry flux of the bands
155 below the touching out of a small sphere around each point, over
156 :math:`2\pi`: :math:`\pm1` for a simple Weyl point.
157 energy : Real ndarray, shape (n,)
158 Energy of the touching.
159 gap : Real ndarray, shape (n,)
160 Remaining gap at the refined point
161 (zero up to the tolerance).
162 '''
164 def __init__(self, k, chirality, energy, gap) -> None:
165 self.k = k
166 self.chirality = chirality
167 self.energy = energy
168 self.gap = gap
170 @property
171 def total_chirality(self) -> int:
172 '''
173 Sum of the chiralities: zero over a whole Brillouin zone
174 (Nielsen-Ninomiya doubling theorem).
175 '''
176 return int(np.sum(self.chirality))
179def _gap(ks: KSpace, kpts: NDArray[np.float64], n: int) -> NDArray[np.float64]:
180 '''
181 Private function. Gap E_{n+1} - E_n at the k-points *kpts* (shape (m, 3)).
182 '''
183 en = ks._eigs(kpts)
184 return en[:, n + 1] - en[:, n]
187def _newton(ks: KSpace, k0: NDArray[np.float64], n: int, h: float, tol: float, max_iter: int):
188 r'''
189 Private function. Newton's method on the touching of bands n and n+1:
190 at the current point, the two bands span a 2D subspace :math:`W`, and
191 :math:`W^\dagger H(\mathbf{k})W - \mathrm{tr}/2 = \mathbf{d}\cdot\boldsymbol\sigma`
192 near it; the step solves :math:`\mathbf{d} + J\delta\mathbf{k} = 0`
193 (central-difference Jacobian). It returns the point and its gap.
194 '''
195 sigma = np.array([[[0., 1.], [1., 0.]], [[0., -1j], [1j, 0.]], [[1., 0.], [0., -1.]]])
196 k = np.array(k0, dtype='f8')
197 for _ in range(max_iter):
198 _, vec = ks._eigs(k[None], eigenvec=True)
199 w = vec[0][:, n:n + 2]
201 def d(kk):
202 pts = np.atleast_2d(kk)
203 heff = w.conj().T @ ks._hams(pts) @ w
204 return np.einsum('aij,kji->ka', sigma, heff).real / 2
205 shifts = h * np.eye(3)
206 dd = d(np.concatenate([k + shifts, k - shifts]))
207 jac = ((dd[:3] - dd[3:]) / (2 * h)).T
208 step = np.linalg.lstsq(jac, -d(k)[0], rcond=None)[0]
209 k = k + step
210 if np.linalg.norm(step) < tol:
211 break
212 return k, float(_gap(ks, k[None], n)[0])
215def _sphere_chirality(ks: KSpace, center: NDArray[np.float64], radius: float,
216 bands: list[int], n_theta: int = 16) -> int:
217 r'''
218 Private function. Berry flux of *bands* out of the sphere of *radius*
219 around *center*, over :math:`2\pi`: the Fukui-Hatsugai-Suzuki
220 plaquettes of a (theta, phi) mesh, poles included, oriented outward
221 (the sign convention of *KSpace.berry_curvature*).
222 '''
223 n_phi = 2 * n_theta
224 theta = np.linspace(0., PI, n_theta + 1)
225 phi = 2 * PI * np.arange(n_phi) / n_phi
226 tt, pp = np.meshgrid(theta, phi, indexing='ij')
227 dirs = np.stack([np.sin(tt) * np.cos(pp), np.sin(tt) * np.sin(pp), np.cos(tt)], axis=-1)
228 pts = center + radius * dirs.reshape(-1, 3)
229 v1 = ks._subspaces(pts, bands).reshape(n_theta + 1, n_phi, ks.norb, len(bands))
230 # each pole is one point: use the same basis for all its phi
231 v1[0] = v1[0, 0]
232 v1[-1] = v1[-1, 0]
234 def overlap(a, b):
235 return np.linalg.det(a.conj().swapaxes(-1, -2) @ b)
236 a, b = v1[:-1], v1[1:]
237 a_next, b_next = np.roll(a, -1, axis=1), np.roll(b, -1, axis=1)
238 # corners (theta, phi) -> (theta + d, phi) -> (theta + d, phi + d) -> (theta, phi + d):
239 # counterclockwise about e_theta x e_phi = e_r, the outward normal
240 link = overlap(a, b) * overlap(b, b_next) * overlap(b_next, a_next) * overlap(a_next, a)
241 return int(np.rint(-np.angle(link).sum() / (2 * PI)))
244def find_weyl_points(
245 ks: KSpace, bands: int | list[int], nk: int = 20, radius: float | None = None,
246 tol: float = 1e-10, max_iter: int = 50,
247) -> WeylPoints:
248 r'''
249 Locate the Weyl points of a 3D model: the points where the highest band
250 of *bands* (band :math:`n`) touches the next one, and their chirality.
252 1. The gap :math:`E_{n+1} - E_n` is computed on an :math:`nk^3` mesh of
253 the Brillouin zone (shifted off the high-symmetry planes); its local
254 minima (over the 26 neighbours, periodically) are the candidates.
255 2. Each candidate is refined by Newton's method on the two-band
256 Hamiltonian :math:`\mathbf{d}(\mathbf{k})\cdot\boldsymbol\sigma` of bands
257 :math:`n, n+1`, whose zeros are the touchings; candidates whose gap
258 does not close (avoided crossings) are dropped, and points equal up
259 to a reciprocal lattice vector merged.
260 3. The chirality of each point is the Berry flux of *bands* out of a
261 small sphere around it, divided by :math:`2\pi` -- the monopole
262 charge, :math:`\pm1` for a linear Weyl point. With the sign
263 convention of *KSpace.berry_curvature*, the Chern number of the
264 planes :math:`k_l = \mathrm{const}` (*KSpace.chern_number*) jumps by
265 the chirality as the plane crosses the point along
266 :math:`+\mathbf{b}_l`.
268 The chiralities add up to zero (*total_chirality*): Weyl points come in
269 pairs of opposite chirality (Nielsen-Ninomiya). The mesh must resolve
270 the gap minima: on too coarse a mesh, two candidates may refine to the
271 same point and miss another one, which a nonzero *total_chirality*
272 reveals (increase *nk*).
273 Nodal lines (degeneracies along curves) are not isolated points and are
274 outside the scope of this finder.
276 Parameters
277 ----------
278 ks : **KSpace** instance with three primitive vectors (Hermitian)
279 bands
280 Band index, or list of band indices (e.g. the occupied
281 bands): the touching is between ``max(bands)`` and ``max(bands) + 1``.
282 nk : Positive integer, default 20
283 Mesh points along each
284 reciprocal vector.
285 radius : Positive real
286 Default value None: a fifth of the mesh
287 step. Radius of the sphere for the chirality (in units of k); it
288 must enclose no other Weyl point.
289 tol : Positive real, default 1e-10
290 Newton's method stops
291 when a step is shorter than *tol*; a point is kept if its gap is
292 below :math:`10^4` *tol* times the bandwidth.
293 max_iter : Positive integer, default 50
294 Newton iterations.
296 Returns
297 -------
298 weyl : **WeylPoints** instance (*k*, *chirality*, *energy*, *gap*)
299 '''
300 error_handling.kspace(ks, KSpace)
301 error_handling.dim_exact(ks.dim, 3)
302 error_handling.hermitian_kspace(ks.is_hermitian(), ks._overlap_hop)
303 if isinstance(bands, int):
304 bands = [bands]
305 error_handling.band_indices(bands, ks.norb)
306 n = max(bands)
307 error_handling.band_below_top(n, ks.norb)
308 error_handling.positive_int(nk, 'nk')
309 step = min(np.linalg.norm(b) for b in ks.rec_vec_k) / nk
310 if radius is None:
311 radius = step / 5
312 error_handling.positive_real(radius, 'radius')
313 error_handling.positive_real(tol, 'tol')
314 error_handling.positive_int(max_iter, 'max_iter')
315 fr = (np.arange(nk)[:, None] + np.array(_MESH_SHIFT)[None, :]) / nk
316 fracs = np.stack(np.meshgrid(fr[:, 0], fr[:, 1], fr[:, 2], indexing='ij'), axis=-1).reshape(-1, 3)
317 kmesh = fracs @ ks.rec_vec_k
318 en = ks._eigs(kmesh)
319 scale = float(en.max() - en.min()) or 1.
320 gap = (en[:, n + 1] - en[:, n]).reshape(nk, nk, nk)
321 is_min = np.ones_like(gap, bool)
322 for shift in np.ndindex(3, 3, 3):
323 if shift != (1, 1, 1):
324 is_min &= gap <= np.roll(gap, np.array(shift) - 1, axis=(0, 1, 2))
325 inv = np.linalg.inv(ks.rec_vec_k)
326 points: list = []
327 gaps: list = []
328 for idx in np.argwhere(is_min):
329 k, g = _newton(ks, kmesh[np.ravel_multi_index(idx, gap.shape)], n,
330 1e-5 * step, tol, max_iter)
331 if not np.isfinite(g) or g > 1e4 * tol * scale:
332 continue
333 # fold into the zone centred on Gamma, and merge points equal up to
334 # a reciprocal vector
335 frac = (k @ inv + 0.5) % 1. - 0.5
336 if any(np.linalg.norm((frac - p + 0.5) % 1. - 0.5) < 1e-6 for p in points):
337 continue
338 points.append(frac)
339 gaps.append(g)
340 k_w = np.array(points).reshape(-1, 3) @ ks.rec_vec_k
341 chirality = np.array([_sphere_chirality(ks, k, radius, bands) for k in k_w], int)
342 energy = np.array([ks._eigs(k[None])[0, n] for k in k_w])
343 return WeylPoints(k_w, chirality, energy.reshape(-1), np.array(gaps))