Coverage for calorine/tools/spectra.py: 99%
162 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-18 13:01 +0000
1r"""
2Tools for computing optical and vibrational spectra from MD trajectories.
4Covers three routes depending on the model and GPUMD output:
6* **Dielectric tensor** (:func:`get_dielectric_function`): full 3×3 tensor
7 from ``dpdt.out`` (qNEP, ``column='dP'`` or ``column='P'``) or ``dipole.out``
8 (TNEP, ``column='mu'``).
9* **IR spectrum** (:func:`get_ir_spectrum`): from ``dpdt.out``
10 (qNEP, ``column='dP'`` or ``column='P'``) or ``dipole.out``
11 (TNEP, ``column='mu'``).
12* **Raman spectrum** (:func:`get_raman_spectrum`): from ``polarizability.out``
13 (TNEP polarizability/susceptibility model).
15Nomenclature
16------------
17- **Molecules / nanoparticles (localized systems)**: dipole moment :math:`\mu`
18 (e·Å) and polarizability :math:`\alpha` (bohr³ or training-data units).
19- **Solids / liquids (extended systems)**: polarization :math:`P` (e·Å, total
20 cell dipole) and susceptibility :math:`\chi` (same unit as training data per
21 cell).
23Both cases use the same GPUMD output files and formulas; the distinction is
24whether :attr:`volume` is supplied for normalization.
26Quantum corrections
27-------------------
28The spectra from classical MD can be corrected for quantum statistics via
29:func:`apply_quantum_correction`. The correction factors differ by scattering
30order (M. Cardona, *Topics in Applied Physics*, Vol. 50 (Springer, 1982);
31Rosander *et al.*, PRB **111**, 064107 (2025)).
32"""
34import numpy as np
35from pandas import DataFrame
36from ase.units import _k as kB_SI, _e as e_SI, _eps0 as eps0_SI, _c as c_SI, _hbar as hbar_SI
38_rad_s_to_thz = 1e-12
39_rad_s_to_invcm = 1 / (2 * np.pi * c_SI * 100)
42def _compute_correlation_function(Z1, Z2):
43 """Normalized cross-correlation; normalization divides by the number of contributing pairs."""
44 N = len(Z1)
45 return np.correlate(Z1, Z2, mode='full')[N - 1:] / np.arange(N, 0, -1)
48def _gaussian_decay(t, t_sigma):
49 r"""Gaussian envelope :math:`\exp(-t^2 / 2\sigma^2)` used as a window applied to the ACF."""
50 return np.exp(-0.5 * (t / t_sigma) ** 2)
53def _psd_from_acf(acf, dt):
54 """Power spectral density via even extension of the ACF followed by rfft."""
55 signal = np.hstack((acf, acf[:0:-1]))
56 fft = dt * np.fft.rfft(signal)
57 freqs = 2 * np.pi * np.fft.rfftfreq(len(signal), dt)
58 return freqs, fft.real
61def _parse_polarization_pair(polarization_in, polarization_out):
62 """Validate and parse a pair of polarization unit vectors, or return (None, None)."""
63 if (polarization_in is None) != (polarization_out is None):
64 raise ValueError('polarization_in and polarization_out must be given together')
65 if polarization_in is None:
66 return None, None
67 n_in = np.asarray(polarization_in, dtype=float)
68 n_out = np.asarray(polarization_out, dtype=float)
69 if n_in.shape != (3,) or n_out.shape != (3,):
70 raise ValueError('polarization_in and polarization_out must have shape (3,)')
71 return n_in, n_out
74# ── Signal column mapping ──────────────────────────────────────────────────────
76_SIGNAL_COLS = {
77 'dP': (['dPx', 'dPy', 'dPz'], True),
78 'P': (['Px', 'Py', 'Pz'], False),
79 'mu': (['mu_x', 'mu_y', 'mu_z'], False),
80}
83def _resolve_signal(signal, column, dt):
84 """Select signal columns from a DataFrame, infer derivative flag, auto-extract dt."""
85 if column not in _SIGNAL_COLS:
86 raise ValueError(
87 f'column must be one of {list(_SIGNAL_COLS)!r}, got {column!r}'
88 )
89 col_names, derivative = _SIGNAL_COLS[column]
90 missing = [c for c in col_names if c not in signal.columns]
91 if missing:
92 raise ValueError(
93 f'signal is missing columns {missing!r} for column={column!r}'
94 )
95 if dt is None: 95 ↛ 102line 95 didn't jump to line 102 because the condition on line 95 was always true
96 if 'time' not in signal.columns:
97 raise ValueError(
98 'dt must be provided when signal has no time column '
99 f'(got columns: {list(signal.columns)!r})'
100 )
101 dt = float(signal['time'].diff().dropna().iloc[0])
102 return signal[col_names].to_numpy(dtype=float), dt, derivative
105# ── IR spectrum ────────────────────────────────────────────────────────────────
107def _acf_to_psd(acf, dt, t_sigma, max_lag):
108 """Truncate ACF, apply optional Gaussian window, and compute PSD. dt and t_sigma in ps."""
109 acf = np.real(acf[:int(len(acf) * max_lag)])
110 t_acf = dt * np.arange(len(acf)) # ps
111 if t_sigma is not None:
112 acf = acf * _gaussian_decay(t_acf, t_sigma=t_sigma) # both ps
113 return _psd_from_acf(acf, dt * 1e-12) # _psd_from_acf uses seconds
116def _ir_spectrum_1d(dt, signal, volume_SI, temperature, t_sigma, max_lag, derivative=False):
117 """Compute IR spectrum for a single dipole/polarization component in SI units."""
118 D = e_SI * 1e-10 / (1e-15 if derivative else 1.0)
119 signal_SI = (signal - np.mean(signal)) * D
120 acf = _compute_correlation_function(signal_SI, signal_SI)
121 w, S = _acf_to_psd(acf, dt, t_sigma, max_lag)
122 # When input is the signal itself, multiply by omega**2 to get PSD of its derivative
123 # (PSD[dx/dt] = omega**2 * PSD[x]). When input is already the derivative, use PSD directly.
124 S_dpdt = S if derivative else w ** 2 * S
126 if volume_SI is not None:
127 beta = 1.0 / (kB_SI * temperature)
128 conductivity = (beta / (3.0 * volume_SI)) * S_dpdt
129 with np.errstate(divide='ignore', invalid='ignore'):
130 eps_imag = conductivity / (eps0_SI * w)
131 return w, eps_imag, conductivity
132 else:
133 return w, S_dpdt
136def get_ir_spectrum(
137 signal: DataFrame,
138 volume: float = None,
139 temperature: float = None,
140 column: str = 'dP',
141 dt: float = None,
142 polarization: np.ndarray = None,
143 t_sigma: float = None,
144 max_lag: float = 0.25,
145) -> DataFrame:
146 r"""Compute the IR spectrum from a time series of dipole moments or polarizations.
148 Parameters
149 ----------
150 signal
151 ``DataFrame`` as returned by :func:`~calorine.gpumd.read_dpdt` or
152 :func:`~calorine.gpumd.read_dipole`. The columns to use are selected
153 via :attr:`column`.
154 volume
155 Simulation cell volume in ų. Required for extended systems; when
156 provided the function returns the imaginary dielectric function and
157 conductivity. Pass ``None`` for molecules to obtain an unnormalized
158 line shape in arbitrary units.
159 temperature
160 Temperature in K. Required when :attr:`volume` is not ``None``.
161 column
162 Selects which columns of :attr:`signal` to use and implies whether
163 the input is a time derivative:
165 ``'dP'`` (default)
166 Uses ``dPx``, ``dPy``, ``dPz`` from :func:`~calorine.gpumd.read_dpdt`
167 (qNEP); signal is :math:`\mathrm{d}P/\mathrm{d}t` in e·Å/fs.
168 ``'P'``
169 Uses ``Px``, ``Py``, ``Pz`` from :func:`~calorine.gpumd.read_dpdt`
170 (qNEP); signal is the polarization :math:`P(t)` in e·Å.
171 ``'mu'``
172 Uses ``mu_x``, ``mu_y``, ``mu_z`` from
173 :func:`~calorine.gpumd.read_dipole` (TNEP); signal is the dipole
174 moment :math:`\mu(t)` in e·Å.
176 dt
177 Time between consecutive frames in ps. Auto-extracted from the
178 ``time`` column when present (i.e. for :func:`~calorine.gpumd.read_dpdt`
179 output). Must be supplied explicitly for step-indexed inputs such as
180 :func:`~calorine.gpumd.read_dipole` output.
181 polarization
182 Unit vector ``(3,)`` defining the electric-field polarization direction.
183 When given, the spectrum is computed from the projected signal
184 :math:`s(t) = \hat{n} \cdot \mathrm{signal}(t)` instead of the
185 isotropic average.
186 t_sigma
187 Width of the Gaussian window applied to the ACF in ps.
188 ``None`` uses no windowing.
189 max_lag
190 Maximum ACF length as a fraction of the trajectory length; must be
191 in (0, 1]. Defaults to 0.25 (25 %).
193 Returns
194 -------
195 DataFrame
196 Always contains ``angular_frequency`` (THz) and
197 ``wavenumber_invcm`` (:math:`\mathrm{cm}^{-1}`).
198 When :attr:`volume` and :attr:`temperature` are given:
199 ``epsilon_imag`` (dimensionless) and ``conductivity`` (S/m).
200 When :attr:`volume` is ``None``: ``ir_intensity`` (arbitrary units,
201 proportional to :math:`\mathrm{PSD}[\dot{\mu}]`).
202 """
203 if volume is not None and temperature is None:
204 raise ValueError('temperature must be provided when volume is given')
206 arr, dt, derivative = _resolve_signal(signal, column, dt)
207 volume_SI = volume * 1e-30 if volume is not None else None
209 if polarization is not None:
210 n = np.asarray(polarization, dtype=float)
211 if n.shape != (3,):
212 raise ValueError('polarization must have shape (3,)')
213 components = [arr @ n]
214 else:
215 components = [arr[:, i] for i in range(3)]
217 results = [_ir_spectrum_1d(dt, c, volume_SI, temperature, t_sigma, max_lag, derivative)
218 for c in components]
220 if volume_SI is not None:
221 w = results[0][0]
222 eps_imag = np.mean([r[1] for r in results], axis=0)
223 conductivity = np.mean([r[2] for r in results], axis=0)
224 mask = np.isfinite(eps_imag)
225 cols = ['angular_frequency', 'wavenumber_invcm', 'epsilon_imag', 'conductivity']
226 data = np.column_stack((w[mask] * _rad_s_to_thz, w[mask] * _rad_s_to_invcm,
227 eps_imag[mask], conductivity[mask]))
228 return DataFrame(data, columns=cols)
229 else:
230 w = results[0][0]
231 ir_intensity = np.mean([r[1] for r in results], axis=0)
232 cols = ['angular_frequency', 'wavenumber_invcm', 'ir_intensity']
233 data = np.column_stack((w * _rad_s_to_thz, w * _rad_s_to_invcm, ir_intensity))
234 return DataFrame(data, columns=cols)
237# ── Raman spectrum ─────────────────────────────────────────────────────────────
239_POL_COLS = ['xx', 'yy', 'zz', 'yz', 'xz', 'xy']
242def get_raman_spectrum(
243 dt: float,
244 polarizability: DataFrame,
245 polarization_in: np.ndarray = None,
246 polarization_out: np.ndarray = None,
247 t_sigma: float = None,
248 max_lag: float = 0.25,
249) -> DataFrame:
250 r"""Compute the Raman spectrum from a time series of polarizability tensors.
252 Parameters
253 ----------
254 dt
255 Time between consecutive frames in ps.
256 polarizability
257 DataFrame with columns ``xx``, ``yy``, ``zz``, ``xy``, ``yz``, ``xz``
258 containing the polarizability :math:`\alpha` (molecules) or susceptibility
259 :math:`\chi` (extended systems), as returned by
260 :func:`~calorine.gpumd.read_polarizability`. Units are those of the
261 TNEP training data (typically bohr³).
262 polarization_in
263 Unit vector ``(3,)`` for the polarization of the incoming light. If
264 both :attr:`polarization_in` and :attr:`polarization_out` are given,
265 the polarization-resolved intensity
266 :math:`I(\omega) \propto \mathrm{FT}[\langle s(0)s(t)\rangle]` with
267 :math:`s(t) = \hat{n}^\mathrm{out} \cdot \alpha(t) \cdot \hat{n}^\mathrm{in}`
268 is added as the column ``raman_polarized``.
269 polarization_out
270 Unit vector ``(3,)`` for the polarization of the outgoing (scattered)
271 light.
272 t_sigma
273 Width of the Gaussian window applied to the ACF in ps. ``None`` uses no
274 windowing.
275 max_lag
276 Maximum ACF length as a fraction of the trajectory length; must be
277 in (0, 1]. Defaults to 0.25 (25 %).
279 Returns
280 -------
281 DataFrame
282 Always contains ``angular_frequency`` (THz),
283 ``wavenumber_invcm`` (:math:`\mathrm{cm}^{-1}`),
284 ``raman_isotropic`` (proportional to
285 :math:`\mathrm{FT}[\langle\gamma(0)\gamma(t)\rangle]`), and
286 ``raman_anisotropic`` (proportional to
287 :math:`\mathrm{FT}[\langle\mathrm{Tr}[\beta(0)\beta(t)]\rangle]`).
288 If both polarization vectors are given, also contains
289 ``raman_polarized``.
290 """
291 missing = [c for c in _POL_COLS if c not in polarizability.columns]
292 if missing:
293 raise ValueError(f'polarizability is missing columns: {missing}')
295 # Remove mean (Rayleigh component) before computing fluctuation spectrum.
296 pol = polarizability[_POL_COLS] - polarizability[_POL_COLS].mean()
298 n_in, n_out = _parse_polarization_pair(polarization_in, polarization_out)
300 # Isotropic component gamma(t) = Tr[alpha(t)] / 3.
301 gamma = (pol['xx'] + pol['yy'] + pol['zz']).to_numpy() / 3.0
303 # Traceless part beta = alpha - gamma*I.
304 beta_xx = pol['xx'].to_numpy() - gamma
305 beta_yy = pol['yy'].to_numpy() - gamma
306 beta_zz = pol['zz'].to_numpy() - gamma
307 # Off-diagonal elements are unchanged: beta_xy = alpha_xy, etc.
309 with np.errstate(invalid='ignore', divide='ignore'):
310 # Isotropic Raman: FT[<gamma(0)*gamma(t)>].
311 acf_iso = _compute_correlation_function(gamma, gamma)
312 w, L_iso = _acf_to_psd(acf_iso, dt, t_sigma, max_lag)
314 # Anisotropic Raman: FT[<Tr[beta(0)*beta(t)]>].
315 # Tr[beta(0)*beta(t)] = sum_ij beta_ij(0)*beta_ij(t) — sum of per-element
316 # ACFs with factor 2 for off-diagonal (since beta_ij = beta_ji).
317 acf_aniso = (
318 _compute_correlation_function(beta_xx, beta_xx)
319 + _compute_correlation_function(beta_yy, beta_yy)
320 + _compute_correlation_function(beta_zz, beta_zz)
321 + 2.0 * _compute_correlation_function(pol['yz'].to_numpy(), pol['yz'].to_numpy())
322 + 2.0 * _compute_correlation_function(pol['xz'].to_numpy(), pol['xz'].to_numpy())
323 + 2.0 * _compute_correlation_function(pol['xy'].to_numpy(), pol['xy'].to_numpy())
324 )
325 _, L_aniso = _acf_to_psd(acf_aniso, dt, t_sigma, max_lag)
327 cols = ['angular_frequency', 'wavenumber_invcm', 'raman_isotropic', 'raman_anisotropic']
328 data = np.column_stack((w * _rad_s_to_thz, w * _rad_s_to_invcm, L_iso, L_aniso))
329 df = DataFrame(data, columns=cols)
331 # Polarization-resolved spectrum: s(t) = n_out . alpha(t) . n_in.
332 # s = sum_ij n_out_i * alpha_ij * n_in_j; use symmetry alpha_ij = alpha_ji.
333 if n_in is not None:
334 s = (n_in[0]*n_out[0] * pol['xx'].to_numpy()
335 + n_in[1]*n_out[1] * pol['yy'].to_numpy()
336 + n_in[2]*n_out[2] * pol['zz'].to_numpy()
337 + (n_in[1]*n_out[2] + n_in[2]*n_out[1]) * pol['yz'].to_numpy()
338 + (n_in[0]*n_out[2] + n_in[2]*n_out[0]) * pol['xz'].to_numpy()
339 + (n_in[0]*n_out[1] + n_in[1]*n_out[0]) * pol['xy'].to_numpy())
340 with np.errstate(invalid='ignore', divide='ignore'):
341 acf_pol = _compute_correlation_function(s, s)
342 _, L_pol = _acf_to_psd(acf_pol, dt, t_sigma, max_lag)
343 df['raman_polarized'] = L_pol
345 return df
348# ── Quantum correction ─────────────────────────────────────────────────────────
350def apply_quantum_correction(
351 df: DataFrame,
352 temperature: float,
353 column: str,
354 order: str = 'first',
355 force: bool = False,
356) -> DataFrame:
357 r"""Apply a quantum correction to a classically computed spectrum.
359 Classical MD underestimates spectral intensities because it samples the
360 classical Boltzmann distribution rather than the Bose-Einstein distribution.
361 The correction factor depends on the scattering order and mode type
362 [Cardona1982]_ [Rosander2025]_.
364 Parameters
365 ----------
366 df
367 DataFrame with an ``angular_frequency`` column (THz) and the column
368 to correct.
369 temperature
370 Temperature in K.
371 column
372 Name of the spectral column to correct (e.g. ``'raman_isotropic'``,
373 ``'ir_intensity'``, ``'epsilon_imag'``).
374 order
375 Correction type:
377 ``'first'``
378 First-order scattering (IR absorption and first-order Raman):
379 :math:`f = \beta\hbar\omega / (1 - e^{-\beta\hbar\omega})`.
380 ``'overtone'``
381 Second-order overtone (same mode, :math:`\omega_1 = \omega/2`):
382 :math:`f = [y/(1 - e^{-y})]^2 (2 - e^{-y})` with
383 :math:`y = \beta\hbar\omega/2`.
384 ``'combination'``
385 Second-order combination band upper bound
386 (:math:`\omega_1 = \omega_2 = \omega/2`):
387 :math:`f = [y/(1 - e^{-y})]^2` with :math:`y = \beta\hbar\omega/2`.
388 force
389 If ``True``, overwrite an existing ``column + '_qm'`` column.
390 Default ``False`` raises :class:`ValueError` to prevent accidental
391 double-correction.
393 Returns
394 -------
395 DataFrame
396 Copy of :attr:`df` with a new column ``column + '_qm'`` containing
397 the quantum-corrected intensities. The input is not mutated.
399 References
400 ----------
401 .. [Cardona1982] M. Cardona, *Resonance phenomena*, in *Topics in Applied
402 Physics*, Vol. 50, edited by M. Cardona and G. Güntherodt
403 (Springer, Berlin, 1982).
404 .. [Rosander2025] Rosander *et al.*, Phys. Rev. B **111**, 064107 (2025),
405 Supp. Eqs. S5, S7, S10.
406 """
407 if order not in ('first', 'overtone', 'combination'):
408 raise ValueError("order must be 'first', 'overtone', or 'combination'")
410 if not force and column + '_qm' in df.columns:
411 raise ValueError(
412 f"column '{column}_qm' already exists; pass force=True to overwrite"
413 )
415 w = df['angular_frequency'].to_numpy() * 1e12 # THz → rad/s
416 beta = 1.0 / (kB_SI * temperature)
417 x = beta * hbar_SI * w # beta*hbar*omega (dimensionless)
419 with np.errstate(divide='ignore', invalid='ignore', over='ignore'):
420 if order == 'first':
421 # f = beta*hbar*omega / (1 - exp(-beta*hbar*omega))
422 factor = np.where(x > 0, x / (1.0 - np.exp(-x)), 1.0)
423 elif order == 'overtone':
424 # f = [y/(1-exp(-y))]**2 * (2-exp(-y)), y = beta*hbar*omega/2
425 y = x / 2.0
426 base = np.where(y > 0, y / (1.0 - np.exp(-y)), 1.0)
427 factor = base ** 2 * np.where(y > 0, 2.0 - np.exp(-y), 1.0)
428 else: # combination
429 # f = [y/(1-exp(-y))]**2, y = beta*hbar*omega/2
430 y = x / 2.0
431 base = np.where(y > 0, y / (1.0 - np.exp(-y)), 1.0)
432 factor = base ** 2
434 result = df.copy()
435 result[column + '_qm'] = df[column].to_numpy() * factor
436 return result
439# ── Dielectric tensor ─────────────────────────────────────────────────────────
441_TENSOR_COMPONENTS = [(0, 0, 'xx'), (1, 1, 'yy'), (2, 2, 'zz'),
442 (1, 2, 'yz'), (0, 2, 'xz'), (0, 1, 'xy')]
445def get_dielectric_function(
446 signal: DataFrame,
447 volume: float,
448 temperature: float,
449 column: str = 'dP',
450 dt: float = None,
451 t_sigma: float = None,
452 max_lag: float = 0.25,
453 return_real_part: bool = True,
454 kk_method: str = 'vectorized',
455) -> DataFrame:
456 r"""Compute the dielectric function from a three-component polarization time series.
458 Computes all six unique Voigt components of :math:`\epsilon_2^{\alpha\beta}(\omega)`
459 (xx, yy, zz, yz, xz, xy) from the three-component time series written by GPUMD.
460 The real part :math:`\epsilon_1^{\alpha\beta}(\omega)` is obtained via the
461 Kramers-Kronig relation when :attr:`return_real_part` is ``True`` (the default).
463 The imaginary part is related to the symmetrized cross-correlation via
465 .. math::
467 \epsilon_2^{\alpha\beta}(\omega) = \frac{\beta}{\epsilon_0 \omega V}
468 \,\mathrm{Re}\!\left[\int_0^\infty
469 \langle \dot{P}_\alpha(0)\,\dot{P}_\beta(t) \rangle
470 e^{-i\omega t}\,dt\right],
472 where :math:`\beta = 1/(k_\mathrm{B}T)` and :math:`V` is the cell volume.
474 Parameters
475 ----------
476 signal
477 ``DataFrame`` as returned by :func:`~calorine.gpumd.read_dpdt` or
478 :func:`~calorine.gpumd.read_dipole`. The columns to use are selected
479 via :attr:`column`.
480 volume
481 Simulation cell volume in ų.
482 temperature
483 Temperature in K.
484 column
485 Selects which columns of :attr:`signal` to use and implies whether
486 the input is a time derivative:
488 ``'dP'`` (default)
489 Uses ``dPx``, ``dPy``, ``dPz`` from :func:`~calorine.gpumd.read_dpdt`
490 (qNEP); input is :math:`\dot{\mathbf{P}}(t)` in e·Å/fs.
491 ``'P'``
492 Uses ``Px``, ``Py``, ``Pz`` from :func:`~calorine.gpumd.read_dpdt`
493 (qNEP); input is the polarization :math:`\mathbf{P}(t)` in e·Å.
494 ``'mu'``
495 Uses ``mu_x``, ``mu_y``, ``mu_z`` from
496 :func:`~calorine.gpumd.read_dipole` (TNEP); input is the dipole
497 moment :math:`\boldsymbol{\mu}(t)` in e·Å.
499 dt
500 Time between consecutive frames in ps. Auto-extracted from the
501 ``time`` column when present (i.e. for :func:`~calorine.gpumd.read_dpdt`
502 output). Must be supplied explicitly for step-indexed inputs such as
503 :func:`~calorine.gpumd.read_dipole` output.
504 t_sigma
505 Width of the Gaussian window applied to the ACF in ps.
506 ``None`` uses no windowing.
507 max_lag
508 Maximum ACF length as a fraction of the trajectory length; must be
509 in (0, 1]. Defaults to 0.25 (25 %).
510 return_real_part
511 When ``True`` (default), the real part :math:`\epsilon_1^{\alpha\beta}(\omega)`
512 is computed via the Kramers-Kronig relation and appended as
513 ``epsilon_real_{xx,yy,zz,yz,xz,xy}`` columns.
514 kk_method
515 Integration method passed to :func:`~calorine.tools.apply_kramers_kronig`
516 when :attr:`return_real_part` is ``True``. ``'vectorized'`` (default) uses
517 an exact trapezoid rule; ``'fft'`` uses a faster Hilbert transform
518 approximation. Ignored when :attr:`return_real_part` is ``False``.
520 Returns
521 -------
522 DataFrame
523 Contains ``angular_frequency`` (THz), ``wavenumber_invcm``
524 (:math:`\mathrm{cm}^{-1}`),
525 ``epsilon_imag_{xx,yy,zz,yz,xz,xy}`` (dimensionless, imaginary dielectric
526 tensor components), and ``conductivity_{xx,yy,zz,yz,xz,xy}`` (S/m).
527 When :attr:`return_real_part` is ``True``, ``epsilon_real_{xx,yy,zz,yz,xz,xy}``
528 columns are appended.
529 """
530 arr, dt, derivative = _resolve_signal(signal, column, dt)
532 D = e_SI * 1e-10 / (1e-15 if derivative else 1.0)
533 volume_SI = volume * 1e-30
534 beta = 1.0 / (kB_SI * temperature)
536 signal_SI = (arr - arr.mean(axis=0)) * D
538 def _psd(z1, z2=None):
539 if z2 is None:
540 acf = _compute_correlation_function(z1, z1)
541 else:
542 acf = (_compute_correlation_function(z1, z2)
543 + _compute_correlation_function(z2, z1)) / 2
544 w, S = _acf_to_psd(acf, dt, t_sigma, max_lag)
545 # PSD[dx/dt] = omega**2 * PSD[x] when input is the signal itself.
546 return w, S if derivative else w ** 2 * S
548 w = None
549 conductivities = {}
550 eps_imags = {}
551 for i, j, lbl in _TENSOR_COMPONENTS:
552 zi = signal_SI[:, i]
553 zj = signal_SI[:, j]
554 w, S_dpdt = _psd(zi) if i == j else _psd(zi, zj)
555 conductivity = (beta / volume_SI) * S_dpdt
556 with np.errstate(divide='ignore', invalid='ignore'):
557 eps_imag = conductivity / (eps0_SI * w)
558 conductivities[lbl] = conductivity
559 eps_imags[lbl] = eps_imag
561 mask = w > 0
562 cols = (['angular_frequency', 'wavenumber_invcm']
563 + [f'epsilon_imag_{lbl}' for _, _, lbl in _TENSOR_COMPONENTS]
564 + [f'conductivity_{lbl}' for _, _, lbl in _TENSOR_COMPONENTS])
565 arrays = ([w[mask] * _rad_s_to_thz, w[mask] * _rad_s_to_invcm]
566 + [eps_imags[lbl][mask] for _, _, lbl in _TENSOR_COMPONENTS]
567 + [conductivities[lbl][mask] for _, _, lbl in _TENSOR_COMPONENTS])
568 result = DataFrame(np.column_stack(arrays), columns=cols)
569 if return_real_part:
570 from .kramers_kronig import apply_kramers_kronig
571 result = apply_kramers_kronig(result, method=kk_method)
572 return result