feat: add SOAP-based symmetry-equivalent site detection (#10)

* feat: add SOAP-based symmetry-equivalent site detection

Add a second method for identifying symmetry non-equivalent adsorption
sites using SOAP (Smooth Overlap of Atomic Positions) descriptors from
the dscribe package with hierarchical clustering, as an alternative to
ASE's SymmetryEquivalenceCheck.

Core implementation (autoadsorbate/utils.py):
- _compute_soap_site_vectors(): compute SOAP descriptors for candidate sites
- _soap_similarity_matrix(): pairwise cosine similarity matrix
- _filter_unique_sites_by_soap(): cluster sites via hierarchical clustering
  or greedy merging, with configurable similarity_threshold (default 0.99)

Unified API (autoadsorbate/autoadsorbate.py):
- Surface.get_nonequivalent_sites() accepts method="ase"|"soap"
- Surface.sym_reduce() accepts method="ase"|"soap"
- New methods: get_soap_similarity_matrix(), get_nonequivalent_sites_soap(),
  soap_reduce()

Tests (tests/test_all.py):
- 13 tests covering SOAP similarity matrix, hcluster/greedy modes,
  soap_reduce, method switching, ASE vs SOAP comparison on Ni/Ru slab,
  timing benchmarks

Documentation (README.md):
- SOAP method usage section with API examples
- Threshold tuning guide
- Dedicated ASE vs SOAP benchmark section with:
  - Method comparison table (ASE: 16 sites/8.4s vs SOAP: 8 sites/0.016s)
  - Detailed breakdown by site type (8 SOAP clusters)
  - Site comparison plot (site_comparison.png)
  - SOAP distance histogram (soap_histogram.png)
  - Discussion of why SOAP outperforms ASE on broken-symmetry surfaces

Scripts:
- scripts/compare_sites.py: detailed ASE vs SOAP comparison
- scripts/plot_sites.py: site visualization with color-coded SOAP clusters
- scripts/site_comparison.png: generated comparison image
- scripts/soap_histogram.png: SOAP distance distribution plot

On Ni(111) 3x3 slab with one Ru dopant:
- ASE finds 16 unique sites in ~8.4 s
- SOAP finds 8 unique sites in ~0.016 s (~525x speedup)

* update version

* updated dependency

---------

Co-authored-by: Sandip De <desa@basfad.basf.net>

Author: sandipde

.gitignore

@@ -401,8 +401,34 @@ plt.show()
 
 
 
-We can reduce the complete list of sites based on symmetry (```ase.utils.structure_comparator.SymmetryEquivalenceCheck```):
+We can reduce the complete list of sites based on symmetry. Two methods are available:
 
+1. **ASE crystallographic symmetry** (default) – uses `ase.utils.structure_comparator.SymmetryEquivalenceCheck`:
+
+```python
+s.sym_reduce()                          # equivalent to s.sym_reduce(method='ase')
+s.site_df
+```
+
+2. **SOAP descriptor similarity** – uses SOAP descriptors ([DScribe](https://singroup.github.io/dscribe/)) and hierarchical clustering, robust for disordered / low-symmetry surfaces:
+
+```python
+# Requires: pip install dscribe scikit-learn scipy
+s.sym_reduce(method='soap')             # default similarity_threshold=0.99
+s.site_df
+```
+
+The SOAP method can be further tuned:
+```python
+s.sym_reduce(
+    method='soap',
+    similarity_threshold=0.95,           # lower = more aggressive merging
+    soap_params={'r_cut': 6.0, 'n_max': 8, 'l_max': 6, 'sigma': 0.1},
+    soap_cluster_method='hcluster',      # or 'greedy'
+)
+```
+
+Using the default (ASE) method:
 
 ```python
 s.sym_reduce()
@@ -496,6 +522,89 @@ plot_atoms(s.view_surface(return_atoms=True))
 
 
 
+## Benchmark: ASE vs SOAP symmetry reduction
+
+The two symmetry-reduction methods differ fundamentally in how they define "equivalent":
+
+| Aspect | ASE (`SymmetryEquivalenceCheck`) | SOAP (dscribe + hierarchical clustering) |
+|---|---|---|
+| **Equivalence criterion** | Strict crystallographic space-group symmetry | Cosine similarity of local SOAP descriptors |
+| **Adjustable threshold** | No — binary match/no-match | Yes — continuous `similarity_threshold` |
+| **Scaling** | O(n²) pairwise structure comparisons | O(n) SOAP evaluations + O(n²) similarity matrix |
+| **Robustness to disorder** | Breaks down when symmetry is broken | Groups by local chemical environment |
+
+### Ni(111) 3×3 slab with 1 Ni → Ru substitution
+
+To illustrate the difference, we replace one surface Ni with Ru on a Ni(111) 3×3 slab (54 candidate adsorption sites) and compare both methods:
+
+```python
+import copy, time
+from ase.build import fcc111
+from autoadsorbate import Surface
+
+slab = fcc111("Ni", (3, 3, 3), periodic=True, vacuum=10)
+top_z = slab.positions[:, 2].max()
+for atom in slab:
+    if abs(atom.position[2] - top_z) < 0.1:
+        atom.symbol = "Ru"
+        break
+
+s = Surface(slab)
+s_ase = copy.deepcopy(s)
+s_soap = copy.deepcopy(s)
+
+s_ase.sym_reduce(method="ase")
+s_soap.sym_reduce(method="soap", similarity_threshold=0.99)
+```
+
+#### Timing results
+
+| Method | Unique sites | Time | Speedup |
+|---|---:|---:|---:|
+| ASE `SymmetryEquivalenceCheck` | 16 | 8.4 s | 1× |
+| SOAP + hierarchical clustering | 8 | 0.016 s | **~525×** |
+
+#### Detailed breakdown by site type
+
+| SOAP Cluster | Site type | Formula | Total sites | ASE reps | SOAP rep | Interpretation |
+|---|---|---|---:|---:|---:|---|
+| 8 | atop | {Ru: 1} | 1 | 1 | 1 | Unique Ru atop — both agree |
+| 1 | atop | {Ni: 1} | 8 | 3 | 1 | ASE splits by distance-to-Ru; SOAP merges (all Ni atop) |
+| 2 | bridge | {Ru:1, Ni:1} | 6 | 1 | 1 | Ru-Ni bridges — both agree |
+| 3 | bridge | {Ni: 2} | 21 | 5 | 1 | ASE splits into 5; SOAP merges all Ni-Ni bridges |
+| 4 | hollow | {Ru:1, Ni:2} | 3 | 1 | 1 | Near-Ru hollows — both agree |
+| 6 | hollow | {Ru:1, Ni:2} | 3 | 1 | 1 | Far-Ru hollows — both agree |
+| 5 | hollow | {Ni: 3} | 6 | 2 | 1 | ASE splits into 2; SOAP merges |
+| 7 | hollow | {Ni: 3} | 6 | 2 | 1 | ASE splits into 2; SOAP merges |
+
+The SOAP method captures the **physically meaningful site diversity** (8 distinct local environments) while the ASE method finds 16 sites that differ only by their distance from the Ru dopant within an otherwise identical coordination shell.
+
+#### Site map
+
+![Site comparison](README_files/site_comparison.png)
+
+*Left: 16 ASE representatives. Right: 8 SOAP representatives. Coloured by SOAP cluster; marker shape = site type (○ atop, □ bridge, △ hollow). Grey dots = all 54 candidate sites.*
+
+#### SOAP distance histogram
+
+```python
+sim = s.get_soap_similarity_matrix()
+dist = 1.0 - sim
+upper = dist[np.triu_indices_from(dist, k=1)]
+
+fig, ax = plt.subplots(figsize=(7, 3.5))
+ax.hist(upper, bins=60, edgecolor="black", linewidth=0.4, color="#4C72B0")
+ax.axvline(0.01, color="red", ls="--", lw=1.5,
+           label="threshold = 0.01\n(similarity = 0.99)")
+ax.set_xlabel("SOAP distance  (1 − cosine similarity)")
+ax.set_ylabel("Number of site pairs")
+ax.legend(fontsize=9)
+```
+
+![SOAP distance histogram](README_files/soap_histogram.png)
+
+The histogram shows a clear separation between intra-cluster pairs (distance ≈ 0) and inter-cluster pairs, confirming that the 0.99 similarity threshold sits in the natural gap between equivalent and non-equivalent site pairs.
+
 ## Making surogate SMILES automatically
 
 Simple methods of brute force SMILES enumeration are implemented as well. For example, only using a few lines of code we can initialize multiple conformers of all reaction intermediates in the nitrogen hydrogenation reaction. A template of the required information can be found here:
@@ -678,7 +787,8 @@ from autoadsorbate import Surface, Fragment
 
 slab = fcc211(symbol = 'Cu', size=(6,3,3), vacuum=10)  # any ase.Atoms object
 s=Surface(slab, touch_sphere_size=2.7)                 # finding all surface atoms
-s.sym_reduce()                                         # keeping only non-identical sites
+s.sym_reduce()                                         # keeping only non-identical sites (default: method='ase')
+# s.sym_reduce(method='soap')                          # alternative: SOAP-descriptor based reduction
 
 fragments = [
     Fragment('S1S[OH+]CC(N)[OH+]1', to_initialize=20), # For each *SMILES we can request a differnet number of conformers 

README.md

@@ -2,7 +2,7 @@
 
 __author__ = """Fakoe Edvin"""
 __email__ = "edvinfako@gmail.com"
-__version__ = "0.2.5"
+__version__ = "0.2.6"
 
 from autoadsorbate.autoadsorbate import Fragment, Surface
 from autoadsorbate.Smile import get_marked_smiles

README_files/site_comparison.png

@@ -18,6 +18,9 @@
 from .utils import (
     get_sorted_by_snap_dist,
     make_site_info_writable,
+    _compute_soap_site_vectors,
+    _soap_similarity_matrix,
+    _filter_unique_sites_by_soap,
 )
 
 from .Particle import get_shrinkwrap_particle_ads_sites
@@ -393,13 +396,53 @@ def compare_sites(self, site_index1: int, site_index2: int, **kwargs) -> bool:
 
         return SEC.compare(self.atoms + site1, self.atoms + site2)
 
-    def get_nonequivalent_sites(self, **kwargs) -> List[int]:
-        """
-        Returns a list of indices for nonequivalent sites.
+    def get_nonequivalent_sites(
+        self,
+        method: str = "ase",
+        similarity_threshold: float = 0.99,
+        soap_params: dict = None,
+        soap_cluster_method: str = "hcluster",
+        **kwargs,
+    ) -> List[int]:
+        """Return a list of indices for nonequivalent sites.
+
+        Args:
+            method (str, optional): Algorithm used to determine equivalence.
+                ``"ase"`` (default) – ASE's crystallographic
+                ``SymmetryEquivalenceCheck`` (pairwise comparison, exact
+                space-group symmetry).  ``"soap"`` – SOAP descriptor
+                similarity with hierarchical clustering (robust for
+                disordered / low-symmetry surfaces).  Any extra
+                ``**kwargs`` are forwarded to
+                ``SymmetryEquivalenceCheck`` when *method="ase"*.
+            similarity_threshold (float, optional): Cosine-similarity
+                threshold (only used when *method="soap"*).  Defaults
+                to ``0.99``.
+            soap_params (dict, optional): SOAP hyper-parameters forwarded
+                to ``dscribe.descriptors.SOAP`` (only used when
+                *method="soap"*).
+            soap_cluster_method (str, optional): ``"hcluster"`` or
+                ``"greedy"`` (only used when *method="soap"*).
 
         Returns:
             List[int]: A list of indices for nonequivalent sites.
         """
+        if method == "ase":
+            return self._get_nonequivalent_sites_ase(**kwargs)
+        elif method == "soap":
+            return self.get_nonequivalent_sites_soap(
+                similarity_threshold=similarity_threshold,
+                soap_params=soap_params,
+                method=soap_cluster_method,
+            )
+        else:
+            raise ValueError(
+                f"Unknown method '{method}'. Choose 'ase' or 'soap'."
+            )
+
+    def _get_nonequivalent_sites_ase(self, **kwargs) -> List[int]:
+        """Return nonequivalent-site indices using ASE's
+        ``SymmetryEquivalenceCheck`` (legacy behaviour)."""
         original = []
         i_s = self.site_df.index.values
         matches = np.array([False for _ in i_s])
@@ -414,11 +457,142 @@ def get_nonequivalent_sites(self, **kwargs) -> List[int]:
                 break
         return original
 
-    def sym_reduce(self, **kwargs):
+    def sym_reduce(
+        self,
+        method: str = "ase",
+        similarity_threshold: float = 0.99,
+        soap_params: dict = None,
+        soap_cluster_method: str = "hcluster",
+        **kwargs,
+    ):
+        """Reduce the site DataFrame to nonequivalent sites.
+
+        Args:
+            method (str, optional): ``"ase"`` (default) for
+                crystallographic symmetry, ``"soap"`` for SOAP-descriptor
+                based clustering.  See :meth:`get_nonequivalent_sites`
+                for full parameter descriptions.
+            similarity_threshold (float, optional): SOAP cosine-similarity
+                threshold (only when *method="soap"*).
+            soap_params (dict, optional): SOAP hyper-parameters (only when
+                *method="soap"*).
+            soap_cluster_method (str, optional): ``"hcluster"`` or
+                ``"greedy"`` (only when *method="soap"*).
+            **kwargs: Extra arguments forwarded to
+                ``SymmetryEquivalenceCheck`` when *method="ase"*.
         """
-        Reduces the site DataFrame to nonequivalent sites.
+        include = self.get_nonequivalent_sites(
+            method=method,
+            similarity_threshold=similarity_threshold,
+            soap_params=soap_params,
+            soap_cluster_method=soap_cluster_method,
+            **kwargs,
+        )
+        include_filter = [i in include for i in self.site_df.index.values]
+        self.site_df = self.site_df[include_filter]
+        self.site_dict = self.site_df.to_dict(orient="list")
+
+    # ------------------------------------------------------------------
+    # SOAP-descriptor based symmetry reduction
+    # ------------------------------------------------------------------
+
+    def get_soap_similarity_matrix(
+        self,
+        soap_params: dict = None,
+        probe_element: str = "X",
+    ) -> np.ndarray:
+        """Return the pairwise cosine-similarity matrix for all sites using
+        SOAP descriptors (via *dscribe*).
+
+        A ghost probe atom is placed at each site position and SOAP descriptors
+        are evaluated within the periodic slab environment.
+
+        Args:
+            soap_params (dict, optional): SOAP hyper-parameters forwarded to
+                ``dscribe.descriptors.SOAP``.  Defaults to
+                ``{r_cut: 5.0, n_max: 8, l_max: 6, sigma: 0.1}``.
+            probe_element (str, optional): Element for the probe atom
+                (must not be present in the slab). Defaults to ``"X"``.
+
+        Returns:
+            np.ndarray: Symmetric similarity matrix of shape
+            ``(n_sites, n_sites)`` with values in ``[-1, 1]``.
+
+        Requires:
+            ``dscribe`` (``pip install dscribe``).
         """
-        include = self.get_nonequivalent_sites(**kwargs)
+        soap_vectors = _compute_soap_site_vectors(
+            self.atoms, self.site_df,
+            soap_params=soap_params,
+            probe_element=probe_element,
+        )
+        return _soap_similarity_matrix(soap_vectors)
+
+    def get_nonequivalent_sites_soap(
+        self,
+        similarity_threshold: float = 0.99,
+        soap_params: dict = None,
+        method: str = "hcluster",
+    ) -> List[int]:
+        """Return indices of non-equivalent sites determined by SOAP
+        descriptor similarity and hierarchical clustering.
+
+        This is an alternative to :meth:`get_nonequivalent_sites` which relies
+        on ASE's crystallographic ``SymmetryEquivalenceCheck``.  The SOAP
+        approach is more robust for disordered or low-symmetry surfaces and
+        allows continuous tuning of the similarity threshold.
+
+        Args:
+            similarity_threshold (float, optional): Cosine similarity above
+                which two sites are deemed equivalent.  Defaults to ``0.99``.
+            soap_params (dict, optional): SOAP hyper-parameters forwarded to
+                ``dscribe.descriptors.SOAP``.  Defaults to
+                ``{r_cut: 5.0, n_max: 8, l_max: 6, sigma: 0.1}``.
+            method (str, optional): Clustering algorithm.  ``"hcluster"``
+                (default) uses agglomerative hierarchical clustering
+                (``scipy.cluster.hierarchy``);  ``"greedy"`` uses the legacy
+                greedy merging approach.
+
+        Returns:
+            List[int]: DataFrame index labels of one representative per
+            equivalence class.
+
+        Requires:
+            ``dscribe``, ``scikit-learn``, and (for *method="hcluster"*)
+            ``scipy``.
+        """
+        reduced = _filter_unique_sites_by_soap(
+            slab=self.atoms,
+            site_df=self.site_df,
+            soap_params=soap_params,
+            similarity_threshold=similarity_threshold,
+            method=method,
+        )
+        return list(reduced.index)
+
+    def soap_reduce(
+        self,
+        similarity_threshold: float = 0.99,
+        soap_params: dict = None,
+        method: str = "hcluster",
+    ):
+        """Reduce the site DataFrame to non-equivalent sites using SOAP
+        descriptors and (hierarchical) clustering.
+
+        This is the SOAP analogue of :meth:`sym_reduce`.
+
+        Args:
+            similarity_threshold (float, optional): Cosine similarity above
+                which two sites are deemed equivalent.  Defaults to ``0.99``.
+            soap_params (dict, optional): SOAP hyper-parameters forwarded to
+                ``dscribe.descriptors.SOAP``.
+            method (str, optional): ``"hcluster"`` (default) or ``"greedy"``.
+        """
+        include = self.get_nonequivalent_sites_soap(
+            similarity_threshold=similarity_threshold,
+            soap_params=soap_params,
+            method=method,
+        )
         include_filter = [i in include for i in self.site_df.index.values]
         self.site_df = self.site_df[include_filter]
         self.site_dict = self.site_df.to_dict(orient="list")