Added implementation of the Kabsch algorithm in ARCSpecies API and as a standalone function in converter. (#819)

### About the algorithm
Kabsch’s algorithm computes the optimal rotation that superimposes two
sets of vectors (represented as matrices) by minimizing the
root-mean-square deviation (RMSD) between them.
For numerical stability and robustness, this implementation relies on
the optimized routine provided by SciPy via
[scipy.spatial.transform.Rotation.align_vectors](https://docs.scipy.org/doc/scipy-1.16.2/reference/generated/scipy.spatial.transform.Rotation.align_vectors.html)

### This PR
This PR introduces a Kabsch-based alignment utility in
`species/converter.py` and integrates it into the `ARCSpecies` API.
The implementation explicitly accounts for key assumptions required for
molecular alignment, including S2S mapping.
To support this workflow, a new helper method, `translate_to_indices`,
was added. This method converts an atom mapping into index-based
ordering, enabling consistent reordering and alignment of species
geometries prior to applying the Kabsch algorithm.

Overall, this addition enables robust, minimal-RMSD alignment of mapped
species geometries while preserving chemical correspondence.

Author: kfir4444

arc/species/converter.py

@@ -8,6 +8,7 @@
 from typing import TYPE_CHECKING, Dict, Iterable, List, Optional, Tuple, Union
 
 from ase import Atoms
+from scipy.spatial.transform import Rotation
 from openbabel import openbabel as ob
 from openbabel import pybel
 from rdkit import Chem
@@ -2446,3 +2447,25 @@ def sorted_distances_of_atom(xyz_dict: dict, atom_index: int) -> List[Tuple[int,
 
     distances = [(i, d_matrix[atom_index, i]) for i in range(d_matrix.shape[0]) if i != atom_index]
     return sorted(distances, key=lambda x: x[1])
+
+
+def kabsch(xyz1: dict, xyz2: dict) -> float:
+    """
+    Return the kabsch similarity score between two sets of Cartesian coordinates in Ångstrom.
+    The algorithm requires the atoms to be ordered the same in both sets of coordinates.
+    This will not be directly useful for comparing two conformers of the same species,
+    but could be used to compare two different species with the same atom types and counts. (e.g., isomers, reactants and products, etc.).
+
+    Args:
+        xyz1 (dict): The first set of Cartesian coordinates.
+        xyz2 (dict): The second set of Cartesian coordinates.
+
+    Returns:
+        float: The Kabsch similarity score in Ångstrom.
+    """
+    if xyz1["symbols"] != xyz2["symbols"]:
+        raise ValueError("The two xyz coordinates must have the same atom symbols to compute Kabsch score.")
+    xyz1, xyz2 = translate_to_center_of_mass(xyz1), translate_to_center_of_mass(xyz2)
+    coords1, coords2 = np.array(xyz1['coords']), np.array(xyz2['coords'])
+    _, score = Rotation.align_vectors(coords1, coords2)
+    return score

arc/species/converter_test.py

@@ -9,6 +9,7 @@
 import os
 
 import numpy as np
+from scipy.spatial.transform import Rotation
 import unittest
 
 from ase import Atoms
@@ -4966,6 +4967,70 @@ def test_sorted_distances_of_atom(self):
         with self.assertRaises(IndexError):
             converter.sorted_distances_of_atom(xyz_dict, 5)
 
+    def test_kabsch(self):
+        """Test the kabsch function"""
+        xyz1 = {'symbols': ('O', 'H', 'H'), 'isotopes': (16, 1, 1),
+                'coords': ((0.0, 0.0, 0.0),
+                           (0.0, 0.757, 0.586),
+                           (0.0, -0.757, 0.586))}
+        xyz2 = converter.translate_xyz(xyz1, (10.0, 0.0, 0.0))
+        score = converter.kabsch(xyz1, xyz2)
+        self.assertAlmostEqual(score, 0.0, places=5)
+
+        r = Rotation.from_quat([0, 0, np.sin(np.pi/4), np.cos(np.pi/4)])
+        self.assertAlmostEqual(converter.kabsch(xyz1, converter.xyz_from_data(coords=r.apply(np.array(xyz1["coords"])), symbols=xyz1["symbols"], isotopes=xyz1["isotopes"])), 0.0, places=5)
+        xyz2 = {i:v for i, v in xyz1.items()}
+        xyz2['symbols'] = ('O', 'H', 'H', 'H')
+        with self.assertRaises(ValueError):
+            converter.kabsch(xyz1, xyz2)
+
+        # Wildtype test: Aspirin (21 atoms)
+        aspirin_xyz = {
+            'symbols': ('O', 'C', 'O', 'C', 'C', 'C', 'C', 'C', 'C', 'H', 'H', 'H', 'H', 'C', 'O', 'O', 'C', 'H', 'H', 'H', 'H'),
+            'isotopes': (16, 12, 16, 12, 12, 12, 12, 12, 12, 1, 1, 1, 1, 12, 16, 16, 12, 1, 1, 1, 1),
+            'coords': ((-2.6373, 1.2580, -0.2078),
+                       (-1.7770, 0.3860, -0.0632),
+                       (-2.1558, -0.8741, 0.1691),
+                       (-0.3592, 0.7788, -0.1627),
+                       (0.6120, -0.2186, -0.0163),
+                       (1.9423, 0.1873, -0.1118),
+                       (2.2874, 1.5173, -0.3475),
+                       (1.3090, 2.4837, -0.4907),
+                       (-0.0249, 2.1150, -0.3989),
+                       (2.6841, -0.5841, 0.0004),
+                       (3.3216, 1.8105, -0.4206),
+                       (1.5971, 3.5230, -0.6740),
+                       (-0.7675, 2.8711, -0.5103),
+                       (0.1843, -1.6369, 0.2443),
+                       (0.8550, -2.5186, 0.7303),
+                       (-1.1095, -1.8841, -0.1121),
+                       (-1.6441, -3.2084, 0.0818),
+                       (-2.6718, -3.1782, -0.2678),
+                       (-1.0850, -3.8967, -0.5484),
+                       (-1.6041, -3.4832, 1.1345),
+                       (-3.5658, 0.9859, -0.1477))
+        }
+
+        # Case 1: Pure rotation (should be 0.0)
+        # Rotate 90 degrees about the z-axis, then 45 degrees about the y-axis
+        r = Rotation.from_euler('zy', [90, 45], degrees=True)
+        rotated_coords = r.apply(np.array(aspirin_xyz["coords"]))
+        aspirin_rotated_xyz = converter.xyz_from_data(coords=rotated_coords,
+                                                      symbols=aspirin_xyz["symbols"],
+                                                      isotopes=aspirin_xyz["isotopes"])
+        score = converter.kabsch(aspirin_xyz, aspirin_rotated_xyz)
+        self.assertAlmostEqual(score, 0.0, places=4)
+
+        # Case 2: Random structural perturbation (should not be 0.0)
+        # Add random noise to coordinates
+        rng = np.random.RandomState(42)
+        perturbed_coords = np.array(aspirin_xyz["coords"]) + 0.1 * rng.rand(*np.array(aspirin_xyz["coords"]).shape)
+        aspirin_perturbed_xyz = converter.xyz_from_data(coords=perturbed_coords,
+                                                        symbols=aspirin_xyz["symbols"],
+                                                        isotopes=aspirin_xyz["isotopes"])
+        score = converter.kabsch(aspirin_xyz, aspirin_perturbed_xyz)
+        self.assertGreater(score, 0.01)
+
     @classmethod
     def tearDownClass(cls):
         """

arc/species/species.py

@@ -46,10 +46,12 @@
                                    order_atoms_in_mol_list,
                                    remove_dummies,
                                    rmg_mol_from_inchi,
+                                   sort_xyz_using_indices,
                                    str_to_xyz,
                                    translate_to_center_of_mass,
                                    xyz_from_data,
                                    xyz_to_str,
+                                   kabsch,
                                    )
 from arc.species.perceive import perceive_molecule_from_xyz, is_mol_valid
 from arc.species.vectors import calculate_angle, calculate_distance, calculate_dihedral_angle
@@ -2142,6 +2144,27 @@ def get_bonds(self) -> List[tuple]:
             for atom2, bond12 in atom1.edges.items():
                 bonds.append(tuple(sorted((self.mol.atoms.index(atom1), self.mol.atoms.index(atom2)))))
         return list(set(bonds))
+    
+    def kabsch(self, other: 'ARCSpecies', map_: list) -> float:
+        """
+        Calculate the Kabsch RMSD between this species and another species.
+
+        Args:
+            other (ARCSpecies): The other species to compare to.
+            map_ (list): A list of atom indices mapping atoms from this species to the other species. (i.e., if
+                         this species has atoms [A, B, C] and the other species has atoms [C, A, B], then map_ would be [1, 2, 0]
+        Returns:
+            float: The Kabsch RMSD value.
+        """
+        if not isinstance(other, ARCSpecies):
+            raise SpeciesError(f'Other must be an ARCSpecies instance, got {type(other)}.\n'
+                               f'If you meant to use the XYZ coordinates directly, use arc.species.converter.kabsch.')
+
+        if len(map_) != self.number_of_atoms:
+            raise SpeciesError(f'The map_ list must have the same length as the number of atoms in {self.label} '
+                               f'({self.number_of_atoms}), got {len(map_)}.')
+
+        return kabsch(self.get_xyz(), sort_xyz_using_indices(other.get_xyz(), map_))
 
 
 class TSGuess(object):

arc/species/species_test.py

@@ -2804,6 +2804,45 @@ def test_species_indexing(self):
         self.assertEqual(spc_b.index, 1)
         self.assertEqual(spc_c.index, 2)
 
+    def test_kabsch(self):
+        """Test the kabsch() method."""
+        # Test with self (RMSD should be 0)
+        rmsd = self.spc1.kabsch(self.spc1, [0, 1, 2, 3, 4, 5])
+        self.assertAlmostEqual(rmsd, 0.0, delta=1e-5)
+
+        # Test with a copy (RMSD should be 0)
+        spc1_copy = self.spc1.copy()
+        rmsd = self.spc1.kabsch(spc1_copy, [0, 1, 2, 3, 4, 5])
+        self.assertAlmostEqual(rmsd, 0.0, delta=1e-5)
+
+        # Test with shuffled atoms
+        # Create a shuffled version of spc1: swap first two C atoms
+        # spc1: C, C, O, H, H, H
+        # shuffled: C(1), C(0), O(2), H(3), H(4), H(5)
+        # xyz of spc1
+        xyz = self.spc1.get_xyz()
+        new_coords = (xyz['coords'][1], xyz['coords'][0]) + xyz['coords'][2:]
+        new_symbols = (xyz['symbols'][1], xyz['symbols'][0]) + xyz['symbols'][2:]
+        new_isotopes = (xyz['isotopes'][1], xyz['isotopes'][0]) + xyz['isotopes'][2:]
+        shuffled_xyz = {'symbols': new_symbols, 'isotopes': new_isotopes, 'coords': new_coords}
+        spc_shuffled = ARCSpecies(label='shuffled', xyz=shuffled_xyz, smiles='C=C[O]')
+        
+        # Map: we want to pull atoms from spc_shuffled to match spc1
+        # spc1[0] is C(0). In spc_shuffled, C(0) is at index 1.
+        # spc1[1] is C(1). In spc_shuffled, C(1) is at index 0.
+        # Rest are same.
+        map_indices = [1, 0, 2, 3, 4, 5]
+        rmsd = self.spc1.kabsch(spc_shuffled, map_indices)
+        self.assertAlmostEqual(rmsd, 0.0, delta=1e-5)
+
+        # Test exception
+        with self.assertRaises(SpeciesError):
+            self.spc1.kabsch("not a species", [])
+
+        # Test incorrect map_ length
+        with self.assertRaises(SpeciesError):
+            self.spc1.kabsch(self.spc1, [0, 1, 2])
+
 
 class TestTSGuess(unittest.TestCase):
     """