Coverage for dynasor/post_processing/electron_scattering_factors.py: 100%

1import json 

2from re import search 

3from importlib.resources import files 

4from typing import Optional 

5from warnings import warn 

6 

7from numpy import exp, abs, pi 

8from pandas import DataFrame, concat 

9from .weights import Weights 

10 

11 

12class ElectronScatteringFactors(Weights): 

13 r"""This class generates sample weights corresponding to electron scattering factors. 

14 The scattering factors are parametrized as sums of exponentials 

15 of the following form 

16 

17 .. math:: 

18 

19 f(q) = \sum_{i=1}^k a_i \exp(-b_i s^2), 

20 

21 with :math:`s=\sin(\theta)/\lambda`, where and :math:`a_i` and :math:`b_i` 

22 are fitting parameters. Note that there are two parametrizations for the elemental 

23 scattering factors; one valid up to :math:`s = 2.0\,\mathrm{Å}^{-1}` 

24 and the other valid up to :math:`s = 6.0\,\mathrm{Å}^{-1}`. 

25 By default, the parametrization valid up to :math:`s = 6.0\,\mathrm{Å}^{-1}` 

26 is used. Ionic scattering factors only have a single parametrization. 

27 

28 The parametrizations are based on the following publications: 

29 

30 * Elemental scattering factors are based on Table 1 and Table 3 

31 from L.-M. Peng, G. Ren, S. L. Dudarev and M. J. Whelan, 

32 Acta Crystallographica Section A **52**, 257-276 (1996); 

33 `doi: 10.1107/S0108767395014371 <https://doi.org/10.1107/S0108767395014371>`_. 

34 * Ionic scattering factors are based on Table 1 from 

35 Lian-Mao. Peng, Acta Crystallographica Section A **54**, 481-485 (1998); 

36 `doi: 10.1107/S0108767398001901 <https://doi.org/10.1107/S0108767398001901>`_. 

37 

38 Parameters 

39 ---------- 

40 atom_types 

41 List of atomic species for which to retrieve scattering lengths. 

42 parametrization 

43 Parametrization to use. 

44 """ 

45 

46 def __init__( 

47 self, 

48 atom_types: list[str], 

49 parametrization: Optional[str] = 'smax6', 

50 ): 

51 self._parametrization = parametrization 

52 scattering_factors = self._read_scattering_factors(self._parametrization) 

53 # Select the relevant species 

54 scattering_factors = scattering_factors[ 

55 scattering_factors.index.isin(atom_types)] # Atom species is index 

56 self._scattering_factors = scattering_factors 

57 

58 # Check if any of the fetched form factors are missing, 

59 # indicating that it is missing in the experimental database. 

60 for s in atom_types: 

61 row = scattering_factors[scattering_factors.index == s] 

62 if row.empty: 

63 raise ValueError('Missing tabulated values ' 

64 f'for requested species {s}.') 

65 

66 weights_coh = scattering_factors.to_dict(orient='index') 

67 supports_currents = False 

68 super().__init__(weights_coh, None, supports_currents=supports_currents) 

69 

70 def get_weight_coh(self, atom_type, q_norm): 

71 """Get the coherent weight for a given atom type and q-vector norm.""" 

72 return self._compute_f(self._weights_coh[atom_type], 

73 q_norm, 

74 self._parametrization, 

75 self._get_charge(atom_type)) 

76 

77 @property 

78 def parameters(self) -> DataFrame: 

79 """Parametrization used to compute the coherent scattering factors 

80 :math:`f(q)` for the selected species. 

81 """ 

82 return self._scattering_factors 

83 

84 def _get_charge(self, atom_type: str) -> float: 

85 """ 

86 Extracts the ionic charge from the atom type. 

87 """ 

88 match = search(r'(\d+)?([+-])', atom_type) 

89 if match: 

90 magnitude = int(match.group(1)) if match.group(1) else 1 

