Coverage for /builds/kinetik161/ase/ase/utils/xrdebye.py: 82.54%

1# flake8: noqa 

2"""Definition of the XrDebye class. 

3 

4This module defines the XrDebye class for calculation 

5of X-ray scattering properties from atomic cluster 

6using Debye formula. 

7Also contains routine for calculation of atomic form factors and 

8X-ray wavelength dict. 

9""" 

10 

11from math import acos, cos, exp, pi, sin, sqrt 

12 

13import numpy as np 

14 

15from ase.data import atomic_numbers 

16 

17# Table (1) of 

18# D. WAASMAIER AND A. KIRFEL, Acta Cryst. (1995). A51, 416-431 

19waasmaier = { 

20 # a1 b1 a2 b2 a3 b3 a4 b4 a5 b5 c 

21 'C': [2.657506, 14.780758, 1.078079, 0.776775, 1.490909, 42.086843, -4.241070, -0.000294, 0.713791, 0.239535, 4.297983], 

22 'N': [11.893780, 0.000158, 3.277479, 10.232723, 1.858092, 30.344690, 0.858927, 0.656065, 0.912985, 0.217287, -11.804902], 

23 'O': [2.960427, 14.182259, 2.5088111, 5.936858, 0.637053, 0.112726, 0.722838, 34.958481, 1.142756, 0.390240, 0.027014], 

24 'P': [1.950541, 0.908139, 4.146930, 27.044953, 1.494560, 0.071280, 1.522042, 67.520190, 5.729711, 1.981173, 0.155233], 

25 'S': [6.372157, 1.514347, 5.154568, 22.092528, 1.473732, 0.061373, 1.635073, 55.445176, 1.209372, 0.646925, 0.154722], 

26 'Cl': [1.446071, 0.052357, 6.870609, 1.193165, 6.151801, 18.343416, 1.750347, 46.398394, 0.634168, 0.401005, 0.146773], 

27 'Ni': [13.521865, 4.077277, 6.947285, 0.286763, 3.866028, 14.622634, 2.135900, 71.966078, 4.284731, 0.004437, -2.762697], 

28 'Cu': [14.014192, 3.738280, 4.784577, 0.003744, 5.056806, 13.034982, 1.457971, 72.554793, 6.932996, 0.265666, -3.774477], 

29 'Pd': [6.121511, 0.062549, 4.784063, 0.784031, 16.631683, 8.751391, 4.318258, 34.489983, 13.246773, 0.784031, 0.883099], 

30 'Ag': [6.073874, 0.055333, 17.155437, 7.896512, 4.173344, 28.443739, 0.852238, 110.376108, 17.988685, 0.716809, 0.756603], 

31 'Pt': [31.273891, 1.316992, 18.445441, 8.797154, 17.063745, 0.124741, 5.555933, 40.177994, 1.575270, 1.316997, 4.050394], 

32 'Au': [16.777389, 0.122737, 19.317156, 8.621570, 32.979682, 1.256902, 5.595453, 38.008821, 10.576854, 0.000601, -6.279078], 

33} 

34 

35wavelengths = { 

36 'CuKa1': 1.5405981, 

37 'CuKa2': 1.54443, 

38 'CuKb1': 1.39225, 

39 'WLa1': 1.47642, 

40 'WLa2': 1.48748 

41} 

42 

43 

44class XrDebye: 

45 """ 

46 Class for calculation of XRD or SAXS patterns. 

47 """ 

48 

49 def __init__(self, atoms, wavelength, damping=0.04, 

50 method='Iwasa', alpha=1.01, warn=True): 

