Coverage for /builds/ase/ase/ase/utils/xrdebye.py: 82.40%

1# fmt: off 

2 

3# flake8: noqa 

4"""Definition of the XrDebye class. 

5 

6This module defines the XrDebye class for calculation 

7of X-ray scattering properties from atomic cluster 

8using Debye formula. 

9Also contains routine for calculation of atomic form factors and 

10X-ray wavelength dict. 

11""" 

12 

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

14 

15import numpy as np 

16 

17from ase.data import atomic_numbers 

18 

19# Table (1) of 

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

21waasmaier = { 

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

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

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

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

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

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

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

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

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

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

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

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

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

35} 

36 

37wavelengths = { 

38 'CuKa1': 1.5405981, 

39 'CuKa2': 1.54443, 

40 'CuKb1': 1.39225, 

41 'WLa1': 1.47642, 

42 'WLa2': 1.48748 

43} 

44 

45 

46class XrDebye: 

47 """ 

48 Class for calculation of XRD or SAXS patterns. 

49 """ 

50 

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

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

53 """ 

54 Initilize the calculation of X-ray diffraction patterns 

55 

56 Parameters: 

57 

58 atoms: ase.Atoms 

59 atoms object for which calculation will be performed. 

60 

61 wavelength: float, Angstrom 

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

63 

64 damping : float, Angstrom**2 

65 thermal damping factor parameter (B-factor). 

66 

67 method: {'Iwasa'} 

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

69 

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

71 atomic factors are used. 

72 

73 For any other string there will be only thermal damping 

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

75 

76 alpha: float 

77 parameter for angular damping of scattering intensity. 

78 Close to 1.0 for unplorized beam. 

79 

80 warn: boolean 

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

82 """ 

83 self.wavelength = wavelength 

84 self.damping = damping 

85 self.mode = '' 

86 self.method = method 

87 self.alpha = alpha 

88 self.warn = warn 

89 

90 self.twotheta_list = [] 

91 self.q_list = [] 

92 self.intensity_list = [] 

93 

94 self.atoms = atoms 

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

96 

97 def set_damping(self, damping): 

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

99 self.damping = damping 

100 

101 def get(self, s): 

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

103 using the Debye-Formula at single point. 

104 

105 Parameters: 

106 

107 s: float, in inverse Angstrom 

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

109 

110 Returns: 

111 Intensity at given scattering vector `s`. 

112 """ 

113 

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

115 

116 if self.method == 'Iwasa': 

117 sinth = self.wavelength * s / 2. 

118 positive = 1. - sinth**2 

119 if positive < 0: 

120 positive = 0 

121 costh = sqrt(positive) 

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

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

124 

125 f = {} 

126 

127 def atomic(symbol): 

128 """ 

129 get atomic factor, using cache. 

130 """ 

131 if symbol not in f: 

132 if self.method == 'Iwasa': 

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

134 else: 

135 f[symbol] = atomic_numbers[symbol] 

136 return f[symbol] 

137 

138 I = 0. 

139 fa = [] # atomic factors list 

140 for a in self.atoms: 

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

142 

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

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

145 

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

147 vr = pos - pos[i] 

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

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

150 

151 return pre * I 

152 

153 def get_waasmaier(self, symbol, s): 

154 r"""Scattering factor for free atoms. 

155 

156 Parameters: 

157 

158 symbol: string 

159 atom element symbol. 

160 

161 s: float, in inverse Angstrom 

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

163 

164 Returns: 

165 Intensity at given scattering vector `s`. 

166 

167 Note: 

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

169 if symbol == 'H': 

170 # XXXX implement analytical H 

171 return 0 

172 elif symbol in waasmaier: 

173 abc = waasmaier[symbol] 

174 f = abc[10] 

175 s2 = s * s 

176 for i in range(5): 

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

178 return f 

179 if self.warn: 

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

181 return 0 

182 

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

184 r""" 

185 Calculate X-ray diffraction pattern or 

186 small angle X-ray scattering pattern. 

187 

188 Parameters: 

189 

190 x: float array 

191 points where intensity will be calculated. 

192 XRD - 2theta values, in degrees; 

193 SAXS - q values in 1/A 

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

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

196 

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

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

199 small-angle scattering (SAXS). 

200 

201 Returns: 

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

203 """ 

204 self.mode = mode.upper() 

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

206 

207 result = [] 

208 if mode == 'XRD': 

209 if x is None: 

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

211 else: 

212 self.twotheta_list = x 

213 self.q_list = [] 

214 if verbose: 

215 print('#2theta\tIntensity') 

216 for twotheta in self.twotheta_list: 

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

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

219 if verbose: 

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

221 elif mode == 'SAXS': 

222 if x is None: 

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

224 else: 

225 self.q_list = x 

226 self.twotheta_list = [] 

227 if verbose: 

228 print('#q\tIntensity') 

229 for q in self.q_list: 

230 s = q / (2 * pi) 

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

232 if verbose: 

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

234 self.intensity_list = np.array(result) 

235 return self.intensity_list 

236 

237 def write_pattern(self, filename): 

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

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

240 self._write_pattern(fd) 

241 

242 def _write_pattern(self, fd): 

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

244 if self.mode == 'XRD': 

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

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

247 elif self.mode == 'SAXS': 

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

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

250 else: 

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

252 

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

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

255 

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

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

258 

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

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

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

262 

263 Returns: 

264 ``matplotlib.axes.Axes`` object.""" 

265 

266 import matplotlib.pyplot as plt 

267 

268 if ax is None: 

269 plt.clf() # clear figure 

270 ax = plt.gca() 

271 

272 if self.mode == 'XRD': 

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

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

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

276 ax.set_ylabel('Intensity') 

277 elif self.mode == 'SAXS': 

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

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

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

281 ax.set_ylabel('Intensity') 

282 else: 

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

284 

285 if show: 

286 plt.show() 

287 if filename is not None: 

288 fig = ax.get_figure() 

289 fig.savefig(filename) 

290 

291 return ax