Coverage for tbkit/optics.py: 100%
149 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +0000
1r"""
2Linear and nonlinear optical response of 2D Bloch models (:class:`~tbkit.kspace.KSpace`)
3over a uniform Brillouin-zone mesh.
5* :func:`optical_conductivity` -- the Kubo-Greenwood conductivity
6 :math:`\sigma_{ab}(\omega)`, interband and intraband (Drude), with a
7 broadening :math:`\eta`.
8* :func:`joint_dos` -- the joint density of states of vertical transitions.
9* :func:`berry_curvature_dipole` -- the Berry curvature dipole of the
10 nonlinear Hall effect (Sodemann and Fu 2015).
11* :func:`shift_current` -- the shift-current conductivity
12 :math:`\sigma^{abc}(0;\omega,-\omega)` (Sipe and Shkrebtii 2000).
14Every function works in units :math:`\hbar = e = 1`, energies (and
15:math:`\hbar\omega`) in the units of the hoppings, lengths in those of
16*prim_vec*, and counts each orbital once: a spinless model gives the
17response *per spin*.
19The velocities are the exact :math:`\partial H/\partial k_a` of
20*KSpace._bloch_derivatives* along the Cartesian axes :math:`x, y` (for a 2D
21lattice embedded in 3D space, the axes of *KSpace.k_basis*). With
22``positions=True`` (the default) the bond vectors include the orbital
23positions, :math:`\mathbf{d} = \mathbf{R} + \boldsymbol\tau_j -
24\boldsymbol\tau_i`, so that the velocity is the physical :math:`i[H,
25\mathbf{r}]` (see *KSpace.hall_conductivity*).
27Examples
28--------
29Example usage::
31 from tbkit.optics import optical_conductivity
32 omega = np.linspace(0.05, 3., 60)
33 sigma = optical_conductivity(graphene, omega, e_fermi=0., nk=300)
34 # Re sigma ~ pi/4 e^2/h per spin for hbar omega << t
35"""
36from __future__ import annotations
38from typing import Any, Iterator
40import numpy as np
41from numpy.typing import ArrayLike, NDArray
43import tbkit.error_handling as error_handling
44import tbkit.occupation as occupation
45from tbkit.kspace import KSpace
48_AXES = {'x': 0, 'y': 1}
49_DEGENERATE = 1e-9
52def _mesh(ks: KSpace, nk: int | tuple[int, int], positions: bool) -> NDArray[np.float64]:
53 '''
54 Private function. Validate the model (a static, Hermitian, 2D KSpace)
55 and return the k-points of a uniform Brillouin-zone mesh.
56 '''
57 error_handling.kspace(ks, KSpace)
58 ks._check_static()
59 error_handling.dim_exact(ks.dim, 2)
60 error_handling.hermitian_model(ks.is_hermitian())
61 error_handling.boolean(positions, 'positions')
62 return ks.mesh_grid(nk)[1]
65def cell_area(ks: KSpace) -> float:
66 r'''
67 Get the area :math:`A_c = |\mathbf{a}_1\times\mathbf{a}_2|` of the unit
68 cell of a 2D model (in the units of *prim_vec*), which turns a
69 Brillouin-zone average into an integral:
70 :math:`\int\frac{d^2k}{(2\pi)^2}g = \frac{1}{N_kA_c}\sum_{\mathbf{k}}g`.
72 Parameters
73 ----------
74 ks : **KSpace** instance, 2D
76 Returns
77 -------
78 area : Positive real
79 '''
80 error_handling.kspace(ks, KSpace)
81 error_handling.dim_exact(ks.dim, 2)
82 a1, a2 = (np.zeros(3) for _ in range(2))
83 a1[:ks.space_dim] = ks.lat.prim_vec[0]
84 a2[:ks.space_dim] = ks.lat.prim_vec[1]
85 return float(np.linalg.norm(np.cross(a1, a2)))
88def _bands(
89 ks: KSpace, kpts: NDArray[np.float64], positions: bool, second: bool = False,
90) -> Iterator[tuple]:
91 r'''
92 Private function. Over chunks of the k-points *kpts*: the band energies
93 (nk, n), the velocity matrices :math:`v^a_{nm} = \langle n|\partial_aH|m\rangle`
94 in the eigenbasis (2, nk, n, n) and, with *second*, the matrices
95 :math:`w^{ab}_{nm} = \langle n|\partial_a\partial_bH|m\rangle` (2, 2, nk, n, n).
96 '''
97 directions = [ks.k_basis[:, 0], ks.k_basis[:, 1]]
98 chunk = max(1, 200000 // ks.norb ** 2)
99 if second:
100 # d^2 H / dk_a dk_b is the first derivative of the Bloch sum of the
101 # hoppings t * i(d.u_a): the same _bloch_sum, with rescaled amplitudes.
102 hops = ks._hop
103 d = np.array([h[2] for h in hops], dtype='f8').reshape(len(hops), ks.space_dim)
104 if positions:
105 tau = ks.orbital_positions()
106 d = d + tau[[h[1] for h in hops]] - tau[[h[0] for h in hops]]
107 scaled = [[(h[0], h[1], h[2], h[3] * 1j * (dd @ u)) for h, dd in zip(hops, d)]
108 for u in directions]
109 for c0 in range(0, len(kpts), chunk):
110 kc = kpts[c0:c0 + chunk]
111 ham, dham = ks._bloch_derivatives(kc, directions, positions)
112 en, vec = np.linalg.eigh(ham)
113 vh = vec.conj().transpose(0, 2, 1)
114 vel = np.array([vh @ dh @ vec for dh in dham])
115 if not second:
116 yield en, vel
117 continue
118 k_cart = kc @ ks.k_basis.T
119 w = np.array([ks._bloch_sum(h, k_cart, directions, positions)[1] for h in scaled])
120 yield en, vel, vh[None, None] @ w @ vec[None, None]
123def _as_output(values: NDArray, like: ArrayLike) -> Any:
124 '''
125 Private function. A scalar for a scalar input, else shaped like it.
126 '''
127 if np.ndim(like) == 0:
128 return values[0]
129 return values.reshape(np.shape(like) + values.shape[1:])
132def optical_conductivity(
133 ks: KSpace, omega: float | ArrayLike, e_fermi: float = 0., temperature: float = 0.,
134 eta: float = 0.05, nk: int | tuple[int, int] = 100, component: str = 'xx',
135 positions: bool = True,
136) -> complex | NDArray[np.complex128]:
137 r'''
138 Get the optical conductivity :math:`\sigma_{ab}(\omega)` of a 2D model
139 from the Kubo-Greenwood formula over a uniform Brillouin-zone mesh,
141 .. math::
143 \sigma_{ab}(\omega) = \frac{i e^2}{\hbar}\int\frac{d^2k}{(2\pi)^2}
144 \sum_{n,m}\frac{f_n - f_m}{E_m - E_n}\,
145 \frac{v^a_{nm}\,v^b_{mn}}{\hbar\omega + E_n - E_m + i\eta}\, ,
147 with :math:`v^a_{nm} = \langle n|\partial H/\partial k_a|m\rangle`,
148 :math:`f_n` the Fermi-Dirac occupations at *e_fermi* and *temperature*,
149 and :math:`\eta` a phenomenological broadening (the inverse lifetime).
150 Pairs of degenerate states, :math:`n = m` included, enter through the
151 limit :math:`(f_n-f_m)/(E_m-E_n)\to-\partial f/\partial E`: the
152 intraband (Drude) term :math:`\frac{i}{\hbar\omega+i\eta}\sum_n
153 (-f'_n)v^a_{nn}v^b_{nn}`. It samples the Fermi surface only through the
154 thermal window, so it needs a *temperature* larger than the level
155 spacing of the mesh; at ``temperature=0`` only the interband part is
156 computed.
158 It is the Ohm's-law tensor, :math:`j_a = \sigma_{ab}E_b`. In a gap its
159 static Hall part is the TKNN conductance:
160 :math:`\sigma_{yx}(0) = -\sigma_{xy}(0)` equals
161 *KSpace.hall_conductivity* (which returns :math:`\sigma^{\mathrm{Ohm}}_{yx}`,
162 see its docstring). The absorptive part obeys the f-sum rule
163 :math:`\int_0^\infty\mathrm{Re}\,\sigma_{aa}\,d\omega =
164 \frac{\pi e^2}{2\hbar^2}\int\frac{d^2k}{(2\pi)^2}\sum_n f_n\langle n|\partial_a^2H|n\rangle`.
165 For graphene at neutrality, :math:`\mathrm{Re}\,\sigma_{xx} \to
166 \frac{e^2}{16\hbar}` per spin and valley for :math:`\eta \ll \hbar\omega
167 \ll t`: the universal :math:`\sigma_0 = e^2/4\hbar` once spin and
168 valleys are counted, :math:`\pi/4` in units of :math:`e^2/h` for a
169 spinless model.
171 Parameters
172 ----------
173 ks : **KSpace** instance, 2D, Hermitian
174 omega : Real number or array of real numbers
175 Photon energies
176 :math:`\hbar\omega`.
177 e_fermi : Real number, default 0
178 Fermi energy.
179 temperature : Positive real or zero
180 Default value 0 (:math:`k_B = 1`).
181 eta : Positive real, default 0.05
182 Broadening, larger than
183 the spacing of the transition energies on the mesh.
184 nk : Positive integer, or tuple of 2, default 100
185 k-mesh.
186 component : 'xx', 'xy', 'yx' or 'yy', default 'xx'
187 positions : Boolean, default True
188 Orbital positions in the
189 bond vectors (see *KSpace.hall_conductivity*).
191 Returns
192 -------
193 sigma
194 Complex number (or ndarray shaped like *omega*), in
195 units of :math:`e^2/h`.
196 '''
197 kpts = _mesh(ks, nk, positions)
198 error_handling.frequencies(omega, 'omega')
199 error_handling.real_number(e_fermi, 'e_fermi')
200 error_handling.positive_real_zero(temperature, 'temperature')
201 error_handling.positive_real(eta, 'eta')
202 a, b = error_handling.tensor_component(component, 2)
203 w_arr = np.atleast_1d(np.asarray(omega, dtype='f8')).ravel()
204 total = np.zeros(len(w_arr), 'c16')
205 for en, vel in _bands(ks, kpts, positions):
206 f = occupation.fermi_dirac(en, e_fermi, temperature)
207 minus_df = f * (1. - f) / temperature if temperature > 0 else np.zeros_like(f)
208 de = en[:, None, :] - en[:, :, None] # E_m - E_n
209 degenerate = np.abs(de) < _DEGENERATE
210 ratio = np.where(degenerate, minus_df[:, :, None],
211 (f[:, :, None] - f[:, None, :]) / np.where(degenerate, 1., de))
212 weight = ratio * vel[a] * vel[b].transpose(0, 2, 1)
213 keep = weight != 0
214 wk, dk = weight[keep], np.where(degenerate, 0., de)[keep]
215 step = max(1, 2000000 // len(w_arr))
216 for c0 in range(0, len(wk), step):
217 total += (wk[None, c0:c0 + step]
218 / (w_arr[:, None] - dk[None, c0:c0 + step] + 1j * eta)).sum(axis=1)
219 sigma = 1j * 2 * np.pi * total / (len(kpts) * cell_area(ks))
220 return _as_output(sigma, omega)
223def joint_dos(
224 ks: KSpace, omega: float | ArrayLike, e_fermi: float = 0., temperature: float = 0.,
225 eta: float = 0.05, nk: int | tuple[int, int] = 100,
226) -> float | NDArray[np.float64]:
227 r'''
228 Get the joint density of states of the vertical (interband) transitions,
229 per unit cell,
231 .. math::
233 J(\omega) = \frac{1}{N_k}\sum_{\mathbf{k}}\sum_{n\neq m}
234 f_n\,(1 - f_m)\,\delta(\hbar\omega - E_m + E_n)\, ,
236 with the delta function a Gaussian of standard deviation :math:`\eta`.
237 Its integral over :math:`\omega` is the number of (occupied, empty)
238 pairs per k-point. For graphene at neutrality,
239 :math:`J = A_c\,\hbar\omega/(4\pi\hbar^2 v_F^2)` near the Dirac points
240 (both valleys, per spin).
242 Parameters
243 ----------
244 ks : **KSpace** instance, 2D, Hermitian
245 omega : Real number or array of real numbers
246 Photon energies.
247 e_fermi : Real number, default 0
248 temperature : Positive real or zero, default 0
249 eta : Positive real, default 0.05
250 Gaussian width.
251 nk : Positive integer, or tuple of 2, default 100
253 Returns
254 -------
255 jdos
256 Real number (or ndarray shaped like *omega*), in states
257 per unit cell and per unit energy.
258 '''
259 kpts = _mesh(ks, nk, True)
260 error_handling.frequencies(omega, 'omega')
261 error_handling.real_number(e_fermi, 'e_fermi')
262 error_handling.positive_real_zero(temperature, 'temperature')
263 error_handling.positive_real(eta, 'eta')
264 w_arr = np.atleast_1d(np.asarray(omega, dtype='f8')).ravel()
265 total = np.zeros(len(w_arr))
266 for en, _ in _bands(ks, kpts, False):
267 f = occupation.fermi_dirac(en, e_fermi, temperature)
268 weight = f[:, :, None] * (1. - f[:, None, :])
269 weight[:, np.arange(ks.norb), np.arange(ks.norb)] = 0.
270 keep = weight > 0
271 wk = weight[keep]
272 dk = (en[:, None, :] - en[:, :, None])[keep]
273 total += _gaussian_sum(w_arr, dk, wk, eta)
274 return _as_output(total / len(kpts), omega)
277def _gaussian_sum(
278 x: NDArray[np.float64], centers: NDArray[np.float64], weights: NDArray[np.float64], eta: float,
279) -> NDArray[np.float64]:
280 r'''
281 Private function. :math:`\sum_i w_i\,g(x - c_i)` with the normalized
282 Gaussian :math:`g` of standard deviation *eta*, in chunks of the centers.
283 '''
284 out = np.zeros(len(x))
285 step = max(1, 2000000 // len(x))
286 for c0 in range(0, len(centers), step):
287 u = (x[:, None] - centers[None, c0:c0 + step]) / eta
288 out += (weights[None, c0:c0 + step] * np.exp(-0.5 * u ** 2)).sum(axis=1)
289 return out / (np.sqrt(2 * np.pi) * eta)
292def berry_curvature_dipole(
293 ks: KSpace, e_fermi: float | ArrayLike, temperature: float,
294 nk: int | tuple[int, int] = 100, positions: bool = True,
295) -> NDArray[np.float64]:
296 r'''
297 Get the Berry curvature dipole of a 2D model (Sodemann and Fu, Phys.
298 Rev. Lett. 115, 216806 (2015)),
300 .. math::
302 D_a = \int\frac{d^2k}{(2\pi)^2}\sum_n f_n\,\partial_a\Omega_n
303 = \int\frac{d^2k}{(2\pi)^2}\sum_n
304 \left(-\frac{\partial f}{\partial E}\right)_{E_n}
305 v^a_{nn}\,\Omega_n\, ,
307 the first moment of the Berry curvature :math:`\Omega_n` (the
308 convention of *KSpace.berry_curvature* and *hall_conductivity*,
309 :math:`\Omega_n = -2\,\mathrm{Im}\sum_{m\neq n}v^x_{nm}v^y_{mn}/(E_n-E_m)^2`)
310 over the occupied states. It is a Fermi-surface property: the second
311 form, used here, samples it through the thermal window
312 :math:`-\partial f/\partial E`, so *temperature* must be positive and
313 larger than the level spacing of the mesh near :math:`E_F`. A finite
314 dipole needs broken inversion symmetry and a single mirror line at
315 most. In a DC field it drives the nonlinear Hall current
316 :math:`j_a = \chi_{abc}E_bE_c`,
317 :math:`\chi_{abc} = -\varepsilon_{adc}\frac{e^3\tau}{2\hbar^2(1+i\omega\tau)}D_{bd}`,
318 with :math:`D_{bd} = D_b` along :math:`d = z` in 2D.
320 For the tilted massive Dirac cone
321 :math:`H = t k_x + v(k_x\sigma_x + k_y\sigma_y) + m\sigma_z` doped into
322 the conduction band (:math:`\mu > |m|`), the first order in the tilt is
323 :math:`D_x = -\frac{3\,t\,m\,(\mu^2 - m^2)}{8\pi\mu^4}`,
324 :math:`D_y = 0`, from :math:`\Omega_c = -mv^2/2\varepsilon^3`
325 (Sodemann and Fu's model, in these conventions).
327 Parameters
328 ----------
329 ks : **KSpace** instance, 2D, Hermitian
330 e_fermi : Real number, or array of real numbers
331 Fermi energies.
332 temperature : Positive real
333 Temperature (:math:`k_B = 1`).
334 nk : Positive integer, or tuple of 2, default 100
335 positions : Boolean, default True
336 See *KSpace.hall_conductivity*.
338 Returns
339 -------
340 dipole
341 Real ndarray of shape (2,): :math:`(D_x, D_y)` (or
342 ``e_fermi.shape + (2,)``), in units of length (of *prim_vec*).
343 '''
344 kpts = _mesh(ks, nk, positions)
345 error_handling.fermi_energies(e_fermi)
346 error_handling.positive_real(temperature, 'temperature')
347 mus = np.atleast_1d(np.asarray(e_fermi, dtype='f8')).ravel()
348 total = np.zeros((len(mus), 2))
349 for en, vel in _bands(ks, kpts, positions):
350 de = en[:, :, None] - en[:, None, :]
351 degenerate = np.abs(de) < _DEGENERATE
352 inv2 = np.where(degenerate, 0., 1. / np.where(degenerate, 1., de) ** 2)
353 omega_n = -2 * np.sum(np.imag(vel[0] * vel[1].transpose(0, 2, 1)) * inv2, axis=2)
354 vdiag = np.real(np.diagonal(vel, axis1=2, axis2=3)) # (2, nk, n)
355 for e, mu in enumerate(mus):
356 f = occupation.fermi_dirac(en, float(mu), temperature)
357 window = f * (1. - f) / temperature
358 total[e] += np.einsum('kn,akn->a', window * omega_n, vdiag)
359 return _as_output(total / (len(kpts) * cell_area(ks)), e_fermi)
362def generalized_derivative(
363 en: NDArray[np.float64], vel: NDArray[np.complex128], w: NDArray[np.complex128],
364) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
365 r'''
366 Get the interband position matrix elements and their generalized
367 derivatives from the band energies and the first and second derivatives
368 of :math:`H` in the eigenbasis (a sum over states, Sipe and Shkrebtii
369 2000; Cook et al. 2017; :math:`\hbar = 1`, :math:`\omega_{nm} = E_n - E_m`,
370 :math:`\Delta^a_{nm} = v^a_{nn} - v^a_{mm}`):
372 .. math::
374 r^b_{nm} = \frac{v^b_{nm}}{i\omega_{nm}}\, ,\qquad
375 r^b_{nm;a} = \frac{i}{\omega_{nm}}\left[
376 \frac{v^b_{nm}\Delta^a_{nm} + v^a_{nm}\Delta^b_{nm}}{\omega_{nm}}
377 - w^{ba}_{nm} + \sum_{p\neq n,m}\left(
378 \frac{v^b_{np}v^a_{pm}}{\omega_{pm}} - \frac{v^a_{np}v^b_{pm}}{\omega_{np}}
379 \right)\right]\, ,
381 :math:`r^b_{nm;a} = \partial_ar^b_{nm} - i(A^a_{nn} - A^a_{mm})r^b_{nm}`,
382 gauge covariant. Degenerate pairs (:math:`|\omega_{nm}| < 10^{-9}`) are
383 set to zero, and dropped from the sum over :math:`p`.
385 Parameters
386 ----------
387 en : Real ndarray, shape (nk, n)
388 Band energies.
389 vel : Complex ndarray, shape (2, nk, n, n)
390 :math:`v^a_{nm}`.
391 w : Complex ndarray, shape (2, 2, nk, n, n)
392 :math:`w^{ab}_{nm} =
393 \langle n|\partial_a\partial_bH|m\rangle`.
395 Returns
396 -------
397 r : Complex ndarray, shape (2, nk, n, n): r[b] is :math:`r^b_{nm}`
398 r_der
399 Complex ndarray, shape (2, 2, nk, n, n): r_der[b, a]
400 is :math:`r^b_{nm;a}`.
401 '''
402 en, vel, w = np.asarray(en), np.asarray(vel), np.asarray(w)
403 error_handling.band_derivatives(en, vel, w)
404 om = en[:, :, None] - en[:, None, :]
405 degenerate = np.abs(om) < _DEGENERATE
406 inv = np.where(degenerate, 0., 1. / np.where(degenerate, 1., om))
407 diag = np.diagonal(vel, axis1=2, axis2=3) # (2, nk, n)
408 delta = diag[:, :, :, None] - diag[:, :, None, :]
409 r = -1j * vel * inv[None]
410 r_der = np.zeros((2, 2) + om.shape, 'c16')
411 for b in range(2):
412 for a in range(2):
413 # the full sums over p include p = n, m: those terms add
414 # Delta^b v^a / omega, removed again below
415 full = vel[b] @ (vel[a] * inv) - (vel[a] * inv) @ vel[b]
416 bracket = ((vel[b] * delta[a] + vel[a] * delta[b]) * inv - w[b, a] + full
417 - vel[a] * delta[b] * inv)
418 r_der[b, a] = 1j * inv * bracket
419 return r, r_der
422def shift_current(
423 ks: KSpace, omega: float | ArrayLike, e_fermi: float = 0., temperature: float = 0.,
424 eta: float = 0.05, nk: int | tuple[int, int] = 100, component: str = 'xxx',
425 positions: bool = True,
426) -> float | NDArray[np.float64]:
427 r'''
428 Get the shift-current conductivity of a 2D model under linearly
429 polarized light (Sipe and Shkrebtii, Phys. Rev. B 61, 5337 (2000)),
430 :math:`j^a = 2\sigma^{abc}(0;\omega,-\omega)E^b(\omega)E^c(-\omega)`,
432 .. math::
434 \sigma^{abc}(\omega) = -\frac{\pi e^3}{2\hbar^2}\int\frac{d^2k}{(2\pi)^2}
435 \sum_{n,m}f_{nm}\,\mathrm{Im}\left[r^b_{mn}r^c_{nm;a}
436 + r^c_{mn}r^b_{nm;a}\right]\delta(\omega_{mn} - \omega)\, ,
438 :math:`f_{nm} = f_n - f_m`, with the position matrix elements
439 :math:`r` and their generalized derivatives from
440 :func:`generalized_derivative` (a sum over states, with the exact second
441 derivatives of :math:`H`). The delta function is a Gaussian of standard
442 deviation :math:`\eta`. For :math:`b = c`,
443 :math:`\mathrm{Im}[r^b_{mn}r^b_{nm;a}] = |r^b_{nm}|^2R^{a,b}_{nm}` with
444 :math:`R^{a,b}_{nm}` the shift vector, the real-space displacement of
445 the electron in the transition. It vanishes with inversion symmetry.
446 Unlike the Hall conductance of full bands, it depends on where the
447 orbitals sit in the cell (Ibanez-Azpiroz, Tsirkin and Souza, Phys. Rev.
448 B 97, 245143 (2018)): ``positions=False``, which puts them all at the
449 cell origin, describes a different crystal.
451 Parameters
452 ----------
453 ks : **KSpace** instance, 2D, Hermitian, without overlap
454 omega : Real number or array of real numbers
455 Photon energies :math:`\hbar\omega`.
456 e_fermi : Real number, default 0
457 temperature : Positive real or zero, default 0
458 eta : Positive real, default 0.05
459 Gaussian width.
460 nk : Positive integer, or tuple of 2, default 100
461 component
462 Three letters among 'x', 'y' (current direction, then
463 the two field directions). Default value 'xxx'.
464 positions : Boolean, default True
465 See *KSpace.hall_conductivity*.
467 Returns
468 -------
469 sigma
470 Real number (or ndarray shaped like *omega*), in units
471 of :math:`e^3/\hbar\times` length / energy (:math:`\hbar = e = 1`),
472 counting each orbital once.
473 '''
474 kpts = _mesh(ks, nk, positions)
475 error_handling.no_overlap(bool(ks._overlap_hop))
476 error_handling.frequencies(omega, 'omega')
477 error_handling.real_number(e_fermi, 'e_fermi')
478 error_handling.positive_real_zero(temperature, 'temperature')
479 error_handling.positive_real(eta, 'eta')
480 a, b, c = error_handling.tensor_component(component, 3)
481 w_arr = np.atleast_1d(np.asarray(omega, dtype='f8')).ravel()
482 total = np.zeros(len(w_arr))
483 for en, vel, w in _bands(ks, kpts, positions, second=True):
484 r, r_der = generalized_derivative(en, vel, w)
485 f = occupation.fermi_dirac(en, e_fermi, temperature)
486 fnm = f[:, :, None] - f[:, None, :]
487 integrand = np.imag(r[b].transpose(0, 2, 1) * r_der[c, a]
488 + r[c].transpose(0, 2, 1) * r_der[b, a]) * fnm
489 keep = integrand != 0
490 wmn = (en[:, None, :] - en[:, :, None])[keep]
491 total += _gaussian_sum(w_arr, wmn, integrand[keep], eta)
492 sigma = -np.pi / 2 * total / (len(kpts) * cell_area(ks))
493 return _as_output(sigma, omega)