51 """ 

52 Initilize the calculation of X-ray diffraction patterns 

53 

54 Parameters: 

55 

56 atoms: ase.Atoms 

57 atoms object for which calculation will be performed. 

58 

59 wavelength: float, Angstrom 

60 X-ray wavelength in Angstrom. Used for XRD and to setup dumpings. 

61 

62 damping : float, Angstrom**2 

63 thermal damping factor parameter (B-factor). 

64 

65 method: {'Iwasa'} 

66 method of calculation (damping and atomic factors affected). 

67 

68 If set to 'Iwasa' than angular damping and q-dependence of 

69 atomic factors are used. 

70 

71 For any other string there will be only thermal damping 

72 and constant atomic factors (`f_a(q) = Z_a`). 

73 

74 alpha: float 

75 parameter for angular damping of scattering intensity. 

76 Close to 1.0 for unplorized beam. 

77 

78 warn: boolean 

79 flag to show warning if atomic factor can't be calculated 

80 """ 

81 self.wavelength = wavelength 

82 self.damping = damping 

83 self.mode = '' 

84 self.method = method 

85 self.alpha = alpha 

86 self.warn = warn 

87 

88 self.twotheta_list = [] 

89 self.q_list = [] 

90 self.intensity_list = [] 

91 

92 self.atoms = atoms 

93 # TODO: setup atomic form factors if method != 'Iwasa' 

94 

95 def set_damping(self, damping): 

96 """ set B-factor for thermal damping """ 

97 self.damping = damping 

98 

99 def get(self, s): 

100 r"""Get the powder x-ray (XRD) scattering intensity 

101 using the Debye-Formula at single point. 

102 

103 Parameters: 

104 

105 s: float, in inverse Angstrom 

106 scattering vector value (`s = q / 2\pi`). 

107 

108 Returns: 

109 Intensity at given scattering vector `s`. 

110 """ 

111 

112 pre = exp(-self.damping * s**2 / 2) 

113 

114 if self.method == 'Iwasa': 

115 sinth = self.wavelength * s / 2. 

116 positive = 1. - sinth**2 

117 if positive < 0: 

118 positive = 0 

119 costh = sqrt(positive) 

120 cos2th = cos(2. * acos(costh)) 

121 pre *= costh / (1. + self.alpha * cos2th**2) 

122 

123 f = {} 

124 

125 def atomic(symbol): 

126 """ 

127 get atomic factor, using cache. 

128 """ 

129 if symbol not in f: 

130 if self.method == 'Iwasa': 

131 f[symbol] = self.get_waasmaier(symbol, s) 

132 else: 

133 f[symbol] = atomic_numbers[symbol] 

134 return f[symbol] 

135 

136 I = 0. 

137 fa = [] # atomic factors list 

138 for a in self.atoms: 

139 fa.append(atomic(a.symbol)) 

140 

141 pos = self.atoms.get_positions() # positions of atoms 

142 fa = np.array(fa) # atomic factors array 

143 

144 for i in range(len(self.atoms)): 

145 vr = pos - pos[i] 

146 I += np.sum(fa[i] * fa * np.sinc(2 * s * 

147 np.sqrt(np.sum(vr * vr, axis=1)))) 

148 

149 return pre * I 

150 

151 def get_waasmaier(self, symbol, s): 

152 r"""Scattering factor for free atoms. 

153 

154 Parameters: 

155 

156 symbol: string 

157 atom element symbol. 

158 

159 s: float, in inverse Angstrom 

160 scattering vector value (`s = q / 2\pi`). 

161 

162 Returns: 

163 Intensity at given scattering vector `s`. 

164 

165 Note: 

166 for hydrogen will be returned zero value.""" 

167 if symbol == 'H': 

168 # XXXX implement analytical H 

169 return 0 

170 elif symbol in waasmaier: 

171 abc = waasmaier[symbol] 

172 f = abc[10] 

173 s2 = s * s 

174 for i in range(5): 

175 f += abc[2 * i] * exp(-abc[2 * i + 1] * s2) 

176 return f 

177 if self.warn: 

178 print('<xrdebye::get_atomic> Element', symbol, 'not available') 

179 return 0 

180 

181 def calc_pattern(self, x=None, mode='XRD', verbose=False): 

