Coverage for strongcoca / response / mlwa.py: 100%

1import logging 

2 

3import numpy as np 

4from scipy.integrate import trapezoid 

5 

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 

12 

13logger = logging.getLogger(__name__) 

14 

15 

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). 

19 

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 

22 

23 .. math:: 

24 

25 \alpha^{(0)}_{\mu\mu}(\omega) = \varepsilon_0 V 

26 \frac{\varepsilon_r(\omega) - 1}{1 + N_\mu (\varepsilon_r(\omega) - 1)} 

27 

28 where :math:`N_\mu` are depolarization factors 

29 

30 .. math:: 

31 

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)}}. 

34 

35 In the MLWA, the quasistatic 

36 polarizability is corrected by dynamic depolarization and radiation 

37 damping. The diagonal polarizability becomes 

38 

39 .. math:: 

40 

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)} 

46 

47 where :math:`k = \omega/c`. 

48 

49 In Hartree atomic units, :math:`\varepsilon_0 = 1/(4\pi)`, so the radiation 

50 damping term becomes :math:`-i\,2k^3/3`. 

51 

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. 

61 

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) 

75 self._set_semiaxes(semiaxes) 

76 self._factor_v = self._calc_depolarization_factors() 

77 self._eps = dielectric_function 

78 

79 def __str__(self) -> str: 

80 fmt = ClassFormatter(self, pad=23) 

81 fmt.append_class_name() 

82 fmt.append_attr('name') 

83 fmt.append_attr('semiaxes', unit='Å') 

84 fmt.append_attr('Depolarization factors', self.depolarization_factors) 

85 fmt.append_attr('Dielectric function', self.dielectric_function.name) 

86 return fmt.to_string() 

87 

88 @property 

89 def semiaxes(self) -> np.ndarray: 

90 """Semi-axes of the ellipsoid in Cartesian coordinates in units of Å.""" 

91 return self._a_v * au_to_A # type: ignore 

92 

93 @property 

94 def depolarization_factors(self) -> np.ndarray: 

95 """Depolarization factors.""" 

96 return self._factor_v 

97 

98 @property 

99 def dielectric_function(self) -> DielectricFunction: 

100 """Dielectric function of the ellipsoid.""" 

101 return self._eps 

102 

103 def _dynamic_correction(self, k: np.ndarray, polarizability: np.ndarray) -> np.ndarray: 

104 return k**2*polarizability/self.semiaxes[:, np.newaxis] 

105 

106 def _radiative_correction(self, k: np.ndarray, polarizability: np.ndarray) -> np.ndarray: 

107 return 1j*2.0/3.0*k**3*polarizability 

108 

109 def _mlwa_correction(self, k, polarizability: Array) -> np.ndarray: 

110 correction = 1 - self._radiative_correction(k, polarizability) \ 

111 - self._dynamic_correction(k, polarizability) 

112 return polarizability/correction 

113 

114 def _get_dynamic_polarizability( 

115 self, frequencies: Array) -> np.ndarray: 

116 

117 freq_w = np.asarray(frequencies) 

118 

119 k_w = au_to_k*freq_w 

120 k_w = k_w[np.newaxis, :] 

121 eps_w = self._eps._eval_at(freq_w) 

122 epsm_w = np.ones_like(eps_w) 

123 deps_w = eps_w - epsm_w 

124 depolarization_factors = self.depolarization_factors 

125 

126 dm_vw = deps_w / (epsm_w + depolarization_factors[:, np.newaxis] * deps_w) 

127 dm_vw *= 1 / 3 * np.prod(self._a_v) 

128 dm_mlwa = self._mlwa_correction(k_w, dm_vw) 

129 dm_mlwav = np.einsum('xn,xy->nxy', dm_mlwa, np.eye(3), optimize=True) 

130 

131 return dm_mlwav 

132 

133 def _get_dynamic_polarizability_imaginary_frequency( 

134 self, frequencies: Array) -> np.ndarray: 

135 raise NotImplementedError() 

136 

137 def _set_semiaxes(self, val: Array) -> None: 

138 self._a_v = np.asarray(val, dtype=float) * A_to_au 

139 

140 if self._a_v.shape in [(), (1,)]: 

141 self._a_v = self._a_v * np.ones(3) 

142 elif self._a_v.shape != (3,): 

143 raise TypeError('semiaxes must be scalar or three-dimensional array, ' 

144 f'has shape {self._a_v.shape}') 

145 

146 if not np.all(self._a_v > 0): 

147 raise ValueError('All components of semiaxes must be strictly positive') 

148 

149 def _calc_depolarization_factors(self, mstep: int = 100, mmax: int = 100) -> np.ndarray: 

150 """Calculates the depolarization factors. 

151 

152 See https://doi.org/10.1155/2007/45090 Eq. (6). 

153 """ 

154 ds = np.min(self.semiaxes**2) / float(mstep) 

155 smax = np.max(self.semiaxes**2) * mmax 

156 s_t = np.arange(0, smax, ds) 

157 sa2_vt = s_t + self.semiaxes[:, np.newaxis]**2 

158 y_vt = 1.0 / (sa2_vt * np.sqrt(np.prod(sa2_vt, axis=0))) 

159 

160 # Integrate 

161 depolarization_factors = trapezoid(y_vt, dx=ds, axis=-1) 

162 # Add remainder 

163 depolarization_factors += 2.0 / 3.0 * s_t[-1]**(-3.0 / 2.0) 

164 

165 depolarization_factors *= np.prod(self.semiaxes) / 2.0 

166 return depolarization_factors # type: ignore