Reuse first termination slabs in CoherentInterfaceBuilder._find_matches

src/pymatgen/analysis/interfaces/coherent_interfaces.py

@@ -11,14 +11,35 @@
 from pymatgen.analysis.elasticity.strain import Deformation
 from pymatgen.analysis.interfaces.zsl import ZSLGenerator, fast_norm
 from pymatgen.core.interface import Interface, label_termination
-from pymatgen.core.surface import SlabGenerator
+from pymatgen.core.surface import Slab, SlabGenerator
 
 if TYPE_CHECKING:
     from collections.abc import Iterator, Sequence
 
     from pymatgen.core import Structure
 
 
+def _check_equal_slab_lattices(slabs: Sequence[Slab], miller: tuple[int, int, int], structure: Structure) -> None:
+    """Verifies that the lattices of a group of slabs are equivalent.
+
+    Args:
+        slabs (Sequence[Slab]): Sequence of Slabs to check. Intended: From the same SlabGenerator.
+        miller (tuple[int, int, int]): Miller index of the slabs (for error messages).
+        structure (Structure): Source Structure of the slabs (for error messages).
+
+    Raises:
+        ValueError: If no slabs are present or the slabs have diverging lattices.
+    """
+    if not slabs:
+        raise ValueError(f"Could not generate Slabs for miller index {miller} of {structure.composition}.")
+    for slab in slabs[1:]:
+        if not np.allclose(slab.lattice.matrix[:2], slabs[0].lattice.matrix[:2]):
+            raise ValueError(
+                f"Not all slabs for miller index {miller} of {structure.composition} have the same lattice. "
+                f"First slab: {slabs[0].lattice.params_dict}, found: {slab.lattice.params_dict}."
+            )
+
+
 class CoherentInterfaceBuilder:
     """
     This class constructs the coherent interfaces between two crystalline slabs
@@ -68,37 +89,17 @@ def __init__(
         self.termination_ftol = termination_ftol
         self.label_index = label_index
         self.filter_out_sym_slabs = filter_out_sym_slabs
-        self._find_matches()
-        self._find_terminations()
-
-    def _find_matches(self) -> None:
-        """Find and stores the ZSL matches."""
-        self.zsl_matches = []
+        film_slab, substrate_slab = self._find_terminations()
+        self._find_matches(film_slab, substrate_slab)
 
-        film_sg = SlabGenerator(
-            self.film_structure,
-            self.film_miller,
-            min_slab_size=1,
-            min_vacuum_size=3,
-            in_unit_planes=True,
-            center_slab=True,
-            primitive=True,
-            reorient_lattice=False,  # This is necessary to not screw up the lattice
-        )
+    def _find_matches(self, film_slab: Slab, sub_slab: Slab) -> None:
+        """Find and stores the ZSL matches.
 
-        sub_sg = SlabGenerator(
-            self.substrate_structure,
-            self.substrate_miller,
-            min_slab_size=1,
-            min_vacuum_size=3,
-            in_unit_planes=True,
-            center_slab=True,
-            primitive=True,
-            reorient_lattice=False,  # This is necessary to not screw up the lattice
-        )
-
-        film_slab = film_sg.get_slab(shift=0)
-        sub_slab = sub_sg.get_slab(shift=0)
+        Args:
+            film_slab (Slab): A slab of the film structure, transformed as used in the other functions.
+            sub_slab (Slab): Associated slab of the substrate structure.
+        """
+        self.zsl_matches = []
 
         film_vectors = film_slab.lattice.matrix
         substrate_vectors = sub_slab.lattice.matrix
@@ -119,8 +120,12 @@ def _find_matches(self) -> None:
                     "Substrate lattice vectors changed during ZSL match, check your ZSL Generator parameters"
                 )
 
-    def _find_terminations(self):
-        """Find all terminations."""
+    def _find_terminations(self) -> tuple[Slab, Slab]:
+        """Find all terminations.
+
+        Returns:
+            tuple[Slab, Slab]: The slab pair of the first termination, useful to check the primitivization.
+        """
         film_sg = SlabGenerator(
             self.film_structure,
             self.film_miller,
@@ -180,6 +185,12 @@ def _find_terminations(self):
         }
         self.terminations = list(self._terminations)
 
+        # Verify slabs: All slabs should have the same lattice (primitivization should be shift-independent)
+        # (If they don't, the assumption of _find_matches [all lattices are slab 0's lattice] is false)
+        _check_equal_slab_lattices(film_slabs, self.film_miller, self.film_structure)
+        _check_equal_slab_lattices(sub_slabs, self.substrate_miller, self.substrate_structure)
+        return film_slabs[0], sub_slabs[0]
+
     def get_interfaces(
         self,
         termination: tuple[str, str],

tests/analysis/interfaces/test_coherent_interface.py

@@ -89,3 +89,19 @@ def test_termination_searching(self):
             ("1_Ge_C2/m_2", "1_Si_P4/mmm_1"),
             ("1_Ge_C2/m_2", "2_Si_P4/mmm_1"),
         ], "cib2 termination results wrong"
+
+    def test_shiftdependent_primitive(self):
+        """Test that a crystal structure which primitivizes correctly at nonzero shift works."""
+        # This is technically unwanted behaviour (ideally shift 0 would also work)
+        # thus we do not check that shift 0 fails, instead only nonzero shift is tested
+        # Approximately calcite (spacegroup R-3c)
+        struct = Structure.from_spacegroup(
+            167, Lattice.hexagonal(5, 17), ("Ca", "C", "O"), ((0, 0, 0), (0, 0, 0.25), (0.25, 0, 0.25))
+        )
+        cib = CoherentInterfaceBuilder(struct, struct, (0, 0, 1), (0, 0, 1))
+        # The names are allowed to change, but if they do, verify that the slabs are the same as before
+        assert cib.terminations[0] == ("CaCO_Pmma_6", "CaCO_Pmma_6")
+        # The lattice should be small (5 and ~8.09 instead of 5 and ~24.27)
+        inter = next(cib.get_interfaces(cib.terminations[0]))
+        assert_allclose(inter.lattice.matrix[0, 0], 5, atol=1e-4)
+        assert_allclose(inter.lattice.matrix[1, 1], 8.089774, atol=1e-4)