Coverage for strongcoca / response / mie_gans.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 

11from .utilities import NoArtificialBroadening 

12 

13logger = logging.getLogger(__name__) 

14 

15 

16class MieGansResponse(BaseResponse): 

17 r"""Objects of this class hold representations of the Mie-Gans response of a 

18 ellipsoid in vacuum. 

19 

20 For an ellipsoid with semi-axes :math:`a_x, a_y, a_z` in the Cartesian directions, 

21 the Mie-Gans polarizability is diagonal 

22 `[Sihvola, J. Nanomater. 2007, 045090] <https://doi.org/10.1155/2007/45090>`_ 

23 

24 .. math:: 

25 

26 \alpha_{\mu\mu}(\omega) = \varepsilon_0 V 

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

28 

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

30 

31 .. math:: 

32 

33 N_\mu = \frac{a_x a_y a_z}{2} \int_0^\infty 

34 \frac{\mathrm{d}s}{(s+a_\mu^2)\sqrt{(s+a_x^2)(s+a_y^2)(s+a_z^2)}}. 

35 

36 It is useful to keep in mind that in Hartree atomic units the vacuum permitivity 

37 

38 .. math:: 

39 

40 \varepsilon_0 = \frac{1}{4\pi}. 

41 

42 Thus the expression is simplified 

43 

44 .. math:: 

45 

46 \alpha_{\mu\mu}(\omega) = \frac{a_x a_y a_z}{3} 

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

48 

49 

50 Parameters 

51 ---------- 

52 semiaxes 

53 Semi-axes of the ellipsoid in Cartesian coordinates in units of Å; shape {3}. 

54 

55 Optionally radius of sphere in units of Å; shape {1}. 

56 dielectric_function 

57 Dielectric function, which can have an analytic form or be sampled on a frequency grid. 

58 name 

59 Name of response. 

60 

61 Raises 

62 ------ 

63 TypeError 

64 If :attr:`semiaxes` is not scalar or three-dimensional array. 

65 ValueError 

66 If any component of :attr:`semiaxes` is zero or smaller. 

67 

68 Examples 

69 -------- 

70 

71 The following snippet illustrates how a Mie-Gans response object can be constructed. 

72 Here, the analytic Drude dielectric function of the free-electron gas is used. 

73 Instead, a numeric dielectric function can be used with the help of 

74 :func:`~strongcoca.response.read_dielec_function_file`. 

75 

76 >>> import numpy as np 

77 >>> from strongcoca.response import MieGansResponse, DrudeDielectricFunction 

78 >>> 

79 >>> semiaxes = 20 

80 >>> omega_p = 3 

81 >>> gamma = 0.5 

82 >>> eps = DrudeDielectricFunction(omega_p, gamma) 

83 >>> response = MieGansResponse(semiaxes, eps, name='Leaves') 

84 >>> print(response) 

85 --------------------- MieGansResponse ---------------------- 

86 name : Leaves 

87 semiaxes (Å) : [20. 20. 20.] 

88 Depolarization factors : [0.33335 0.33335 0.33335] 

89 Dielectric function : DrudeDielectricFunction 

90 

91 Different shapes of the ellipsoid can be set by supplying three cartesian semiaxes. Here a 

92 sphere, a prolate ellipsoid and an oblate ellipsoid are created. 

93 

94 >>> response = MieGansResponse([30, 30, 30], eps, name='Sphere') 

95 >>> print(response) 

96 --------------------- MieGansResponse ---------------------- 

97 name : Sphere 

98 semiaxes (Å) : [30. 30. 30.] 

99 Depolarization factors : [0.33335 0.33335 0.33335] 

100 Dielectric function : DrudeDielectricFunction 

101 

102 >>> response = MieGansResponse([25, 25, 35], eps, name='Prolate ellipsoid') 

103 >>> print(response) 

104 --------------------- MieGansResponse ---------------------- 

105 name : Prolate ellipsoid 

106 semiaxes (Å) : [25. 25. 35.] 

107 Depolarization factors : [0.37561 0.37561 0.24881] 

108 Dielectric function : DrudeDielectricFunction 

109 

110 >>> response = MieGansResponse([30, 30, 20], eps, name='Oblate ellipsoid') 

111 >>> print(response) 

112 --------------------- MieGansResponse ---------------------- 

113 name : Oblate ellipsoid 

114 semiaxes (Å) : [30. 30. 20.] 

115 Depolarization factors : [0.27705 0.27705 0.44592] 

116 Dielectric function : DrudeDielectricFunction 

117 """ 

118 def __init__(self, 

119 semiaxes: Array, 

120 dielectric_function: DielectricFunction, 

121 name: str = 'MieGansResponse') -> None: 

122 logger.debug(f'Entering {self.__class__.__name__}.__init__') 

123 super().__init__(pbc=False, broadening=NoArtificialBroadening(), name=name) 

124 

125 self._set_semiaxes(semiaxes) 

126 self._factor_v = self._calc_depolarization_factors() 

127 self._eps = dielectric_function 

128 

129 def __str__(self) -> str: 

130 fmt = ClassFormatter(self, pad=23) 

131 fmt.append_class_name() 

132 fmt.append_attr('name') 

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

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

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

136 return fmt.to_string() 

137 

138 @property 

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

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

141 return self._a_v * au_to_A # type: ignore 

142 

143 @property 

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

145 """Depolarization factors.""" 

146 return self._factor_v 

147 

148 @property 

149 def dielectric_function(self) -> DielectricFunction: 

150 """Dielectric function of the Mie-Gans ellipsoid.""" 

151 return self._eps 

152 

153 def _get_dynamic_polarizability( 

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

155 freq_w = np.asarray(frequencies) 

156 

157 eps_w = self._eps._eval_at(freq_w) 

158 epsm_w = np.ones_like(eps_w) 

159 deps_w = eps_w - epsm_w 

160 depolarization_factors = self.depolarization_factors 

161 

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

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

164 

165 # Transform into (wvv) matrix where the dm_wvv[i, :, :] is a diagonal matrix for any i 

166 dm_wvv = np.einsum('xn,xy->nxy', dm_vw, np.eye(3), optimize=True) 

167 return dm_wvv # type: ignore 

168 

169 def _get_dynamic_polarizability_imaginary_frequency( 

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

171 raise NotImplementedError() 

172 

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

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

175 

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

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

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

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

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

181 

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

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

184 

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

186 """Calculates the depolarization factors. 

187 

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

189 """ 

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

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

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

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

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

195 

196 # Integrate 

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

198 # Add remainder 

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

200 

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

202 return depolarization_factors # type: ignore