Coverage for strongcoca / response / mlwa.py: 100%
74 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-07-25 16:26 +0000
1import logging
3import numpy as np
4from scipy.integrate import trapezoid
6from .base import BaseResponse
7from .dielectric import DielectricFunction
8from ..types import Array
9from ..utilities import ClassFormatter
10from ..units import au_to_A, A_to_au, au_to_k
11from .utilities import NoArtificialBroadening
13logger = logging.getLogger(__name__)
16class MLWAResponse(BaseResponse):
17 r"""Objects of this class hold representations of the Mie-Gans response of an
18 ellipsoid in vacuum, corrected using the modified long wavelength approximation (MLWA).
20 For an ellipsoid with semi-axes :math:`a_x, a_y, a_z` in the Cartesian directions,
21 the quasistatic Mie-Gans polarizability is diagonal
23 .. math::
25 \alpha^{(0)}_{\mu\mu}(\omega) = \varepsilon_0 V
26 \frac{\varepsilon_r(\omega) - 1}{1 + N_\mu (\varepsilon_r(\omega) - 1)}
28 where :math:`N_\mu` are depolarization factors
30 .. math::
32 N_\mu = \frac{a_x a_y a_z}{2} \int_0^\infty
33 \frac{\mathrm{d}s}{(s+a_\mu^2)\sqrt{(s+a_x^2)(s+a_y^2)(s+a_z^2)}}.
35 In the MLWA, the quasistatic
36 polarizability is corrected by dynamic depolarization and radiation
37 damping. The diagonal polarizability becomes
39 .. math::
41 \alpha_{\mu\mu}^{\mathrm{MLWA}}(\omega) =
42 \frac{\alpha_{\mu\mu}^{(0)}(\omega)}
43 {1
44 - \dfrac{k^2}{a_\mu}\alpha_{\mu\mu}^{(0)}(\omega)
45 - i\dfrac{k^3}{6\pi\varepsilon_0}\alpha_{\mu\mu}^{(0)}(\omega)}
47 where :math:`k = \omega/c`.
49 In Hartree atomic units, :math:`\varepsilon_0 = 1/(4\pi)`, so the radiation
50 damping term becomes :math:`-i\,2k^3/3`.
52 Parameters
53 ----------
54 semiaxes
55 Semi-axes of the ellipsoid in Cartesian coordinates in units of Å; shape {3}.
56 Optionally radius of sphere in units of Å; shape {1}.
57 dielectric_function
58 Dielectric function, which can have an analytic form or be sampled on a frequency grid.
59 name
60 Name of response object.
62 Raises
63 ------
64 TypeError
65 If :attr:`semiaxes` is not scalar or a three-dimensional array.
66 ValueError
67 If any component of :attr:`semiaxes` is zero or smaller.
68 """
69 def __init__(self,
70 semiaxes: Array,
71 dielectric_function: DielectricFunction,
72 name: str = 'MLWAResponse') -> None:
73 logger.debug(f'Entering {self.__class__.__name__}.__init__')
74 super().__init__(pbc=False, broadening=NoArtificialBroadening(), name=name)
76 self._set_semiaxes(semiaxes)
77 self._factor_v = self._calc_depolarization_factors()
78 self._eps = dielectric_function
80 def __str__(self) -> str:
81 fmt = ClassFormatter(self, pad=23)
82 fmt.append_class_name()
83 fmt.append_attr('name')
84 fmt.append_attr('semiaxes', unit='Å')
85 fmt.append_attr('Depolarization factors', self.depolarization_factors)
86 fmt.append_attr('Dielectric function', self.dielectric_function.name)
87 return fmt.to_string()
89 @property
90 def semiaxes(self) -> np.ndarray:
91 """Semi-axes of the ellipsoid in Cartesian coordinates in units of Å."""
92 return self._a_v * au_to_A
94 @property
95 def depolarization_factors(self) -> np.ndarray:
96 """Depolarization factors."""
97 return self._factor_v
99 @property
100 def dielectric_function(self) -> DielectricFunction:
101 """Dielectric function of the ellipsoid."""
102 return self._eps
104 def _dynamic_correction(self, k: np.ndarray, polarizability: np.ndarray) -> np.ndarray:
105 return k**2*polarizability/self._a_v[:, np.newaxis] # type: ignore
107 def _radiative_correction(self, k: np.ndarray, polarizability: np.ndarray) -> np.ndarray:
108 return 1j*2.0/3.0*k**3*polarizability # type: ignore
110 def _mlwa_correction(self, k, polarizability: Array) -> np.ndarray:
111 polarizability = np.asarray(polarizability)
112 correction = (1 - self._radiative_correction(k, polarizability)
113 - self._dynamic_correction(k, polarizability))
114 return polarizability/correction # type: ignore
116 def _get_dynamic_polarizability(
117 self, frequencies: Array) -> np.ndarray:
119 freq_w = np.asarray(frequencies)
121 k_w = au_to_k * freq_w
122 k_w = k_w[np.newaxis, :]
123 eps_w = self._eps._eval_at(freq_w)
124 epsm_w = np.ones_like(eps_w)
125 deps_w = eps_w - epsm_w
126 depolarization_factors = self.depolarization_factors
128 dm_vw = deps_w / (epsm_w + depolarization_factors[:, np.newaxis] * deps_w)
129 dm_vw *= 1 / 3 * np.prod(self._a_v)
130 dm_mlwa = self._mlwa_correction(k_w, dm_vw)
131 dm_mlwav = np.einsum('xn,xy->nxy', dm_mlwa, np.eye(3), optimize=True)
133 return dm_mlwav # type: ignore
135 def _get_dynamic_polarizability_imaginary_frequency(
136 self, frequencies: Array) -> np.ndarray:
137 raise NotImplementedError()
139 def _set_semiaxes(self, val: Array) -> None:
140 self._a_v = np.asarray(val, dtype=float) * A_to_au
142 if self._a_v.shape in [(), (1,)]:
143 self._a_v = self._a_v * np.ones(3)
144 elif self._a_v.shape != (3,):
145 raise TypeError('semiaxes must be scalar or three-dimensional array, '
146 f'has shape {self._a_v.shape}')
148 if not np.all(self._a_v > 0):
149 raise ValueError('All components of semiaxes must be strictly positive')
151 def _calc_depolarization_factors(self, mstep: int = 100, mmax: int = 100) -> np.ndarray:
152 """Calculates the depolarization factors.
154 See https://doi.org/10.1155/2007/45090 Eq. (6).
155 """
156 ds = np.min(self.semiaxes**2) / float(mstep)
157 smax = np.max(self.semiaxes**2) * mmax
158 s_t = np.arange(0, smax, ds)
159 sa2_vt = s_t + self.semiaxes[:, np.newaxis]**2
160 y_vt = 1.0 / (sa2_vt * np.sqrt(np.prod(sa2_vt, axis=0)))
162 # Integrate
163 depolarization_factors = trapezoid(y_vt, dx=ds, axis=-1)
164 # Add remainder
165 depolarization_factors += 2.0 / 3.0 * s_t[-1]**(-3.0 / 2.0)
167 depolarization_factors *= np.prod(self.semiaxes) / 2.0
168 return depolarization_factors # type: ignore