182 r""" 

183 Calculate X-ray diffraction pattern or 

184 small angle X-ray scattering pattern. 

185 

186 Parameters: 

187 

188 x: float array 

189 points where intensity will be calculated. 

190 XRD - 2theta values, in degrees; 

191 SAXS - q values in 1/A 

192 (`q = 2 \pi \cdot s = 4 \pi \sin( \theta) / \lambda`). 

193 If ``x`` is ``None`` then default values will be used. 

194 

195 mode: {'XRD', 'SAXS'} 

196 the mode of calculation: X-ray diffraction (XRD) or 

197 small-angle scattering (SAXS). 

198 

199 Returns: 

200 list of intensities calculated for values given in ``x``. 

201 """ 

202 self.mode = mode.upper() 

203 assert (mode in ['XRD', 'SAXS']) 

204 

205 result = [] 

206 if mode == 'XRD': 

207 if x is None: 

208 self.twotheta_list = np.linspace(15, 55, 100) 

209 else: 

210 self.twotheta_list = x 

211 self.q_list = [] 

212 if verbose: 

213 print('#2theta\tIntensity') 

214 for twotheta in self.twotheta_list: 

215 s = 2 * sin(twotheta * pi / 180 / 2.0) / self.wavelength 

216 result.append(self.get(s)) 

217 if verbose: 

218 print(f'{twotheta:.3f}\t{result[-1]:f}') 

219 elif mode == 'SAXS': 

220 if x is None: 

221 self.twotheta_list = np.logspace(-3, -0.3, 100) 

222 else: 

223 self.q_list = x 

224 self.twotheta_list = [] 

225 if verbose: 

226 print('#q\tIntensity') 

227 for q in self.q_list: 

228 s = q / (2 * pi) 

229 result.append(self.get(s)) 

230 if verbose: 

231 print(f'{q:.4f}\t{result[-1]:f}') 

232 self.intensity_list = np.array(result) 

233 return self.intensity_list 

234 

235 def write_pattern(self, filename): 

236 """ Save calculated data to file specified by ``filename`` string.""" 

237 with open(filename, 'w') as fd: 

238 self._write_pattern(fd) 

239 

240 def _write_pattern(self, fd): 

241 fd.write('# Wavelength = %f\n' % self.wavelength) 

242 if self.mode == 'XRD': 

243 x, y = self.twotheta_list, self.intensity_list 

244 fd.write('# 2theta \t Intesity\n') 

245 elif self.mode == 'SAXS': 

246 x, y = self.q_list, self.intensity_list 

247 fd.write('# q(1/A)\tIntesity\n') 

248 else: 

249 raise Exception('No data available, call calc_pattern() first.') 

250 

251 for i in range(len(x)): 

252 fd.write(f' {x[i]:f}\t{y[i]:f}\n') 

253 

254 def plot_pattern(self, filename=None, show=False, ax=None): 

255 """ Plot XRD or SAXS depending on filled data 

256 

257 Uses Matplotlib to plot pattern. Use *show=True* to 

258 show the figure and *filename='abc.png'* or 

259 *filename='abc.eps'* to save the figure to a file. 

260 

261 Returns: 

262 ``matplotlib.axes.Axes`` object.""" 

263 

264 import matplotlib.pyplot as plt 

265 

266 if ax is None: 

267 plt.clf() # clear figure 

268 ax = plt.gca() 

269 

270 if self.mode == 'XRD': 

271 x, y = np.array(self.twotheta_list), np.array(self.intensity_list) 

272 ax.plot(x, y / np.max(y), '.-') 

273 ax.set_xlabel('2$\\theta$') 

274 ax.set_ylabel('Intensity') 

275 elif self.mode == 'SAXS': 

276 x, y = np.array(self.q_list), np.array(self.intensity_list) 

277 ax.loglog(x, y / np.max(y), '.-') 

278 ax.set_xlabel('q, 1/Angstr.') 

279 ax.set_ylabel('Intensity') 

280 else: 

281 raise Exception('No data available, call calc_pattern() first') 

282 

283 if show: 

284 plt.show() 

285 if filename is not None: 

286 fig = ax.get_figure() 

287 fig.savefig(filename) 

288 

289 return ax