Use Phase.expand_asymmetric_unit

Author: viljarjf

diffsims/crystallography/_diffracting_vector.py

@@ -37,7 +37,10 @@ class RotatedPhase(Phase):
     """
 
     def __init__(self, phase: Phase, rotation: Rotation):
+        # Copy structure to not override the input phase.
+        # Use private members to avoid re-computing
         self._structure = Structure(phase.structure)
+        self._diffpy_lattice = phase._diffpy_lattice
         self.name = phase.name
         self.space_group = phase.space_group
         self.point_group = phase.point_group

diffsims/crystallography/reciprocal_lattice_vector.py

@@ -19,8 +19,6 @@
 from collections import defaultdict
 from copy import deepcopy
 
-from diffpy.structure.symmetryutilities import expandPosition
-from diffpy.structure import Structure, Lattice
 import numba as nb
 import numpy as np
 from orix.vector import Miller, Vector3d
@@ -919,22 +917,12 @@ def sanitise_phase(self):
 
         self._raise_if_no_space_group()
 
-        space_group = self.phase.space_group
-        structure = self.phase.structure
+        self.phase.expand_asymmetric_unit()
 
-        # Lattice basis transform for hexagonal crystals can sometimes move atoms below 0
-        for atom in structure:
-            for i in range(3):
-                if -1e-6 < atom.xyz[i] < 0:
-                    atom.xyz[i] = 0
-
-        new_structure = _expand_unit_cell(space_group, structure)
-        for atom in new_structure:
+        for atom in self.phase.structure:
             if np.issubdtype(type(atom.element), np.integer):
                 atom.element = _get_string_from_element_id(atom.element)
 
-        self.phase.structure = new_structure
-
     def symmetrise(self, return_multiplicity=False, return_index=False):
         """Unique vectors symmetrically equivalent to the vectors.
 
@@ -1528,7 +1516,9 @@ def unique(self, use_symmetry=False, return_index=False, return_inverse=False):
 
         # TODO: Reduce floating point precision in orix!
         def miller_unique(miller, use_symmetry=False):
-            v, idx, inv = Vector3d(miller).unique(return_index=True, return_inverse=True)
+            v, idx, inv = Vector3d(miller).unique(
+                return_index=True, return_inverse=True, ignore_zero=False
+            )
 
             if use_symmetry:
                 operations = miller.phase.point_group
@@ -1540,8 +1530,10 @@ def miller_unique(miller, use_symmetry=False):
                 for i in range(n_v):
                     a = data[i]
                     order = np.lexsort(a.T)  # Sort by column 1, 2, then 3
-                    data_sorted[i] = a[order[::-1]] # Reverse to preserve order
-                _, idx, inv = np.unique(data_sorted, return_index=True, return_inverse=True, axis=0)
+                    data_sorted[i] = a[order[::-1]]  # Reverse to preserve order
+                _, idx, inv = np.unique(
+                    data_sorted, return_index=True, return_inverse=True, axis=0
+                )
                 v = v[idx]
 
             m = miller.__class__(xyz=v.data, phase=miller.phase)
@@ -1598,74 +1590,6 @@ def stack(cls, sequence):
         return new
 
 
