Fix lll_reduce for slab generation (materialsproject#3927)

* Update surface and interface generating functions

1. fixing problem for the lll_reduce process when making slabs, doing mapping before updating the structure
2. allow to set ftol of the termination distances for hierarchical cluster so that some non-identical terminations close to each other can be identified
3. allow to add index for terminations so that terminations with the same space group can be distinguished

Interfaces made by identical slabs can be non-identical because the relative transformation of the misorientation and the termination variation do not ensure symmetry, especially when the film and substrate have different point groups. Therefore, the termination finding function should allow to generate all the possible terminations. This can help others to develop more robust algorithm to group the equivalent interfaces made by different terminations.

---------

Signed-off-by: Jason Xie <[email protected]>
Signed-off-by: Jason Xie <[email protected]>
Co-authored-by: Janosh Riebesell <[email protected]>

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]])