Coverage for calorine/tools/refractive_index.py: 100%

1r""" 

2Conversions between the complex dielectric function :math:`\epsilon = \epsilon_1 + 

3i\epsilon_2` and the complex index of refraction :math:`N = n + i\kappa`, and between 

4:math:`\kappa` and the absorption coefficient :math:`\alpha`. 

5 

6Public entry points: :func:`get_refractive_index` (:math:`\epsilon \to n,\kappa`), 

7:func:`get_absorption_coefficient` (:math:`\kappa \to \alpha`), and 

8:func:`get_dielectric_function_from_refractive_index` (the inverse of the above, 

9accepting either :math:`n,\kappa` or :math:`n,\alpha`). All three operate on the same 

10DataFrame convention as :func:`~calorine.tools.get_dielectric_function` and 

11:func:`~calorine.tools.apply_kramers_kronig`: components are matched by column-name 

12suffix (e.g. ``epsilon_real_xx``/``epsilon_imag_xx``), so both the full Voigt tensor and 

13the bare unsuffixed (isotropic/molecular) case are handled uniformly. 

14""" 

15 

16import numpy as np 

17from pandas import DataFrame 

18from ase.units import _c as c_SI 

19 

20from .spectra import _rad_s_to_thz 

21 

22 

23def _eps_to_nk(eps_real, eps_imag): 

24 """Complex dielectric function -> (n, kappa). 

25 

26 n is taken non-negative (the physical branch of the square root); kappa then takes 

27 the sign of eps_imag, since eps_imag = 2*n*kappa with n >= 0. np.clip guards 

28 against tiny negative values from floating-point roundoff before taking the square 

29 root; mathematically ``eps_abs +/- eps_real`` is always non-negative. 

30 """ 

31 eps_abs = np.sqrt(eps_real ** 2 + eps_imag ** 2) 

32 n = np.sqrt(np.clip((eps_abs + eps_real) / 2, 0.0, None)) 

33 kappa = np.sqrt(np.clip((eps_abs - eps_real) / 2, 0.0, None)) 

34 kappa = np.copysign(kappa, eps_imag) 

35 return n, kappa 

36 

37 

38def get_refractive_index(df: DataFrame) -> DataFrame: 

39 r"""Compute the complex index of refraction from the complex dielectric function. 

40 

41 Given the complex dielectric function :math:`\epsilon = \epsilon_1 + i\epsilon_2`, 

42 the complex index of refraction :math:`N = n + i\kappa` is 

43 

44 .. math:: 

45 

46 n = \sqrt{\frac{|\epsilon| + \epsilon_1}{2}}, \qquad 

47 \kappa = \sqrt{\frac{|\epsilon| - \epsilon_1}{2}}, 

48 

49 where :math:`|\epsilon| = \sqrt{\epsilon_1^2 + \epsilon_2^2}`, taking the physical 

50 branch :math:`n \geq 0` and giving :math:`\kappa` the sign of :math:`\epsilon_2` 

51 (since :math:`\epsilon_2 = 2n\kappa` with :math:`n \geq 0`). For a passive medium 

52 at positive frequency :math:`\epsilon_2 \geq 0`, so :math:`\kappa \geq 0` as usual; 

53 a negative :math:`\epsilon_2` (e.g. numerical noise near zero) carries through as a 

54 negative :math:`\kappa` rather than being silently discarded. Both :math:`n` and 

55 :math:`\kappa` are dimensionless. 

56 

57 Each ``epsilon_real{suffix}``/``epsilon_imag{suffix}`` column pair in :attr:`df` 

58 (e.g. ``epsilon_real_xx``/``epsilon_imag_xx``, or the bare 

59 ``epsilon_real``/``epsilon_imag``) produces a corresponding 

60 ``refractive_index_real{suffix}``/``refractive_index_imag{suffix}`` column pair. 

61 

62 Parameters 

63 ---------- 

64 df 

65 DataFrame as returned by :func:`~calorine.tools.get_dielectric_function` (with 

66 ``return_real_part=True``) or :func:`~calorine.tools.apply_kramers_kronig`; 

67 must contain at least one ``epsilon_real{suffix}``/``epsilon_imag{suffix}`` 

68 column pair sharing the same suffix. 

69 

70 Returns 

71 ------- 

72 DataFrame 

73 Input DataFrame with additional ``refractive_index_real{suffix}`` (:math:`n`, 

74 dimensionless) and ``refractive_index_imag{suffix}`` (:math:`\kappa`, 

75 dimensionless) columns appended. 

76 

77 Raises 

78 ------ 

79 ValueError 

80 If :attr:`df` contains no ``epsilon_imag*`` columns, or an 

81 ``epsilon_imag{suffix}`` column without a matching ``epsilon_real{suffix}`` 

82 column (e.g. output of :func:`~calorine.tools.get_dielectric_function` called 

83 with ``return_real_part=False``). 

84 """ 

