From 54ae57be81bd62a17cb2bb6aec3b2d36751b1a6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20Wiik=20=C3=85nes?= Date: Sat, 22 Aug 2020 17:37:30 +0200 Subject: [PATCH 1/8] Add reciprocal_lattice_point module w/ rlp class, myself to credits MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Håkon Wiik Ånes --- diffsims/reciprocal_lattice_point/__init__.py | 47 +++ .../atomic_scattering_parameters.py | 208 ++++++++++ .../reciprocal_lattice_point.py | 358 ++++++++++++++++++ .../structure_factor.py | 128 +++++++ diffsims/release_info.py | 8 +- 5 files changed, 748 insertions(+), 1 deletion(-) create mode 100644 diffsims/reciprocal_lattice_point/__init__.py create mode 100644 diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py create mode 100644 diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py create mode 100644 diffsims/reciprocal_lattice_point/structure_factor.py diff --git a/diffsims/reciprocal_lattice_point/__init__.py b/diffsims/reciprocal_lattice_point/__init__.py new file mode 100644 index 00000000..5208e597 --- /dev/null +++ b/diffsims/reciprocal_lattice_point/__init__.py @@ -0,0 +1,47 @@ +# -*- 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 . + +"""Generation of HKL reflectors for a crystal structure.""" + +from diffsims.reciprocal_lattice_point.atomic_scattering_parameters import ( + get_atomic_scattering_parameters, + get_element_id_from_string, +) +from diffsims.reciprocal_lattice_point.reciprocal_lattice_point import ( + ReciprocalLatticePoint, + get_equivalent_hkl, + get_highest_hkl, + get_hkl, +) +from diffsims.reciprocal_lattice_point.structure_factor import ( + find_asymmetric_positions, + get_atomic_scattering_factor, + get_xray_structure_factor, +) + +__all__ = [ + "get_atomic_scattering_parameters", + "get_element_id_from_string", + "ReciprocalLatticePoint", + "get_equivalent_hkl", + "get_highest_hkl", + "get_hkl", + "find_asymmetric_positions", + "get_atomic_scattering_factor", + "get_xray_structure_factor", +] diff --git a/diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py b/diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py new file mode 100644 index 00000000..bda75c93 --- /dev/null +++ b/diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py @@ -0,0 +1,208 @@ +# -*- 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 . + +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/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py b/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py new file mode 100644 index 00000000..bb6f4fe7 --- /dev/null +++ b/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py @@ -0,0 +1,358 @@ +# -*- 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 collections import defaultdict +from itertools import product + +import numpy as np +from orix.vector import Vector3d + +from diffsims.reciprocal_lattice_point.structure_factor import ( + get_xray_structure_factor +) + + +class ReciprocalLatticePoint: + """Reciprocal lattice points (reflectors) g 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 reciprocal lattice points (reflectors) g. + + 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 + + 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_rlp = self.__class__(self.phase, self.hkl[key]) + new_rlp._structure_factor = self.structure_factor[key] + return new_rlp + + @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 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 + + @classmethod + def from_min_dspacing(cls, phase, min_dspacing=0.5): + """Create a ReciprocalLatticePoint 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 ReciprocalLatticePoint 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() + + @classmethod + def from_nfamilies(cls, phase, nfamilies=5): + raise NotImplementedError + + def calculate_structure_factor(self, method=None): + """Populate `self.structure_factor` with the structure factor F + for each point. + + Parameters + ---------- + method + Either "xray" for the X-ray structure factor or + "doyleturner" for the structure factor using Doyle-Turner + atomic scattering factors. + """ + structure_factors = np.zeros(self.size) + hkls = self._hkldata + scattering_parameters = self.scattering_parameter + for i, (hkl, s) in enumerate(zip(hkls, scattering_parameters)): + structure_factors[i] = get_xray_structure_factor( + phase=self.phase, hkl=hkl, scattering_parameter=s + ) + self._structure_factor = structure_factors + + def unique(self, use_symmetry=True): + """Return reciprocal lattice points with unique Miller indices. + + Parameters + ---------- + use_symmetry : bool, optional + Whether to use symmetry to remove the 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() + return self.__class__(phase=self.phase, hkl=unique_hkl) + + def symmetrise( + self, antipodal=True, unique=True, return_multiplicity=False, + ): + """Return reciprocal lattice points 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 the indices. This + option is only available if `unique` is True. Default is + False. + + Returns + ------- + ReciprocalLatticePoint + Reciprocal lattice points with Miller indices symmetrically + equivalent to the original lattice points. + 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, + ) + + # 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 get_highest_hkl(lattice, min_dspacing=0.5): + """Return the highest Miller indices hkl of the reciprocal + lattice point with a direct space interplanar spacing greater + than but closest to a lower threshold. + + Parameters + ---------- + lattice : diffpy.structure.Lattice + Crystal structure 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 reciprocal lattice points 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 + An array of reciprocal lattice points. + """ + 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 hkl. 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/diffsims/reciprocal_lattice_point/structure_factor.py b/diffsims/reciprocal_lattice_point/structure_factor.py new file mode 100644 index 00000000..869887c8 --- /dev/null +++ b/diffsims/reciprocal_lattice_point/structure_factor.py @@ -0,0 +1,128 @@ +# -*- 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 diffsims.reciprocal_lattice_point.atomic_scattering_parameters import ( + get_atomic_scattering_parameters, + get_element_id_from_string, +) + + +def get_atomic_scattering_factor(atom, scattering_parameter): + """Return the atomic scattering factor f for a certain atom and + scattering parameter. + + 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 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_xray_structure_factor(phase, hkl, scattering_parameter): + """Assumes structure's lattice parameters and Debye-Waller factors + are expressed in Angstroms. + + 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_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 diff --git a/diffsims/release_info.py b/diffsims/release_info.py index 5c35c942..bc9b1382 100644 --- a/diffsims/release_info.py +++ b/diffsims/release_info.py @@ -2,7 +2,13 @@ version = "0.2.3" author = "Duncan Johnstone, Phillip Crout" copyright = "Copyright 2017-2020, The pyXem Developers" -credits = ["Duncan Johnstone", "Phillip Crout", "Ben Martineau", "Simon Hogas"] +credits = [ + "Duncan Johnstone", + "Phillip Crout", + "Ben Martineau", + "Simon Hogas", + "Håkon Wiik Ånes", +] license = "GPLv3" maintainer = "Duncan Johnstone, Phillip Crout" email = "pyxem.team@gmail.com" From a5d776e2fca7645f053949a269ff454cf2627ecb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20Wiik=20=C3=85nes?= Date: Sun, 23 Aug 2020 18:00:15 +0200 Subject: [PATCH 2/8] Add atomic scattering factor and structure factor Doyle-Turner calc MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Håkon Wiik Ånes --- diffsims/reciprocal_lattice_point/__init__.py | 14 ++- .../atomic_scattering_factor.py | 101 +++++++++++++++ .../atomic_scattering_parameters.py | 4 +- .../reciprocal_lattice_point.py | 48 +++++-- .../structure_factor.py | 118 +++++++++++------- 5 files changed, 222 insertions(+), 63 deletions(-) create mode 100644 diffsims/reciprocal_lattice_point/atomic_scattering_factor.py diff --git a/diffsims/reciprocal_lattice_point/__init__.py b/diffsims/reciprocal_lattice_point/__init__.py index 5208e597..9e09680f 100644 --- a/diffsims/reciprocal_lattice_point/__init__.py +++ b/diffsims/reciprocal_lattice_point/__init__.py @@ -30,8 +30,12 @@ ) from diffsims.reciprocal_lattice_point.structure_factor import ( find_asymmetric_positions, - get_atomic_scattering_factor, - get_xray_structure_factor, + get_kinematical_structure_factor, + get_doyleturner_structure_factor, +) +from diffsims.reciprocal_lattice_point.atomic_scattering_factor import ( + get_kinematical_atomic_scattering_factor, + get_doyleturner_atomic_scattering_factor, ) __all__ = [ @@ -42,6 +46,8 @@ "get_highest_hkl", "get_hkl", "find_asymmetric_positions", - "get_atomic_scattering_factor", - "get_xray_structure_factor", + "get_kinematical_atomic_scattering_factor", + "get_doyleturner_atomic_scattering_factor", + "get_kinematical_structure_factor", + "get_doyleturner_structure_factor", ] diff --git a/diffsims/reciprocal_lattice_point/atomic_scattering_factor.py b/diffsims/reciprocal_lattice_point/atomic_scattering_factor.py new file mode 100644 index 00000000..58f5a76d --- /dev/null +++ b/diffsims/reciprocal_lattice_point/atomic_scattering_factor.py @@ -0,0 +1,101 @@ +# -*- 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 . + +import numpy as np + +from diffsims.reciprocal_lattice_point.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/diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py b/diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py index bda75c93..a27d6ef8 100644 --- a/diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py +++ b/diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py @@ -201,8 +201,6 @@ def get_element_id_from_string(element_str): ] # fmt: on n_elements = len(elements) - element2periodic = dict( - zip(elements[:n_elements], np.arange(1, n_elements)) - ) + element2periodic = dict(zip(elements[:n_elements], np.arange(1, n_elements))) element_id = element2periodic[element_str.lower()] return element_id diff --git a/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py b/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py index bb6f4fe7..05648feb 100644 --- a/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py +++ b/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py @@ -23,10 +23,14 @@ from orix.vector import Vector3d from diffsims.reciprocal_lattice_point.structure_factor import ( - get_xray_structure_factor + get_kinematical_structure_factor, + get_doyleturner_structure_factor, ) +_FLOAT_EPS = np.finfo(np.float).eps # Used to round values below 1e-16 to zero + + class ReciprocalLatticePoint: """Reciprocal lattice points (reflectors) g with Miller indices, length of the reciprocal lattice vectors and other relevant @@ -150,25 +154,45 @@ def from_highest_hkl(cls, phase, highest_hkl): def from_nfamilies(cls, phase, nfamilies=5): raise NotImplementedError - def calculate_structure_factor(self, method=None): + def calculate_structure_factor(self, method=None, voltage=None): """Populate `self.structure_factor` with the structure factor F for each point. Parameters ---------- - method - Either "xray" for the X-ray structure factor or - "doyleturner" for the structure factor using Doyle-Turner - atomic scattering factors. + 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 voltage 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)): - structure_factors[i] = get_xray_structure_factor( - phase=self.phase, hkl=hkl, scattering_parameter=s - ) - self._structure_factor = structure_factors + 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 unique(self, use_symmetry=True): """Return reciprocal lattice points with unique Miller indices. @@ -304,9 +328,7 @@ def get_hkl(highest_hkl): return np.asarray(list(product(*index_ranges))) -def get_equivalent_hkl( - hkl, operations, unique=False, return_multiplicity=False -): +def get_equivalent_hkl(hkl, operations, unique=False, return_multiplicity=False): """Return symmetrically equivalent Miller indices. Parameters diff --git a/diffsims/reciprocal_lattice_point/structure_factor.py b/diffsims/reciprocal_lattice_point/structure_factor.py index 869887c8..ee0fe486 100644 --- a/diffsims/reciprocal_lattice_point/structure_factor.py +++ b/diffsims/reciprocal_lattice_point/structure_factor.py @@ -16,50 +16,17 @@ # 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, -) +from diffpy.structure.symmetryutilities import expandPosition, SymmetryConstraints import numpy as np +from scipy.constants import c, e, h, physical_constants -from diffsims.reciprocal_lattice_point.atomic_scattering_parameters import ( - get_atomic_scattering_parameters, - get_element_id_from_string, +from diffsims.reciprocal_lattice_point.atomic_scattering_factor import ( + get_kinematical_atomic_scattering_factor, + get_doyleturner_atomic_scattering_factor, ) -def get_atomic_scattering_factor(atom, scattering_parameter): - """Return the atomic scattering factor f for a certain atom and - scattering parameter. - - 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 +rest_mass = physical_constants["atomic unit of mass"][0] def find_asymmetric_positions(positions, space_group): @@ -85,9 +52,12 @@ def find_asymmetric_positions(positions, space_group): ][0] -def get_xray_structure_factor(phase, hkl, scattering_parameter): - """Assumes structure's lattice parameters and Debye-Waller factors - are expressed in Angstroms. +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 ---------- @@ -117,7 +87,7 @@ def get_xray_structure_factor(phase, hkl, scattering_parameter): continue # Get atomic scattering factor for this atom - f = get_atomic_scattering_factor(atom, scattering_parameter) + 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] @@ -126,3 +96,65 @@ def get_xray_structure_factor(phase, hkl, scattering_parameter): 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 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. + + Returns + ------- + structure_factor : float + Structure factor F. + """ + 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 factors + # 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 + + return (pre * v_g).real From 8c92bdf5fb55a4076846f42795d16d5062e0c625 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20Wiik=20=C3=85nes?= Date: Sun, 23 Aug 2020 19:29:53 +0200 Subject: [PATCH 3/8] Add allowed hkl without hexagonal rules, quick acces to h, k and l MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Håkon Wiik Ånes --- .../reciprocal_lattice_point.py | 46 +++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py b/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py index 05648feb..095f4f63 100644 --- a/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py +++ b/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py @@ -75,6 +75,21 @@ 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`.""" @@ -114,6 +129,37 @@ def structure_factor(self): """Return :class:`np.ndarray` of structure factors F or None.""" return self._structure_factor + @property + def allowed(self): + """Return whether reciprocal lattice points 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 + @classmethod def from_min_dspacing(cls, phase, min_dspacing=0.5): """Create a ReciprocalLatticePoint object populated by unique From 09098c85d036e45ce1dcb3a0f8f4aecb85db322a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20Wiik=20=C3=85nes?= Date: Mon, 24 Aug 2020 15:04:37 +0200 Subject: [PATCH 4/8] Rename RLP class to CrystalPlane, restructure into two modules MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Håkon Wiik Ånes --- diffsims/crystallography/__init__.py | 33 ++++ .../crystal_plane.py} | 174 ++++++++++-------- .../__init__.py | 34 ++-- .../atomic_scattering_factor.py | 2 +- .../atomic_scattering_parameters.py | 0 .../structure_factor.py | 70 ++++++- 6 files changed, 208 insertions(+), 105 deletions(-) create mode 100644 diffsims/crystallography/__init__.py rename diffsims/{reciprocal_lattice_point/reciprocal_lattice_point.py => crystallography/crystal_plane.py} (83%) rename diffsims/{reciprocal_lattice_point => diffraction}/__init__.py (71%) rename diffsims/{reciprocal_lattice_point => diffraction}/atomic_scattering_factor.py (97%) rename diffsims/{reciprocal_lattice_point => diffraction}/atomic_scattering_parameters.py (100%) rename diffsims/{reciprocal_lattice_point => diffraction}/structure_factor.py (71%) diff --git a/diffsims/crystallography/__init__.py b/diffsims/crystallography/__init__.py new file mode 100644 index 00000000..c5d5b7a5 --- /dev/null +++ b/diffsims/crystallography/__init__.py @@ -0,0 +1,33 @@ +# -*- 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 . + +"""Generation of crystal planes hkl for a crystal structure.""" + +from diffsims.crystallography.crystal_plane import ( + CrystalPlane, + get_equivalent_hkl, + get_highest_hkl, + get_hkl, +) + +__all__ = [ + "CrystalPlane", + "get_equivalent_hkl", + "get_highest_hkl", + "get_hkl", +] diff --git a/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py b/diffsims/crystallography/crystal_plane.py similarity index 83% rename from diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py rename to diffsims/crystallography/crystal_plane.py index 095f4f63..090f48a4 100644 --- a/diffsims/reciprocal_lattice_point/reciprocal_lattice_point.py +++ b/diffsims/crystallography/crystal_plane.py @@ -22,24 +22,26 @@ import numpy as np from orix.vector import Vector3d -from diffsims.reciprocal_lattice_point.structure_factor import ( +from diffsims.diffraction.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 points (reflectors) g with Miller indices, - length of the reciprocal lattice vectors and other relevant - diffraction parameters. +class CrystalPlane: + """Crystal plane (or reciprocal lattice point, reflectors, 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 reciprocal lattice points (reflectors) g. + parameters for crystal planes (reciprocal lattice points, + reflectors, g, etc.). Parameters ---------- @@ -52,6 +54,7 @@ def __init__(self, phase, hkl): self._hkl = Vector3d(hkl) self.phase = phase self._structure_factor = [None] * self.size + self._theta = [None] * self.size def __repr__(self): return ( @@ -61,9 +64,9 @@ def __repr__(self): ) def __getitem__(self, key): - new_rlp = self.__class__(self.phase, self.hkl[key]) - new_rlp._structure_factor = self.structure_factor[key] - return new_rlp + new_cp = self.__class__(self.phase, self.hkl[key]) + new_cp._structure_factor = self.structure_factor[key] + return new_cp @property def hkl(self): @@ -131,9 +134,8 @@ def structure_factor(self): @property def allowed(self): - """Return whether reciprocal lattice points diffract according to - diffraction selection rules assuming kinematical scattering - theory. + """Return whether planes diffract according to diffraction + selection rules assuming kinematical scattering theory. """ # Translational symmetry centering = self.phase.space_group.short_name[0] @@ -160,11 +162,16 @@ def allowed(self): 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 ReciprocalLatticePoint object populated by unique - Miller indices with a direct space interplanar spacing greater - than a lower threshold. + """Create a CrystalPlane object populated by unique Miller indices + with a direct space interplanar spacing greater than a lower + threshold. Parameters ---------- @@ -182,8 +189,8 @@ def from_min_dspacing(cls, phase, min_dspacing=0.5): @classmethod def from_highest_hkl(cls, phase, highest_hkl): - """Create a ReciprocalLatticePoint object populated by unique - Miller indices below, but including, a set of higher indices. + """Create a CrystalPlane object populated by unique Miller indices + below, but including, a set of higher indices. Parameters ---------- @@ -200,58 +207,18 @@ def from_highest_hkl(cls, phase, highest_hkl): def from_nfamilies(cls, phase, nfamilies=5): raise NotImplementedError - def calculate_structure_factor(self, method=None, voltage=None): - """Populate `self.structure_factor` with the structure factor F - for each point. - - 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 voltage 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 unique(self, use_symmetry=True): - """Return reciprocal lattice points with unique Miller indices. + """Return planes with unique Miller indices. Parameters ---------- use_symmetry : bool, optional - Whether to use symmetry to remove the indices symmetrically - equivalent to another set of indices. + Whether to use symmetry to remove the planes with indices + symmetrically equivalent to another set of indices. Returns ------- - ReciprocalLatticePoint + CrystalPlane """ if use_symmetry: all_hkl = self._hkldata @@ -276,8 +243,7 @@ def unique(self, use_symmetry=True): def symmetrise( self, antipodal=True, unique=True, return_multiplicity=False, ): - """Return reciprocal lattice points with symmetrically equivalent - Miller indices. + """Return planes with symmetrically equivalent Miller indices. Parameters ---------- @@ -289,15 +255,14 @@ def symmetrise( If true, zero entries which are assumed to be degenerate are removed. return_multiplicity : bool, optional - Whether to return the multiplicity of the indices. This - option is only available if `unique` is True. Default is - False. + Whether to return the multiplicity of indices. This option is + only available if `unique` is True. Default is False. Returns ------- - ReciprocalLatticePoint - Reciprocal lattice points with Miller indices symmetrically - equivalent to the original lattice points. + CrystalPlane + 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. @@ -327,16 +292,68 @@ def symmetrise( 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 reciprocal - lattice point with a direct space interplanar spacing greater - than but closest to a lower threshold. + """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 structure lattice. + Crystal lattice. min_dspacing : float, optional Smallest interplanar spacing to consider. Default is 0.5 Å. @@ -357,8 +374,7 @@ def get_highest_hkl(lattice, min_dspacing=0.5): def get_hkl(highest_hkl): - """Return a list of reciprocal lattice points from a set of highest - Miller indices. + """Return a list of planes from a set of highest Miller indices. Parameters ---------- @@ -367,8 +383,8 @@ def get_hkl(highest_hkl): Returns ------- - hkl - An array of reciprocal lattice points. + 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))) @@ -386,8 +402,8 @@ def get_equivalent_hkl(hkl, operations, unique=False, return_multiplicity=False) 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 hkl. Default is - False. + Whether to return the multiplicity of the input indices. Default + is False. Returns ------- diff --git a/diffsims/reciprocal_lattice_point/__init__.py b/diffsims/diffraction/__init__.py similarity index 71% rename from diffsims/reciprocal_lattice_point/__init__.py rename to diffsims/diffraction/__init__.py index 9e09680f..31df42bc 100644 --- a/diffsims/reciprocal_lattice_point/__init__.py +++ b/diffsims/diffraction/__init__.py @@ -16,38 +16,30 @@ # You should have received a copy of the GNU General Public License # along with diffsims. If not, see . -"""Generation of HKL reflectors for a crystal structure.""" +"""Calculation of scattering factors and structure factors.""" -from diffsims.reciprocal_lattice_point.atomic_scattering_parameters import ( +from diffsims.diffraction.atomic_scattering_factor import ( + get_doyleturner_atomic_scattering_factor, + get_kinematical_atomic_scattering_factor, +) +from diffsims.diffraction.atomic_scattering_parameters import ( get_atomic_scattering_parameters, get_element_id_from_string, ) -from diffsims.reciprocal_lattice_point.reciprocal_lattice_point import ( - ReciprocalLatticePoint, - get_equivalent_hkl, - get_highest_hkl, - get_hkl, -) -from diffsims.reciprocal_lattice_point.structure_factor import ( +from diffsims.diffraction.structure_factor import ( find_asymmetric_positions, - get_kinematical_structure_factor, get_doyleturner_structure_factor, -) -from diffsims.reciprocal_lattice_point.atomic_scattering_factor import ( - get_kinematical_atomic_scattering_factor, - get_doyleturner_atomic_scattering_factor, + get_kinematical_structure_factor, + get_refraction_corrected_wavelength, ) __all__ = [ + "get_doyleturner_atomic_scattering_factor", + "get_kinematical_atomic_scattering_factor", "get_atomic_scattering_parameters", "get_element_id_from_string", - "ReciprocalLatticePoint", - "get_equivalent_hkl", - "get_highest_hkl", - "get_hkl", "find_asymmetric_positions", - "get_kinematical_atomic_scattering_factor", - "get_doyleturner_atomic_scattering_factor", - "get_kinematical_structure_factor", "get_doyleturner_structure_factor", + "get_kinematical_structure_factor", + "get_refraction_corrected_wavelength", ] diff --git a/diffsims/reciprocal_lattice_point/atomic_scattering_factor.py b/diffsims/diffraction/atomic_scattering_factor.py similarity index 97% rename from diffsims/reciprocal_lattice_point/atomic_scattering_factor.py rename to diffsims/diffraction/atomic_scattering_factor.py index 58f5a76d..d4b56352 100644 --- a/diffsims/reciprocal_lattice_point/atomic_scattering_factor.py +++ b/diffsims/diffraction/atomic_scattering_factor.py @@ -18,7 +18,7 @@ import numpy as np -from diffsims.reciprocal_lattice_point.atomic_scattering_parameters import ( +from diffsims.diffraction.atomic_scattering_parameters import ( get_atomic_scattering_parameters, get_element_id_from_string, ) diff --git a/diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py b/diffsims/diffraction/atomic_scattering_parameters.py similarity index 100% rename from diffsims/reciprocal_lattice_point/atomic_scattering_parameters.py rename to diffsims/diffraction/atomic_scattering_parameters.py diff --git a/diffsims/reciprocal_lattice_point/structure_factor.py b/diffsims/diffraction/structure_factor.py similarity index 71% rename from diffsims/reciprocal_lattice_point/structure_factor.py rename to diffsims/diffraction/structure_factor.py index ee0fe486..54ab35d8 100644 --- a/diffsims/reciprocal_lattice_point/structure_factor.py +++ b/diffsims/diffraction/structure_factor.py @@ -20,7 +20,7 @@ import numpy as np from scipy.constants import c, e, h, physical_constants -from diffsims.reciprocal_lattice_point.atomic_scattering_factor import ( +from diffsims.diffraction.atomic_scattering_factor import ( get_kinematical_atomic_scattering_factor, get_doyleturner_atomic_scattering_factor, ) @@ -98,7 +98,9 @@ def get_kinematical_structure_factor(phase, hkl, scattering_parameter): return structure_factor.real -def get_doyleturner_structure_factor(phase, hkl, scattering_parameter, voltage): +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]_. @@ -116,11 +118,17 @@ def get_doyleturner_structure_factor(phase, hkl, scattering_parameter, voltage): 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 @@ -148,7 +156,7 @@ def get_doyleturner_structure_factor(phase, hkl, scattering_parameter, voltage): arg = 2 * np.pi * np.sum(hkl * xyz) structure_factor += f * np.exp(-arg * 1j) - # Derived factors + # 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)) @@ -157,4 +165,58 @@ def get_doyleturner_structure_factor(phase, hkl, scattering_parameter, voltage): v_g = v_mod * np.exp(v_phase * 1j) pre = 2 * (rest_mass * e / h ** 2) * 1e-18 - return (pre * v_g).real + 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 From 3f656c71c0c4f71ffef3c63e367c79bfc2dfa88f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20Wiik=20=C3=85nes?= Date: Tue, 8 Sep 2020 19:39:20 +0200 Subject: [PATCH 5/8] Rename to ReciprocalLatticePoint, other minor changes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Håkon Wiik Ånes --- diffsims/__init__.py | 6 +-- diffsims/crystallography/__init__.py | 10 +++-- ...l_plane.py => reciprocal_lattice_point.py} | 43 +++++++++++-------- .../__init__.py | 6 +-- .../atomic_scattering_factor.py | 6 ++- .../atomic_scattering_parameters.py | 0 .../structure_factor.py | 14 +++++- 7 files changed, 54 insertions(+), 31 deletions(-) rename diffsims/crystallography/{crystal_plane.py => reciprocal_lattice_point.py} (92%) rename diffsims/{diffraction => structure_factor}/__init__.py (88%) rename diffsims/{diffraction => structure_factor}/atomic_scattering_factor.py (95%) rename diffsims/{diffraction => structure_factor}/atomic_scattering_parameters.py (100%) rename diffsims/{diffraction => structure_factor}/structure_factor.py (96%) diff --git a/diffsims/__init__.py b/diffsims/__init__.py index c741ab77..d55827f9 100644 --- a/diffsims/__init__.py +++ b/diffsims/__init__.py @@ -17,10 +17,6 @@ # along with diffsims. If not, see . import logging -import os -import warnings - -import numpy as np from .generators.diffraction_generator import ( DiffractionGenerator, @@ -44,6 +40,8 @@ from .utils.atomic_diffraction_generator_support.discretise_utils import ( get_discretisation, ) +from .crystallography import * # What's imported is specified in the modules' init +from .structure_factor import * # What's imported is specified in the modules' init from . import release_info diff --git a/diffsims/crystallography/__init__.py b/diffsims/crystallography/__init__.py index c5d5b7a5..46372805 100644 --- a/diffsims/crystallography/__init__.py +++ b/diffsims/crystallography/__init__.py @@ -16,17 +16,19 @@ # You should have received a copy of the GNU General Public License # along with diffsims. If not, see . -"""Generation of crystal planes hkl for a crystal structure.""" +"""Generation of reciprocal lattice points (crystal plane, reflector, g, +hkl) for a crystal structure. +""" -from diffsims.crystallography.crystal_plane import ( - CrystalPlane, +from diffsims.crystallography.reciprocal_lattice_point import ( + ReciprocalLatticePoint, get_equivalent_hkl, get_highest_hkl, get_hkl, ) __all__ = [ - "CrystalPlane", + "ReciprocalLatticePoint", "get_equivalent_hkl", "get_highest_hkl", "get_hkl", diff --git a/diffsims/crystallography/crystal_plane.py b/diffsims/crystallography/reciprocal_lattice_point.py similarity index 92% rename from diffsims/crystallography/crystal_plane.py rename to diffsims/crystallography/reciprocal_lattice_point.py index 090f48a4..8e6459d0 100644 --- a/diffsims/crystallography/crystal_plane.py +++ b/diffsims/crystallography/reciprocal_lattice_point.py @@ -22,7 +22,7 @@ import numpy as np from orix.vector import Vector3d -from diffsims.diffraction.structure_factor import ( +from diffsims.structure_factor.structure_factor import ( get_kinematical_structure_factor, get_doyleturner_structure_factor, get_refraction_corrected_wavelength, @@ -32,10 +32,10 @@ _FLOAT_EPS = np.finfo(np.float).eps # Used to round values below 1e-16 to zero -class CrystalPlane: - """Crystal plane (or reciprocal lattice point, reflectors, g, etc.) +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. + other relevant structure_factor parameters. """ def __init__(self, phase, hkl): @@ -65,7 +65,10 @@ def __repr__(self): def __getitem__(self, key): new_cp = self.__class__(self.phase, self.hkl[key]) - new_cp._structure_factor = self.structure_factor[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 @@ -134,7 +137,7 @@ def structure_factor(self): @property def allowed(self): - """Return whether planes diffract according to diffraction + """Return whether planes diffract according to structure_factor selection rules assuming kinematical scattering theory. """ # Translational symmetry @@ -203,10 +206,6 @@ def from_highest_hkl(cls, phase, highest_hkl): hkl = get_hkl(highest_hkl=highest_hkl) return cls(phase=phase, hkl=hkl).unique() - @classmethod - def from_nfamilies(cls, phase, nfamilies=5): - raise NotImplementedError - def unique(self, use_symmetry=True): """Return planes with unique Miller indices. @@ -218,7 +217,7 @@ def unique(self, use_symmetry=True): Returns ------- - CrystalPlane + ReciprocalLatticePoint """ if use_symmetry: all_hkl = self._hkldata @@ -226,7 +225,7 @@ def unique(self, use_symmetry=True): families = defaultdict(list) for this_hkl in all_hkl.tolist(): for that_hkl in families.keys(): - if is_equivalent(this_hkl, that_hkl): + if _is_equivalent(this_hkl, that_hkl): families[tuple(that_hkl)].append(this_hkl) break else: @@ -238,10 +237,14 @@ def unique(self, use_symmetry=True): 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, + self, + antipodal=True, + unique=True, + return_multiplicity=False, ): """Return planes with symmetrically equivalent Miller indices. @@ -260,7 +263,7 @@ def symmetrise( Returns ------- - CrystalPlane + ReciprocalLatticePoint Planes with Miller indices symmetrically equivalent to the original planes. multiplicity : np.ndarray @@ -283,6 +286,7 @@ def symmetrise( 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] @@ -322,11 +326,16 @@ def calculate_structure_factor(self, method=None, voltage=None): 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, + phase=phase, + hkl=hkl, + scattering_parameter=s, ) else: structure_factors[i] = get_doyleturner_structure_factor( - phase=phase, hkl=hkl, scattering_parameter=s, voltage=voltage, + phase=phase, + hkl=hkl, + scattering_parameter=s, + voltage=voltage, ) self._structure_factor = np.where( structure_factors < _FLOAT_EPS, 0, structure_factors @@ -438,5 +447,5 @@ def get_equivalent_hkl(hkl, operations, unique=False, return_multiplicity=False) return new_hkl -def is_equivalent(this_hkl: list, that_hkl: list) -> bool: +def _is_equivalent(this_hkl: list, that_hkl: list) -> bool: return sorted(np.abs(this_hkl)) == sorted(np.abs(that_hkl)) diff --git a/diffsims/diffraction/__init__.py b/diffsims/structure_factor/__init__.py similarity index 88% rename from diffsims/diffraction/__init__.py rename to diffsims/structure_factor/__init__.py index 31df42bc..52befcb9 100644 --- a/diffsims/diffraction/__init__.py +++ b/diffsims/structure_factor/__init__.py @@ -18,15 +18,15 @@ """Calculation of scattering factors and structure factors.""" -from diffsims.diffraction.atomic_scattering_factor import ( +from diffsims.structure_factor.atomic_scattering_factor import ( get_doyleturner_atomic_scattering_factor, get_kinematical_atomic_scattering_factor, ) -from diffsims.diffraction.atomic_scattering_parameters import ( +from diffsims.structure_factor.atomic_scattering_parameters import ( get_atomic_scattering_parameters, get_element_id_from_string, ) -from diffsims.diffraction.structure_factor import ( +from diffsims.structure_factor.structure_factor import ( find_asymmetric_positions, get_doyleturner_structure_factor, get_kinematical_structure_factor, diff --git a/diffsims/diffraction/atomic_scattering_factor.py b/diffsims/structure_factor/atomic_scattering_factor.py similarity index 95% rename from diffsims/diffraction/atomic_scattering_factor.py rename to diffsims/structure_factor/atomic_scattering_factor.py index d4b56352..9b73cadf 100644 --- a/diffsims/diffraction/atomic_scattering_factor.py +++ b/diffsims/structure_factor/atomic_scattering_factor.py @@ -18,7 +18,7 @@ import numpy as np -from diffsims.diffraction.atomic_scattering_parameters import ( +from diffsims.structure_factor.atomic_scattering_parameters import ( get_atomic_scattering_parameters, get_element_id_from_string, ) @@ -30,6 +30,8 @@ def get_kinematical_atomic_scattering_factor(atom, scattering_parameter): Assumes structure's Debye-Waller factors are expressed in Ångströms. + This function is adapted from EMsoft. + Parameters ---------- atom : diffpy.structure.Atom @@ -69,6 +71,8 @@ def get_doyleturner_atomic_scattering_factor( Assumes structure's Debye-Waller factors are expressed in Ångströms. + This function is adapted from EMsoft. + Parameters ---------- atom : diffpy.structure.Atom diff --git a/diffsims/diffraction/atomic_scattering_parameters.py b/diffsims/structure_factor/atomic_scattering_parameters.py similarity index 100% rename from diffsims/diffraction/atomic_scattering_parameters.py rename to diffsims/structure_factor/atomic_scattering_parameters.py diff --git a/diffsims/diffraction/structure_factor.py b/diffsims/structure_factor/structure_factor.py similarity index 96% rename from diffsims/diffraction/structure_factor.py rename to diffsims/structure_factor/structure_factor.py index 54ab35d8..d15ce004 100644 --- a/diffsims/diffraction/structure_factor.py +++ b/diffsims/structure_factor/structure_factor.py @@ -20,7 +20,7 @@ import numpy as np from scipy.constants import c, e, h, physical_constants -from diffsims.diffraction.atomic_scattering_factor import ( +from diffsims.structure_factor.atomic_scattering_factor import ( get_kinematical_atomic_scattering_factor, get_doyleturner_atomic_scattering_factor, ) @@ -59,6 +59,8 @@ def get_kinematical_structure_factor(phase, hkl, scattering_parameter): Assumes structure's lattice parameters and Debye-Waller factors are expressed in Ångströms. + This function is adapted from EMsoft. + Parameters ---------- phase : orix.crystal_map.phase_list.Phase @@ -99,7 +101,11 @@ def get_kinematical_structure_factor(phase, hkl, scattering_parameter): def get_doyleturner_structure_factor( - phase, hkl, scattering_parameter, voltage, return_parameters=False, + 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]_. @@ -107,6 +113,8 @@ def get_doyleturner_structure_factor( Assumes structure's lattice parameters and Debye-Waller factors are expressed in Ångströms. + This function is adapted from EMsoft. + Parameters ---------- phase : orix.crystal_map.phase_list.Phase @@ -183,6 +191,8 @@ 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. + This function is adapted from EMsoft. + Parameters ---------- phase : orix.crystal_map.phase_list.Phase From def86c45a7c473c7a2939ea31cde059010dbe0e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20Wiik=20=C3=85nes?= Date: Tue, 8 Sep 2020 23:01:31 +0200 Subject: [PATCH 6/8] Test getting Doyle Turner scattering params, start structure factor MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Håkon Wiik Ånes --- .../atomic_scattering_parameters.py | 250 +++++++++--------- diffsims/tests/conftest.py | 54 +++- .../tests/test_crystallography/__init__.py | 17 ++ .../test_reciprocal_lattice_point.py | 17 ++ .../tests/test_structure_factor/__init__.py | 17 ++ .../test_atomic_scattering_factor.py | 96 +++++++ .../test_atomic_scattering_parameter.py | 170 ++++++++++++ .../test_structure_factor.py | 56 ++++ 8 files changed, 546 insertions(+), 131 deletions(-) create mode 100644 diffsims/tests/test_crystallography/__init__.py create mode 100644 diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py create mode 100644 diffsims/tests/test_structure_factor/__init__.py create mode 100644 diffsims/tests/test_structure_factor/test_atomic_scattering_factor.py create mode 100644 diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py create mode 100644 diffsims/tests/test_structure_factor/test_structure_factor.py diff --git a/diffsims/structure_factor/atomic_scattering_parameters.py b/diffsims/structure_factor/atomic_scattering_parameters.py index a27d6ef8..019b8a8c 100644 --- a/diffsims/structure_factor/atomic_scattering_parameters.py +++ b/diffsims/structure_factor/atomic_scattering_parameters.py @@ -18,6 +18,127 @@ import numpy as np +# List of elements Z = 1-118 +# 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) + +# fmt: off +ATOMIC_SCATTERING_PARAMETERS_DOYLETURNER = np.array([ + # a1 a2 a3 a4 b1 b2 b3 b4 + 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 + def get_atomic_scattering_parameters(element, unit=None): """Return the eight atomic scattering parameters a_1-4, b_1-4 for @@ -43,8 +164,8 @@ def get_atomic_scattering_parameters(element, unit=None): References ---------- - .. [DeGraef2007] M. De Graef, M. E. McHenry, "Structure of - \Materials," Cambridge University Press (2007). + .. [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. @@ -52,109 +173,6 @@ def get_atomic_scattering_parameters(element, unit=None): 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: @@ -165,8 +183,8 @@ def get_atomic_scattering_parameters(element, unit=None): if unit.lower() == "nm": factor = 1e-2 - a = all_parameters[element_id - 1, :4] * factor - b = all_parameters[element_id - 1, 4:] * factor + a = ATOMIC_SCATTERING_PARAMETERS_DOYLETURNER[element_id - 1, :4] * factor + b = ATOMIC_SCATTERING_PARAMETERS_DOYLETURNER[element_id - 1, 4:] * factor return a, b @@ -185,22 +203,6 @@ def get_element_id_from_string(element_str): 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))) + element2periodic = dict(zip(ELEMENTS[:N_ELEMENTS], np.arange(1, N_ELEMENTS))) element_id = element2periodic[element_str.lower()] return element_id diff --git a/diffsims/tests/conftest.py b/diffsims/tests/conftest.py index 57a43c63..de12095e 100644 --- a/diffsims/tests/conftest.py +++ b/diffsims/tests/conftest.py @@ -17,20 +17,18 @@ # along with diffsims. If not, see . import pytest -import diffpy.structure -import numpy as np -from transforms3d.euler import euler2mat +from diffpy.structure import Atom, Lattice, Structure +from orix.crystal_map import Phase -from diffsims.libraries.vector_library import DiffractionVectorLibrary from diffsims.generators.diffraction_generator import DiffractionGenerator @pytest.fixture def default_structure(): """An atomic structure represented using diffpy """ - latt = diffpy.structure.lattice.Lattice(3, 3, 5, 90, 90, 120) - atom = diffpy.structure.atom.Atom(atype="Ni", xyz=[0, 0, 0], lattice=latt) - hexagonal_structure = diffpy.structure.Structure(atoms=[atom], lattice=latt) + latt = Lattice(3, 3, 5, 90, 90, 120) + atom = Atom(atype="Ni", xyz=[0, 0, 0], lattice=latt) + hexagonal_structure = Structure(atoms=[atom], lattice=latt) return hexagonal_structure @@ -38,3 +36,45 @@ def default_structure(): def default_simulator(): accelerating_voltage = 300 return DiffractionGenerator(accelerating_voltage) + + +@pytest.fixture +def nickel_phase(): + return Phase( + space_group=225, + structure=Structure( + lattice=Lattice(3.5236, 3.5236, 3.5236, 90, 90, 90), + atoms=[Atom(xyz=[0, 0, 0], atype="Ni")] + ) + ) + + +@pytest.fixture +def ferrite_phase(): + return Phase( + space_group=229, + structure=Structure( + lattice=Lattice(2.8665, 2.8665, 2.8665, 90, 90, 90), + atoms=[ + Atom(xyz=[0, 0, 0], atype="Fe"), + Atom(xyz=[0.5, 0.5, 0.5], atype="Fe"), + ] + ) + ) + + +@pytest.fixture +def silicon_carbide_phase(): + """Silicon Carbide 4H polytype (hexagonal, space group 186).""" + return Phase( + space_group=186, + structure=Structure( + lattice=Lattice(3.073, 3.073, 10.053, 90, 90, 120), + atoms=[ + Atom(atype="Si", xyz=[0, 0, 0]), + Atom(atype="Si", xyz=[0.33, 0.667, 0.25]), + Atom(atype="C", xyz=[0, 0, 0.188]), + Atom(atype="C", xyz=[0.333, 0.667, 0.438]), + ], + ) + ) diff --git a/diffsims/tests/test_crystallography/__init__.py b/diffsims/tests/test_crystallography/__init__.py new file mode 100644 index 00000000..063cffcc --- /dev/null +++ b/diffsims/tests/test_crystallography/__init__.py @@ -0,0 +1,17 @@ +# -*- 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 . diff --git a/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py b/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py new file mode 100644 index 00000000..063cffcc --- /dev/null +++ b/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py @@ -0,0 +1,17 @@ +# -*- 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 . diff --git a/diffsims/tests/test_structure_factor/__init__.py b/diffsims/tests/test_structure_factor/__init__.py new file mode 100644 index 00000000..063cffcc --- /dev/null +++ b/diffsims/tests/test_structure_factor/__init__.py @@ -0,0 +1,17 @@ +# -*- 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 . diff --git a/diffsims/tests/test_structure_factor/test_atomic_scattering_factor.py b/diffsims/tests/test_structure_factor/test_atomic_scattering_factor.py new file mode 100644 index 00000000..61861c9b --- /dev/null +++ b/diffsims/tests/test_structure_factor/test_atomic_scattering_factor.py @@ -0,0 +1,96 @@ +# -*- 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 import Atom, Structure, Lattice +import numpy as np +import pytest + +from diffsims.structure_factor import ( + get_kinematical_atomic_scattering_factor, + get_doyleturner_atomic_scattering_factor, +) + + +@pytest.mark.parametrize( + "element, occupancy, displacement_factor, scattering_parameter, desired_factor", + [ + ("fe", 1, 0.5, 0.17442875, 6.306987), + ("fe", 1, 0.5, 0.90635835, 5.911219e-14), + ("ni", 1, 0.5, 0.14190033, 10.964046), + ("ni", 0.5, 0.5, 0.14190033, 5.482023), + ], +) +def test_get_kinematical_atomic_scattering_factor( + element, occupancy, displacement_factor, scattering_parameter, desired_factor +): + atom = Atom(atype=element, occupancy=occupancy, Uisoequiv=displacement_factor) + assert np.allclose( + get_kinematical_atomic_scattering_factor( + atom=atom, + scattering_parameter=scattering_parameter, + ), + desired_factor, + ) + + +nickel = Structure( + lattice=Lattice(3.5236, 3.5236, 3.5236, 90, 90, 90), + atoms=[Atom(xyz=[0, 0, 0], atype="Ni")], +) +ferrite = Structure( + lattice=Lattice(2.8665, 2.8665, 2.8665, 90, 90, 90), + atoms=[Atom(xyz=[0, 0, 0], atype="Fe")], +) +# Silicon Carbide 4H polytype (hexagonal, space group 186) +sic4h = Structure( + lattice=Lattice(3.073, 3.073, 10.053, 90, 90, 120), + atoms=[ + Atom(atype="Si", xyz=[0, 0, 0]), + Atom(atype="Si", xyz=[0.33, 0.667, 0.25]), + Atom(atype="C", xyz=[0, 0, 0.188]), + Atom(atype="C", xyz=[0.333, 0.667, 0.438]), + ], +) + + +@pytest.mark.parametrize( + "structure, displacement_factor, scattering_parameter, desired_factor", + [ + (ferrite, 0.5, 0.17442875, 2.422658), + (ferrite, 0.5, 0.90635835, 9.172283e-15), + (ferrite, 0, 0.90635835, 1.114444), + (nickel, 0.5, 0.14190033, 2.186904), + (sic4h, 0, 0.19894559, 1.512963), + ], +) +def test_get_doyleturner_atomic_scattering_factor( + structure, + displacement_factor, + scattering_parameter, + desired_factor, +): + atom = structure[0] + atom.Uisoequiv = displacement_factor + assert np.allclose( + get_doyleturner_atomic_scattering_factor( + atom=atom, + scattering_parameter=scattering_parameter, + unit_cell_volume=structure.lattice.volume, + ), + desired_factor, + ) diff --git a/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py b/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py new file mode 100644 index 00000000..3c0b4784 --- /dev/null +++ b/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py @@ -0,0 +1,170 @@ +# -*- 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 . + +import numpy as np +import pytest + +from diffsims.structure_factor import ( + get_atomic_scattering_parameters, + get_element_id_from_string, +) + + +@pytest.mark.parametrize( + "element, unit, desired_parameters", + [ + (10, "m", [0.303, 0.720, 0.475, 0.153, 17.640, 5.860, 1.762, 0.266]), + ("Fe", None, [2.544, 2.343, 1.759, 0.506, 64.424, 14.880, 2.854, 0.350]), + ( + "Fe", + "NM", + [0.02544, 0.02343, 0.01759, 0.00506, 0.64424, 0.14880, 0.02854, 0.00350], + ), + ("bk", "Å", [6.502, 5.478, 2.510, 0.000, 28.375, 4.975, 0.561, 0.0]), + ], +) +def test_get_atomic_scattering_parameters(element, unit, desired_parameters): + a, b = get_atomic_scattering_parameters(element, unit=unit) + assert np.allclose(a, desired_parameters[:4]) + assert np.allclose(b, desired_parameters[4:]) + + +def test_get_element_id_from_string(): + 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", + ] + elements_id = np.arange(98) + 1 + for i, element in zip(elements_id, elements): + assert get_element_id_from_string(element) == i diff --git a/diffsims/tests/test_structure_factor/test_structure_factor.py b/diffsims/tests/test_structure_factor/test_structure_factor.py new file mode 100644 index 00000000..e83071b9 --- /dev/null +++ b/diffsims/tests/test_structure_factor/test_structure_factor.py @@ -0,0 +1,56 @@ +# -*- 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.spacegroups import GetSpaceGroup +import numpy as np +import pytest + +from diffsims.structure_factor import ( + find_asymmetric_positions, + get_kinematical_structure_factor, + get_doyleturner_structure_factor, + get_refraction_corrected_wavelength, +) + + +@pytest.mark.parametrize( + "positions, space_group, desired_asymmetric", + [ + ([[0, 0, 0], [0, 0.5, 0.5]], 229, [True, False]), + ([[0.0349, 0.3106, 0.2607], [0.0349, 0.3106, 0.2607]], 229, [True, True]), + ], +) +def test_find_asymmetric_positions(positions, space_group, desired_asymmetric): + assert np.allclose( + find_asymmetric_positions( + positions=positions, space_group=GetSpaceGroup(space_group) + ), + desired_asymmetric, + ) + + +def test_get_kinematical_structure_factor(): + pass + + +def test_get_doyleturner_structure_factor(): + pass + + +def test_get_refraction_corrected_wavelength(): + pass From 2f0383c905609394b796b9be505e0b9130015ff2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20Wiik=20=C3=85nes?= Date: Tue, 8 Sep 2020 23:05:51 +0200 Subject: [PATCH 7/8] Restrict black formatting MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Håkon Wiik Ånes --- .../test_atomic_scattering_parameter.py | 134 ++---------------- 1 file changed, 14 insertions(+), 120 deletions(-) diff --git a/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py b/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py index 3c0b4784..c25f7d8b 100644 --- a/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py +++ b/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py @@ -33,10 +33,10 @@ ( "Fe", "NM", - [0.02544, 0.02343, 0.01759, 0.00506, 0.64424, 0.14880, 0.02854, 0.00350], + [0.02544, 0.02343, 0.01759, 0.00506, 0.64424, 0.14880, 0.02854, 0.00350] ), ("bk", "Å", [6.502, 5.478, 2.510, 0.000, 28.375, 4.975, 0.561, 0.0]), - ], + ] ) def test_get_atomic_scattering_parameters(element, unit, desired_parameters): a, b = get_atomic_scattering_parameters(element, unit=unit) @@ -45,126 +45,20 @@ def test_get_atomic_scattering_parameters(element, unit, desired_parameters): def test_get_element_id_from_string(): + # 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", + "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 elements_id = np.arange(98) + 1 for i, element in zip(elements_id, elements): assert get_element_id_from_string(element) == i From d9a6dafb55d451ea2940feb1efe6bb8184b3a360 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20Wiik=20=C3=85nes?= Date: Wed, 9 Sep 2020 22:09:34 +0200 Subject: [PATCH 8/8] Structure factor tested against EMsoft, almost finished rlp tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Håkon Wiik Ånes --- diffsims/structure_factor/structure_factor.py | 16 +- diffsims/tests/conftest.py | 8 +- .../test_reciprocal_lattice_point.py | 142 ++++++++++++++++++ .../test_atomic_scattering_parameter.py | 4 +- .../test_structure_factor.py | 117 ++++++++++++++- 5 files changed, 268 insertions(+), 19 deletions(-) diff --git a/diffsims/structure_factor/structure_factor.py b/diffsims/structure_factor/structure_factor.py index d15ce004..946c5098 100644 --- a/diffsims/structure_factor/structure_factor.py +++ b/diffsims/structure_factor/structure_factor.py @@ -27,6 +27,7 @@ rest_mass = physical_constants["atomic unit of mass"][0] +# rest_mass = 9.109383709e-31 # Used by EMsoft def find_asymmetric_positions(positions, space_group): @@ -66,7 +67,7 @@ def get_kinematical_structure_factor(phase, hkl, scattering_parameter): 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 + hkl : np.ndarray or list Miller indices. scattering_parameter : float Scattering parameter for these Miller indices. @@ -120,7 +121,7 @@ def get_doyleturner_structure_factor( 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 + hkl : np.ndarray or list Miller indices. scattering_parameter : float Scattering parameter for these Miller indices. @@ -164,15 +165,13 @@ def get_doyleturner_structure_factor( 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 + # Relativistic correction factor of wavelength gamma_relcor = 1 + (2 * e * 0.5 * voltage / rest_mass / (c ** 2)) + # Mean inner potential 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: @@ -210,7 +209,8 @@ def get_refraction_corrected_wavelength(phase, voltage): temp1 = 1e9 * h / np.sqrt(2 * rest_mass * e) temp2 = e * 0.5 * voltage / rest_mass / (c ** 2) - # Relativistic correction factor (known as gamma) + # Relativistic correction factor (known as gamma). (This is used by EMsoft but not + # here, for now.) # gamma_relcor = 1 + (2 * temp2) # Relativistic acceleration voltage @@ -226,7 +226,7 @@ def get_refraction_corrected_wavelength(phase, voltage): psi_hat += v_mod wavelength = temp1 / np.sqrt(psi_hat) - # Interaction constant sigma + # Interaction constant sigma (this is used by EMsoft but not here, for now) # sigma = 1e-18 * (2 * np.pi * rest_mass * gamma_relcor * e * wavelength) / h ** 2 return wavelength diff --git a/diffsims/tests/conftest.py b/diffsims/tests/conftest.py index de12095e..70168058 100644 --- a/diffsims/tests/conftest.py +++ b/diffsims/tests/conftest.py @@ -41,10 +41,11 @@ def default_simulator(): @pytest.fixture def nickel_phase(): return Phase( + name="nickel", space_group=225, structure=Structure( lattice=Lattice(3.5236, 3.5236, 3.5236, 90, 90, 90), - atoms=[Atom(xyz=[0, 0, 0], atype="Ni")] + atoms=[Atom(xyz=[0, 0, 0], atype="Ni", Uisoequiv=0.006332)] ) ) @@ -52,12 +53,13 @@ def nickel_phase(): @pytest.fixture def ferrite_phase(): return Phase( + name="ferrite", space_group=229, structure=Structure( lattice=Lattice(2.8665, 2.8665, 2.8665, 90, 90, 90), atoms=[ - Atom(xyz=[0, 0, 0], atype="Fe"), - Atom(xyz=[0.5, 0.5, 0.5], atype="Fe"), + Atom(xyz=[0, 0, 0], atype="Fe", Uisoequiv=0.006332), + Atom(xyz=[0.5, 0.5, 0.5], atype="Fe", Uisoequiv=0.006332), ] ) ) diff --git a/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py b/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py index 063cffcc..b095aa91 100644 --- a/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py +++ b/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py @@ -15,3 +15,145 @@ # # You should have received a copy of the GNU General Public License # along with diffsims. If not, see . + +from diffsims.crystallography import ReciprocalLatticePoint +import numpy as np +from orix.vector import Vector3d +import pytest + + +class TestReciprocalLatticePoint: + @pytest.mark.parametrize( + "hkl", [[[1, 1, 1], [2, 0, 0]], np.array([[1, 1, 1], [2, 0, 0]])] + ) + def test_init_rlp(self, nickel_phase, hkl): + rlp = ReciprocalLatticePoint(phase=nickel_phase, hkl=hkl) + assert rlp.phase.name == nickel_phase.name + assert isinstance(rlp.hkl, Vector3d) + assert rlp.structure_factor[0] is None + assert rlp.theta[0] is None + assert rlp.size == 2 + assert rlp.shape == (2, 3) + assert rlp.hkl[0].shape == (1,) + assert rlp._hkldata[0].shape == (3,) + assert np.issubdtype(rlp.hkl.data.dtype, int) + + @pytest.mark.parametrize("min_dspacing, desired_size", [(2, 9), (1, 19), (0.5, 83)]) + def test_init_from_min_dspacing(self, ferrite_phase, min_dspacing, desired_size): + assert ReciprocalLatticePoint.from_min_dspacing( + phase=ferrite_phase, min_dspacing=min_dspacing + ).size == desired_size + + @pytest.mark.parametrize( + "highest_hkl, desired_highest_hkl, desired_lowest_hkl, desired_size", + [ + ([3, 3, 3], [3, 3, 3], [1, 0, 0], 19), + ([3, 4, 0], [3, 4, 0], [0, 4, 0], 13), + ([4, 3, 0], [4, 3, 0], [1, 0, 0], 13) + ] + ) + def test_init_from_highest_hkl( + self, + silicon_carbide_phase, + highest_hkl, + desired_highest_hkl, + desired_lowest_hkl, + desired_size + ): + rlp = ReciprocalLatticePoint.from_highest_hkl( + phase=silicon_carbide_phase, highest_hkl=highest_hkl + ) + assert np.allclose(rlp[0]._hkldata, desired_highest_hkl) + assert np.allclose(rlp[-1]._hkldata, desired_lowest_hkl) + assert rlp.size == desired_size + + def test_repr(self, ferrite_phase): + rlp = ReciprocalLatticePoint.from_min_dspacing(ferrite_phase, min_dspacing=2) + assert repr(rlp) == ( + f"ReciprocalLatticePoint (9,)\n" + f"Phase: ferrite (m-3m)\n" + "[[2 2 2]\n [2 2 1]\n [2 2 0]\n [2 1 1]\n [2 1 0]\n [2 0 0]\n [1 1 1]\n" + " [1 1 0]\n [1 0 0]]" + ) + + def test_get_item(self): + pass + + def test_get_hkl(self, silicon_carbide_phase): + rlp = ReciprocalLatticePoint.from_min_dspacing( + silicon_carbide_phase, min_dspacing=3 + ) + assert np.allclose(rlp.h, [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]) + assert np.allclose(rlp.k, [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0]) + assert np.allclose(rlp.l, [4, 3, 2, 1, 0, 4, 3, 2, 0, 4, 3, 2]) + + def test_multiplicity(self, nickel_phase, silicon_carbide_phase): + assert np.allclose( + ReciprocalLatticePoint.from_min_dspacing( + phase=nickel_phase, min_dspacing=1 + ).multiplicity, + # fmt: off + np.array([ + 8, 24, 24, 24, 12, 24, 48, 48, 24, 24, 48, 24, 24, 24, 6, 8, 24, + 24, 12, 24, 48, 24, 24, 24, 6, 8, 24, 12, 24, 24, 6, 8, 12, 6 + ]) + # fmt: on + ) + assert np.allclose( + ReciprocalLatticePoint.from_min_dspacing( + phase=silicon_carbide_phase, min_dspacing=1 + ).multiplicity, + # fmt: off + np.array([ + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, + 6, 6, 6, 6, 6, 6, 6, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 6, 6, 6, 6, 6, 6, 6, 6, 12, + 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 6, 6, 6, 6, 6, 6, 6, 6, 1, 1, + 1, 1, 1, 1, 1, 1 + ]) + # fmt: on + ) + + def test_gspacing_dspacing_scattering_parameter(self, ferrite_phase): + rlp = ReciprocalLatticePoint.from_min_dspacing( + phase=ferrite_phase, min_dspacing=2 + ) + # fmt: off + assert np.allclose( + rlp.gspacing, + np.array([ + 1.2084778, 1.04657248, 0.98671799, 0.85452285, 0.78006907, 0.69771498, + 0.604238, 0.493359, 0.34885749 + ]) + ) + assert np.allclose( + rlp.dspacing, + np.array([ + 0.82748727, 0.9555, 1.01346079, 1.17024372, 1.28193777, 1.43325, + 1.65497455, 2.02692159, 2.8665 + ]) + ) + assert np.allclose( + rlp.scattering_parameter, + np.array([ + 0.6042389, 0.52328624, 0.493359, 0.42726142, 0.39003453, 0.34885749, + 0.30211945, 0.2466795, 0.17442875 + ]) + ) + # fmt: on + + def test_allowed(self): + pass + + def test_unique(self): + pass + + def test_symmetrise(self): + pass + + def test_calculate_structure_factor(self): + pass + + def test_calculate_theta(self): + pass diff --git a/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py b/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py index c25f7d8b..b425a4b1 100644 --- a/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py +++ b/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py @@ -33,10 +33,10 @@ ( "Fe", "NM", - [0.02544, 0.02343, 0.01759, 0.00506, 0.64424, 0.14880, 0.02854, 0.00350] + [0.02544, 0.02343, 0.01759, 0.00506, 0.64424, 0.14880, 0.02854, 0.00350], ), ("bk", "Å", [6.502, 5.478, 2.510, 0.000, 28.375, 4.975, 0.561, 0.0]), - ] + ], ) def test_get_atomic_scattering_parameters(element, unit, desired_parameters): a, b = get_atomic_scattering_parameters(element, unit=unit) diff --git a/diffsims/tests/test_structure_factor/test_structure_factor.py b/diffsims/tests/test_structure_factor/test_structure_factor.py index e83071b9..618b32a3 100644 --- a/diffsims/tests/test_structure_factor/test_structure_factor.py +++ b/diffsims/tests/test_structure_factor/test_structure_factor.py @@ -17,7 +17,9 @@ # along with diffsims. If not, see . from diffpy.structure.spacegroups import GetSpaceGroup +from diffpy.structure import Atom, Lattice, Structure import numpy as np +from orix.crystal_map import Phase import pytest from diffsims.structure_factor import ( @@ -27,6 +29,37 @@ get_refraction_corrected_wavelength, ) +# Debye-Waller factor in Å^-2: Biosequiv = 8 * np.pi ** 2 * Uisoequiv +nickel = Phase( + space_group=225, + structure=Structure( + lattice=Lattice(3.5236, 3.5236, 3.5236, 90, 90, 90), + atoms=[Atom(xyz=[0, 0, 0], atype="Ni", Uisoequiv=0.006332)], + ), +) +ferrite = Phase( + space_group=229, + structure=Structure( + lattice=Lattice(2.8665, 2.8665, 2.8665, 90, 90, 90), + atoms=[ + Atom(xyz=[0, 0, 0], atype="Fe", Uisoequiv=0.006332), + Atom(xyz=[0.5, 0.5, 0.5], atype="Fe", Uisoequiv=0.006332), + ], + ), +) +sic4h = Phase( + space_group=186, + structure=Structure( + lattice=Lattice(3.073, 3.073, 10.053, 90, 90, 120), + atoms=[ + Atom(atype="Si", xyz=[0, 0, 0], Uisoequiv=0.006332), + Atom(atype="Si", xyz=[0.33, 0.667, 0.25], Uisoequiv=0.006332), + Atom(atype="C", xyz=[0, 0, 0.188], Uisoequiv=0.006332), + Atom(atype="C", xyz=[0.333, 0.667, 0.438], Uisoequiv=0.006332), + ], + ), +) + @pytest.mark.parametrize( "positions, space_group, desired_asymmetric", @@ -44,13 +77,85 @@ def test_find_asymmetric_positions(positions, space_group, desired_asymmetric): ) -def test_get_kinematical_structure_factor(): - pass +@pytest.mark.parametrize( + "phase, hkl, scattering_parameter, desired_factor", + [ + (nickel, [1, 1, 1], 0.245778, 79.665427), + (ferrite, [1, 1, 0], 0.2466795, 35.783295), + (ferrite, [1, 1, 1], 0.2466795, 0), + (sic4h, [0, 0, 4], 0.198945, 19.027004), + ], +) +def test_get_kinematical_structure_factor( + phase, hkl, scattering_parameter, desired_factor +): + assert np.allclose( + get_kinematical_structure_factor( + phase=phase, hkl=hkl, scattering_parameter=scattering_parameter + ), + desired_factor, + ) -def test_get_doyleturner_structure_factor(): - pass +@pytest.mark.parametrize( + "phase, hkl, scattering_parameter, voltage, desired_factor", + [ + (nickel, [1, 1, 1], 0.245778, 15e3, 8.606363), + (ferrite, [1, 1, 0], 0.2466795, 20e3, 8.096651), + (ferrite, [1, 1, 1], 0.2466795, 20e3, 0), + (sic4h, [0, 0, 4], 0.198945, 200e3, 2.744304), + ], +) +def test_get_doyleturner_structure_factor( + phase, hkl, scattering_parameter, voltage, desired_factor +): + """Tested against EMsoft v5.""" + assert np.allclose( + get_doyleturner_structure_factor( + phase=phase, + hkl=hkl, + scattering_parameter=scattering_parameter, + voltage=voltage, + ), + desired_factor, + ) -def test_get_refraction_corrected_wavelength(): - pass +@pytest.mark.parametrize( + "hkl, scattering_parameter, desired_factor, desired_params", + [ + ([1, 1, 0], 0.2466795, 8.096628, [1.03, 12.17, 1.22e-16, 12.17 + 2.45e-15j]), + ([2, 0, 0], 0.348857, 5.581301, [1.03, 8.39, 1.22e-16, 8.39 + 1.02e-15j]), + ], +) +def test_get_doyleturner_structure_factor_returns( + ferrite_phase, hkl, scattering_parameter, desired_factor, desired_params +): + """Tested against EMsoft v5. Only the imaginary part of v_g is off.""" + sf, params = get_doyleturner_structure_factor( + phase=ferrite_phase, + hkl=hkl, + scattering_parameter=scattering_parameter, + voltage=20e3, + return_parameters=True, + ) + assert np.allclose(sf, desired_factor) + assert np.allclose(list(params.values()), desired_params, atol=1e-2) + + +@pytest.mark.parametrize( + "phase, voltage, desired_wavelength", + [ + (nickel, 20e3, 0.008582), + (nickel, 200e3, 0.002507), + (ferrite, 10e3, 0.012186), + (sic4h, 100e3, 0.003701), + ], +) +def test_get_refraction_corrected_wavelength(phase, voltage, desired_wavelength): + """Tested against EMsoft v5.""" + assert np.allclose( + get_refraction_corrected_wavelength(phase=phase, voltage=voltage), + desired_wavelength, + atol=1e-6, + )