Source code for matador.fingerprints.pxrd

# coding: utf-8
# Distributed under the terms of the MIT License.

""" This file implements the PXRD class for simulating powder XRD pattern
of a crystal.

"""

import itertools
import os
from typing import Tuple

import numpy as np

from matador.fingerprints.fingerprint import Fingerprint, FingerprintFactory
from matador.crystal import Crystal
from matador.utils.cell_utils import standardize_doc_cell, real2recip
from matador.utils.chem_utils import get_formula_from_stoich

THETA_TOL = 1e-5


[docs]class PXRD(Fingerprint): """ This class for computes powder X-ray diffraction patterns of a given crystal for a certain incident wavelength. The cell is standardised with spglib before computing PXRD. This calculation takes into account atomic scattering factors, Lorentz polarisation and thermal broadening (with Debye-Waller factors set to 1). Note: this class does not perform any q-dependent peak broadening, and instead uses a simple Lorentzian broadening. The default width of 0.03 provides good agreement with e.g. GSAS-II's default CuKa setup. Only one wavelength can be used at a time, but multiple patterns could be combined post hoc. Attributes: self.peak_positions (numpy.ndarray): sorted peak positions as values in 2θ self.hkls (numpy.ndarray): Miller indices correspnding to peaks, sorted by peak angle. self.peak_intensities (numpy.ndarray): intensity of each peak. self.pattern (numpy.ndarray): Lorentzian-broadened pattern at values of self.two_thetas. self.two_thetas (numpy.ndarray): continuous space of 2θ values corresponding to sample points of self.pattern. """ def __init__( self, doc, wavelength: float = 1.5406, lorentzian_width: float = 0.03, two_theta_resolution: float = 0.01, two_theta_bounds: Tuple[float, float] = (0, 90), theta_m: float = 0.0, scattering_factors: str = "RASPA", lazy=False, plot=False, progress=False, *args, **kwargs ): """ Set up the PXRD, and compute it, if lazy is False. Parameters: doc (dict/Crystal): matador document to compute PXRD for. Keyword arguments: lorentzian_width (float): width of Lorentzians for broadening (DEFAULT: 0.03) wavelength (float): incident X-ray wavelength (DEFAULT: CuKa, 1.5406). theta_m (float): the monochromator angle in degrees (DEFAULT: 0) two_theta_resolution (float): resolution of grid 2θ used for plotting. two_theta_bounds (tuple of float): values between which to compute the PXRD pattern. scattering_factors (str): either "GSAS" or "RASPA" (default), which set of atomic scattering factors to use. lazy (bool): whether to compute PXRD or just set it up. plot (bool): whether to display PXRD as a plot. """ self.wavelength = wavelength self.lorentzian_width = lorentzian_width self.two_theta_resolution = two_theta_resolution if two_theta_bounds is not None: self.two_theta_bounds = list(two_theta_bounds) else: self.two_theta_bounds = [0, 90] self.theta_m = theta_m self.scattering_factors = scattering_factors self.progress = progress if self.two_theta_bounds[0] < THETA_TOL: self.two_theta_bounds[0] = THETA_TOL if np.min(doc.get('site_occupancy', [1.0])) < 1.0: print("System has partial occupancy, not refining with spglib.") self.doc = Crystal(doc) else: self.doc = Crystal(standardize_doc_cell(doc, primitive=True)) self.formula = get_formula_from_stoich(self.doc['stoichiometry'], tex=True) self.spg = self.doc['space_group'] species = list(set(self.doc['atom_types'])) # this could be cached across PXRD objects but is much faster than the XRD calculation itself if self.scattering_factors == "GSAS": from matador.data.atomic_scattering import GSAS_ATOMIC_SCATTERING_COEFFS self.atomic_scattering_coeffs = {spec: GSAS_ATOMIC_SCATTERING_COEFFS[spec] for spec in species} elif self.scattering_factors == "RASPA": from matador.data.atomic_scattering import RASPA_ATOMIC_SCATTERING_COEFFS self.atomic_scattering_coeffs = {spec: RASPA_ATOMIC_SCATTERING_COEFFS[spec] for spec in species} else: raise RuntimeError( "No set of scattering factors matched: {}. Please use 'GSAS' or 'RASPA'." .format(self.scattering_factors) ) if not lazy: self.calculate() if plot: self.plot()
[docs] def calc_pxrd(self): """ Calculate the PXRD pattern. """ # set crystallographic data lattice_abc = np.asarray(self.doc.lattice_abc) lattice_cart = np.asarray(self.doc.lattice_cart) positions_abs = np.asarray(self.doc.positions_abs) site_occupancies = np.asarray(self.doc.site_occupancies) # find allowed reciprocal lattice points within limiting sphere min_r, max_r = [2 / self.wavelength * np.sin(np.pi / 180 * t / 2) for t in self.two_theta_bounds] Ns = np.floor(max_r * lattice_abc[0, :]).astype(int) recip = np.asarray(real2recip(lattice_cart)).T qs = np.zeros((2*sum(Ns), 3), dtype=np.float64) hkls = np.asarray( list(itertools.product( range(-Ns[0], Ns[0] + 1), range(-Ns[1], Ns[1] + 1), range(-Ns[2], Ns[2] + 1) )), dtype=np.float64 ) hkls = hkls[np.argsort(np.linalg.norm(hkls, axis=-1))] qs = np.dot(recip, hkls.T).T # filter out by theta bounds q_mags = np.linalg.norm(qs, axis=-1) allowed = np.where(np.logical_and(q_mags <= max_r * 2 * np.pi, q_mags >= min_r * 2 * np.pi)) qs = qs[allowed] hkls = hkls[allowed] # compute Bragg condition to find peak locations sin_tau = np.linalg.norm(qs, axis=1) * self.wavelength / (4 * np.pi) sin_tau[sin_tau > 1] = 0 taus = 2 * np.arcsin(sin_tau) # compute structure factor S(q) as sum of atomic scattering factors S_q = np.zeros_like(taus) if self.progress: import tqdm bar = tqdm.tqdm else: def bar(x): return x for ind, q_vector in bar(enumerate(qs)): # accumulate atomic scattering factors atomic_factor = {} for species in set(self.doc.atom_types): atomic_factor[species] = self.atomic_scattering_factor(q_mags[ind], species) factors = np.array([atomic_factor[species] for species in self.doc.atom_types]) F_s = np.sum(np.exp(1j * positions_abs @ q_vector) * factors * site_occupancies) S_q[ind] = np.abs(F_s)**2 # apply Lorentz correction for polarisation and finite size effects S_q *= 2 * (1 + np.cos(taus) ** 2 * np.cos(2 * self.theta_m) ** 2) / (np.sin(taus) * np.sin(0.5 * taus)) # thermal correction assuming no Debye-Waller factor S_q *= np.exp(-np.sin(taus)**2 / self.wavelength**2)**2 S_q /= np.max(S_q) # create histogram and broaden onto 2 theta space in degrees self.peak_positions = (180 / np.pi) * taus self.hkls = hkls self.two_thetas = np.arange(self.two_theta_bounds[0], self.two_theta_bounds[1] + self.two_theta_resolution, self.two_theta_resolution) self.pattern, bins = np.histogram(self.peak_positions, bins=self.two_thetas, weights=S_q) if self.lorentzian_width > 0: self.pattern = self._broadening_unrolled( self.pattern, self.two_thetas, self.lorentzian_width, broadening_type='lorentzian' ) else: # shift and clip the last two theta value if we didnt do broadening self.two_thetas = self.two_thetas[:-1] + self.two_theta_resolution / 2 self.pattern /= np.max(self.pattern) order = np.argsort(self.peak_positions) self.hkls = hkls[order] self.peak_intensities = S_q self.peak_positions = self.peak_positions[order] # alias old name for compatibility... self.spectrum = self.pattern
[docs] def calculate(self): """ Alias for calculating the PXRD pattern. """ self.calc_pxrd()
[docs] def atomic_scattering_factor(self, q_mag, species): """ Return fit for particular atom at given q-vector. Parameters: q_mag (float): magnitude of the q_vector. species (str): the element label. Returns: float: the atomic scattering factor. """ a = self.atomic_scattering_coeffs[species][0] b = self.atomic_scattering_coeffs[species][1] c = self.atomic_scattering_coeffs[species][2] return c + np.sum(a * np.exp(-b * (q_mag / (4 * np.pi))**2))
[docs] def plot(self, **kwargs): """ Wrapper function to plot the PXRD pattern. """ from matador.plotting.pxrd_plotting import plot_pxrd plot_pxrd(self, **kwargs)
[docs] def save_pattern(pxrd, fname): """ Write a file to `fname` that contains the xy coordinates of the PXRD pattern. """ if os.path.isfile(fname): raise RuntimeError(f"Requested filename {fname} already exists!") from matador import __version__ header = f""" PXRD pattern computed with matador {__version__}. Input file: {pxrd.doc.source[0]} Structure: {pxrd.doc} Settings: wavelength = {pxrd.wavelength} Å theta_m = {pxrd.theta_m} degrees lorentzian_width = {pxrd.lorentzian_width} degrees 2θ (degrees),\t\t\tRelative intensity""" np.savetxt(fname, np.vstack([pxrd.two_thetas, pxrd.pattern]).T, header=header, fmt='%.14e', delimiter='\t')
[docs] def save_peaks(pxrd, fname): """ Write a file to `fname` that contains the peak list. """ if os.path.isfile(fname): raise RuntimeError(f"Requested filename {fname} already exists!") from matador import __version__ header = f""" PXRD peaks computed with matador {__version__}. Input file: {pxrd.doc.source[0]} Structure: {pxrd.doc} Settings: wavelength = {pxrd.wavelength} Å theta_m = {pxrd.theta_m} degrees lorentzian_width = {pxrd.lorentzian_width} degrees <hkl>,\t\tPeak position (degrees)""" np.savetxt( fname, np.vstack( [pxrd.hkls[:, 0], pxrd.hkls[:, 1], pxrd.hkls[:, 2], pxrd.peak_positions] ).T, header=header, fmt=['% d', '% d', '% d', '%.14e'], delimiter='\t' )
[docs]class PXRDFactory(FingerprintFactory): fingerprint = PXRD default_key = 'pxrd'