85 imag_cols = [c for c in df.columns if c.startswith('epsilon_imag')] 

86 if not imag_cols: 

87 raise ValueError( 

88 "df must contain at least one 'epsilon_imag*' column; " 

89 f'got columns: {list(df.columns)}' 

90 ) 

91 

92 missing = [ 

93 'epsilon_real' + col[len('epsilon_imag'):] 

94 for col in imag_cols 

95 if 'epsilon_real' + col[len('epsilon_imag'):] not in df.columns 

96 ] 

97 if missing: 

98 raise ValueError( 

99 f'df is missing columns {missing!r} required to compute the refractive ' 

100 'index; call get_dielectric_function with return_real_part=True, or run ' 

101 'apply_kramers_kronig first.' 

102 ) 

103 

104 df = df.copy() 

105 for col in imag_cols: 

106 suffix = col[len('epsilon_imag'):] 

107 n, kappa = _eps_to_nk(df['epsilon_real' + suffix].to_numpy(), df[col].to_numpy()) 

108 df['refractive_index_real' + suffix] = n 

109 df['refractive_index_imag' + suffix] = kappa 

110 return df 

111 

112 

113def get_absorption_coefficient(df: DataFrame) -> DataFrame: 

114 r"""Compute the absorption coefficient from the extinction coefficient. 

115 

116 .. math:: 

117 

118 \alpha(\omega) = \frac{2\kappa\omega}{c}, 

119 

120 where :math:`\kappa` is the extinction coefficient (imaginary part of the complex 

121 index of refraction) and :math:`\omega` is the angular frequency. With 

122 :math:`\omega` in rad/s and :math:`c` in m/s this gives :math:`\alpha` in 

123 m\ :sup:`-1`; :math:`\alpha` is returned here in cm\ :sup:`-1` (matching the 

124 ``wavenumber_invcm`` convention used elsewhere in this module, though note that 

125 :math:`\alpha` is a physically distinct quantity from the spectroscopic wavenumber 

126 :math:`\omega/2\pi c`): 

127 

128 .. math:: 

129 

130 \alpha\,[\mathrm{cm}^{-1}] = \frac{2\kappa\omega}{100\,c}. 

131 

132 Each ``refractive_index_imag{suffix}`` column in :attr:`df` (as returned by 

133 :func:`~calorine.tools.get_refractive_index`) produces a corresponding 

134 ``absorption_coefficient{suffix}`` column. As a convenience, if :attr:`df` lacks 

135 ``refractive_index_imag{suffix}`` but contains a matching ``epsilon_real{suffix}``/ 

136 ``epsilon_imag{suffix}`` pair, :math:`\kappa` is derived internally using the same 

137 relation as :func:`~calorine.tools.get_refractive_index` (without adding the 

138 corresponding ``refractive_index_real{suffix}``/``refractive_index_imag{suffix}`` 

139 columns to the output). 

140 

141 Parameters 

142 ---------- 

143 df 

144 DataFrame containing ``angular_frequency`` (THz) and at least one 

145 ``refractive_index_imag{suffix}`` column (as returned by 

146 :func:`~calorine.tools.get_refractive_index`), or, alternatively, a matching 

147 ``epsilon_real{suffix}``/``epsilon_imag{suffix}`` column pair (as returned by 

148 :func:`~calorine.tools.get_dielectric_function` with ``return_real_part=True``). 

149 

150 Returns 

151 ------- 

152 DataFrame 

153 Input DataFrame with additional ``absorption_coefficient{suffix}`` columns 

154 (cm\ :sup:`-1`) appended. 

155 

156 Raises 

157 ------ 

158 ValueError 

159 If :attr:`df` lacks ``angular_frequency``, or contains neither 

160 ``refractive_index_imag*`` columns nor a matching ``epsilon_real*``/ 

161 ``epsilon_imag*`` column pair for a given component. 

162 """ 

163 if 'angular_frequency' not in df.columns: 

164 raise ValueError("df must contain an 'angular_frequency' column") 

165 omega = df['angular_frequency'].to_numpy() / _rad_s_to_thz 

166 

167 kappas = {} 

168 for col in df.columns: 

169 if col.startswith('refractive_index_imag'): 

170 kappas[col[len('refractive_index_imag'):]] = df[col].to_numpy() 

171 

172 missing = [] 

173 for col in df.columns: 

174 if not col.startswith('epsilon_imag'): 

175 continue 

176 suffix = col[len('epsilon_imag'):] 

177 if suffix in kappas: 

178 continue 

179 real_col = 'epsilon_real' + suffix 

180 if real_col not in df.columns: 

181 missing.append(real_col) 

182 continue 

