diff --git a/doc/changelog.rst b/doc/changelog.rst index d79cd306..1c818f05 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -34,6 +34,9 @@ Added detector, which can be added to an EBSD signal as markers. (`#204 `_, `#219 `_) +- Dependency on the diffsims package (https://github.com/pyxem/diffsims/) for + handling of electron scattering and diffraction. + (`#220 `_) Changed ------- diff --git a/kikuchipy/conftest.py b/kikuchipy/conftest.py index 566846ad..65659f4c 100644 --- a/kikuchipy/conftest.py +++ b/kikuchipy/conftest.py @@ -21,13 +21,13 @@ import tempfile from diffpy.structure import Atom, Lattice, Structure +from diffsims.crystallography import ReciprocalLatticePoint import numpy as np from orix.crystal_map import Phase from orix.quaternion.rotation import Rotation from orix.vector import Vector3d, neo_euler import pytest -from kikuchipy.crystallography import ReciprocalLatticePoint from kikuchipy.detectors import EBSDDetector from kikuchipy.generators import EBSDSimulationGenerator from kikuchipy.signals import EBSD diff --git a/kikuchipy/crystallography/__init__.py b/kikuchipy/crystallography/__init__.py index ea82c49e..f58d299c 100644 --- a/kikuchipy/crystallography/__init__.py +++ b/kikuchipy/crystallography/__init__.py @@ -21,12 +21,7 @@ from diffpy.structure import Lattice import numpy as np -from kikuchipy.crystallography.reciprocal_lattice_point import ( - ReciprocalLatticePoint, -) - __all__ = [ - "ReciprocalLatticePoint", "get_direct_structure_matrix", "get_reciprocal_structure_matrix", "get_reciprocal_metric_tensor", diff --git a/kikuchipy/crystallography/atomic_scattering_factor.py b/kikuchipy/crystallography/atomic_scattering_factor.py deleted file mode 100644 index cce2d1a4..00000000 --- a/kikuchipy/crystallography/atomic_scattering_factor.py +++ /dev/null @@ -1,103 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright 2017-2020 The diffsims developers -# -# This file is part of diffsims. -# -# diffsims is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# diffsims is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with diffsims. If not, see . - -# TODO: This file will be moved to diffsims - -import numpy as np - -from kikuchipy.crystallography.atomic_scattering_parameters import ( - get_atomic_scattering_parameters, - get_element_id_from_string, -) - - -def get_kinematical_atomic_scattering_factor(atom, scattering_parameter): - """Return the kinematical (X-ray) atomic scattering factor f for a - certain atom and scattering parameter. - - Assumes structure's Debye-Waller factors are expressed in Ångströms. - - Parameters - ---------- - atom : diffpy.structure.Atom - Atom with element type, Debye-Waller factor and occupancy number. - scattering_parameter : float - The scattering parameter s for these Miller indices describing - the crystal plane in which the atom lies. - - Returns - ------- - f : float - Scattering factor for this atom on this plane. - """ - # Get the atomic scattering parameters - element_id = get_element_id_from_string(atom.element) - a, b = get_atomic_scattering_parameters(element_id) - - # Get the scattering parameter squared - s2 = scattering_parameter ** 2 - - # Get the atomic scattering factor - f = element_id - (41.78214 * s2 * np.sum(a * np.exp(-b * s2))) - - # Correct for occupancy and the Debye-Waller factor - dw_factor = np.exp(-atom.Bisoequiv * s2) - f *= atom.occupancy * dw_factor - - return f - - -def get_doyleturner_atomic_scattering_factor( - atom, scattering_parameter, unit_cell_volume -): - """Return the atomic scattering factor f for a certain atom and - scattering parameter using Doyle-Turner atomic scattering parameters - [Doyle1968]_. - - Assumes structure's Debye-Waller factors are expressed in Ångströms. - - Parameters - ---------- - atom : diffpy.structure.Atom - Atom with element type, Debye-Waller factor and occupancy number. - scattering_parameter : float - The scattering parameter s for these Miller indices describing - the crystal plane in which the atom lies. - unit_cell_volume : float - Volume of the unit cell. - - Returns - ------- - f : float - Scattering factor for this atom on this plane. - """ - # Get the atomic scattering parameters - element_id = get_element_id_from_string(atom.element) - a, b = get_atomic_scattering_parameters(element_id) - - # Get the scattering parameter squared - s2 = scattering_parameter ** 2 - - # Get the atomic scattering factor - f = (47.87801 / unit_cell_volume) * np.sum(a * np.exp(-b * s2)) - - # Correct for occupancy and the Debye-Waller factor - dw_factor = np.exp(-atom.Bisoequiv * s2) - f *= atom.occupancy * dw_factor - - return f diff --git a/kikuchipy/crystallography/atomic_scattering_parameters.py b/kikuchipy/crystallography/atomic_scattering_parameters.py deleted file mode 100644 index aef03497..00000000 --- a/kikuchipy/crystallography/atomic_scattering_parameters.py +++ /dev/null @@ -1,210 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright 2017-2020 The diffsims developers -# -# This file is part of diffsims. -# -# diffsims is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# diffsims is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with diffsims. If not, see . - -# TODO: This file will be moved to diffsims - -import numpy as np - - -def get_atomic_scattering_parameters(element, unit=None): - """Return the eight atomic scattering parameters a_1-4, b_1-4 for - elements with atomic numbers Z = 1-98 from Table 12.1 in - [DeGraef2007]_, which are themselves from [Doyle1968]_ and - [Smith1962]_. - - Parameters - ---------- - element : int or str - Element to return scattering parameters for. Either one-two - letter string or integer atomic number. - unit : str, optional - Either "nm" or "Å"/"A". Whether to return parameters in terms - of Å^-2 or nm^-2. If None (default), Å^-2 is used. - - Returns - ------- - a : np.ndarray - The four atomic scattering parameters a_1-4. - b : np.ndarray - The four atomic scattering parameters b_1-4. - - References - ---------- - .. [DeGraef2007] M. De Graef, M. E. McHenry, "Structure of - \Materials," Cambridge University Press (2007). - .. [Doyle1968] P. A. Doyle, P. S. Turner, "Relativistic Hartree-Fock - X-ray and electron scattering factors," *Acta Cryst.* **24** - (1968), doi: https://doi.org/10.1107/S0567739468000756. - .. [Smith1962] G. Smith, R. Burge, "The analytical representation - of atomic scattering amplitudes for electrons," *Acta Cryst.* - **A15** (1962), doi: https://doi.org/10.1107/S0365110X62000481. - """ - # fmt: off - all_parameters = np.array([ - 0.202, 0.244, 0.082, 0.000, 30.868, 8.544, 1.273, 0.000, # H - 0.091, 0.181, 0.110, 0.036, 18.183, 6.212, 1.803, 0.284, # He - 1.611, 1.246, 0.326, 0.099, 107.638, 30.480, 4.533, 0.495, # etc. - 1.250, 1.334, 0.360, 0.106, 60.804, 18.591, 3.653, 0.416, - 0.945, 1.312, 0.419, 0.116, 46.444, 14.178, 3.223, 0.377, - 0.731, 1.195, 0.456, 0.125, 36.995, 11.297, 2.814, 0.346, - 0.572, 1.043, 0.465, 0.131, 28.847, 9.054, 2.421, 0.317, - 0.455, 0.917, 0.472, 0.138, 23.780, 7.622, 2.144, 0.296, - 0.387, 0.811, 0.475, 0.146, 20.239, 6.609, 1.931, 0.279, - 0.303, 0.720, 0.475, 0.153, 17.640, 5.860, 1.762, 0.266, - 2.241, 1.333, 0.907, 0.286, 108.004, 24.505, 3.391, 0.435, - 2.268, 1.803, 0.839, 0.289, 73.670, 20.175, 3.013, 0.405, - 2.276, 2.428, 0.858, 0.317, 72.322, 19.773, 3.080, 0.408, - 2.129, 2.533, 0.835, 0.322, 57.775, 16.476, 2.880, 0.386, - 1.888, 2.469, 0.805, 0.320, 44.876, 13.538, 2.642, 0.361, - 1.659, 2.386, 0.790, 0.321, 36.650, 11.488, 2.469, 0.340, - 1.452, 2.292, 0.787, 0.322, 30.935, 9.980, 2.234, 0.323, - 1.274, 2.190, 0.793, 0.326, 26.682, 8.813, 2.219, 0.307, - 3.951, 2.545, 1.980, 0.482, 137.075, 22.402, 4.532, 0.434, - 4.470, 2.971, 1.970, 0.482, 99.523, 22.696, 4.195, 0.417, - 3.966, 2.917, 1.925, 0.480, 88.960, 20.606, 3.856, 0.399, - 3.565, 2.818, 1.893, 0.483, 81.982, 19.049, 3.590, 0.386, - 3.245, 2.698, 1.860, 0.486, 76.379, 17.726, 3.363, 0.374, - 2.307, 2.334, 1.823, 0.490, 78.405, 15.785, 3.157, 0.364, - 2.747, 2.456, 1.792, 0.498, 67.786, 15.674, 3.000, 0.357, - 2.544, 2.343, 1.759, 0.506, 64.424, 14.880, 2.854, 0.350, - 2.367, 2.236, 1.724, 0.515, 61.431, 14.180, 2.725, 0.344, - 2.210, 2.134, 1.689, 0.524, 58.727, 13.553, 2.609, 0.339, - 1.579, 1.820, 1.658, 0.532, 62.940, 12.453, 2.504, 0.333, - 1.942, 1.950, 1.619, 0.543, 54.162, 12.518, 2.416, 0.330, - 2.321, 2.486, 1.688, 0.599, 65.602, 15.458, 2.581, 0.351, - 2.447, 2.702, 1.616, 0.601, 55.893, 14.393, 2.446, 0.342, - 2.399, 2.790, 1.529, 0.594, 45.718, 12.817, 2.280, 0.328, - 2.298, 2.854, 1.456, 0.590, 38.830, 11.536, 2.146, 0.316, - 2.166, 2.904, 1.395, 0.589, 33.899, 10.497, 2.041, 0.307, - 2.034, 2.927, 1.342, 0.589, 29.999, 9.598, 1.952, 0.299, - 4.776, 3.859, 2.234, 0.868, 140.782, 18.991, 3.701, 0.419, - 5.848, 4.003, 2.342, 0.880, 104.972, 19.367, 3.737, 0.414, - 4.129, 3.012, 1.179, 0.000, 27.548, 5.088, 0.591, 0.0, - 4.105, 3.144, 1.229, 0.000, 28.492, 5.277, 0.601, 0.0, - 4.237, 3.105, 1.234, 0.000, 27.415, 5.074, 0.593, 0.0, - 3.120, 3.906, 2.361, 0.850, 72.464, 14.642, 3.237, 0.366, - 4.318, 3.270, 1.287, 0.000, 28.246, 5.148, 0.590, 0.0, - 4.358, 3.298, 1.323, 0.000, 27.881, 5.179, 0.594, 0.0, - 4.431, 3.343, 1.345, 0.000, 27.911, 5.153, 0.592, 0.0, - 4.436, 3.454, 1.383, 0.000, 28.670, 5.269, 0.595, 0.0, - 2.036, 3.272, 2.511, 0.837, 61.497, 11.824, 2.846, 0.327, - 2.574, 3.259, 2.547, 0.838, 55.675, 11.838, 2.784, 0.322, - 3.153, 3.557, 2.818, 0.884, 66.649, 14.449, 2.976, 0.335, - 3.450, 3.735, 2.118, 0.877, 59.104, 14.179, 2.855, 0.327, - 3.564, 3.844, 2.687, 0.864, 50.487, 13.316, 2.691, 0.316, - 4.785, 3.688, 1.500, 0.000, 27.999, 5.083, 0.581, 0.0, - 3.473, 4.060, 2.522, 0.840, 39.441, 11.816, 2.415, 0.298, - 3.366, 4.147, 2.443, 0.829, 35.509, 11.117, 2.294, 0.289, - 6.062, 5.986, 3.303, 1.096, 155.837, 19.695, 3.335, 0.379, - 7.821, 6.004, 3.280, 1.103, 117.657, 18.778, 3.263, 0.376, - 4.940, 3.968, 1.663, 0.000, 28.716, 5.245, 0.594, 0.0, - 5.007, 3.980, 1.678, 0.000, 28.283, 5.183, 0.589, 0.0, - 5.085, 4.043, 1.684, 0.000, 28.588, 5.143, 0.581, 0.0, - 5.151, 4.075, 1.683, 0.000, 28.304, 5.073, 0.571, 0.0, - 5.201, 4.094, 1.719, 0.000, 28.079, 5.081, 0.576, 0.0, - 5.255, 4.113, 1.743, 0.000, 28.016, 5.037, 0.577, 0.0, - 6.267, 4.844, 3.202, 1.200, 100.298, 16.066, 2.980, 0.367, - 5.225, 4.314, 1.827, 0.000, 29.158, 5.259, 0.586, 0.0, - 5.272, 4.347, 1.844, 0.000, 29.046, 5.226, 0.585, 0.0, - 5.332, 4.370, 1.863, 0.000, 28.888, 5.198, 0.581, 0.0, - 5.376, 4.403, 1.884, 0.000, 28.773, 5.174, 0.582, 0.0, - 5.436, 4.437, 1.891, 0.000, 28.655, 5.117, 0.577, 0.0, - 5.441, 4.510, 1.956, 0.000, 29.149, 5.264, 0.590, 0.0, - 5.529, 4.533, 1.945, 0.000, 28.927, 5.144, 0.578, 0.0, - 5.553, 4.580, 1.969, 0.000, 28.907, 5.160, 0.577, 0.0, - 5.588, 4.619, 1.997, 0.000, 29.001, 5.164, 0.579, 0.0, - 5.659, 4.630, 2.014, 0.000, 28.807, 5.114, 0.578, 0.0, - 5.709, 4.677, 2.019, 0.000, 28.782, 5.084, 0.572, 0.0, - 5.695, 4.740, 2.064, 0.000, 28.968, 5.156, 0.575, 0.0, - 5.750, 4.773, 2.079, 0.000, 28.933, 5.139, 0.573, 0.0, - 5.754, 4.851, 2.096, 0.000, 29.159, 5.152, 0.570, 0.0, - 5.803, 4.870, 2.127, 0.000, 29.016, 5.150, 0.572, 0.0, - 2.388, 4.226, 2.689, 1.255, 42.866, 9.743, 2.264, 0.307, - 2.682, 4.241, 2.755, 1.270, 42.822, 9.856, 2.295, 0.307, - 5.932, 4.972, 2.195, 0.000, 29.086, 5.126, 0.572, 0.0, - 3.510, 4.552, 3.154, 1.359, 52.914, 11.884, 2.571, 0.321, - 3.841, 4.679, 3.192, 1.363, 50.261, 11.999, 2.560, 0.318, - 6.070, 4.997, 2.232, 0.000, 28.075, 4.999, 0.563, 0.0, - 6.133, 5.031, 2.239, 0.000, 28.047, 4.957, 0.558, 0.0, - 4.078, 4.978, 3.096, 1.326, 38.406, 11.020, 2.355, 0.299, - 6.201, 5.121, 2.275, 0.000, 28.200, 4.954, 0.556, 0.0, - 6.215, 5.170, 2.316, 0.000, 28.382, 5.002, 0.562, 0.0, - 6.278, 5.195, 2.321, 0.000, 28.323, 4.949, 0.557, 0.0, - 6.264, 5.263, 2.367, 0.000, 28.651, 5.030, 0.563, 0.0, - 6.306, 5.303, 2.386, 0.000, 28.688, 5.026, 0.561, 0.0, - 6.767, 6.729, 4.014, 1.561, 85.951, 15.642, 2.936, 0.335, - 6.323, 5.414, 2.453, 0.000, 29.142, 5.096, 0.568, 0.0, - 6.415, 5.419, 2.449, 0.000, 28.836, 5.022, 0.561, 0.0, - 6.378, 5.495, 2.495, 0.000, 29.156, 5.102, 0.565, 0.0, - 6.460, 5.469, 2.471, 0.000, 28.396, 4.970, 0.554, 0.0, - 6.502, 5.478, 2.510, 0.000, 28.375, 4.975, 0.561, 0.0, - 6.548, 5.526, 2.520, 0.000, 28.461, 4.965, 0.557, 0.0, - ]).reshape(98, 8) # 1 / Å^2 - # fmt: on - - if isinstance(element, str): - element_id = get_element_id_from_string(element) - else: - element_id = int(element) - - factor = 1 # Ångstrøm - if unit is not None: - if unit.lower() == "nm": - factor = 1e-2 - - a = all_parameters[element_id - 1, :4] * factor - b = all_parameters[element_id - 1, 4:] * factor - - return a, b - - -def get_element_id_from_string(element_str): - """Get periodic element ID for elements Z = 1-98 from one-two letter - string. - - Parameters - ---------- - element_str : str - One-two letter string. - - Returns - ------- - element_id : int - Integer ID in the periodic table of elements. - """ - # List of elements Z = 1-98 - # fmt: off - elements = [ - "h", "he", "li", "be", "b", "c", "n", "o", "f", "ne", "na", "mg", "al", - "si", "p", "s", "cl", "ar", "k", "ca", "sc", "ti", "v", "cr", "mn", - "fe", "co", "ni", "cu", "zn", "ga", "ge", "as", "se", "br", "kr", "rb", - "sr", "y", "zr", "nb", "mo", "tc", "ru", "rh", "pd", "ag", "cd", "in", - "sn", "sb", "te", "i", "xe", "cs", "ba", "la", "ce", "pr", "nd", "pm", - "sm", "eu", "gd", "tb", "dy", "ho", "er", "tm", "yb", "lu", "hf", "ta", - "w", "re", "os", "ir", "pt", "au", "hg", "tl", "pb", "bi", "po", "at", - "rn", "fr", "ra", "ac", "th", "pa", "u", "np", "pu", "am", "cm", "bk", - "cf", "es", "fm", "md", "no", "lr", "rf", "db", "sg", "bh", "hs", "mt", - "ds", "rg", "cn", "nh", "fl", "mc", "lv", "ts", "og" - ] - # fmt: on - n_elements = len(elements) - element2periodic = dict( - zip(elements[:n_elements], np.arange(1, n_elements)) - ) - element_id = element2periodic[element_str.lower()] - return element_id diff --git a/kikuchipy/crystallography/reciprocal_lattice_point.py b/kikuchipy/crystallography/reciprocal_lattice_point.py deleted file mode 100644 index 80f8f112..00000000 --- a/kikuchipy/crystallography/reciprocal_lattice_point.py +++ /dev/null @@ -1,452 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright 2017-2020 The diffsims developers -# -# This file is part of diffsims. -# -# diffsims is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# diffsims is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with diffsims. If not, see . - -# TODO: This file will be moved to diffsims - -from collections import defaultdict -from itertools import product - -import numpy as np -from orix.vector import Vector3d - -from kikuchipy.crystallography.structure_factor import ( - get_kinematical_structure_factor, - get_doyleturner_structure_factor, - get_refraction_corrected_wavelength, -) - - -_FLOAT_EPS = np.finfo(np.float).eps # Used to round values below 1e-16 to zero - - -class ReciprocalLatticePoint: - """Reciprocal lattice point (or crystal plane, reflector, g, etc.) - with Miller indices, length of the reciprocal lattice vectors and - other relevant diffraction parameters. - """ - - def __init__(self, phase, hkl): - """A container for Miller indices, structure factors and related - parameters for crystal planes (reciprocal lattice points, - reflectors, g, etc.). - - Parameters - ---------- - phase : orix.crystal_map.phase_list.Phase - A phase container with a crystal structure and a space and - point group describing the allowed symmetry operations. - hkl : orix.vector.Vector3d, np.ndarray, list, or tuple - Miller indices. - """ - self._hkl = Vector3d(hkl) - self.phase = phase - self._structure_factor = [None] * self.size - self._theta = [None] * self.size - - def __repr__(self): - return ( - f"{self.__class__.__name__} {self.hkl.shape}\n" - f"Phase: {self.phase.name} ({self.phase.point_group.name})\n" - f"{np.array_str(self.hkl.data, precision=4, suppress_small=True)}" - ) - - def __getitem__(self, key): - new_cp = self.__class__(self.phase, self.hkl[key]) - if self.structure_factor[0] is None: - new_cp._structure_factor = [None] * new_cp.size - else: - new_cp._structure_factor = self.structure_factor[key] - return new_cp - - @property - def hkl(self): - """Return :class:`~orix.vector.Vector3d` of Miller indices.""" - return Vector3d(self._hkl.data.astype(int)) - - @property - def _hkldata(self): - """Return :class:`np.ndarray` without 1-dimensions.""" - return np.squeeze(self.hkl.data) - - @property - def h(self): - """Return :class:`np.ndarray` of Miller index h.""" - return self._hkldata[..., 0] - - @property - def k(self): - """Return :class:`np.ndarray` of Miller index k.""" - return self._hkldata[..., 1] - - @property - def l(self): - """Return :class:`np.ndarray` of Miller index l.""" - return self._hkldata[..., 2] - - @property - def size(self): - """Return `int`.""" - return self.hkl.size - - @property - def shape(self): - """Return `tuple`.""" - return self._hkldata.shape - - @property - def multiplicity(self): - """Return either `int` or :class:`np.ndarray` of `int`.""" - return self.symmetrise(antipodal=True, return_multiplicity=True)[1] - - @property - def gspacing(self): - """Return :class:`np.ndarray` of reciprocal lattice point - spacings. - """ - return self.phase.structure.lattice.rnorm(self._hkldata) - - @property - def dspacing(self): - """Return :class:`np.ndarray` of direct lattice interplanar - spacings. - """ - return 1 / self.gspacing - - @property - def scattering_parameter(self): - """Return :class:`np.ndarray` of scattering parameters s.""" - return 0.5 * self.gspacing - - @property - def structure_factor(self): - """Return :class:`np.ndarray` of structure factors F or None.""" - return self._structure_factor - - @property - def allowed(self): - """Return whether planes diffract according to diffraction - selection rules assuming kinematical scattering theory. - """ - # Translational symmetry - centering = self.phase.space_group.short_name[0] - - if centering == "P": # Primitive - if self.phase.space_group.crystal_system == "HEXAGONAL": - # TODO: See rules in e.g. - # https://mcl1.ncifcrf.gov/dauter_pubs/284.pdf, Table 4 - # http://xrayweb.chem.ou.edu/notes/symmetry.html, Systematic Absences - raise NotImplementedError - else: # Any hkl - return np.ones(self.size, dtype=bool) - elif centering == "F": # Face-centred, hkl all odd/even - selection = np.sum(np.mod(self._hkldata, 2), axis=1) - return np.array([i not in [1, 2] for i in selection], dtype=bool) - elif centering == "I": # Body-centred, h + k + l = 2n (even) - return np.mod(np.sum(self._hkldata, axis=1), 2) == 0 - elif centering == "A": # Centred on A faces only - return np.mod(self.k + self.l, 2) == 0 - elif centering == "B": # Centred on B faces only - return np.mod(self.h + self.l, 2) == 0 - elif centering == "C": # Centred on C faces only - return np.mod(self.h + self.k, 2) == 0 - elif centering == "R": # Rhombohedral - return np.mod(-self.h + self.k + self.l, 3) == 0 - - @property - def theta(self): - """Return :class:`np.ndarray` of twice the Bragg angle.""" - return self._theta - - @classmethod - def from_min_dspacing(cls, phase, min_dspacing=0.5): - """Create a CrystalPlane object populated by unique Miller indices - with a direct space interplanar spacing greater than a lower - threshold. - - Parameters - ---------- - phase : orix.crystal_map.phase_list.Phase - A phase container with a crystal structure and a space and - point group describing the allowed symmetry operations. - min_dspacing : float, optional - Smallest interplanar spacing to consider. Default is 0.5 Å. - """ - highest_hkl = get_highest_hkl( - lattice=phase.structure.lattice, min_dspacing=min_dspacing - ) - hkl = get_hkl(highest_hkl=highest_hkl) - return cls(phase=phase, hkl=hkl).unique() - - @classmethod - def from_highest_hkl(cls, phase, highest_hkl): - """Create a CrystalPlane object populated by unique Miller indices - below, but including, a set of higher indices. - - Parameters - ---------- - phase : orix.crystal_map.phase_list.Phase - A phase container with a crystal structure and a space and - point group describing the allowed symmetry operations. - highest_hkl : np.ndarray, list, or tuple of int - Highest Miller indices to consider (including). - """ - hkl = get_hkl(highest_hkl=highest_hkl) - return cls(phase=phase, hkl=hkl).unique() - - def unique(self, use_symmetry=True): - """Return planes with unique Miller indices. - - Parameters - ---------- - use_symmetry : bool, optional - Whether to use symmetry to remove the planes with indices - symmetrically equivalent to another set of indices. - - Returns - ------- - ReciprocalLatticePoint - """ - if use_symmetry: - all_hkl = self._hkldata - all_hkl = all_hkl[~np.all(np.isclose(all_hkl, 0), axis=1)] - families = defaultdict(list) - for this_hkl in all_hkl.tolist(): - for that_hkl in families.keys(): - if is_equivalent(this_hkl, that_hkl): - families[tuple(that_hkl)].append(this_hkl) - break - else: - families[tuple(this_hkl)].append(this_hkl) - - n_families = len(families) - unique_hkl = np.zeros((n_families, 3)) - for i, all_hkl_in_family in enumerate(families.values()): - unique_hkl[i] = sorted(all_hkl_in_family)[-1] - else: - unique_hkl = self.hkl.unique() - # TODO: Enable inheriting classes pass on their properties in this new object - return self.__class__(phase=self.phase, hkl=unique_hkl) - - def symmetrise( - self, antipodal=True, unique=True, return_multiplicity=False, - ): - """Return planes with symmetrically equivalent Miller indices. - - Parameters - ---------- - antipodal : bool, optional - Whether to include antipodal symmetry operations. Default is - True. - unique : bool, optional - Whether to return only distinct indices. Default is True. - If true, zero entries which are assumed to be degenerate are - removed. - return_multiplicity : bool, optional - Whether to return the multiplicity of indices. This option is - only available if `unique` is True. Default is False. - - Returns - ------- - ReciprocalLatticePoint - Planes with Miller indices symmetrically equivalent to the - original planes. - multiplicity : np.ndarray - Multiplicity of the original Miller indices. Only returned if - `return_multiplicity` is True. - - Notes - ----- - Should be the same as EMsoft's CalcFamily in their symmetry.f90 - module. - """ - # Get symmetry operations - pg = self.phase.point_group - operations = pg[~pg.improper] if not antipodal else pg - - out = get_equivalent_hkl( - hkl=self.hkl, - operations=operations, - unique=unique, - return_multiplicity=return_multiplicity, - ) - - # TODO: Enable inheriting classes pass on their properties in this new object - # Format output and return - if unique and return_multiplicity: - multiplicity = out[1] - if multiplicity.size == 1: - multiplicity = multiplicity[0] - return self.__class__(phase=self.phase, hkl=out[0]), multiplicity - else: - return self.__class__(phase=self.phase, hkl=out) - - def calculate_structure_factor(self, method=None, voltage=None): - """Populate `self.structure_factor` with the structure factor F - for each plane. - - Parameters - ---------- - method : str, optional - Either "kinematical" for kinematical X-ray structure factors - or "doyleturner" for structure factors using Doyle-Turner - atomic scattering factors. If None (default), kinematical - structure factors are calculated. - voltage : float, optional - Beam energy in V used when `method=doyleturner`. - """ - if method is None: - method = "kinematical" - methods = ["kinematical", "doyleturner"] - if method not in methods: - raise ValueError(f"method={method} must be among {methods}") - elif method == "doyleturner" and voltage is None: - raise ValueError( - "'voltage' parameter must set when method='doyleturner'" - ) - - structure_factors = np.zeros(self.size) - hkls = self._hkldata - scattering_parameters = self.scattering_parameter - phase = self.phase - # TODO: Find a better way to call different methods in the loop - for i, (hkl, s) in enumerate(zip(hkls, scattering_parameters)): - if method == "kinematical": - structure_factors[i] = get_kinematical_structure_factor( - phase=phase, hkl=hkl, scattering_parameter=s, - ) - else: - structure_factors[i] = get_doyleturner_structure_factor( - phase=phase, - hkl=hkl, - scattering_parameter=s, - voltage=voltage, - ) - self._structure_factor = np.where( - structure_factors < _FLOAT_EPS, 0, structure_factors - ) - - def calculate_theta(self, voltage): - """Populate `self.theta` with the Bragg angle :math:`theta_B` for - each plane. - - Parameters - ---------- - voltage : float - Beam energy in V. - """ - wavelength = get_refraction_corrected_wavelength(self.phase, voltage) - self._theta = np.arcsin(0.5 * wavelength * self.gspacing) - - -def get_highest_hkl(lattice, min_dspacing=0.5): - """Return the highest Miller indices hkl of the plane with a direct - space interplanar spacing greater than but closest to a lower - threshold. - - Parameters - ---------- - lattice : diffpy.structure.Lattice - Crystal lattice. - min_dspacing : float, optional - Smallest interplanar spacing to consider. Default is 0.5 Å. - - Returns - ------- - highest_hkl : np.ndarray - Highest Miller indices. - """ - highest_hkl = np.ones(3, dtype=int) - for i in range(3): - hkl = np.zeros(3) - d = min_dspacing + 1 - while d > min_dspacing: - hkl[i] += 1 - d = 1 / lattice.rnorm(hkl) - highest_hkl[i] = hkl[i] - return highest_hkl - - -def get_hkl(highest_hkl): - """Return a list of planes from a set of highest Miller indices. - - Parameters - ---------- - highest_hkl : orix.vector.Vector3d, np.ndarray, list, or tuple of int - Highest Miller indices to consider. - - Returns - ------- - hkl : np.ndarray - An array of Miller indices. - """ - index_ranges = [np.arange(-i, i + 1) for i in highest_hkl] - return np.asarray(list(product(*index_ranges))) - - -def get_equivalent_hkl( - hkl, operations, unique=False, return_multiplicity=False -): - """Return symmetrically equivalent Miller indices. - - Parameters - ---------- - hkl : orix.vector.Vector3d, np.ndarray, list or tuple of int - Miller indices. - operations : orix.quaternion.symmetry.Symmetry - Point group describing allowed symmetry operations. - unique : bool, optional - Whether to return only unique Miller indices. Default is False. - return_multiplicity : bool, optional - Whether to return the multiplicity of the input indices. Default - is False. - - Returns - ------- - new_hkl : orix.vector.Vector3d - The symmetrically equivalent Miller indices. - multiplicity : np.ndarray - Number of symmetrically equivalent indices. Only returned if - `return_multiplicity` is True. - """ - new_hkl = operations.outer(Vector3d(hkl)) - new_hkl = new_hkl.flatten().reshape(*new_hkl.shape[::-1]) - - multiplicity = None - if unique: - n_families = new_hkl.shape[0] - multiplicity = np.zeros(n_families, dtype=int) - temp_hkl = new_hkl[0].unique().data - multiplicity[0] = temp_hkl.shape[0] - if n_families > 1: - for i, hkl in enumerate(new_hkl[1:]): - temp_hkl2 = hkl.unique() - multiplicity[i + 1] = temp_hkl2.size - temp_hkl = np.append(temp_hkl, temp_hkl2.data, axis=0) - new_hkl = Vector3d(temp_hkl[: multiplicity.sum()]) - - # Remove 1-dimensions - new_hkl = new_hkl.squeeze() - - if unique and return_multiplicity: - return new_hkl, multiplicity - else: - return new_hkl - - -def is_equivalent(this_hkl: list, that_hkl: list) -> bool: - return sorted(np.abs(this_hkl)) == sorted(np.abs(that_hkl)) diff --git a/kikuchipy/crystallography/structure_factor.py b/kikuchipy/crystallography/structure_factor.py deleted file mode 100644 index d9ff010e..00000000 --- a/kikuchipy/crystallography/structure_factor.py +++ /dev/null @@ -1,225 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright 2017-2020 The diffsims developers -# -# This file is part of diffsims. -# -# diffsims is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# diffsims is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with diffsims. If not, see . - -from diffpy.structure.symmetryutilities import ( - expandPosition, - SymmetryConstraints, -) -import numpy as np -from scipy.constants import c, e, h, physical_constants - -from kikuchipy.crystallography.atomic_scattering_factor import ( - get_kinematical_atomic_scattering_factor, - get_doyleturner_atomic_scattering_factor, -) - - -rest_mass = physical_constants["atomic unit of mass"][0] - - -def find_asymmetric_positions(positions, space_group): - """Return the asymmetric atom positions among a set of positions - when considering symmetry operations defined by a space group. - - Parameters - ---------- - positions : list - A list of cartesian atom positions. - space_group : diffpy.structure.spacegroups.SpaceGroup - Space group describing the symmetry operations. - - Returns - ------- - np.ndarray - Asymmetric atom positions. - """ - asymmetric_positions = SymmetryConstraints(space_group, positions).corepos - return [ - np.array([np.allclose(xyz, asym_xyz) for xyz in positions]) - for asym_xyz in asymmetric_positions - ][0] - - -def get_kinematical_structure_factor(phase, hkl, scattering_parameter): - """Return the kinematical (X-ray) structure factor for a given family - of Miller indices. - - Assumes structure's lattice parameters and Debye-Waller factors are - expressed in Ångströms. - - Parameters - ---------- - phase : orix.crystal_map.phase_list.Phase - A phase container with a crystal structure and a space and point - group describing the allowed symmetry operations. - hkl : np.ndarray - Miller indices. - scattering_parameter : float - Scattering parameter for these Miller indices. - - Returns - ------- - structure_factor : float - Structure factor F. - """ - # Initialize real and imaginary parts of the structure factor - structure_factor = 0 + 0j - - structure = phase.structure - space_group = phase.space_group - - # Loop over asymmetric unit - asymmetric_positions = find_asymmetric_positions(structure.xyz, space_group) - for is_asymmetric, atom in zip(asymmetric_positions, structure): - if not is_asymmetric: - continue - - # Get atomic scattering factor for this atom - f = get_kinematical_atomic_scattering_factor(atom, scattering_parameter) - - # Loop over all atoms in the orbit - equiv_pos = expandPosition(spacegroup=space_group, xyz=atom.xyz)[0] - for xyz in equiv_pos: - arg = 2 * np.pi * np.sum(hkl * xyz) - structure_factor += f * (np.cos(arg) - (np.sin(arg) * 1j)) - - return structure_factor.real - - -def get_doyleturner_structure_factor( - phase, hkl, scattering_parameter, voltage, return_parameters=False, -): - """Return the structure factor for a given family of Miller indices - using Doyle-Turner atomic scattering parameters [Doyle1968]_. - - Assumes structure's lattice parameters and Debye-Waller factors are - expressed in Ångströms. - - Parameters - ---------- - phase : orix.crystal_map.phase_list.Phase - A phase container with a crystal structure and a space and point - group describing the allowed symmetry operations. - hkl : np.ndarray - Miller indices. - scattering_parameter : float - Scattering parameter for these Miller indices. - voltage : float - Beam energy in V. - return_parameters : bool, optional - Whether to return a set of parameters derived from the - calculation as a dictionary. Default is False. - - Returns - ------- - structure_factor : float - Structure factor F. - params : dict - A dictionary with (key, item) (str, float) of parameters derived - from the calculation. Only returned if `return_parameters=True`. - """ - structure = phase.structure - space_group = phase.space_group - - # Initialize real and imaginary parts of the structure factor - structure_factor = 0 + 0j - - # Get unit cell volume for the atomic scattering factor calculation - unit_cell_volume = structure.lattice.volume - - # Loop over all atoms in the asymmetric unit - asymmetric_positions = find_asymmetric_positions(structure.xyz, space_group) - for is_asymmetric, atom in zip(asymmetric_positions, structure): - if not is_asymmetric: - continue - - # Get atomic scattering factor for this atom - f = get_doyleturner_atomic_scattering_factor( - atom, scattering_parameter, unit_cell_volume - ) - - # Loop over all atoms in the orbit - equiv_pos = expandPosition(spacegroup=space_group, xyz=atom.xyz)[0] - for xyz in equiv_pos: - arg = 2 * np.pi * np.sum(hkl * xyz) - structure_factor += f * np.exp(-arg * 1j) - - # Derived parameters - # TODO: Comment these factors with stuff from Structure of Materials by De Graef - # and McHenry - gamma_relcor = 1 + (2 * e * 0.5 * voltage / rest_mass / (c ** 2)) - v_mod = abs(structure_factor) * gamma_relcor - v_phase = np.arctan2(structure_factor.imag, structure_factor.real) - v_g = v_mod * np.exp(v_phase * 1j) - pre = 2 * (rest_mass * e / h ** 2) * 1e-18 - - structure_factor = (pre * v_g).real - - if return_parameters: - params = { - "gamma_relcor": gamma_relcor, - "v_mod": v_mod, - "v_phase": v_phase, - "v_g": v_g, - } - return structure_factor, params - else: - return structure_factor - - -def get_refraction_corrected_wavelength(phase, voltage): - """Return the refraction corrected relativistic electron wavelength - in Ångströms for a given crystal structure and beam energy in V. - - Parameters - ---------- - phase : orix.crystal_map.phase_list.Phase - A phase container with a crystal structure and a space and point - group describing the allowed symmetry operations. - voltage : float - Beam energy in V. - - Returns - ------- - wavelength : float - Refraction corrected relativistic electron wavelength in - Ångströms. - """ - temp1 = 1e9 * h / np.sqrt(2 * rest_mass * e) - temp2 = e * 0.5 * voltage / rest_mass / (c ** 2) - - # Relativistic correction factor (known as gamma) - # gamma_relcor = 1 + (2 * temp2) - - # Relativistic acceleration voltage - psi_hat = voltage * (1 + temp2) - - # Compute the electron wavelength in nm - hkl = np.zeros(3, dtype=int) - scattering_parameter = 0 - _, params = get_doyleturner_structure_factor( - phase, hkl, scattering_parameter, voltage, return_parameters=True - ) - v_mod = params["v_mod"] - psi_hat += v_mod - wavelength = temp1 / np.sqrt(psi_hat) - - # Interaction constant sigma - # sigma = 1e-18 * (2 * np.pi * rest_mass * gamma_relcor * e * wavelength) / h ** 2 - - return wavelength diff --git a/kikuchipy/generators/ebsd_simulation_generator.py b/kikuchipy/generators/ebsd_simulation_generator.py index e6679b02..8a7d7fbf 100644 --- a/kikuchipy/generators/ebsd_simulation_generator.py +++ b/kikuchipy/generators/ebsd_simulation_generator.py @@ -19,11 +19,11 @@ from copy import deepcopy from typing import Optional +from diffsims.crystallography import ReciprocalLatticePoint import numpy as np from orix.crystal_map import Phase from orix.quaternion.rotation import Rotation -from kikuchipy.crystallography import ReciprocalLatticePoint from kikuchipy.detectors import EBSDDetector from kikuchipy.projections.ebsd_projections import ( detector2reciprocal_lattice, diff --git a/kikuchipy/simulations/features.py b/kikuchipy/simulations/features.py index 5499b6d5..0f1da7ee 100644 --- a/kikuchipy/simulations/features.py +++ b/kikuchipy/simulations/features.py @@ -18,12 +18,11 @@ from typing import Union +from diffsims.crystallography import ReciprocalLatticePoint import numpy as np from orix.crystal_map import Phase from orix.vector import Vector3d -from kikuchipy.crystallography import ReciprocalLatticePoint - class KikuchiBand(ReciprocalLatticePoint): def __init__( diff --git a/setup.py b/setup.py index d1797585..1c957c86 100644 --- a/setup.py +++ b/setup.py @@ -114,6 +114,7 @@ extras_require=extra_feature_requirements, install_requires=[ "dask[array] >= 2.14", + "diffsims >= 0.3", "hyperspy >= 1.5.2", "h5py >= 2.10", "matplotlib >= 3.2",