91 charge = magnitude * (1 if match.group(2) == '+' else -1) 

92 return charge 

93 return None # If no charge is found 

94 

95 def _compute_f( 

96 self, 

97 coefficients: dict, 

98 q_norm: float, 

99 parametrization: str, 

100 charge: Optional[int] = None): 

101 r"""Compute electronic scattering factors :math:`f(q)`. 

102 The scattering factors are parametrized as a sum 

103 of exponentials of the form 

104 

105 .. math:: 

106 

107 f(q) = \sum_{i=1}^k a_i * exp(-b_i * s**2) 

108 

109 Ions are also offset by 

110 

111 .. math:: 

112 

113 (m_0 e^2) / (8 \pi^2 \hbar^2) C / s^2 \approx 0.023934 C / s^2 

114 

115 where :math:`C` is the ionic charge. 

116 

117 Parameters 

118 ---------- 

119 coefficients 

120 Parametrization parameters; read from the corresponding source file. 

121 q_norm 

122 The :math:`|q|`-value at which to evaluate the form factor. 

123 parametrization 

124 If the Peng parametrization valid up to :math:`s = 6.0\,\mathrm{Å}^{-1}` 

125 or :math:`s = 2.0\,\mathrm{Å}^{-1}` should be used. 

126 charge 

127 Integer representing the ionic charge. `None` if no charge to avoid division by zero. 

128 """ 

129 s_max = 6.0 if parametrization == 'smax6' else 2.0 

130 

131 s = q_norm / (4 * pi) # q in dynasor is q = 4 pi sin(theta) / lambda. 

132 s_squared = s*s # s = sin(theta) / lambda in the Waasmaier paper. 

133 if abs(s) > s_max: 

134 warn(f'Peng parametrization is not reliable for q' 

135 f' above {(s_max * 4 * pi):.2f} rad/Å' 

136 f' (corresponding to s={s_max} 1/Å).' 

137 ' Parametrisations for ions are accurate up to s=6.0 1/Å') 

138 nmax = 5 

139 

140 if charge is not None: 

141 if s_squared == 0: 

142 warn('Electron scattering factor for ionic species diverges at q=0; ' 

143 'returning infinity.') 

144 return float('inf') if charge > 0 else float('-inf') 

145 f = 0.023934 * charge / s_squared 

146 else: 

147 f = 0.0 

148 for i in range(1, nmax+1): 

149 f += coefficients[f'a{i}'] * exp(-coefficients[f'b{i}'] * s_squared) 

150 return f 

151 

152 def _read_scattering_factors(self, parametrization: str) -> DataFrame: 

153 r""" 

154 Extracts the parametrization for the form factors :math:`f(q)`, 

155 for both elemental systems and ionic species. 

156 

157 Parameters 

158 ---------- 

159 parametrization 

160 If the Peng parametrization valid up to :math:`s = 6.0\,\mathrm{Å}^{-1}` 

161 or :math:`s = 2.0\,\mathrm{Å}^{-1}` should be used. 

162 """ 

163 if parametrization == 'smax2': 

164 data_file = files(__package__) / \ 

165 'form-factors/electron-parameters-kmax2-peng-1996.json' 

166 elif parametrization == 'smax6': 

167 data_file = files(__package__) / \ 

168 'form-factors/electron-parameters-kmax6-peng-1996.json' 

169 else: 

170 raise ValueError(f'Unknown parametrization {parametrization}') 

171 

172 data_file_ions = files(__package__) / \ 

173 'form-factors/electron-parameters-ions-peng-1998.json' 

174 

175 frames = [] 

176 for df in [data_file, data_file_ions]: 

177 with open(df) as fp: 

178 coefficients = json.load(fp) 

179 frames.append(DataFrame.from_dict(coefficients)) 

180 

181 scattering_factors = concat(frames) 

182 scattering_factors.index.names = ['species'] 

183 return scattering_factors