import logging
import numpy as np
from scipy.integrate import trapezoid
from .base import BaseResponse
from .dielectric import DielectricFunction
from ..types import Array
from ..utilities import ClassFormatter
from ..units import au_to_A, A_to_au
from .utilities import NoArtificialBroadening
logger = logging.getLogger(__name__)
[docs]
class MieGansResponse(BaseResponse):
r"""Objects of this class hold representations of the Mie-Gans response of a
ellipsoid in vacuum.
For an ellipsoid with semi-axes :math:`a_x, a_y, a_z` in the Cartesian directions,
the Mie-Gans polarizability is diagonal
`[Sihvola, J. Nanomater. 2007, 045090] <https://doi.org/10.1155/2007/45090>`_
.. math::
\alpha_{\mu\mu}(\omega) = \varepsilon_0 V
\frac{\varepsilon_r(\omega) - 1}{1 + N_\mu (\varepsilon_r(\omega) - 1)}
where :math:`N_\mu` are depolarization factors
.. math::
N_\mu = \frac{a_x a_y a_z}{2} \int_0^\infty
\frac{\mathrm{d}s}{(s+a_\mu^2)\sqrt{(s+a_x^2)(s+a_y^2)(s+a_z^2)}}.
It is useful to keep in mind that in Hartree atomic units the vacuum permitivity
.. math::
\varepsilon_0 = \frac{1}{4\pi}.
Thus the expression is simplified
.. math::
\alpha_{\mu\mu}(\omega) = \frac{a_x a_y a_z}{3}
\frac{\varepsilon_r(\omega) - 1}{1 + N_\mu (\varepsilon_r(\omega) - 1)}
Parameters
----------
semiaxes
Semi-axes of the ellipsoid in Cartesian coordinates in units of Å; shape {3}.
Optionally radius of sphere in units of Å; shape {1}.
dielectric_function
Dielectric function, which can have an analytic form or be sampled on a frequency grid.
name
Name of response.
Raises
------
TypeError
If :attr:`semiaxes` is not scalar or three-dimensional array.
ValueError
If any component of :attr:`semiaxes` is zero or smaller.
Examples
--------
The following snippet illustrates how a Mie-Gans response object can be constructed.
Here, the analytic Drude dielectric function of the free-electron gas is used.
Instead, a numeric dielectric function can be used with the help of
:func:`~strongcoca.response.read_dielec_function_file`.
>>> import numpy as np
>>> from strongcoca.response import MieGansResponse, DrudeDielectricFunction
>>>
>>> semiaxes = 20
>>> omega_p = 3
>>> gamma = 0.5
>>> eps = DrudeDielectricFunction(omega_p, gamma)
>>> response = MieGansResponse(semiaxes, eps, name='Leaves')
>>> print(response)
--------------------- MieGansResponse ----------------------
name : Leaves
semiaxes (Å) : [20. 20. 20.]
depol. factors : [0.33335 0.33335 0.33335]
Dielectric function : DrudeDielectricFunction
Different shapes of the ellipsoid can be set by supplying three cartesian semiaxes. Here a
sphere, a prolate ellipsoid and an oblate ellipsoid are created.
>>> response = MieGansResponse([30, 30, 30], eps, name='Sphere')
>>> print(response)
--------------------- MieGansResponse ----------------------
name : Sphere
semiaxes (Å) : [30. 30. 30.]
depol. factors : [0.33335 0.33335 0.33335]
Dielectric function : DrudeDielectricFunction
>>> response = MieGansResponse([25, 25, 35], eps, name='Prolate ellipsoid')
>>> print(response)
--------------------- MieGansResponse ----------------------
name : Prolate ellipsoid
semiaxes (Å) : [25. 25. 35.]
depol. factors : [0.37561 0.37561 0.24881]
Dielectric function : DrudeDielectricFunction
>>> response = MieGansResponse([30, 30, 20], eps, name='Oblate ellipsoid')
>>> print(response)
--------------------- MieGansResponse ----------------------
name : Oblate ellipsoid
semiaxes (Å) : [30. 30. 20.]
depol. factors : [0.27705 0.27705 0.44592]
Dielectric function : DrudeDielectricFunction
"""
def __init__(self,
semiaxes: Array,
dielectric_function: DielectricFunction,
name: str = 'MieGansResponse') -> None:
logger.debug(f'Entering {self.__class__.__name__}.__init__')
super().__init__(pbc=False, broadening=NoArtificialBroadening(), name=name)
self._set_semiaxes(semiaxes)
self._factor_v = self._calc_depolarization_factors()
self._eps = dielectric_function
def __str__(self) -> str:
fmt = ClassFormatter(self, pad=20)
fmt.append_class_name()
fmt.append_attr('name')
fmt.append_attr('semiaxes', unit='Å')
fmt.append_attr('depol. factors', self.depolarization_factors)
fmt.append_attr('Dielectric function', self.dielectric_function.name)
return fmt.to_string()
@property
def semiaxes(self) -> np.ndarray:
"""Semi-axes of the ellipsoid in Cartesian coordinates in units of Å."""
return self._a_v * au_to_A # type: ignore
@property
def depolarization_factors(self) -> np.ndarray:
"""Depolarization factors."""
return self._factor_v
@property
def dielectric_function(self) -> DielectricFunction:
"""Dielectric function of the Mie-Gans ellipsoid."""
return self._eps
def _get_dynamic_polarizability(
self, frequencies: Array) -> np.ndarray:
freq_w = np.asarray(frequencies)
eps_w = self._eps._eval_at(freq_w)
epsm_w = np.ones_like(eps_w)
deps_w = eps_w - epsm_w
depolarization_factors = self.depolarization_factors
dm_vw = deps_w / (epsm_w + depolarization_factors[:, np.newaxis] * deps_w)
dm_vw *= 1 / 3 * np.prod(self._a_v)
# Transform into (wvv) matrix where the dm_wvv[i, :, :] is a diagonal matrix for any i
dm_wvv = np.einsum('xn,xy->nxy', dm_vw, np.eye(3), optimize=True)
return dm_wvv # type: ignore
def _get_dynamic_polarizability_imaginary_frequency(
self, frequencies: Array) -> np.ndarray:
raise NotImplementedError()
def _set_semiaxes(self, val: Array) -> None:
self._a_v = np.asarray(val, dtype=float) * A_to_au
if self._a_v.shape in [(), (1,)]:
self._a_v = self._a_v * np.ones(3)
elif self._a_v.shape != (3,):
raise TypeError('semiaxes must be scalar or three-dimensional array, '
f'has shape {self._a_v.shape}')
if not np.all(self._a_v > 0):
raise ValueError('All components of semiaxes must be strictly positive')
def _calc_depolarization_factors(self, mstep: int = 100, mmax: int = 100) -> np.ndarray:
"""Calculates the depolarization factors.
See https://doi.org/10.1155/2007/45090 Eq. (6).
"""
ds = np.min(self.semiaxes**2) / float(mstep)
smax = np.max(self.semiaxes**2) * mmax
s_t = np.arange(0, smax, ds)
sa2_vt = s_t + self.semiaxes[:, np.newaxis]**2
y_vt = 1.0 / (sa2_vt * np.sqrt(np.prod(sa2_vt, axis=0)))
# Integrate
depolarization_factors = trapezoid(y_vt, dx=ds, axis=-1)
# Add remainder
depolarization_factors += 2.0 / 3.0 * s_t[-1]**(-3.0 / 2.0)
depolarization_factors *= np.prod(self.semiaxes) / 2.0
return depolarization_factors # type: ignore
