Coverage for strongcoca/response/dielectric.py: 100%
124 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-26 18:44 +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 @property
61 def name(self) -> str:
62 """Name of dielectric function."""
63 return self._name
66class AnalyticDielectricFunction(DielectricFunction):
68 """Objects of this class hold different representations of dielectric functions that
69 have analytic form. The dielectric function can thus be evaluated at any frequency
70 without interpolation.
72 Parameters
73 ----------
74 name
75 Name of dielectric function.
76 """
78 def __init__(self,
79 name: str = 'AnalyticDielectricFunction') -> None:
81 logger.debug(f'Entering {self.__class__.__name__}.__init__')
82 super().__init__(name=name)
85class NumericDielectricFunction(DielectricFunction):
87 """Objects of this class hold different representations of dielectric functions sampled
88 on a discrete grid of frequencies. The dielectric function can be evaluated at any frequency
89 within the supplied grid by linear interpolation.
91 Parameters
92 ----------
93 frequencies
94 List of frequencies at which the dielectric function :attr:`dielectric_function`
95 is sampled; units are eV by default, optionally atomic units
96 (see :attr:`units`).
98 The list of frequencies and the corresponding dielectric function are sorted if the former
99 is not monotonically increasing.
100 dielectric_function
101 Complex dielectric function sampled at the frequencies :attr:`frequencies`.
102 units
103 `eV` to specify :attr:`freq_w` in eV, or `au` to specify in atomic units.
105 This parameter determines whether conversion should be performed during
106 initialization and has no effect on instance methods and variables.
107 name
108 Name of dielectric function.
110 Raises
111 ------
112 ValueError
113 If :attr:`frequencies` or :attr:`dielectric_function` have zero length.
114 ValueError
115 If there is a mismatch in the length of :attr:`frequencies` and :attr:`dielectric_function`.
116 ValueError
117 If :attr:`units` is not one of the allowed values `eV` or `au`.
119 Examples
120 --------
122 The following snippet illustrates how a NumericDielectricFunction object can be constructed.
123 Here, the analytic Drude dielectric function of the free-electron gas is used.
124 To construct a dielectric function object from a numeric dielectric function the function
125 :func:`~strongcoca.response.read_dielec_function_file` can be used.
127 >>> import numpy as np
128 >>> from strongcoca.response.dielectric import NumericDielectricFunction
129 >>>
130 >>> omega_p = 3
131 >>> gamma = 0.5
132 >>> freq_w = np.linspace(0.1, 10)
133 >>> eps_w = 1 - omega_p**2 / (freq_w**2 + 1j*gamma*freq_w)
134 >>> dielec = NumericDielectricFunction(freq_w, eps_w, name='Drude-function')
135 >>> print(dielec)
136 ---------------- NumericDielectricFunction -----------------
137 name : Drude-function
138 Frequency grid : 50 points between 0.10 and 10.00eV
139 frequencies (eV) : [ 0.1 0.30204 0.50408 ...
140 9.59592 9.79796 10. ]
141 dielectric_function : [-33.61538+1.73077e+02j
142 -25.37528+4.36618e+01j
143 -16.85366+1.77091e+01j ...
144 0.90253+5.07897e-03j
145 0.90649+4.77173e-03j
146 0.91022+4.48878e-03j]
147 >>> dense_freq_w = np.linspace(3, 4)
148 >>> print(dielec(dense_freq_w))
149 [0.02425453+0.16310249j 0.0367996 +0.15996402j 0.04934466+0.15682555j
150 0.06188972+0.15368708j 0.07443478+0.15054861j 0.08697984+0.14741014j
151 ...
152 """
154 def __init__(self,
155 frequencies: Union[Array, float],
156 dielectric_function: Union[ComplexArray, complex],
157 units: str = 'eV',
158 name: str = 'NumericDielectricFunction') -> None:
160 logger.debug(f'Entering {self.__class__.__name__}.__init__')
161 super().__init__(name=name)
163 if isinstance(frequencies, float):
164 frequencies = [frequencies]
166 if isinstance(dielectric_function, complex) or isinstance(dielectric_function, float):
167 dielectric_function = [dielectric_function]
169 freq_w = np.array(frequencies)
170 eps_w = np.array(dielectric_function, dtype=complex)
172 if units == 'eV':
173 freq_w = freq_w * eV_to_au
174 elif units != 'au':
175 raise ValueError(f"units has to be 'eV' or 'au', not '{units}'")
177 if len(freq_w) != len(eps_w):
178 raise ValueError(f'length of frequencies {len(freq_w)} must match length of '
179 f'dielectric function {len(eps_w)}')
181 if len(freq_w) == 0:
182 raise ValueError('length of frequencies and dielectric function must not be zero')
184 # np.diff works if input has at least length 2
185 if len(freq_w) > 1 and not np.all(np.diff(freq_w) > 0):
186 sort_w = np.argsort(freq_w)
187 freq_w = freq_w[sort_w]
188 eps_w = eps_w[sort_w]
190 self._freq_w = freq_w
191 self._eps_w = eps_w
193 def _eval_at(self, freq_w: np.ndarray) -> np.ndarray:
194 """Returns the dielectric function at the given frequencies by interpolating the numeric
195 dielectric function.
197 Parameters
198 ----------
199 frequencies
200 List of frequencies at which the dielectric function is returned;
201 in units of eV.
203 Raises
204 ------
205 ValueError
206 If any component of :attr:`frequencies` is outside of the interval of the numeric
207 dielectric function.
208 """
209 if np.min(freq_w) < self._freq_w[0]:
210 raise ValueError('Minimum of array frequencies is outside the lower bound of the'
211 'numeric dielectric function')
212 if np.max(freq_w) > self._freq_w[-1]:
213 raise ValueError('Maximum of array frequencies is outside the upper bound of the'
214 'numeric dielectric function')
216 eps_w = np.interp(freq_w, self._freq_w, self._eps_w)
218 return eps_w # type: ignore
220 def __str__(self) -> str:
221 fmt = ClassFormatter(self, pad=20)
222 fmt.append_class_name()
224 fmt.append_attr('name')
226 fmt.append_attr('Frequency grid', f'{self.n} points between {self.minfreq:.2f} and '
227 f'{self.maxfreq:.2f}eV')
229 fmt.append_attr('frequencies', unit='eV')
230 fmt.append_attr('dielectric_function')
232 return fmt.to_string()
234 @property
235 def n(self) -> int:
236 """number of samples for the dielectric function
237 :attr:`~strongcoca.response.NumericDielectricFunction.dielectric_function`."""
238 return len(self._eps_w)
240 @property
241 def minfreq(self) -> float:
242 """Lowest frequency :attr:`~strongcoca.response.NumericDielectricFunction.frequencies`
243 for which there is a value of the dielectric function; in units of eV."""
244 m = self._freq_w[0] * au_to_eV
245 assert isinstance(m, float)
246 return m
248 @property
249 def maxfreq(self) -> float:
250 """Highest frequency :attr:`~strongcoca.response.NumericDielectricFunction.frequencies`
251 for which there is a value of the dielectric function; in units of eV."""
252 m = self._freq_w[-1] * au_to_eV
253 assert isinstance(m, float)
254 return m
256 @property
257 def frequencies(self) -> np.ndarray:
258 """Frequencies corresponding to dielectric function
259 :attr:`~strongcoca.response.NumericDielectricFunction.dielectric_function` in units of eV;
260 array with shape {n}."""
261 return self._freq_w * au_to_eV # type: ignore
263 @property
264 def dielectric_function(self) -> np.ndarray:
265 """Complex unitless dielectric function at frequencies
266 :attr:`~strongcoca.response.MieGansResponse.frequencies`; array with shape {n}."""
267 return self._eps_w
270class DrudeDielectricFunction(AnalyticDielectricFunction):
272 r"""Objects of this class represent the Drude dielectric function, which is defined for the
273 entire positive frequency domain. The function is implemented as
275 .. math::
277 \varepsilon_r(\omega) = 1 - \frac{\omega_p^2}{\omega^2 + i\gamma\omega}
279 where the plasma frequency :math:`\omega_p` determines where the real part of the dielectric
280 function is zero and :math:`\gamma` is a broadening.
282 Parameters
283 ----------
284 plasma_frequency
285 The plasma frequency :math:`\omega_p`; units are eV by default, optionally atomic units
286 (see :attr:`units`).
287 broadening
288 Broadening of the dielectric function :math:`\gamma` units are eV by default, optionally
289 atomic units (see :attr:`units`).
290 name
291 Name of the dielectric function.
292 units
293 `eV` to specify :attr:`plasma_frequency` and :attr:`broadening` in eV, or `au` to specify
294 in atomic units.
296 This parameter determines whether conversion should be performed during
297 initialization and has no effect on instance methods and variables.
298 """
300 def __init__(self,
301 plasma_frequency: float,
302 broadening: float,
303 units: str = 'eV',
304 name: str = 'DrudeDielectricFunction') -> None:
306 logger.debug(f'Entering {self.__class__.__name__}.__init__')
307 super().__init__(name=name)
309 omega_pl = plasma_frequency
310 gamma = broadening
311 if units == 'eV':
312 omega_pl *= eV_to_au
313 gamma *= eV_to_au
314 elif units != 'au':
315 raise ValueError(f"units has to be 'eV' or 'au', not '{units}'")
317 self._omega_pl = omega_pl
318 self._gamma = gamma
320 def _eval_at(self, freq_w: np.ndarray) -> np.ndarray:
321 eps_w = 1 - self._omega_pl**2 / (freq_w**2 + 1j*self._gamma*freq_w)
322 return eps_w # type: ignore
324 def __str__(self) -> str:
325 fmt = ClassFormatter(self, pad=15)
326 fmt.append_class_name()
328 fmt.append_attr('name')
330 fmt.append_attr('plasma_frequency', unit='eV')
331 fmt.append_attr('broadening', unit='eV')
333 return fmt.to_string()
335 @property
336 def plasma_frequency(self) -> float:
337 """plasma frequency."""
338 return self._omega_pl * au_to_eV
340 @property
341 def broadening(self) -> float:
342 """Broadening parameter."""
343 return self._gamma * au_to_eV
346def read_dielec_function_file(dielectric_function_file: Path_str,
347 frequency_column: int = 0,
348 dielectric_function_real_column: int = 1,
349 dielectric_function_imag_column: int = 2,
350 units: str = 'eV',
351 name: str = 'NumericDielectricFunction from file',
352 **np_loadtxt_kwargs,
353 ) -> NumericDielectricFunction:
354 """Parses a textfile with columns describing the dielectric function of a material.
356 Parameters
357 ----------
358 dielectric_function_file
359 Path to the file containing the dielectic function.
360 frequency_column
361 Index of column containing the frequency vector in units of eV.
362 dielectric_function_real_column
363 Index of column containing the real part of the unitless dielectric function.
364 dielectric_function_imag_column
365 Index of column containing the imaginary part of the unitless dielectric function.
366 units
367 Units of the frequency column;
368 passed to :class:`strongcoca.response.NumericDielectricFunction`.
369 name
370 Passed to :class:`~strongcoca.response.NumericDielectricFunction`.
371 np_loadtxt_kwargs
372 Parameters passed to :func:`numpy.loadtxt`; useful parameters are `delimiter` or `comments`.
373 """
375 freq_w, df_real_w, df_imag_w = np.loadtxt(dielectric_function_file,
376 usecols=[frequency_column,
377 dielectric_function_real_column,
378 dielectric_function_imag_column],
379 unpack=True,
380 **np_loadtxt_kwargs)
381 df_w = df_real_w + 1j*df_imag_w
383 df = NumericDielectricFunction(freq_w, df_w, units=units, name=name)
385 return df
388def read_dielec_function_from_refractive_index_file(
389 dielectric_function_file: Path_str,
390 frequency_column: int = 0,
391 refractive_index_column: int = 1,
392 extinction_coefficient_column: int = 2,
393 units: str = 'eV',
394 name: str = 'NumericDielectricFunction from file',
395 **np_loadtxt_kwargs,
396 ) -> NumericDielectricFunction:
397 r"""Parses a textfile with columns describing the refractive index and extinction coefficient
398 of a material.
400 The refractive index :math:`n(\omega)` and extinction coefficient :math:`k(\omega)` are
401 related to the complex dielectric function :math:`\varepsilon_r(\omega)` by
403 .. math::
404 \varepsilon_r(\omega) = n^2(\omega) - k^2(\omega) + 2jn(\omega)k(\omega)
406 Parameters
407 ----------
408 dielectric_function_file
409 Path to the file containing the dielectic function.
410 frequency_column
411 Index of column containing the frequency vector in units of eV.
412 refractive_index_column
413 Index of column containing the unitless refractive index.
414 extinction_coefficient_column
415 Index of column containing the unitless extinction coefficient.
416 units
417 Units of the frequency column;
418 passed to :class:`~strongcoca.response.NumericDielectricFunction`.
419 name
420 Passed to :class:`~strongcoca.response.NumericDielectricFunction`.
421 np_loadtxt_kwargs
422 Parameters passed to :func:`numpy.loadtxt`; useful parameters are `delimiter` or `comments`.
423 """
425 freq_w, n_w, k_w = np.loadtxt(dielectric_function_file,
426 usecols=[frequency_column,
427 refractive_index_column,
428 extinction_coefficient_column],
429 unpack=True,
430 **np_loadtxt_kwargs)
431 df_w = n_w**2 - k_w**2 + 2j*n_w*k_w
433 df = NumericDielectricFunction(freq_w, df_w, units=units, name=name)
435 return df