Coverage for tbkit/kpm.py: 100%
142 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +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 Parameters
47 ----------
48 ham : Square sparse (or dense) Hermitian matrix
50 Returns
51 -------
52 e_min, e_max : Real numbers
53 '''
54 ham = sparse.csr_matrix(ham)
55 diag = ham.diagonal().real
56 radii = np.asarray(abs(ham).sum(axis=1)).ravel() - np.abs(ham.diagonal())
57 return float(np.min(diag - radii)), float(np.max(diag + radii))
60def jackson_kernel(n_moments: int) -> NDArray[np.float64]:
61 r'''
62 Get the Jackson kernel
64 .. math::
66 g_m = \frac{(M-m+1)\cos\frac{\pi m}{M+1}
67 + \sin\frac{\pi m}{M+1}\cot\frac{\pi}{M+1}}{M+1}\, ,
69 the optimal choice for spectral densities (positive, and resolution
70 about :math:`\pi/M` in the scaled energy).
72 Parameters
73 ----------
74 n_moments : Positive integer :math:`M`
76 Returns
77 -------
78 g : Real ndarray, shape (n_moments,)
79 '''
80 error_handling.positive_int(n_moments, 'n_moments')
81 m = np.arange(n_moments)
82 q = PI / (n_moments + 1)
83 return ((n_moments - m + 1) * np.cos(q * m) + np.sin(q * m) / np.tan(q)) / (n_moments + 1)
86def lorentz_kernel(n_moments: int, lam: float = 4.) -> NDArray[np.float64]:
87 r'''
88 Get the Lorentz kernel :math:`g_m = \sinh[\lambda(1 - m/M)]/\sinh\lambda`,
89 which gives Lorentzian broadening (the better choice for Green's
90 functions).
92 Parameters
93 ----------
94 n_moments : Positive integer :math:`M`
95 lam : Positive real, default 4
97 Returns
98 -------
99 g : Real ndarray, shape (n_moments,)
100 '''
101 error_handling.positive_int(n_moments, 'n_moments')
102 error_handling.positive_real(lam, 'lam')
103 return np.sinh(lam * (1. - np.arange(n_moments) / n_moments)) / np.sinh(lam)
106def _setup(
107 ham, n_moments, kernel, bounds, e_grid, n_grid=1001,
108) -> tuple[sparse.csr_matrix, NDArray[np.float64], NDArray[np.float64], float, float]:
109 '''
110 Private function. Validate, rescale the Hamiltonian, and build the
111 coefficients c_m(E) (shape (n_moments, len(e_grid))).
112 '''
113 ham = sparse.csr_matrix(ham, dtype='c16')
114 error_handling.square_matrix(ham, 'ham')
115 error_handling.hermitian(ham)
116 error_handling.positive_int(n_moments, 'n_moments')
117 error_handling.kpm_kernel(kernel)
118 if bounds is None:
119 bounds = spectral_bounds(ham)
120 error_handling.lims(bounds)
121 e_min, e_max = bounds
122 a = (e_max - e_min) / (2. - 0.02) # keep the scaled spectrum inside (-0.99, 0.99)
123 b = (e_max + e_min) / 2.
124 if e_grid is None:
125 e_grid = b + a * np.cos(PI * (np.arange(n_grid) + 0.5) / n_grid)[::-1]
126 e_grid = np.atleast_1d(np.asarray(e_grid, dtype='f8'))
127 error_handling.ndarray_empty(e_grid, 'e_grid')
128 eps = (e_grid - b) / a
129 error_handling.kpm_grid(eps)
130 g = jackson_kernel(n_moments) if kernel == 'jackson' else lorentz_kernel(n_moments)
131 tm = np.cos(np.arange(n_moments)[:, None] * np.arccos(eps)[None, :])
132 weight = np.where(np.arange(n_moments) == 0, 1., 2.)
133 coef = (weight * g)[:, None] * tm / (a * PI * np.sqrt(1. - eps ** 2))[None, :]
134 ham_s = (ham - b * sparse.identity(ham.shape[0], format='csr')) / a
135 return ham_s, e_grid, coef, a, b
138def _vectors(n: int, n_random: int | None, seed) -> NDArray[np.complex128]:
139 '''
140 Private function. Random-phase vectors (columns, |r_i| = 1), or, if
141 n_random is None, the full basis (the exact trace).
142 '''
143 if n_random is None:
144 return np.eye(n, dtype='c16')
145 error_handling.positive_int(n_random, 'n_random')
146 rng = np.random.default_rng(seed)
147 return np.exp(2j * PI * rng.random((n, n_random)))
150def _chebyshev(ham_s, vec: NDArray[np.complex128], n_moments: int):
151 '''
152 Private function. Yield T_m(H~) vec for m < n_moments, keeping only two
153 vectors in memory.
154 '''
155 prev, cur = None, vec
156 for m in range(n_moments):
157 yield cur
158 prev, cur = cur, (ham_s @ cur if m == 0 else 2 * (ham_s @ cur) - prev)
161def dos(
162 ham, n_moments: int = 256, n_random: int | None = 10, e_grid: ArrayLike | None = None,
163 kernel: str = 'jackson', seed=None, bounds: tuple[float, float] | None = None,
164) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
165 r'''
166 Get the density of states :math:`\rho(E) = \mathrm{Tr}\,\delta(E - H)`
167 (normalized, like *tbkit.dos.density_of_states*, to integrate to the
168 number of states).
170 Parameters
171 ----------
172 ham : Square sparse Hermitian matrix (e.g. *System.ham*)
173 n_moments : Positive integer, default 256
174 Number of
175 Chebyshev moments (the energy resolution is about
176 :math:`\pi a/M`, :math:`2a` the spectral width).
177 n_random : Positive integer, or None, default 10
178 Number of
179 random vectors for the trace; None for the exact trace (N vectors:
180 small systems only).
181 e_grid : Real array
182 Default value None: 1001 Chebyshev nodes over
183 the spectrum. Energies, inside the spectral bounds.
184 kernel : String, default 'jackson'
185 'jackson' or 'lorentz'.
186 seed
187 Default value None. Seed of the random vectors.
188 bounds : Pair of reals
189 Default value None (Gershgorin bounds, see
190 *spectral_bounds*). Bounds enclosing the whole spectrum.
192 Returns
193 -------
194 e_grid : Real ndarray
195 rho : Real ndarray, same shape
196 '''
197 ham_s, e_grid, coef, _, _ = _setup(ham, n_moments, kernel, bounds, e_grid)
198 vec = _vectors(ham_s.shape[0], n_random, seed)
199 mu = np.array([np.vdot(vec, t).real for t in _chebyshev(ham_s, vec, n_moments)])
200 if n_random is not None:
201 mu = mu / n_random # E<r|A|r> = Tr A for random phases
202 return e_grid, mu @ coef
205def dos_from_levels(
206 levels: ArrayLike, n_moments: int, bounds: tuple[float, float], e_grid: ArrayLike,
207 kernel: str = 'jackson',
208) -> NDArray[np.float64]:
209 r'''
210 Get the KPM density of states of a *known* spectrum: every level
211 broadened by exactly the truncated, kernel-damped Chebyshev series that
212 *dos* uses. It is the reference *dos* converges to (up to the random
213 trace), e.g. for the bands of a periodic model on a k-mesh.
215 Parameters
216 ----------
217 levels : Real array
218 Energy levels.
219 n_moments : Positive integer
220 Number of Chebyshev moments.
221 bounds : Pair of reals
222 Spectral bounds (those of *dos*, e.g.
223 *spectral_bounds* of the Hamiltonian it is compared with).
224 e_grid : Real array
225 Energies.
226 kernel : String, default 'jackson'
228 Returns
229 -------
230 rho
231 Real ndarray, same shape as *e_grid* (integrating to the number of levels).
232 '''
233 levels = np.asarray(levels, dtype='f8').ravel()
234 error_handling.ndarray_empty(levels, 'levels')
235 _, e_grid, coef, a, b = _setup(sparse.identity(1, format='csr'), n_moments, kernel, bounds, e_grid)
236 error_handling.kpm_grid((levels - b) / a)
237 tm = np.cos(np.arange(n_moments)[:, None] * np.arccos((levels - b) / a)[None, :])
238 return coef.T @ tm.sum(axis=1)
241def ldos(
242 ham, site: int, n_moments: int = 256, e_grid: ArrayLike | None = None,
243 kernel: str = 'jackson', bounds: tuple[float, float] | None = None,
244) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
245 r'''
246 Get the local density of states on one site (or orbital),
247 :math:`\rho_i(E) = \langle i|\delta(E-H)|i\rangle`, exactly (no random vectors).
249 Parameters
250 ----------
251 ham : Square sparse Hermitian matrix
252 site : Integer
253 Site (row) index.
254 n_moments : Positive integer, default 256
255 e_grid
256 See *dos*.
257 kernel
258 See *dos*.
259 bounds
260 See *dos*.
262 Returns
263 -------
264 e_grid : Real ndarray
265 rho : Real ndarray, same shape
266 '''
267 ham_s, e_grid, coef, _, _ = _setup(ham, n_moments, kernel, bounds, e_grid)
268 error_handling.site_index(site, ham_s.shape[0])
269 vec = np.zeros(ham_s.shape[0], 'c16')
270 vec[site] = 1.
271 mu = np.array([t[site].real for t in _chebyshev(ham_s, vec, n_moments)])
272 return e_grid, mu @ coef
275def conductivity(
276 ham, x: ArrayLike, n_moments: int = 128, n_random: int | None = 10,
277 e_grid: ArrayLike | None = None, kernel: str = 'jackson', seed=None,
278 bounds: tuple[float, float] | None = None, area: float | None = None,
279) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
280 r'''
281 Get the zero-temperature DC conductivity along :math:`x` from the
282 Kubo-Greenwood formula,
284 .. math::
286 \sigma_{xx}(E) = \frac{\pi}{\Omega}\,
287 \mathrm{Tr}\left[v_x\,\delta(E-H)\,v_x\,\delta(E-H)\right]\, ,\qquad
288 v_x = i[H, x]\, ,
290 in units of :math:`e^2/\hbar` (:math:`\hbar = e = 1`), by the
291 two-dimensional Chebyshev expansion of Weisse et al. (the moments
292 :math:`\mu_{mn} = \mathrm{Tr}[v_x T_m(\tilde H) v_x T_n(\tilde H)]`,
293 :math:`O(M N)` memory, the DOS needing only :math:`O(N)`).
295 Parameters
296 ----------
297 ham : Square sparse Hermitian matrix
298 x : Real array, one coordinate per site (e.g. *lat.coor['x']*)
299 n_moments : Positive integer, default 128
300 n_random : Positive integer, or None, default 10
301 See *dos*.
302 e_grid
303 See *dos*.
304 kernel
305 See *dos*.
306 seed
307 See *dos*.
308 bounds
309 See *dos*.
310 area : Positive real
311 Default value None (the number of sites):
312 the area (volume) :math:`\Omega` of the sample.
314 Returns
315 -------
316 e_grid : Real ndarray
317 sigma : Real ndarray, same shape
318 '''
319 ham_s, e_grid, coef, a, b = _setup(ham, n_moments, kernel, bounds, e_grid)
320 n = ham_s.shape[0]
321 x = np.asarray(x, dtype='f8')
322 error_handling.ndarray(x, 'x', n)
323 if area is None:
324 area = float(n)
325 error_handling.positive_real(area, 'area')
326 # v = i [H, x]: v_ij = i H_ij (x_j - x_i); the scale a comes back below
327 ham_c = sparse.csr_matrix(ham, dtype='c16')
328 vel = 1j * (ham_c @ sparse.diags(x) - sparse.diags(x) @ ham_c)
329 vec = _vectors(n, n_random, seed)
330 mu = np.zeros((n_moments, n_moments), 'c16')
331 for r in vec.T:
332 left = np.array(list(_chebyshev(ham_s, vel @ r, n_moments))) # T_m v |r>
333 right = np.array([vel @ t for t in _chebyshev(ham_s, r, n_moments)]) # v T_n |r>
334 mu += left.conj() @ right.T
335 if n_random is not None:
336 mu = mu / n_random
337 sigma = PI / area * np.einsum('me,mn,ne->e', coef, mu, coef).real
338 return e_grid, sigma
341def hall_conductivity(
342 ham, vx, vy, n_moments: int = 256, n_random: int | None = 10,
343 e_grid: ArrayLike | None = None, kernel: str = 'jackson', seed=None,
344 bounds: tuple[float, float] | None = None, area: float | None = None,
345 temperature: float = 0.,
346) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
347 r'''
348 Get the Hall conductivity :math:`\sigma_{xy}(\mu)` of a large sample, at
349 any Fermi level and temperature, from the Kubo-Bastin formula expanded
350 in Chebyshev polynomials (Garcia, Covaci and Rappoport, Phys. Rev. Lett.
351 114, 116602 (2015), Eqs. (1)-(5)). Bastin's formula,
353 .. math::
355 \sigma^{B}_{\alpha\beta}(\mu) = \frac{ie^2\hbar}{\Omega}\int d\varepsilon\,f(\varepsilon)\,
356 \mathrm{Tr}\left[v_\alpha\,\delta(\varepsilon-H)\,v_\beta\,\frac{dG^+}{d\varepsilon}
357 - v_\alpha\,\frac{dG^-}{d\varepsilon}\,v_\beta\,\delta(\varepsilon-H)\right]\, ,
359 :math:`G^\pm = (\varepsilon - H \pm i0)^{-1}`, becomes, with the spectrum
360 rescaled into :math:`(-1, 1)`, :math:`\tilde H = (H-b)/a`,
362 .. math::
364 \sigma^{B}_{\alpha\beta}(\mu) = \frac{4e^2\hbar}{\pi\Omega a^2}\int_{-1}^{1}
365 d\tilde\varepsilon\,\frac{f(\tilde\varepsilon)}{(1-\tilde\varepsilon^2)^2}
366 \sum_{m,n}\Gamma_{nm}(\tilde\varepsilon)\,\mu^{\alpha\beta}_{nm}\, ,
368 .. math::
370 \mu^{\alpha\beta}_{mn} = \frac{g_mg_n}{(1+\delta_{m0})(1+\delta_{n0})}
371 \mathrm{Tr}\left[v_\alpha T_m(\tilde H)\,v_\beta T_n(\tilde H)\right]\, ,
373 .. math::
375 \Gamma_{mn}(\tilde\varepsilon) = (\tilde\varepsilon - in\sqrt{1-\tilde\varepsilon^2})\,
376 e^{in\theta}\,T_m(\tilde\varepsilon)
377 + (\tilde\varepsilon + im\sqrt{1-\tilde\varepsilon^2})\,e^{-im\theta}\,T_n(\tilde\varepsilon)\, ,
379 :math:`\theta = \arccos\tilde\varepsilon`, :math:`g_m` the kernel. The
380 moments :math:`\mu_{mn}` (the expensive part: :math:`5MR` sparse
381 products, :math:`O(MN)` memory) serve every Fermi level and temperature.
382 The energy integral runs over :math:`|\tilde\varepsilon|\le 0.995`
383 (the spectrum lies within :math:`\pm0.99`): at the ends,
384 :math:`(1-\tilde\varepsilon^2)^{-2}` would only amplify rounding errors.
386 **Hall part.** The value is computed from the antisymmetric moments
387 :math:`(\mu^{xy}_{mn} - \mu^{yx}_{mn})/2`, i.e. it is the Hall
388 conductivity proper, :math:`(\sigma^B_{xy} - \sigma^B_{yx})/2`, the part
389 that the Berry curvature gives (*KSpace.hall_conductivity*). The
390 symmetric part :math:`(\sigma^B_{xy} + \sigma^B_{yx})/2` is a
391 Fermi-surface term, like :math:`\sigma_{xx}` (see *conductivity*): it
392 vanishes by symmetry in many lattices (e.g. with a :math:`C_3` axis),
393 but in a clean metal of lower symmetry it grows as the kernel's
394 broadening :math:`\pi a/M` shrinks (the Drude weight), and with random
395 vectors it is the main source of noise. Leaving it out removes that
396 noise too. In a clean metal the random-vector estimate is still noisy
397 inside the bands (increase *n_random*, or the size of the sample);
398 disorder, which makes the Fermi-surface terms finite, and gaps, where
399 they vanish, are much quieter.
401 **Sign and units.** The value returned is :math:`-(\sigma^B_{xy}-\sigma^B_{yx})/2` in units
402 of :math:`e^2/h` (:math:`\hbar = e = 1`), the sign of
403 *KSpace.hall_conductivity*: :math:`C` in the gap of a Chern insulator
404 of Chern number :math:`C` (*KSpace.chern_number*). Bastin's
405 :math:`\sigma^B_{xy}` is the Ohm's-law tensor, :math:`j_x = \sigma^B_{xy}E_y`
406 (see *KSpace.hall_conductivity*).
408 **Sample.** The trace must run over a sample without edges -- a torus
409 (*KSpace.finite_ham* and *KSpace.finite_velocity* with
410 ``periodic=True``): over a whole open flake, the edge currents cancel
411 the bulk and the trace vanishes in a gap.
413 Parameters
414 ----------
415 ham : Square sparse Hermitian matrix
416 vx : Square sparse (or dense) Hermitian matrix, like *ham*
417 Velocity
418 :math:`v_x = i[H, x]` (see *KSpace.finite_velocity*). For a spin
419 Hall conductivity, pass the spin current :math:`\{s, v_x\}/2`
420 instead (units :math:`e/2\pi`, see *KSpace.spin_hall_conductivity*).
421 vy : Square sparse (or dense) Hermitian matrix
422 Velocity :math:`v_y`.
423 n_moments : Positive integer, default 256
424 Number of
425 Chebyshev moments :math:`M` (energy resolution about :math:`\pi a/M`).
426 n_random : Positive integer, or None, default 10
427 See *dos*.
428 e_grid : Real array
429 Default value None: 1001 Chebyshev nodes over
430 the spectrum. Fermi energies :math:`\mu`.
431 kernel : String, default 'jackson'
432 See *dos*.
433 seed
434 See *dos*.
435 bounds
436 See *dos*.
437 area : Positive real
438 Default value None (the number of orbitals).
439 Area :math:`\Omega` of the sample (e.g. number of cells times the
440 cell area).
441 temperature : Positive real or zero, default 0
442 In energy units.
444 Returns
445 -------
446 e_grid : Real ndarray
447 sigma : Real ndarray, same shape, in units of :math:`e^2/h`
448 '''
449 ham_s, e_grid, _, a, b = _setup(ham, n_moments, kernel, bounds, e_grid)
450 n = ham_s.shape[0]
451 error_handling.velocity(vx, n, 'vx')
452 error_handling.velocity(vy, n, 'vy')
453 vx = sparse.csr_matrix(vx, dtype='c16')
454 vy = sparse.csr_matrix(vy, dtype='c16')
455 error_handling.hermitian(vx)
456 error_handling.hermitian(vy)
457 if area is None:
458 area = float(n)
459 error_handling.positive_real(area, 'area')
460 error_handling.positive_real_zero(temperature, 'temperature')
461 g = jackson_kernel(n_moments) if kernel == 'jackson' else lorentz_kernel(n_moments)
462 w = g / np.where(np.arange(n_moments) == 0, 2., 1.)
463 mu = np.zeros((n_moments, n_moments), 'c16')
464 for r in _vectors(n, n_random, seed).T:
465 # <r| v_x T_m = (T_m v_x |r>)^dagger and v_y T_n |r>, for mu^xy;
466 # the same with x and y swapped, for mu^yx
467 t_r = list(_chebyshev(ham_s, r, n_moments))
468 left_x = np.array(list(_chebyshev(ham_s, vx @ r, n_moments)))
469 left_y = np.array(list(_chebyshev(ham_s, vy @ r, n_moments)))
470 right_y = np.array([vy @ t for t in t_r])
471 right_x = np.array([vx @ t for t in t_r])
472 mu += (left_x.conj() @ right_y.T - left_y.conj() @ right_x.T) / 2
473 if n_random is not None:
474 mu = mu / n_random
475 mu = w[:, None] * mu * w[None, :]
476 # the energy integrand on a grid uniform in theta = arccos(eps)
477 n_int = max(2001, 8 * n_moments + 1)
478 theta = np.linspace(np.arccos(-0.995), np.arccos(0.995), n_int)
479 eps = np.cos(theta)
480 root = np.sin(theta)
481 m = np.arange(n_moments)[:, None]
482 cheb = np.cos(m * theta[None, :])
483 plus = (eps[None, :] - 1j * m * root[None, :]) * np.exp(1j * m * theta[None, :])
484 minus = (eps[None, :] + 1j * m * root[None, :]) * np.exp(-1j * m * theta[None, :])
485 # sum_{n,m} Gamma_nm mu_nm, Gamma_nm = plus_m T_n + minus_n T_m
486 gamma_mu = np.sum(cheb * (mu @ plus), axis=0) + np.sum(minus * (mu @ cheb), axis=0)
487 integrand = (gamma_mu / (1. - eps ** 2) ** 2).real
488 prefactor = -2 * PI * 4. / (PI * area * a ** 2) # e^2/hbar = 2 pi e^2/h, TKNN sign
489 energies = a * eps + b
490 if temperature == 0:
491 cum = np.concatenate([[0.], np.cumsum(0.5 * (integrand[1:] + integrand[:-1]) * np.diff(eps))])
492 sigma = prefactor * np.interp(e_grid, energies, cum)
493 else:
494 f = np.array([occupation.fermi_dirac(energies, float(e), temperature) for e in e_grid])
495 sigma = prefactor * trapezoid(f * integrand[None, :], eps, axis=1)
496 return e_grid, sigma