Coverage for tbkit/topology.py: 100%
120 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 06:13 +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 :param ham: Square Hermitian matrix (dense or sparse), shape (N, N).
73 :param positions: Real array, shape (N, d). Position of the orbital of each row.
74 :param cell: Real array, shape (2, d). The two vectors
75 :math:`\mathbf{L}_1, \mathbf{L}_2` spanning the torus (e.g.
76 :math:`N_1\mathbf{a}_1` and :math:`N_2\mathbf{a}_2`).
77 :param e_fermi: Real number. Default value 0. Fermi energy.
79 :returns:
80 * **bott** -- Real number, close to an integer.
82 Example usage::
84 n = 12
85 ham = ks.finite_ham((n, n), periodic=True)
86 pos = flake_positions(ks, (n, n))
87 cell = [n * np.array(a) for a in ks.lat.prim_vec]
88 bott = bott_index(ham, pos, cell)
89 '''
90 occ = _occupied(ham, e_fermi)
91 positions = np.asarray(positions, dtype='f8')
92 cell = np.asarray(cell, dtype='f8')
93 error_handling.site_positions(positions, len(occ))
94 error_handling.torus_cell(cell, positions.shape[1])
95 fracs = positions @ np.linalg.pinv(cell)
96 u = occ.conj().T @ (np.exp(2j * PI * fracs[:, 0])[:, None] * occ)
97 v = occ.conj().T @ (np.exp(2j * PI * fracs[:, 1])[:, None] * occ)
98 w = v @ u @ v.conj().T @ u.conj().T
99 return float(np.sum(np.angle(np.linalg.eigvals(w))) / (2 * PI))
102def entanglement_spectrum(ham, region: list[int], e_fermi: float = 0.) -> NDArray[np.float64]:
103 r'''
104 Get the entanglement spectrum of the states below *e_fermi* of a finite
105 sample: the eigenvalues :math:`\xi_n \in [0, 1]` of the correlation
106 matrix :math:`C_{ij} = \langle c_i^\dagger c_j\rangle` restricted to the
107 orbitals :math:`i, j` of *region* (Peschel, J. Phys. A 36, L205 (2003)).
108 The entanglement Hamiltonian of the region is quadratic, with
109 single-particle energies :math:`\epsilon_n = \ln[(1-\xi_n)/\xi_n]`, and
110 the entanglement entropy is
111 :math:`S = -\sum_n[\xi_n\ln\xi_n + (1-\xi_n)\ln(1-\xi_n)]`.
112 Its low-lying part mirrors the edge spectrum at the cut (Li and Haldane,
113 Phys. Rev. Lett. 101, 010504 (2008)): a topological phase leaves
114 entanglement modes inside (0, 1), e.g. a mode at exactly 1/2 at each cut
115 of a chiral-symmetric chain -- without any physical edge, so it is a
116 property of the ground state alone. See
117 *KSpace.entanglement_spectrum* for the k-resolved version.
119 :param ham: Square Hermitian matrix (dense or sparse), shape (N, N).
120 :param region: List of row indices (the region :math:`A`).
121 :param e_fermi: Real number. Default value 0. Fermi energy.
123 :returns:
124 * **xi** -- Real ndarray, shape (len(region),), sorted ascending.
125 '''
126 occ = _occupied(ham, e_fermi)
127 error_handling.region(region, len(occ))
128 occ = occ[region]
129 return np.linalg.eigvalsh(occ @ occ.conj().T)
132class WeylPoints():
133 r'''
134 Weyl points of a 3D model, see *find_weyl_points*.
136 :ivar k: Real ndarray, shape (n, 3). Positions, in the k coordinates of
137 *KSpace.get_ham*, folded to fractional coordinates in [-1/2, 1/2).
138 :ivar chirality: Integer ndarray, shape (n,). Berry flux of the bands
139 below the touching out of a small sphere around each point, over
140 :math:`2\pi`: :math:`\pm1` for a simple Weyl point.
141 :ivar energy: Real ndarray, shape (n,). Energy of the touching.
142 :ivar gap: Real ndarray, shape (n,). Remaining gap at the refined point
143 (zero up to the tolerance).
144 '''
146 def __init__(self, k, chirality, energy, gap) -> None:
147 self.k = k
148 self.chirality = chirality
149 self.energy = energy
150 self.gap = gap
152 @property
153 def total_chirality(self) -> int:
154 '''
155 Sum of the chiralities: zero over a whole Brillouin zone
156 (Nielsen-Ninomiya doubling theorem).
157 '''
158 return int(np.sum(self.chirality))
161def _gap(ks: KSpace, kpts: NDArray[np.float64], n: int) -> NDArray[np.float64]:
162 '''
163 Private function. Gap E_{n+1} - E_n at the k-points *kpts* (shape (m, 3)).
164 '''
165 en = ks._eigs(kpts)
166 return en[:, n + 1] - en[:, n]
169def _newton(ks: KSpace, k0: NDArray[np.float64], n: int, h: float, tol: float, max_iter: int):
170 r'''
171 Private function. Newton's method on the touching of bands n and n+1:
172 at the current point, the two bands span a 2D subspace :math:`W`, and
173 :math:`W^\dagger H(\mathbf{k})W - \mathrm{tr}/2 = \mathbf{d}\cdot\boldsymbol\sigma`
174 near it; the step solves :math:`\mathbf{d} + J\delta\mathbf{k} = 0`
175 (central-difference Jacobian). It returns the point and its gap.
176 '''
177 sigma = np.array([[[0., 1.], [1., 0.]], [[0., -1j], [1j, 0.]], [[1., 0.], [0., -1.]]])
178 k = np.array(k0, dtype='f8')
179 for _ in range(max_iter):
180 _, vec = ks._eigs(k[None], eigenvec=True)
181 w = vec[0][:, n:n + 2]
183 def d(kk):
184 pts = np.atleast_2d(kk)
185 heff = w.conj().T @ ks._hams(pts) @ w
186 return np.einsum('aij,kji->ka', sigma, heff).real / 2
187 shifts = h * np.eye(3)
188 dd = d(np.concatenate([k + shifts, k - shifts]))
189 jac = ((dd[:3] - dd[3:]) / (2 * h)).T
190 step = np.linalg.lstsq(jac, -d(k)[0], rcond=None)[0]
191 k = k + step
192 if np.linalg.norm(step) < tol:
193 break
194 return k, float(_gap(ks, k[None], n)[0])
197def _sphere_chirality(ks: KSpace, center: NDArray[np.float64], radius: float,
198 bands: list[int], n_theta: int = 16) -> int:
199 r'''
200 Private function. Berry flux of *bands* out of the sphere of *radius*
201 around *center*, over :math:`2\pi`: the Fukui-Hatsugai-Suzuki
202 plaquettes of a (theta, phi) mesh, poles included, oriented outward
203 (the sign convention of *KSpace.berry_curvature*).
204 '''
205 n_phi = 2 * n_theta
206 theta = np.linspace(0., PI, n_theta + 1)
207 phi = 2 * PI * np.arange(n_phi) / n_phi
208 tt, pp = np.meshgrid(theta, phi, indexing='ij')
209 dirs = np.stack([np.sin(tt) * np.cos(pp), np.sin(tt) * np.sin(pp), np.cos(tt)], axis=-1)
210 pts = center + radius * dirs.reshape(-1, 3)
211 v1 = ks._subspaces(pts, bands).reshape(n_theta + 1, n_phi, ks.norb, len(bands))
212 # each pole is one point: use the same basis for all its phi
213 v1[0] = v1[0, 0]
214 v1[-1] = v1[-1, 0]
216 def overlap(a, b):
217 return np.linalg.det(a.conj().swapaxes(-1, -2) @ b)
218 a, b = v1[:-1], v1[1:]
219 a_next, b_next = np.roll(a, -1, axis=1), np.roll(b, -1, axis=1)
220 # corners (theta, phi) -> (theta + d, phi) -> (theta + d, phi + d) -> (theta, phi + d):
221 # counterclockwise about e_theta x e_phi = e_r, the outward normal
222 link = overlap(a, b) * overlap(b, b_next) * overlap(b_next, a_next) * overlap(a_next, a)
223 return int(np.rint(-np.angle(link).sum() / (2 * PI)))
226def find_weyl_points(
227 ks: KSpace, bands: int | list[int], nk: int = 20, radius: float | None = None,
228 tol: float = 1e-10, max_iter: int = 50,
229) -> WeylPoints:
230 r'''
231 Locate the Weyl points of a 3D model: the points where the highest band
232 of *bands* (band :math:`n`) touches the next one, and their chirality.
234 1. The gap :math:`E_{n+1} - E_n` is computed on an :math:`nk^3` mesh of
235 the Brillouin zone (shifted off the high-symmetry planes); its local
236 minima (over the 26 neighbours, periodically) are the candidates.
237 2. Each candidate is refined by Newton's method on the two-band
238 Hamiltonian :math:`\mathbf{d}(\mathbf{k})\cdot\boldsymbol\sigma` of bands
239 :math:`n, n+1`, whose zeros are the touchings; candidates whose gap
240 does not close (avoided crossings) are dropped, and points equal up
241 to a reciprocal lattice vector merged.
242 3. The chirality of each point is the Berry flux of *bands* out of a
243 small sphere around it, divided by :math:`2\pi` -- the monopole
244 charge, :math:`\pm1` for a linear Weyl point. With the sign
245 convention of *KSpace.berry_curvature*, the Chern number of the
246 planes :math:`k_l = \mathrm{const}` (*KSpace.chern_number*) jumps by
247 the chirality as the plane crosses the point along
248 :math:`+\mathbf{b}_l`.
250 The chiralities add up to zero (*total_chirality*): Weyl points come in
251 pairs of opposite chirality (Nielsen-Ninomiya). The mesh must resolve
252 the gap minima: on too coarse a mesh, two candidates may refine to the
253 same point and miss another one, which a nonzero *total_chirality*
254 reveals (increase *nk*).
255 Nodal lines (degeneracies along curves) are not isolated points and are
256 outside the scope of this finder.
258 :param ks: **KSpace** instance with three primitive vectors (Hermitian).
259 :param bands: Band index, or list of band indices (e.g. the occupied
260 bands): the touching is between ``max(bands)`` and ``max(bands) + 1``.
261 :param nk: Positive integer. Default value 20. Mesh points along each
262 reciprocal vector.
263 :param radius: Positive real. Default value None: a fifth of the mesh
264 step. Radius of the sphere for the chirality (in units of k); it
265 must enclose no other Weyl point.
266 :param tol: Positive real. Default value 1e-10. Newton's method stops
267 when a step is shorter than *tol*; a point is kept if its gap is
268 below :math:`10^4` *tol* times the bandwidth.
269 :param max_iter: Positive integer. Default value 50. Newton iterations.
271 :returns:
272 * **weyl** -- **WeylPoints** instance (*k*, *chirality*, *energy*, *gap*).
273 '''
274 error_handling.kspace(ks, KSpace)
275 error_handling.dim_exact(ks.dim, 3)
276 error_handling.hermitian_kspace(ks.is_hermitian(), ks._overlap_hop)
277 if isinstance(bands, int):
278 bands = [bands]
279 error_handling.band_indices(bands, ks.norb)
280 n = max(bands)
281 error_handling.band_below_top(n, ks.norb)
282 error_handling.positive_int(nk, 'nk')
283 step = min(np.linalg.norm(b) for b in ks.rec_vec_k) / nk
284 if radius is None:
285 radius = step / 5
286 error_handling.positive_real(radius, 'radius')
287 error_handling.positive_real(tol, 'tol')
288 error_handling.positive_int(max_iter, 'max_iter')
289 fr = (np.arange(nk)[:, None] + np.array(_MESH_SHIFT)[None, :]) / nk
290 fracs = np.stack(np.meshgrid(fr[:, 0], fr[:, 1], fr[:, 2], indexing='ij'), axis=-1).reshape(-1, 3)
291 kmesh = fracs @ ks.rec_vec_k
292 en = ks._eigs(kmesh)
293 scale = float(en.max() - en.min()) or 1.
294 gap = (en[:, n + 1] - en[:, n]).reshape(nk, nk, nk)
295 is_min = np.ones_like(gap, bool)
296 for shift in np.ndindex(3, 3, 3):
297 if shift != (1, 1, 1):
298 is_min &= gap <= np.roll(gap, np.array(shift) - 1, axis=(0, 1, 2))
299 inv = np.linalg.inv(ks.rec_vec_k)
300 points: list = []
301 gaps: list = []
302 for idx in np.argwhere(is_min):
303 k, g = _newton(ks, kmesh[np.ravel_multi_index(idx, gap.shape)], n,
304 1e-5 * step, tol, max_iter)
305 if not np.isfinite(g) or g > 1e4 * tol * scale:
306 continue
307 # fold into the zone centred on Gamma, and merge points equal up to
308 # a reciprocal vector
309 frac = (k @ inv + 0.5) % 1. - 0.5
310 if any(np.linalg.norm((frac - p + 0.5) % 1. - 0.5) < 1e-6 for p in points):
311 continue
312 points.append(frac)
313 gaps.append(g)
314 k_w = np.array(points).reshape(-1, 3) @ ks.rec_vec_k
315 chirality = np.array([_sphere_chirality(ks, k, radius, bands) for k in k_w], int)
316 energy = np.array([ks._eigs(k[None])[0, n] for k in k_w])
317 return WeylPoints(k_w, chirality, energy.reshape(-1), np.array(gaps))