Fix lll_reduce for slab generation

src/pymatgen/analysis/interfaces/coherent_interfaces.py

@@ -34,22 +34,32 @@ def __init__(
         film_miller: Tuple3Ints,
         substrate_miller: Tuple3Ints,
         zslgen: ZSLGenerator | None = None,
+        termination_ftol: float = 0.25,
+        label_index: bool = False,  # necessary to add index to termination
+        filter_out_sym_slabs: bool = True,
     ):
         """
         Args:
-            substrate_structure: structure of substrate
-            film_structure: structure of film
-            film_miller: miller index of the film layer
-            substrate_miller: miller index for the substrate layer
-            zslgen: BiDirectionalZSL if you want custom lattice matching tolerances for coherency.
+            substrate_structure (Structure): substrate structure
+            film_structure (Structure): film structure
+            film_miller (tuple[int, int, int]): miller index for the film layer
+            substrate_miller (tuple[int, int, int]): miller index for the substrate layer
+            zslgen (ZSLGenerator | None): BiDirectionalZSL if you want custom lattice matching tolerances for coherency.
+            termination_ftol (float): tolerance to distinguish different terminating atomic planes.
+            label_index (bool): If True add an extra index at the beginning of the termination label.
+            filter_out_sym_slabs (bool): If True filter out identical slabs with different terminations.
+                This might need to be set as False to find more non-identical terminations because slab
+                identity separately does not mean combinational identity.
         """
         # Bulk structures
         self.substrate_structure = substrate_structure
         self.film_structure = film_structure
         self.film_miller = film_miller
         self.substrate_miller = substrate_miller
         self.zslgen = zslgen or ZSLGenerator(bidirectional=True)
-
+        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()
 
@@ -131,14 +141,24 @@ def _find_terminations(self):
             reorient_lattice=False,  # This is necessary to not screw up the lattice
         )
 
-        film_slabs = film_sg.get_slabs()
-        sub_slabs = sub_sg.get_slabs()
-
+        film_slabs = film_sg.get_slabs(ftol=self.termination_ftol, filter_out_sym_slabs=self.filter_out_sym_slabs)
+        sub_slabs = sub_sg.get_slabs(ftol=self.termination_ftol, filter_out_sym_slabs=self.filter_out_sym_slabs)
         film_shifts = [slab.shift for slab in film_slabs]
-        film_terminations = [label_termination(slab) for slab in film_slabs]
+
+        if self.label_index:
+            film_terminations = [
+                label_termination(slab, self.termination_ftol, t_idx) for t_idx, slab in enumerate(film_slabs, start=1)
+            ]
+        else:
+            film_terminations = [label_termination(slab, self.termination_ftol) for slab in film_slabs]
 
         sub_shifts = [slab.shift for slab in sub_slabs]