183 _, kappa = _eps_to_nk(df[real_col].to_numpy(), df[col].to_numpy()) 

184 kappas[suffix] = kappa 

185 

186 if missing: 

187 raise ValueError( 

188 f'df is missing columns {missing!r} required to compute the absorption ' 

189 "coefficient (no matching 'refractive_index_imag*' column and no matching " 

190 "'epsilon_real*' column to derive kappa from)." 

191 ) 

192 if not kappas: 

193 raise ValueError( 

194 "df must contain at least one 'refractive_index_imag*' column (as returned " 

195 "by get_refractive_index) or a matching 'epsilon_real*'/'epsilon_imag*' " 

196 'column pair' 

197 ) 

198 

199 df = df.copy() 

200 for suffix, kappa in kappas.items(): 

201 df['absorption_coefficient' + suffix] = 2.0 * kappa * omega / (100.0 * c_SI) 

202 return df 

203 

204 

205def get_dielectric_function_from_refractive_index(df: DataFrame) -> DataFrame: 

206 r"""Compute the complex dielectric function from the complex index of refraction. 

207 

208 Given the complex index of refraction :math:`N = n + i\kappa`, 

209 

210 .. math:: 

211 

212 \epsilon_1 = n^2 - \kappa^2, \qquad \epsilon_2 = 2 n \kappa. 

213 

214 This is the exact inverse of :func:`~calorine.tools.get_refractive_index`. 

215 

216 Each ``refractive_index_real{suffix}`` column (:math:`n`) must be paired with 

217 either a ``refractive_index_imag{suffix}`` column (:math:`\kappa`, as returned by 

218 :func:`~calorine.tools.get_refractive_index`) or an ``absorption_coefficient{suffix}`` 

219 column (:math:`\alpha`, in cm\ :sup:`-1`, as returned by 

220 :func:`~calorine.tools.get_absorption_coefficient`), from which :math:`\kappa` is 

221 recovered via the inverse of :math:`\alpha = 2\kappa\omega/(100 c)`, 

222 

223 .. math:: 

224 

225 \kappa = \frac{100\,c\,\alpha}{2\omega}, 

226 

227 which additionally requires the ``angular_frequency`` column. If both a 

228 ``refractive_index_imag{suffix}`` and an ``absorption_coefficient{suffix}`` column 

229 are present for the same suffix, the former takes precedence. 

230 

231 Parameters 

232 ---------- 

233 df 

234 DataFrame containing at least one ``refractive_index_real{suffix}`` column 

235 together with a matching ``refractive_index_imag{suffix}`` or 

236 ``absorption_coefficient{suffix}`` column (in the latter case 

237 ``angular_frequency`` in THz must also be present). 

238 

239 Returns 

240 ------- 

241 DataFrame 

242 Input DataFrame with additional ``epsilon_real{suffix}`` and 

243 ``epsilon_imag{suffix}`` columns appended. 

244 

245 Raises 

246 ------ 

247 ValueError 

248 If a ``refractive_index_real{suffix}`` column has no matching 

249 ``refractive_index_imag{suffix}`` or ``absorption_coefficient{suffix}`` 

250 column, or if the latter is used but ``angular_frequency`` is missing. 

251 """ 

252 real_cols = [c for c in df.columns if c.startswith('refractive_index_real')] 

253 if not real_cols: 

254 raise ValueError("df must contain at least one 'refractive_index_real*' column") 

255 

256 have_omega = 'angular_frequency' in df.columns 

257 omega = df['angular_frequency'].to_numpy() / _rad_s_to_thz if have_omega else None 

258 

259 kappas = {} 

260 missing = [] 

261 for col in real_cols: 

262 suffix = col[len('refractive_index_real'):] 

263 kappa_col = 'refractive_index_imag' + suffix 

264 alpha_col = 'absorption_coefficient' + suffix 

265 if kappa_col in df.columns: 

266 kappas[suffix] = df[kappa_col].to_numpy() 

267 elif alpha_col in df.columns: 

268 if not have_omega: 

269 missing.append( 

270 f"'angular_frequency' (required to convert {alpha_col!r} to kappa)" 

271 ) 

272 continue 

273 kappas[suffix] = 100.0 * c_SI * df[alpha_col].to_numpy() / (2.0 * omega) 

274 else: 

275 missing.append(f'{kappa_col!r} or {alpha_col!r}') 

276 

277 if missing: 

278 raise ValueError( 

279 f'df is missing required columns to invert refractive_index_real*: {missing}' 

280 ) 

281 

282 df = df.copy() 

283 for suffix, kappa in kappas.items(): 

284 n = df['refractive_index_real' + suffix].to_numpy() 

285 df['epsilon_real' + suffix] = n ** 2 - kappa ** 2 

286 df['epsilon_imag' + suffix] = 2.0 * n * kappa 

287 return df