Coverage for strongcoca / response / dielectric.py: 100%
134 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-04-15 18:15 +0000
1import logging
2from abc import ABC, abstractmethod
3from typing import Union
5import numpy as np
7from ..types import Array, ComplexArray, Path_str
8from ..units import au_to_eV, eV_to_au
9from ..utilities import ClassFormatter
12logger = logging.getLogger(__name__)
15class DielectricFunction(ABC):
17 """Objects of this class hold different representations of dielectric functions.
19 Parameters
20 ----------
21 name
22 Name of dielectric function.
23 """
25 def __init__(self,
26 name: str = 'DielectricFunction') -> None:
27 logger.debug(f'Entering {self.__class__.__name__}.__init__')
28 self._name = name
30 def __call__(self, frequencies: Array) -> np.ndarray:
31 """Returns the dielectric function at the given frequencies
33 Parameters
34 ----------
35 frequencies
36 List of frequencies at which the dielectric function is returned;
37 in units of eV.
38 """
39 freq_w = np.array(frequencies, dtype=float)
40 freq_w *= eV_to_au
42 return self._eval_at(freq_w)
44 @abstractmethod
45 def _eval_at(self, freq_w: np.ndarray) -> np.ndarray:
46 """Returns the dielectric function at the given frequencies
48 Parameters
49 ----------
50 frequencies
51 List of frequencies at which the dielectric function is returned;
52 in atomic units.
53 """
54 raise NotImplementedError()
56 @abstractmethod
57 def __str__(self) -> str:
58 raise NotImplementedError()
60 def _repr_html_(self) -> str:
61 """HTML representation for Jupyter notebooks."""
62 rows = [f'<tr><td style="text-align: left;">name</td><td>{self._name}</td></tr>']
63 return (
64 f'<h4>{self.__class__.__name__}</h4>'
65 '<table>'
66 '<thead><tr>'
67 '<th style="text-align: left;">field</th>'
68 '<th style="text-align: left;">value</th>'
69 '</tr></thead>'
70 '<tbody>' + ''.join(rows) + '</tbody>'
71 '</table>'
72 )
74 @property
75 def name(self) -> str:
76 """Name of dielectric function."""
77 return self._name
80class AnalyticDielectricFunction(DielectricFunction):
82 """Objects of this class hold different representations of dielectric functions that
83 have analytic form. The dielectric function can thus be evaluated at any frequency
84 without interpolation.
86 Parameters
87 ----------
88 name
89 Name of dielectric function.
90 """
92 def __init__(self,
93 name: str = 'AnalyticDielectricFunction') -> None:
95 logger.debug(f'Entering {self.__class__.__name__}.__init__')
96 super().__init__(name=name)
99class NumericDielectricFunction(DielectricFunction):
101 """Objects of this class hold different representations of dielectric functions sampled
102 on a discrete grid of frequencies. The dielectric function can be evaluated at any frequency
103 within the supplied grid by linear interpolation.
105 Parameters
106 ----------
107 frequencies
108 List of frequencies at which the dielectric function :attr:`dielectric_function`
109 is sampled; units are eV by default, optionally atomic units
110 (see :attr:`units`).
112 The list of frequencies and the corresponding dielectric function are sorted if the former
113 is not monotonically increasing.
114 dielectric_function
115 Complex dielectric function sampled at the frequencies :attr:`frequencies`.
116 units
117 `eV` to specify :attr:`freq_w` in eV, or `au` to specify in atomic units.
119 This parameter determines whether conversion should be performed during
120 initialization and has no effect on instance methods and variables.
121 name
122 Name of dielectric function.
124 Raises
125 ------
126 ValueError
127 If :attr:`frequencies` or :attr:`dielectric_function` have zero length.
128 ValueError
129 If there is a mismatch in the length of :attr:`frequencies` and :attr:`dielectric_function`.
130 ValueError
131 If :attr:`units` is not one of the allowed values `eV` or `au`.
133 Examples
134 --------
136 The following snippet illustrates how a NumericDielectricFunction object can be constructed.
137 Here, the analytic Drude dielectric function of the free-electron gas is used.
138 To construct a dielectric function object from a numeric dielectric function the function
139 :func:`~strongcoca.response.read_dielec_function_file` can be used.
141 >>> import numpy as np
142 >>> from strongcoca.response.dielectric import NumericDielectricFunction
143 >>>
144 >>> omega_p = 3
145 >>> gamma = 0.5
146 >>> freq_w = np.linspace(0.1, 10)
147 >>> eps_w = 1 - omega_p**2 / (freq_w**2 + 1j*gamma*freq_w)
148 >>> dielec = NumericDielectricFunction(freq_w, eps_w, name='Drude-function')
149 >>> print(dielec)
150 ---------------- NumericDielectricFunction -----------------
151 name : Drude-function
152 Frequency grid : 50 points between 0.10 and 10.00eV
153 frequencies (eV) : [ 0.1 0.30204 0.50408 ...
154 9.59592 9.79796 10. ]
155 dielectric_function : [-33.61538+1.73077e+02j
156 -25.37528+4.36618e+01j
157 -16.85366+1.77091e+01j ...
158 0.90253+5.07897e-03j
159 0.90649+4.77173e-03j
160 0.91022+4.48878e-03j]
161 >>> dense_freq_w = np.linspace(3, 4)
162 >>> print(dielec(dense_freq_w))
163 [0.02425453+0.16310249j 0.0367996 +0.15996402j 0.04934466+0.15682555j
164 0.06188972+0.15368708j 0.07443478+0.15054861j 0.08697984+0.14741014j
165 ...
166 """
168 def __init__(self,
169 frequencies: Union[Array, float],
170 dielectric_function: Union[ComplexArray, complex],
171 units: str = 'eV',
172 name: str = 'NumericDielectricFunction') -> None:
174 logger.debug(f'Entering {self.__class__.__name__}.__init__')
175 super().__init__(name=name)
177 if isinstance(frequencies, float):
178 frequencies = [frequencies]
180 if isinstance(dielectric_function, complex) or isinstance(dielectric_function, float):
181 dielectric_function = [dielectric_function]
183 freq_w = np.array(frequencies)
184 eps_w = np.array(dielectric_function, dtype=complex)
186 if units == 'eV':
187 freq_w = freq_w * eV_to_au
188 elif units != 'au':
189 raise ValueError(f"units has to be 'eV' or 'au', not '{units}'")
191 if len(freq_w) != len(eps_w):
192 raise ValueError(f'length of frequencies {len(freq_w)} must match length of '
193 f'dielectric function {len(eps_w)}')
195 if len(freq_w) == 0:
196 raise ValueError('length of frequencies and dielectric function must not be zero')
198 # np.diff works if input has at least length 2
199 if len(freq_w) > 1 and not np.all(np.diff(freq_w) > 0):
200 sort_w = np.argsort(freq_w)
201 freq_w = freq_w[sort_w]
202 eps_w = eps_w[sort_w]
204 self._freq_w = freq_w
205 self._eps_w = eps_w
207 def _eval_at(self, freq_w: np.ndarray) -> np.ndarray:
208 """Returns the dielectric function at the given frequencies by interpolating the numeric
209 dielectric function.
211 Parameters
212 ----------
213 frequencies
214 List of frequencies at which the dielectric function is returned;
215 in units of eV.
217 Raises
218 ------
219 ValueError
220 If any component of :attr:`frequencies` is outside of the interval of the numeric
221 dielectric function.
222 """
223 atol = np.finfo(float).eps * max(abs(self._freq_w[0]), abs(self._freq_w[-1]))
224 if np.min(freq_w) < self._freq_w[0] - atol:
225 raise ValueError('Minimum of array frequencies is outside the lower bound of the'
226 'numeric dielectric function')
227 if np.max(freq_w) > self._freq_w[-1] + atol:
228 raise ValueError('Maximum of array frequencies is outside the upper bound of the'
229 'numeric dielectric function')
231 eps_w = np.interp(freq_w, self._freq_w, self._eps_w)
233 return eps_w # type: ignore
235 def _repr_html_(self) -> str:
236 """HTML representation for Jupyter notebooks."""
237 rows = [
238 f'<tr><td style="text-align: left;">name</td><td>{self._name}</td></tr>',
239 f'<tr><td style="text-align: left;">frequency grid</td>'
240 f'<td>{self.n} points between {self.minfreq:.2f} and {self.maxfreq:.2f} eV</td></tr>',
241 ]
242 return (
243 f'<h4>{self.__class__.__name__}</h4>'
244 '<table>'
245 '<thead><tr>'
246 '<th style="text-align: left;">field</th>'
247 '<th style="text-align: left;">value</th>'
248 '</tr></thead>'
249 '<tbody>' + ''.join(rows) + '</tbody>'
250 '</table>'
251 )
253 def __str__(self) -> str:
254 fmt = ClassFormatter(self, pad=20)
255 fmt.append_class_name()
257 fmt.append_attr('name')
259 fmt.append_attr('Frequency grid', f'{self.n} points between {self.minfreq:.2f} and '
260 f'{self.maxfreq:.2f}eV')
262 fmt.append_attr('frequencies', unit='eV')
263 fmt.append_attr('dielectric_function')
265 return fmt.to_string()
267 @property
268 def n(self) -> int:
269 """number of samples for the dielectric function
270 :attr:`~strongcoca.response.NumericDielectricFunction.dielectric_function`."""
271 return len(self._eps_w)
273 @property
274 def minfreq(self) -> float:
275 """Lowest frequency :attr:`~strongcoca.response.NumericDielectricFunction.frequencies`
276 for which there is a value of the dielectric function; in units of eV."""
277 m = self._freq_w[0] * au_to_eV
278 assert isinstance(m, float)
279 return m
281 @property
282 def maxfreq(self) -> float:
283 """Highest frequency :attr:`~strongcoca.response.NumericDielectricFunction.frequencies`
284 for which there is a value of the dielectric function; in units of eV."""
285 m = self._freq_w[-1] * au_to_eV
286 assert isinstance(m, float)
287 return m
289 @property
290 def frequencies(self) -> np.ndarray:
291 """Frequencies corresponding to dielectric function
292 :attr:`~strongcoca.response.NumericDielectricFunction.dielectric_function` in units of eV;
293 array with shape {n}."""
294 return self._freq_w * au_to_eV # type: ignore
296 @property
297 def dielectric_function(self) -> np.ndarray:
298 """Complex unitless dielectric function at frequencies
299 :attr:`~strongcoca.response.MieGansResponse.frequencies`; array with shape {n}."""
300 return self._eps_w
303class DrudeDielectricFunction(AnalyticDielectricFunction):
305 r"""Objects of this class represent the Drude dielectric function, which is defined for the
306 entire positive frequency domain. The function is implemented as
308 .. math::
310 \varepsilon_r(\omega) = 1 - \frac{\omega_p^2}{\omega^2 + i\gamma\omega}
312 where the plasma frequency :math:`\omega_p` determines where the real part of the dielectric
313 function is zero and :math:`\gamma` is a broadening.
315 Parameters
316 ----------
317 plasma_frequency
318 The plasma frequency :math:`\omega_p`; units are eV by default, optionally atomic units
319 (see :attr:`units`).
320 broadening
321 Broadening of the dielectric function :math:`\gamma` units are eV by default, optionally
322 atomic units (see :attr:`units`).
323 name
324 Name of the dielectric function.
325 units
326 `eV` to specify :attr:`plasma_frequency` and :attr:`broadening` in eV, or `au` to specify
327 in atomic units.
329 This parameter determines whether conversion should be performed during
330 initialization and has no effect on instance methods and variables.
331 """
333 def __init__(self,
334 plasma_frequency: float,
335 broadening: float,
336 units: str = 'eV',
337 name: str = 'DrudeDielectricFunction') -> None:
339 logger.debug(f'Entering {self.__class__.__name__}.__init__')
340 super().__init__(name=name)
342 omega_pl = plasma_frequency
343 gamma = broadening
344 if units == 'eV':
345 omega_pl *= eV_to_au
346 gamma *= eV_to_au
347 elif units != 'au':
348 raise ValueError(f"units has to be 'eV' or 'au', not '{units}'")
350 self._omega_pl = omega_pl
351 self._gamma = gamma
353 def _eval_at(self, freq_w: np.ndarray) -> np.ndarray:
354 eps_w = 1 - self._omega_pl**2 / (freq_w**2 + 1j*self._gamma*freq_w)
355 return eps_w # type: ignore
357 def _repr_html_(self) -> str:
358 """HTML representation for Jupyter notebooks."""
359 rows = [
360 f'<tr><td style="text-align: left;">name</td><td>{self._name}</td></tr>',
361 f'<tr><td style="text-align: left;">plasma_frequency (eV)</td>'
362 f'<td>{self.plasma_frequency:.4f}</td></tr>',
363 f'<tr><td style="text-align: left;">broadening (eV)</td>'
364 f'<td>{self.broadening:.4f}</td></tr>',
365 ]
366 return (
367 f'<h4>{self.__class__.__name__}</h4>'
368 '<table>'
369 '<thead><tr>'
370 '<th style="text-align: left;">field</th>'
371 '<th style="text-align: left;">value</th>'
372 '</tr></thead>'
373 '<tbody>' + ''.join(rows) + '</tbody>'
374 '</table>'
375 )
377 def __str__(self) -> str:
378 fmt = ClassFormatter(self, pad=15)
379 fmt.append_class_name()
381 fmt.append_attr('name')
383 fmt.append_attr('plasma_frequency', unit='eV')
384 fmt.append_attr('broadening', unit='eV')
386 return fmt.to_string()
388 @property
389 def plasma_frequency(self) -> float:
390 """plasma frequency."""
391 return self._omega_pl * au_to_eV
393 @property
394 def broadening(self) -> float:
395 """Broadening parameter."""
396 return self._gamma * au_to_eV
399def read_dielec_function_file(dielectric_function_file: Path_str,
400 frequency_column: int = 0,
401 dielectric_function_real_column: int = 1,
402 dielectric_function_imag_column: int = 2,
403 units: str = 'eV',
404 name: str = 'NumericDielectricFunction from file',
405 **np_loadtxt_kwargs,
406 ) -> NumericDielectricFunction:
407 """Parses a textfile with columns describing the dielectric function of a material.
409 Parameters
410 ----------
411 dielectric_function_file
412 Path to the file containing the dielectic function.
413 frequency_column
414 Index of column containing the frequency vector in units of eV.
415 dielectric_function_real_column
416 Index of column containing the real part of the unitless dielectric function.
417 dielectric_function_imag_column
418 Index of column containing the imaginary part of the unitless dielectric function.
419 units
420 Units of the frequency column;
421 passed to :class:`strongcoca.response.NumericDielectricFunction`.
422 name
423 Passed to :class:`~strongcoca.response.NumericDielectricFunction`.
424 np_loadtxt_kwargs
425 Parameters passed to :func:`numpy.loadtxt`; useful parameters are `delimiter` or `comments`.
426 """
428 freq_w, df_real_w, df_imag_w = np.loadtxt(dielectric_function_file,
429 usecols=[frequency_column,
430 dielectric_function_real_column,
431 dielectric_function_imag_column],
432 unpack=True,
433 **np_loadtxt_kwargs)
434 df_w = df_real_w + 1j*df_imag_w
436 df = NumericDielectricFunction(freq_w, df_w, units=units, name=name)
438 return df
441def read_dielec_function_from_refractive_index_file(
442 dielectric_function_file: Path_str,
443 frequency_column: int = 0,
444 refractive_index_column: int = 1,
445 extinction_coefficient_column: int = 2,
446 units: str = 'eV',
447 name: str = 'NumericDielectricFunction from file',
448 **np_loadtxt_kwargs,
449 ) -> NumericDielectricFunction:
450 r"""Parses a textfile with columns describing the refractive index and extinction coefficient
451 of a material.
453 The refractive index :math:`n(\omega)` and extinction coefficient :math:`k(\omega)` are
454 related to the complex dielectric function :math:`\varepsilon_r(\omega)` by
456 .. math::
457 \varepsilon_r(\omega) = n^2(\omega) - k^2(\omega) + 2jn(\omega)k(\omega)
459 Parameters
460 ----------
461 dielectric_function_file
462 Path to the file containing the dielectic function.
463 frequency_column
464 Index of column containing the frequency vector in units of eV.
465 refractive_index_column
466 Index of column containing the unitless refractive index.
467 extinction_coefficient_column
468 Index of column containing the unitless extinction coefficient.
469 units
470 Units of the frequency column;
471 passed to :class:`~strongcoca.response.NumericDielectricFunction`.
472 name
473 Passed to :class:`~strongcoca.response.NumericDielectricFunction`.
474 np_loadtxt_kwargs
475 Parameters passed to :func:`numpy.loadtxt`; useful parameters are `delimiter` or `comments`.
476 """
478 freq_w, n_w, k_w = np.loadtxt(dielectric_function_file,
479 usecols=[frequency_column,
480 refractive_index_column,
481 extinction_coefficient_column],
482 unpack=True,
483 **np_loadtxt_kwargs)
484 df_w = n_w**2 - k_w**2 + 2j*n_w*k_w
486 df = NumericDielectricFunction(freq_w, df_w, units=units, name=name)
488 return df