-# TODO: Upstream to diffpy.structure.Atom.__eq__()
-def _atom_eq(atom1, atom2):
-    """Determine whether two atoms are equal.
-
-    Parameters
-    ----------
-    atom1, atom2 : diffpy.structure.Atom
-
-    Returns
-    -------
-    bool
-
-    """
-
-    return (
-        atom1.element == atom2.element
-        and np.allclose(atom1.xyz_cartn, atom2.xyz_cartn, atol=1e-7)
-        and atom1.occupancy == atom2.occupancy
-        and np.allclose(atom1.U, atom2.U, atol=1e-7)
-        and np.allclose(atom1.Uisoequiv, atom2.Uisoequiv, atol=1e-7)
-    )
-
-
-# TODO: Upstream to orix.crystal_map.Phase.expand_structure()
-def _expand_unit_cell(space_group, structure):
-    """Expand a unit cell with symmetrically equivalent atoms.
-
-    Parameters
-    ----------
-    space_group : diffpy.structure.spacegroupmod.SpaceGroup
-        Space group describing the symmetry operations of the unit cell.
-    structure : diffpy.structure.Structure
-        Initial structure with atoms.
-
-    Returns
-    -------
-    new_structure : diffpy.structure.Structure
-
-    """
-    # Transform to diffpy axis conventions
-    structure_matrix = structure.lattice.base
-    pos = structure.xyz_cartn
-    structure = Structure(
-        atoms=list(structure),
-        lattice=Lattice(*structure.lattice.abcABG()),
-    )
-    new_structure = Structure(lattice=structure.lattice)
-    structure.xyz_cartn = pos
-
-    # Perform symmetry expansion
-    for atom in structure:
-        equal = []
-        for atom2 in new_structure:
-            equal.append(_atom_eq(atom, atom2))
-        if not any(equal):
-            new_positions = expandPosition(space_group, atom.xyz)[0]
-            for new_position in new_positions:
-                new_atom = deepcopy(atom)
-                new_atom.xyz = new_position
-                new_structure.append(new_atom)
-
-    # Transform back to orix convention
-    old_xyz_cartn = new_structure.xyz_cartn
-    new_structure.lattice.setLatBase(structure_matrix)
-    new_structure.xyz_cartn = old_xyz_cartn
-    return new_structure
-
-
 @nb.njit(
     "int64[:](float64[:, :], float64[:, :], int64[:])",
     cache=True,

diffsims/tests/generators/test_simulation_generator.py

@@ -197,7 +197,9 @@ def test_appropriate_scaling(self, diffraction_calculator: SimulationGenerator):
             phase=big_silicon, reciprocal_radius=5.0
         )
         indices = [tuple(i.tolist()) for i in diffraction.coordinates.hkl.astype(int)]
-        big_indices = [tuple(i.tolist()) for i in big_diffraction.coordinates.hkl.astype(int)]
+        big_indices = [
+            tuple(i.tolist()) for i in big_diffraction.coordinates.hkl.astype(int)
+        ]
         target = (2, 2, 0)
         assert target in indices
         assert target in big_indices
@@ -413,3 +415,41 @@ def test_calculate_diffraction2d_progressbar_multi_phase(capsys):
     expected2 = "B: 100\%\|█+\| 10\/10 "
     assert re.findall(expected1, captured.err)
     assert re.findall(expected2, captured.err)
+def test_space_group_is_accounted_for():
+    structure = diffpy.structure.Structure(
+        atoms=[diffpy.structure.Atom("Au", xyz=(0, 0, 0))],
+        lattice=diffpy.structure.Lattice(4, 4, 4, 90, 90, 90),
+    )
+    P1 = Phase(space_group=1, structure=structure)
+    Fd3m = Phase(space_group=227, structure=structure)
+
+    gen = SimulationGenerator()
+    P1_sim = gen.calculate_diffraction2d(P1)
+    Fd3m_sim = gen.calculate_diffraction2d(Fd3m)
+
+    # Use ReciprocalLatticeVector.sanitize_phase
+    Fd3m_expanded = Fd3m.deepcopy()
+    Fd3m_expanded.expand_asymmetric_unit()
+    Fd3m_expanded_sim = gen.calculate_diffraction2d(Fd3m_expanded)
+
+    # Check coordinates. There are many extinct reflections in Fd3m, so look for those in P1
+    P1_hkl = [tuple(hkl.round().astype(int).tolist()) for hkl in P1_sim.coordinates.hkl]
+    Fd3m_hkl = [
+        tuple(hkl.round().astype(int).tolist()) for hkl in Fd3m_sim.coordinates.hkl
+    ]
+    Fd3m_expanded_hkl = [
+        tuple(hkl.round().astype(int).tolist())
+        for hkl in Fd3m_expanded_sim.coordinates.hkl
+    ]
+    assert set(P1_hkl).issuperset(set(Fd3m_hkl))
+    assert set(P1_hkl).issuperset(set(Fd3m_expanded_hkl))
+    assert set(Fd3m_hkl) == set(Fd3m_expanded_hkl)
+
+    # Check intensities, should be different
+    order = [P1_hkl.index(hkl) for hkl in Fd3m_hkl]
+    assert not np.allclose(
+        P1_sim.coordinates.intensity[order], Fd3m_sim.coordinates.intensity
+    )
+    assert np.allclose(
+        Fd3m_sim.coordinates.intensity, Fd3m_expanded_sim.coordinates.intensity
+    )