1# coding: utf-8 2# Copyright (c) Pymatgen Development Team. 3# Distributed under the terms of the MIT License. 4 5""" 6This module implements an XRD pattern calculator. 7""" 8 9import json 10import os 11from math import asin, cos, degrees, pi, radians, sin 12 13import numpy as np 14 15from pymatgen.symmetry.analyzer import SpacegroupAnalyzer 16 17from .core import ( 18 AbstractDiffractionPatternCalculator, 19 DiffractionPattern, 20 get_unique_families, 21) 22 23# XRD wavelengths in angstroms 24WAVELENGTHS = { 25 "CuKa": 1.54184, 26 "CuKa2": 1.54439, 27 "CuKa1": 1.54056, 28 "CuKb1": 1.39222, 29 "MoKa": 0.71073, 30 "MoKa2": 0.71359, 31 "MoKa1": 0.70930, 32 "MoKb1": 0.63229, 33 "CrKa": 2.29100, 34 "CrKa2": 2.29361, 35 "CrKa1": 2.28970, 36 "CrKb1": 2.08487, 37 "FeKa": 1.93735, 38 "FeKa2": 1.93998, 39 "FeKa1": 1.93604, 40 "FeKb1": 1.75661, 41 "CoKa": 1.79026, 42 "CoKa2": 1.79285, 43 "CoKa1": 1.78896, 44 "CoKb1": 1.63079, 45 "AgKa": 0.560885, 46 "AgKa2": 0.563813, 47 "AgKa1": 0.559421, 48 "AgKb1": 0.497082, 49} 50 51with open(os.path.join(os.path.dirname(__file__), "atomic_scattering_params.json")) as f: 52 ATOMIC_SCATTERING_PARAMS = json.load(f) 53 54 55class XRDCalculator(AbstractDiffractionPatternCalculator): 56 r""" 57 Computes the XRD pattern of a crystal structure. 58 59 This code is implemented by Shyue Ping Ong as part of UCSD's NANO106 - 60 Crystallography of Materials. The formalism for this code is based on 61 that given in Chapters 11 and 12 of Structure of Materials by Marc De 62 Graef and Michael E. McHenry. This takes into account the atomic 63 scattering factors and the Lorentz polarization factor, but not 64 the Debye-Waller (temperature) factor (for which data is typically not 65 available). Note that the multiplicity correction is not needed since 66 this code simply goes through all reciprocal points within the limiting 67 sphere, which includes all symmetrically equivalent facets. The algorithm 68 is as follows 69 70 1. Calculate reciprocal lattice of structure. Find all reciprocal points 71 within the limiting sphere given by :math:`\\frac{2}{\\lambda}`. 72 73 2. For each reciprocal point :math:`\\mathbf{g_{hkl}}` corresponding to 74 lattice plane :math:`(hkl)`, compute the Bragg condition 75 :math:`\\sin(\\theta) = \\frac{\\lambda}{2d_{hkl}}` 76 77 3. Compute the structure factor as the sum of the atomic scattering 78 factors. The atomic scattering factors are given by 79 80 .. math:: 81 82 f(s) = Z - 41.78214 \\times s^2 \\times \\sum\\limits_{i=1}^n a_i \ 83 \\exp(-b_is^2) 84 85 where :math:`s = \\frac{\\sin(\\theta)}{\\lambda}` and :math:`a_i` 86 and :math:`b_i` are the fitted parameters for each element. The 87 structure factor is then given by 88 89 .. math:: 90 91 F_{hkl} = \\sum\\limits_{j=1}^N f_j \\exp(2\\pi i \\mathbf{g_{hkl}} 92 \\cdot \\mathbf{r}) 93 94 4. The intensity is then given by the modulus square of the structure 95 factor. 96 97 .. math:: 98 99 I_{hkl} = F_{hkl}F_{hkl}^* 100 101 5. Finally, the Lorentz polarization correction factor is applied. This 102 factor is given by: 103 104 .. math:: 105 106 P(\\theta) = \\frac{1 + \\cos^2(2\\theta)} 107 {\\sin^2(\\theta)\\cos(\\theta)} 108 """ 109 110 # Tuple of available radiation keywords. 111 AVAILABLE_RADIATION = tuple(WAVELENGTHS.keys()) 112 113 def __init__(self, wavelength="CuKa", symprec=0, debye_waller_factors=None): 114 """ 115 Initializes the XRD calculator with a given radiation. 116 117 Args: 118 wavelength (str/float): The wavelength can be specified as either a 119 float or a string. If it is a string, it must be one of the 120 supported definitions in the AVAILABLE_RADIATION class 121 variable, which provides useful commonly used wavelengths. 122 If it is a float, it is interpreted as a wavelength in 123 angstroms. Defaults to "CuKa", i.e, Cu K_alpha radiation. 124 symprec (float): Symmetry precision for structure refinement. If 125 set to 0, no refinement is done. Otherwise, refinement is 126 performed using spglib with provided precision. 127 debye_waller_factors ({element symbol: float}): Allows the 128 specification of Debye-Waller factors. Note that these 129 factors are temperature dependent. 130 """ 131 if isinstance(wavelength, (float, int)): 132 self.wavelength = wavelength 133 elif isinstance(wavelength, str): 134 self.radiation = wavelength 135 self.wavelength = WAVELENGTHS[wavelength] 136 else: 137 raise TypeError("'wavelength' must be either of: float, int or str") 138 self.symprec = symprec 139 self.debye_waller_factors = debye_waller_factors or {} 140 141 def get_pattern(self, structure, scaled=True, two_theta_range=(0, 90)): 142 """ 143 Calculates the diffraction pattern for a structure. 144 145 Args: 146 structure (Structure): Input structure 147 scaled (bool): Whether to return scaled intensities. The maximum 148 peak is set to a value of 100. Defaults to True. Use False if 149 you need the absolute values to combine XRD plots. 150 two_theta_range ([float of length 2]): Tuple for range of 151 two_thetas to calculate in degrees. Defaults to (0, 90). Set to 152 None if you want all diffracted beams within the limiting 153 sphere of radius 2 / wavelength. 154 155 Returns: 156 (XRDPattern) 157 """ 158 if self.symprec: 159 finder = SpacegroupAnalyzer(structure, symprec=self.symprec) 160 structure = finder.get_refined_structure() 161 162 wavelength = self.wavelength 163 latt = structure.lattice 164 is_hex = latt.is_hexagonal() 165 166 # Obtained from Bragg condition. Note that reciprocal lattice 167 # vector length is 1 / d_hkl. 168 min_r, max_r = ( 169 (0, 2 / wavelength) 170 if two_theta_range is None 171 else [2 * sin(radians(t / 2)) / wavelength for t in two_theta_range] 172 ) 173 174 # Obtain crystallographic reciprocal lattice points within range 175 recip_latt = latt.reciprocal_lattice_crystallographic 176 recip_pts = recip_latt.get_points_in_sphere([[0, 0, 0]], [0, 0, 0], max_r) 177 if min_r: 178 recip_pts = [pt for pt in recip_pts if pt[1] >= min_r] 179 180 # Create a flattened array of zs, coeffs, fcoords and occus. This is 181 # used to perform vectorized computation of atomic scattering factors 182 # later. Note that these are not necessarily the same size as the 183 # structure as each partially occupied specie occupies its own 184 # position in the flattened array. 185 zs = [] 186 coeffs = [] 187 fcoords = [] 188 occus = [] 189 dwfactors = [] 190 191 for site in structure: 192 for sp, occu in site.species.items(): 193 zs.append(sp.Z) 194 try: 195 c = ATOMIC_SCATTERING_PARAMS[sp.symbol] 196 except KeyError: 197 raise ValueError( 198 "Unable to calculate XRD pattern as " 199 "there is no scattering coefficients for" 200 " %s." % sp.symbol 201 ) 202 coeffs.append(c) 203 dwfactors.append(self.debye_waller_factors.get(sp.symbol, 0)) 204 fcoords.append(site.frac_coords) 205 occus.append(occu) 206 207 zs = np.array(zs) 208 coeffs = np.array(coeffs) 209 fcoords = np.array(fcoords) 210 occus = np.array(occus) 211 dwfactors = np.array(dwfactors) 212 peaks = {} 213 two_thetas = [] 214 215 for hkl, g_hkl, ind, _ in sorted(recip_pts, key=lambda i: (i[1], -i[0][0], -i[0][1], -i[0][2])): 216 # Force miller indices to be integers. 217 hkl = [int(round(i)) for i in hkl] 218 if g_hkl != 0: 219 220 d_hkl = 1 / g_hkl 221 222 # Bragg condition 223 theta = asin(wavelength * g_hkl / 2) 224 225 # s = sin(theta) / wavelength = 1 / 2d = |ghkl| / 2 (d = 226 # 1/|ghkl|) 227 s = g_hkl / 2 228 229 # Store s^2 since we are using it a few times. 230 s2 = s ** 2 231 232 # Vectorized computation of g.r for all fractional coords and 233 # hkl. 234 g_dot_r = np.dot(fcoords, np.transpose([hkl])).T[0] 235 236 # Highly vectorized computation of atomic scattering factors. 237 # Equivalent non-vectorized code is:: 238 # 239 # for site in structure: 240 # el = site.specie 241 # coeff = ATOMIC_SCATTERING_PARAMS[el.symbol] 242 # fs = el.Z - 41.78214 * s2 * sum( 243 # [d[0] * exp(-d[1] * s2) for d in coeff]) 244 fs = zs - 41.78214 * s2 * np.sum(coeffs[:, :, 0] * np.exp(-coeffs[:, :, 1] * s2), axis=1) 245 246 dw_correction = np.exp(-dwfactors * s2) 247 248 # Structure factor = sum of atomic scattering factors (with 249 # position factor exp(2j * pi * g.r and occupancies). 250 # Vectorized computation. 251 f_hkl = np.sum(fs * occus * np.exp(2j * pi * g_dot_r) * dw_correction) 252 253 # Lorentz polarization correction for hkl 254 lorentz_factor = (1 + cos(2 * theta) ** 2) / (sin(theta) ** 2 * cos(theta)) 255 256 # Intensity for hkl is modulus square of structure factor. 257 i_hkl = (f_hkl * f_hkl.conjugate()).real 258 259 two_theta = degrees(2 * theta) 260 261 if is_hex: 262 # Use Miller-Bravais indices for hexagonal lattices. 263 hkl = (hkl[0], hkl[1], -hkl[0] - hkl[1], hkl[2]) 264 # Deal with floating point precision issues. 265 ind = np.where( 266 np.abs(np.subtract(two_thetas, two_theta)) < AbstractDiffractionPatternCalculator.TWO_THETA_TOL 267 ) 268 if len(ind[0]) > 0: 269 peaks[two_thetas[ind[0][0]]][0] += i_hkl * lorentz_factor 270 peaks[two_thetas[ind[0][0]]][1].append(tuple(hkl)) 271 else: 272 peaks[two_theta] = [i_hkl * lorentz_factor, [tuple(hkl)], d_hkl] 273 two_thetas.append(two_theta) 274 275 # Scale intensities so that the max intensity is 100. 276 max_intensity = max([v[0] for v in peaks.values()]) 277 x = [] 278 y = [] 279 hkls = [] 280 d_hkls = [] 281 for k in sorted(peaks.keys()): 282 v = peaks[k] 283 fam = get_unique_families(v[1]) 284 if v[0] / max_intensity * 100 > AbstractDiffractionPatternCalculator.SCALED_INTENSITY_TOL: 285 x.append(k) 286 y.append(v[0]) 287 hkls.append([{"hkl": hkl, "multiplicity": mult} for hkl, mult in fam.items()]) 288 d_hkls.append(v[2]) 289 xrd = DiffractionPattern(x, y, hkls, d_hkls) 290 if scaled: 291 xrd.normalize(mode="max", value=100) 292 return xrd 293