Coverage for tbkit/kpm.py: 100%
142 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 06:13 +0000
1r"""
2The kernel polynomial method (KPM; Silver and Roder 1994, Weisse,
3Wellein, Alvermann and Fehske, Rev. Mod. Phys. 78, 275 (2006)): spectral
4quantities of a large sparse Hamiltonian from its Chebyshev moments,
5without diagonalizing it. The spectrum is first mapped into
6:math:`(-1, 1)`, :math:`\tilde H = (H - b)/a`; then
8.. math::
10 \delta(E - H) = \sum_{m=0}^{M-1} c_m(E)\, T_m(\tilde H)\, ,\qquad
11 c_m(E) = \frac{(2-\delta_{m0})\,g_m\,T_m(\epsilon)}{a\pi\sqrt{1-\epsilon^2}}\, ,
13with :math:`\epsilon = (E-b)/a`, :math:`T_m` the Chebyshev polynomials
14(computed by the recursion :math:`T_{m+1}(\tilde H)|r\rangle =
152\tilde H\,T_m(\tilde H)|r\rangle - T_{m-1}(\tilde H)|r\rangle`, one sparse
16matrix-vector product each), and :math:`g_m` a kernel that damps the Gibbs
17oscillations of the truncated series (Jackson: an almost Gaussian
18broadening of width about :math:`\pi a/M`). Traces are estimated with a
19few random-phase vectors, :math:`\mathrm{Tr}\,A \approx
20\frac1R\sum_r\langle r|A|r\rangle`, whose error falls as :math:`1/\sqrt{RN}`:
21the cost is :math:`O(MRN)` for :math:`N` sites, so millions of sites are
22within reach.
24The Hamiltonian must be Hermitian (*System.ham*, a sparse matrix).
25"""
26from __future__ import annotations
28import numpy as np
29from numpy.typing import ArrayLike, NDArray
30import scipy.sparse as sparse
31from scipy.integrate import trapezoid
33import tbkit.error_handling as error_handling
34import tbkit.occupation as occupation
37PI = np.pi
40def spectral_bounds(ham) -> tuple[float, float]:
41 r'''
42 Get bounds on the spectrum of a Hermitian matrix from Gershgorin's
43 circle theorem: every eigenvalue lies within
44 :math:`H_{ii} \pm \sum_{j\neq i}|H_{ij}|` for some :math:`i`.
46 :param ham: Square sparse (or dense) Hermitian matrix.
48 :returns:
49 * **e_min**, **e_max** -- Real numbers.
50 '''
51 ham = sparse.csr_matrix(ham)
52 diag = ham.diagonal().real
53 radii = np.asarray(abs(ham).sum(axis=1)).ravel() - np.abs(ham.diagonal())
54 return float(np.min(diag - radii)), float(np.max(diag + radii))
57def jackson_kernel(n_moments: int) -> NDArray[np.float64]:
58 r'''
59 Get the Jackson kernel
61 .. math::
63 g_m = \frac{(M-m+1)\cos\frac{\pi m}{M+1}
64 + \sin\frac{\pi m}{M+1}\cot\frac{\pi}{M+1}}{M+1}\, ,
66 the optimal choice for spectral densities (positive, and resolution
67 about :math:`\pi/M` in the scaled energy).
69 :param n_moments: Positive integer :math:`M`.
71 :returns:
72 * **g** -- Real ndarray, shape (n_moments,).
73 '''
74 error_handling.positive_int(n_moments, 'n_moments')
75 m = np.arange(n_moments)
76 q = PI / (n_moments + 1)
77 return ((n_moments - m + 1) * np.cos(q * m) + np.sin(q * m) / np.tan(q)) / (n_moments + 1)
80def lorentz_kernel(n_moments: int, lam: float = 4.) -> NDArray[np.float64]:
81 r'''
82 Get the Lorentz kernel :math:`g_m = \sinh[\lambda(1 - m/M)]/\sinh\lambda`,
83 which gives Lorentzian broadening (the better choice for Green's
84 functions).
86 :param n_moments: Positive integer :math:`M`.
87 :param lam: Positive real. Default value 4.
89 :returns:
90 * **g** -- Real ndarray, shape (n_moments,).
91 '''
92 error_handling.positive_int(n_moments, 'n_moments')
93 error_handling.positive_real(lam, 'lam')
94 return np.sinh(lam * (1. - np.arange(n_moments) / n_moments)) / np.sinh(lam)
97def _setup(
98 ham, n_moments, kernel, bounds, e_grid, n_grid=1001,
99) -> tuple[sparse.csr_matrix, NDArray[np.float64], NDArray[np.float64], float, float]:
100 '''
101 Private function. Validate, rescale the Hamiltonian, and build the
102 coefficients c_m(E) (shape (n_moments, len(e_grid))).
103 '''
104 ham = sparse.csr_matrix(ham, dtype='c16')
105 error_handling.square_matrix(ham, 'ham')
106 error_handling.hermitian(ham)
107 error_handling.positive_int(n_moments, 'n_moments')
108 error_handling.kpm_kernel(kernel)
109 if bounds is None:
110 bounds = spectral_bounds(ham)
111 error_handling.lims(bounds)
112 e_min, e_max = bounds
113 a = (e_max - e_min) / (2. - 0.02) # keep the scaled spectrum inside (-0.99, 0.99)
114 b = (e_max + e_min) / 2.
115 if e_grid is None:
116 e_grid = b + a * np.cos(PI * (np.arange(n_grid) + 0.5) / n_grid)[::-1]
117 e_grid = np.atleast_1d(np.asarray(e_grid, dtype='f8'))
118 error_handling.ndarray_empty(e_grid, 'e_grid')
119 eps = (e_grid - b) / a
120 error_handling.kpm_grid(eps)
121 g = jackson_kernel(n_moments) if kernel == 'jackson' else lorentz_kernel(n_moments)
122 tm = np.cos(np.arange(n_moments)[:, None] * np.arccos(eps)[None, :])
123 weight = np.where(np.arange(n_moments) == 0, 1., 2.)
124 coef = (weight * g)[:, None] * tm / (a * PI * np.sqrt(1. - eps ** 2))[None, :]
125 ham_s = (ham - b * sparse.identity(ham.shape[0], format='csr')) / a
126 return ham_s, e_grid, coef, a, b
129def _vectors(n: int, n_random: int | None, seed) -> NDArray[np.complex128]:
130 '''
131 Private function. Random-phase vectors (columns, |r_i| = 1), or, if
132 n_random is None, the full basis (the exact trace).
133 '''
134 if n_random is None:
135 return np.eye(n, dtype='c16')
136 error_handling.positive_int(n_random, 'n_random')
137 rng = np.random.default_rng(seed)
138 return np.exp(2j * PI * rng.random((n, n_random)))
141def _chebyshev(ham_s, vec: NDArray[np.complex128], n_moments: int):
142 '''
143 Private function. Yield T_m(H~) vec for m < n_moments, keeping only two
144 vectors in memory.
145 '''
146 prev, cur = None, vec
147 for m in range(n_moments):
148 yield cur
149 prev, cur = cur, (ham_s @ cur if m == 0 else 2 * (ham_s @ cur) - prev)
152def dos(
153 ham, n_moments: int = 256, n_random: int | None = 10, e_grid: ArrayLike | None = None,
154 kernel: str = 'jackson', seed=None, bounds: tuple[float, float] | None = None,
155) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
156 r'''
157 Get the density of states :math:`\rho(E) = \mathrm{Tr}\,\delta(E - H)`
158 (normalized, like *tbkit.dos.density_of_states*, to integrate to the
159 number of states).
161 :param ham: Square sparse Hermitian matrix (e.g. *System.ham*).
162 :param n_moments: Positive integer. Default value 256. Number of
163 Chebyshev moments (the energy resolution is about
164 :math:`\pi a/M`, :math:`2a` the spectral width).
165 :param n_random: Positive integer, or None. Default value 10. Number of
166 random vectors for the trace; None for the exact trace (N vectors:
167 small systems only).
168 :param e_grid: Real array. Default value None: 1001 Chebyshev nodes over
169 the spectrum. Energies, inside the spectral bounds.
170 :param kernel: String. Default value 'jackson'. 'jackson' or 'lorentz'.
171 :param seed: Default value None. Seed of the random vectors.
172 :param bounds: Pair of reals. Default value None (Gershgorin bounds, see
173 *spectral_bounds*). Bounds enclosing the whole spectrum.
175 :returns:
176 * **e_grid** -- Real ndarray.
177 * **rho** -- Real ndarray, same shape.
178 '''
179 ham_s, e_grid, coef, _, _ = _setup(ham, n_moments, kernel, bounds, e_grid)
180 vec = _vectors(ham_s.shape[0], n_random, seed)
181 mu = np.array([np.vdot(vec, t).real for t in _chebyshev(ham_s, vec, n_moments)])
182 if n_random is not None:
183 mu = mu / n_random # E<r|A|r> = Tr A for random phases
184 return e_grid, mu @ coef
187def dos_from_levels(
188 levels: ArrayLike, n_moments: int, bounds: tuple[float, float], e_grid: ArrayLike,
189 kernel: str = 'jackson',
190) -> NDArray[np.float64]:
191 r'''
192 Get the KPM density of states of a *known* spectrum: every level
193 broadened by exactly the truncated, kernel-damped Chebyshev series that
194 *dos* uses. It is the reference *dos* converges to (up to the random
195 trace), e.g. for the bands of a periodic model on a k-mesh.
197 :param levels: Real array. Energy levels.
198 :param n_moments: Positive integer. Number of Chebyshev moments.
199 :param bounds: Pair of reals. Spectral bounds (those of *dos*, e.g.
200 *spectral_bounds* of the Hamiltonian it is compared with).
201 :param e_grid: Real array. Energies.
202 :param kernel: String. Default value 'jackson'.
204 :returns:
205 * **rho** -- Real ndarray, same shape as *e_grid* (integrating to the number of levels).
206 '''
207 levels = np.asarray(levels, dtype='f8').ravel()
208 error_handling.ndarray_empty(levels, 'levels')
209 _, e_grid, coef, a, b = _setup(sparse.identity(1, format='csr'), n_moments, kernel, bounds, e_grid)
210 error_handling.kpm_grid((levels - b) / a)
211 tm = np.cos(np.arange(n_moments)[:, None] * np.arccos((levels - b) / a)[None, :])
212 return coef.T @ tm.sum(axis=1)
215def ldos(
216 ham, site: int, n_moments: int = 256, e_grid: ArrayLike | None = None,
217 kernel: str = 'jackson', bounds: tuple[float, float] | None = None,
218) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
219 r'''
220 Get the local density of states on one site (or orbital),
221 :math:`\rho_i(E) = \langle i|\delta(E-H)|i\rangle`, exactly (no random vectors).
223 :param ham: Square sparse Hermitian matrix.
224 :param site: Integer. Site (row) index.
225 :param n_moments: Positive integer. Default value 256.
226 :param e_grid: See *dos*.
227 :param kernel: See *dos*.
228 :param bounds: See *dos*.
230 :returns:
231 * **e_grid** -- Real ndarray.
232 * **rho** -- Real ndarray, same shape.
233 '''
234 ham_s, e_grid, coef, _, _ = _setup(ham, n_moments, kernel, bounds, e_grid)
235 error_handling.site_index(site, ham_s.shape[0])
236 vec = np.zeros(ham_s.shape[0], 'c16')
237 vec[site] = 1.
238 mu = np.array([t[site].real for t in _chebyshev(ham_s, vec, n_moments)])
239 return e_grid, mu @ coef
242def conductivity(
243 ham, x: ArrayLike, n_moments: int = 128, n_random: int | None = 10,
244 e_grid: ArrayLike | None = None, kernel: str = 'jackson', seed=None,
245 bounds: tuple[float, float] | None = None, area: float | None = None,
246) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
247 r'''
248 Get the zero-temperature DC conductivity along :math:`x` from the
249 Kubo-Greenwood formula,
251 .. math::
253 \sigma_{xx}(E) = \frac{\pi}{\Omega}\,
254 \mathrm{Tr}\left[v_x\,\delta(E-H)\,v_x\,\delta(E-H)\right]\, ,\qquad
255 v_x = i[H, x]\, ,
257 in units of :math:`e^2/\hbar` (:math:`\hbar = e = 1`), by the
258 two-dimensional Chebyshev expansion of Weisse et al. (the moments
259 :math:`\mu_{mn} = \mathrm{Tr}[v_x T_m(\tilde H) v_x T_n(\tilde H)]`,
260 :math:`O(M N)` memory, the DOS needing only :math:`O(N)`).
262 :param ham: Square sparse Hermitian matrix.
263 :param x: Real array, one coordinate per site (e.g. *lat.coor['x']*).
264 :param n_moments: Positive integer. Default value 128.
265 :param n_random: Positive integer, or None. Default value 10. See *dos*.
266 :param e_grid: See *dos*.
267 :param kernel: See *dos*.
268 :param seed: See *dos*.
269 :param bounds: See *dos*.
270 :param area: Positive real. Default value None (the number of sites):
271 the area (volume) :math:`\Omega` of the sample.
273 :returns:
274 * **e_grid** -- Real ndarray.
275 * **sigma** -- Real ndarray, same shape.
276 '''
277 ham_s, e_grid, coef, a, b = _setup(ham, n_moments, kernel, bounds, e_grid)
278 n = ham_s.shape[0]
279 x = np.asarray(x, dtype='f8')
280 error_handling.ndarray(x, 'x', n)
281 if area is None:
282 area = float(n)
283 error_handling.positive_real(area, 'area')
284 # v = i [H, x]: v_ij = i H_ij (x_j - x_i); the scale a comes back below
285 ham_c = sparse.csr_matrix(ham, dtype='c16')
286 vel = 1j * (ham_c @ sparse.diags(x) - sparse.diags(x) @ ham_c)
287 vec = _vectors(n, n_random, seed)
288 mu = np.zeros((n_moments, n_moments), 'c16')
289 for r in vec.T:
290 left = np.array(list(_chebyshev(ham_s, vel @ r, n_moments))) # T_m v |r>
291 right = np.array([vel @ t for t in _chebyshev(ham_s, r, n_moments)]) # v T_n |r>
292 mu += left.conj() @ right.T
293 if n_random is not None:
294 mu = mu / n_random
295 sigma = PI / area * np.einsum('me,mn,ne->e', coef, mu, coef).real
296 return e_grid, sigma
299def hall_conductivity(
300 ham, vx, vy, n_moments: int = 256, n_random: int | None = 10,
301 e_grid: ArrayLike | None = None, kernel: str = 'jackson', seed=None,
302 bounds: tuple[float, float] | None = None, area: float | None = None,
303 temperature: float = 0.,
304) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
305 r'''
306 Get the Hall conductivity :math:`\sigma_{xy}(\mu)` of a large sample, at
307 any Fermi level and temperature, from the Kubo-Bastin formula expanded
308 in Chebyshev polynomials (Garcia, Covaci and Rappoport, Phys. Rev. Lett.
309 114, 116602 (2015), Eqs. (1)-(5)). Bastin's formula,
311 .. math::
313 \sigma^{B}_{\alpha\beta}(\mu) = \frac{ie^2\hbar}{\Omega}\int d\varepsilon\,f(\varepsilon)\,
314 \mathrm{Tr}\left[v_\alpha\,\delta(\varepsilon-H)\,v_\beta\,\frac{dG^+}{d\varepsilon}
315 - v_\alpha\,\frac{dG^-}{d\varepsilon}\,v_\beta\,\delta(\varepsilon-H)\right]\, ,
317 :math:`G^\pm = (\varepsilon - H \pm i0)^{-1}`, becomes, with the spectrum
318 rescaled into :math:`(-1, 1)`, :math:`\tilde H = (H-b)/a`,
320 .. math::
322 \sigma^{B}_{\alpha\beta}(\mu) = \frac{4e^2\hbar}{\pi\Omega a^2}\int_{-1}^{1}
323 d\tilde\varepsilon\,\frac{f(\tilde\varepsilon)}{(1-\tilde\varepsilon^2)^2}
324 \sum_{m,n}\Gamma_{nm}(\tilde\varepsilon)\,\mu^{\alpha\beta}_{nm}\, ,
326 .. math::
328 \mu^{\alpha\beta}_{mn} = \frac{g_mg_n}{(1+\delta_{m0})(1+\delta_{n0})}
329 \mathrm{Tr}\left[v_\alpha T_m(\tilde H)\,v_\beta T_n(\tilde H)\right]\, ,
331 .. math::
333 \Gamma_{mn}(\tilde\varepsilon) = (\tilde\varepsilon - in\sqrt{1-\tilde\varepsilon^2})\,
334 e^{in\theta}\,T_m(\tilde\varepsilon)
335 + (\tilde\varepsilon + im\sqrt{1-\tilde\varepsilon^2})\,e^{-im\theta}\,T_n(\tilde\varepsilon)\, ,
337 :math:`\theta = \arccos\tilde\varepsilon`, :math:`g_m` the kernel. The
338 moments :math:`\mu_{mn}` (the expensive part: :math:`5MR` sparse
339 products, :math:`O(MN)` memory) serve every Fermi level and temperature.
340 The energy integral runs over :math:`|\tilde\varepsilon|\le 0.995`
341 (the spectrum lies within :math:`\pm0.99`): at the ends,
342 :math:`(1-\tilde\varepsilon^2)^{-2}` would only amplify rounding errors.
344 **Hall part.** The value is computed from the antisymmetric moments
345 :math:`(\mu^{xy}_{mn} - \mu^{yx}_{mn})/2`, i.e. it is the Hall
346 conductivity proper, :math:`(\sigma^B_{xy} - \sigma^B_{yx})/2`, the part
347 that the Berry curvature gives (*KSpace.hall_conductivity*). The
348 symmetric part :math:`(\sigma^B_{xy} + \sigma^B_{yx})/2` is a
349 Fermi-surface term, like :math:`\sigma_{xx}` (see *conductivity*): it
350 vanishes by symmetry in many lattices (e.g. with a :math:`C_3` axis),
351 but in a clean metal of lower symmetry it grows as the kernel's
352 broadening :math:`\pi a/M` shrinks (the Drude weight), and with random
353 vectors it is the main source of noise. Leaving it out removes that
354 noise too. In a clean metal the random-vector estimate is still noisy
355 inside the bands (increase *n_random*, or the size of the sample);
356 disorder, which makes the Fermi-surface terms finite, and gaps, where
357 they vanish, are much quieter.
359 **Sign and units.** The value returned is :math:`-(\sigma^B_{xy}-\sigma^B_{yx})/2` in units
360 of :math:`e^2/h` (:math:`\hbar = e = 1`), the sign of
361 *KSpace.hall_conductivity*: :math:`C` in the gap of a Chern insulator
362 of Chern number :math:`C` (*KSpace.chern_number*). Bastin's
363 :math:`\sigma^B_{xy}` is the Ohm's-law tensor, :math:`j_x = \sigma^B_{xy}E_y`
364 (see *KSpace.hall_conductivity*).
366 **Sample.** The trace must run over a sample without edges -- a torus
367 (*KSpace.finite_ham* and *KSpace.finite_velocity* with
368 ``periodic=True``): over a whole open flake, the edge currents cancel
369 the bulk and the trace vanishes in a gap.
371 :param ham: Square sparse Hermitian matrix.
372 :param vx: Square sparse (or dense) Hermitian matrix, like *ham*. Velocity
373 :math:`v_x = i[H, x]` (see *KSpace.finite_velocity*). For a spin
374 Hall conductivity, pass the spin current :math:`\{s, v_x\}/2`
375 instead (units :math:`e/2\pi`, see *KSpace.spin_hall_conductivity*).
376 :param vy: Square sparse (or dense) Hermitian matrix. Velocity :math:`v_y`.
377 :param n_moments: Positive integer. Default value 256. Number of
378 Chebyshev moments :math:`M` (energy resolution about :math:`\pi a/M`).
379 :param n_random: Positive integer, or None. Default value 10. See *dos*.
380 :param e_grid: Real array. Default value None: 1001 Chebyshev nodes over
381 the spectrum. Fermi energies :math:`\mu`.
382 :param kernel: String. Default value 'jackson'. See *dos*.
383 :param seed: See *dos*.
384 :param bounds: See *dos*.
385 :param area: Positive real. Default value None (the number of orbitals).
386 Area :math:`\Omega` of the sample (e.g. number of cells times the
387 cell area).
388 :param temperature: Positive real or zero. Default value 0. In energy units.
390 :returns:
391 * **e_grid** -- Real ndarray.
392 * **sigma** -- Real ndarray, same shape, in units of :math:`e^2/h`.
393 '''
394 ham_s, e_grid, _, a, b = _setup(ham, n_moments, kernel, bounds, e_grid)
395 n = ham_s.shape[0]
396 error_handling.velocity(vx, n, 'vx')
397 error_handling.velocity(vy, n, 'vy')
398 vx = sparse.csr_matrix(vx, dtype='c16')
399 vy = sparse.csr_matrix(vy, dtype='c16')
400 error_handling.hermitian(vx)
401 error_handling.hermitian(vy)
402 if area is None:
403 area = float(n)
404 error_handling.positive_real(area, 'area')
405 error_handling.positive_real_zero(temperature, 'temperature')
406 g = jackson_kernel(n_moments) if kernel == 'jackson' else lorentz_kernel(n_moments)
407 w = g / np.where(np.arange(n_moments) == 0, 2., 1.)
408 mu = np.zeros((n_moments, n_moments), 'c16')
409 for r in _vectors(n, n_random, seed).T:
410 # <r| v_x T_m = (T_m v_x |r>)^dagger and v_y T_n |r>, for mu^xy;
411 # the same with x and y swapped, for mu^yx
412 t_r = list(_chebyshev(ham_s, r, n_moments))
413 left_x = np.array(list(_chebyshev(ham_s, vx @ r, n_moments)))
414 left_y = np.array(list(_chebyshev(ham_s, vy @ r, n_moments)))
415 right_y = np.array([vy @ t for t in t_r])
416 right_x = np.array([vx @ t for t in t_r])
417 mu += (left_x.conj() @ right_y.T - left_y.conj() @ right_x.T) / 2
418 if n_random is not None:
419 mu = mu / n_random
420 mu = w[:, None] * mu * w[None, :]
421 # the energy integrand on a grid uniform in theta = arccos(eps)
422 n_int = max(2001, 8 * n_moments + 1)
423 theta = np.linspace(np.arccos(-0.995), np.arccos(0.995), n_int)
424 eps = np.cos(theta)
425 root = np.sin(theta)
426 m = np.arange(n_moments)[:, None]
427 cheb = np.cos(m * theta[None, :])
428 plus = (eps[None, :] - 1j * m * root[None, :]) * np.exp(1j * m * theta[None, :])
429 minus = (eps[None, :] + 1j * m * root[None, :]) * np.exp(-1j * m * theta[None, :])
430 # sum_{n,m} Gamma_nm mu_nm, Gamma_nm = plus_m T_n + minus_n T_m
431 gamma_mu = np.sum(cheb * (mu @ plus), axis=0) + np.sum(minus * (mu @ cheb), axis=0)
432 integrand = (gamma_mu / (1. - eps ** 2) ** 2).real
433 prefactor = -2 * PI * 4. / (PI * area * a ** 2) # e^2/hbar = 2 pi e^2/h, TKNN sign
434 energies = a * eps + b
435 if temperature == 0:
436 cum = np.concatenate([[0.], np.cumsum(0.5 * (integrand[1:] + integrand[:-1]) * np.diff(eps))])
437 sigma = prefactor * np.interp(e_grid, energies, cum)
438 else:
439 f = np.array([occupation.fermi_dirac(energies, float(e), temperature) for e in e_grid])
440 sigma = prefactor * trapezoid(f * integrand[None, :], eps, axis=1)
441 return e_grid, sigma