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