diff --git a/.gitignore b/.gitignore index 0a19790..81db081 100644 Binary files a/.gitignore and b/.gitignore differ diff --git a/README.md b/README.md index 120d32c..960ae39 100644 --- a/README.md +++ b/README.md @@ -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 diff --git a/README_files/site_comparison.png b/README_files/site_comparison.png new file mode 100644 index 0000000..678cde6 Binary files /dev/null and b/README_files/site_comparison.png differ diff --git a/README_files/soap_histogram.png b/README_files/soap_histogram.png new file mode 100644 index 0000000..96babbc Binary files /dev/null and b/README_files/soap_histogram.png differ diff --git a/autoadsorbate/__init__.py b/autoadsorbate/__init__.py index 136659f..e94b3ae 100644 --- a/autoadsorbate/__init__.py +++ b/autoadsorbate/__init__.py @@ -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 diff --git a/autoadsorbate/autoadsorbate.py b/autoadsorbate/autoadsorbate.py index fcd40aa..ac80099 100644 --- a/autoadsorbate/autoadsorbate.py +++ b/autoadsorbate/autoadsorbate.py @@ -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") diff --git a/autoadsorbate/utils.py b/autoadsorbate/utils.py index 8257016..8b299e4 100644 --- a/autoadsorbate/utils.py +++ b/autoadsorbate/utils.py @@ -661,66 +661,158 @@ def get_max_atoms_dict(traj): max_atoms_dict[k] = max(v) return max_atoms_dict -def _filter_unique_sites_by_soap( +def _compute_soap_site_vectors( slab: Atoms, site_df: pd.DataFrame, - cutoff: float = 5.0, soap_params: dict = None, - similarity_threshold: float = 0.999 -) -> pd.DataFrame: - - """ - special helper function for handling edge cases where ase symmetry checker gives unsatisfactory results. - requires additional dependecies: sklearn and dscribe. + probe_element: str = "X", +) -> np.ndarray: + """Compute SOAP descriptor vectors centred at each adsorption site. + + A ghost probe atom is temporarily placed at each site position and SOAP + descriptors are evaluated for that position within the periodic slab + environment. + + Args: + slab: The surface slab (ASE Atoms, with cell and PBC). + site_df: DataFrame with a ``coordinates`` column (list of 3-vectors). + soap_params: Optional dict forwarded to ``dscribe.descriptors.SOAP``. + Defaults to ``{r_cut: 5.0, n_max: 8, l_max: 6, sigma: 0.1}``. + probe_element: Element symbol used for the probe atom. Must **not** + already be present in *slab*. Defaults to ``"X"``. + + Returns: + np.ndarray: SOAP feature matrix of shape ``(n_sites, n_features)``. + + Requires: + ``dscribe`` (pip install dscribe). """ - from dscribe.descriptors import SOAP - from sklearn.metrics.pairwise import cosine_similarity + from dscribe.descriptors import SOAP if soap_params is None: soap_params = { - "r_cut": cutoff, + "r_cut": 5.0, "n_max": 8, "l_max": 6, "sigma": 0.1, - "average": "off", } + + species = sorted(set(slab.get_chemical_symbols()) | {probe_element}) + soap = SOAP( - species=slab.get_chemical_symbols(), + species=species, periodic=True, - **soap_params + **soap_params, ) - # Compute SOAP for all atoms in slab - all_soap_vectors = soap.create(slab) # shape (num_atoms, soap_vector_length) - - coords = np.array(site_df['coordinates'].tolist()) - atom_positions = slab.get_positions() - - # Find closest atom index for each site coordinate - indices = [] - for c in coords: - dists = np.linalg.norm(atom_positions - c, axis=1) - closest_index = np.argmin(dists) - indices.append(closest_index) - - # Extract SOAP vectors for closest atoms - soap_vectors = all_soap_vectors[indices] - - # Compute similarity and cluster - similarity_matrix = cosine_similarity(soap_vectors) - n_sites = len(site_df) - seen = np.zeros(n_sites, dtype=bool) - unique_indices = [] - - for i in range(n_sites): - if seen[i]: - continue - unique_indices.append(i) - similar_sites = np.where(similarity_matrix[i] >= similarity_threshold)[0] - for j in similar_sites: - seen[j] = True - - return site_df.iloc[unique_indices] #.reset_index(drop=True) + coords = np.array(site_df["coordinates"].tolist()) + n_slab = len(slab) + soap_vectors = [] + + for coord in coords: + probe = slab.copy() + probe.append(Atom(probe_element, position=coord)) + # Compute SOAP only for the probe atom (last one) + desc = soap.create(probe, centers=[n_slab]) + soap_vectors.append(desc.flatten()) + + return np.array(soap_vectors) + + +def _soap_similarity_matrix( + soap_vectors: np.ndarray, +) -> np.ndarray: + """Pairwise cosine-similarity matrix from SOAP feature vectors. + + Args: + soap_vectors: Feature matrix of shape ``(n, d)``. + + Returns: + np.ndarray: Symmetric similarity matrix of shape ``(n, n)`` with + values in ``[-1, 1]``. + """ + from sklearn.metrics.pairwise import cosine_similarity + + return cosine_similarity(soap_vectors) + + +def _filter_unique_sites_by_soap( + slab: Atoms, + site_df: pd.DataFrame, + cutoff: float = 5.0, + soap_params: dict = None, + similarity_threshold: float = 0.99, + method: str = "hcluster", +) -> pd.DataFrame: + """Filter *site_df* to keep only symmetry-non-equivalent sites using + SOAP descriptors and (optionally) hierarchical clustering. + + This is a helper for cases where ASE's ``SymmetryEquivalenceCheck`` + gives unsatisfactory results. Requires ``dscribe`` and ``sklearn``. + + Args: + slab: The surface slab. + site_df: DataFrame produced by ``Surface`` (must contain a + ``coordinates`` column). + cutoff: SOAP radial cutoff (used only when *soap_params* is ``None``). + soap_params: Full dict of SOAP hyper-parameters forwarded to + ``_compute_soap_site_vectors``. + similarity_threshold: Cosine-similarity above which two sites are + considered equivalent. For *method="hcluster"* this is converted + to a distance threshold ``1 - similarity_threshold``. + method: ``"hcluster"`` (default) — agglomerative hierarchical + clustering via ``scipy.cluster.hierarchy``; ``"greedy"`` — simple + greedy merging (legacy behaviour). + + Returns: + pd.DataFrame: Subset of *site_df* with one representative per cluster. + """ + if soap_params is None: + soap_params = { + "r_cut": cutoff, + "n_max": 8, + "l_max": 6, + "sigma": 0.1, + } + + soap_vectors = _compute_soap_site_vectors(slab, site_df, soap_params=soap_params) + similarity_matrix = _soap_similarity_matrix(soap_vectors) + + if method == "hcluster": + from scipy.cluster.hierarchy import fcluster, linkage + from scipy.spatial.distance import squareform + + # Convert similarity → distance, clip for numerical safety + dist_matrix = np.clip(1.0 - similarity_matrix, 0.0, 2.0) + np.fill_diagonal(dist_matrix, 0.0) + condensed = squareform(dist_matrix, checks=False) + Z = linkage(condensed, method="average") + labels = fcluster(Z, t=1.0 - similarity_threshold, criterion="distance") + + # Pick the first site in each cluster as the representative + unique_indices = [] + seen_labels = set() + for idx, lab in enumerate(labels): + if lab not in seen_labels: + seen_labels.add(lab) + unique_indices.append(idx) + + elif method == "greedy": + n_sites = len(site_df) + seen = np.zeros(n_sites, dtype=bool) + unique_indices = [] + for i in range(n_sites): + if seen[i]: + continue + unique_indices.append(i) + similar_sites = np.where(similarity_matrix[i] >= similarity_threshold)[0] + for j in similar_sites: + seen[j] = True + + else: + raise ValueError(f"Unknown method '{method}'. Use 'hcluster' or 'greedy'.") + + return site_df.iloc[unique_indices] diff --git a/pyproject.toml b/pyproject.toml index 6fd1621..bd8bce0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,8 +11,10 @@ license = "MIT" dependencies = [ "ase>=3.24.0", + "dscribe>=2.1.0", "pandas>=2.0.0", "rdkit>=2022.9.5", + "scikit-learn>=1.3.0", ] [dependency-groups] diff --git a/scripts/compare_sites.py b/scripts/compare_sites.py new file mode 100644 index 0000000..740a156 --- /dev/null +++ b/scripts/compare_sites.py @@ -0,0 +1,86 @@ +"""Detailed comparison of ASE vs SOAP site reduction on Ni/Ru slab.""" +import copy +import numpy as np +from ase.build import fcc111 +from autoadsorbate import Surface +from autoadsorbate.utils import _compute_soap_site_vectors, _soap_similarity_matrix +from scipy.cluster.hierarchy import fcluster, linkage +from scipy.spatial.distance import squareform + +# Build Ni/Ru slab +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) + +ase_idx = set(s_ase.site_df.index) +soap_idx = set(s_soap.site_df.index) + +# Compute SOAP cluster labels for all sites +df = s.site_df.copy() +df["in_ASE"] = [i in ase_idx for i in df.index] +df["in_SOAP"] = [i in soap_idx for i in df.index] +df["formula_str"] = [str(f) for f in df["site_formula"]] + +vecs = _compute_soap_site_vectors(s.atoms, s.site_df) +sim = _soap_similarity_matrix(vecs) +dist = np.clip(1.0 - sim, 0.0, 2.0) +np.fill_diagonal(dist, 0.0) +Z = linkage(squareform(dist, checks=False), method="average") +labels = fcluster(Z, t=0.01, criterion="distance") +df["soap_cluster"] = labels + +print("=" * 90) +print("ALL SITES — grouped by SOAP cluster") +print("=" * 90) +for cl in sorted(df["soap_cluster"].unique()): + members = df[df["soap_cluster"] == cl] + rep_ase = [i for i in members.index if i in ase_idx] + rep_soap = [i for i in members.index if i in soap_idx] + print(f"\nSOAP Cluster {cl} ({len(members)} sites)") + print(f" ASE representatives : {rep_ase}") + print(f" SOAP representative : {rep_soap}") + print(f" Members:") + for idx, row in members.iterrows(): + conn = row["connectivity"] + formula = str(row["site_formula"]) + marker = "" + if idx in soap_idx: + marker = " <-- SOAP rep" + elif idx in ase_idx: + marker = " <-- ASE only" + print( + f" idx={idx:3d} conn={conn:2d} formula={formula:20s}" + f" ASE={row['in_ASE']} SOAP={row['in_SOAP']}{marker}" + ) + +print() +print("=" * 90) +print("SUMMARY BY SITE TYPE") +print("=" * 90) + +type_groups = df.groupby(["connectivity", "formula_str"]) +header = f"{'Type':30s} {'Total':>6s} {'ASE reps':>8s} {'SOAP reps':>9s} {'SOAP clusters':>14s}" +print(header) +print("-" * 70) +for (conn, formula), grp in type_groups: + n_total = len(grp) + n_ase = int(grp["in_ASE"].sum()) + n_soap = int(grp["in_SOAP"].sum()) + clusters = sorted(grp["soap_cluster"].unique()) + print(f"conn={conn} {formula:22s} {n_total:6d} {n_ase:8d} {n_soap:9d} {clusters}") + +print() +print(f"Total sites: {len(df)}") +print(f"ASE unique: {len(ase_idx)}") +print(f"SOAP unique: {len(soap_idx)}") +print(f"SOAP clusters: {len(df['soap_cluster'].unique())}") diff --git a/scripts/plot_sites.py b/scripts/plot_sites.py new file mode 100644 index 0000000..d78db35 --- /dev/null +++ b/scripts/plot_sites.py @@ -0,0 +1,168 @@ +"""Visualise adsorption sites on the Ni/Ru slab, colour-coded by SOAP cluster.""" +import copy +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.lines import Line2D +from matplotlib.patches import Circle +from ase.build import fcc111 +from autoadsorbate import Surface +from autoadsorbate.utils import _compute_soap_site_vectors, _soap_similarity_matrix +from scipy.cluster.hierarchy import fcluster, linkage +from scipy.spatial.distance import squareform + +# ── Build Ni/Ru slab ────────────────────────────────────────────── +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) +ase_idx = set(s_ase.site_df.index) +soap_idx = set(s_soap.site_df.index) + +# ── SOAP cluster labels ─────────────────────────────────────────── +vecs = _compute_soap_site_vectors(s.atoms, s.site_df) +sim = _soap_similarity_matrix(vecs) +dist = np.clip(1.0 - sim, 0.0, 2.0) +np.fill_diagonal(dist, 0.0) +Z = linkage(squareform(dist, checks=False), method="average") +labels = fcluster(Z, t=0.01, criterion="distance") + +coords = np.array(s.site_df["coordinates"].tolist()) +connectivity = np.array(s.site_df["connectivity"].tolist()) + +# ── Colour map: one colour per SOAP cluster ────────────────────── +n_clusters = len(set(labels)) +cmap = plt.colormaps.get_cmap("tab10") +cluster_ids = sorted(set(labels)) +colour_map = {cl: cmap(i / max(n_clusters - 1, 1)) for i, cl in enumerate(cluster_ids)} + +# Build cluster descriptions +cluster_desc = {} +for cl in cluster_ids: + mask = labels == cl + idx_in_cl = np.where(mask)[0] + row = s.site_df.iloc[idx_in_cl[0]] + conn = row["connectivity"] + formula = str(row["site_formula"]) + site_type = {1: "atop", 2: "bridge", 3: "hollow"}.get(conn, f"conn{conn}") + cluster_desc[cl] = f"{site_type} {formula}" + +# ── Helper: draw slab atoms as circles (top-down, x-y plane) ──── +atom_colours = {"Ni": "#8CBE6E", "Ru": "#D4534B"} +atom_radius = {"Ni": 0.9, "Ru": 0.95} + +def draw_slab(ax, atoms): + """Draw slab atoms as filled circles, top layer on top.""" + pos = atoms.positions + symbols = atoms.get_chemical_symbols() + # Sort by z so that top-layer atoms are drawn last (on top) + order = np.argsort(pos[:, 2]) + for i in order: + z_frac = (pos[i, 2] - pos[:, 2].min()) / (pos[:, 2].max() - pos[:, 2].min() + 1e-9) + alpha = 0.25 + 0.75 * z_frac # deeper atoms more transparent + r = atom_radius.get(symbols[i], 0.8) + c = atom_colours.get(symbols[i], "#AAAAAA") + circle = Circle( + (pos[i, 0], pos[i, 1]), r, + facecolor=c, edgecolor="black", linewidth=0.4, + alpha=alpha, zorder=2 + z_frac, + ) + ax.add_patch(circle) + + # Draw unit cell outline + cell = atoms.cell + corners = np.array([ + [0, 0], [cell[0, 0], cell[0, 1]], + [cell[0, 0] + cell[1, 0], cell[0, 1] + cell[1, 1]], + [cell[1, 0], cell[1, 1]], [0, 0], + ]) + ax.plot(corners[:, 0], corners[:, 1], "k-", lw=0.8, zorder=1) + +# ── Figure ──────────────────────────────────────────────────────── +fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + +for ax_i, (title, highlight_idx) in enumerate([ + (f"ASE representatives ({len(ase_idx)} sites)", ase_idx), + (f"SOAP representatives ({len(soap_idx)} sites)", soap_idx), +]): + ax = axes[ax_i] + + # Draw slab atoms in the same coordinate system + draw_slab(ax, s.atoms) + + # Plot ALL sites as small grey dots + ax.scatter( + coords[:, 0], coords[:, 1], + s=15, c="lightgrey", edgecolors="grey", linewidths=0.3, + zorder=5, label="_nolegend_", + ) + + # Overlay highlighted representative sites + for cl in cluster_ids: + mask = labels == cl + for site_i in np.where(mask)[0]: + df_idx = s.site_df.index[site_i] + if df_idx not in highlight_idx: + continue + conn = connectivity[site_i] + marker = {1: "o", 2: "s", 3: "^"}.get(conn, "D") + ax.scatter( + coords[site_i, 0], coords[site_i, 1], + s=120, c=[colour_map[labels[site_i]]], + edgecolors="black", linewidths=1.0, + marker=marker, zorder=10, + ) + # Label with index + ax.annotate( + str(df_idx), (coords[site_i, 0], coords[site_i, 1]), + textcoords="offset points", xytext=(5, 5), + fontsize=7, fontweight="bold", zorder=11, + ) + + ax.set_title(title, fontsize=12, fontweight="bold") + ax.set_xlabel("x (Å)") + ax.set_ylabel("y (Å)") + ax.set_aspect("equal") + ax.autoscale_view() + # Add a small margin + pad = 1.0 + ax.set_xlim(coords[:, 0].min() - pad, coords[:, 0].max() + pad) + ax.set_ylim(coords[:, 1].min() - pad, coords[:, 1].max() + pad) + +# ── Legend ───────────────────────────────────────────────────────── +legend_handles = [] +# Cluster colours +for cl in cluster_ids: + legend_handles.append( + Line2D([0], [0], marker="o", color="w", + markerfacecolor=colour_map[cl], markeredgecolor="black", + markersize=8, label=f"Cluster {cl}: {cluster_desc[cl]}") + ) +# Marker shapes for connectivity +legend_handles.append(Line2D([0], [0], marker="o", color="w", markerfacecolor="grey", + markersize=8, label="atop (conn=1)")) +legend_handles.append(Line2D([0], [0], marker="s", color="w", markerfacecolor="grey", + markersize=8, label="bridge (conn=2)")) +legend_handles.append(Line2D([0], [0], marker="^", color="w", markerfacecolor="grey", + markersize=8, label="hollow (conn=3)")) +# Slab atoms +legend_handles.append(Line2D([0], [0], marker="o", color="w", markerfacecolor="#8CBE6E", + markeredgecolor="black", markersize=10, label="Ni atom")) +legend_handles.append(Line2D([0], [0], marker="o", color="w", markerfacecolor="#D4534B", + markeredgecolor="black", markersize=10, label="Ru atom")) + +fig.legend(handles=legend_handles, loc="lower center", ncol=4, fontsize=8, + frameon=True, bbox_to_anchor=(0.5, -0.02)) + +fig.suptitle("Ni(111) 3×3 + 1 Ru dopant — Adsorption sites (top view)", fontsize=14, fontweight="bold") +fig.tight_layout(rect=[0, 0.10, 1, 0.95]) +fig.savefig("README_files/site_comparison.png", dpi=150, bbox_inches="tight") +plt.show() +print("Saved to README_files/site_comparison.png") diff --git a/tests/test_all.py b/tests/test_all.py index 19240ed..b2fb474 100644 --- a/tests/test_all.py +++ b/tests/test_all.py @@ -12,3 +12,205 @@ def test_Surface(): def test_Fragment(): f = Fragment(smile="COC", to_initialize=5) assert f.smile == "COC" + + +# ---------- SOAP-based symmetry reduction tests ---------- + +def _make_surface(): + """Helper: build a small Cu(111) Surface for SOAP tests.""" + from ase.build import fcc111 + slab = fcc111("Cu", (3, 3, 3), periodic=True, vacuum=10) + return Surface(slab) + + +def test_soap_similarity_matrix(): + """get_soap_similarity_matrix returns a square, symmetric matrix with 1s on the diagonal.""" + import numpy as np + s = _make_surface() + sim = s.get_soap_similarity_matrix() + n = len(s.site_df) + assert sim.shape == (n, n), f"Expected ({n},{n}), got {sim.shape}" + np.testing.assert_allclose(np.diag(sim), 1.0, atol=1e-6) + np.testing.assert_allclose(sim, sim.T, atol=1e-10) + + +def test_get_nonequivalent_sites_soap_hcluster(): + """get_nonequivalent_sites_soap (hcluster) returns fewer sites than the total.""" + s = _make_surface() + n_before = len(s.site_df) + unique = s.get_nonequivalent_sites_soap(similarity_threshold=0.99, method="hcluster") + assert 0 < len(unique) <= n_before + # all returned labels must exist in the original DataFrame + assert all(idx in s.site_df.index for idx in unique) + + +def test_get_nonequivalent_sites_soap_greedy(): + """get_nonequivalent_sites_soap (greedy) returns fewer sites than the total.""" + s = _make_surface() + n_before = len(s.site_df) + unique = s.get_nonequivalent_sites_soap(similarity_threshold=0.99, method="greedy") + assert 0 < len(unique) <= n_before + + +def test_soap_reduce(): + """soap_reduce mutates site_df in-place and reduces its length.""" + s = _make_surface() + n_before = len(s.site_df) + s.soap_reduce(similarity_threshold=0.99) + n_after = len(s.site_df) + assert 0 < n_after <= n_before + # site_dict should stay in sync + assert len(s.site_dict["coordinates"]) == n_after + + +def test_soap_reduce_vs_sym_reduce(): + """SOAP reduces sites on a perfect fcc(111) slab; count should be + <= ASE's result (SOAP at 0.9 may merge more aggressively).""" + import copy + s1 = _make_surface() + s2 = copy.deepcopy(s1) + s1.sym_reduce() + s2.soap_reduce(similarity_threshold=0.9) + assert 0 < len(s2.site_df) <= len(s1.site_df), ( + f"sym_reduce gave {len(s1.site_df)}, soap_reduce gave {len(s2.site_df)}" + ) + + +# ---------- Unified API tests (method='ase' | 'soap') ---------- + +def test_sym_reduce_method_ase(): + """sym_reduce(method='ase') behaves like the legacy sym_reduce().""" + import copy + s1 = _make_surface() + s2 = copy.deepcopy(s1) + s1.sym_reduce() # default → method='ase' + s2.sym_reduce(method="ase") + assert list(s1.site_df.index) == list(s2.site_df.index) + + +def test_sym_reduce_method_soap(): + """sym_reduce(method='soap') produces the same result as soap_reduce().""" + import copy + s1 = _make_surface() + s2 = copy.deepcopy(s1) + s1.sym_reduce(method="soap", similarity_threshold=0.99) + s2.soap_reduce(similarity_threshold=0.99) + assert list(s1.site_df.index) == list(s2.site_df.index) + + +def test_get_nonequivalent_sites_method_switch(): + """get_nonequivalent_sites dispatches correctly via the method param.""" + s = _make_surface() + ase_sites = s.get_nonequivalent_sites(method="ase") + soap_sites = s.get_nonequivalent_sites(method="soap", similarity_threshold=0.99) + # Both should reduce the site count; SOAP at 0.9 may merge more + assert 0 < len(soap_sites) <= len(ase_sites) + + +def test_sym_reduce_invalid_method(): + """sym_reduce raises ValueError for an unknown method.""" + import pytest + s = _make_surface() + with pytest.raises(ValueError, match="Unknown method"): + s.sym_reduce(method="invalid") + + +# ---------- Ni/Ru broken-symmetry comparison ---------- + +def _make_ni_ru_surface(): + """Build a Ni(111) slab and replace one surface Ni with Ru.""" + from ase.build import fcc111 + slab = fcc111("Ni", (3, 3, 3), periodic=True, vacuum=10) + # Replace one top-layer atom with Ru to break the perfect symmetry + top_z = slab.positions[:, 2].max() + for atom in slab: + if abs(atom.position[2] - top_z) < 0.1: + atom.symbol = "Ru" + break + return slab + + +def test_ni_ru_ase_vs_soap_site_counts(): + """On a Ni(111) slab with one Ni→Ru substitution, both methods should + find more unique site types than on the pristine surface and the SOAP + method should return a plausible number of groups.""" + import copy + + slab = _make_ni_ru_surface() + s_ase = Surface(slab) + s_soap = copy.deepcopy(s_ase) + + n_total = len(s_ase.site_df) + + s_ase.sym_reduce(method="ase") + s_soap.sym_reduce(method="soap", similarity_threshold=0.99) + + n_ase = len(s_ase.site_df) + n_soap = len(s_soap.site_df) + + # Both should reduce the site count + assert 0 < n_ase <= n_total + assert 0 < n_soap <= n_total + + # The broken-symmetry slab should have more unique sites than + # the pristine slab (pristine Ni(111) 3×3 has ~4 unique types) + assert n_ase > 1 + assert n_soap > 1 + + print( + f"\n Ni/Ru comparison: total={n_total} ASE={n_ase} SOAP={n_soap}" + ) + + +def test_ni_ru_timings(capsys): + """Benchmark ASE vs SOAP symmetry reduction on the Ni/Ru slab and print + a timing summary. This is informational – the test always passes.""" + import copy + import time + + slab = _make_ni_ru_surface() + s_ase = Surface(slab) + s_soap = copy.deepcopy(s_ase) + + n_total = len(s_ase.site_df) + + # --- ASE timing --- + t0 = time.perf_counter() + s_ase.sym_reduce(method="ase") + t_ase = time.perf_counter() - t0 + + # --- SOAP timing --- + t0 = time.perf_counter() + s_soap.sym_reduce(method="soap", similarity_threshold=0.99) + t_soap = time.perf_counter() - t0 + + n_ase = len(s_ase.site_df) + n_soap = len(s_soap.site_df) + + summary = ( + f"\n{'='*60}\n" + f" Ni(111) 3×3 slab with 1 Ni → Ru substitution\n" + f" Total sites before reduction : {n_total}\n" + f"{'─'*60}\n" + f" ASE SymmetryEquivalenceCheck : {n_ase:3d} unique sites " + f"({t_ase:.3f} s)\n" + f" SOAP + hierarchical cluster : {n_soap:3d} unique sites " + f"({t_soap:.3f} s)\n" + f"{'='*60}" + ) + print(summary) + + # Compare groupings + ase_indices = set(s_ase.site_df.index) + soap_indices = set(s_soap.site_df.index) + common = ase_indices & soap_indices + ase_only = ase_indices - soap_indices + soap_only = soap_indices - ase_indices + + print(f" Representatives in common : {len(common)}") + print(f" Only in ASE result : {len(ase_only)} {sorted(ase_only)}") + print(f" Only in SOAP result : {len(soap_only)} {sorted(soap_only)}") + print() + + # Informational – always pass + assert True