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