-        sub_terminations = [label_termination(slab) for slab in sub_slabs]
+        if self.label_index:
+            sub_terminations = [
+                label_termination(slab, self.termination_ftol, t_idx) for t_idx, slab in enumerate(sub_slabs, start=1)
+            ]
+        else:
+            sub_terminations = [label_termination(slab, self.termination_ftol) for slab in sub_slabs]
 
         self._terminations = {
             (film_label, sub_label): (film_shift, sub_shift)

src/pymatgen/core/interface.py

@@ -2843,8 +2843,14 @@ def from_slabs(
         return iface
 
 
-def label_termination(slab: Structure) -> str:
-    """Label the slab surface termination."""
+def label_termination(slab: Structure, ftol: float = 0.25, t_idx: int | None = None) -> str:
+    """Label the slab surface termination.
+
+    Args:
+        slab (Slab): film or substrate slab to label termination for
+        ftol (float): tolerance for terminating position hierarchical clustering
+        t_idx (None | int): if not None, adding an extra index to the termination label output
+    """
     frac_coords = slab.frac_coords
     n = len(frac_coords)
 
@@ -2867,7 +2873,7 @@ def label_termination(slab: Structure) -> str:
 
     condensed_m = squareform(dist_matrix)
     z = linkage(condensed_m)
-    clusters = fcluster(z, 0.25, criterion="distance")
+    clusters = fcluster(z, ftol, criterion="distance")
 
     clustered_sites: dict[int, list[Site]] = {c: [] for c in clusters}
     for idx, cluster in enumerate(clusters):
@@ -2880,7 +2886,11 @@ def label_termination(slab: Structure) -> str:
 
     sp_symbol = SpacegroupAnalyzer(top_plane, symprec=0.1).get_space_group_symbol()
     form = top_plane.reduced_formula
-    return f"{form}_{sp_symbol}_{len(top_plane)}"
+
+    if t_idx is None:
+        return f"{form}_{sp_symbol}_{len(top_plane)}"
+
+    return f"{t_idx}_{form}_{sp_symbol}_{len(top_plane)}"
 
 
 def count_layers(struct: Structure, el: Element | None = None) -> int:

src/pymatgen/core/surface.py

@@ -1134,10 +1134,10 @@ def get_slab(
         if self.lll_reduce:
             # Sanitize Slab (LLL reduction + site sorting + map frac_coords)
             lll_slab = struct.copy(sanitize=True)
-            struct = lll_slab
 
             # Apply reduction on the scaling factor
             mapping = lll_slab.lattice.find_mapping(struct.lattice)
+            struct = lll_slab
             if mapping is None:
                 raise RuntimeError("LLL reduction has failed")
             scale_factor = np.dot(mapping[2], scale_factor)
@@ -1194,6 +1194,7 @@ def get_slabs(
         symmetrize: bool = False,
         repair: bool = False,
         ztol: float = 0,
+        filter_out_sym_slabs: bool = True,
     ) -> list[Slab]:
         """Generate slabs with shift values calculated from the internal
         gen_possible_terminations func. If the user decide to avoid breaking
@@ -1217,6 +1218,7 @@ def get_slabs(
                 can lead to many more possible slabs.
             ztol (float): Fractional tolerance for determine overlapping z-ranges,
                 smaller ztol might result in more possible Slabs.
+            filter_out_sym_slabs (bool): If True filter out identical slabs with different terminations.
 
         Returns:
             list[Slab]: All possible Slabs of a particular surface,
@@ -1342,22 +1344,25 @@ def get_z_ranges(
                 slabs.append(self.repair_broken_bonds(slab=slab, bonds=bonds))
 
         # Filter out surfaces that might be the same
-        matcher = StructureMatcher(ltol=tol, stol=tol, primitive_cell=False, scale=False)
+        if filter_out_sym_slabs:
+            matcher = StructureMatcher(ltol=tol, stol=tol, primitive_cell=False, scale=False)
+
+            final_slabs: list[Slab] = []
+            for group in matcher.group_structures(slabs):
+                # For each unique slab, symmetrize the
+                # surfaces by removing sites from the bottom
+                if symmetrize:
+                    sym_slabs = self.nonstoichiometric_symmetrized_slab(group[0])
+                    final_slabs.extend(sym_slabs)
+                else:
+                    final_slabs.append(group[0])
 
-        final_slabs: list[Slab] = []
-        for group in matcher.group_structures(slabs):
-            # For each unique slab, symmetrize the
-            # surfaces by removing sites from the bottom
+            # Filter out similar surfaces generated by symmetrization
             if symmetrize:
-                sym_slabs = self.nonstoichiometric_symmetrized_slab(group[0])
-                final_slabs.extend(sym_slabs)
-            else:
-                final_slabs.append(group[0])
-
-        # Filter out similar surfaces generated by symmetrization
-        if symmetrize:
-            matcher_sym = StructureMatcher(ltol=tol, stol=tol, primitive_cell=False, scale=False)
-            final_slabs = [group[0] for group in matcher_sym.group_structures(final_slabs)]
+                matcher_sym = StructureMatcher(ltol=tol, stol=tol, primitive_cell=False, scale=False)
+                final_slabs = [group[0] for group in matcher_sym.group_structures(final_slabs)]
+        else:
+            final_slabs = slabs
 
         return cast(list[Slab], sorted(final_slabs, key=lambda slab: slab.energy))
 

tests/analysis/interfaces/test_coherent_interface.py

@@ -1,5 +1,7 @@
 from __future__ import annotations
 
+import unittest
+
 from numpy.testing import assert_allclose
 
 from pymatgen.analysis.interfaces.coherent_interfaces import (
@@ -8,6 +10,9 @@
     get_2d_transform,
     get_rot_3d_for_2d,
 )
+from pymatgen.analysis.interfaces.substrate_analyzer import SubstrateAnalyzer
+from pymatgen.core.lattice import Lattice
+from pymatgen.core.structure import Structure
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 from pymatgen.util.testing import PymatgenTest
 
@@ -44,3 +49,31 @@ def test_coherent_interface_builder(self):
         # SP: this test is super fragile and the result fluctuates between 6, 30 and 42 for
         # no apparent reason. The author should fix this.
         assert len(list(builder.get_interfaces(termination=("O2_Pmmm_1", "Si_R-3m_1")))) >= 6
+
+
+class TestCoherentInterfaceBuilder(unittest.TestCase):
+    def setUp(self):
+        # build substrate & film structure
+        basis = [[0, 0, 0], [0.25, 0.25, 0.25]]
+        self.substrate = Structure(Lattice.cubic(a=5.431), ["Si", "Si"], basis)
+        self.film = Structure(Lattice.cubic(a=5.658), ["Ge", "Ge"], basis)
+
+    def test_termination_searching(self):
+        sub_analyzer = SubstrateAnalyzer()
+        matches = list(sub_analyzer.calculate(substrate=self.substrate, film=self.film))
+        cib = CoherentInterfaceBuilder(
+            film_structure=self.film,
+            substrate_structure=self.substrate,
+            film_miller=matches[0].film_miller,
+            substrate_miller=matches[0].substrate_miller,
+            zslgen=sub_analyzer,
+            termination_ftol=1e-4,
+            label_index=True,
+            filter_out_sym_slabs=False,
+        )
+        assert cib.terminations == [
+            ("1_Ge_P4/mmm_1", "1_Si_P4/mmm_1"),
+            ("1_Ge_P4/mmm_1", "2_Si_P4/mmm_1"),
+            ("2_Ge_P4/mmm_1", "1_Si_P4/mmm_1"),
+            ("2_Ge_P4/mmm_1", "2_Si_P4/mmm_1"),
+        ], "termination results wrong"

tests/core/test_surface.py

@@ -686,7 +686,7 @@ def test_previous_reconstructions(self):
             assert any(len(match.group_structures([struct, slab])) == 1 for slab in slabs)
 
 
-class MillerIndexFinderTests(PymatgenTest):
+class TestMillerIndexFinder(PymatgenTest):
     def setUp(self):
         self.cscl = Structure.from_spacegroup("Pm-3m", Lattice.cubic(4.2), ["Cs", "Cl"], [[0, 0, 0], [0.5, 0.5, 0.5]])
         self.Fe = Structure.from_spacegroup("Im-3m", Lattice.cubic(2.82), ["Fe"], [[0, 0, 0]])