Small cluster issue (#26)

* Trying out a fix for the small cluster problem.

* Updated formatting.

* Fixed formatting.

* Simplified logic.

* Improvement.

* Fixed linting.

* Updated workflows.

Author: lauri-codes

.github/workflows/wheels.yml

@@ -14,7 +14,7 @@ jobs:
       - uses: actions/checkout@v4
 
       - name: Build wheels
-        uses: pypa/cibuildwheel@v2.18.0
+        uses: pypa/cibuildwheel@v2.20.0
         env:
           CIBW_SKIP: cp36-*
 

matid/core/periodicfinder.py

@@ -124,7 +124,7 @@ def get_region(
             neighbour_factors,
             search_mask,
         ) = self._find_possible_bases(system, seed_index)
-        proto_cell, offset, dim = self._find_proto_cell(
+        proto_cell, offset, dim, n_periodic_spans_selected = self._find_proto_cell(
             system,
             seed_index,
             possible_spans,
@@ -153,7 +153,11 @@ def get_region(
 
             i_indices = unit_collection.get_basis_indices()
 
-            if len(i_indices) > 0:
+            # The region is accepted if it either has more than 1+dim atoms in
+            # it or one of the spans is periodic. In many cases the prototype
+            # cell can be found (seed + number of spans), but it cannot really
+            # be extended to cover any significant portion of the system.
+            if len(i_indices) > 1 + dim or n_periodic_spans_selected >= 1:
                 region = unit_collection
                 region._pos_tol = pos_tol
 
@@ -223,6 +227,8 @@ def _find_proto_cell(
         Returns:
             ase.Atoms: A system representing the best cell that was found
             np.ndarray: Position of the seed atom in the cell
+            int: Dimensionality of the found system
+            int: Number of periodic unit cell vectors used
         """
         positions = system.get_positions()
         numbers = system.get_atomic_numbers()
@@ -342,18 +348,26 @@ def _find_proto_cell(
             | (metric >= 0.75 * n_neighbours)
         )[0]
 
-        if len(valid_span_indices) == 0:
-            return None, None, None
+        total_valid_spans = len(valid_span_indices)
+        if total_valid_spans == 0:
+            return None, None, None, None
 
         # Find the best basis
         valid_span_metrics = metric[valid_span_indices]
         valid_spans = possible_spans[valid_span_indices]
         best_combo = self._find_best_basis(valid_spans, valid_span_metrics)
         dim = len(best_combo)
 
+        # Check how many of the periodic spans are still selected as prototype
+        # unit cell vectors
+        selected_spans = range(total_valid_spans - n_periodic_spans, total_valid_spans)
+        n_periodic_spans_selected = sum(
+            span_index in best_combo for span_index in selected_spans
+        )
+
         # Currently 1D is not handled
         if dim == 1:
-            return None, None, 1
+            return None, None, 1, None
 
         best_spans = valid_spans[best_combo]
         n_spans = len(best_spans)
@@ -382,7 +396,7 @@ def _find_proto_cell(
 
         # If the seed atom is not in a valid graph, no region could be found.
         if seed_group_index is None:
-            return None, None, None
+            return None, None, None, None
 
         # Notice that the seed group index can get updated by the cell search if
         # an atom is dropped out of the cell due to appearing too infrequently.
@@ -409,7 +423,7 @@ def _find_proto_cell(
                 pos_tol,
             )
         if proto_cell is None:
-            return None, None, None
+            return None, None, None, None
 
         two_valid_spans = n_spans == 2
         if n_spans == 3:
@@ -424,7 +438,7 @@ def _find_proto_cell(
                     proto_cell, bond_threshold
                 )
             except MatIDError:
-                return None, None, None
+                return None, None, None, None
             if dimensionality != 3:
                 # Try if the cell can be "reduced" to a 2D material
                 if dimensionality == 2:
@@ -458,20 +472,18 @@ def _find_proto_cell(
                     two_valid_spans = True
                     n_spans = 2
                 else:
-                    return None, None, None
+                    return None, None, None, None
 
         if two_valid_spans:
             # If the best 2D vectors consists only of the simulation basis cell
             # vectors, check that these vectors are below a predefined size.
             # Otherwise the cell cannot be accepted because there is not enough
             # statistics about the cell contents to distinguish outliers.
             if n_periodic_spans > 0:
-                periodic_span_indices = valid_span_indices[-n_periodic_spans:]
-                best_span_ind = valid_span_indices[best_combo]
-                if set(best_span_ind).issubset(set(periodic_span_indices)):
+                if n_periodic_spans_selected == 2:
                     cell_lens = np.linalg.norm(best_spans, axis=1)
                     if np.any(cell_lens > self.max_2d_single_cell_size):
-                        return None, None, None
+                        return None, None, None, None
 
             # Check the dimensionality
             dimensionality, cluster_labels = matid.geometry.get_dimensionality(
@@ -482,7 +494,7 @@ def _find_proto_cell(
                 # has multiple stacked 2D sheets with identical periodicity. In
                 # this case the unit cell should only comprise of atoms in the
                 # cluster where the seed atom is in.
-                for i_index, i_cluster in enumerate(cluster_labels):
+                for i_cluster in cluster_labels:
                     try:
                         seed_group_index = i_cluster.index(seed_group_index)
                     except ValueError:
@@ -498,13 +510,13 @@ def _find_proto_cell(
                     proto_cell, bond_threshold
                 )
                 if dimensionality is None:
-                    return None, None, None
+                    return None, None, None, None
                 else:
                     if dimensionality != 2:
-                        return None, None, None
+                        return None, None, None, None
             else:
                 if dimensionality != 2:
-                    return None, None, None
+                    return None, None, None, None
 
             # Check the cell thickness. 2D materials that are thicker than a
             # specified threshold are not accepted.
@@ -514,16 +526,16 @@ def _find_proto_cell(
             offset = proto_cell.get_positions()[seed_group_index]
             thickness = matid.geometry.get_thickness(proto_cell, 2)
             if thickness > self.max_2d_cell_height:
-                return None, None, None
+                return None, None, None, None
 
         # Check that the final proto cell atoms don't overlap
         if proto_cell is not None:
             dist_proto_cell = matid.geometry.get_distances(proto_cell).dist_matrix_mic
             dist_proto_cell = dist_proto_cell[np.triu_indices(dist_proto_cell.shape[0])]
             if dist_proto_cell.min() < overlap_threshold:
-                return None, None, None
+                return None, None, None, None
 
-        return proto_cell, offset, n_spans
+        return proto_cell, offset, n_spans, n_periodic_spans_selected
 
     def _find_graphs(
         self,

matid/symmetry/symmetryanalyzer.py

@@ -1139,8 +1139,7 @@ def _find_wyckoff_ground_state(
                     i_n = r["wyckoff_positions"].get((w, z))
                     if i_n is not None:
                         n_atoms_map[i_n].append(r)
-                        if i_n > n_atoms_max:
-                            n_atoms_max = i_n
+                        n_atoms_max = max(i_n, n_atoms_max)
                 if n_atoms_max != 0:
                     representations = n_atoms_map[n_atoms_max]
                 if len(representations) == 1:

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = 'matid'
-version = '2.1.0'
+version = '2.1.1'
 description = 'MatID is a Python package for identifying and analyzing atomistic systems based on their structure.'
 readme = "README.md"
 authors = [{ name = "Lauri Himanen" }]

tests/clustering/test_sbc.py

@@ -5,8 +5,9 @@
 import ase.io
 from ase import Atoms
 from ase.build import bulk
+from ase.visualize import view
 
-from conftest import surface, stack, rattle, assert_topology
+from conftest import surface, stack, rattle, assert_topology, create_graphene, create_fe
 from matid.clustering import SBC, Cluster
 
 
@@ -302,6 +303,53 @@ def test_clustering_coincidental_patterns(system, clusters_expected):
     assert_topology(results, clusters_expected)
 
 
+@pytest.mark.parametrize(
+    "system, clusters_expected",
+    [
+        pytest.param(
+            create_fe(),
+            [
+                Cluster(
+                    [0, 1],
+                    dimensionality=3,
+                ),
+            ],
+            id="clusters that consist only of the 3D prototype cell but contain a periodic span should be accepted",
+        ),
+        pytest.param(
+            create_graphene(),
+            [
+                Cluster(
+                    [0, 1],
+                    dimensionality=2,
+                ),
+            ],
+            id="clusters that consist only of the 2D prototype cell but contain a periodic span should be accepted",
+        ),
+        pytest.param(
+            ase.build.surface(
+                ase.build.bulk(
+                    "MoC", crystalstructure="rocksalt", a=4.38, b=4.38, c=4.38, alpha=90
+                ),
+                (0, 0, 1),
+                5,
+                vacuum=5,
+                tol=1e-10,
+                periodic=True,
+            )
+            * [3, 3, 1],
+            [],
+            id="clusters that consist only of the prototype cell but don't contain a periodic span should be ignored",
+        ),
+    ],
+)
+def test_small_clusters(system, clusters_expected):
+    """Tests that small clusters that do not cover"""
+    system = rattle(system, 0.1)
+    results = SBC().get_clusters(system, pos_tol=0.01)
+    assert_topology(results, clusters_expected)
+
+
 fcc_clusters = [Cluster(range(len(surface_fcc)), dimensionality=0)]
 single_atom_clusters = []
 

tests/performance/performance.py

@@ -198,14 +198,10 @@ def plot(show):
         i_timemin = times.min()
         i_nmax = n_atoms.max()
         i_nmin = n_atoms.min()
-        if i_timemax > timemax:
-            timemax = i_timemax
-        if i_timemin < timemin:
-            timemin = i_timemin
-        if i_nmax > nmax:
-            nmax = i_nmax
-        if i_nmin < nmin:
-            nmin = i_nmin
+        timemax = max(i_timemax, timemax)
+        timemin = min(i_timemin, timemin)
+        nmax = max(i_nmax, nmax)
+        nmin = min(i_nmin, nmin)
         # ax1.fill_between(n_atoms, times_mean - times_std, times_mean + times_std, color=color, alpha=0.3)
         ax1.plot(
             n_atoms,