Structure factor tested against EMsoft, almost finished rlp tests

Signed-off-by: Håkon Wiik Ånes <hwaanes@gmail.com>

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

diffsims/tests/conftest.py

@@ -41,23 +41,25 @@ 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)]
         )
     )
 
 
 @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),
             ]
         )
     )

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 <http://www.gnu.org/licenses/>.
+
+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

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)

diffsims/tests/test_structure_factor/test_structure_factor.py

@@ -17,7 +17,9 @@
 # along with diffsims.  If not, see <http://www.gnu.org/licenses/>.
 
 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,
+    )