Skip to content

Commit 793ada9

Browse files
committed
Added kabsch algorithm to converter
The goal is to be able to match the rmsd of two atom mapped species xyzs. Good for TS, verification, atom mapping, etc.
1 parent a504aa7 commit 793ada9

File tree

1 file changed

+23
-0
lines changed

1 file changed

+23
-0
lines changed

arc/species/converter.py

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
from typing import TYPE_CHECKING, Dict, Iterable, List, Optional, Tuple, Union
99

1010
from ase import Atoms
11+
from scipy.spatial.transform import Rotation
1112
from openbabel import openbabel as ob
1213
from openbabel import pybel
1314
from rdkit import Chem
@@ -2446,3 +2447,25 @@ def sorted_distances_of_atom(xyz_dict: dict, atom_index: int) -> List[Tuple[int,
24462447

24472448
distances = [(i, d_matrix[atom_index, i]) for i in range(d_matrix.shape[0]) if i != atom_index]
24482449
return sorted(distances, key=lambda x: x[1])
2450+
2451+
2452+
def kabsch(xyz1: dict, xyz2: dict) -> float:
2453+
"""
2454+
Return the kabsch similarity score between two sets of Cartesian coordinates in Ångstrom.
2455+
The algorithm requires the atoms to be ordered the same in both sets of coordinates.
2456+
This will not be directly useful for comparing two conformers of the same species,
2457+
but could be used to compare two different species with the same atom types and counts. (e.g., isomers, reactants and products, etc.).
2458+
2459+
Args:
2460+
xyz1 (dict): The first set of Cartesian coordinates.
2461+
xyz2 (dict): The second set of Cartesian coordinates.
2462+
2463+
Returns:
2464+
float: The Kabsch similarity score in Ångstrom.
2465+
"""
2466+
if xyz1["symbols"] != xyz2["symbols"]:
2467+
raise ValueError("The two xyz coordinates must have the same atom symbols to compute Kabsch score.")
2468+
xyz1, xyz2 = translate_to_center_of_mass(xyz1), translate_to_center_of_mass(xyz2)
2469+
coords1, coords2 = np.array(xyz1['coords']), np.array(xyz2['coords'])
2470+
_, score = Rotation.align_vectors(coords1, coords2)
2471+
return score

0 commit comments

Comments
 (0)