From b588f4870b04faa0b1af0a5e406579af92bbdebb Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Wed, 17 Sep 2025 14:11:21 +0300 Subject: [PATCH 1/6] Added ic_params classes These classes unify the concept of internal coordinates and their parameters for adding to the zmat. --- arc/species/__init__.py | 1 + arc/species/ic_params.py | 133 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 134 insertions(+) create mode 100644 arc/species/ic_params.py diff --git a/arc/species/__init__.py b/arc/species/__init__.py index 9baf4ba33e..1c0302bb8d 100644 --- a/arc/species/__init__.py +++ b/arc/species/__init__.py @@ -2,4 +2,5 @@ import arc.species.converter import arc.species.perceive import arc.species.species +import arc.species.ic_params from arc.species.species import ARCSpecies, TSGuess diff --git a/arc/species/ic_params.py b/arc/species/ic_params.py new file mode 100644 index 0000000000..a3ff09e6ef --- /dev/null +++ b/arc/species/ic_params.py @@ -0,0 +1,133 @@ +from dataclasses import dataclass +from enum import Enum +from typing import List, Optional + + +class ParamKey(str, Enum): + """Internal-coordinate parameter type.""" + R = "R" # bond length + A = "A" # bond angle + D = "D" # dihedral + + +@dataclass +class AnchorSpec: + """ + Indices that define this parameter relative to the two molecules: + - m2i: indices in xyz2 (local to the second structure being added) + - m1i: indices in xyz1 (anchors in the first structure) + """ + m2i: List[int] + m1i: List[int] + + def ensure_primary(self, required_first_m2i: int) -> None: + """ + Enforce that the first entry in m2i is `required_first_m2i`. + If absent, insert it; if present but not first, move it to front. + """ + if not isinstance(self.m2i, list): + self.m2i = [required_first_m2i] + return + self.m2i = [int(x) for x in self.m2i] + if self.m2i and self.m2i[0] == required_first_m2i: + return + if required_first_m2i in self.m2i: + self.m2i.remove(required_first_m2i) + self.m2i.insert(0, required_first_m2i) + + +@dataclass +class PlacementParam: + """ + One IC parameter (R/A/D) for a single atom placement. + + Attributes + ---------- + key : ParamKey + Type of parameter (R/A/D). + val : float + Numeric value of the parameter (Å, degrees). + anchors : AnchorSpec + Anchor indices in xyz2 (m2i) and xyz1 (m1i). + index : Optional[int] + Local atom index in xyz2 this parameter primarily refers to. If not + given, it is inferred as anchors.m2i[0] after normalization. + final_index : Optional[int] + The shifted index in the combined structure (n + index), set by + calling `set_base_n(n)` where n = len(xyz1). + """ + key: ParamKey + val: float + anchors: AnchorSpec + index: Optional[int] = None + final_index: Optional[int] = None + + def ensure_primary_m2i(self, required_first_m2i: int) -> None: + self.anchors.ensure_primary(required_first_m2i) + if self.index is None: + self.index = self.anchors.m2i[0] + + def set_base_n(self, base_n: int) -> None: + if self.index is not None: + self.final_index = base_n + int(self.index) + + def build_key(self, base_n: int) -> str: + """ + Build zmat variable name like 'A__...' + m2 indices are shifted by base_n; m1 indices are unchanged. + """ + m2_part = [str(int(i) + base_n) for i in self.anchors.m2i] + m1_part = [str(int(i)) for i in self.anchors.m1i] + return "_".join([self.key.value] + m2_part + m1_part) + + +@dataclass +class Atom1Params: + """Placement parameters for the FIRST atom of xyz2 (must target m2i[0] == 0).""" + R: PlacementParam + A: PlacementParam + D: PlacementParam + + def __post_init__(self) -> None: + self.R.ensure_primary_m2i(0) + self.A.ensure_primary_m2i(0) + self.D.ensure_primary_m2i(0) + # Key sanity: + if self.R.key is not ParamKey.R or self.A.key is not ParamKey.A or self.D.key is not ParamKey.D: + raise ValueError("Atom1Params expects keys R, A, D respectively.") + + def set_base_n(self, base_n: int) -> None: + self.R.set_base_n(base_n) + self.A.set_base_n(base_n) + self.D.set_base_n(base_n) + + +@dataclass +class Atom2Params: + """Placement parameters for the SECOND atom of xyz2 (must target m2i[0] == 1).""" + A: PlacementParam + D: PlacementParam + + def __post_init__(self) -> None: + self.A.ensure_primary_m2i(1) + self.D.ensure_primary_m2i(1) + if self.A.key is not ParamKey.A or self.D.key is not ParamKey.D: + raise ValueError("Atom2Params expects keys A, D respectively.") + + def set_base_n(self, base_n: int) -> None: + self.A.set_base_n(base_n) + self.D.set_base_n(base_n) + + +@dataclass +class Atom3Params: + """Placement parameters for the THIRD atom of xyz2 (must target m2i[0] == 2).""" + D: PlacementParam + + def __post_init__(self) -> None: + self.D.ensure_primary_m2i(2) + if self.D.key is not ParamKey.D: + raise ValueError("Atom3Params expects key D.") + + def set_base_n(self, base_n: int) -> None: + self.D.set_base_n(base_n) From f814f1b9ea89fd5d2f566f5f6cb581f25e8cd031 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Wed, 17 Sep 2025 14:11:30 +0300 Subject: [PATCH 2/6] test: ic_params --- arc/species/ic_params_test.py | 173 ++++++++++++++++++++++++++++++++++ 1 file changed, 173 insertions(+) create mode 100644 arc/species/ic_params_test.py diff --git a/arc/species/ic_params_test.py b/arc/species/ic_params_test.py new file mode 100644 index 0000000000..b716f3f6cb --- /dev/null +++ b/arc/species/ic_params_test.py @@ -0,0 +1,173 @@ +import unittest + +from arc.common import ARC_PATH +from arc.species import ic_params + +class TestPerceive(unittest.TestCase): + @classmethod + def setUpClass(cls): + pass + + def test_AnchorSpec_ensure_primary_when_not_list(self): + """test the method AnchorSpec.ensure_primary""" + a = ic_params.AnchorSpec(m2i="not-a-list", m1i=[7, 8]) + a.ensure_primary(required_first_m2i=0) + self.assertEqual(a.m2i, [0]) + + def test_AnchorSpec_ensure_primary_move_to_front(self): + """test the method AnchorSpec.ensure_primary""" + a = ic_params.AnchorSpec(m2i=[2, 0, 3], m1i=[7, 8]) + a.ensure_primary(0) + self.assertEqual(a.m2i, [0, 2, 3]) + + def test_AnchorSpec_ensure_primary_noop_if_already_first(self): + """test the method AnchorSpec.ensure_primary""" + a = ic_params.AnchorSpec(m2i=[1, 4, 5], m1i=[7, 8]) + a.ensure_primary(1) + self.assertEqual(a.m2i, [1, 4, 5]) + + def test_AnchorSpec_ensure_primary_casts_to_int(self): + """test the method AnchorSpec.ensure_primary""" + a = ic_params.AnchorSpec(m2i=["1", "3"], m1i=[7, 8]) + a.ensure_primary(1) + self.assertEqual(a.m2i, [1, 3]) + + def test_PlacementParam_ensure_primary_m2i_sets_index_if_none(self): + """test the method ensure_primary_m2i""" + p = ic_params.PlacementParam( + key=ic_params.ParamKey.R, + val=1.23, + anchors=ic_params.AnchorSpec(m2i=[3, 0, 2], m1i=[9]), + index=None, + ) + p.ensure_primary_m2i(0) + self.assertEqual(p.anchors.m2i[0], 0) + self.assertEqual(p.index, 0) + + def test_PlacementParam_set_base_n_sets_final_index(self): + """test the method set_base_n""" + p = ic_params.PlacementParam( + key=ic_params.ParamKey.A, + val=109.5, + anchors=ic_params.AnchorSpec(m2i=[1, 0], m1i=[5]), + index=1, + ) + p.set_base_n(10) + self.assertEqual(p.final_index, 11) + + def test_PlacementParam_set_base_n_noop_if_index_none(self): + """test the method set_base_n""" + p = ic_params.PlacementParam( + key=ic_params.ParamKey.D, + val=180.0, + anchors=ic_params.AnchorSpec(m2i=[2, 1, 0], m1i=[4, 5]), + index=None, + ) + p.set_base_n(7) + self.assertIsNone(p.final_index) + + def test_PlacementParam_build_key_shifts_only_m2(self): + """test the method build_key""" + p = ic_params.PlacementParam( + key=ic_params.ParamKey.A, + val=109.5, + anchors=ic_params.AnchorSpec(m2i=[0, 1, 2], m1i=[5, 6]), + index=0, + ) + key = p.build_key(base_n=10) + self.assertEqual(key, "A_10_11_12_5_6") + + def test_Atom1Params_post_init_enforces_m2i0_and_keys(self): + """test the class Atom1Params """ + R = ic_params.PlacementParam(ic_params.ParamKey.R, 1.0, ic_params.AnchorSpec([3, 0], [100]), index=None) + A = ic_params.PlacementParam(ic_params.ParamKey.A, 90.0, ic_params.AnchorSpec([2, 0], [101, 102]), index=None) + D = ic_params.PlacementParam(ic_params.ParamKey.D, 180.0, ic_params.AnchorSpec([5, 0, 4], [103, 104, 105]), index=None) + atom1 = ic_params.Atom1Params(R=R, A=A, D=D) + self.assertEqual(atom1.R.anchors.m2i[0], 0) + self.assertEqual(atom1.A.anchors.m2i[0], 0) + self.assertEqual(atom1.D.anchors.m2i[0], 0) + self.assertEqual(atom1.R.index, 0) + self.assertEqual(atom1.A.index, 0) + self.assertEqual(atom1.D.index, 0) + + def test_Atom1Params_set_base_n_propagates_final_index(self): + """test the method set_base_n for Atom1Params""" + atom1 = ic_params.Atom1Params( + R=ic_params.PlacementParam(ic_params.ParamKey.R, 1.1, ic_params.AnchorSpec([0], [1]), index=None), + A=ic_params.PlacementParam(ic_params.ParamKey.A, 109.5, ic_params.AnchorSpec([0, 1], [2]), index=None), + D=ic_params.PlacementParam(ic_params.ParamKey.D, 180.0, ic_params.AnchorSpec([0, 1, 2], [3]), index=None), + ) + atom1.set_base_n(7) + self.assertEqual(atom1.R.final_index, 7) + self.assertEqual(atom1.A.final_index, 7) + self.assertEqual(atom1.D.final_index, 7) + + def test_Atom1Params_invalid_keys_raise(self): + """test the method Atom1Params raises ValueError on invalid keys""" + with self.assertRaises(ValueError): + ic_params.Atom1Params( + R=ic_params.PlacementParam(ic_params.ParamKey.A, 1.0, ic_params.AnchorSpec([0], [1])), + A=ic_params.PlacementParam(ic_params.ParamKey.R, 90.0, ic_params.AnchorSpec([0, 1], [2])), + D=ic_params.PlacementParam(ic_params.ParamKey.D, 180.0, ic_params.AnchorSpec([0, 1, 2], [3])), + ) + + def test_Atom2Params_post_init_enforces_m2i1_and_keys(self): + """test that the method enforces m2i[0] == 1 and correct keys""" + A = ic_params.PlacementParam(ic_params.ParamKey.A, 120.0, ic_params.AnchorSpec([4, 1, 0], [9, 10]), index=None) + D = ic_params.PlacementParam(ic_params.ParamKey.D, 60.0, ic_params.AnchorSpec([2, 1, 0], [11, 12, 13]), index=None) + atom2 = ic_params.Atom2Params(A=A, D=D) + self.assertEqual(atom2.A.anchors.m2i[0], 1) + self.assertEqual(atom2.D.anchors.m2i[0], 1) + self.assertEqual(atom2.A.index, 1) + self.assertEqual(atom2.D.index, 1) + + def test_Atom2Params_set_base_n_propagates_final_index(self): + """test the method Atom2Params.set_base_n""" + atom2 = ic_params.Atom2Params( + A=ic_params.PlacementParam(ic_params.ParamKey.A, 100.0, ic_params.AnchorSpec([1, 0], [2]), index=None), + D=ic_params.PlacementParam(ic_params.ParamKey.D, 180.0, ic_params.AnchorSpec([1, 0, 2], [3]), index=None), + ) + atom2.set_base_n(20) + self.assertEqual(atom2.A.final_index, 21) + self.assertEqual(atom2.D.final_index, 21) + + def test_Atom2Params_invalid_keys_raise(self): + """test that the method Atom2Params raises ValueError on invalid keys""" + with self.assertRaises(ValueError): + ic_params.Atom2Params( + A=ic_params.PlacementParam(ic_params.ParamKey.R, 1.0, ic_params.AnchorSpec([1, 0], [2])), + D=ic_params.PlacementParam(ic_params.ParamKey.D, 180.0, ic_params.AnchorSpec([1, 0, 2], [3])), + ) + + def test_Atom3Params_post_init_enforces_m2i2_and_key(self): + """test that the class Atom3Params enforces m2i[0] == 2 and correct key""" + D = ic_params.PlacementParam(ic_params.ParamKey.D, 179.9, ic_params.AnchorSpec([5, 2, 1, 0], [7, 8, 9]), index=None) + atom3 = ic_params.Atom3Params(D=D) + self.assertEqual(atom3.D.anchors.m2i[0], 2) + self.assertEqual(atom3.D.index, 2) + + def test_Atom3Params_set_base_n_propagates_final_index(self): + """test the method Atom3Params.set_base_n""" + atom3 = ic_params.Atom3Params( + D=ic_params.PlacementParam(ic_params.ParamKey.D, 120.0, ic_params.AnchorSpec([2, 1, 0], [4, 5, 6]), index=None) + ) + atom3.set_base_n(30) + self.assertEqual(atom3.D.final_index, 32) + + def test_Atom3Params_invalid_key_raises(self): + """test the class Atom3Params""" + with self.assertRaises(ValueError): + ic_params.Atom3Params(D=ic_params.PlacementParam(ic_params.ParamKey.A, 90.0, ic_params.AnchorSpec([2, 1, 0], [4, 5]))) + + def test_build_key_with_string_indices_and_shift(self): + """test the method build_key""" + p = ic_params.PlacementParam( + key=ic_params.ParamKey.R, + val=1.0, + anchors=ic_params.AnchorSpec(m2i=["2", "0"], m1i=["10", "11"]), + index=2, + ) + self.assertEqual(p.build_key(5), "R_7_5_10_11") + +if __name__ == '__main__': + unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) \ No newline at end of file From 7d171b3869882c19e5fd7bd8edd590c49f51c9ba Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Wed, 17 Sep 2025 14:12:19 +0300 Subject: [PATCH 3/6] Added add_two_xyzs logic --- arc/species/converter.py | 131 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 129 insertions(+), 2 deletions(-) diff --git a/arc/species/converter.py b/arc/species/converter.py index b8ec5e7d06..4253fe799f 100644 --- a/arc/species/converter.py +++ b/arc/species/converter.py @@ -5,7 +5,9 @@ import math import numpy as np import os -from typing import TYPE_CHECKING, Dict, Iterable, List, Optional, Tuple, Union +from typing import TYPE_CHECKING, Any, Dict, Iterable, List, Optional, Tuple, Union +from copy import deepcopy +from arc.species.ic_params import Atom1Params, Atom2Params, Atom3Params, ParamKey from ase import Atoms from openbabel import openbabel as ob @@ -35,7 +37,8 @@ get_atom_indices_from_zmat_parameter, get_parameter_from_atom_indices, zmat_to_coords, - xyz_to_zmat) + xyz_to_zmat, + check_zmat_vs_coords,) if TYPE_CHECKING: from arc.species.species import ARCSpecies @@ -2353,6 +2356,7 @@ def generate_initial_guess_r_a(atom_r_coord: tuple, X_position = np.array(atom_r_coord) + r_value * X_direction return X_position + def generate_midpoint_initial_guess(atom_r_coord, r_value, atom_a_coord, atom_b_coord, a_value): """Generate an initial guess midway between the two reference atoms.""" midpoint = (np.array(atom_a_coord) + np.array(atom_b_coord)) / 2.0 @@ -2360,6 +2364,7 @@ def generate_midpoint_initial_guess(atom_r_coord, r_value, atom_a_coord, atom_b_ direction /= np.linalg.norm(direction) return midpoint + r_value * direction + def generate_perpendicular_initial_guess(atom_r_coord, r_value, atom_a_coord, atom_b_coord, a_value): """Generate an initial guess that is perpendicular to the plane defined by the reference atoms.""" BA = np.array(atom_a_coord) - np.array(atom_b_coord) @@ -2372,18 +2377,140 @@ def generate_perpendicular_initial_guess(atom_r_coord, r_value, atom_a_coord, at perpendicular_vector /= np.linalg.norm(perpendicular_vector) return np.array(atom_r_coord) + r_value * perpendicular_vector + def generate_random_initial_guess(atom_r_coord, r_value, atom_a_coord, atom_b_coord, a_value): perturbation = np.random.uniform(-0.1, 0.1, 3) base_guess = generate_initial_guess_r_a(atom_r_coord, r_value, atom_a_coord, atom_b_coord, a_value) return base_guess + perturbation + def generate_shifted_initial_guess(atom_r_coord, r_value, atom_a_coord, atom_b_coord, a_value): shift = np.array([0.1, -0.1, 0.1]) # A deterministic shift base_guess = generate_initial_guess_r_a(atom_r_coord, r_value, atom_a_coord, atom_b_coord, a_value) return base_guess + shift + def generate_bond_length_initial_guess(atom_r_coord, r_value, atom_a_coord, atom_b_coord, a_value): """Generate an initial guess considering only the bond length to the reference atom.""" direction = np.random.uniform(-1.0, 1.0, 3) # Random direction direction /= np.linalg.norm(direction) # Normalize to unit vector return np.array(atom_r_coord) + r_value * direction + + +def update_zmat(zmat: dict, + element: str, + coords: tuple, + vars: dict, + n: int, + ) -> dict: + """ + Update the coordinates and internal coordinates of a zmat dictionary. + Args: + zmat (dict): The zmat dictionary to update. + coords (tuple): The new Cartesian coordinates. + vars (dict): The new internal coordinates. + n (int): The index of the new atom in the zmat. + Returns: + dict: The updated zmat dictionary. + """ + zmat = deepcopy(zmat) + zmat["symbols"] += (element,) + zmat["coords"] += (coords,) + zmat["vars"].update(vars) + zmat["map"].update({n:n}) + return zmat + + +def add_two_xyzs( + xyz1: Dict[str, Any], + xyz2: Dict[str, Any], + atom1: Atom1Params, + atom2: Atom2Params, + atom3: Atom3Params, + return_zmat: bool = False, +) -> Optional[Union[Dict[str, Any], Tuple[Dict[str, Any], Dict[str, Any]]]]: + """ + Combine two xyz dictionaries into one, based on internal-coordinate parameters + (objects) for the first three atoms in xyz2. The remaining atoms in xyz2 are + added based on their internal coordinates relative to the first three. + + Parameters + ---------- + xyz1, xyz2 : dict + Your existing xyz dicts. + atom1 : Atom1Params + R/A/D for placing xyz2's first atom relative to anchors in xyz1. + atom2 : Atom2Params + A/D for placing xyz2's second atom (its R is taken from xyz2's own zmat). + atom3 : Atom3Params + D for placing xyz2's third atom (its R/A are taken from xyz2's own zmat). + return_zmat : bool + If True, return (zmat, combined_xyz); otherwise return combined xyz dict. + + Returns + ------- + dict or (dict, dict) or None + Combined xyz, or (zmat, xyz) if requested. Returns None on failure. + """ + zmat1 = xyz_to_zmat(xyz1, consolidate=False) + zmat2 = xyz_to_zmat(xyz2, consolidate=False) + + if zmat1 is None or zmat2 is None: + logger.error('Cannot add two xyzs if one of them cannot be converted to a zmat.') + return None + + n = len(zmat1['symbols']) + z2_vars: Dict[str, float] = zmat2['vars'] + + atom1.set_base_n(n) + atom2.set_base_n(n) + atom3.set_base_n(n) + + shift_piece = lambda piece: str(int(piece) + n) if piece.isnumeric() else piece + shift_token = lambda token: '_'.join(shift_piece(p) for p in token.split('_')) + + for idx, (symbol, coords) in enumerate(zip(zmat2['symbols'], zmat2['coords'])): + # coords is a tuple like (), (R), (R, A), or (R, A, D) from z2 + new_coords = tuple(shift_token(c) for c in coords if c) + vmap: Dict[str, float] = {} + + if len(new_coords) == 0: + # First atom of xyz2 relative to anchors in xyz1 + R_key = atom1.R.build_key(n) + A_key = atom1.A.build_key(n) + D_key = atom1.D.build_key(n) + coords_out = (R_key, A_key, D_key) + vmap[R_key] = atom1.R.val + vmap[A_key] = atom1.A.val + vmap[D_key] = atom1.D.val + + elif len(new_coords) == 1: + # Second atom of xyz2: use its own R from z2, A/D from anchors in xyz1 + A_key = atom2.A.build_key(n) + D_key = atom2.D.build_key(n) + coords_out = (new_coords[0], A_key, D_key) + vmap[new_coords[0]] = z2_vars[coords[0]] # carry R from z2 + vmap[A_key] = atom2.A.val + vmap[D_key] = atom2.D.val + + elif len(new_coords) == 2: + # Third atom of xyz2: use R/A from z2, D from anchors in xyz1 + D_key = atom3.D.build_key(n) + coords_out = (new_coords[0], new_coords[1], D_key) + vmap[new_coords[0]] = z2_vars[coords[0]] # carry R + vmap[new_coords[1]] = z2_vars[coords[1]] # carry A + vmap[D_key] = atom3.D.val + + else: + # Remaining atoms: just shift tokens and copy the values from z2 + coords_out = new_coords + vmap.update({nc: z2_vars[oc] for nc, oc in zip(new_coords, coords)}) + + zmat1 = update_zmat(zmat1, symbol, coords_out, vmap, n + idx) + if zmat1 is None or xyz1 is None: + raise ValueError("The combined zmat or xyz is invalid. Please check the input parameters and xyzs.") + if not check_zmat_vs_coords(zmat1, xyz_to_np_array(xyz1)): + raise ValueError("The combined zmat is invalid. Please check the input parameters and xyzs.") + if return_zmat: + return zmat1, zmat_to_xyz(zmat1) + return zmat_to_xyz(zmat1) From 8c9fc48d5e9321aef45b452219bfd1adaf90eec3 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Wed, 17 Sep 2025 14:12:41 +0300 Subject: [PATCH 4/6] Test: updating zmat and add combining two xyz --- arc/species/converter_test.py | 109 ++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) diff --git a/arc/species/converter_test.py b/arc/species/converter_test.py index 8203d2f56c..2612fd06d8 100644 --- a/arc/species/converter_test.py +++ b/arc/species/converter_test.py @@ -4751,6 +4751,115 @@ def test_add_atom_to_xyz_using_internal_coords(self): self.assertAlmostEqual(calculate_param(coords=new_xyz_4['coords'], atoms=[3, 1, 14]), 77.4, places=1) self.assertAlmostEqual(calculate_param(coords=new_xyz_4['coords'], atoms=[2, 0, 1, 14]), 140, places=1) + def test_update_zmat(self): + """Test the update_zmat() function.""" + zmat1 = {'coords': ((None, None, None), + ('R_1_0', None, None), + ('R_2|3|4|5|6|7_0|0|0|1|1|1', + 'A_2|3|4|5|6|7_0|0|0|1|1|1_1|1|1|0|0|0', + None), + ('R_2|3|4|5|6|7_0|0|0|1|1|1', + 'A_2|3|4|5|6|7_0|0|0|1|1|1_1|1|1|0|0|0', + 'D_3|4|6|7_0|0|1|1_1|1|0|0_2|3|5|6'), + ('R_2|3|4|5|6|7_0|0|0|1|1|1', + 'A_2|3|4|5|6|7_0|0|0|1|1|1_1|1|1|0|0|0', + 'D_3|4|6|7_0|0|1|1_1|1|0|0_2|3|5|6'), + ('R_2|3|4|5|6|7_0|0|0|1|1|1', + 'A_2|3|4|5|6|7_0|0|0|1|1|1_1|1|1|0|0|0', + 'D_5_1_0_4'), + ('R_2|3|4|5|6|7_0|0|0|1|1|1', + 'A_2|3|4|5|6|7_0|0|0|1|1|1_1|1|1|0|0|0', + 'D_3|4|6|7_0|0|1|1_1|1|0|0_2|3|5|6'), + ('R_2|3|4|5|6|7_0|0|0|1|1|1', + 'A_2|3|4|5|6|7_0|0|0|1|1|1_1|1|1|0|0|0', + 'D_3|4|6|7_0|0|1|1_1|1|0|0_2|3|5|6')), + 'map': {0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6, 7: 7}, + 'symbols': ('C', 'C', 'H', 'H', 'H', 'H', 'H', 'H'), + 'vars': {'A_2|3|4|5|6|7_0|0|0|1|1|1_1|1|1|0|0|0': 110.56796391805013, + 'D_3|4|6|7_0|0|1|1_1|1|0|0_2|3|5|6': 240.00000225241752, + 'D_5_1_0_4': 59.99998551792418, + 'R_1_0': 1.512057900428772, + 'R_2|3|4|5|6|7_0|0|0|1|1|1': 1.0940840641657512}} + add_atom = "C" + coords1 = ("R_8_0", "A_8_1_0", "D_8_2_1_0",) + vars = {"R_8_0": 5.0, "A_8_1_0": 109.5, "D_8_2_1_0": 60.0} + n = len(zmat1['symbols']) + new_zmat1 = converter.update_zmat(zmat=zmat1, + element=add_atom, + coords=coords1, + vars=vars, + n=n, + ) + self.assertIsNot(zmat1, new_zmat1) # Test that a copy was made + self.assertEqual(len(new_zmat1['symbols']), len(zmat1['symbols']) + 1) # exactly one atom added + self.assertEqual(new_zmat1['symbols'][-1], add_atom) # the right atom was added at the end + self.assertEqual(new_zmat1['coords'][-1], coords1) + self.assertNotIn(vars.keys(), zmat1['vars'].keys()) + self.assertIn(vars.keys(), new_zmat1['vars'].keys()) + + zmat2 = {'coords': ((None, None, None),), 'map': {0: 0}, 'symbols': ('H',), 'vars': {}} # lone H atom + add_atom = "H" + coords2 = ("R_1_0", None, None,) + vars = {"R_1_0": 0.74, } + n = len(zmat2['symbols']) + new_zmat2 = converter.update_zmat(zmat=zmat2, + element=add_atom, + coords=coords2, + vars=vars, + n=n, + ) + self.assertIsNot(zmat2, new_zmat2) # Test that a copy was made + self.assertEqual(len(new_zmat2['symbols']), len(zmat2['symbols']) + 1) # exactly one atom added + self.assertEqual(new_zmat2['symbols'][-1], add_atom) # the right atom was added at the end + self.assertEqual(new_zmat2['coords'][-1], coords2) + self.assertNotIn(vars.keys(), zmat2['vars'].keys()) + self.assertIn(vars.keys(), new_zmat2['vars'].keys()) + + def test_add_two_xyzs(self): + """Test the add_two_xyzs() function.""" + xyz1 = {'symbols': ('C', 'H', 'H', 'H', 'H'), + 'isotopes': (12, 1, 1, 1, 1), + 'coords': ((0.0, 0.0, 0.0), + (0.0, 0.0, 1.089), + (1.026719, 0.0, -0.363), + (-0.51336, -0.889165, -0.363), + (-0.51336, 0.889165, -0.363))} + xyz2 = {'symbols': ('O', 'H', 'H'), + 'isotopes': (16, 1, 1), + 'coords': ((0.0, 0.0, 0.0), + (0.758602, 0.584882, 0.0), + (-0.758602, 0.584882, 0.0))} + combined_xyz = converter.add_two_xyzs(xyz1=xyz1, + xyz2=xyz2, + atom1_params= {"R": {"val": 5.0, + "m2i": [0], + "m1i": [0], + }, + "A": {"val": 90.0, + "m2i": [0], + "m1i": [0, 1], + }, + "D": {"val": 90.0, + "m2i": [0], + "m1i": [0, 1, 2], + } + }, + atom2_params={"A": {"val": 90.0, + "m2i": [1], + "m1i": [0, 1],}, + "D" : {"val": 90.0, + "m2i": [1], + "m1i": [0, 1, 2],} + }, + atom3_params={"D": {"val": 90.0, + "m2i": [2], + "m1i": [0, 1, 2],} + }, + return_zmat=False, + ) + + + # self.assertTrue(almost_equal_coords_lists(combined_xyz, expected_combined_xyz)) @classmethod def tearDownClass(cls): From cc1366bb58ac2c934e4d26ddbe3b5a27ab9c5262 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Wed, 17 Sep 2025 14:13:39 +0300 Subject: [PATCH 5/6] Added check_zmat_vs_coords This function ensures that the coordinates are consistent with zmat requirements. --- arc/species/zmat.py | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/arc/species/zmat.py b/arc/species/zmat.py index cbeb75ee8b..c5f11355d4 100644 --- a/arc/species/zmat.py +++ b/arc/species/zmat.py @@ -2338,3 +2338,46 @@ def map_index_to_int(index: Union[int, str]) -> int: if isinstance(index, str) and all(char.isdigit() for char in index[1:]): return int(index[1:]) raise TypeError(f'Expected either an int or a string on the format "X15", got {index}') + + +def check_zmat_vs_coords(zmat: dict, + coords: np.array, + rtol: float = 0.05, + atol: float = 1e-4 + ) -> bool: + """ + Check the zmat for consistency vs coordinates. + Args: + zmat (dict): The zmat. + coords (np.array): The cartesian coordinates. + rtol (float, optional): Tolerance to the radius (in ångström). + atol (float, optional): Tolerance to the angle (in degrees). + Returns: + bool: Whether the zmat is consistent with the coordinates. + """ + nones = 0 + for zmat_coords in zmat["coords"]: + for param in zmat_coords: + if param is not None: + if param not in zmat["vars"]: + return False + else: + indices = list(get_atom_indices_from_zmat_parameter(param)[0]) + calculated_param = calculate_param(coords=coords.tolist(), atoms=indices) + zmat_param = zmat["vars"][param] + if abs(calculated_param - zmat_param) > (rtol if param.startswith('R') else atol): + return False + else: + nones += 1 + if nones < 6: + n = len(zmat["symbols"]) + if n == 1: + return nones == 3 # Exactly 3 Nones for one atom + elif n == 2: + return nones == 5 # Exactly 5 Nones for two atoms + else: + return False + elif nones > 6: + return False + return True + \ No newline at end of file From 2d34c22c29bc98904bc03e18e47b9a31c39389b4 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Wed, 17 Sep 2025 14:13:51 +0300 Subject: [PATCH 6/6] test: check_zmat_vs_coords --- arc/species/zmat_test.py | 100 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 98 insertions(+), 2 deletions(-) diff --git a/arc/species/zmat_test.py b/arc/species/zmat_test.py index 2e36ee47d6..aaa5c25ee4 100644 --- a/arc/species/zmat_test.py +++ b/arc/species/zmat_test.py @@ -7,6 +7,8 @@ import unittest +import numpy as np + import arc.species.zmat as zmat from arc.exceptions import ZMatError from arc.species.species import ARCSpecies @@ -2081,7 +2083,101 @@ def test_map_index_to_int(self): self.assertEqual(zmat.map_index_to_int('X5486'), 5486) with self.assertRaises(TypeError): zmat.map_index_to_int('XY5486') - - + + def test_check_zmat_vs_coords(self): + zmat1 = {'coords': ((None, None, None), + ('R_1_0', None, None), + ('R_2_1', 'A_2_1_0', None), + ('R_3_2', 'A_3_2_0', 'D_3_2_0_1'), + ('R_4_3', 'A_4_3_2', 'D_4_3_2_1'), + ('R_5_0', 'A_5_0_1', 'D_5_0_1_2'), + ('R_6_5', 'A_6_0_1', 'D_6_0_1_2'), + ('R_7_6', 'A_7_6_5', 'D_7_0_1_2')), + 'map': {0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6, 7: 7}, + 'symbols': ('C', 'H', 'H', 'H', 'H', 'O', 'H', 'H'), + 'vars': {'A_2_1_0': 35.2643928527832, + 'A_3_2_0': 35.264400482177734, + 'A_4_3_2': 60.000003814697266, + 'A_5_0_1': 90.0, + 'A_6_0_1': 90.0, + 'A_7_6_5': 37.6322135925293, + 'D_3_2_0_1': 120.00000298873425, + 'D_4_3_2_1': 70.52878150493633, + 'D_5_0_1_2': 90.0, + 'D_6_0_1_2': 90.0, + 'D_7_0_1_2': 90.0, + 'R_1_0': 1.0889999866485596, + 'R_2_1': 1.778329610824585, + 'R_3_2': 1.7783300876617432, + 'R_4_3': 1.7783299684524536, + 'R_5_0': 5.0, + 'R_6_5': 0.9578955769538879, + 'R_7_6': 1.5172040462493896}} + coords1 = np.array([[-2.4050857575832123, 2.082149685245476e-08, -0.03557225452565956], + [-2.4050857575832123, 2.082149685245476e-08, 1.0534277321229], + [-2.4050857575832123, 1.0267191806743907, -0.3985722580721713], + [-1.5159205035161651, -0.5133597301749281, -0.3985725441317558], + [-3.2942504719685823, -0.5133600913342189, -0.3985725821758469], + [2.5949142424167877, 2.082149715861646e-08, -0.03557225452565926], + [-1.4471901806293244, 2.0821496911108946e-08, -0.03557225452565951], + [-1.478695351514872, 2.0821496909179812e-08, 1.1659721406175476]]) + self.assertFalse(zmat.check_zmat_vs_coords(zmat1, coords1, rtol=0.05, atol=1e-4)) + + zmat2 = {'coords': ((None, None, None), + ('R_1_0', None, None), + ('R_2_1', 'A_2_1_0', None), + ('R_3_2', 'A_3_2_0', 'D_3_2_0_1'), + ('R_4_3', 'A_4_3_2', 'D_4_3_2_1'), + ('R_5_0', 'A_5_0_1', 'D_5_0_1_2'), + ('R_6_5', 'A_6_5_0', 'D_6_5_0_1'), + ('R_7_6', 'A_7_6_5', 'D_7_6_5_0')), + 'map': {0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6, 7: 7}, + 'symbols': ('C', 'H', 'H', 'H', 'H', 'O', 'H', 'H'), + 'vars': {'A_2_1_0': 35.2643928527832, + 'A_3_2_0': 35.264400482177734, + 'A_4_3_2': 60.000003814697266, + 'A_5_0_1': 90.0, + 'A_6_5_0': 90.0, + 'A_7_6_5': 37.6322135925293, + 'D_3_2_0_1': 120.00000298873425, + 'D_4_3_2_1': 70.52878150493633, + 'D_5_0_1_2': 90.0, + 'D_6_5_0_1': 90.0, + 'D_7_6_5_0': 90.0, + 'R_1_0': 1.0889999866485596, + 'R_2_1': 1.778329610824585, + 'R_3_2': 1.7783300876617432, + 'R_4_3': 1.7783299684524536, + 'R_5_0': 5.0, + 'R_6_5': 0.9578955769538879, + 'R_7_6': 1.5172040462493896} + } + + coords2 = np.array(((-2.645355195986168, 0.021145623937931058, 0.027426231938455496), + (-2.645355195986168, 0.021145623937931058, 1.116426218587015), + (-2.645355195986168, 1.047864783790825, -0.33557377160805624), + (-1.756189941919121, -0.4922141270584939, -0.33557405766764076), + (-3.5345199103715377, -0.4922144882177847, -0.33557409571173186), + (2.354644804013832, 0.021145623937931363, 0.0274262319384558), + (2.354644804013832, -0.9367499530159565, 0.02742623193845586), + (2.354644804013832, 0.2647944421272507, -0.8989641741298846))) + self.assertTrue(zmat.check_zmat_vs_coords(zmat2, coords2, rtol=0.05, atol=1e-4)) + + change_to_none = np.array(zmat2["coords"]) + change_to_none[np.random.randint(3, 7, size=1)[0]][np.random.randint(0, 2, size=1)[0]] = None + zmat2['coords'] = tuple(map(tuple, change_to_none)) + self.assertFalse(zmat.check_zmat_vs_coords(zmat2, coords2, rtol=0.05, atol=1e-4)) + + zmat3 = {'coords': ((None, None, None), ('R_1_0', None, None)), + 'map': {0: 0, 1: 1}, + 'symbols': ('C', 'C'), + 'vars': {'R_1_0': 1.512057900428772}} + coords3 = np.array(((-2.5000000225118915, -2.2511891158743052e-08, 0.0), + (-2.5000000225118915, -2.2511891158743052e-08, 1.5120579004287718),) + ) + self.assertTrue(zmat.check_zmat_vs_coords(zmat3, coords3, rtol=0.05, atol=1e-4)) + zmat3['vars']['R_1_0'] = 1.5 + self.assertFalse(zmat.check_zmat_vs_coords(zmat3, coords3, rtol=1e-3, atol=1e-4)) + if __name__ == '__main__': unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))