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 exp, pi, sin, sqrt, cos, acos 

12import numpy as np 

13 

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 def __init__(self, atoms, wavelength, damping=0.04, 

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

50 """ 

51 Initilize the calculation of X-ray diffraction patterns 

52 

53 Parameters: 

54 

55 atoms: ase.Atoms 

56 atoms object for which calculation will be performed. 

57 

58 wavelength: float, Angstrom 

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

60 

61 damping : float, Angstrom**2 

62 thermal damping factor parameter (B-factor). 

63 

64 method: {'Iwasa'} 

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

66 

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

68 atomic factors are used. 

69 

70 For any other string there will be only thermal damping 

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

72 

73 alpha: float 

74 parameter for angular damping of scattering intensity. 

75 Close to 1.0 for unplorized beam. 

76 

77 warn: boolean 

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

79 """ 

80 self.wavelength = wavelength 

81 self.damping = damping 

82 self.mode = '' 

83 self.method = method 

84 self.alpha = alpha 

85 self.warn = warn 

86 

87 self.twotheta_list = [] 

88 self.q_list = [] 

89 self.intensity_list = [] 

90 

91 self.atoms = atoms 

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

93 

94 def set_damping(self, damping): 

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

96 self.damping = damping 

97 

98 def get(self, s): 

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

100 using the Debye-Formula at single point. 

101 

102 Parameters: 

103 

104 s: float, in inverse Angstrom 

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

106 

107 Returns: 

108 Intensity at given scattering vector `s`. 

109 """ 

110 

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

112 

113 if self.method == 'Iwasa': 

114 sinth = self.wavelength * s / 2. 

115 positive = 1. - sinth**2 

116 if positive < 0: 

117 positive = 0 

118 costh = sqrt(positive) 

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

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

121 

122 f = {} 

123 def atomic(symbol): 

124 """ 

125 get atomic factor, using cache. 

126 """ 

127 if symbol not in f: 

128 if self.method == 'Iwasa': 

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

130 else: 

131 f[symbol] = atomic_numbers[symbol] 

132 return f[symbol] 

133 

134 I = 0. 

135 fa = [] # atomic factors list 

136 for a in self.atoms: 

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

138 

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

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

141 

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

143 vr = pos - pos[i] 

144 I += np.sum(fa[i] * fa * np.sinc(2 * s * np.sqrt(np.sum(vr * vr, axis=1)))) 

145 

146 return pre * I 

147 

148 def get_waasmaier(self, symbol, s): 

149 r"""Scattering factor for free atoms. 

150 

151 Parameters: 

152 

153 symbol: string 

154 atom element symbol. 

155 

156 s: float, in inverse Angstrom 

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

158 

159 Returns: 

160 Intensity at given scattering vector `s`. 

161 

162 Note: 

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

164 if symbol == 'H': 

165 # XXXX implement analytical H 

166 return 0 

167 elif symbol in waasmaier: 

168 abc = waasmaier[symbol] 

169 f = abc[10] 

170 s2 = s * s 

171 for i in range(5): 

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

173 return f 

174 if self.warn: 

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

176 return 0 

177 

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

179 r""" 

180 Calculate X-ray diffraction pattern or 

181 small angle X-ray scattering pattern. 

182 

183 Parameters: 

184 

185 x: float array 

186 points where intensity will be calculated. 

187 XRD - 2theta values, in degrees; 

188 SAXS - q values in 1/A 

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

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

191 

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

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

194 small-angle scattering (SAXS). 

195 

196 Returns: 

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

198 """ 

199 self.mode = mode.upper() 

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

201 

202 result = [] 

203 if mode == 'XRD': 

204 if x is None: 

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

206 else: 

207 self.twotheta_list = x 

208 self.q_list = [] 

209 if verbose: 

210 print('#2theta\tIntensity') 

211 for twotheta in self.twotheta_list: 

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

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

214 if verbose: 

215 print('%.3f\t%f' % (twotheta, result[-1])) 

216 elif mode == 'SAXS': 

217 if x is None: 

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

219 else: 

220 self.q_list = x 

221 self.twotheta_list = [] 

222 if verbose: 

223 print('#q\tIntensity') 

224 for q in self.q_list: 

225 s = q / (2 * pi) 

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

227 if verbose: 

228 print('%.4f\t%f' % (q, result[-1])) 

229 self.intensity_list = np.array(result) 

230 return self.intensity_list 

231 

232 def write_pattern(self, filename): 

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

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

235 self._write_pattern(fd) 

236 

237 def _write_pattern(self, fd): 

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

239 if self.mode == 'XRD': 

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

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

242 elif self.mode == 'SAXS': 

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

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

245 else: 

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

247 

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

249 fd.write(' %f\t%f\n' % (x[i], y[i])) 

250 

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

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

253 

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

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

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

257 

258 Returns: 

259 ``matplotlib.axes.Axes`` object.""" 

260 

261 import matplotlib.pyplot as plt 

262 

263 if ax is None: 

264 plt.clf() # clear figure 

265 ax = plt.gca() 

266 

267 if self.mode == 'XRD': 

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

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

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

271 ax.set_ylabel('Intensity') 

272 elif self.mode == 'SAXS': 

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

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

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

276 ax.set_ylabel('Intensity') 

277 else: 

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

279 

280 if show: 

281 plt.show() 

282 if filename is not None: 

283 fig = ax.get_figure() 

284 fig.savefig(filename) 

285 

286 return ax