From 29e1915dfaddec05e345c4f5d1e0a6508581642b Mon Sep 17 00:00:00 2001 From: crjacinto Date: Fri, 17 Apr 2026 13:59:30 +0200 Subject: [PATCH 01/27] RMSD and Rotationl Constants Modules --- src/opi/utils/rmsd.py | 132 +++++++++++++++ src/opi/utils/rotconst.py | 345 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 477 insertions(+) create mode 100644 src/opi/utils/rmsd.py create mode 100644 src/opi/utils/rotconst.py diff --git a/src/opi/utils/rmsd.py b/src/opi/utils/rmsd.py new file mode 100644 index 00000000..221546d1 --- /dev/null +++ b/src/opi/utils/rmsd.py @@ -0,0 +1,132 @@ +import numpy as np +import os +from typing import Tuple, List, Union + + +def read_xyz( + data: Union[str, Tuple[List[str], np.ndarray]] +) -> Tuple[List[str], np.ndarray]: + """ + Read geometry from: + - XYZ file path + - XYZ block string + - (symbols, coords) tuple + + Returns + ------- + symbols : list[str] + coords : (N, 3) ndarray + """ + + # ------------------------- + # Case 1: string input + # ------------------------- + if isinstance(data, str): + + # Case 1a: a file path + if os.path.isfile(data): + with open(data, "r") as f: + lines = f.readlines() + + # Case 1b: an XYZ block string + else: + lines = data.strip().splitlines() + + n_atoms = int(lines[0].strip()) + symbols = [] + coords = [] + + for line in lines[2:2 + n_atoms]: + parts = line.split() + sym = parts[0] + xyz = [float(x) for x in parts[1:4]] + + symbols.append(sym) + coords.append(xyz) + + return symbols, np.array(coords, dtype=np.float64) + + # ------------------------- + # Case 2: in-memory tuple + # ------------------------- + elif isinstance(data, tuple) and len(data) == 2: + symbols, coords = data + + coords = np.asarray(coords, dtype=np.float64) + + if coords.ndim != 2 or coords.shape[1] != 3: + raise ValueError("Coordinates must be an (N,3) array") + + if len(symbols) != len(coords): + raise ValueError("Symbols and coordinates must have same length") + + return list(symbols), coords + + else: + raise TypeError( + "Input must be: file path, XYZ string, or (symbols, coords)" + ) + + +def _validate_geometries(symA: List[str], symB: List[str]) -> None: + if len(symA) != len(symB): + raise ValueError("Geometries have different number of atoms") + + for i, (a, b) in enumerate(zip(symA, symB)): + if a != b: + raise ValueError(f"Atom mismatch at index {i}: {a} != {b}") + + +def kabsch_rmsd( + ref_xyz: str, + target_xyz: str, + *, + align: bool = True, +) -> float: + """Compute RMSD between two XYZ geometries. + + Parameters + ---------- + ref_xyz : str + Path to reference geometry. + target_xyz : str + Path to target geometry. + align : bool, default True + Whether to perform optimal alignment (Kabsch). + + Returns + ------- + float + RMSD (Å) + """ + + symA, A = read_xyz(ref_xyz) + symB, B = read_xyz(target_xyz) + + _validate_geometries(symA, symB) + + # Center using centroid (simple average) + A_cent = A - A.mean(axis=0) + B_cent = B - B.mean(axis=0) + + if not align: + diff = A_cent - B_cent + return float(np.sqrt(np.sum(diff**2) / len(A_cent))) + + # Standard Kabsch covariance matrix + H = B_cent.T @ A_cent + + U, _, Vt = np.linalg.svd(H) + + d = np.linalg.det(Vt.T @ U.T) + D = np.diag([1.0, 1.0, d]) + + R = Vt.T @ D @ U.T + + B_rot = B_cent @ R + + diff = A_cent - B_rot + + rmsd = np.sqrt(np.sum(diff**2) / len(A_cent)) + + return float(rmsd) diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py new file mode 100644 index 00000000..d43ad295 --- /dev/null +++ b/src/opi/utils/rotconst.py @@ -0,0 +1,345 @@ +import numpy as np +import os +from dataclasses import dataclass +import warnings +from typing import Iterable, Union, Tuple, List + + +# ============================================================ +# Physical constants +# ============================================================ +_AMU_TO_KG = 1.66053906660e-27 +_ANGSTROM_TO_M = 1.0e-10 +_H_PLANCK = 6.62607015e-34 +_C_CM = 2.99792458e10 + + +# ============================================================ +# Atomic masses +# ============================================================ +ATOMIC_MASSES = { + "X": 0.0, + "PointCharge": 0.0, + + "H": 1.008, "He": 4.003, + "Li": 6.941, "Be": 9.012, "B": 10.810, "C": 12.011, "N": 14.007, + "O": 15.999, "F": 18.998, "Ne": 20.179, + + "Na": 22.990, "Mg": 24.305, "Al": 26.982, "Si": 28.086, + "P": 30.974, "S": 32.060, "Cl": 35.453, "Ar": 39.948, + + "K": 39.100, "Ca": 40.080, + + "Sc": 44.960, "Ti": 47.900, "V": 50.940, "Cr": 52.000, + "Mn": 54.940, "Fe": 55.850, "Co": 58.930, "Ni": 58.700, + "Cu": 63.550, "Zn": 65.380, + + "Ga": 69.720, "Ge": 72.590, "As": 74.920, "Se": 78.960, + "Br": 79.900, "Kr": 83.800, + + "Rb": 85.479, "Sr": 87.620, + + "Y": 88.910, "Zr": 91.220, "Nb": 92.910, "Mo": 95.940, + "Tc": 97.000, "Ru": 101.070, "Rh": 102.910, "Pd": 106.400, + "Ag": 107.870, "Cd": 112.410, + + "In": 114.820, "Sn": 118.690, "Sb": 121.750, "Te": 127.600, + "I": 126.900, "Xe": 131.300, + + "Cs": 132.9054, "Ba": 137.3300, + + "La": 138.9055, "Ce": 140.1200, "Pr": 140.9077, "Nd": 144.2400, + "Pm": 145.0000, "Sm": 150.4000, "Eu": 151.9600, "Gd": 157.2500, + "Tb": 158.9254, "Dy": 162.5000, "Ho": 164.9304, "Er": 167.2600, + "Tm": 168.9342, "Yb": 173.0400, "Lu": 174.9670, + + "Hf": 178.4900, "Ta": 180.9479, "W": 183.8500, "Re": 186.2070, + "Os": 190.2000, "Ir": 192.2200, "Pt": 195.0900, "Au": 196.9665, + "Hg": 200.5900, + + "Tl": 204.3700, "Pb": 207.2000, "Bi": 208.9804, "Po": 209.0000, + "At": 210.0000, "Rn": 222.0000, + + "Fr": 223.0000, "Ra": 226.0254, + + "Ac": 227.0278, "Th": 232.0381, "Pa": 231.0359, "U": 238.0290, + "Np": 237.0482, "Pu": 244.0000, "Am": 243.0000, "Cm": 247.0000, + "Bk": 247.0000, "Cf": 251.0000, "Es": 252.0000, "Fm": 257.0000, + "Md": 258.0000, "No": 259.0000, "Lr": 262.0000, + + "Rf": 267.0000, "Db": 268.0000, "Sg": 269.0000, "Bh": 270.0000, + "Hs": 269.0000, "Mt": 278.0000, "Ds": 281.0000, "Rg": 281.0000, + "Cn": 285.0000, + + "Nh": 283.0000, "Fl": 289.0000, "Mc": 288.0000, + "Lv": 293.0000, "Ts": 294.0000, "Og": 294.0000, +} + +# ============================================================ +# Dataclass +# ============================================================ +@dataclass +class RotationalConstants: + A: float | None + B: float | None + C: float | None + A_cm: float | None + B_cm: float | None + C_cm: float | None + moments: tuple[float, float, float] + rotor_type: str + + def __str__(self) -> str: + def fmt(x, unit=""): + return f"{x:.6f} {unit}" if x is not None else "None" + + return ( + # "Rotational Spectrum\n" + # "--------------------\n" + f"Rotor type : {self.rotor_type}\n\n" + "Moments of inertia (amu·Å²):\n" + f" Ia = {self.moments[0]:.6f}\n" + f" Ib = {self.moments[1]:.6f}\n" + f" Ic = {self.moments[2]:.6f}\n\n" + "Rotational constants:\n" + f" A = {fmt(self.A, 'MHz')} ({fmt(self.A_cm, 'cm⁻¹')})\n" + f" B = {fmt(self.B, 'MHz')} ({fmt(self.B_cm, 'cm⁻¹')})\n" + f" C = {fmt(self.C, 'MHz')} ({fmt(self.C_cm, 'cm⁻¹')})" + ) + + +# ============================================================ +# Utilities +# ============================================================ +def _normalize_symbol(s: str) -> str: + if s in ("X", "PointCharge"): + return s + return s.capitalize() + + +def _read_xyz( + data: Union[str, Tuple[List[str], np.ndarray]] +) -> Tuple[List[str], np.ndarray]: + """ + Read geometry from: + - XYZ file path + - XYZ block string + - (symbols, coords) tuple + + Returns + ------- + symbols : list[str] + coords : (N, 3) ndarray + """ + + # ------------------------- + # Case 1: string input + # ------------------------- + if isinstance(data, str): + + # Case 1a: file path + if os.path.isfile(data): + with open(data, "r") as f: + lines = f.readlines() + else: + # Case 1b: XYZ block string + lines = data.strip().splitlines() + + if len(lines) < 3: + raise ValueError("Invalid XYZ format: too few lines") + + try: + n_atoms = int(lines[0].strip()) + except ValueError: + raise ValueError("First line must contain number of atoms") + + symbols = [] + coords = [] + + for i, line in enumerate(lines[2:2 + n_atoms], start=3): + parts = line.split() + + if len(parts) < 4: + raise ValueError(f"Line {i} malformed: '{line}'") + + # --- clean symbol --- + raw_sym = parts[0].strip() + + # Remove numeric labels (e.g. C1 → C) + sym = ''.join(filter(str.isalpha, raw_sym)) + + # Capitalization (cap insensitive) + sym = sym.capitalize() + + try: + xyz = [float(x) for x in parts[1:4]] + except ValueError: + raise ValueError(f"Invalid coordinates at line {i}: '{line}'") + + symbols.append(sym) + coords.append(xyz) + + return symbols, np.array(coords, dtype=np.float64) + + # ------------------------- + # Case 2: in-memory tuple + # ------------------------- + elif isinstance(data, tuple) and len(data) == 2: + symbols, coords = data + + coords = np.asarray(coords, dtype=np.float64) + + if coords.ndim != 2 or coords.shape[1] != 3: + raise ValueError("Coordinates must be an (N,3) array") + + if len(symbols) != len(coords): + raise ValueError("Symbols and coordinates must have same length") + + # Clean symbols here as well (consistency!) + clean_symbols = [] + for s in symbols: + s = ''.join(filter(str.isalpha, str(s))) + clean_symbols.append(s.capitalize()) + + return clean_symbols, coords + + else: + raise TypeError( + "Input must be: file path, XYZ string, or (symbols, coords)" + ) + + +# ============================================================ +# Main function +# ============================================================ +def rotational_constants( + symbols: list[str] | None = None, + coords: np.ndarray | None = None, + xyz: str | None = None, + masses: np.ndarray | None = None, + weights: dict[str, float] | None = None, + atom_weights: dict[int, float] | None = None, +) -> RotationalConstants | None: + """ + Flexible rotational constant calculator. + + Input options + ------------- + - symbols + coords + - xyz (file path, string, or lines) + + Mass priority + ------------- + masses > atom_weights > weights > default + + Unknown atoms + ------------- + Assigned mass = 0 with warning (unless overridden). + """ + + # --- Input parsing --- + if xyz is not None: + symbols, coords = _read_xyz(xyz) + + if symbols is None or coords is None: + raise ValueError("Provide either (symbols, coords) or xyz input.") + + coords = np.asarray(coords, dtype=np.float64) + + # --- Normalize symbols --- + symbols = [_normalize_symbol(s) for s in symbols] + + # --- Prepare weights --- + weights = { _normalize_symbol(k): v for k, v in (weights or {}).items() } + atom_weights = atom_weights or {} + + # --- Assign masses --- + if masses is not None: + masses = np.asarray(masses, dtype=np.float64) + + else: + masses_list = [] + for i, s in enumerate(symbols): + + if i in atom_weights: + m = atom_weights[i] + + elif s in weights: + m = weights[s] + + elif s in ATOMIC_MASSES: + m = ATOMIC_MASSES[s] + + else: + warnings.warn(f"Unknown element '{s}' → mass set to 0.0") + m = 0.0 + + masses_list.append(m) + + masses = np.array(masses_list, dtype=np.float64) + + # --- Filter zero-mass atoms --- + mask = masses > 0.0 + if not np.any(mask): + return None + + masses = masses[mask] + coords = coords[mask] + + total_mass = masses.sum() + + # --- Center of mass --- + com = (masses[:, None] * coords).sum(axis=0) / total_mass + coords -= com + + # --- Inertia tensor --- + inertia = np.zeros((3, 3), dtype=np.float64) + for m, r in zip(masses, coords): + inertia += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r)) + + # --- Diagonalize --- + moments_raw, _ = np.linalg.eigh(inertia) + moments_raw = np.maximum(moments_raw, 0.0) + Ia, Ib, Ic = moments_raw + + # --- Convert to rotational constants --- + def _moment_to_mhz(I): + if I < 1e-6: + return None + I_si = I * _AMU_TO_KG * (_ANGSTROM_TO_M ** 2) + return _H_PLANCK / (8.0 * np.pi**2 * I_si) / 1e6 + + def _mhz_to_cm(mhz): + return None if mhz is None else mhz * 1e6 / _C_CM + + A = _moment_to_mhz(Ia) + B = _moment_to_mhz(Ib) + C = _moment_to_mhz(Ic) + + # --- Rotor classification --- + tol = 1e-3 + n_zero = sum(m < 1e-6 for m in (Ia, Ib, Ic)) + + if n_zero == 3: + rotor = "monoatomic" + elif n_zero == 2: + rotor = "linear" + elif abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: + rotor = "spherical top" + elif abs(Ia - Ib) < tol: + rotor = "symmetric top (oblate)" + elif abs(Ib - Ic) < tol: + rotor = "symmetric top (prolate)" + else: + rotor = "asymmetric top" + + return RotationalConstants( + A=A, + B=B, + C=C, + A_cm=_mhz_to_cm(A), + B_cm=_mhz_to_cm(B), + C_cm=_mhz_to_cm(C), + moments=(Ia, Ib, Ic), + rotor_type=rotor, + ) \ No newline at end of file From 21df17b25038b8bb36bdc8bff3c3fdb2b6be3c97 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Fri, 17 Apr 2026 15:58:58 +0200 Subject: [PATCH 02/27] Fixed missing type annotations --- src/opi/utils/rmsd.py | 35 +----- src/opi/utils/rotconst.py | 235 ++++++++++++++++++++++---------------- 2 files changed, 143 insertions(+), 127 deletions(-) diff --git a/src/opi/utils/rmsd.py b/src/opi/utils/rmsd.py index 221546d1..2c75ff68 100644 --- a/src/opi/utils/rmsd.py +++ b/src/opi/utils/rmsd.py @@ -1,11 +1,10 @@ -import numpy as np import os -from typing import Tuple, List, Union +from typing import List, Tuple, Union + +import numpy as np -def read_xyz( - data: Union[str, Tuple[List[str], np.ndarray]] -) -> Tuple[List[str], np.ndarray]: +def read_xyz(data: Union[str, Tuple[List[str], np.ndarray]]) -> Tuple[List[str], np.ndarray]: """ Read geometry from: - XYZ file path @@ -18,11 +17,7 @@ def read_xyz( coords : (N, 3) ndarray """ - # ------------------------- - # Case 1: string input - # ------------------------- if isinstance(data, str): - # Case 1a: a file path if os.path.isfile(data): with open(data, "r") as f: @@ -36,7 +31,7 @@ def read_xyz( symbols = [] coords = [] - for line in lines[2:2 + n_atoms]: + for line in lines[2 : 2 + n_atoms]: parts = line.split() sym = parts[0] xyz = [float(x) for x in parts[1:4]] @@ -46,26 +41,8 @@ def read_xyz( return symbols, np.array(coords, dtype=np.float64) - # ------------------------- - # Case 2: in-memory tuple - # ------------------------- - elif isinstance(data, tuple) and len(data) == 2: - symbols, coords = data - - coords = np.asarray(coords, dtype=np.float64) - - if coords.ndim != 2 or coords.shape[1] != 3: - raise ValueError("Coordinates must be an (N,3) array") - - if len(symbols) != len(coords): - raise ValueError("Symbols and coordinates must have same length") - - return list(symbols), coords - else: - raise TypeError( - "Input must be: file path, XYZ string, or (symbols, coords)" - ) + raise TypeError("Input must be: file path or XYZ string") def _validate_geometries(symA: List[str], symB: List[str]) -> None: diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index d43ad295..e82e9e78 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -1,17 +1,17 @@ -import numpy as np import os -from dataclasses import dataclass import warnings -from typing import Iterable, Union, Tuple, List +from dataclasses import dataclass +from typing import List, Tuple, Union +import numpy as np # ============================================================ # Physical constants # ============================================================ -_AMU_TO_KG = 1.66053906660e-27 -_ANGSTROM_TO_M = 1.0e-10 -_H_PLANCK = 6.62607015e-34 -_C_CM = 2.99792458e10 +_AMU_TO_KG = 1.66053906660e-27 +_ANGSTROM_TO_M = 1.0e-10 +_H_PLANCK = 6.62607015e-34 +_C_CM = 2.99792458e10 # ============================================================ @@ -20,61 +20,127 @@ ATOMIC_MASSES = { "X": 0.0, "PointCharge": 0.0, - - "H": 1.008, "He": 4.003, - "Li": 6.941, "Be": 9.012, "B": 10.810, "C": 12.011, "N": 14.007, - "O": 15.999, "F": 18.998, "Ne": 20.179, - - "Na": 22.990, "Mg": 24.305, "Al": 26.982, "Si": 28.086, - "P": 30.974, "S": 32.060, "Cl": 35.453, "Ar": 39.948, - - "K": 39.100, "Ca": 40.080, - - "Sc": 44.960, "Ti": 47.900, "V": 50.940, "Cr": 52.000, - "Mn": 54.940, "Fe": 55.850, "Co": 58.930, "Ni": 58.700, - "Cu": 63.550, "Zn": 65.380, - - "Ga": 69.720, "Ge": 72.590, "As": 74.920, "Se": 78.960, - "Br": 79.900, "Kr": 83.800, - - "Rb": 85.479, "Sr": 87.620, - - "Y": 88.910, "Zr": 91.220, "Nb": 92.910, "Mo": 95.940, - "Tc": 97.000, "Ru": 101.070, "Rh": 102.910, "Pd": 106.400, - "Ag": 107.870, "Cd": 112.410, - - "In": 114.820, "Sn": 118.690, "Sb": 121.750, "Te": 127.600, - "I": 126.900, "Xe": 131.300, - - "Cs": 132.9054, "Ba": 137.3300, - - "La": 138.9055, "Ce": 140.1200, "Pr": 140.9077, "Nd": 144.2400, - "Pm": 145.0000, "Sm": 150.4000, "Eu": 151.9600, "Gd": 157.2500, - "Tb": 158.9254, "Dy": 162.5000, "Ho": 164.9304, "Er": 167.2600, - "Tm": 168.9342, "Yb": 173.0400, "Lu": 174.9670, - - "Hf": 178.4900, "Ta": 180.9479, "W": 183.8500, "Re": 186.2070, - "Os": 190.2000, "Ir": 192.2200, "Pt": 195.0900, "Au": 196.9665, + "H": 1.008, + "He": 4.003, + "Li": 6.941, + "Be": 9.012, + "B": 10.810, + "C": 12.011, + "N": 14.007, + "O": 15.999, + "F": 18.998, + "Ne": 20.179, + "Na": 22.990, + "Mg": 24.305, + "Al": 26.982, + "Si": 28.086, + "P": 30.974, + "S": 32.060, + "Cl": 35.453, + "Ar": 39.948, + "K": 39.100, + "Ca": 40.080, + "Sc": 44.960, + "Ti": 47.900, + "V": 50.940, + "Cr": 52.000, + "Mn": 54.940, + "Fe": 55.850, + "Co": 58.930, + "Ni": 58.700, + "Cu": 63.550, + "Zn": 65.380, + "Ga": 69.720, + "Ge": 72.590, + "As": 74.920, + "Se": 78.960, + "Br": 79.900, + "Kr": 83.800, + "Rb": 85.479, + "Sr": 87.620, + "Y": 88.910, + "Zr": 91.220, + "Nb": 92.910, + "Mo": 95.940, + "Tc": 97.000, + "Ru": 101.070, + "Rh": 102.910, + "Pd": 106.400, + "Ag": 107.870, + "Cd": 112.410, + "In": 114.820, + "Sn": 118.690, + "Sb": 121.750, + "Te": 127.600, + "I": 126.900, + "Xe": 131.300, + "Cs": 132.9054, + "Ba": 137.3300, + "La": 138.9055, + "Ce": 140.1200, + "Pr": 140.9077, + "Nd": 144.2400, + "Pm": 145.0000, + "Sm": 150.4000, + "Eu": 151.9600, + "Gd": 157.2500, + "Tb": 158.9254, + "Dy": 162.5000, + "Ho": 164.9304, + "Er": 167.2600, + "Tm": 168.9342, + "Yb": 173.0400, + "Lu": 174.9670, + "Hf": 178.4900, + "Ta": 180.9479, + "W": 183.8500, + "Re": 186.2070, + "Os": 190.2000, + "Ir": 192.2200, + "Pt": 195.0900, + "Au": 196.9665, "Hg": 200.5900, - - "Tl": 204.3700, "Pb": 207.2000, "Bi": 208.9804, "Po": 209.0000, - "At": 210.0000, "Rn": 222.0000, - - "Fr": 223.0000, "Ra": 226.0254, - - "Ac": 227.0278, "Th": 232.0381, "Pa": 231.0359, "U": 238.0290, - "Np": 237.0482, "Pu": 244.0000, "Am": 243.0000, "Cm": 247.0000, - "Bk": 247.0000, "Cf": 251.0000, "Es": 252.0000, "Fm": 257.0000, - "Md": 258.0000, "No": 259.0000, "Lr": 262.0000, - - "Rf": 267.0000, "Db": 268.0000, "Sg": 269.0000, "Bh": 270.0000, - "Hs": 269.0000, "Mt": 278.0000, "Ds": 281.0000, "Rg": 281.0000, + "Tl": 204.3700, + "Pb": 207.2000, + "Bi": 208.9804, + "Po": 209.0000, + "At": 210.0000, + "Rn": 222.0000, + "Fr": 223.0000, + "Ra": 226.0254, + "Ac": 227.0278, + "Th": 232.0381, + "Pa": 231.0359, + "U": 238.0290, + "Np": 237.0482, + "Pu": 244.0000, + "Am": 243.0000, + "Cm": 247.0000, + "Bk": 247.0000, + "Cf": 251.0000, + "Es": 252.0000, + "Fm": 257.0000, + "Md": 258.0000, + "No": 259.0000, + "Lr": 262.0000, + "Rf": 267.0000, + "Db": 268.0000, + "Sg": 269.0000, + "Bh": 270.0000, + "Hs": 269.0000, + "Mt": 278.0000, + "Ds": 281.0000, + "Rg": 281.0000, "Cn": 285.0000, - - "Nh": 283.0000, "Fl": 289.0000, "Mc": 288.0000, - "Lv": 293.0000, "Ts": 294.0000, "Og": 294.0000, + "Nh": 283.0000, + "Fl": 289.0000, + "Mc": 288.0000, + "Lv": 293.0000, + "Ts": 294.0000, + "Og": 294.0000, } + # ============================================================ # Dataclass # ============================================================ @@ -90,7 +156,7 @@ class RotationalConstants: rotor_type: str def __str__(self) -> str: - def fmt(x, unit=""): + def fmt(x: float | None, unit: str = "") -> str: return f"{x:.6f} {unit}" if x is not None else "None" return ( @@ -117,9 +183,7 @@ def _normalize_symbol(s: str) -> str: return s.capitalize() -def _read_xyz( - data: Union[str, Tuple[List[str], np.ndarray]] -) -> Tuple[List[str], np.ndarray]: +def _read_xyz(data: Union[str, Tuple[List[str], np.ndarray]]) -> Tuple[List[str], np.ndarray]: """ Read geometry from: - XYZ file path @@ -136,7 +200,6 @@ def _read_xyz( # Case 1: string input # ------------------------- if isinstance(data, str): - # Case 1a: file path if os.path.isfile(data): with open(data, "r") as f: @@ -156,7 +219,7 @@ def _read_xyz( symbols = [] coords = [] - for i, line in enumerate(lines[2:2 + n_atoms], start=3): + for i, line in enumerate(lines[2 : 2 + n_atoms], start=3): parts = line.split() if len(parts) < 4: @@ -166,7 +229,7 @@ def _read_xyz( raw_sym = parts[0].strip() # Remove numeric labels (e.g. C1 → C) - sym = ''.join(filter(str.isalpha, raw_sym)) + sym = "".join(filter(str.isalpha, raw_sym)) # Capitalization (cap insensitive) sym = sym.capitalize() @@ -181,32 +244,8 @@ def _read_xyz( return symbols, np.array(coords, dtype=np.float64) - # ------------------------- - # Case 2: in-memory tuple - # ------------------------- - elif isinstance(data, tuple) and len(data) == 2: - symbols, coords = data - - coords = np.asarray(coords, dtype=np.float64) - - if coords.ndim != 2 or coords.shape[1] != 3: - raise ValueError("Coordinates must be an (N,3) array") - - if len(symbols) != len(coords): - raise ValueError("Symbols and coordinates must have same length") - - # Clean symbols here as well (consistency!) - clean_symbols = [] - for s in symbols: - s = ''.join(filter(str.isalpha, str(s))) - clean_symbols.append(s.capitalize()) - - return clean_symbols, coords - else: - raise TypeError( - "Input must be: file path, XYZ string, or (symbols, coords)" - ) + raise TypeError("Input must be: file path or XYZ string") # ============================================================ @@ -215,7 +254,7 @@ def _read_xyz( def rotational_constants( symbols: list[str] | None = None, coords: np.ndarray | None = None, - xyz: str | None = None, + xyz: str | None = None, masses: np.ndarray | None = None, weights: dict[str, float] | None = None, atom_weights: dict[int, float] | None = None, @@ -250,7 +289,7 @@ def rotational_constants( symbols = [_normalize_symbol(s) for s in symbols] # --- Prepare weights --- - weights = { _normalize_symbol(k): v for k, v in (weights or {}).items() } + weights = {_normalize_symbol(k): v for k, v in (weights or {}).items()} atom_weights = atom_weights or {} # --- Assign masses --- @@ -260,7 +299,6 @@ def rotational_constants( else: masses_list = [] for i, s in enumerate(symbols): - if i in atom_weights: m = atom_weights[i] @@ -294,6 +332,7 @@ def rotational_constants( # --- Inertia tensor --- inertia = np.zeros((3, 3), dtype=np.float64) + assert coords is not None for m, r in zip(masses, coords): inertia += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r)) @@ -303,13 +342,13 @@ def rotational_constants( Ia, Ib, Ic = moments_raw # --- Convert to rotational constants --- - def _moment_to_mhz(I): - if I < 1e-6: + def _moment_to_mhz(inertia: float) -> float | None: + if inertia < 1e-6: return None - I_si = I * _AMU_TO_KG * (_ANGSTROM_TO_M ** 2) + I_si = inertia * _AMU_TO_KG * (_ANGSTROM_TO_M**2) return _H_PLANCK / (8.0 * np.pi**2 * I_si) / 1e6 - def _mhz_to_cm(mhz): + def _mhz_to_cm(mhz: float | None) -> float | None: return None if mhz is None else mhz * 1e6 / _C_CM A = _moment_to_mhz(Ia) @@ -342,4 +381,4 @@ def _mhz_to_cm(mhz): C_cm=_mhz_to_cm(C), moments=(Ia, Ib, Ic), rotor_type=rotor, - ) \ No newline at end of file + ) From 02f9514d5ce21110cc3d3e7f723c61d00a9b5134 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Mon, 20 Apr 2026 10:14:48 +0200 Subject: [PATCH 03/27] Fixed error in documentation --- src/opi/utils/rmsd.py | 5 ++--- src/opi/utils/rotconst.py | 1 - 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/opi/utils/rmsd.py b/src/opi/utils/rmsd.py index 2c75ff68..3adf9bdd 100644 --- a/src/opi/utils/rmsd.py +++ b/src/opi/utils/rmsd.py @@ -9,7 +9,6 @@ def read_xyz(data: Union[str, Tuple[List[str], np.ndarray]]) -> Tuple[List[str], Read geometry from: - XYZ file path - XYZ block string - - (symbols, coords) tuple Returns ------- @@ -65,9 +64,9 @@ def kabsch_rmsd( Parameters ---------- ref_xyz : str - Path to reference geometry. + Reference geometry. target_xyz : str - Path to target geometry. + Target geometry. align : bool, default True Whether to perform optimal alignment (Kabsch). diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index e82e9e78..d7c06cbd 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -188,7 +188,6 @@ def _read_xyz(data: Union[str, Tuple[List[str], np.ndarray]]) -> Tuple[List[str] Read geometry from: - XYZ file path - XYZ block string - - (symbols, coords) tuple Returns ------- From a901cb2fe135445252636b65d825707ff8038e5f Mon Sep 17 00:00:00 2001 From: crjacinto Date: Wed, 29 Apr 2026 10:22:00 +0200 Subject: [PATCH 04/27] Fixed RMSD and RotConst implementation. Modified structure.py and element.py. Added units.py and rotconts.py --- src/opi/input/structures/structure.py | 324 +++++++++++++++++++- src/opi/utils/constants.py | 7 + src/opi/utils/element.py | 255 ++++++++++++++++ src/opi/utils/rotconst.py | 406 +++----------------------- src/opi/utils/units.py | 4 + 5 files changed, 634 insertions(+), 362 deletions(-) create mode 100644 src/opi/utils/constants.py diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 8dccf329..95bae8fc 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -19,7 +19,9 @@ PointCharge, ) from opi.input.structures.coordinates import Coordinates -from opi.utils.element import Element +from opi.utils.element import Element, ATOMIC_MASSES_FROM_ELEMENT +from opi.utils import units, constants +from opi.utils.rotconst import * from opi.utils.tracking_text_io import TrackingTextIO __all__ = ("Structure",) @@ -32,7 +34,6 @@ RGX_FRAG_ID = re.compile(r"(?<=\()\d+(?=\))") RGX_ATOM_SYMBOL_FRAG_ID = re.compile(r"(?P[A-Za-z]{1,2})(\((?P\d+)\))?") - class Structure: """ Class to model internal structure for ORCA calculations. @@ -1066,3 +1067,322 @@ def _iter_xyz_structures( n_struc += 1 if n_struc_limit and n_struc >= n_struc_limit: break + + def centered_structure(self) -> "Structure": + """ + Return a new Structure centered at its centroid. + + Returns + ------- + Structure + New Structure with centered coordinates. + """ + import copy + + coords = np.array([a.coordinates.coordinates for a in self.atoms], dtype=np.float64) + centered_coords = coords - coords.mean(axis=0) + + new_structure = copy.deepcopy(self) + for atom, new_coord in zip(new_structure.atoms, centered_coords): + atom.coordinates = new_coord + + return new_structure + + def _filtered_atoms( + self, + only_atoms: Sequence[int], + ignore_hs: bool, + ) -> list[Atom]: + """ + Return Atom instances after applying only_atoms and ignore_hs filters. + """ + atom_list = self.atoms + if only_atoms: + candidates = [atom_list[i] for i in only_atoms] + elif ignore_hs: + candidates = [a for a in atom_list if isinstance(a, Atom) and a.element != Element.H] + else: + candidates = [a for a in atom_list if isinstance(a, Atom)] + return [a for a in candidates if isinstance(a, Atom)] + + @staticmethod + def _validate_rmsd_compatibility(atoms1: list[Atom], atoms2: list[Atom]) -> None: + """ + Raise ValueError if atoms1 and atoms2 differ in count or element order. + """ + if len(atoms1) != len(atoms2): + raise ValueError( + f"Structures have different number of atoms: {len(atoms1)} vs {len(atoms2)}" + ) + for i, (a, b) in enumerate(zip(atoms1, atoms2)): + if a.element != b.element: + raise ValueError(f"Atom mismatch at index {i}: {a.element!r} != {b.element!r}") + + @staticmethod + def _rmsd_coords( + coords1: npt.NDArray[np.float64], + coords2: npt.NDArray[np.float64], + ) -> float: + """ + Compute RMSD between two aligned (N, 3) coordinate arrays. + """ + diff = coords1 - coords2 + return float(np.sqrt(np.sum(diff**2) / len(coords1))) + + def rmsd( + self, + other: "Structure", + /, + only_atoms: Sequence[int] = (), + *, + ignore_hs: bool = False, + ) -> float: + """ + Compute RMSD between this structure and *other* (no rotational alignment). + + Both structures are translated to their centroid before comparison. + + Parameters + ---------- + other : Structure + Structure to compare against. + only_atoms : Sequence[int], default () + Atom indices to include. If non-empty, ``ignore_hs`` is ignored. + ignore_hs : bool, default False + Exclude hydrogen atoms from the RMSD computation. + + Returns + ------- + float + RMSD in Ångström. + + Raises + ------ + ValueError + If the filtered atom sets differ in size or element order. + """ + atoms1 = self._filtered_atoms(only_atoms, ignore_hs) + atoms2 = other._filtered_atoms(only_atoms, ignore_hs) + self._validate_rmsd_compatibility(atoms1, atoms2) + + coords1 = np.array([a.coordinates.coordinates for a in atoms1], dtype=np.float64) + coords2 = np.array([a.coordinates.coordinates for a in atoms2], dtype=np.float64) + + coords1 -= coords1.mean(axis=0) + coords2 -= coords2.mean(axis=0) + + return self._rmsd_coords(coords1, coords2) + + def rmsd_kabsch( + self, + other: "Structure", + /, + only_atoms: Sequence[int] = (), + *, + ignore_hs: bool = False, + ) -> float: + """ + Compute RMSD between this structure and *other* using the Kabsch algorithm. + + Translates both structures to their centroid, then finds the optimal + rotation matrix via SVD before computing RMSD. + + Parameters + ---------- + other : Structure + Structure to compare against. + only_atoms : Sequence[int], default () + Atom indices to include. If non-empty, ``ignore_hs`` is ignored. + ignore_hs : bool, default False + Exclude hydrogen atoms from the RMSD computation. + + Returns + ------- + float + RMSD in Ångström after optimal rotational alignment. + + Raises + ------ + ValueError + If the filtered atom sets differ in size or element order. + """ + atoms1 = self._filtered_atoms(only_atoms, ignore_hs) + atoms2 = other._filtered_atoms(only_atoms, ignore_hs) + self._validate_rmsd_compatibility(atoms1, atoms2) + + coords1 = np.array([a.coordinates.coordinates for a in atoms1], dtype=np.float64) + coords2 = np.array([a.coordinates.coordinates for a in atoms2], dtype=np.float64) + + # Translate to centroid + coords1 -= coords1.mean(axis=0) + coords2_centered = coords2 - coords2.mean(axis=0) + + # Kabsch algorithm: find optimal rotation matrix via SVD + # https://doi.org/10.1107/S0567739476001873 + + # ------------------------------------------------------------------ + # Build covariance matrix + # ------------------------------------------------------------------ + # H = B^T A + # This matrix captures how coordinates from B map onto A + H = coords2_centered.T @ coords1 + + # ------------------------------------------------------------------ + # Singular Value Decomposition (SVD) + # ------------------------------------------------------------------ + # H = U S V^T + # This decomposes the transformation into rotations + scaling + U, _, Vt = np.linalg.svd(H) + + # ------------------------------------------------------------------ + # Optimal rotation matrix + # ------------------------------------------------------------------ + # R = V U^T + # Reflection correction: ensure det(R) = +1 + d = np.linalg.det(Vt.T @ U.T) + D = np.diag([1.0, 1.0, d]) + + R = Vt.T @ D @ U.T + + coords2_rotated = coords2_centered @ R + + return self._rmsd_coords(coords1, coords2_rotated) + + def rotational_constants( + self, + masses: npt.NDArray[np.float64] | None = None, + weights: dict[str, float] | None = None, + atom_weights: dict[int, float] | None = None, + ) ->RotationalConstants | None: + """ + Compute rotational constants for this structure. + + Only Atom instances contribute; GhostAtom, PointCharge, and + EmbeddingPotential atoms are silently ignored. + + Mass priority + ------------- + masses > atom_weights > weights > default (ATOMIC_MASSES_FROM_ELEMENT) + + Parameters + ---------- + masses : npt.NDArray[np.float64] | None, default None + Per-atom masses (amu) for all Atom instances, overriding every + other mass source. + weights : dict[str, float] | None, default None + Per-element mass overrides keyed by element symbol string + (e.g. ``{"C": 13.003}``) + atom_weights : dict[int, float] | None, default None + Per-atom mass overrides keyed by index within the Atom list. + + Returns + ------- + RotationalConstants | None + None if no Atom instances are present or all masses are zero. + + Raises + ------ + ValueError + If the length of *masses* does not match the number of Atom + instances in the structure. + """ + # --- Collect Atom instances only --- + atom_list = [a for a in self.atoms if isinstance(a, Atom)] + + if not atom_list: + return None + + coords = np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) + + # --- Prepare weight overrides --- + weights = {k: v for k, v in (weights or {}).items()} + atom_weights = atom_weights or {} + + # --- Assign masses --- + if masses is not None: + masses = np.asarray(masses, dtype=np.float64) + if len(masses) != len(atom_list): + raise ValueError( + f"masses length ({len(masses)}) does not match number of atoms ({len(atom_list)})" + ) + else: + masses_list: list[float] = [] + for i, atom in enumerate(atom_list): + if i in atom_weights: + m = atom_weights[i] + elif atom.element.value in weights: + m = weights[atom.element.value] + elif atom.element in ATOMIC_MASSES_FROM_ELEMENT: + m = ATOMIC_MASSES_FROM_ELEMENT[atom.element] + else: + warn(f"Unknown element '{atom.element.value}' → mass set to 0.0") + m = 0.0 + masses_list.append(m) + masses = np.array(masses_list, dtype=np.float64) + + # --- Filter zero-mass atoms --- + mask = masses > 0.0 + if not np.any(mask): + return None + + masses = masses[mask] + coords = coords[mask] + total_mass = float(masses.sum()) + + # --- Center of mass --- + com = (masses[:, None] * coords).sum(axis=0) / total_mass + coords -= com + + # --- Inertia tensor --- + inertia = np.zeros((3, 3), dtype=np.float64) + for m, r in zip(masses, coords): + inertia += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r)) + + # --- Diagonalize: eigenvalues are the principal moments --- + moments_raw, _ = np.linalg.eigh(inertia) + moments_raw = np.maximum(moments_raw, 0.0) + Ia, Ib, Ic = moments_raw + + # --- Convert moments (amu·Å²) to rotational constants (MHz) --- + def _moment_to_mhz(I: float) -> float | None: + if I < 1e-6: + return None + I_si = I * units.AMU_TO_KG * (units.ANGST_TO_M ** 2) + return constants.H_PLANCK / (8.0 * np.pi ** 2 * I_si) / 1e6 + + def _mhz_to_cm(mhz: float | None) -> float | None: + # MHz → cm⁻¹ : divide by speed of light in cm/s + return None if mhz is None else mhz * 1e6 / constants.C + + A = _moment_to_mhz(Ia) + B = _moment_to_mhz(Ib) + C = _moment_to_mhz(Ic) + + # --- Rotor classification --- + tol = 1e-3 + n_zero = sum(moment < 1e-6 for moment in (Ia, Ib, Ic)) + + if n_zero == 3: + rotor = RotorType.MONOATOMIC + elif n_zero == 2: + rotor = RotorType.LINEAR + elif abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: + rotor = RotorType.SPHERICAL_TOP + elif abs(Ia - Ib) < tol: + rotor = RotorType.OBLATE_TOP + elif abs(Ib - Ic) < tol: + rotor = RotorType.PROLATE_TOP + else: + rotor = RotorType.ASYMMETRIC_TOP + + return RotationalConstants( + A=A, + B=B, + C=C, + A_cm=_mhz_to_cm(A), + B_cm=_mhz_to_cm(B), + C_cm=_mhz_to_cm(C), + moments=(float(Ia), float(Ib), float(Ic)), + rotor_type=rotor, + ) + diff --git a/src/opi/utils/constants.py b/src/opi/utils/constants.py new file mode 100644 index 00000000..585359ae --- /dev/null +++ b/src/opi/utils/constants.py @@ -0,0 +1,7 @@ +""" +Contains universal constants. +""" +# Planck's Constant [m^2*kg/s] +H_PLANCK = 6.62607015e-34 +# Speed of light in vacuum [m/s] +C = 2.99792458e10 diff --git a/src/opi/utils/element.py b/src/opi/utils/element.py index a8d00c6d..0d907621 100644 --- a/src/opi/utils/element.py +++ b/src/opi/utils/element.py @@ -900,3 +900,258 @@ def _missing_(cls, value: object, /) -> "Element": Element.OG: 118, Element.OGANESSON: 118, } + +# ============================================================ +# Atomic masses +# ============================================================ +ATOMIC_MASSES_FROM_ELEMENT: dict[Element, float] = { + Element.X: 0.0, + + Element.H: 1.008, + Element.HYDROGEN: 1.008, + Element.HE: 4.003, + Element.HELIUM: 4.003, + Element.LI: 6.941, + Element.LITHIUM: 6.941, + Element.BE: 9.012, + Element.BERYLLIUM: 9.012, + Element.B: 10.810, + Element.BORON: 10.810, + Element.C: 12.011, + Element.CARBON: 12.011, + Element.N: 14.007, + Element.NITROGEN: 14.007, + Element.O: 15.999, + Element.OXYGEN: 15.999, + Element.F: 18.998, + Element.FLUORINE: 18.998, + Element.NE: 20.179, + Element.NEON: 20.179, + + Element.NA: 22.990, + Element.SODIUM: 22.990, + Element.MG: 24.305, + Element.MAGNESIUM: 24.305, + Element.AL: 26.982, + Element.ALUMINUM: 26.982, + Element.SI: 28.086, + Element.SILICON: 28.086, + Element.P: 30.974, + Element.PHOSPHORUS: 30.974, + Element.S: 32.060, + Element.SULFUR: 32.060, + Element.CL: 35.453, + Element.CHLORINE: 35.453, + Element.AR: 39.948, + Element.ARGON: 39.948, + Element.K: 39.100, + Element.POTASSIUM: 39.100, + Element.CA: 40.080, + Element.CALCIUM: 40.080, + + Element.SC: 44.960, + Element.SCANDIUM: 44.960, + Element.TI: 47.900, + Element.TITANIUM: 47.900, + Element.V: 50.940, + Element.VANADIUM: 50.940, + Element.CR: 52.000, + Element.CHROMIUM: 52.000, + Element.MN: 54.940, + Element.MANGANESE: 54.940, + Element.FE: 55.850, + Element.IRON: 55.850, + Element.CO: 58.930, + Element.COBALT: 58.930, + Element.NI: 58.700, + Element.NICKEL: 58.700, + Element.CU: 63.550, + Element.COPPER: 63.550, + Element.ZN: 65.380, + Element.ZINC: 65.380, + + Element.GA: 69.720, + Element.GALLIUM: 69.720, + Element.GE: 72.590, + Element.GERMANIUM: 72.590, + Element.AS: 74.920, + Element.ARSENIC: 74.920, + Element.SE: 78.960, + Element.SELENIUM: 78.960, + Element.BR: 79.900, + Element.BROMINE: 79.900, + Element.KR: 83.800, + Element.KRYPTON: 83.800, + Element.RB: 85.479, + Element.RUBIDIUM: 85.479, + Element.SR: 87.620, + Element.STRONTIUM: 87.620, + Element.Y: 88.910, + Element.YTTRIUM: 88.910, + Element.ZR: 91.220, + Element.ZIRCONIUM: 91.220, + + Element.NB: 92.910, + Element.NIOBIUM: 92.910, + Element.MO: 95.940, + Element.MOLYBDENUM: 95.940, + Element.TC: 97.000, + Element.TECHNETIUM: 97.000, + Element.RU: 101.070, + Element.RUTHENIUM: 101.070, + Element.RH: 102.910, + Element.RHODIUM: 102.910, + Element.PD: 106.400, + Element.PALLADIUM: 106.400, + Element.AG: 107.870, + Element.SILVER: 107.870, + Element.CD: 112.410, + Element.CADMIUM: 112.410, + Element.IN: 114.820, + Element.INDIUM: 114.820, + Element.SN: 118.690, + Element.TIN: 118.690, + + Element.SB: 121.750, + Element.ANTIMONY: 121.750, + Element.TE: 127.600, + Element.TELLURIUM: 127.600, + Element.I: 126.900, + Element.IODINE: 126.900, + Element.XE: 131.300, + Element.XENON: 131.300, + Element.CS: 132.9054, + Element.CESIUM: 132.9054, + Element.BA: 137.3300, + Element.BARIUM: 137.3300, + Element.LA: 138.9055, + Element.LANTHANUM: 138.9055, + Element.CE: 140.1200, + Element.CERIUM: 140.1200, + Element.PR: 140.9077, + Element.PRASEODYMIUM: 140.9077, + Element.ND: 144.2400, + Element.NEODYMIUM: 144.2400, + + Element.PM: 145.0000, + Element.PROMETHIUM: 145.0000, + Element.SM: 150.4000, + Element.SAMARIUM: 150.4000, + Element.EU: 151.9600, + Element.EUROPIUM: 151.9600, + Element.GD: 157.2500, + Element.GADOLINIUM: 157.2500, + Element.TB: 158.9254, + Element.TERBIUM: 158.9254, + Element.DY: 162.5000, + Element.DYSPROSIUM: 162.5000, + Element.HO: 164.9304, + Element.HOLMIUM: 164.9304, + Element.ER: 167.2600, + Element.ERBIUM: 167.2600, + Element.TM: 168.9342, + Element.THULIUM: 168.9342, + Element.YB: 173.0400, + Element.YTTERBIUM: 173.0400, + + Element.LU: 174.9670, + Element.LUTETIUM: 174.9670, + Element.HF: 178.4900, + Element.HAFNIUM: 178.4900, + Element.TA: 180.9479, + Element.TANTALUM: 180.9479, + Element.W: 183.8500, + Element.WOLFRAM: 183.8500, + Element.RE: 186.2070, + Element.RHENIUM: 186.2070, + Element.OS: 190.2000, + Element.OSMIUM: 190.2000, + Element.IR: 192.2200, + Element.IRIDIUM: 192.2200, + Element.PT: 195.0900, + Element.PLATINUM: 195.0900, + Element.AU: 196.9665, + Element.GOLD: 196.9665, + Element.HG: 200.5900, + Element.MERCURY: 200.5900, + + Element.TL: 204.3700, + Element.THALLIUM: 204.3700, + Element.PB: 207.2000, + Element.LEAD: 207.2000, + Element.BI: 208.9804, + Element.BISMUTH: 208.9804, + Element.PO: 209.0000, + Element.POLONIUM: 209.0000, + Element.AT: 210.0000, + Element.ASTATINE: 210.0000, + Element.RN: 222.0000, + Element.RADON: 222.0000, + Element.FR: 223.0000, + Element.FRANCIUM: 223.0000, + Element.RA: 226.0254, + Element.RADIUM: 226.0254, + Element.AC: 227.0278, + Element.ACTINIUM: 227.0278, + Element.TH: 232.0381, + Element.THORIUM: 232.0381, + + Element.PA: 231.0359, + Element.PROTACTINIUM: 231.0359, + Element.U: 238.0290, + Element.URANIUM: 238.0290, + Element.NP: 237.0482, + Element.NEPTUNIUM: 237.0482, + Element.PU: 244.0000, + Element.PLUTONIUM: 244.0000, + Element.AM: 243.0000, + Element.AMERICIUM: 243.0000, + Element.CM: 247.0000, + Element.CURIUM: 247.0000, + Element.BK: 247.0000, + Element.BERKELIUM: 247.0000, + Element.CF: 251.0000, + Element.CALIFORNIUM: 251.0000, + Element.ES: 252.0000, + Element.EINSTEINIUM: 252.0000, + Element.FM: 257.0000, + Element.FERMIUM: 257.0000, + + Element.MD: 258.0000, + Element.MENDELEVIUM: 258.0000, + Element.NO: 259.0000, + Element.NOBELIUM: 259.0000, + Element.LR: 262.0000, + Element.LAWRENCIUM: 262.0000, + Element.RF: 267.0000, + Element.RUTHERFORDIUM: 267.0000, + Element.DB: 268.0000, + Element.DUBNIUM: 268.0000, + Element.SG: 269.0000, + Element.SEABORGIUM: 269.0000, + Element.BH: 270.0000, + Element.BOHRIUM: 270.0000, + Element.HS: 269.0000, + Element.HASSIUM: 269.0000, + Element.MT: 278.0000, + Element.MEITNERIUM: 278.0000, + Element.DS: 281.0000, + Element.DARMSTADTIUM: 281.0000, + + Element.RG: 281.0000, + Element.ROENTGENIUM: 281.0000, + Element.CN: 285.0000, + Element.COPERNICIUM: 285.0000, + Element.NH: 283.0000, + Element.NIHONIUM: 283.0000, + Element.FL: 289.0000, + Element.FLEROVIUM: 289.0000, + Element.MC: 288.0000, + Element.MOSCOVIUM: 288.0000, + Element.LV: 293.0000, + Element.LIVERMORIUM: 293.0000, + Element.TS: 294.0000, + Element.TENNESSINE: 294.0000, + Element.OG: 294.0000, + Element.OGANESSON: 294.0000, +} \ No newline at end of file diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index d7c06cbd..6b5f02f2 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -1,383 +1,69 @@ -import os -import warnings -from dataclasses import dataclass -from typing import List, Tuple, Union +from opi.models.string_enum import StringEnum -import numpy as np +__all__ = ("RotationalConstants", "RotorType",) # ============================================================ -# Physical constants +# Rotor type classification # ============================================================ -_AMU_TO_KG = 1.66053906660e-27 -_ANGSTROM_TO_M = 1.0e-10 -_H_PLANCK = 6.62607015e-34 -_C_CM = 2.99792458e10 +class RotorType(StringEnum): + MONOATOMIC = "monoatomic" + LINEAR = "linear" + SPHERICAL_TOP = "spherical top" + OBLATE_TOP = "symmetric top (oblate)" + PROLATE_TOP = "symmetric top (prolate)" + ASYMMETRIC_TOP = "asymmetric top" # ============================================================ -# Atomic masses +# Rotational constants result container # ============================================================ -ATOMIC_MASSES = { - "X": 0.0, - "PointCharge": 0.0, - "H": 1.008, - "He": 4.003, - "Li": 6.941, - "Be": 9.012, - "B": 10.810, - "C": 12.011, - "N": 14.007, - "O": 15.999, - "F": 18.998, - "Ne": 20.179, - "Na": 22.990, - "Mg": 24.305, - "Al": 26.982, - "Si": 28.086, - "P": 30.974, - "S": 32.060, - "Cl": 35.453, - "Ar": 39.948, - "K": 39.100, - "Ca": 40.080, - "Sc": 44.960, - "Ti": 47.900, - "V": 50.940, - "Cr": 52.000, - "Mn": 54.940, - "Fe": 55.850, - "Co": 58.930, - "Ni": 58.700, - "Cu": 63.550, - "Zn": 65.380, - "Ga": 69.720, - "Ge": 72.590, - "As": 74.920, - "Se": 78.960, - "Br": 79.900, - "Kr": 83.800, - "Rb": 85.479, - "Sr": 87.620, - "Y": 88.910, - "Zr": 91.220, - "Nb": 92.910, - "Mo": 95.940, - "Tc": 97.000, - "Ru": 101.070, - "Rh": 102.910, - "Pd": 106.400, - "Ag": 107.870, - "Cd": 112.410, - "In": 114.820, - "Sn": 118.690, - "Sb": 121.750, - "Te": 127.600, - "I": 126.900, - "Xe": 131.300, - "Cs": 132.9054, - "Ba": 137.3300, - "La": 138.9055, - "Ce": 140.1200, - "Pr": 140.9077, - "Nd": 144.2400, - "Pm": 145.0000, - "Sm": 150.4000, - "Eu": 151.9600, - "Gd": 157.2500, - "Tb": 158.9254, - "Dy": 162.5000, - "Ho": 164.9304, - "Er": 167.2600, - "Tm": 168.9342, - "Yb": 173.0400, - "Lu": 174.9670, - "Hf": 178.4900, - "Ta": 180.9479, - "W": 183.8500, - "Re": 186.2070, - "Os": 190.2000, - "Ir": 192.2200, - "Pt": 195.0900, - "Au": 196.9665, - "Hg": 200.5900, - "Tl": 204.3700, - "Pb": 207.2000, - "Bi": 208.9804, - "Po": 209.0000, - "At": 210.0000, - "Rn": 222.0000, - "Fr": 223.0000, - "Ra": 226.0254, - "Ac": 227.0278, - "Th": 232.0381, - "Pa": 231.0359, - "U": 238.0290, - "Np": 237.0482, - "Pu": 244.0000, - "Am": 243.0000, - "Cm": 247.0000, - "Bk": 247.0000, - "Cf": 251.0000, - "Es": 252.0000, - "Fm": 257.0000, - "Md": 258.0000, - "No": 259.0000, - "Lr": 262.0000, - "Rf": 267.0000, - "Db": 268.0000, - "Sg": 269.0000, - "Bh": 270.0000, - "Hs": 269.0000, - "Mt": 278.0000, - "Ds": 281.0000, - "Rg": 281.0000, - "Cn": 285.0000, - "Nh": 283.0000, - "Fl": 289.0000, - "Mc": 288.0000, - "Lv": 293.0000, - "Ts": 294.0000, - "Og": 294.0000, -} - - -# ============================================================ -# Dataclass -# ============================================================ -@dataclass class RotationalConstants: - A: float | None - B: float | None - C: float | None - A_cm: float | None - B_cm: float | None - C_cm: float | None - moments: tuple[float, float, float] - rotor_type: str + """ + Stores rotational constants and molecular moments of inertia. + + Attributes + ---------- + A, B, C : float | None + Rotational constants in MHz (None for a degenerate axis). + A_cm, B_cm, C_cm : float | None + Rotational constants in cm⁻¹. + moments : tuple[float, float, float] + Principal moments of inertia (Ia, Ib, Ic) in amu·Å². + rotor_type : RotorType + Molecular rotor classification. + """ + + def __init__( + self, + A: float | None, + B: float | None, + C: float | None, + A_cm: float | None, + B_cm: float | None, + C_cm: float | None, + moments: tuple[float, float, float], + rotor_type: RotorType, + ) -> None: + self.A = A + self.B = B + self.C = C + self.A_cm = A_cm + self.B_cm = B_cm + self.C_cm = C_cm + self.moments = moments + self.rotor_type = rotor_type def __str__(self) -> str: def fmt(x: float | None, unit: str = "") -> str: return f"{x:.6f} {unit}" if x is not None else "None" return ( - # "Rotational Spectrum\n" - # "--------------------\n" f"Rotor type : {self.rotor_type}\n\n" "Moments of inertia (amu·Å²):\n" f" Ia = {self.moments[0]:.6f}\n" f" Ib = {self.moments[1]:.6f}\n" f" Ic = {self.moments[2]:.6f}\n\n" - "Rotational constants:\n" + "Rotational constants (yeih):\n" f" A = {fmt(self.A, 'MHz')} ({fmt(self.A_cm, 'cm⁻¹')})\n" f" B = {fmt(self.B, 'MHz')} ({fmt(self.B_cm, 'cm⁻¹')})\n" f" C = {fmt(self.C, 'MHz')} ({fmt(self.C_cm, 'cm⁻¹')})" - ) - - -# ============================================================ -# Utilities -# ============================================================ -def _normalize_symbol(s: str) -> str: - if s in ("X", "PointCharge"): - return s - return s.capitalize() - - -def _read_xyz(data: Union[str, Tuple[List[str], np.ndarray]]) -> Tuple[List[str], np.ndarray]: - """ - Read geometry from: - - XYZ file path - - XYZ block string - - Returns - ------- - symbols : list[str] - coords : (N, 3) ndarray - """ - - # ------------------------- - # Case 1: string input - # ------------------------- - if isinstance(data, str): - # Case 1a: file path - if os.path.isfile(data): - with open(data, "r") as f: - lines = f.readlines() - else: - # Case 1b: XYZ block string - lines = data.strip().splitlines() - - if len(lines) < 3: - raise ValueError("Invalid XYZ format: too few lines") - - try: - n_atoms = int(lines[0].strip()) - except ValueError: - raise ValueError("First line must contain number of atoms") - - symbols = [] - coords = [] - - for i, line in enumerate(lines[2 : 2 + n_atoms], start=3): - parts = line.split() - - if len(parts) < 4: - raise ValueError(f"Line {i} malformed: '{line}'") - - # --- clean symbol --- - raw_sym = parts[0].strip() - - # Remove numeric labels (e.g. C1 → C) - sym = "".join(filter(str.isalpha, raw_sym)) - - # Capitalization (cap insensitive) - sym = sym.capitalize() - - try: - xyz = [float(x) for x in parts[1:4]] - except ValueError: - raise ValueError(f"Invalid coordinates at line {i}: '{line}'") - - symbols.append(sym) - coords.append(xyz) - - return symbols, np.array(coords, dtype=np.float64) - - else: - raise TypeError("Input must be: file path or XYZ string") - - -# ============================================================ -# Main function -# ============================================================ -def rotational_constants( - symbols: list[str] | None = None, - coords: np.ndarray | None = None, - xyz: str | None = None, - masses: np.ndarray | None = None, - weights: dict[str, float] | None = None, - atom_weights: dict[int, float] | None = None, -) -> RotationalConstants | None: - """ - Flexible rotational constant calculator. - - Input options - ------------- - - symbols + coords - - xyz (file path, string, or lines) - - Mass priority - ------------- - masses > atom_weights > weights > default - - Unknown atoms - ------------- - Assigned mass = 0 with warning (unless overridden). - """ - - # --- Input parsing --- - if xyz is not None: - symbols, coords = _read_xyz(xyz) - - if symbols is None or coords is None: - raise ValueError("Provide either (symbols, coords) or xyz input.") - - coords = np.asarray(coords, dtype=np.float64) - - # --- Normalize symbols --- - symbols = [_normalize_symbol(s) for s in symbols] - - # --- Prepare weights --- - weights = {_normalize_symbol(k): v for k, v in (weights or {}).items()} - atom_weights = atom_weights or {} - - # --- Assign masses --- - if masses is not None: - masses = np.asarray(masses, dtype=np.float64) - - else: - masses_list = [] - for i, s in enumerate(symbols): - if i in atom_weights: - m = atom_weights[i] - - elif s in weights: - m = weights[s] - - elif s in ATOMIC_MASSES: - m = ATOMIC_MASSES[s] - - else: - warnings.warn(f"Unknown element '{s}' → mass set to 0.0") - m = 0.0 - - masses_list.append(m) - - masses = np.array(masses_list, dtype=np.float64) - - # --- Filter zero-mass atoms --- - mask = masses > 0.0 - if not np.any(mask): - return None - - masses = masses[mask] - coords = coords[mask] - - total_mass = masses.sum() - - # --- Center of mass --- - com = (masses[:, None] * coords).sum(axis=0) / total_mass - coords -= com - - # --- Inertia tensor --- - inertia = np.zeros((3, 3), dtype=np.float64) - assert coords is not None - for m, r in zip(masses, coords): - inertia += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r)) - - # --- Diagonalize --- - moments_raw, _ = np.linalg.eigh(inertia) - moments_raw = np.maximum(moments_raw, 0.0) - Ia, Ib, Ic = moments_raw - - # --- Convert to rotational constants --- - def _moment_to_mhz(inertia: float) -> float | None: - if inertia < 1e-6: - return None - I_si = inertia * _AMU_TO_KG * (_ANGSTROM_TO_M**2) - return _H_PLANCK / (8.0 * np.pi**2 * I_si) / 1e6 - - def _mhz_to_cm(mhz: float | None) -> float | None: - return None if mhz is None else mhz * 1e6 / _C_CM - - A = _moment_to_mhz(Ia) - B = _moment_to_mhz(Ib) - C = _moment_to_mhz(Ic) - - # --- Rotor classification --- - tol = 1e-3 - n_zero = sum(m < 1e-6 for m in (Ia, Ib, Ic)) - - if n_zero == 3: - rotor = "monoatomic" - elif n_zero == 2: - rotor = "linear" - elif abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: - rotor = "spherical top" - elif abs(Ia - Ib) < tol: - rotor = "symmetric top (oblate)" - elif abs(Ib - Ic) < tol: - rotor = "symmetric top (prolate)" - else: - rotor = "asymmetric top" - - return RotationalConstants( - A=A, - B=B, - C=C, - A_cm=_mhz_to_cm(A), - B_cm=_mhz_to_cm(B), - C_cm=_mhz_to_cm(C), - moments=(Ia, Ib, Ic), - rotor_type=rotor, - ) diff --git a/src/opi/utils/units.py b/src/opi/utils/units.py index 8262558f..4cc8d638 100644 --- a/src/opi/utils/units.py +++ b/src/opi/utils/units.py @@ -10,3 +10,7 @@ AU_TO_KCAL: float = 627.509 # Atomic unit of energy (Hartree) to eV AU_TO_EV: float = 27.2107 +# Angstrom to meter +ANGST_TO_M = 1.0e-10 +# Atomic unit of mass (Dalton) to kilogram (kg) +AMU_TO_KG = 1.66053906660e-27 From 72c82772c3a2b53cccd0fc626c6bdea3c0f912c0 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Wed, 29 Apr 2026 10:35:22 +0200 Subject: [PATCH 05/27] Removed rmsd.py external module --- src/opi/utils/rmsd.py | 108 ------------------------------------------ 1 file changed, 108 deletions(-) delete mode 100644 src/opi/utils/rmsd.py diff --git a/src/opi/utils/rmsd.py b/src/opi/utils/rmsd.py deleted file mode 100644 index 3adf9bdd..00000000 --- a/src/opi/utils/rmsd.py +++ /dev/null @@ -1,108 +0,0 @@ -import os -from typing import List, Tuple, Union - -import numpy as np - - -def read_xyz(data: Union[str, Tuple[List[str], np.ndarray]]) -> Tuple[List[str], np.ndarray]: - """ - Read geometry from: - - XYZ file path - - XYZ block string - - Returns - ------- - symbols : list[str] - coords : (N, 3) ndarray - """ - - if isinstance(data, str): - # Case 1a: a file path - if os.path.isfile(data): - with open(data, "r") as f: - lines = f.readlines() - - # Case 1b: an XYZ block string - else: - lines = data.strip().splitlines() - - n_atoms = int(lines[0].strip()) - symbols = [] - coords = [] - - for line in lines[2 : 2 + n_atoms]: - parts = line.split() - sym = parts[0] - xyz = [float(x) for x in parts[1:4]] - - symbols.append(sym) - coords.append(xyz) - - return symbols, np.array(coords, dtype=np.float64) - - else: - raise TypeError("Input must be: file path or XYZ string") - - -def _validate_geometries(symA: List[str], symB: List[str]) -> None: - if len(symA) != len(symB): - raise ValueError("Geometries have different number of atoms") - - for i, (a, b) in enumerate(zip(symA, symB)): - if a != b: - raise ValueError(f"Atom mismatch at index {i}: {a} != {b}") - - -def kabsch_rmsd( - ref_xyz: str, - target_xyz: str, - *, - align: bool = True, -) -> float: - """Compute RMSD between two XYZ geometries. - - Parameters - ---------- - ref_xyz : str - Reference geometry. - target_xyz : str - Target geometry. - align : bool, default True - Whether to perform optimal alignment (Kabsch). - - Returns - ------- - float - RMSD (Å) - """ - - symA, A = read_xyz(ref_xyz) - symB, B = read_xyz(target_xyz) - - _validate_geometries(symA, symB) - - # Center using centroid (simple average) - A_cent = A - A.mean(axis=0) - B_cent = B - B.mean(axis=0) - - if not align: - diff = A_cent - B_cent - return float(np.sqrt(np.sum(diff**2) / len(A_cent))) - - # Standard Kabsch covariance matrix - H = B_cent.T @ A_cent - - U, _, Vt = np.linalg.svd(H) - - d = np.linalg.det(Vt.T @ U.T) - D = np.diag([1.0, 1.0, d]) - - R = Vt.T @ D @ U.T - - B_rot = B_cent @ R - - diff = A_cent - B_rot - - rmsd = np.sqrt(np.sum(diff**2) / len(A_cent)) - - return float(rmsd) From 3b3644c32a5212db28816a8ba77baa3bd4a26ccc Mon Sep 17 00:00:00 2001 From: crjacinto Date: Wed, 29 Apr 2026 10:55:52 +0200 Subject: [PATCH 06/27] Fixed typos and ambiguous variable names --- src/opi/input/structures/structure.py | 20 ++++++++++---------- src/opi/utils/constants.py | 1 + src/opi/utils/element.py | 14 +------------- src/opi/utils/rotconst.py | 7 ++++++- 4 files changed, 18 insertions(+), 24 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 95bae8fc..53da8497 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -19,9 +19,9 @@ PointCharge, ) from opi.input.structures.coordinates import Coordinates -from opi.utils.element import Element, ATOMIC_MASSES_FROM_ELEMENT -from opi.utils import units, constants -from opi.utils.rotconst import * +from opi.utils import constants, units +from opi.utils.element import ATOMIC_MASSES_FROM_ELEMENT, Element +from opi.utils.rotconst import RotationalConstants, RotorType from opi.utils.tracking_text_io import TrackingTextIO __all__ = ("Structure",) @@ -34,6 +34,7 @@ RGX_FRAG_ID = re.compile(r"(?<=\()\d+(?=\))") RGX_ATOM_SYMBOL_FRAG_ID = re.compile(r"(?P[A-Za-z]{1,2})(\((?P\d+)\))?") + class Structure: """ Class to model internal structure for ORCA calculations. @@ -1231,7 +1232,7 @@ def rmsd_kabsch( # Singular Value Decomposition (SVD) # ------------------------------------------------------------------ # H = U S V^T - # This decomposes the transformation into rotations + scaling + # This decomposes the transformation into rotations + scaling U, _, Vt = np.linalg.svd(H) # ------------------------------------------------------------------ @@ -1253,7 +1254,7 @@ def rotational_constants( masses: npt.NDArray[np.float64] | None = None, weights: dict[str, float] | None = None, atom_weights: dict[int, float] | None = None, - ) ->RotationalConstants | None: + ) -> RotationalConstants | None: """ Compute rotational constants for this structure. @@ -1344,11 +1345,11 @@ def rotational_constants( Ia, Ib, Ic = moments_raw # --- Convert moments (amu·Å²) to rotational constants (MHz) --- - def _moment_to_mhz(I: float) -> float | None: - if I < 1e-6: + def _moment_to_mhz(Inertia: float) -> float | None: + if Inertia < 1e-6: return None - I_si = I * units.AMU_TO_KG * (units.ANGST_TO_M ** 2) - return constants.H_PLANCK / (8.0 * np.pi ** 2 * I_si) / 1e6 + I_si = Inertia * units.AMU_TO_KG * (units.ANGST_TO_M**2) + return constants.H_PLANCK / (8.0 * np.pi**2 * I_si) / 1e6 def _mhz_to_cm(mhz: float | None) -> float | None: # MHz → cm⁻¹ : divide by speed of light in cm/s @@ -1385,4 +1386,3 @@ def _mhz_to_cm(mhz: float | None) -> float | None: moments=(float(Ia), float(Ib), float(Ic)), rotor_type=rotor, ) - diff --git a/src/opi/utils/constants.py b/src/opi/utils/constants.py index 585359ae..033e6bbe 100644 --- a/src/opi/utils/constants.py +++ b/src/opi/utils/constants.py @@ -1,6 +1,7 @@ """ Contains universal constants. """ + # Planck's Constant [m^2*kg/s] H_PLANCK = 6.62607015e-34 # Speed of light in vacuum [m/s] diff --git a/src/opi/utils/element.py b/src/opi/utils/element.py index 0d907621..5983c513 100644 --- a/src/opi/utils/element.py +++ b/src/opi/utils/element.py @@ -906,7 +906,6 @@ def _missing_(cls, value: object, /) -> "Element": # ============================================================ ATOMIC_MASSES_FROM_ELEMENT: dict[Element, float] = { Element.X: 0.0, - Element.H: 1.008, Element.HYDROGEN: 1.008, Element.HE: 4.003, @@ -927,7 +926,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.FLUORINE: 18.998, Element.NE: 20.179, Element.NEON: 20.179, - Element.NA: 22.990, Element.SODIUM: 22.990, Element.MG: 24.305, @@ -948,7 +946,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.POTASSIUM: 39.100, Element.CA: 40.080, Element.CALCIUM: 40.080, - Element.SC: 44.960, Element.SCANDIUM: 44.960, Element.TI: 47.900, @@ -969,7 +966,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.COPPER: 63.550, Element.ZN: 65.380, Element.ZINC: 65.380, - Element.GA: 69.720, Element.GALLIUM: 69.720, Element.GE: 72.590, @@ -990,7 +986,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.YTTRIUM: 88.910, Element.ZR: 91.220, Element.ZIRCONIUM: 91.220, - Element.NB: 92.910, Element.NIOBIUM: 92.910, Element.MO: 95.940, @@ -1011,7 +1006,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.INDIUM: 114.820, Element.SN: 118.690, Element.TIN: 118.690, - Element.SB: 121.750, Element.ANTIMONY: 121.750, Element.TE: 127.600, @@ -1032,7 +1026,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.PRASEODYMIUM: 140.9077, Element.ND: 144.2400, Element.NEODYMIUM: 144.2400, - Element.PM: 145.0000, Element.PROMETHIUM: 145.0000, Element.SM: 150.4000, @@ -1053,7 +1046,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.THULIUM: 168.9342, Element.YB: 173.0400, Element.YTTERBIUM: 173.0400, - Element.LU: 174.9670, Element.LUTETIUM: 174.9670, Element.HF: 178.4900, @@ -1074,7 +1066,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.GOLD: 196.9665, Element.HG: 200.5900, Element.MERCURY: 200.5900, - Element.TL: 204.3700, Element.THALLIUM: 204.3700, Element.PB: 207.2000, @@ -1095,7 +1086,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.ACTINIUM: 227.0278, Element.TH: 232.0381, Element.THORIUM: 232.0381, - Element.PA: 231.0359, Element.PROTACTINIUM: 231.0359, Element.U: 238.0290, @@ -1116,7 +1106,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.EINSTEINIUM: 252.0000, Element.FM: 257.0000, Element.FERMIUM: 257.0000, - Element.MD: 258.0000, Element.MENDELEVIUM: 258.0000, Element.NO: 259.0000, @@ -1137,7 +1126,6 @@ def _missing_(cls, value: object, /) -> "Element": Element.MEITNERIUM: 278.0000, Element.DS: 281.0000, Element.DARMSTADTIUM: 281.0000, - Element.RG: 281.0000, Element.ROENTGENIUM: 281.0000, Element.CN: 285.0000, @@ -1154,4 +1142,4 @@ def _missing_(cls, value: object, /) -> "Element": Element.TENNESSINE: 294.0000, Element.OG: 294.0000, Element.OGANESSON: 294.0000, -} \ No newline at end of file +} diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index 6b5f02f2..700ada4b 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -1,6 +1,10 @@ from opi.models.string_enum import StringEnum -__all__ = ("RotationalConstants", "RotorType",) +__all__ = ( + "RotationalConstants", + "RotorType", +) + # ============================================================ # Rotor type classification @@ -67,3 +71,4 @@ def fmt(x: float | None, unit: str = "") -> str: f" A = {fmt(self.A, 'MHz')} ({fmt(self.A_cm, 'cm⁻¹')})\n" f" B = {fmt(self.B, 'MHz')} ({fmt(self.B_cm, 'cm⁻¹')})\n" f" C = {fmt(self.C, 'MHz')} ({fmt(self.C_cm, 'cm⁻¹')})" + ) From e688f5465ddb49a2ea14bc96aebf261323fd5198 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Wed, 6 May 2026 16:09:10 +0200 Subject: [PATCH 07/27] Splitted RotationalConstatns() function in 3 to get separately Rot. Constatns, Inertia Moments, and Rotor Type --- src/opi/input/structures/structure.py | 179 ++++++++++++++++---------- src/opi/utils/rotconst.py | 115 ++++++++++++++--- 2 files changed, 212 insertions(+), 82 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 53da8497..2981803a 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -19,9 +19,15 @@ PointCharge, ) from opi.input.structures.coordinates import Coordinates -from opi.utils import constants, units from opi.utils.element import ATOMIC_MASSES_FROM_ELEMENT, Element -from opi.utils.rotconst import RotationalConstants, RotorType +from opi.utils.rotconst import ( + PrincipalMoments, + RotationalConstants, + RotorType, + classify_rotor_type, + moment_to_mhz, + mhz_to_wavenumber, +) from opi.utils.tracking_text_io import TrackingTextIO __all__ = ("Structure",) @@ -1249,17 +1255,20 @@ def rmsd_kabsch( return self._rmsd_coords(coords1, coords2_rotated) - def rotational_constants( + # ------------------------------------------------------------------ # + # Moment of inertia # + # ------------------------------------------------------------------ # + def calc_moment_of_inertia( self, masses: npt.NDArray[np.float64] | None = None, weights: dict[str, float] | None = None, atom_weights: dict[int, float] | None = None, - ) -> RotationalConstants | None: + ) -> PrincipalMoments | None: """ - Compute rotational constants for this structure. + Compute the principal axes and moments of inertia for this structure. - Only Atom instances contribute; GhostAtom, PointCharge, and - EmbeddingPotential atoms are silently ignored. + Only ``Atom`` instances contribute; ``GhostAtom``, ``PointCharge``, + and ``EmbeddingPotential`` atoms are silently ignored. Mass priority ------------- @@ -1270,41 +1279,44 @@ def rotational_constants( masses : npt.NDArray[np.float64] | None, default None Per-atom masses (amu) for all Atom instances, overriding every other mass source. - weights : dict[str, float] | None, default None - Per-element mass overrides keyed by element symbol string - (e.g. ``{"C": 13.003}``) atom_weights : dict[int, float] | None, default None Per-atom mass overrides keyed by index within the Atom list. + weights : dict[str, float] | None, default None + Per-element mass overrides keyed by element symbol string + (e.g. ``{"C": 13.003}``). Returns ------- - RotationalConstants | None - None if no Atom instances are present or all masses are zero. + tuple[ndarray shape (3, 3), ndarray shape (3,)] | None + ``(principal_axes, moments)`` where *moments* are in amu·Å² + sorted in ascending order, and each column of *principal_axes* + is the corresponding eigenvector. + Returns ``None`` if no ``Atom`` instances are present or all + masses are zero. Raises ------ ValueError - If the length of *masses* does not match the number of Atom + If the length of *masses* does not match the number of ``Atom`` instances in the structure. """ - # --- Collect Atom instances only --- atom_list = [a for a in self.atoms if isinstance(a, Atom)] - if not atom_list: return None - coords = np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) + coords = np.array( + [a.coordinates.coordinates for a in atom_list], dtype=np.float64 + ) - # --- Prepare weight overrides --- weights = {k: v for k, v in (weights or {}).items()} atom_weights = atom_weights or {} - # --- Assign masses --- if masses is not None: masses = np.asarray(masses, dtype=np.float64) if len(masses) != len(atom_list): raise ValueError( - f"masses length ({len(masses)}) does not match number of atoms ({len(atom_list)})" + f"masses length ({len(masses)}) does not match " + f"number of atoms ({len(atom_list)})" ) else: masses_list: list[float] = [] @@ -1321,68 +1333,103 @@ def rotational_constants( masses_list.append(m) masses = np.array(masses_list, dtype=np.float64) - # --- Filter zero-mass atoms --- mask = masses > 0.0 if not np.any(mask): return None - masses = masses[mask] coords = coords[mask] - total_mass = float(masses.sum()) - # --- Center of mass --- + total_mass = float(masses.sum()) com = (masses[:, None] * coords).sum(axis=0) / total_mass coords -= com - # --- Inertia tensor --- inertia = np.zeros((3, 3), dtype=np.float64) for m, r in zip(masses, coords): inertia += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r)) - # --- Diagonalize: eigenvalues are the principal moments --- - moments_raw, _ = np.linalg.eigh(inertia) + moments_raw, axes = np.linalg.eigh(inertia) # ascending order guaranteed moments_raw = np.maximum(moments_raw, 0.0) - Ia, Ib, Ic = moments_raw - # --- Convert moments (amu·Å²) to rotational constants (MHz) --- - def _moment_to_mhz(Inertia: float) -> float | None: - if Inertia < 1e-6: + return PrincipalMoments(Ia=float(moments_raw[0]), Ib=float(moments_raw[1]), Ic=float(moments_raw[2]), axes=axes) + + # ------------------------------------------------------------------ # + # Rotor classification # + # ------------------------------------------------------------------ # + def calc_rotor_type( + self, + moments: PrincipalMoments | None = None, + **mass_kwargs: Any, + ) -> RotorType | None: + """ + Classify the molecular rotor type. + + Parameters + ---------- + moments : ndarray shape (3,) | None, default None + Pre-computed principal moments (amu·Å², ascending). When + ``None`` the moments are computed via + :meth:`calc_moment_of_inertia` using *mass_kwargs*. + **mass_kwargs + Forwarded to :meth:`calc_moment_of_inertia` when *moments* is + ``None`` (i.e. ``masses``, ``weights``, ``atom_weights``). + + Returns + ------- + RotorType | None + ``None`` if the structure has no real atoms or all masses vanish. + """ + if moments is None: + result = self.calc_moment_of_inertia(**mass_kwargs) + if result is None: return None - I_si = Inertia * units.AMU_TO_KG * (units.ANGST_TO_M**2) - return constants.H_PLANCK / (8.0 * np.pi**2 * I_si) / 1e6 - - def _mhz_to_cm(mhz: float | None) -> float | None: - # MHz → cm⁻¹ : divide by speed of light in cm/s - return None if mhz is None else mhz * 1e6 / constants.C - - A = _moment_to_mhz(Ia) - B = _moment_to_mhz(Ib) - C = _moment_to_mhz(Ic) - - # --- Rotor classification --- - tol = 1e-3 - n_zero = sum(moment < 1e-6 for moment in (Ia, Ib, Ic)) - - if n_zero == 3: - rotor = RotorType.MONOATOMIC - elif n_zero == 2: - rotor = RotorType.LINEAR - elif abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: - rotor = RotorType.SPHERICAL_TOP - elif abs(Ia - Ib) < tol: - rotor = RotorType.OBLATE_TOP - elif abs(Ib - Ic) < tol: - rotor = RotorType.PROLATE_TOP - else: - rotor = RotorType.ASYMMETRIC_TOP + moments = result + + return classify_rotor_type(np.array([moments.Ia, moments.Ib, moments.Ic])) + + # ------------------------------------------------------------------ # + # Rotational constants # + # ------------------------------------------------------------------ # + def calc_rotational_constants( + self, + masses: npt.NDArray[np.float64] | None = None, + weights: dict[str, float] | None = None, + atom_weights: dict[int, float] | None = None, + ) -> RotationalConstants | None: + """ + Compute rotational constants for this structure. + + Only ``Atom`` instances contribute; ``GhostAtom``, ``PointCharge``, + and ``EmbeddingPotential`` atoms are silently ignored. + + Mass priority + ------------- + masses > atom_weights > weights > default (ATOMIC_MASSES_FROM_ELEMENT) + + Parameters + ---------- + masses : npt.NDArray[np.float64] | None, default None + Per-atom masses (amu) for all Atom instances. + atom_weights : dict[int, float] | None, default None + Per-atom mass overrides keyed by index within the Atom list. + weights : dict[str, float] | None, default None + Per-element mass overrides (e.g. ``{"C": 13.003}``). + + Returns + ------- + RotationalConstants | None + ``None`` if no ``Atom`` instances are present or all masses are + zero. + """ + pm = self.calc_moment_of_inertia(masses=masses, weights=weights, atom_weights=atom_weights) + if pm is None: + return None + A = moment_to_mhz(pm.Ia) + B = moment_to_mhz(pm.Ib) + C = moment_to_mhz(pm.Ic) return RotationalConstants( - A=A, - B=B, - C=C, - A_cm=_mhz_to_cm(A), - B_cm=_mhz_to_cm(B), - C_cm=_mhz_to_cm(C), - moments=(float(Ia), float(Ib), float(Ic)), - rotor_type=rotor, - ) + A=A, B=B, C=C, + A_cm=mhz_to_wavenumber(A), + B_cm=mhz_to_wavenumber(B), + C_cm=mhz_to_wavenumber(C), + ) \ No newline at end of file diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index 700ada4b..bcc0edcf 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -1,8 +1,17 @@ +from __future__ import annotations +from dataclasses import dataclass +import numpy as np + from opi.models.string_enum import StringEnum +from opi.utils import constants, units __all__ = ( + "PrincipalMoments", "RotationalConstants", "RotorType", + "classify_rotor_type", + "moment_to_mhz", + "mhz_to_wavenumber", ) @@ -10,13 +19,46 @@ # Rotor type classification # ============================================================ class RotorType(StringEnum): - MONOATOMIC = "monoatomic" - LINEAR = "linear" - SPHERICAL_TOP = "spherical top" - OBLATE_TOP = "symmetric top (oblate)" - PROLATE_TOP = "symmetric top (prolate)" - ASYMMETRIC_TOP = "asymmetric top" + MONOATOMIC = "monoatomic" + LINEAR = "linear" + SPHERICAL_TOP = "spherical top" + OBLATE_TOP = "symmetric top (oblate)" + PROLATE_TOP = "symmetric top (prolate)" + ASYMMETRIC_TOP = "asymmetric top" + + + def __str__(self) -> str: + return f"Rotor type : {self.value}" + +# ============================================================ +# Principal Moments +# ============================================================ + +@dataclass +class PrincipalMoments: + """ + Principal moments of inertia (amu·Å²), sorted ascending. + + Attributes + ---------- + Ia, Ib, Ic : float + Principal moments in amu·Å². + axes : np.ndarray, shape (3, 3) + Corresponding eigenvectors (columns). + """ + Ia: float + Ib: float + Ic: float + axes: np.ndarray + + def __str__(self) -> str: + return ( + "Moments of inertia (amu·Å²):\n" + f" Ia = {self.Ia:.6f}\n" + f" Ib = {self.Ib:.6f}\n" + f" Ic = {self.Ic:.6f}" + ) # ============================================================ # Rotational constants result container @@ -45,8 +87,6 @@ def __init__( A_cm: float | None, B_cm: float | None, C_cm: float | None, - moments: tuple[float, float, float], - rotor_type: RotorType, ) -> None: self.A = A self.B = B @@ -54,21 +94,64 @@ def __init__( self.A_cm = A_cm self.B_cm = B_cm self.C_cm = C_cm - self.moments = moments - self.rotor_type = rotor_type def __str__(self) -> str: def fmt(x: float | None, unit: str = "") -> str: return f"{x:.6f} {unit}" if x is not None else "None" return ( - f"Rotor type : {self.rotor_type}\n\n" - "Moments of inertia (amu·Å²):\n" - f" Ia = {self.moments[0]:.6f}\n" - f" Ib = {self.moments[1]:.6f}\n" - f" Ic = {self.moments[2]:.6f}\n\n" - "Rotational constants (yeih):\n" + "Rotational constants:\n" f" A = {fmt(self.A, 'MHz')} ({fmt(self.A_cm, 'cm⁻¹')})\n" f" B = {fmt(self.B, 'MHz')} ({fmt(self.B_cm, 'cm⁻¹')})\n" f" C = {fmt(self.C, 'MHz')} ({fmt(self.C_cm, 'cm⁻¹')})" ) + +# ============================================================ +# Helper functions (used by Structure methods) +# ============================================================ + +def moment_to_mhz(inertia: float) -> float | None: + """ + Convert a principal moment of inertia (amu·Å²) to a rotational + constant in MHz. Returns None when the moment is effectively zero + (degenerate / linear axis). + """ + if inertia < 1e-6: + return None + I_si = inertia * units.AMU_TO_KG * (units.ANGST_TO_M ** 2) + return constants.H_PLANCK / (8.0 * np.pi ** 2 * I_si) / 1e6 + + +def mhz_to_wavenumber(mhz: float | None) -> float | None: + """Convert a rotational constant from MHz to cm⁻¹.""" + if mhz is None: + return None + return mhz * 1e6 / constants.C + + +def classify_rotor_type(moments: np.ndarray, tol: float = 1e-3) -> RotorType: + """ + Classify the molecular rotor from three principal moments of inertia + (amu·Å², sorted ascending). + + Parameters + ---------- + moments : array-like, shape (3,) + Principal moments Ia ≤ Ib ≤ Ic. + tol : float + Absolute tolerance for treating two moments as equal. + """ + Ia, Ib, Ic = moments + n_zero = sum(m < 1e-6 for m in (Ia, Ib, Ic)) + + if n_zero == 3: + return RotorType.MONOATOMIC + if n_zero == 2: + return RotorType.LINEAR + if abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: + return RotorType.SPHERICAL_TOP + if abs(Ia - Ib) < tol: + return RotorType.OBLATE_TOP + if abs(Ib - Ic) < tol: + return RotorType.PROLATE_TOP + return RotorType.ASYMMETRIC_TOP \ No newline at end of file From 4c1c49838344e0970fd528badbb3271bde4aba37 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Wed, 6 May 2026 16:10:53 +0200 Subject: [PATCH 08/27] Ran nox session --- src/opi/input/structures/structure.py | 19 +++++++++-------- src/opi/utils/rotconst.py | 30 ++++++++++++++++----------- 2 files changed, 28 insertions(+), 21 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 2981803a..5d9f77f8 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -25,8 +25,8 @@ RotationalConstants, RotorType, classify_rotor_type, - moment_to_mhz, mhz_to_wavenumber, + moment_to_mhz, ) from opi.utils.tracking_text_io import TrackingTextIO @@ -1304,9 +1304,7 @@ def calc_moment_of_inertia( if not atom_list: return None - coords = np.array( - [a.coordinates.coordinates for a in atom_list], dtype=np.float64 - ) + coords = np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) weights = {k: v for k, v in (weights or {}).items()} atom_weights = atom_weights or {} @@ -1347,10 +1345,12 @@ def calc_moment_of_inertia( for m, r in zip(masses, coords): inertia += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r)) - moments_raw, axes = np.linalg.eigh(inertia) # ascending order guaranteed + moments_raw, axes = np.linalg.eigh(inertia) # ascending order guaranteed moments_raw = np.maximum(moments_raw, 0.0) - return PrincipalMoments(Ia=float(moments_raw[0]), Ib=float(moments_raw[1]), Ic=float(moments_raw[2]), axes=axes) + return PrincipalMoments( + Ia=float(moments_raw[0]), Ib=float(moments_raw[1]), Ic=float(moments_raw[2]), axes=axes + ) # ------------------------------------------------------------------ # # Rotor classification # @@ -1386,7 +1386,6 @@ def calc_rotor_type( return classify_rotor_type(np.array([moments.Ia, moments.Ib, moments.Ic])) - # ------------------------------------------------------------------ # # Rotational constants # # ------------------------------------------------------------------ # @@ -1428,8 +1427,10 @@ def calc_rotational_constants( B = moment_to_mhz(pm.Ib) C = moment_to_mhz(pm.Ic) return RotationalConstants( - A=A, B=B, C=C, + A=A, + B=B, + C=C, A_cm=mhz_to_wavenumber(A), B_cm=mhz_to_wavenumber(B), C_cm=mhz_to_wavenumber(C), - ) \ No newline at end of file + ) diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index bcc0edcf..03c0f27b 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -1,9 +1,11 @@ from __future__ import annotations + from dataclasses import dataclass + import numpy as np from opi.models.string_enum import StringEnum -from opi.utils import constants, units +from opi.utils import constants, units __all__ = ( "PrincipalMoments", @@ -19,22 +21,22 @@ # Rotor type classification # ============================================================ class RotorType(StringEnum): - MONOATOMIC = "monoatomic" - LINEAR = "linear" - SPHERICAL_TOP = "spherical top" - OBLATE_TOP = "symmetric top (oblate)" - PROLATE_TOP = "symmetric top (prolate)" - ASYMMETRIC_TOP = "asymmetric top" - + MONOATOMIC = "monoatomic" + LINEAR = "linear" + SPHERICAL_TOP = "spherical top" + OBLATE_TOP = "symmetric top (oblate)" + PROLATE_TOP = "symmetric top (prolate)" + ASYMMETRIC_TOP = "asymmetric top" def __str__(self) -> str: return f"Rotor type : {self.value}" # ============================================================ -# Principal Moments +# Principal Moments # ============================================================ + @dataclass class PrincipalMoments: """ @@ -47,6 +49,7 @@ class PrincipalMoments: axes : np.ndarray, shape (3, 3) Corresponding eigenvectors (columns). """ + Ia: float Ib: float Ic: float @@ -60,6 +63,7 @@ def __str__(self) -> str: f" Ic = {self.Ic:.6f}" ) + # ============================================================ # Rotational constants result container # ============================================================ @@ -106,10 +110,12 @@ def fmt(x: float | None, unit: str = "") -> str: f" C = {fmt(self.C, 'MHz')} ({fmt(self.C_cm, 'cm⁻¹')})" ) + # ============================================================ # Helper functions (used by Structure methods) # ============================================================ + def moment_to_mhz(inertia: float) -> float | None: """ Convert a principal moment of inertia (amu·Å²) to a rotational @@ -118,8 +124,8 @@ def moment_to_mhz(inertia: float) -> float | None: """ if inertia < 1e-6: return None - I_si = inertia * units.AMU_TO_KG * (units.ANGST_TO_M ** 2) - return constants.H_PLANCK / (8.0 * np.pi ** 2 * I_si) / 1e6 + I_si = inertia * units.AMU_TO_KG * (units.ANGST_TO_M**2) + return constants.H_PLANCK / (8.0 * np.pi**2 * I_si) / 1e6 def mhz_to_wavenumber(mhz: float | None) -> float | None: @@ -154,4 +160,4 @@ def classify_rotor_type(moments: np.ndarray, tol: float = 1e-3) -> RotorType: return RotorType.OBLATE_TOP if abs(Ib - Ic) < tol: return RotorType.PROLATE_TOP - return RotorType.ASYMMETRIC_TOP \ No newline at end of file + return RotorType.ASYMMETRIC_TOP From 0ba7dc736fc92b002c7ce60daec20cc8c25319f0 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Tue, 19 May 2026 19:42:52 +0200 Subject: [PATCH 09/27] Fixed errors and added unit tests for RMSD and RotConst implementations --- src/opi/input/structures/structure.py | 197 +++++++---- src/opi/utils/rotconst.py | 100 +++--- tests/unit/test_rmsd.py | 430 +++++++++++++++++++++++ tests/unit/test_rotconst.py | 474 ++++++++++++++++++++++++++ 4 files changed, 1080 insertions(+), 121 deletions(-) create mode 100644 tests/unit/test_rmsd.py create mode 100644 tests/unit/test_rotconst.py diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 5d9f77f8..4a7af74e 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1,3 +1,4 @@ +import copy import re from collections.abc import Iterable, Iterator, Sequence from os import PathLike @@ -24,8 +25,6 @@ PrincipalMoments, RotationalConstants, RotorType, - classify_rotor_type, - mhz_to_wavenumber, moment_to_mhz, ) from opi.utils.tracking_text_io import TrackingTextIO @@ -1075,42 +1074,94 @@ def _iter_xyz_structures( if n_struc_limit and n_struc >= n_struc_limit: break - def centered_structure(self) -> "Structure": + def get_coordinates( + self, + atoms: Sequence[Atom | EmbeddingPotential | GhostAtom | PointCharge] | None = None, + ) -> npt.NDArray[np.float64]: """ - Return a new Structure centered at its centroid. + Return the coordinates of the given atoms as a numpy array. + + Parameters + ---------- + atoms : sequence of atom-like objects | None, default None + The atoms to extract coordinates from. If `None`, all entries + in `self.atoms` are used. Returns ------- - Structure - New Structure with centered coordinates. + npt.NDArray[np.float64], shape (N, 3) + Coordinates in the same order as the input list. """ - import copy + atom_list = atoms if atoms is not None else self.atoms + return np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) + + def set_coordinates(self, coords: npt.NDArray[np.float64]) -> "Structure": + """ + Return a new Structure with updated coordinates. + + Parameters + ---------- + coords : npt.NDArray[np.float64], shape (N, 3) + New coordinates for all atoms in the same order as `self.atoms`. + + Returns + ------- + Structure + New `Structure` instance with replaced coordinates; all other + attributes (charge, multiplicity, origin) are preserved. - coords = np.array([a.coordinates.coordinates for a in self.atoms], dtype=np.float64) - centered_coords = coords - coords.mean(axis=0) + Raises + ------ + ValueError + If *coords* shape does not match the number of atoms. + """ + coords = np.asarray(coords, dtype=np.float64) + if coords.shape != (len(self.atoms), 3): + raise ValueError( + f"coords shape {coords.shape} does not match expected ({len(self.atoms)}, 3)" + ) new_structure = copy.deepcopy(self) - for atom, new_coord in zip(new_structure.atoms, centered_coords): + for atom, new_coord in zip(new_structure.atoms, coords): atom.coordinates = new_coord - return new_structure + def centered_structure(self) -> "Structure": + """ + Return a new `Structure` centered at its centroid. + Only `Atom` instances contribute to the centroid calculation; + `GhostAtom`, `PointCharge`, and `EmbeddingPotential` are ignored. + + Returns + ------- + Structure + New `Structure` with centered coordinates. + """ + real_atoms = [a for a in self.atoms if type(a) is Atom] + centroid = self.get_coordinates(atoms=real_atoms).mean(axis=0) + return self.set_coordinates(self.get_coordinates() - centroid) + def _filtered_atoms( self, only_atoms: Sequence[int], ignore_hs: bool, ) -> list[Atom]: """ - Return Atom instances after applying only_atoms and ignore_hs filters. + Return real `Atom` instances after applying `only_atoms` and `ignore_hs` filters. + + The final `type(a) is Atom` check ensures that explicitly indexed atoms + (via `only_atoms`) cannot introduce `GhostAtom`, `PointCharge`, or + `EmbeddingPotential` instances into the returned list. Note that `GhostAtom` + is a subclass of `Atom`, so `isinstance` is intentionally avoided here. """ atom_list = self.atoms if only_atoms: candidates = [atom_list[i] for i in only_atoms] elif ignore_hs: - candidates = [a for a in atom_list if isinstance(a, Atom) and a.element != Element.H] + candidates = [a for a in atom_list if type(a) is Atom and a.element != Element.H] else: - candidates = [a for a in atom_list if isinstance(a, Atom)] - return [a for a in candidates if isinstance(a, Atom)] + candidates = [a for a in atom_list if type(a) is Atom] + return [a for a in candidates if type(a) is Atom] @staticmethod def _validate_rmsd_compatibility(atoms1: list[Atom], atoms2: list[Atom]) -> None: @@ -1121,9 +1172,9 @@ def _validate_rmsd_compatibility(atoms1: list[Atom], atoms2: list[Atom]) -> None raise ValueError( f"Structures have different number of atoms: {len(atoms1)} vs {len(atoms2)}" ) - for i, (a, b) in enumerate(zip(atoms1, atoms2)): + for i, (a, b) in enumerate(zip(atoms1, atoms2), start=1): if a.element != b.element: - raise ValueError(f"Atom mismatch at index {i}: {a.element!r} != {b.element!r}") + raise ValueError(f"Atom mismatch at position {i}: {a.element!r} != {b.element!r}") @staticmethod def _rmsd_coords( @@ -1136,6 +1187,27 @@ def _rmsd_coords( diff = coords1 - coords2 return float(np.sqrt(np.sum(diff**2) / len(coords1))) + def get_coordinates_at_centroid( + self, + atoms: Sequence[Atom | EmbeddingPotential | GhostAtom | PointCharge] | None = None, + ) -> npt.NDArray[np.float64]: + """ + Return coordinates translated to the centroid. + + Parameters + ---------- + atoms : sequence of atom-like objects | None, default None + The atoms to extract coordinates from. If `None`, all entries + in `self.atoms` are used. + + Returns + ------- + npt.NDArray[np.float64], shape (N, 3) + Centroid-translated coordinates in the same order as the input list. + """ + coords = self.get_coordinates(atoms=atoms) + return coords - np.mean(coords, axis=0, dtype=np.float64) + def rmsd( self, other: "Structure", @@ -1153,9 +1225,9 @@ def rmsd( ---------- other : Structure Structure to compare against. - only_atoms : Sequence[int], default () - Atom indices to include. If non-empty, ``ignore_hs`` is ignored. - ignore_hs : bool, default False + only_atoms : Sequence[int], default: () + Atom indices to include. If non-empty, `ignore_hs` is ignored. + ignore_hs : bool, default: False Exclude hydrogen atoms from the RMSD computation. Returns @@ -1172,11 +1244,8 @@ def rmsd( atoms2 = other._filtered_atoms(only_atoms, ignore_hs) self._validate_rmsd_compatibility(atoms1, atoms2) - coords1 = np.array([a.coordinates.coordinates for a in atoms1], dtype=np.float64) - coords2 = np.array([a.coordinates.coordinates for a in atoms2], dtype=np.float64) - - coords1 -= coords1.mean(axis=0) - coords2 -= coords2.mean(axis=0) + coords1 = self.get_coordinates_at_centroid(atoms=atoms1) + coords2 = other.get_coordinates_at_centroid(atoms=atoms2) return self._rmsd_coords(coords1, coords2) @@ -1199,7 +1268,7 @@ def rmsd_kabsch( other : Structure Structure to compare against. only_atoms : Sequence[int], default () - Atom indices to include. If non-empty, ``ignore_hs`` is ignored. + Atom indices to include. If non-empty, `ignore_hs` is ignored. ignore_hs : bool, default False Exclude hydrogen atoms from the RMSD computation. @@ -1217,22 +1286,19 @@ def rmsd_kabsch( atoms2 = other._filtered_atoms(only_atoms, ignore_hs) self._validate_rmsd_compatibility(atoms1, atoms2) - coords1 = np.array([a.coordinates.coordinates for a in atoms1], dtype=np.float64) - coords2 = np.array([a.coordinates.coordinates for a in atoms2], dtype=np.float64) - - # Translate to centroid - coords1 -= coords1.mean(axis=0) - coords2_centered = coords2 - coords2.mean(axis=0) + # Getting centered coordinates + coords1 = self.get_coordinates_at_centroid(atoms=atoms1) + coords2 = other.get_coordinates_at_centroid(atoms=atoms2) # Kabsch algorithm: find optimal rotation matrix via SVD - # https://doi.org/10.1107/S0567739476001873 + # # ------------------------------------------------------------------ # Build covariance matrix # ------------------------------------------------------------------ - # H = B^T A - # This matrix captures how coordinates from B map onto A - H = coords2_centered.T @ coords1 + # H = A^T B + # This matrix captures how coordinates from B (coords2) map onto A (coords1) + H = coords1.T @ coords2 # ------------------------------------------------------------------ # Singular Value Decomposition (SVD) @@ -1251,14 +1317,12 @@ def rmsd_kabsch( R = Vt.T @ D @ U.T - coords2_rotated = coords2_centered @ R - - return self._rmsd_coords(coords1, coords2_rotated) + return self._rmsd_coords(coords1, coords2 @ R) # ------------------------------------------------------------------ # # Moment of inertia # # ------------------------------------------------------------------ # - def calc_moment_of_inertia( + def calc_moments_of_inertia( self, masses: npt.NDArray[np.float64] | None = None, weights: dict[str, float] | None = None, @@ -1267,9 +1331,6 @@ def calc_moment_of_inertia( """ Compute the principal axes and moments of inertia for this structure. - Only ``Atom`` instances contribute; ``GhostAtom``, ``PointCharge``, - and ``EmbeddingPotential`` atoms are silently ignored. - Mass priority ------------- masses > atom_weights > weights > default (ATOMIC_MASSES_FROM_ELEMENT) @@ -1283,29 +1344,32 @@ def calc_moment_of_inertia( Per-atom mass overrides keyed by index within the Atom list. weights : dict[str, float] | None, default None Per-element mass overrides keyed by element symbol string - (e.g. ``{"C": 13.003}``). + (e.g. `{"C": 13.003}`). + + Only `Atom` instances contribute; `GhostAtom`, `PointCharge`, + and `EmbeddingPotential` atoms are silently ignored. Returns ------- - tuple[ndarray shape (3, 3), ndarray shape (3,)] | None - ``(principal_axes, moments)`` where *moments* are in amu·Å² - sorted in ascending order, and each column of *principal_axes* - is the corresponding eigenvector. - Returns ``None`` if no ``Atom`` instances are present or all - masses are zero. + PrincipalMoments | None + Principal axes and moments of inertia in amu·Å², sorted ascending. + Returns `None` if no `Atom` instances are present or all masses are zero. Raises ------ ValueError - If the length of *masses* does not match the number of ``Atom`` + If the length of *masses* does not match the number of `Atom` instances in the structure. """ - atom_list = [a for a in self.atoms if isinstance(a, Atom)] + atom_list = [a for a in self.atoms if type(a) is Atom] if not atom_list: return None - coords = np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) + # Extract coordinates from real atoms only + coords = self.get_coordinates(atoms=atom_list) + # --- Mass assignment priority: masses > atom_weights > weights > default --- + # Normalise the optional dicts to avoid None checks below weights = {k: v for k, v in (weights or {}).items()} atom_weights = atom_weights or {} @@ -1334,10 +1398,12 @@ def calc_moment_of_inertia( mask = masses > 0.0 if not np.any(mask): return None + # --- Filter out zero-mass atoms before computing the inertia tensor --- + # This handles dummy atoms or unknown elements that were assigned 0.0 masses = masses[mask] coords = coords[mask] - total_mass = float(masses.sum()) + com = (masses[:, None] * coords).sum(axis=0) / total_mass coords -= com @@ -1379,12 +1445,12 @@ def calc_rotor_type( ``None`` if the structure has no real atoms or all masses vanish. """ if moments is None: - result = self.calc_moment_of_inertia(**mass_kwargs) + result = self.calc_moments_of_inertia(**mass_kwargs) if result is None: return None moments = result - return classify_rotor_type(np.array([moments.Ia, moments.Ib, moments.Ic])) + return moments.rotor_type() # ------------------------------------------------------------------ # # Rotational constants # @@ -1398,8 +1464,8 @@ def calc_rotational_constants( """ Compute rotational constants for this structure. - Only ``Atom`` instances contribute; ``GhostAtom``, ``PointCharge``, - and ``EmbeddingPotential`` atoms are silently ignored. + Only `Atom`` instances contribute; `GhostAtom`, `PointCharge`, + and `EmbeddingPotential` atoms are silently ignored. Mass priority ------------- @@ -1412,25 +1478,18 @@ def calc_rotational_constants( atom_weights : dict[int, float] | None, default None Per-atom mass overrides keyed by index within the Atom list. weights : dict[str, float] | None, default None - Per-element mass overrides (e.g. ``{"C": 13.003}``). + Per-element mass overrides (e.g. `{"C": 13.003}`). Returns ------- RotationalConstants | None - ``None`` if no ``Atom`` instances are present or all masses are + `None` if no `Atom` instances are present or all masses are zero. """ - pm = self.calc_moment_of_inertia(masses=masses, weights=weights, atom_weights=atom_weights) + pm = self.calc_moments_of_inertia(masses=masses, weights=weights, atom_weights=atom_weights) if pm is None: return None A = moment_to_mhz(pm.Ia) B = moment_to_mhz(pm.Ib) C = moment_to_mhz(pm.Ic) - return RotationalConstants( - A=A, - B=B, - C=C, - A_cm=mhz_to_wavenumber(A), - B_cm=mhz_to_wavenumber(B), - C_cm=mhz_to_wavenumber(C), - ) + return RotationalConstants(A=A, B=B, C=C) diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index 03c0f27b..05575330 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -11,7 +11,6 @@ "PrincipalMoments", "RotationalConstants", "RotorType", - "classify_rotor_type", "moment_to_mhz", "mhz_to_wavenumber", ) @@ -55,6 +54,22 @@ class PrincipalMoments: Ic: float axes: np.ndarray + def rotor_type(self, tol: float = 1e-3) -> RotorType: + Ia, Ib, Ic = self.Ia, self.Ib, self.Ic + n_zero = sum(m < tol for m in (Ia, Ib, Ic)) + + if n_zero == 3: + return RotorType.MONOATOMIC + if n_zero == 1 and abs(Ib - Ic) < tol: + return RotorType.LINEAR + if abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: + return RotorType.SPHERICAL_TOP + if abs(Ib - Ic) < tol: # Ia < Ib == Ic → prolate + return RotorType.PROLATE_TOP + if abs(Ia - Ib) < tol: # Ia == Ib < Ic → oblate + return RotorType.OBLATE_TOP + return RotorType.ASYMMETRIC_TOP + def __str__(self) -> str: return ( "Moments of inertia (amu·Å²):\n" @@ -67,20 +82,16 @@ def __str__(self) -> str: # ============================================================ # Rotational constants result container # ============================================================ + + class RotationalConstants: """ - Stores rotational constants and molecular moments of inertia. + Stores rotational constants in MHz. Attributes ---------- A, B, C : float | None - Rotational constants in MHz (None for a degenerate axis). - A_cm, B_cm, C_cm : float | None - Rotational constants in cm⁻¹. - moments : tuple[float, float, float] - Principal moments of inertia (Ia, Ib, Ic) in amu·Å². - rotor_type : RotorType - Molecular rotor classification. + Rotational constants in MHz (`None` for a degenerate axis). """ def __init__( @@ -88,26 +99,36 @@ def __init__( A: float | None, B: float | None, C: float | None, - A_cm: float | None, - B_cm: float | None, - C_cm: float | None, ) -> None: self.A = A self.B = B self.C = C - self.A_cm = A_cm - self.B_cm = B_cm - self.C_cm = C_cm + + def get_in_wavenumbers(self) -> tuple[float | None, float | None, float | None]: + """ + Return rotational constants A, B, C converted to cm⁻¹. + + Returns + ------- + tuple[float | None, float | None, float | None] + A, B, C in cm⁻¹. `None` for any degenerate axis. + """ + return ( + mhz_to_wavenumber(self.A), + mhz_to_wavenumber(self.B), + mhz_to_wavenumber(self.C), + ) def __str__(self) -> str: def fmt(x: float | None, unit: str = "") -> str: return f"{x:.6f} {unit}" if x is not None else "None" + A_cm, B_cm, C_cm = self.get_in_wavenumbers() return ( - "Rotational constants:\n" - f" A = {fmt(self.A, 'MHz')} ({fmt(self.A_cm, 'cm⁻¹')})\n" - f" B = {fmt(self.B, 'MHz')} ({fmt(self.B_cm, 'cm⁻¹')})\n" - f" C = {fmt(self.C, 'MHz')} ({fmt(self.C_cm, 'cm⁻¹')})" + "Rotational constants (yeih):\n" + f" A = {fmt(self.A, 'MHz')} ({fmt(A_cm, 'cm⁻¹')})\n" + f" B = {fmt(self.B, 'MHz')} ({fmt(B_cm, 'cm⁻¹')})\n" + f" C = {fmt(self.C, 'MHz')} ({fmt(C_cm, 'cm⁻¹')})" ) @@ -116,48 +137,23 @@ def fmt(x: float | None, unit: str = "") -> str: # ============================================================ -def moment_to_mhz(inertia: float) -> float | None: +def moment_to_mhz(inertia: float | None) -> float | None: """ Convert a principal moment of inertia (amu·Å²) to a rotational - constant in MHz. Returns None when the moment is effectively zero - (degenerate / linear axis). + constant in MHz. Returns `None` when the moment is `None` or + effectively zero (degenerate / linear axis). """ - if inertia < 1e-6: + if inertia is None or inertia < 1e-6: return None I_si = inertia * units.AMU_TO_KG * (units.ANGST_TO_M**2) return constants.H_PLANCK / (8.0 * np.pi**2 * I_si) / 1e6 def mhz_to_wavenumber(mhz: float | None) -> float | None: - """Convert a rotational constant from MHz to cm⁻¹.""" + """ + Convert a rotational constant from MHz to cm⁻¹. + Returns `None` if `mhz` is `None`. + """ if mhz is None: return None return mhz * 1e6 / constants.C - - -def classify_rotor_type(moments: np.ndarray, tol: float = 1e-3) -> RotorType: - """ - Classify the molecular rotor from three principal moments of inertia - (amu·Å², sorted ascending). - - Parameters - ---------- - moments : array-like, shape (3,) - Principal moments Ia ≤ Ib ≤ Ic. - tol : float - Absolute tolerance for treating two moments as equal. - """ - Ia, Ib, Ic = moments - n_zero = sum(m < 1e-6 for m in (Ia, Ib, Ic)) - - if n_zero == 3: - return RotorType.MONOATOMIC - if n_zero == 2: - return RotorType.LINEAR - if abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: - return RotorType.SPHERICAL_TOP - if abs(Ia - Ib) < tol: - return RotorType.OBLATE_TOP - if abs(Ib - Ic) < tol: - return RotorType.PROLATE_TOP - return RotorType.ASYMMETRIC_TOP diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_rmsd.py new file mode 100644 index 00000000..7eb1093b --- /dev/null +++ b/tests/unit/test_rmsd.py @@ -0,0 +1,430 @@ +""" +Unit tests for RMSD and coordinate utilities. + +Covers: +- get_coordinates() +- get_coordinates_at_centroid() +- set_coordinates() +- centered_structure() +- _filtered_atoms() +- _validate_rmsd_compatibility() +- rmsd() +- rmsd_kabsch() + +Edge cases: ignore_hs, only_atoms, mismatched structures, +mixed atom types, wrong shapes, identical structures. +""" + +import numpy as np +import pytest + +from opi.input.structures.atom import Atom, PointCharge +from opi.input.structures.coordinates import Coordinates +from opi.input.structures.structure import Structure +from opi.utils.element import Element + +# ============================================================ +# Helpers +# ============================================================ + + +def make_atom(element: str, x: float, y: float, z: float) -> Atom: + return Atom( + element=Element(element), + coordinates=Coordinates(coordinates=(x, y, z)), + ) + + +def make_structure(*atoms: Atom) -> Structure: + return Structure(atoms=list(atoms)) + + +# ============================================================ +# Fixtures +# ============================================================ + + +@pytest.fixture() +def water() -> Structure: + """H2O at a known geometry.""" + return make_structure( + make_atom("O", 0.000000, 0.000000, 0.119748), + make_atom("H", 0.000000, 0.756950, -0.478993), + make_atom("H", 0.000000, -0.756950, -0.478993), + ) + + +@pytest.fixture() +def water_translated() -> Structure: + """H2O shifted by (1, 2, 3) — RMSD vs water should be 0 after centring.""" + return make_structure( + make_atom("O", 1.000000, 2.000000, 3.119748), + make_atom("H", 1.000000, 2.756950, 2.521007), + make_atom("H", 1.000000, 1.243050, 2.521007), + ) + + +@pytest.fixture() +def water_rotated(water) -> Structure: + """ + H2O rotated 90° around Z after centring. + Kabsch RMSD vs the original centred water should be ~0. + """ + centered = water.centered_structure() + coords = centered.get_coordinates() + theta = np.pi / 2 + R = np.array( + [ + [np.cos(theta), -np.sin(theta), 0], + [np.sin(theta), np.cos(theta), 0], + [0, 0, 1], + ] + ) + return centered.set_coordinates(coords @ R.T) + + +@pytest.fixture() +def ethanol() -> Structure: + """ + Simple ethanol-like structure (C2H5OH) for ignore_hs tests. + Coordinates are approximate — correctness of geometry is not critical here. + """ + return make_structure( + make_atom("C", 0.000, 0.000, 0.000), + make_atom("C", 1.540, 0.000, 0.000), + make_atom("O", 2.060, 1.190, 0.000), + make_atom("H", -0.390, 1.020, 0.000), + make_atom("H", -0.390, -0.510, 0.890), + make_atom("H", -0.390, -0.510, -0.890), + make_atom("H", 1.930, -0.510, 0.890), + make_atom("H", 1.930, -0.510, -0.890), + make_atom("H", 2.980, 1.190, 0.000), + ) + + +@pytest.fixture() +def mixed_structure() -> Structure: + """Structure with a real Atom and a PointCharge (not an Atom subclass).""" + return Structure( + atoms=[ + Atom(element=Element("O"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), + PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), + ] + ) + + +# ============================================================ +# get_coordinates +# ============================================================ + + +class TestGetCoordinates: + def test_shape(self, water): + coords = water.get_coordinates() + assert coords.shape == (3, 3) + + def test_dtype(self, water): + coords = water.get_coordinates() + assert coords.dtype == np.float64 + + def test_values(self, water): + coords = water.get_coordinates() + np.testing.assert_allclose(coords[0], [0.0, 0.0, 0.119748]) + + def test_atoms_kwarg_filters(self, mixed_structure): + """Passing only real Atom instances should exclude PointCharge.""" + real = [a for a in mixed_structure.atoms if type(a) is Atom] + coords = mixed_structure.get_coordinates(atoms=real) + assert coords.shape == (1, 3) + + def test_default_includes_all(self, mixed_structure): + """Default call includes all atom types.""" + coords = mixed_structure.get_coordinates() + assert coords.shape == (2, 3) + + def test_single_atom(self): + s = make_structure(make_atom("C", 1.0, 2.0, 3.0)) + coords = s.get_coordinates() + np.testing.assert_allclose(coords[0], [1.0, 2.0, 3.0]) + + +# ============================================================ +# get_coordinates_at_centroid +# ============================================================ + + +class TestGetCoordinatesAtCentroid: + def test_centroid_is_zero(self, water): + coords = water.get_coordinates_at_centroid() + np.testing.assert_allclose(coords.mean(axis=0), [0.0, 0.0, 0.0], atol=1e-12) + + def test_shape_preserved(self, water): + coords = water.get_coordinates_at_centroid() + assert coords.shape == (3, 3) + + def test_with_atoms_kwarg(self, water): + atoms = water.atoms[:2] + coords = water.get_coordinates_at_centroid(atoms=atoms) + assert coords.shape == (2, 3) + np.testing.assert_allclose(coords.mean(axis=0), [0.0, 0.0, 0.0], atol=1e-12) + + def test_translation_invariant(self, water, water_translated): + c1 = water.get_coordinates_at_centroid() + c2 = water_translated.get_coordinates_at_centroid() + np.testing.assert_allclose(c1, c2, atol=1e-6) + + +# ============================================================ +# set_coordinates +# ============================================================ + + +class TestSetCoordinates: + def test_returns_new_structure(self, water): + new_coords = water.get_coordinates() + 1.0 + new_s = water.set_coordinates(new_coords) + assert new_s is not water + + def test_original_unchanged(self, water): + original_coords = water.get_coordinates().copy() + water.set_coordinates(water.get_coordinates() + 5.0) + np.testing.assert_allclose(water.get_coordinates(), original_coords) + + def test_new_coords_applied(self, water): + new_coords = np.zeros((3, 3)) + new_s = water.set_coordinates(new_coords) + np.testing.assert_allclose(new_s.get_coordinates(), new_coords) + + def test_raises_on_wrong_shape(self, water): + with pytest.raises(ValueError, match="coords shape"): + water.set_coordinates(np.zeros((2, 3))) + + def test_charge_preserved(self, water): + water_charged = Structure(atoms=water.atoms, charge=1) + new_s = water_charged.set_coordinates(water.get_coordinates()) + assert new_s.charge == 1 + + def test_multiplicity_preserved(self, water): + water_mult = Structure(atoms=water.atoms, multiplicity=3) + new_s = water_mult.set_coordinates(water.get_coordinates()) + assert new_s.multiplicity == 3 + + +# ============================================================ +# centered_structure +# ============================================================ + + +class TestCenteredStructure: + def test_centroid_at_origin(self, water): + centered = water.centered_structure() + real_atoms = [a for a in centered.atoms if isinstance(a, Atom)] + coords = centered.get_coordinates(atoms=real_atoms) + np.testing.assert_allclose(coords.mean(axis=0), [0.0, 0.0, 0.0], atol=1e-12) + + def test_returns_new_structure(self, water): + assert water.centered_structure() is not water + + def test_original_unchanged(self, water): + original = water.get_coordinates().copy() + water.centered_structure() + np.testing.assert_allclose(water.get_coordinates(), original) + + def test_pointcharge_excluded_from_centroid(self, mixed_structure): + """ + The PointCharge at (1, 0, 0) should not pull the centroid away + from the real Atom at (0, 0, 0). + """ + centered = mixed_structure.centered_structure() + real_atoms = [a for a in centered.atoms if type(a) is Atom] + coords = centered.get_coordinates(atoms=real_atoms) + # Real atom was already at origin, so it should stay at origin + np.testing.assert_allclose(coords[0], [0.0, 0.0, 0.0], atol=1e-12) + + def test_translated_centered_equals_original_centered(self, water, water_translated): + c1 = water.centered_structure().get_coordinates() + c2 = water_translated.centered_structure().get_coordinates() + np.testing.assert_allclose(c1, c2, atol=1e-6) + + +# ============================================================ +# _filtered_atoms +# ============================================================ + + +class TestFilteredAtoms: + def test_default_returns_only_atoms(self, mixed_structure): + """PointCharge should be excluded since it is not an Atom instance.""" + filtered = mixed_structure._filtered_atoms((), False) + assert all(isinstance(a, Atom) for a in filtered) + assert len(filtered) == 1 + + def test_ignore_hs(self, ethanol): + filtered = ethanol._filtered_atoms((), True) + elements = [a.element for a in filtered] + assert Element("H") not in elements + + def test_only_atoms_indices(self, ethanol): + filtered = ethanol._filtered_atoms([0, 1, 2], False) + assert len(filtered) == 3 + + def test_only_atoms_excludes_pointcharge(self): + """Explicitly indexed PointCharge must be excluded by the final isinstance check.""" + structure = Structure( + atoms=[ + Atom(element=Element("C"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), + PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), + ] + ) + filtered = structure._filtered_atoms([0, 1], False) + assert all(isinstance(a, Atom) for a in filtered) + assert len(filtered) == 1 + + def test_only_atoms_takes_priority_over_ignore_hs(self, ethanol): + """When only_atoms is set, ignore_hs should be ignored.""" + filtered_with = ethanol._filtered_atoms([3, 4], True) + filtered_without = ethanol._filtered_atoms([3, 4], False) + assert len(filtered_with) == len(filtered_without) + + +# ============================================================ +# _validate_rmsd_compatibility +# ============================================================ + + +class TestValidateRmsdCompatibility: + def test_compatible_structures(self, water): + atoms = [a for a in water.atoms if isinstance(a, Atom)] + Structure._validate_rmsd_compatibility(atoms, atoms) # should not raise + + def test_raises_on_different_count(self, water, ethanol): + a1 = [a for a in water.atoms if isinstance(a, Atom)] + a2 = [a for a in ethanol.atoms if isinstance(a, Atom)] + with pytest.raises(ValueError, match="different number of atoms"): + Structure._validate_rmsd_compatibility(a1, a2) + + def test_raises_on_element_mismatch(self): + a1 = [make_atom("C", 0, 0, 0), make_atom("H", 1, 0, 0)] + a2 = [make_atom("C", 0, 0, 0), make_atom("O", 1, 0, 0)] + with pytest.raises(ValueError, match="position 2"): + Structure._validate_rmsd_compatibility(a1, a2) + + def test_error_message_uses_natural_counting(self): + """Position in error message should start at 1, not 0.""" + a1 = [make_atom("C", 0, 0, 0)] + a2 = [make_atom("O", 0, 0, 0)] + with pytest.raises(ValueError, match="position 1"): + Structure._validate_rmsd_compatibility(a1, a2) + + +# ============================================================ +# rmsd +# ============================================================ + + +class TestRmsd: + def test_identical_structures_zero_rmsd(self, water): + assert pytest.approx(water.rmsd(water), abs=1e-10) == 0.0 + + def test_translated_structure_zero_rmsd(self, water, water_translated): + """Pure translation should give 0 RMSD after centring.""" + assert pytest.approx(water.rmsd(water_translated), abs=1e-6) == 0.0 + + def test_symmetry(self, water, water_translated): + assert pytest.approx(water.rmsd(water_translated)) == water_translated.rmsd(water) + + def test_raises_on_incompatible_structures(self, water, ethanol): + with pytest.raises(ValueError): + water.rmsd(ethanol) + + def test_ignore_hs(self, ethanol): + """ + Displace only heavy atoms — RMSD with and without H should differ + since the centroid shift differs between the two filtered sets. + """ + coords = ethanol.get_coordinates() + new_coords = coords.copy() + new_coords[0] += np.array([0.5, 0.0, 0.0]) # shift only C0 + shifted = ethanol.set_coordinates(new_coords) + rmsd_all = ethanol.rmsd(shifted) + rmsd_no_h = ethanol.rmsd(shifted, ignore_hs=True) + assert rmsd_all != pytest.approx(rmsd_no_h, abs=1e-6) + + def test_only_atoms_subset(self, ethanol): + """RMSD over a subset of atoms should differ from the full-molecule RMSD.""" + coords = ethanol.get_coordinates() + new_coords = coords.copy() + new_coords[0] += np.array([0.5, 0.0, 0.0]) # shift only C0 + shifted = ethanol.set_coordinates(new_coords) + rmsd_all = ethanol.rmsd(shifted) + rmsd_subset = ethanol.rmsd(shifted, only_atoms=[0, 1, 2]) + assert rmsd_all != pytest.approx(rmsd_subset, abs=1e-6) + + def test_nonzero_rmsd_for_different_structures(self, water): + other = water.set_coordinates(water.get_coordinates() + np.array([0, 0, 1.0])) + # After centring the z-shift cancels, but individual atom positions differ + # We just confirm it's finite and non-negative + result = water.rmsd(other) + assert result >= 0.0 + + def test_result_is_float(self, water): + assert isinstance(water.rmsd(water), float) + + +# ============================================================ +# rmsd_kabsch +# ============================================================ + + +class TestRmsdKabsch: + def test_identical_structures_zero_rmsd(self, water): + assert pytest.approx(water.rmsd_kabsch(water), abs=1e-10) == 0.0 + + def test_translated_zero_rmsd(self, water, water_translated): + assert pytest.approx(water.rmsd_kabsch(water_translated), abs=1e-6) == 0.0 + + def test_rotated_zero_rmsd(self, water, water_rotated): + """Kabsch should align the rotation and give ~0 RMSD.""" + centered = water.centered_structure() + assert pytest.approx(centered.rmsd_kabsch(water_rotated), abs=1e-6) == 0.0 + + def test_kabsch_le_rmsd(self, water): + """Kabsch RMSD ≤ plain RMSD (optimal rotation can only help or be neutral).""" + centered = water.centered_structure() + coords = centered.get_coordinates() + R = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]], dtype=float) + other = centered.set_coordinates(coords @ R.T) + assert centered.rmsd_kabsch(other) <= centered.rmsd(other) + 1e-10 + + def test_symmetry(self, water, water_rotated): + assert pytest.approx( + water.rmsd_kabsch(water_rotated), abs=1e-6 + ) == water_rotated.rmsd_kabsch(water) + + def test_raises_on_incompatible_structures(self, water, ethanol): + with pytest.raises(ValueError): + water.rmsd_kabsch(ethanol) + + def test_ignore_hs(self, ethanol): + coords = ethanol.get_coordinates() + new_coords = coords.copy() + new_coords[0] += np.array([0.5, 0.0, 0.0]) + shifted = ethanol.set_coordinates(new_coords) + rmsd_all = ethanol.rmsd_kabsch(shifted) + rmsd_no_h = ethanol.rmsd_kabsch(shifted, ignore_hs=True) + assert rmsd_all != pytest.approx(rmsd_no_h, abs=1e-6) + + def test_only_atoms_subset(self, ethanol): + coords = ethanol.get_coordinates() + new_coords = coords.copy() + new_coords[0] += np.array([0.5, 0.0, 0.0]) + shifted = ethanol.set_coordinates(new_coords) + rmsd_all = ethanol.rmsd_kabsch(shifted) + rmsd_subset = ethanol.rmsd_kabsch(shifted, only_atoms=[0, 1, 2]) + assert rmsd_all != pytest.approx(rmsd_subset, abs=1e-6) + + def test_result_is_float(self, water): + assert isinstance(water.rmsd_kabsch(water), float) + + def test_result_nonnegative(self, water, water_rotated): + assert water.rmsd_kabsch(water_rotated) >= 0.0 diff --git a/tests/unit/test_rotconst.py b/tests/unit/test_rotconst.py new file mode 100644 index 00000000..657c38eb --- /dev/null +++ b/tests/unit/test_rotconst.py @@ -0,0 +1,474 @@ +""" +Unit tests for rotational constants functionality. + +Covers: +- calc_moments_of_inertia() +- calc_rotor_type() +- calc_rotational_constants() +- PrincipalMoments.rotor_type() +- PrincipalMoments.__str__() +- RotationalConstants.get_in_wavenumbers() +- RotationalConstants.__str__() +- moment_to_mhz() and mhz_to_wavenumber() +- All RotorType classifications +- Edge cases: empty structure, all-zero masses, mass overrides, unknown elements +""" + +import math + +import numpy as np +import pytest + +from opi.input.structures.atom import Atom, PointCharge +from opi.input.structures.coordinates import Coordinates +from opi.input.structures.structure import Structure +from opi.utils.element import Element +from opi.utils.rotconst import ( + RotationalConstants, + RotorType, + mhz_to_wavenumber, + moment_to_mhz, +) + +# ============================================================ +# Helpers +# ============================================================ + + +def make_atom(element: str, x: float, y: float, z: float) -> Atom: + return Atom( + element=Element(element), + coordinates=Coordinates(coordinates=(x, y, z)), + ) + + +def make_no_real_atoms_structure() -> Structure: + """Structure containing only a PointCharge — no real Atom instances.""" + return Structure( + atoms=[PointCharge(coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0)), charge=1.0)] + ) + + +def make_structure(*atoms: Atom) -> Structure: + return Structure(atoms=list(atoms)) + + +# ============================================================ +# Fixtures +# ============================================================ + + +@pytest.fixture() +def water() -> Structure: + """H2O — asymmetric top.""" + return make_structure( + make_atom("O", 0.000000, 0.000000, 0.119748), + make_atom("H", 0.000000, 0.756950, -0.478993), + make_atom("H", 0.000000, -0.756950, -0.478993), + ) + + +@pytest.fixture() +def co2() -> Structure: + """ + CO2 — linear molecule aligned exactly on Z. + After COM shift, both O atoms lie on the Z axis so Ia = 0 exactly. + Using exact masses (C=12, O=15.999) the COM is at z=0 by symmetry. + """ + return make_structure( + make_atom("C", 0.000000, 0.000000, 0.000000), + make_atom("O", 0.000000, 0.000000, 1.160000), + make_atom("O", 0.000000, 0.000000, -1.160000), + ) + + +@pytest.fixture() +def hcl() -> Structure: + """HCl — guaranteed linear (diatomic), Ia = 0 exactly on Z axis.""" + return make_structure( + make_atom("H", 0.000000, 0.000000, 0.000000), + make_atom("Cl", 0.000000, 0.000000, 1.274500), + ) + + +@pytest.fixture() +def methane() -> Structure: + """CH4 — spherical top.""" + d = 0.6276 # C-H bond / sqrt(3) + return make_structure( + make_atom("C", 0.000, 0.000, 0.000), + make_atom("H", d, d, d), + make_atom("H", -d, -d, d), + make_atom("H", -d, d, -d), + make_atom("H", d, -d, -d), + ) + + +@pytest.fixture() +def benzene() -> Structure: + """C6H6 — oblate symmetric top.""" + r_c, r_h = 1.3970, 2.4832 + atoms = [] + for i in range(6): + angle = math.radians(i * 60) + atoms.append(make_atom("C", r_c * math.cos(angle), r_c * math.sin(angle), 0.0)) + atoms.append(make_atom("H", r_h * math.cos(angle), r_h * math.sin(angle), 0.0)) + return make_structure(*atoms) + + +@pytest.fixture() +def single_atom() -> Structure: + """Single atom — monoatomic.""" + return make_structure(make_atom("C", 0.0, 0.0, 0.0)) + + +@pytest.fixture() +def ammonia() -> Structure: + """NH3 — oblate symmetric top (C3v). XTB2 optimized geometry.""" + return make_structure( + make_atom("N", -0.000000, -0.000066, 0.100407), + make_atom("H", 0.000000, 0.943825, -0.266468), + make_atom("H", 0.817493, -0.471879, -0.266439), + make_atom("H", -0.817493, -0.471879, -0.266439), + ) + + +@pytest.fixture() +def ch3cl() -> Structure: + """CH3Cl — prolate symmetric top (C3v). XTB2 optimized geometry.""" + return make_structure( + make_atom("Cl", 0.9754830000, 0.0921220000, -0.0239260000), + make_atom("C", 2.7424550000, 0.0921300000, -0.0239350000), + make_atom("H", 3.0994960000, 0.3238830000, 0.9818970000), + make_atom("H", 3.0994920000, -0.8948230000, -0.3261570000), + make_atom("H", 3.0994750000, 0.8473380000, -0.7275610000), + ) + + +# ============================================================ +# moment_to_mhz and mhz_to_wavenumber +# ============================================================ + + +class TestHelperFunctions: + def test_moment_to_mhz_none_input(self): + assert moment_to_mhz(None) is None + + def test_moment_to_mhz_zero(self): + assert moment_to_mhz(0.0) is None + + def test_moment_to_mhz_below_threshold(self): + assert moment_to_mhz(1e-7) is None + + def test_moment_to_mhz_positive(self): + result = moment_to_mhz(100.0) + assert result is not None + assert result > 0.0 + + def test_moment_to_mhz_known_value(self): + # For water Ib ≈ 1.022 amu·Å² → B ≈ 435000 MHz (rough check) + result = moment_to_mhz(1.022) + assert result is not None + assert 400_000 < result < 500_000 + + def test_mhz_to_wavenumber_none(self): + assert mhz_to_wavenumber(None) is None + + def test_mhz_to_wavenumber_positive(self): + result = mhz_to_wavenumber(100_000.0) + assert result is not None + assert result > 0.0 + + def test_mhz_to_wavenumber_proportional(self): + a = mhz_to_wavenumber(100.0) + b = mhz_to_wavenumber(200.0) + assert a is not None and b is not None + assert pytest.approx(b, rel=1e-9) == 2 * a + + +# ============================================================ +# PrincipalMoments +# ============================================================ + + +class TestPrincipalMoments: + def test_rotor_type_monoatomic(self, single_atom): + pm = single_atom.calc_moments_of_inertia() + assert pm is not None + assert pm.rotor_type() == RotorType.MONOATOMIC + + def test_rotor_type_linear(self, hcl): + pm = hcl.calc_moments_of_inertia() + assert pm is not None + assert pm.rotor_type() == RotorType.LINEAR + + def test_rotor_type_spherical_top(self, methane): + pm = methane.calc_moments_of_inertia() + assert pm is not None + assert pm.rotor_type() == RotorType.SPHERICAL_TOP + + def test_rotor_type_oblate_top(self, benzene): + pm = benzene.calc_moments_of_inertia() + assert pm is not None + assert pm.rotor_type() == RotorType.OBLATE_TOP + + def test_rotor_type_oblate_top_ammonia(self, ammonia): + pm = ammonia.calc_moments_of_inertia() + assert pm is not None + assert pm.rotor_type() == RotorType.OBLATE_TOP + + def test_rotor_type_prolate_top(self, ch3cl): + pm = ch3cl.calc_moments_of_inertia() + assert pm is not None + assert pm.rotor_type() == RotorType.PROLATE_TOP + + def test_rotor_type_asymmetric_top(self, water): + pm = water.calc_moments_of_inertia() + assert pm is not None + assert pm.rotor_type() == RotorType.ASYMMETRIC_TOP + + def test_moments_sorted_ascending(self, water): + pm = water.calc_moments_of_inertia() + assert pm is not None + assert pm.Ia <= pm.Ib <= pm.Ic + + def test_axes_shape(self, water): + pm = water.calc_moments_of_inertia() + assert pm is not None + assert pm.axes.shape == (3, 3) + + def test_axes_orthonormal(self, water): + pm = water.calc_moments_of_inertia() + assert pm is not None + product = pm.axes.T @ pm.axes + np.testing.assert_allclose(product, np.eye(3), atol=1e-10) + + def test_str_output(self, water): + pm = water.calc_moments_of_inertia() + assert pm is not None + s = str(pm) + assert "Moments of inertia (amu·Å²):" in s + assert "Ia" in s + assert "Ib" in s + assert "Ic" in s + + +# ============================================================ +# calc_moments_of_inertia +# ============================================================ + + +class TestCalcMomentsOfInertia: + def test_returns_none_for_empty_structure(self): + assert make_no_real_atoms_structure().calc_moments_of_inertia() is None + + def test_returns_none_for_all_zero_masses(self, water): + masses = np.zeros(3) + assert water.calc_moments_of_inertia(masses=masses) is None + + def test_raises_on_wrong_masses_length(self, water): + with pytest.raises(ValueError, match="masses length"): + water.calc_moments_of_inertia(masses=np.array([1.0, 2.0])) + + def test_ghost_atoms_ignored(self): + """PointCharge should not contribute to moments.""" + real = make_structure( + make_atom("O", 0.0, 0.0, 0.119748), + make_atom("H", 0.0, 0.756950, -0.478993), + make_atom("H", 0.0, -0.756950, -0.478993), + ) + with_pc = Structure( + atoms=[ + Atom( + element=Element("O"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.119748)) + ), + Atom( + element=Element("H"), + coordinates=Coordinates(coordinates=(0.0, 0.756950, -0.478993)), + ), + Atom( + element=Element("H"), + coordinates=Coordinates(coordinates=(0.0, -0.756950, -0.478993)), + ), + PointCharge(coordinates=Coordinates(coordinates=(5.0, 5.0, 5.0)), charge=1.0), + ] + ) + pm_real = real.calc_moments_of_inertia() + pm_pc = with_pc.calc_moments_of_inertia() + assert pm_real is not None and pm_pc is not None + assert pytest.approx(pm_real.Ia, rel=1e-6) == pm_pc.Ia + assert pytest.approx(pm_real.Ib, rel=1e-6) == pm_pc.Ib + assert pytest.approx(pm_real.Ic, rel=1e-6) == pm_pc.Ic + + def test_custom_masses_override(self, water): + """Passing explicit masses should change the moments.""" + pm_default = water.calc_moments_of_inertia() + pm_custom = water.calc_moments_of_inertia(masses=np.array([18.0, 2.0, 2.0])) + assert pm_default is not None and pm_custom is not None + assert pm_default.Ia != pytest.approx(pm_custom.Ia, rel=1e-3) + + def test_weights_per_element(self, water): + """Per-element weight override should change the moments.""" + pm_default = water.calc_moments_of_inertia() + pm_weights = water.calc_moments_of_inertia(weights={"H": 2.014}) + assert pm_default is not None and pm_weights is not None + assert pm_default.Ic != pytest.approx(pm_weights.Ic, rel=1e-3) + + def test_atom_weights_override(self, water): + """Per-atom weight override (index 0) should change the moments.""" + pm_default = water.calc_moments_of_inertia() + pm_aw = water.calc_moments_of_inertia(atom_weights={0: 17.999}) + assert pm_default is not None and pm_aw is not None + assert pm_default.Ic != pytest.approx(pm_aw.Ic, rel=1e-3) + + def test_atom_weights_priority_over_weights(self, water): + """atom_weights should take priority over weights for the same atom.""" + pm_aw = water.calc_moments_of_inertia(atom_weights={0: 17.999}) + pm_both = water.calc_moments_of_inertia(weights={"O": 16.0}, atom_weights={0: 17.999}) + assert pm_aw is not None and pm_both is not None + assert pytest.approx(pm_aw.Ic, rel=1e-9) == pm_both.Ic + + def test_unknown_element_warns_and_excludes(self): + """Atoms with unknown elements should warn and be assigned mass 0.""" + structure = Structure( + atoms=[ + Atom(element=Element("C"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), + Atom(element=Element("H"), coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0))), + ] + ) + # Patch via weights to simulate unknown: pass mass 0 explicitly for C + pm = structure.calc_moments_of_inertia(atom_weights={0: 0.0}) + # Only H contributes → single non-zero mass → linear or monoatomic + assert pm is not None + + def test_zero_mass_atoms_filtered(self, water): + """Atoms with zero mass are excluded; remaining atoms determine the rotor type.""" + # Zero out the two H atoms → only O remains → monoatomic + masses = np.array([15.999, 0.0, 0.0]) + pm = water.calc_moments_of_inertia(masses=masses) + assert pm is not None + assert pm.rotor_type() == RotorType.MONOATOMIC + + def test_water_moments_known_values(self, water): + """ + Check moments are self-consistent: Ia <= Ib <= Ic and all positive. + The geometry used here gives Ia ≈ 0.6418 amu·Å² (not the experimental + equilibrium value of 0.5791) because the O-H distance and angle differ + slightly. We verify internal consistency rather than an absolute value. + """ + pm = water.calc_moments_of_inertia() + assert pm is not None + assert pm.Ia > 0.0 + assert pm.Ia <= pm.Ib <= pm.Ic + # Rough sanity check: moments should be in a physically reasonable range + assert 0.1 < pm.Ia < 5.0 + assert 0.1 < pm.Ib < 5.0 + assert 0.1 < pm.Ic < 5.0 + + +# ============================================================ +# calc_rotor_type +# ============================================================ + + +class TestCalcRotorType: + def test_with_precomputed_moments(self, water): + pm = water.calc_moments_of_inertia() + assert water.calc_rotor_type(moments=pm) == RotorType.ASYMMETRIC_TOP + + def test_without_precomputed_moments(self, water): + assert water.calc_rotor_type() == RotorType.ASYMMETRIC_TOP + + def test_returns_none_for_empty_structure(self): + assert make_no_real_atoms_structure().calc_rotor_type() is None + + def test_all_rotor_types(self, single_atom, hcl, methane, benzene, ammonia, ch3cl, water): + assert single_atom.calc_rotor_type() == RotorType.MONOATOMIC + assert hcl.calc_rotor_type() == RotorType.LINEAR + assert methane.calc_rotor_type() == RotorType.SPHERICAL_TOP + assert benzene.calc_rotor_type() == RotorType.OBLATE_TOP + assert ammonia.calc_rotor_type() == RotorType.OBLATE_TOP + assert ch3cl.calc_rotor_type() == RotorType.PROLATE_TOP + assert water.calc_rotor_type() == RotorType.ASYMMETRIC_TOP + + +# ============================================================ +# calc_rotational_constants +# ============================================================ + + +class TestCalcRotationalConstants: + def test_returns_none_for_empty_structure(self): + assert make_no_real_atoms_structure().calc_rotational_constants() is None + + def test_returns_rotational_constants(self, water): + rc = water.calc_rotational_constants() + assert rc is not None + assert isinstance(rc, RotationalConstants) + + def test_abc_ordering(self, water): + """A ≥ B ≥ C by convention (ascending moments → descending constants).""" + rc = water.calc_rotational_constants() + assert rc is not None + assert rc.A is not None and rc.B is not None and rc.C is not None + assert rc.A >= rc.B >= rc.C + + def test_linear_molecule_has_none_constants(self, hcl): + """HCl is linear: Ia = 0, so A should be None.""" + rc = hcl.calc_rotational_constants() + assert rc is not None + assert rc.A is None + + def test_monoatomic_all_none(self, single_atom): + rc = single_atom.calc_rotational_constants() + assert rc is not None + assert rc.A is None + assert rc.B is None + assert rc.C is None + + def test_get_in_wavenumbers_returns_tuple(self, water): + rc = water.calc_rotational_constants() + assert rc is not None + wn = rc.get_in_wavenumbers() + assert isinstance(wn, tuple) + assert len(wn) == 3 + + def test_get_in_wavenumbers_none_propagates(self, hcl): + rc = hcl.calc_rotational_constants() + assert rc is not None + wn = rc.get_in_wavenumbers() + assert wn[0] is None # A is None for linear + + def test_get_in_wavenumbers_values_positive(self, water): + rc = water.calc_rotational_constants() + assert rc is not None + wn = rc.get_in_wavenumbers() + assert all(v is not None and v > 0 for v in wn) + + def test_get_in_wavenumbers_consistent_with_mhz(self, water): + """Manual conversion should match get_in_wavenumbers.""" + rc = water.calc_rotational_constants() + assert rc is not None + wn = rc.get_in_wavenumbers() + assert wn[0] is not None and rc.A is not None + assert pytest.approx(wn[0], rel=1e-9) == mhz_to_wavenumber(rc.A) + + def test_str_output(self, water): + rc = water.calc_rotational_constants() + assert rc is not None + s = str(rc) + assert "Rotational constants" in s + assert "MHz" in s + assert "cm⁻¹" in s + + def test_a_not_in_sync_problem(self, water): + """ + Mutating rc.A should NOT affect get_in_wavenumbers, + since wavenumbers are computed on the fly. + """ + rc = water.calc_rotational_constants() + assert rc is not None + wn_before = rc.get_in_wavenumbers()[0] + rc.A = rc.A * 5 if rc.A is not None else None # type: ignore[operator] + wn_after = rc.get_in_wavenumbers()[0] + assert wn_before != wn_after # confirms on-the-fly computation From bd8da5793e28d8b8a415f2d4914d31a2e0d2eb0c Mon Sep 17 00:00:00 2001 From: crjacinto Date: Tue, 19 May 2026 19:47:32 +0200 Subject: [PATCH 10/27] Corrected double ticks --- src/opi/input/structures/structure.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 4a7af74e..4706bee1 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1433,16 +1433,16 @@ def calc_rotor_type( ---------- moments : ndarray shape (3,) | None, default None Pre-computed principal moments (amu·Å², ascending). When - ``None`` the moments are computed via + `None` the moments are computed via :meth:`calc_moment_of_inertia` using *mass_kwargs*. **mass_kwargs Forwarded to :meth:`calc_moment_of_inertia` when *moments* is - ``None`` (i.e. ``masses``, ``weights``, ``atom_weights``). + `None` (i.e. `masses`, `weights`, `atom_weights`). Returns ------- RotorType | None - ``None`` if the structure has no real atoms or all masses vanish. + `None` if the structure has no real atoms or all masses vanish. """ if moments is None: result = self.calc_moments_of_inertia(**mass_kwargs) @@ -1464,7 +1464,7 @@ def calc_rotational_constants( """ Compute rotational constants for this structure. - Only `Atom`` instances contribute; `GhostAtom`, `PointCharge`, + Only `Atom` instances contribute; `GhostAtom`, `PointCharge`, and `EmbeddingPotential` atoms are silently ignored. Mass priority From 03ba66df62181154ba4a8195239628a901521cb1 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Fri, 22 May 2026 15:43:17 +0200 Subject: [PATCH 11/27] Changed Atom lists to Indeces list. Added property for real_atoms --- src/opi/input/structures/structure.py | 116 +++++++++++++++----------- 1 file changed, 65 insertions(+), 51 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 4706bee1..40985793 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1074,26 +1074,35 @@ def _iter_xyz_structures( if n_struc_limit and n_struc >= n_struc_limit: break + @property + def real_atoms(self) -> list[Atom]: + """ + Return only real `Atom` instances, excluding `GhostAtom`, `PointCharge`, + and `EmbeddingPotential`. Note that these are subclasses of `Atom`, + so `type(a) is Atom` is used rather than `isinstance`. + """ + return [a for a in self.atoms if type(a) is Atom] + def get_coordinates( self, - atoms: Sequence[Atom | EmbeddingPotential | GhostAtom | PointCharge] | None = None, + only_atoms: Sequence[int] = (), ) -> npt.NDArray[np.float64]: """ - Return the coordinates of the given atoms as a numpy array. + Return the coordinates of atoms as a numpy array. Parameters ---------- - atoms : sequence of atom-like objects | None, default None - The atoms to extract coordinates from. If `None`, all entries - in `self.atoms` are used. + only_atoms : Sequence[int], default: () + Indices of atoms in `self.atoms` to include. If empty, all + atoms are included. Returns ------- npt.NDArray[np.float64], shape (N, 3) - Coordinates in the same order as the input list. + Coordinates in the same order as the selected atoms. """ - atom_list = atoms if atoms is not None else self.atoms - return np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) + atom_list = [self.atoms[i] for i in only_atoms] if only_atoms else self.atoms + return np.array([a.coordinates.coordinates for a in atom_list]) def set_coordinates(self, coords: npt.NDArray[np.float64]) -> "Structure": """ @@ -1137,31 +1146,33 @@ def centered_structure(self) -> "Structure": Structure New `Structure` with centered coordinates. """ - real_atoms = [a for a in self.atoms if type(a) is Atom] - centroid = self.get_coordinates(atoms=real_atoms).mean(axis=0) + real_indices = [i for i, a in enumerate(self.atoms) if type(a) is Atom] + centroid = self.get_coordinates(only_atoms=real_indices).mean(axis=0) return self.set_coordinates(self.get_coordinates() - centroid) def _filtered_atoms( self, only_atoms: Sequence[int], ignore_hs: bool, - ) -> list[Atom]: + ) -> list[int]: """ - Return real `Atom` instances after applying `only_atoms` and `ignore_hs` filters. + Return indices of real `Atom` instances after applying `only_atoms` + and `ignore_hs` filters. The final `type(a) is Atom` check ensures that explicitly indexed atoms (via `only_atoms`) cannot introduce `GhostAtom`, `PointCharge`, or `EmbeddingPotential` instances into the returned list. Note that `GhostAtom` is a subclass of `Atom`, so `isinstance` is intentionally avoided here. """ - atom_list = self.atoms if only_atoms: - candidates = [atom_list[i] for i in only_atoms] + candidates = list(only_atoms) elif ignore_hs: - candidates = [a for a in atom_list if type(a) is Atom and a.element != Element.H] + candidates = [ + i for i, a in enumerate(self.atoms) if type(a) is Atom and a.element != Element.H + ] else: - candidates = [a for a in atom_list if type(a) is Atom] - return [a for a in candidates if type(a) is Atom] + candidates = [i for i, a in enumerate(self.atoms) if type(a) is Atom] + return [i for i in candidates if type(self.atoms[i]) is Atom] @staticmethod def _validate_rmsd_compatibility(atoms1: list[Atom], atoms2: list[Atom]) -> None: @@ -1189,23 +1200,23 @@ def _rmsd_coords( def get_coordinates_at_centroid( self, - atoms: Sequence[Atom | EmbeddingPotential | GhostAtom | PointCharge] | None = None, + only_atoms: Sequence[int] = (), ) -> npt.NDArray[np.float64]: """ Return coordinates translated to the centroid. Parameters ---------- - atoms : sequence of atom-like objects | None, default None - The atoms to extract coordinates from. If `None`, all entries - in `self.atoms` are used. + only_atoms : Sequence[int], default: () + Indices of atoms in `self.atoms` to include. If empty, all + atoms are included. Returns ------- npt.NDArray[np.float64], shape (N, 3) - Centroid-translated coordinates in the same order as the input list. + Centroid-translated coordinates in the same order as the selected atoms. """ - coords = self.get_coordinates(atoms=atoms) + coords = self.get_coordinates(only_atoms=only_atoms) return coords - np.mean(coords, axis=0, dtype=np.float64) def rmsd( @@ -1240,12 +1251,14 @@ def rmsd( ValueError If the filtered atom sets differ in size or element order. """ - atoms1 = self._filtered_atoms(only_atoms, ignore_hs) - atoms2 = other._filtered_atoms(only_atoms, ignore_hs) - self._validate_rmsd_compatibility(atoms1, atoms2) - - coords1 = self.get_coordinates_at_centroid(atoms=atoms1) - coords2 = other.get_coordinates_at_centroid(atoms=atoms2) + indices1 = self._filtered_atoms(only_atoms, ignore_hs) + indices2 = other._filtered_atoms(only_atoms, ignore_hs) + self._validate_rmsd_compatibility( + cast(list[Atom], [self.atoms[i] for i in indices1]), + cast(list[Atom], [other.atoms[i] for i in indices2]), + ) + coords1 = self.get_coordinates_at_centroid(only_atoms=indices1) + coords2 = other.get_coordinates_at_centroid(only_atoms=indices2) return self._rmsd_coords(coords1, coords2) @@ -1267,9 +1280,9 @@ def rmsd_kabsch( ---------- other : Structure Structure to compare against. - only_atoms : Sequence[int], default () + only_atoms : Sequence[int], default: () Atom indices to include. If non-empty, `ignore_hs` is ignored. - ignore_hs : bool, default False + ignore_hs : bool, default: False Exclude hydrogen atoms from the RMSD computation. Returns @@ -1282,13 +1295,16 @@ def rmsd_kabsch( ValueError If the filtered atom sets differ in size or element order. """ - atoms1 = self._filtered_atoms(only_atoms, ignore_hs) - atoms2 = other._filtered_atoms(only_atoms, ignore_hs) - self._validate_rmsd_compatibility(atoms1, atoms2) + indices1 = self._filtered_atoms(only_atoms, ignore_hs) + indices2 = other._filtered_atoms(only_atoms, ignore_hs) + self._validate_rmsd_compatibility( + cast(list[Atom], [self.atoms[i] for i in indices1]), + cast(list[Atom], [other.atoms[i] for i in indices2]), + ) # Getting centered coordinates - coords1 = self.get_coordinates_at_centroid(atoms=atoms1) - coords2 = other.get_coordinates_at_centroid(atoms=atoms2) + coords1 = self.get_coordinates_at_centroid(only_atoms=indices1) + coords2 = other.get_coordinates_at_centroid(only_atoms=indices2) # Kabsch algorithm: find optimal rotation matrix via SVD # @@ -1331,24 +1347,24 @@ def calc_moments_of_inertia( """ Compute the principal axes and moments of inertia for this structure. - Mass priority - ------------- - masses > atom_weights > weights > default (ATOMIC_MASSES_FROM_ELEMENT) - Parameters ---------- - masses : npt.NDArray[np.float64] | None, default None + masses : npt.NDArray[np.float64] | None, default: None Per-atom masses (amu) for all Atom instances, overriding every other mass source. - atom_weights : dict[int, float] | None, default None + atom_weights : dict[int, float] | None, default: None Per-atom mass overrides keyed by index within the Atom list. - weights : dict[str, float] | None, default None + weights : dict[str, float] | None, default: None Per-element mass overrides keyed by element symbol string (e.g. `{"C": 13.003}`). Only `Atom` instances contribute; `GhostAtom`, `PointCharge`, and `EmbeddingPotential` atoms are silently ignored. + Mass priority + ------------- + masses > atom_weights > weights > default: (ATOMIC_MASSES_FROM_ELEMENT) + Returns ------- PrincipalMoments | None @@ -1361,12 +1377,10 @@ def calc_moments_of_inertia( If the length of *masses* does not match the number of `Atom` instances in the structure. """ - atom_list = [a for a in self.atoms if type(a) is Atom] - if not atom_list: + atom_indices = [i for i, a in enumerate(self.atoms) if type(a) is Atom] + if not atom_indices: return None - - # Extract coordinates from real atoms only - coords = self.get_coordinates(atoms=atom_list) + coords = self.get_coordinates(only_atoms=atom_indices) # --- Mass assignment priority: masses > atom_weights > weights > default --- # Normalise the optional dicts to avoid None checks below @@ -1375,14 +1389,14 @@ def calc_moments_of_inertia( if masses is not None: masses = np.asarray(masses, dtype=np.float64) - if len(masses) != len(atom_list): + if len(masses) != len(self.real_atoms): raise ValueError( f"masses length ({len(masses)}) does not match " - f"number of atoms ({len(atom_list)})" + f"number of atoms ({len(self.real_atoms)})" ) else: masses_list: list[float] = [] - for i, atom in enumerate(atom_list): + for i, atom in enumerate(self.real_atoms): if i in atom_weights: m = atom_weights[i] elif atom.element.value in weights: From 999b4af1a46315cb388a2ef9d5dc16efd8b2bd23 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Fri, 22 May 2026 15:47:40 +0200 Subject: [PATCH 12/27] Updated test_rmsd to use new functions in structure.py --- tests/unit/test_rmsd.py | 176 ++++++++++++++++++++++------------------ 1 file changed, 99 insertions(+), 77 deletions(-) diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_rmsd.py index 7eb1093b..442c12b9 100644 --- a/tests/unit/test_rmsd.py +++ b/tests/unit/test_rmsd.py @@ -2,6 +2,7 @@ Unit tests for RMSD and coordinate utilities. Covers: +- real_atoms property - get_coordinates() - get_coordinates_at_centroid() - set_coordinates() @@ -18,7 +19,7 @@ import numpy as np import pytest -from opi.input.structures.atom import Atom, PointCharge +from opi.input.structures.atom import Atom, GhostAtom, PointCharge from opi.input.structures.coordinates import Coordinates from opi.input.structures.structure import Structure from opi.utils.element import Element @@ -27,7 +28,6 @@ # Helpers # ============================================================ - def make_atom(element: str, x: float, y: float, z: float) -> Atom: return Atom( element=Element(element), @@ -43,14 +43,13 @@ def make_structure(*atoms: Atom) -> Structure: # Fixtures # ============================================================ - @pytest.fixture() def water() -> Structure: """H2O at a known geometry.""" return make_structure( - make_atom("O", 0.000000, 0.000000, 0.119748), - make_atom("H", 0.000000, 0.756950, -0.478993), - make_atom("H", 0.000000, -0.756950, -0.478993), + make_atom("O", 0.000000, 0.000000, 0.119748), + make_atom("H", 0.000000, 0.756950, -0.478993), + make_atom("H", 0.000000, -0.756950, -0.478993), ) @@ -58,9 +57,9 @@ def water() -> Structure: def water_translated() -> Structure: """H2O shifted by (1, 2, 3) — RMSD vs water should be 0 after centring.""" return make_structure( - make_atom("O", 1.000000, 2.000000, 3.119748), - make_atom("H", 1.000000, 2.756950, 2.521007), - make_atom("H", 1.000000, 1.243050, 2.521007), + make_atom("O", 1.000000, 2.000000, 3.119748), + make_atom("H", 1.000000, 2.756950, 2.521007), + make_atom("H", 1.000000, 1.243050, 2.521007), ) @@ -73,13 +72,11 @@ def water_rotated(water) -> Structure: centered = water.centered_structure() coords = centered.get_coordinates() theta = np.pi / 2 - R = np.array( - [ - [np.cos(theta), -np.sin(theta), 0], - [np.sin(theta), np.cos(theta), 0], - [0, 0, 1], - ] - ) + R = np.array([ + [np.cos(theta), -np.sin(theta), 0], + [np.sin(theta), np.cos(theta), 0], + [0, 0, 1], + ]) return centered.set_coordinates(coords @ R.T) @@ -90,34 +87,65 @@ def ethanol() -> Structure: Coordinates are approximate — correctness of geometry is not critical here. """ return make_structure( - make_atom("C", 0.000, 0.000, 0.000), - make_atom("C", 1.540, 0.000, 0.000), - make_atom("O", 2.060, 1.190, 0.000), - make_atom("H", -0.390, 1.020, 0.000), - make_atom("H", -0.390, -0.510, 0.890), + make_atom("C", 0.000, 0.000, 0.000), + make_atom("C", 1.540, 0.000, 0.000), + make_atom("O", 2.060, 1.190, 0.000), + make_atom("H", -0.390, 1.020, 0.000), + make_atom("H", -0.390, -0.510, 0.890), make_atom("H", -0.390, -0.510, -0.890), - make_atom("H", 1.930, -0.510, 0.890), - make_atom("H", 1.930, -0.510, -0.890), - make_atom("H", 2.980, 1.190, 0.000), + make_atom("H", 1.930, -0.510, 0.890), + make_atom("H", 1.930, -0.510, -0.890), + make_atom("H", 2.980, 1.190, 0.000), ) @pytest.fixture() def mixed_structure() -> Structure: """Structure with a real Atom and a PointCharge (not an Atom subclass).""" - return Structure( - atoms=[ + return Structure(atoms=[ + Atom(element=Element("O"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), + PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), + ]) + + +# ============================================================ +# real_atoms +# ============================================================ + +class TestRealAtoms: + def test_returns_only_real_atoms(self, mixed_structure): + """PointCharge should be excluded.""" + assert len(mixed_structure.real_atoms) == 1 + assert all(type(a) is Atom for a in mixed_structure.real_atoms) + + def test_ghost_atom_excluded(self): + """GhostAtom is a subclass of Atom but should be excluded.""" + structure = Structure(atoms=[ Atom(element=Element("O"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), - PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), - ] - ) + GhostAtom(element=Element("H"), coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0))), + ]) + assert len(structure.real_atoms) == 1 + assert type(structure.real_atoms[0]) is Atom + + def test_all_real_atoms(self, water): + """All atoms in water are real.""" + assert len(water.real_atoms) == 3 + + def test_return_type(self, water): + assert isinstance(water.real_atoms, list) + assert all(type(a) is Atom for a in water.real_atoms) + + def test_empty_structure(self): + structure = Structure(atoms=[ + PointCharge(coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0)), charge=1.0) + ]) + assert water.real_atoms == [] if False else structure.real_atoms == [] # ============================================================ # get_coordinates # ============================================================ - class TestGetCoordinates: def test_shape(self, water): coords = water.get_coordinates() @@ -131,10 +159,9 @@ def test_values(self, water): coords = water.get_coordinates() np.testing.assert_allclose(coords[0], [0.0, 0.0, 0.119748]) - def test_atoms_kwarg_filters(self, mixed_structure): - """Passing only real Atom instances should exclude PointCharge.""" - real = [a for a in mixed_structure.atoms if type(a) is Atom] - coords = mixed_structure.get_coordinates(atoms=real) + def test_only_atoms_filters(self, mixed_structure): + """Passing index 0 only should return one row.""" + coords = mixed_structure.get_coordinates(only_atoms=[0]) assert coords.shape == (1, 3) def test_default_includes_all(self, mixed_structure): @@ -142,6 +169,11 @@ def test_default_includes_all(self, mixed_structure): coords = mixed_structure.get_coordinates() assert coords.shape == (2, 3) + def test_only_atoms_subset(self, water): + """Passing indices [0, 1] should return only those two atoms.""" + coords = water.get_coordinates(only_atoms=[0, 1]) + assert coords.shape == (2, 3) + def test_single_atom(self): s = make_structure(make_atom("C", 1.0, 2.0, 3.0)) coords = s.get_coordinates() @@ -152,7 +184,6 @@ def test_single_atom(self): # get_coordinates_at_centroid # ============================================================ - class TestGetCoordinatesAtCentroid: def test_centroid_is_zero(self, water): coords = water.get_coordinates_at_centroid() @@ -162,9 +193,8 @@ def test_shape_preserved(self, water): coords = water.get_coordinates_at_centroid() assert coords.shape == (3, 3) - def test_with_atoms_kwarg(self, water): - atoms = water.atoms[:2] - coords = water.get_coordinates_at_centroid(atoms=atoms) + def test_with_only_atoms(self, water): + coords = water.get_coordinates_at_centroid(only_atoms=[0, 1]) assert coords.shape == (2, 3) np.testing.assert_allclose(coords.mean(axis=0), [0.0, 0.0, 0.0], atol=1e-12) @@ -178,7 +208,6 @@ def test_translation_invariant(self, water, water_translated): # set_coordinates # ============================================================ - class TestSetCoordinates: def test_returns_new_structure(self, water): new_coords = water.get_coordinates() + 1.0 @@ -214,12 +243,11 @@ def test_multiplicity_preserved(self, water): # centered_structure # ============================================================ - class TestCenteredStructure: def test_centroid_at_origin(self, water): centered = water.centered_structure() - real_atoms = [a for a in centered.atoms if isinstance(a, Atom)] - coords = centered.get_coordinates(atoms=real_atoms) + real_indices = [i for i, a in enumerate(centered.atoms) if type(a) is Atom] + coords = centered.get_coordinates(only_atoms=real_indices) np.testing.assert_allclose(coords.mean(axis=0), [0.0, 0.0, 0.0], atol=1e-12) def test_returns_new_structure(self, water): @@ -236,9 +264,8 @@ def test_pointcharge_excluded_from_centroid(self, mixed_structure): from the real Atom at (0, 0, 0). """ centered = mixed_structure.centered_structure() - real_atoms = [a for a in centered.atoms if type(a) is Atom] - coords = centered.get_coordinates(atoms=real_atoms) - # Real atom was already at origin, so it should stay at origin + real_indices = [i for i, a in enumerate(centered.atoms) if type(a) is Atom] + coords = centered.get_coordinates(only_atoms=real_indices) np.testing.assert_allclose(coords[0], [0.0, 0.0, 0.0], atol=1e-12) def test_translated_centered_equals_original_centered(self, water, water_translated): @@ -251,57 +278,56 @@ def test_translated_centered_equals_original_centered(self, water, water_transla # _filtered_atoms # ============================================================ - class TestFilteredAtoms: - def test_default_returns_only_atoms(self, mixed_structure): - """PointCharge should be excluded since it is not an Atom instance.""" - filtered = mixed_structure._filtered_atoms((), False) - assert all(isinstance(a, Atom) for a in filtered) - assert len(filtered) == 1 + def test_default_returns_only_atom_indices(self, mixed_structure): + """PointCharge should be excluded — only index 0 (the Atom) returned.""" + indices = mixed_structure._filtered_atoms((), False) + assert indices == [0] + assert all(type(mixed_structure.atoms[i]) is Atom for i in indices) + + def test_returns_list_of_int(self, water): + indices = water._filtered_atoms((), False) + assert isinstance(indices, list) + assert all(isinstance(i, int) for i in indices) def test_ignore_hs(self, ethanol): - filtered = ethanol._filtered_atoms((), True) - elements = [a.element for a in filtered] + indices = ethanol._filtered_atoms((), True) + elements = [ethanol.atoms[i].element for i in indices] assert Element("H") not in elements def test_only_atoms_indices(self, ethanol): - filtered = ethanol._filtered_atoms([0, 1, 2], False) - assert len(filtered) == 3 + indices = ethanol._filtered_atoms([0, 1, 2], False) + assert indices == [0, 1, 2] def test_only_atoms_excludes_pointcharge(self): - """Explicitly indexed PointCharge must be excluded by the final isinstance check.""" - structure = Structure( - atoms=[ - Atom(element=Element("C"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), - PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), - ] - ) - filtered = structure._filtered_atoms([0, 1], False) - assert all(isinstance(a, Atom) for a in filtered) - assert len(filtered) == 1 + """Explicitly indexed PointCharge must be excluded by the final type check.""" + structure = Structure(atoms=[ + Atom(element=Element("C"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), + PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), + ]) + indices = structure._filtered_atoms([0, 1], False) + assert indices == [0] + assert all(type(structure.atoms[i]) is Atom for i in indices) def test_only_atoms_takes_priority_over_ignore_hs(self, ethanol): """When only_atoms is set, ignore_hs should be ignored.""" - filtered_with = ethanol._filtered_atoms([3, 4], True) - filtered_without = ethanol._filtered_atoms([3, 4], False) - assert len(filtered_with) == len(filtered_without) + indices_with = ethanol._filtered_atoms([3, 4], True) + indices_without = ethanol._filtered_atoms([3, 4], False) + assert indices_with == indices_without # ============================================================ # _validate_rmsd_compatibility # ============================================================ - class TestValidateRmsdCompatibility: def test_compatible_structures(self, water): - atoms = [a for a in water.atoms if isinstance(a, Atom)] + atoms = water.real_atoms Structure._validate_rmsd_compatibility(atoms, atoms) # should not raise def test_raises_on_different_count(self, water, ethanol): - a1 = [a for a in water.atoms if isinstance(a, Atom)] - a2 = [a for a in ethanol.atoms if isinstance(a, Atom)] with pytest.raises(ValueError, match="different number of atoms"): - Structure._validate_rmsd_compatibility(a1, a2) + Structure._validate_rmsd_compatibility(water.real_atoms, ethanol.real_atoms) def test_raises_on_element_mismatch(self): a1 = [make_atom("C", 0, 0, 0), make_atom("H", 1, 0, 0)] @@ -321,7 +347,6 @@ def test_error_message_uses_natural_counting(self): # rmsd # ============================================================ - class TestRmsd: def test_identical_structures_zero_rmsd(self, water): assert pytest.approx(water.rmsd(water), abs=1e-10) == 0.0 @@ -362,8 +387,6 @@ def test_only_atoms_subset(self, ethanol): def test_nonzero_rmsd_for_different_structures(self, water): other = water.set_coordinates(water.get_coordinates() + np.array([0, 0, 1.0])) - # After centring the z-shift cancels, but individual atom positions differ - # We just confirm it's finite and non-negative result = water.rmsd(other) assert result >= 0.0 @@ -375,7 +398,6 @@ def test_result_is_float(self, water): # rmsd_kabsch # ============================================================ - class TestRmsdKabsch: def test_identical_structures_zero_rmsd(self, water): assert pytest.approx(water.rmsd_kabsch(water), abs=1e-10) == 0.0 From 93d55d00539374d95378bb12e435c9e3112b56df Mon Sep 17 00:00:00 2001 From: crjacinto Date: Fri, 22 May 2026 16:10:28 +0200 Subject: [PATCH 13/27] Fixed a typo --- src/opi/utils/rotconst.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index 05575330..62839fea 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -125,7 +125,7 @@ def fmt(x: float | None, unit: str = "") -> str: A_cm, B_cm, C_cm = self.get_in_wavenumbers() return ( - "Rotational constants (yeih):\n" + "Rotational constants:\n" f" A = {fmt(self.A, 'MHz')} ({fmt(A_cm, 'cm⁻¹')})\n" f" B = {fmt(self.B, 'MHz')} ({fmt(B_cm, 'cm⁻¹')})\n" f" C = {fmt(self.C, 'MHz')} ({fmt(C_cm, 'cm⁻¹')})" From 45bdfd0e5cad4f7760190bbce87b31fe07e77528 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Fri, 22 May 2026 16:16:16 +0200 Subject: [PATCH 14/27] Ran nox session --- tests/unit/test_rmsd.py | 89 +++++++++++++++++++++++++---------------- 1 file changed, 55 insertions(+), 34 deletions(-) diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_rmsd.py index 442c12b9..338429bd 100644 --- a/tests/unit/test_rmsd.py +++ b/tests/unit/test_rmsd.py @@ -28,6 +28,7 @@ # Helpers # ============================================================ + def make_atom(element: str, x: float, y: float, z: float) -> Atom: return Atom( element=Element(element), @@ -43,13 +44,14 @@ def make_structure(*atoms: Atom) -> Structure: # Fixtures # ============================================================ + @pytest.fixture() def water() -> Structure: """H2O at a known geometry.""" return make_structure( - make_atom("O", 0.000000, 0.000000, 0.119748), - make_atom("H", 0.000000, 0.756950, -0.478993), - make_atom("H", 0.000000, -0.756950, -0.478993), + make_atom("O", 0.000000, 0.000000, 0.119748), + make_atom("H", 0.000000, 0.756950, -0.478993), + make_atom("H", 0.000000, -0.756950, -0.478993), ) @@ -57,9 +59,9 @@ def water() -> Structure: def water_translated() -> Structure: """H2O shifted by (1, 2, 3) — RMSD vs water should be 0 after centring.""" return make_structure( - make_atom("O", 1.000000, 2.000000, 3.119748), - make_atom("H", 1.000000, 2.756950, 2.521007), - make_atom("H", 1.000000, 1.243050, 2.521007), + make_atom("O", 1.000000, 2.000000, 3.119748), + make_atom("H", 1.000000, 2.756950, 2.521007), + make_atom("H", 1.000000, 1.243050, 2.521007), ) @@ -72,11 +74,13 @@ def water_rotated(water) -> Structure: centered = water.centered_structure() coords = centered.get_coordinates() theta = np.pi / 2 - R = np.array([ - [np.cos(theta), -np.sin(theta), 0], - [np.sin(theta), np.cos(theta), 0], - [0, 0, 1], - ]) + R = np.array( + [ + [np.cos(theta), -np.sin(theta), 0], + [np.sin(theta), np.cos(theta), 0], + [0, 0, 1], + ] + ) return centered.set_coordinates(coords @ R.T) @@ -87,31 +91,34 @@ def ethanol() -> Structure: Coordinates are approximate — correctness of geometry is not critical here. """ return make_structure( - make_atom("C", 0.000, 0.000, 0.000), - make_atom("C", 1.540, 0.000, 0.000), - make_atom("O", 2.060, 1.190, 0.000), - make_atom("H", -0.390, 1.020, 0.000), - make_atom("H", -0.390, -0.510, 0.890), + make_atom("C", 0.000, 0.000, 0.000), + make_atom("C", 1.540, 0.000, 0.000), + make_atom("O", 2.060, 1.190, 0.000), + make_atom("H", -0.390, 1.020, 0.000), + make_atom("H", -0.390, -0.510, 0.890), make_atom("H", -0.390, -0.510, -0.890), - make_atom("H", 1.930, -0.510, 0.890), - make_atom("H", 1.930, -0.510, -0.890), - make_atom("H", 2.980, 1.190, 0.000), + make_atom("H", 1.930, -0.510, 0.890), + make_atom("H", 1.930, -0.510, -0.890), + make_atom("H", 2.980, 1.190, 0.000), ) @pytest.fixture() def mixed_structure() -> Structure: """Structure with a real Atom and a PointCharge (not an Atom subclass).""" - return Structure(atoms=[ - Atom(element=Element("O"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), - PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), - ]) + return Structure( + atoms=[ + Atom(element=Element("O"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), + PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), + ] + ) # ============================================================ # real_atoms # ============================================================ + class TestRealAtoms: def test_returns_only_real_atoms(self, mixed_structure): """PointCharge should be excluded.""" @@ -120,10 +127,14 @@ def test_returns_only_real_atoms(self, mixed_structure): def test_ghost_atom_excluded(self): """GhostAtom is a subclass of Atom but should be excluded.""" - structure = Structure(atoms=[ - Atom(element=Element("O"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), - GhostAtom(element=Element("H"), coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0))), - ]) + structure = Structure( + atoms=[ + Atom(element=Element("O"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), + GhostAtom( + element=Element("H"), coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)) + ), + ] + ) assert len(structure.real_atoms) == 1 assert type(structure.real_atoms[0]) is Atom @@ -136,9 +147,9 @@ def test_return_type(self, water): assert all(type(a) is Atom for a in water.real_atoms) def test_empty_structure(self): - structure = Structure(atoms=[ - PointCharge(coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0)), charge=1.0) - ]) + structure = Structure( + atoms=[PointCharge(coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0)), charge=1.0)] + ) assert water.real_atoms == [] if False else structure.real_atoms == [] @@ -146,6 +157,7 @@ def test_empty_structure(self): # get_coordinates # ============================================================ + class TestGetCoordinates: def test_shape(self, water): coords = water.get_coordinates() @@ -184,6 +196,7 @@ def test_single_atom(self): # get_coordinates_at_centroid # ============================================================ + class TestGetCoordinatesAtCentroid: def test_centroid_is_zero(self, water): coords = water.get_coordinates_at_centroid() @@ -208,6 +221,7 @@ def test_translation_invariant(self, water, water_translated): # set_coordinates # ============================================================ + class TestSetCoordinates: def test_returns_new_structure(self, water): new_coords = water.get_coordinates() + 1.0 @@ -243,6 +257,7 @@ def test_multiplicity_preserved(self, water): # centered_structure # ============================================================ + class TestCenteredStructure: def test_centroid_at_origin(self, water): centered = water.centered_structure() @@ -278,6 +293,7 @@ def test_translated_centered_equals_original_centered(self, water, water_transla # _filtered_atoms # ============================================================ + class TestFilteredAtoms: def test_default_returns_only_atom_indices(self, mixed_structure): """PointCharge should be excluded — only index 0 (the Atom) returned.""" @@ -301,10 +317,12 @@ def test_only_atoms_indices(self, ethanol): def test_only_atoms_excludes_pointcharge(self): """Explicitly indexed PointCharge must be excluded by the final type check.""" - structure = Structure(atoms=[ - Atom(element=Element("C"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), - PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), - ]) + structure = Structure( + atoms=[ + Atom(element=Element("C"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), + PointCharge(coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0)), charge=1.0), + ] + ) indices = structure._filtered_atoms([0, 1], False) assert indices == [0] assert all(type(structure.atoms[i]) is Atom for i in indices) @@ -320,6 +338,7 @@ def test_only_atoms_takes_priority_over_ignore_hs(self, ethanol): # _validate_rmsd_compatibility # ============================================================ + class TestValidateRmsdCompatibility: def test_compatible_structures(self, water): atoms = water.real_atoms @@ -347,6 +366,7 @@ def test_error_message_uses_natural_counting(self): # rmsd # ============================================================ + class TestRmsd: def test_identical_structures_zero_rmsd(self, water): assert pytest.approx(water.rmsd(water), abs=1e-10) == 0.0 @@ -398,6 +418,7 @@ def test_result_is_float(self, water): # rmsd_kabsch # ============================================================ + class TestRmsdKabsch: def test_identical_structures_zero_rmsd(self, water): assert pytest.approx(water.rmsd_kabsch(water), abs=1e-10) == 0.0 From 32236cedf1c4abf8c5bd5d4636b7947b7011ecab Mon Sep 17 00:00:00 2001 From: crjacinto Date: Wed, 27 May 2026 14:25:18 +0200 Subject: [PATCH 15/27] Modified set_coordinates() and centered_structure(). Added relative tolerance for rotor classification. Improved unit tests --- src/opi/input/structures/structure.py | 23 ++--- src/opi/utils/rotconst.py | 49 +++++++++- tests/unit/test_rmsd.py | 131 +++++++++++++------------- tests/unit/test_rotconst.py | 115 +++++++++++----------- 4 files changed, 179 insertions(+), 139 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 40985793..3007b330 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1104,36 +1104,27 @@ def get_coordinates( atom_list = [self.atoms[i] for i in only_atoms] if only_atoms else self.atoms return np.array([a.coordinates.coordinates for a in atom_list]) - def set_coordinates(self, coords: npt.NDArray[np.float64]) -> "Structure": + def set_coordinates(self, coords: npt.NDArray[np.float64]) -> None: """ - Return a new Structure with updated coordinates. + Update the coordinates of all atoms in place. Parameters ---------- coords : npt.NDArray[np.float64], shape (N, 3) New coordinates for all atoms in the same order as `self.atoms`. - Returns - ------- - Structure - New `Structure` instance with replaced coordinates; all other - attributes (charge, multiplicity, origin) are preserved. - Raises ------ ValueError If *coords* shape does not match the number of atoms. """ - coords = np.asarray(coords, dtype=np.float64) if coords.shape != (len(self.atoms), 3): raise ValueError( f"coords shape {coords.shape} does not match expected ({len(self.atoms)}, 3)" ) - new_structure = copy.deepcopy(self) - for atom, new_coord in zip(new_structure.atoms, coords): + for atom, new_coord in zip(self.atoms, coords): atom.coordinates = new_coord - return new_structure def centered_structure(self) -> "Structure": """ @@ -1146,9 +1137,11 @@ def centered_structure(self) -> "Structure": Structure New `Structure` with centered coordinates. """ - real_indices = [i for i, a in enumerate(self.atoms) if type(a) is Atom] - centroid = self.get_coordinates(only_atoms=real_indices).mean(axis=0) - return self.set_coordinates(self.get_coordinates() - centroid) + new_structure = copy.deepcopy(self) + real_indices = [i for i, a in enumerate(new_structure.atoms) if type(a) is Atom] + centroid = new_structure.get_coordinates(only_atoms=real_indices).mean(axis=0) + new_structure.set_coordinates(new_structure.get_coordinates() - centroid) + return new_structure def _filtered_atoms( self, diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index 62839fea..bc110777 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -54,7 +54,42 @@ class PrincipalMoments: Ic: float axes: np.ndarray - def rotor_type(self, tol: float = 1e-3) -> RotorType: + # def rotor_type(self, tol: float = 1e-3) -> RotorType: + # Ia, Ib, Ic = self.Ia, self.Ib, self.Ic + # n_zero = sum(m < tol for m in (Ia, Ib, Ic)) + + # if n_zero == 3: + # return RotorType.MONOATOMIC + # if n_zero == 1 and abs(Ib - Ic) < tol: + # return RotorType.LINEAR + # if abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: + # return RotorType.SPHERICAL_TOP + # if abs(Ib - Ic) < tol: # Ia < Ib == Ic → prolate + # return RotorType.PROLATE_TOP + # if abs(Ia - Ib) < tol: # Ia == Ib < Ic → oblate + # return RotorType.OBLATE_TOP + # return RotorType.ASYMMETRIC_TOP + + def rotor_type(self, tol: float = 1e-3, kappa_tol: float = 0.01) -> RotorType: + """ + Classify the molecular rotor from the principal moments of inertia. + + Parameters + ---------- + tol : float, default 1e-3 + Absolute tolerance (amu·Å²) for treating a moment as zero, + used for `MONOATOMIC`, `LINEAR`, and `SPHERICAL_TOP` detection. + kappa_tol : float, default 0.01 + Tolerance for the asymmetry parameter κ (Ray's asymmetry parameter). + κ = -1 → perfect prolate top, κ = +1 → perfect oblate top. + Values within `kappa_tol` of ±1 are classified as symmetric tops. + Reference: Gordy & Cook, *Microwave Molecular Spectra* (1984). + + Returns + ------- + RotorType + Molecular rotor classification. + """ Ia, Ib, Ic = self.Ia, self.Ib, self.Ic n_zero = sum(m < tol for m in (Ia, Ib, Ic)) @@ -64,9 +99,17 @@ def rotor_type(self, tol: float = 1e-3) -> RotorType: return RotorType.LINEAR if abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: return RotorType.SPHERICAL_TOP - if abs(Ib - Ic) < tol: # Ia < Ib == Ic → prolate + + # --- Ray's asymmetry parameter κ --- + # κ = (2B - A - C) / (A - C), where A ≥ B ≥ C (A = 1/Ia, etc.) + # κ = -1 → prolate symmetric top + # κ = +1 → oblate symmetric top + A, B, C = 1.0 / Ia, 1.0 / Ib, 1.0 / Ic + kappa = (2.0 * B - A - C) / (A - C) + + if abs(kappa - (-1.0)) < kappa_tol: return RotorType.PROLATE_TOP - if abs(Ia - Ib) < tol: # Ia == Ib < Ic → oblate + if abs(kappa - 1.0) < kappa_tol: return RotorType.OBLATE_TOP return RotorType.ASYMMETRIC_TOP diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_rmsd.py index 338429bd..e266063c 100644 --- a/tests/unit/test_rmsd.py +++ b/tests/unit/test_rmsd.py @@ -16,6 +16,8 @@ mixed atom types, wrong shapes, identical structures. """ +import copy + import numpy as np import pytest @@ -48,20 +50,22 @@ def make_structure(*atoms: Atom) -> Structure: @pytest.fixture() def water() -> Structure: """H2O at a known geometry.""" - return make_structure( - make_atom("O", 0.000000, 0.000000, 0.119748), - make_atom("H", 0.000000, 0.756950, -0.478993), - make_atom("H", 0.000000, -0.756950, -0.478993), + return Structure.from_xyz_block( + "3\n\n" + "O 0.000000 0.000000 0.119748\n" + "H 0.000000 0.756950 -0.478993\n" + "H 0.000000 -0.756950 -0.478993" ) @pytest.fixture() def water_translated() -> Structure: """H2O shifted by (1, 2, 3) — RMSD vs water should be 0 after centring.""" - return make_structure( - make_atom("O", 1.000000, 2.000000, 3.119748), - make_atom("H", 1.000000, 2.756950, 2.521007), - make_atom("H", 1.000000, 1.243050, 2.521007), + return Structure.from_xyz_block( + "3\n\n" + "O 1.000000 2.000000 3.119748\n" + "H 1.000000 2.756950 2.521007\n" + "H 1.000000 1.243050 2.521007" ) @@ -81,25 +85,28 @@ def water_rotated(water) -> Structure: [0, 0, 1], ] ) - return centered.set_coordinates(coords @ R.T) + rotated = copy.deepcopy(centered) + rotated.set_coordinates(coords @ R.T) + return rotated @pytest.fixture() def ethanol() -> Structure: """ - Simple ethanol-like structure (C2H5OH) for ignore_hs tests. + Ethanol (C2H5OH) for ignore_hs tests. Coordinates are approximate — correctness of geometry is not critical here. """ - return make_structure( - make_atom("C", 0.000, 0.000, 0.000), - make_atom("C", 1.540, 0.000, 0.000), - make_atom("O", 2.060, 1.190, 0.000), - make_atom("H", -0.390, 1.020, 0.000), - make_atom("H", -0.390, -0.510, 0.890), - make_atom("H", -0.390, -0.510, -0.890), - make_atom("H", 1.930, -0.510, 0.890), - make_atom("H", 1.930, -0.510, -0.890), - make_atom("H", 2.980, 1.190, 0.000), + return Structure.from_xyz_block( + "9\n\n" + "C 0.000 0.000 0.000\n" + "C 1.540 0.000 0.000\n" + "O 2.060 1.190 0.000\n" + "H -0.390 1.020 0.000\n" + "H -0.390 -0.510 0.890\n" + "H -0.390 -0.510 -0.890\n" + "H 1.930 -0.510 0.890\n" + "H 1.930 -0.510 -0.890\n" + "H 2.980 1.190 0.000" ) @@ -150,7 +157,7 @@ def test_empty_structure(self): structure = Structure( atoms=[PointCharge(coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0)), charge=1.0)] ) - assert water.real_atoms == [] if False else structure.real_atoms == [] + assert structure.real_atoms == [] # ============================================================ @@ -187,7 +194,7 @@ def test_only_atoms_subset(self, water): assert coords.shape == (2, 3) def test_single_atom(self): - s = make_structure(make_atom("C", 1.0, 2.0, 3.0)) + s = Structure.from_xyz_block("1\n\nC 1.0 2.0 3.0") coords = s.get_coordinates() np.testing.assert_allclose(coords[0], [1.0, 2.0, 3.0]) @@ -223,34 +230,27 @@ def test_translation_invariant(self, water, water_translated): class TestSetCoordinates: - def test_returns_new_structure(self, water): - new_coords = water.get_coordinates() + 1.0 - new_s = water.set_coordinates(new_coords) - assert new_s is not water - - def test_original_unchanged(self, water): - original_coords = water.get_coordinates().copy() - water.set_coordinates(water.get_coordinates() + 5.0) - np.testing.assert_allclose(water.get_coordinates(), original_coords) - - def test_new_coords_applied(self, water): + def test_mutates_in_place(self, water): + """set_coordinates should modify the structure in place.""" new_coords = np.zeros((3, 3)) - new_s = water.set_coordinates(new_coords) - np.testing.assert_allclose(new_s.get_coordinates(), new_coords) + water.set_coordinates(new_coords) + np.testing.assert_allclose(water.get_coordinates(), new_coords) + + def test_returns_none(self, water): + """set_coordinates should return None.""" + result = water.set_coordinates(water.get_coordinates()) + assert result is None def test_raises_on_wrong_shape(self, water): with pytest.raises(ValueError, match="coords shape"): water.set_coordinates(np.zeros((2, 3))) - def test_charge_preserved(self, water): - water_charged = Structure(atoms=water.atoms, charge=1) - new_s = water_charged.set_coordinates(water.get_coordinates()) - assert new_s.charge == 1 - - def test_multiplicity_preserved(self, water): - water_mult = Structure(atoms=water.atoms, multiplicity=3) - new_s = water_mult.set_coordinates(water.get_coordinates()) - assert new_s.multiplicity == 3 + def test_deepcopy_is_independent(self, water): + """Changes to a deepcopy should not affect the original.""" + original_coords = water.get_coordinates().copy() + new_structure = copy.deepcopy(water) + new_structure.set_coordinates(np.zeros((3, 3))) + np.testing.assert_allclose(water.get_coordinates(), original_coords) # ============================================================ @@ -384,29 +384,31 @@ def test_raises_on_incompatible_structures(self, water, ethanol): def test_ignore_hs(self, ethanol): """ - Displace only heavy atoms — RMSD with and without H should differ - since the centroid shift differs between the two filtered sets. + Displace only one heavy atom — RMSD with and without H should differ. """ - coords = ethanol.get_coordinates() - new_coords = coords.copy() - new_coords[0] += np.array([0.5, 0.0, 0.0]) # shift only C0 - shifted = ethanol.set_coordinates(new_coords) + shifted = copy.deepcopy(ethanol) + coords = shifted.get_coordinates() + coords[0] += np.array([0.5, 0.0, 0.0]) + shifted.set_coordinates(coords) rmsd_all = ethanol.rmsd(shifted) rmsd_no_h = ethanol.rmsd(shifted, ignore_hs=True) assert rmsd_all != pytest.approx(rmsd_no_h, abs=1e-6) def test_only_atoms_subset(self, ethanol): """RMSD over a subset of atoms should differ from the full-molecule RMSD.""" - coords = ethanol.get_coordinates() - new_coords = coords.copy() - new_coords[0] += np.array([0.5, 0.0, 0.0]) # shift only C0 - shifted = ethanol.set_coordinates(new_coords) + shifted = copy.deepcopy(ethanol) + coords = shifted.get_coordinates() + coords[0] += np.array([0.5, 0.0, 0.0]) + shifted.set_coordinates(coords) rmsd_all = ethanol.rmsd(shifted) rmsd_subset = ethanol.rmsd(shifted, only_atoms=[0, 1, 2]) assert rmsd_all != pytest.approx(rmsd_subset, abs=1e-6) def test_nonzero_rmsd_for_different_structures(self, water): - other = water.set_coordinates(water.get_coordinates() + np.array([0, 0, 1.0])) + other = copy.deepcopy(water) + coords = other.get_coordinates() + coords += np.array([0, 0, 1.0]) + other.set_coordinates(coords) result = water.rmsd(other) assert result >= 0.0 @@ -436,7 +438,8 @@ def test_kabsch_le_rmsd(self, water): centered = water.centered_structure() coords = centered.get_coordinates() R = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]], dtype=float) - other = centered.set_coordinates(coords @ R.T) + other = copy.deepcopy(centered) + other.set_coordinates(coords @ R.T) assert centered.rmsd_kabsch(other) <= centered.rmsd(other) + 1e-10 def test_symmetry(self, water, water_rotated): @@ -449,19 +452,19 @@ def test_raises_on_incompatible_structures(self, water, ethanol): water.rmsd_kabsch(ethanol) def test_ignore_hs(self, ethanol): - coords = ethanol.get_coordinates() - new_coords = coords.copy() - new_coords[0] += np.array([0.5, 0.0, 0.0]) - shifted = ethanol.set_coordinates(new_coords) + shifted = copy.deepcopy(ethanol) + coords = shifted.get_coordinates() + coords[0] += np.array([0.5, 0.0, 0.0]) + shifted.set_coordinates(coords) rmsd_all = ethanol.rmsd_kabsch(shifted) rmsd_no_h = ethanol.rmsd_kabsch(shifted, ignore_hs=True) assert rmsd_all != pytest.approx(rmsd_no_h, abs=1e-6) def test_only_atoms_subset(self, ethanol): - coords = ethanol.get_coordinates() - new_coords = coords.copy() - new_coords[0] += np.array([0.5, 0.0, 0.0]) - shifted = ethanol.set_coordinates(new_coords) + shifted = copy.deepcopy(ethanol) + coords = shifted.get_coordinates() + coords[0] += np.array([0.5, 0.0, 0.0]) + shifted.set_coordinates(coords) rmsd_all = ethanol.rmsd_kabsch(shifted) rmsd_subset = ethanol.rmsd_kabsch(shifted, only_atoms=[0, 1, 2]) assert rmsd_all != pytest.approx(rmsd_subset, abs=1e-6) diff --git a/tests/unit/test_rotconst.py b/tests/unit/test_rotconst.py index 657c38eb..9f4db8a9 100644 --- a/tests/unit/test_rotconst.py +++ b/tests/unit/test_rotconst.py @@ -14,8 +14,6 @@ - Edge cases: empty structure, all-zero masses, mass overrides, unknown elements """ -import math - import numpy as np import pytest @@ -49,10 +47,6 @@ def make_no_real_atoms_structure() -> Structure: ) -def make_structure(*atoms: Atom) -> Structure: - return Structure(atoms=list(atoms)) - - # ============================================================ # Fixtures # ============================================================ @@ -61,10 +55,11 @@ def make_structure(*atoms: Atom) -> Structure: @pytest.fixture() def water() -> Structure: """H2O — asymmetric top.""" - return make_structure( - make_atom("O", 0.000000, 0.000000, 0.119748), - make_atom("H", 0.000000, 0.756950, -0.478993), - make_atom("H", 0.000000, -0.756950, -0.478993), + return Structure.from_xyz_block( + "3\n\n" + "O 0.000000 0.000000 0.119748\n" + "H 0.000000 0.756950 -0.478993\n" + "H 0.000000 -0.756950 -0.478993" ) @@ -73,75 +68,85 @@ def co2() -> Structure: """ CO2 — linear molecule aligned exactly on Z. After COM shift, both O atoms lie on the Z axis so Ia = 0 exactly. - Using exact masses (C=12, O=15.999) the COM is at z=0 by symmetry. """ - return make_structure( - make_atom("C", 0.000000, 0.000000, 0.000000), - make_atom("O", 0.000000, 0.000000, 1.160000), - make_atom("O", 0.000000, 0.000000, -1.160000), + return Structure.from_xyz_block( + "3\n\n" + "C 0.000000 0.000000 0.000000\n" + "O 0.000000 0.000000 1.160000\n" + "O 0.000000 0.000000 -1.160000" ) @pytest.fixture() def hcl() -> Structure: """HCl — guaranteed linear (diatomic), Ia = 0 exactly on Z axis.""" - return make_structure( - make_atom("H", 0.000000, 0.000000, 0.000000), - make_atom("Cl", 0.000000, 0.000000, 1.274500), + return Structure.from_xyz_block( + "2\n\nH 0.000000 0.000000 0.000000\nCl 0.000000 0.000000 1.274500" ) @pytest.fixture() def methane() -> Structure: """CH4 — spherical top.""" - d = 0.6276 # C-H bond / sqrt(3) - return make_structure( - make_atom("C", 0.000, 0.000, 0.000), - make_atom("H", d, d, d), - make_atom("H", -d, -d, d), - make_atom("H", -d, d, -d), - make_atom("H", d, -d, -d), + d = 0.6276 + return Structure.from_xyz_block( + f"5\n\n" + f"C 0.000 0.000 0.000\n" + f"H {d} {d} {d}\n" + f"H -{d} -{d} {d}\n" + f"H -{d} {d} -{d}\n" + f"H {d} -{d} -{d}" ) @pytest.fixture() def benzene() -> Structure: - """C6H6 — oblate symmetric top.""" - r_c, r_h = 1.3970, 2.4832 - atoms = [] - for i in range(6): - angle = math.radians(i * 60) - atoms.append(make_atom("C", r_c * math.cos(angle), r_c * math.sin(angle), 0.0)) - atoms.append(make_atom("H", r_h * math.cos(angle), r_h * math.sin(angle), 0.0)) - return make_structure(*atoms) + """C6H6 — oblate symmetric top. Avogadro optimized geometry.""" + return Structure.from_xyz_block( + "12\n\n" + "C -4.0132333739 3.1874756422 0.0000029092\n" + "C -3.7419050401 1.8193724296 0.0000182442\n" + "H -0.3432332788 1.9395893312 0.0002037711\n" + "C -2.9641082474 4.1065582543 0.0000491840\n" + "C -2.4213601440 1.3702446325 0.0000859970\n" + "H -5.0421713998 3.5372885217 -0.0000468671\n" + "C -1.6435439604 3.6575023283 0.0001222947\n" + "H -3.1755373879 5.1725212691 0.0000321721\n" + "C -1.3721905089 2.2894069399 0.0001431310\n" + "H -4.5593921780 1.1033282205 -0.0000191331\n" + "H -0.8260945740 4.3735825219 0.0001644962\n" + "H -2.2100493877 0.3043042675 0.0000939114" + ) @pytest.fixture() def single_atom() -> Structure: """Single atom — monoatomic.""" - return make_structure(make_atom("C", 0.0, 0.0, 0.0)) + return Structure.from_xyz_block("1\n\nC 0.000000 0.000000 0.000000") @pytest.fixture() def ammonia() -> Structure: """NH3 — oblate symmetric top (C3v). XTB2 optimized geometry.""" - return make_structure( - make_atom("N", -0.000000, -0.000066, 0.100407), - make_atom("H", 0.000000, 0.943825, -0.266468), - make_atom("H", 0.817493, -0.471879, -0.266439), - make_atom("H", -0.817493, -0.471879, -0.266439), + return Structure.from_xyz_block( + "4\n\n" + "N -0.000000 -0.000066 0.100407\n" + "H 0.000000 0.943825 -0.266468\n" + "H 0.817493 -0.471879 -0.266439\n" + "H -0.817493 -0.471879 -0.266439" ) @pytest.fixture() def ch3cl() -> Structure: """CH3Cl — prolate symmetric top (C3v). XTB2 optimized geometry.""" - return make_structure( - make_atom("Cl", 0.9754830000, 0.0921220000, -0.0239260000), - make_atom("C", 2.7424550000, 0.0921300000, -0.0239350000), - make_atom("H", 3.0994960000, 0.3238830000, 0.9818970000), - make_atom("H", 3.0994920000, -0.8948230000, -0.3261570000), - make_atom("H", 3.0994750000, 0.8473380000, -0.7275610000), + return Structure.from_xyz_block( + "5\n\n" + "Cl 0.975483 0.092122 -0.023926\n" + "C 2.742455 0.092130 -0.023935\n" + "H 3.099496 0.323883 0.981897\n" + "H 3.099492 -0.894823 -0.326157\n" + "H 3.099475 0.847338 -0.727561" ) @@ -272,10 +277,11 @@ def test_raises_on_wrong_masses_length(self, water): def test_ghost_atoms_ignored(self): """PointCharge should not contribute to moments.""" - real = make_structure( - make_atom("O", 0.0, 0.0, 0.119748), - make_atom("H", 0.0, 0.756950, -0.478993), - make_atom("H", 0.0, -0.756950, -0.478993), + real = Structure.from_xyz_block( + "3\n\n" + "O 0.000000 0.000000 0.119748\n" + "H 0.000000 0.756950 -0.478993\n" + "H 0.000000 -0.756950 -0.478993" ) with_pc = Structure( atoms=[ @@ -336,14 +342,11 @@ def test_unknown_element_warns_and_excludes(self): Atom(element=Element("H"), coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0))), ] ) - # Patch via weights to simulate unknown: pass mass 0 explicitly for C pm = structure.calc_moments_of_inertia(atom_weights={0: 0.0}) - # Only H contributes → single non-zero mass → linear or monoatomic assert pm is not None def test_zero_mass_atoms_filtered(self, water): """Atoms with zero mass are excluded; remaining atoms determine the rotor type.""" - # Zero out the two H atoms → only O remains → monoatomic masses = np.array([15.999, 0.0, 0.0]) pm = water.calc_moments_of_inertia(masses=masses) assert pm is not None @@ -352,15 +355,13 @@ def test_zero_mass_atoms_filtered(self, water): def test_water_moments_known_values(self, water): """ Check moments are self-consistent: Ia <= Ib <= Ic and all positive. - The geometry used here gives Ia ≈ 0.6418 amu·Å² (not the experimental - equilibrium value of 0.5791) because the O-H distance and angle differ - slightly. We verify internal consistency rather than an absolute value. + We verify internal consistency and a physically reasonable range + rather than absolute values. """ pm = water.calc_moments_of_inertia() assert pm is not None assert pm.Ia > 0.0 assert pm.Ia <= pm.Ib <= pm.Ic - # Rough sanity check: moments should be in a physically reasonable range assert 0.1 < pm.Ia < 5.0 assert 0.1 < pm.Ib < 5.0 assert 0.1 < pm.Ic < 5.0 @@ -471,4 +472,4 @@ def test_a_not_in_sync_problem(self, water): wn_before = rc.get_in_wavenumbers()[0] rc.A = rc.A * 5 if rc.A is not None else None # type: ignore[operator] wn_after = rc.get_in_wavenumbers()[0] - assert wn_before != wn_after # confirms on-the-fly computation + assert wn_before != wn_after From 723dfc9b451bcc62ab13ed7e3988254dab8db388 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Wed, 10 Jun 2026 14:30:00 +0200 Subject: [PATCH 16/27] Modified calc_rotataional_constants(), added private function _resolve_masses(), and updated unit tests --- src/opi/input/structures/structure.py | 147 +++++++++++++------------- src/opi/utils/rotconst.py | 16 --- tests/unit/test_rotconst.py | 75 +++++++------ 3 files changed, 112 insertions(+), 126 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 3007b330..edf06872 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1331,86 +1331,91 @@ def rmsd_kabsch( # ------------------------------------------------------------------ # # Moment of inertia # # ------------------------------------------------------------------ # - def calc_moments_of_inertia( + def _resolve_masses( self, - masses: npt.NDArray[np.float64] | None = None, - weights: dict[str, float] | None = None, - atom_weights: dict[int, float] | None = None, - ) -> PrincipalMoments | None: + elem_masses: dict[str, float] | None = None, + ) -> npt.NDArray[np.float64]: """ - Compute the principal axes and moments of inertia for this structure. + Resolve masses for all real `Atom` instances without mutating any atom. + + Mass priority + ------------- + atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) Parameters ---------- - masses : npt.NDArray[np.float64] | None, default: None - Per-atom masses (amu) for all Atom instances, overriding every - other mass source. - atom_weights : dict[int, float] | None, default: None - Per-atom mass overrides keyed by index within the Atom list. - weights : dict[str, float] | None, default: None + elem_masses : dict[str, float] | None, default None Per-element mass overrides keyed by element symbol string - (e.g. `{"C": 13.003}`). + (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. + + Returns + ------- + npt.NDArray[np.float64], shape (N,) + Masses in amu in the same order as `self.real_atoms`. + """ + elem_masses = elem_masses or {} + masses_list: list[float] = [] + for atom in self.real_atoms: + if atom.mass is not None: + # Per-atom mass set directly on the atom takes highest priority + m = atom.mass + elif atom.element.value in elem_masses: + # Per-element override takes second priority + m = elem_masses[atom.element.value] + elif atom.element in ATOMIC_MASSES_FROM_ELEMENT: + # Fall back to default atomic masses + m = ATOMIC_MASSES_FROM_ELEMENT[atom.element] + else: + # Unknown element — warn and exclude from computation + warn(f"Unknown element '{atom.element.value}' → mass set to 0.0") + m = 0.0 + masses_list.append(m) + return np.array(masses_list, dtype=np.float64) + + def calc_moments_of_inertia( + self, + elem_masses: dict[str, float] | None = None, + ) -> PrincipalMoments | None: + """ + Compute the principal axes and moments of inertia for this structure. Only `Atom` instances contribute; `GhostAtom`, `PointCharge`, - and `EmbeddingPotential` atoms are silently ignored. + and `EmbeddingPotential` atoms are silently ignored. Note that + `GhostAtom` is a subclass of `Atom`, so an exact type check + (`type(a) is Atom`) is used rather than `isinstance`. Mass priority ------------- - masses > atom_weights > weights > default: (ATOMIC_MASSES_FROM_ELEMENT) + atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) + + Parameters + ---------- + elem_masses : dict[str, float] | None, default None + Per-element mass overrides keyed by element symbol string + (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. Returns ------- PrincipalMoments | None Principal axes and moments of inertia in amu·Å², sorted ascending. Returns `None` if no `Atom` instances are present or all masses are zero. - - Raises - ------ - ValueError - If the length of *masses* does not match the number of `Atom` - instances in the structure. """ atom_indices = [i for i, a in enumerate(self.atoms) if type(a) is Atom] if not atom_indices: return None coords = self.get_coordinates(only_atoms=atom_indices) - # --- Mass assignment priority: masses > atom_weights > weights > default --- - # Normalise the optional dicts to avoid None checks below - weights = {k: v for k, v in (weights or {}).items()} - atom_weights = atom_weights or {} - - if masses is not None: - masses = np.asarray(masses, dtype=np.float64) - if len(masses) != len(self.real_atoms): - raise ValueError( - f"masses length ({len(masses)}) does not match " - f"number of atoms ({len(self.real_atoms)})" - ) - else: - masses_list: list[float] = [] - for i, atom in enumerate(self.real_atoms): - if i in atom_weights: - m = atom_weights[i] - elif atom.element.value in weights: - m = weights[atom.element.value] - elif atom.element in ATOMIC_MASSES_FROM_ELEMENT: - m = ATOMIC_MASSES_FROM_ELEMENT[atom.element] - else: - warn(f"Unknown element '{atom.element.value}' → mass set to 0.0") - m = 0.0 - masses_list.append(m) - masses = np.array(masses_list, dtype=np.float64) + masses = self._resolve_masses(elem_masses) + # --- Filter out zero-mass atoms before computing the inertia tensor --- + # This handles unknown elements that were assigned 0.0 mask = masses > 0.0 if not np.any(mask): return None - # --- Filter out zero-mass atoms before computing the inertia tensor --- - # This handles dummy atoms or unknown elements that were assigned 0.0 masses = masses[mask] coords = coords[mask] - total_mass = float(masses.sum()) + total_mass = float(masses.sum()) com = (masses[:, None] * coords).sum(axis=0) / total_mass coords -= com @@ -1422,12 +1427,12 @@ def calc_moments_of_inertia( moments_raw = np.maximum(moments_raw, 0.0) return PrincipalMoments( - Ia=float(moments_raw[0]), Ib=float(moments_raw[1]), Ic=float(moments_raw[2]), axes=axes + Ia=float(moments_raw[0]), + Ib=float(moments_raw[1]), + Ic=float(moments_raw[2]), + axes=axes, ) - # ------------------------------------------------------------------ # - # Rotor classification # - # ------------------------------------------------------------------ # def calc_rotor_type( self, moments: PrincipalMoments | None = None, @@ -1438,13 +1443,13 @@ def calc_rotor_type( Parameters ---------- - moments : ndarray shape (3,) | None, default None - Pre-computed principal moments (amu·Å², ascending). When + moments : PrincipalMoments | None, default None + Pre-computed principal moments (amu·Å², ascending). When `None` the moments are computed via - :meth:`calc_moment_of_inertia` using *mass_kwargs*. + :meth:`calc_moments_of_inertia` using *mass_kwargs*. **mass_kwargs - Forwarded to :meth:`calc_moment_of_inertia` when *moments* is - `None` (i.e. `masses`, `weights`, `atom_weights`). + Forwarded to :meth:`calc_moments_of_inertia` when *moments* is + `None` (i.e. `elem_masses`). Returns ------- @@ -1456,17 +1461,11 @@ def calc_rotor_type( if result is None: return None moments = result - return moments.rotor_type() - # ------------------------------------------------------------------ # - # Rotational constants # - # ------------------------------------------------------------------ # def calc_rotational_constants( self, - masses: npt.NDArray[np.float64] | None = None, - weights: dict[str, float] | None = None, - atom_weights: dict[int, float] | None = None, + elem_masses: dict[str, float] | None = None, ) -> RotationalConstants | None: """ Compute rotational constants for this structure. @@ -1476,24 +1475,20 @@ def calc_rotational_constants( Mass priority ------------- - masses > atom_weights > weights > default (ATOMIC_MASSES_FROM_ELEMENT) + elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) Parameters ---------- - masses : npt.NDArray[np.float64] | None, default None - Per-atom masses (amu) for all Atom instances. - atom_weights : dict[int, float] | None, default None - Per-atom mass overrides keyed by index within the Atom list. - weights : dict[str, float] | None, default None - Per-element mass overrides (e.g. `{"C": 13.003}`). + elem_masses : dict[str, float] | None, default None + Per-element mass overrides keyed by element symbol string + (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. Returns ------- RotationalConstants | None - `None` if no `Atom` instances are present or all masses are - zero. + `None` if no `Atom` instances are present or all masses are zero. """ - pm = self.calc_moments_of_inertia(masses=masses, weights=weights, atom_weights=atom_weights) + pm = self.calc_moments_of_inertia(elem_masses=elem_masses) if pm is None: return None A = moment_to_mhz(pm.Ia) diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index bc110777..cb92a59b 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -54,22 +54,6 @@ class PrincipalMoments: Ic: float axes: np.ndarray - # def rotor_type(self, tol: float = 1e-3) -> RotorType: - # Ia, Ib, Ic = self.Ia, self.Ib, self.Ic - # n_zero = sum(m < tol for m in (Ia, Ib, Ic)) - - # if n_zero == 3: - # return RotorType.MONOATOMIC - # if n_zero == 1 and abs(Ib - Ic) < tol: - # return RotorType.LINEAR - # if abs(Ia - Ib) < tol and abs(Ib - Ic) < tol: - # return RotorType.SPHERICAL_TOP - # if abs(Ib - Ic) < tol: # Ia < Ib == Ic → prolate - # return RotorType.PROLATE_TOP - # if abs(Ia - Ib) < tol: # Ia == Ib < Ic → oblate - # return RotorType.OBLATE_TOP - # return RotorType.ASYMMETRIC_TOP - def rotor_type(self, tol: float = 1e-3, kappa_tol: float = 0.01) -> RotorType: """ Classify the molecular rotor from the principal moments of inertia. diff --git a/tests/unit/test_rotconst.py b/tests/unit/test_rotconst.py index 9f4db8a9..5a43601b 100644 --- a/tests/unit/test_rotconst.py +++ b/tests/unit/test_rotconst.py @@ -33,10 +33,11 @@ # ============================================================ -def make_atom(element: str, x: float, y: float, z: float) -> Atom: +def make_atom(element: str, x: float, y: float, z: float, mass: float | None = None) -> Atom: return Atom( element=Element(element), coordinates=Coordinates(coordinates=(x, y, z)), + mass=mass, ) @@ -103,7 +104,8 @@ def methane() -> Structure: def benzene() -> Structure: """C6H6 — oblate symmetric top. Avogadro optimized geometry.""" return Structure.from_xyz_block( - "12\n\n" + "12\n" + "XYZ file generated by Avogadro.\n" "C -4.0132333739 3.1874756422 0.0000029092\n" "C -3.7419050401 1.8193724296 0.0000182442\n" "H -0.3432332788 1.9395893312 0.0002037711\n" @@ -268,12 +270,10 @@ def test_returns_none_for_empty_structure(self): assert make_no_real_atoms_structure().calc_moments_of_inertia() is None def test_returns_none_for_all_zero_masses(self, water): - masses = np.zeros(3) - assert water.calc_moments_of_inertia(masses=masses) is None - - def test_raises_on_wrong_masses_length(self, water): - with pytest.raises(ValueError, match="masses length"): - water.calc_moments_of_inertia(masses=np.array([1.0, 2.0])) + """Setting mass=0.0 on all atoms should return None.""" + for atom in water.real_atoms: + atom.mass = 0.0 + assert water.calc_moments_of_inertia() is None def test_ghost_atoms_ignored(self): """PointCharge should not contribute to moments.""" @@ -306,33 +306,29 @@ def test_ghost_atoms_ignored(self): assert pytest.approx(pm_real.Ib, rel=1e-6) == pm_pc.Ib assert pytest.approx(pm_real.Ic, rel=1e-6) == pm_pc.Ic - def test_custom_masses_override(self, water): - """Passing explicit masses should change the moments.""" + def test_atom_mass_override(self, water): + """Setting `atom.mass` directly should change the moments.""" pm_default = water.calc_moments_of_inertia() - pm_custom = water.calc_moments_of_inertia(masses=np.array([18.0, 2.0, 2.0])) + water.real_atoms[0].mass = 17.999 # override O mass + pm_custom = water.calc_moments_of_inertia() assert pm_default is not None and pm_custom is not None - assert pm_default.Ia != pytest.approx(pm_custom.Ia, rel=1e-3) - - def test_weights_per_element(self, water): - """Per-element weight override should change the moments.""" - pm_default = water.calc_moments_of_inertia() - pm_weights = water.calc_moments_of_inertia(weights={"H": 2.014}) - assert pm_default is not None and pm_weights is not None - assert pm_default.Ic != pytest.approx(pm_weights.Ic, rel=1e-3) + assert pm_default.Ic != pytest.approx(pm_custom.Ic, rel=1e-3) - def test_atom_weights_override(self, water): - """Per-atom weight override (index 0) should change the moments.""" + def test_elem_masses_override(self, water): + """Per-element override via `elem_masses` should change the moments.""" pm_default = water.calc_moments_of_inertia() - pm_aw = water.calc_moments_of_inertia(atom_weights={0: 17.999}) - assert pm_default is not None and pm_aw is not None - assert pm_default.Ic != pytest.approx(pm_aw.Ic, rel=1e-3) - - def test_atom_weights_priority_over_weights(self, water): - """atom_weights should take priority over weights for the same atom.""" - pm_aw = water.calc_moments_of_inertia(atom_weights={0: 17.999}) - pm_both = water.calc_moments_of_inertia(weights={"O": 16.0}, atom_weights={0: 17.999}) - assert pm_aw is not None and pm_both is not None - assert pytest.approx(pm_aw.Ic, rel=1e-9) == pm_both.Ic + pm_elem = water.calc_moments_of_inertia(elem_masses={"H": 2.014}) + assert pm_default is not None and pm_elem is not None + assert pm_default.Ic != pytest.approx(pm_elem.Ic, rel=1e-3) + + def test_atom_mass_priority_over_elem_masses(self, water): + """`atom.mass` should take priority over `elem_masses` for the same atom.""" + water.real_atoms[0].mass = 17.999 # set O mass directly + pm_atom = water.calc_moments_of_inertia() + pm_both = water.calc_moments_of_inertia(elem_masses={"O": 16.0}) + assert pm_atom is not None and pm_both is not None + # elem_masses={"O": 16.0} should be ignored since atom.mass is set + assert pytest.approx(pm_atom.Ic, rel=1e-9) == pm_both.Ic def test_unknown_element_warns_and_excludes(self): """Atoms with unknown elements should warn and be assigned mass 0.""" @@ -342,13 +338,17 @@ def test_unknown_element_warns_and_excludes(self): Atom(element=Element("H"), coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0))), ] ) - pm = structure.calc_moments_of_inertia(atom_weights={0: 0.0}) + # Force C mass to 0 via atom.mass + structure.real_atoms[0].mass = 0.0 + pm = structure.calc_moments_of_inertia() assert pm is not None def test_zero_mass_atoms_filtered(self, water): """Atoms with zero mass are excluded; remaining atoms determine the rotor type.""" - masses = np.array([15.999, 0.0, 0.0]) - pm = water.calc_moments_of_inertia(masses=masses) + # Set H masses to 0 — only O remains → monoatomic + water.real_atoms[1].mass = 0.0 + water.real_atoms[2].mass = 0.0 + pm = water.calc_moments_of_inertia() assert pm is not None assert pm.rotor_type() == RotorType.MONOATOMIC @@ -473,3 +473,10 @@ def test_a_not_in_sync_problem(self, water): rc.A = rc.A * 5 if rc.A is not None else None # type: ignore[operator] wn_after = rc.get_in_wavenumbers()[0] assert wn_before != wn_after + + def test_elem_masses_forwarded(self, water): + """elem_masses passed to calc_rotational_constants should affect the result.""" + rc_default = water.calc_rotational_constants() + rc_deuterated = water.calc_rotational_constants(elem_masses={"H": 2.014}) + assert rc_default is not None and rc_deuterated is not None + assert rc_default.C != pytest.approx(rc_deuterated.C, rel=1e-3) From 1c8fe240d299865289a104411a9747567562ba35 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Wed, 10 Jun 2026 14:43:19 +0200 Subject: [PATCH 17/27] Modified mass priority: elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) --- src/opi/input/structures/structure.py | 14 +++++++------- tests/unit/test_rotconst.py | 15 +++++++-------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index edf06872..a2d47d0a 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1340,7 +1340,7 @@ def _resolve_masses( Mass priority ------------- - atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) + elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) Parameters ---------- @@ -1356,12 +1356,12 @@ def _resolve_masses( elem_masses = elem_masses or {} masses_list: list[float] = [] for atom in self.real_atoms: - if atom.mass is not None: - # Per-atom mass set directly on the atom takes highest priority - m = atom.mass - elif atom.element.value in elem_masses: - # Per-element override takes second priority + if atom.element.value in elem_masses: + # Per-element override takes highest priority m = elem_masses[atom.element.value] + elif atom.mass is not None: + # Mass set directly on the atom takes second priority + m = atom.mass elif atom.element in ATOMIC_MASSES_FROM_ELEMENT: # Fall back to default atomic masses m = ATOMIC_MASSES_FROM_ELEMENT[atom.element] @@ -1386,7 +1386,7 @@ def calc_moments_of_inertia( Mass priority ------------- - atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) + elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) Parameters ---------- diff --git a/tests/unit/test_rotconst.py b/tests/unit/test_rotconst.py index 5a43601b..7e48d50b 100644 --- a/tests/unit/test_rotconst.py +++ b/tests/unit/test_rotconst.py @@ -321,14 +321,13 @@ def test_elem_masses_override(self, water): assert pm_default is not None and pm_elem is not None assert pm_default.Ic != pytest.approx(pm_elem.Ic, rel=1e-3) - def test_atom_mass_priority_over_elem_masses(self, water): - """`atom.mass` should take priority over `elem_masses` for the same atom.""" - water.real_atoms[0].mass = 17.999 # set O mass directly - pm_atom = water.calc_moments_of_inertia() - pm_both = water.calc_moments_of_inertia(elem_masses={"O": 16.0}) - assert pm_atom is not None and pm_both is not None - # elem_masses={"O": 16.0} should be ignored since atom.mass is set - assert pytest.approx(pm_atom.Ic, rel=1e-9) == pm_both.Ic + def test_elem_masses_priority_over_atom_mass(self, water): + """`elem_masses` should take priority over `atom.mass` for the same atom.""" + water.real_atoms[0].mass = 17.999 + pm_atom_mass = water.calc_moments_of_inertia() + pm_elem = water.calc_moments_of_inertia(elem_masses={"O": 16.0}) + assert pm_atom_mass is not None and pm_elem is not None + assert pm_atom_mass.Ic != pytest.approx(pm_elem.Ic, rel=1e-3) def test_unknown_element_warns_and_excludes(self): """Atoms with unknown elements should warn and be assigned mass 0.""" From 293c7c8a14576eaac5fc06d67f35ee5987aa1168 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Mon, 15 Jun 2026 11:53:54 +0200 Subject: [PATCH 18/27] Modified structure.py modules: RMSD and RotConst. Updated corresponding unit tests. --- src/opi/input/structures/structure.py | 74 +++++++++++++++++---------- src/opi/utils/rotconst.py | 2 +- tests/unit/test_rmsd.py | 1 + tests/unit/test_rotconst.py | 20 ++++---- 4 files changed, 60 insertions(+), 37 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index a2d47d0a..0cafd4bd 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1074,11 +1074,15 @@ def _iter_xyz_structures( if n_struc_limit and n_struc >= n_struc_limit: break + # ------------------------------------------------------------------ # + # RMSD # + # ------------------------------------------------------------------ # + @property def real_atoms(self) -> list[Atom]: """ Return only real `Atom` instances, excluding `GhostAtom`, `PointCharge`, - and `EmbeddingPotential`. Note that these are subclasses of `Atom`, + and `EmbeddingPotential`. Note that `GhostAtom` is a subclass of `Atom`, so `type(a) is Atom` is used rather than `isinstance`. """ return [a for a in self.atoms if type(a) is Atom] @@ -1102,7 +1106,7 @@ def get_coordinates( Coordinates in the same order as the selected atoms. """ atom_list = [self.atoms[i] for i in only_atoms] if only_atoms else self.atoms - return np.array([a.coordinates.coordinates for a in atom_list]) + return np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) def set_coordinates(self, coords: npt.NDArray[np.float64]) -> None: """ @@ -1170,7 +1174,7 @@ def _filtered_atoms( @staticmethod def _validate_rmsd_compatibility(atoms1: list[Atom], atoms2: list[Atom]) -> None: """ - Raise ValueError if atoms1 and atoms2 differ in count or element order. + Raise `ValueError` if `atoms1` and `atoms2` differ in count or element order. """ if len(atoms1) != len(atoms2): raise ValueError( @@ -1224,6 +1228,8 @@ def rmsd( Compute RMSD between this structure and *other* (no rotational alignment). Both structures are translated to their centroid before comparison. + Neither `self` nor *other* is modified; all coordinate operations are + performed on temporary arrays. Parameters ---------- @@ -1268,6 +1274,8 @@ def rmsd_kabsch( Translates both structures to their centroid, then finds the optimal rotation matrix via SVD before computing RMSD. + Neither `self` nor *other* is modified; all coordinate operations are + performed on temporary arrays. Parameters ---------- @@ -1331,9 +1339,10 @@ def rmsd_kabsch( # ------------------------------------------------------------------ # # Moment of inertia # # ------------------------------------------------------------------ # + def _resolve_masses( self, - elem_masses: dict[str, float] | None = None, + elem_masses: dict[str | Element, float] | None = None, ) -> npt.NDArray[np.float64]: """ Resolve masses for all real `Atom` instances without mutating any atom. @@ -1344,9 +1353,10 @@ def _resolve_masses( Parameters ---------- - elem_masses : dict[str, float] | None, default None - Per-element mass overrides keyed by element symbol string - (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. + elem_masses : dict[str | Element, float] | None, default None + Per-element mass overrides keyed by element symbol string or + `Element` instance (e.g. `{"C": 13.003}`). If `None`, default + atomic masses are used. Returns ------- @@ -1354,11 +1364,16 @@ def _resolve_masses( Masses in amu in the same order as `self.real_atoms`. """ elem_masses = elem_masses or {} + # Normalise keys to Element for consistent lookup + normalised: dict[Element, float] = {} + for k, v in elem_masses.items(): + normalised[Element(k) if isinstance(k, str) else k] = v + masses_list: list[float] = [] for atom in self.real_atoms: - if atom.element.value in elem_masses: + if atom.element in normalised: # Per-element override takes highest priority - m = elem_masses[atom.element.value] + m = normalised[atom.element] elif atom.mass is not None: # Mass set directly on the atom takes second priority m = atom.mass @@ -1366,15 +1381,16 @@ def _resolve_masses( # Fall back to default atomic masses m = ATOMIC_MASSES_FROM_ELEMENT[atom.element] else: - # Unknown element — warn and exclude from computation - warn(f"Unknown element '{atom.element.value}' → mass set to 0.0") - m = 0.0 + raise ValueError( + f"Unknown element '{atom.element.value}': no mass available. " + "Provide a mass via `elem_masses` or set `atom.mass` directly." + ) masses_list.append(m) return np.array(masses_list, dtype=np.float64) def calc_moments_of_inertia( self, - elem_masses: dict[str, float] | None = None, + elem_masses: dict[str | Element, float] | None = None, ) -> PrincipalMoments | None: """ Compute the principal axes and moments of inertia for this structure. @@ -1390,9 +1406,10 @@ def calc_moments_of_inertia( Parameters ---------- - elem_masses : dict[str, float] | None, default None - Per-element mass overrides keyed by element symbol string - (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. + elem_masses : dict[str | Element, float] | None, default None + Per-element mass overrides keyed by element symbol string or + `Element` instance (e.g. `{"C": 13.003}`). If `None`, default + atomic masses are used. Returns ------- @@ -1408,7 +1425,7 @@ def calc_moments_of_inertia( masses = self._resolve_masses(elem_masses) # --- Filter out zero-mass atoms before computing the inertia tensor --- - # This handles unknown elements that were assigned 0.0 + # This handles atoms explicitly assigned a mass of 0.0 via elem_masses mask = masses > 0.0 if not np.any(mask): return None @@ -1424,6 +1441,9 @@ def calc_moments_of_inertia( inertia += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r)) moments_raw, axes = np.linalg.eigh(inertia) # ascending order guaranteed + # Clamp to zero: floating point arithmetic can produce very small + # negative eigenvalues (e.g. -1e-15) for moments that are physically + # zero, due to numerical noise in the inertia tensor. moments_raw = np.maximum(moments_raw, 0.0) return PrincipalMoments( @@ -1446,9 +1466,9 @@ def calc_rotor_type( moments : PrincipalMoments | None, default None Pre-computed principal moments (amu·Å², ascending). When `None` the moments are computed via - :meth:`calc_moments_of_inertia` using *mass_kwargs*. + :meth:`calc_moments_of_inertia()` using *mass_kwargs*. **mass_kwargs - Forwarded to :meth:`calc_moments_of_inertia` when *moments* is + Forwarded to :meth:`calc_moments_of_inertia()` when `moments` is `None` (i.e. `elem_masses`). Returns @@ -1457,15 +1477,14 @@ def calc_rotor_type( `None` if the structure has no real atoms or all masses vanish. """ if moments is None: - result = self.calc_moments_of_inertia(**mass_kwargs) - if result is None: + moments = self.calc_moments_of_inertia(**mass_kwargs) + if moments is None: return None - moments = result return moments.rotor_type() def calc_rotational_constants( self, - elem_masses: dict[str, float] | None = None, + elem_masses: dict[str | Element, float] | None = None, ) -> RotationalConstants | None: """ Compute rotational constants for this structure. @@ -1475,13 +1494,14 @@ def calc_rotational_constants( Mass priority ------------- - elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) + elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) Parameters ---------- - elem_masses : dict[str, float] | None, default None - Per-element mass overrides keyed by element symbol string - (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. + elem_masses : dict[str | Element, float] | None, default None + Per-element mass overrides keyed by element symbol string or + `Element` instance (e.g. `{"C": 13.003}`). If `None`, default + atomic masses are used. Returns ------- diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index cb92a59b..6ffac7b1 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -170,7 +170,7 @@ def moment_to_mhz(inertia: float | None) -> float | None: constant in MHz. Returns `None` when the moment is `None` or effectively zero (degenerate / linear axis). """ - if inertia is None or inertia < 1e-6: + if inertia is None or inertia < 1e-3: return None I_si = inertia * units.AMU_TO_KG * (units.ANGST_TO_M**2) return constants.H_PLANCK / (8.0 * np.pi**2 * I_si) / 1e6 diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_rmsd.py index e266063c..83d54667 100644 --- a/tests/unit/test_rmsd.py +++ b/tests/unit/test_rmsd.py @@ -431,6 +431,7 @@ def test_translated_zero_rmsd(self, water, water_translated): def test_rotated_zero_rmsd(self, water, water_rotated): """Kabsch should align the rotation and give ~0 RMSD.""" centered = water.centered_structure() + print(centered.get_coordinates()) assert pytest.approx(centered.rmsd_kabsch(water_rotated), abs=1e-6) == 0.0 def test_kabsch_le_rmsd(self, water): diff --git a/tests/unit/test_rotconst.py b/tests/unit/test_rotconst.py index 7e48d50b..3d3260e0 100644 --- a/tests/unit/test_rotconst.py +++ b/tests/unit/test_rotconst.py @@ -251,13 +251,12 @@ def test_axes_orthonormal(self, water): np.testing.assert_allclose(product, np.eye(3), atol=1e-10) def test_str_output(self, water): + """str(PrincipalMoments) should match the expected formatted output exactly.""" pm = water.calc_moments_of_inertia() assert pm is not None - s = str(pm) - assert "Moments of inertia (amu·Å²):" in s - assert "Ia" in s - assert "Ib" in s - assert "Ic" in s + assert str(pm) == ( + "Moments of inertia (amu·Å²):\n Ia = 0.641840\n Ib = 1.155114\n Ic = 1.796955" + ) # ============================================================ @@ -454,12 +453,15 @@ def test_get_in_wavenumbers_consistent_with_mhz(self, water): assert pytest.approx(wn[0], rel=1e-9) == mhz_to_wavenumber(rc.A) def test_str_output(self, water): + """str(RotationalConstants) should match the expected formatted output exactly.""" rc = water.calc_rotational_constants() assert rc is not None - s = str(rc) - assert "Rotational constants" in s - assert "MHz" in s - assert "cm⁻¹" in s + assert str(rc) == ( + "Rotational constants:\n" + " A = 787390.377635 MHz (26.264516 cm⁻¹)\n" + " B = 437514.333074 MHz (14.593907 cm⁻¹)\n" + " C = 281241.939008 MHz (9.381221 cm⁻¹)" + ) def test_a_not_in_sync_problem(self, water): """ From 53b9c4bff4f33135b2eb1c2961813b618b7b5ed2 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Mon, 15 Jun 2026 12:20:39 +0200 Subject: [PATCH 19/27] Fixed some typos on structure.py and added some docstrings in test_rmsd.py --- src/opi/input/structures/structure.py | 10 +++++----- tests/unit/test_rmsd.py | 4 ++++ 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 0cafd4bd..595ba68e 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1110,7 +1110,7 @@ def get_coordinates( def set_coordinates(self, coords: npt.NDArray[np.float64]) -> None: """ - Update the coordinates of all atoms in place. + Update the coordinates of all atoms in-place. Parameters ---------- @@ -1353,7 +1353,7 @@ def _resolve_masses( Parameters ---------- - elem_masses : dict[str | Element, float] | None, default None + elem_masses : dict[str | Element, float] | None, default: None Per-element mass overrides keyed by element symbol string or `Element` instance (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. @@ -1406,7 +1406,7 @@ def calc_moments_of_inertia( Parameters ---------- - elem_masses : dict[str | Element, float] | None, default None + elem_masses : dict[str | Element, float] | None, default: None Per-element mass overrides keyed by element symbol string or `Element` instance (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. @@ -1463,7 +1463,7 @@ def calc_rotor_type( Parameters ---------- - moments : PrincipalMoments | None, default None + moments : PrincipalMoments | None, default: None Pre-computed principal moments (amu·Å², ascending). When `None` the moments are computed via :meth:`calc_moments_of_inertia()` using *mass_kwargs*. @@ -1498,7 +1498,7 @@ def calc_rotational_constants( Parameters ---------- - elem_masses : dict[str | Element, float] | None, default None + elem_masses : dict[str | Element, float] | None, default: None Per-element mass overrides keyed by element symbol string or `Element` instance (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_rmsd.py index 83d54667..fab38f70 100644 --- a/tests/unit/test_rmsd.py +++ b/tests/unit/test_rmsd.py @@ -167,14 +167,17 @@ def test_empty_structure(self): class TestGetCoordinates: def test_shape(self, water): + """get_coordinates() should return an array of shape (N, 3).""" coords = water.get_coordinates() assert coords.shape == (3, 3) def test_dtype(self, water): + """get_coordinates() should always return a float64 array.""" coords = water.get_coordinates() assert coords.dtype == np.float64 def test_values(self, water): + """get_coordinates() should return the correct coordinate values.""" coords = water.get_coordinates() np.testing.assert_allclose(coords[0], [0.0, 0.0, 0.119748]) @@ -194,6 +197,7 @@ def test_only_atoms_subset(self, water): assert coords.shape == (2, 3) def test_single_atom(self): + """get_coordinates() on a single-atom structure should return shape (1, 3).""" s = Structure.from_xyz_block("1\n\nC 1.0 2.0 3.0") coords = s.get_coordinates() np.testing.assert_allclose(coords[0], [1.0, 2.0, 3.0]) From 144d9593575cb0259d6648f80521b53d374da655 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Thu, 2 Jul 2026 15:53:25 +0200 Subject: [PATCH 20/27] Fixed typos, added explanatory docstrings and polished the code --- src/opi/input/structures/structure.py | 19 +++++++++++-------- src/opi/utils/rotconst.py | 4 ++++ tests/unit/test_rmsd.py | 1 - 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 595ba68e..97b8d326 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1104,6 +1104,9 @@ def get_coordinates( ------- npt.NDArray[np.float64], shape (N, 3) Coordinates in the same order as the selected atoms. + + The `dtype` is explicitly set to `float64` to guarantee the return type + regardless of how coordinates are stored internally. """ atom_list = [self.atoms[i] for i in only_atoms] if only_atoms else self.atoms return np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) @@ -1347,10 +1350,6 @@ def _resolve_masses( """ Resolve masses for all real `Atom` instances without mutating any atom. - Mass priority - ------------- - elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) - Parameters ---------- elem_masses : dict[str | Element, float] | None, default: None @@ -1358,6 +1357,10 @@ def _resolve_masses( `Element` instance (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. + Mass priority + ------------- + elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) + Returns ------- npt.NDArray[np.float64], shape (N,) @@ -1400,10 +1403,6 @@ def calc_moments_of_inertia( `GhostAtom` is a subclass of `Atom`, so an exact type check (`type(a) is Atom`) is used rather than `isinstance`. - Mass priority - ------------- - elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) - Parameters ---------- elem_masses : dict[str | Element, float] | None, default: None @@ -1411,6 +1410,10 @@ def calc_moments_of_inertia( `Element` instance (e.g. `{"C": 13.003}`). If `None`, default atomic masses are used. + Mass priority + ------------- + elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) + Returns ------- PrincipalMoments | None diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index 6ffac7b1..98857657 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -169,6 +169,10 @@ def moment_to_mhz(inertia: float | None) -> float | None: Convert a principal moment of inertia (amu·Å²) to a rotational constant in MHz. Returns `None` when the moment is `None` or effectively zero (degenerate / linear axis). + + The tolerance value (1e-3) aligns with the default in `rotor_type()`, so a moment + that's considered "zero" by the classifier also returns `None` from the conversion + rather than producing a large MHz value. """ if inertia is None or inertia < 1e-3: return None diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_rmsd.py index fab38f70..cbdef1b2 100644 --- a/tests/unit/test_rmsd.py +++ b/tests/unit/test_rmsd.py @@ -435,7 +435,6 @@ def test_translated_zero_rmsd(self, water, water_translated): def test_rotated_zero_rmsd(self, water, water_rotated): """Kabsch should align the rotation and give ~0 RMSD.""" centered = water.centered_structure() - print(centered.get_coordinates()) assert pytest.approx(centered.rmsd_kabsch(water_rotated), abs=1e-6) == 0.0 def test_kabsch_le_rmsd(self, water): From 18ca2ce84abe38aec152ae8a66c93aed22eb0a8e Mon Sep 17 00:00:00 2001 From: crjacinto Date: Fri, 10 Jul 2026 16:20:59 +0200 Subject: [PATCH 21/27] Modified change log --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f0ec98bd..d58a5f57 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,9 @@ - Add the keyword `dumpactints` to `BlockOutput` (#245). - Added fallback keyword argument for `get_final_energy`, `get_gradient`, and `get_structure` that allows to parse these properties from the `.out` file if no JSON output is available (#237). - Added `get_frequencies`, `get_imaginary_frequencies`, `is_pes_minimum`, and `is_pes_transition_state` to the `Output` class (#247). +- Added `rmsd` and `rmsd_kabsch` to calculate RMSD without and with rotational alignment, respectively (#230). +- Added `calc_rotational_constants` to calculate molecular rotational constants (#230). +- Added `calc_rotor_type` to classify a molecule's rotor type (#230). ### Changed - Refactored methods from Runner into BaseRunner (#193) From cfd0a616d113cf4253b04136e8d39b47c36bd651 Mon Sep 17 00:00:00 2001 From: crjacinto Date: Fri, 10 Jul 2026 17:08:07 +0200 Subject: [PATCH 22/27] Modified mass priority for rotational constants calculation. Corrected speed of light units. Cleaned up unit tests for RMSD and RotConst. Modifed changelog --- src/opi/input/structures/structure.py | 21 +++++++++++---------- src/opi/utils/constants.py | 2 +- src/opi/utils/rotconst.py | 2 +- tests/unit/test_rmsd.py | 9 +++++++++ tests/unit/test_rotconst.py | 25 +++++++++---------------- 5 files changed, 31 insertions(+), 28 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 97b8d326..08e6cfa8 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1359,7 +1359,7 @@ def _resolve_masses( Mass priority ------------- - elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) + atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) Returns ------- @@ -1374,12 +1374,12 @@ def _resolve_masses( masses_list: list[float] = [] for atom in self.real_atoms: - if atom.element in normalised: - # Per-element override takes highest priority - m = normalised[atom.element] - elif atom.mass is not None: - # Mass set directly on the atom takes second priority + if atom.mass is not None: + # Mass set directly on the atom is the most specific → highest priority m = atom.mass + elif atom.element in normalised: + # Per-element override is next + m = normalised[atom.element] elif atom.element in ATOMIC_MASSES_FROM_ELEMENT: # Fall back to default atomic masses m = ATOMIC_MASSES_FROM_ELEMENT[atom.element] @@ -1412,7 +1412,7 @@ def calc_moments_of_inertia( Mass priority ------------- - elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) + atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) Returns ------- @@ -1497,14 +1497,15 @@ def calc_rotational_constants( Mass priority ------------- - elem_masses > atom.mass > default (ATOMIC_MASSES_FROM_ELEMENT) + atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) Parameters ---------- elem_masses : dict[str | Element, float] | None, default: None Per-element mass overrides keyed by element symbol string or - `Element` instance (e.g. `{"C": 13.003}`). If `None`, default - atomic masses are used. + `Element` instance (e.g. `{"C": 13.003}`). Only applied to atoms + that do not have a mass set directly via `atom.mass`. If `None`, + masses fall back to `atom.mass` or the default atomic masses. Returns ------- diff --git a/src/opi/utils/constants.py b/src/opi/utils/constants.py index 033e6bbe..20b0980c 100644 --- a/src/opi/utils/constants.py +++ b/src/opi/utils/constants.py @@ -5,4 +5,4 @@ # Planck's Constant [m^2*kg/s] H_PLANCK = 6.62607015e-34 # Speed of light in vacuum [m/s] -C = 2.99792458e10 +C = 2.99792458e8 diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index 98857657..62b49b0b 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -187,4 +187,4 @@ def mhz_to_wavenumber(mhz: float | None) -> float | None: """ if mhz is None: return None - return mhz * 1e6 / constants.C + return mhz * 1e6 / constants.C / 100.0 diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_rmsd.py index cbdef1b2..c5b8df7e 100644 --- a/tests/unit/test_rmsd.py +++ b/tests/unit/test_rmsd.py @@ -126,6 +126,7 @@ def mixed_structure() -> Structure: # ============================================================ +@pytest.mark.unit class TestRealAtoms: def test_returns_only_real_atoms(self, mixed_structure): """PointCharge should be excluded.""" @@ -165,6 +166,7 @@ def test_empty_structure(self): # ============================================================ +@pytest.mark.unit class TestGetCoordinates: def test_shape(self, water): """get_coordinates() should return an array of shape (N, 3).""" @@ -208,6 +210,7 @@ def test_single_atom(self): # ============================================================ +@pytest.mark.unit class TestGetCoordinatesAtCentroid: def test_centroid_is_zero(self, water): coords = water.get_coordinates_at_centroid() @@ -233,6 +236,7 @@ def test_translation_invariant(self, water, water_translated): # ============================================================ +@pytest.mark.unit class TestSetCoordinates: def test_mutates_in_place(self, water): """set_coordinates should modify the structure in place.""" @@ -262,6 +266,7 @@ def test_deepcopy_is_independent(self, water): # ============================================================ +@pytest.mark.unit class TestCenteredStructure: def test_centroid_at_origin(self, water): centered = water.centered_structure() @@ -298,6 +303,7 @@ def test_translated_centered_equals_original_centered(self, water, water_transla # ============================================================ +@pytest.mark.unit class TestFilteredAtoms: def test_default_returns_only_atom_indices(self, mixed_structure): """PointCharge should be excluded — only index 0 (the Atom) returned.""" @@ -343,6 +349,7 @@ def test_only_atoms_takes_priority_over_ignore_hs(self, ethanol): # ============================================================ +@pytest.mark.unit class TestValidateRmsdCompatibility: def test_compatible_structures(self, water): atoms = water.real_atoms @@ -371,6 +378,7 @@ def test_error_message_uses_natural_counting(self): # ============================================================ +@pytest.mark.unit class TestRmsd: def test_identical_structures_zero_rmsd(self, water): assert pytest.approx(water.rmsd(water), abs=1e-10) == 0.0 @@ -425,6 +433,7 @@ def test_result_is_float(self, water): # ============================================================ +@pytest.mark.unit class TestRmsdKabsch: def test_identical_structures_zero_rmsd(self, water): assert pytest.approx(water.rmsd_kabsch(water), abs=1e-10) == 0.0 diff --git a/tests/unit/test_rotconst.py b/tests/unit/test_rotconst.py index 3d3260e0..efb371cb 100644 --- a/tests/unit/test_rotconst.py +++ b/tests/unit/test_rotconst.py @@ -157,6 +157,7 @@ def ch3cl() -> Structure: # ============================================================ +@pytest.mark.unit class TestHelperFunctions: def test_moment_to_mhz_none_input(self): assert moment_to_mhz(None) is None @@ -198,6 +199,7 @@ def test_mhz_to_wavenumber_proportional(self): # ============================================================ +@pytest.mark.unit class TestPrincipalMoments: def test_rotor_type_monoatomic(self, single_atom): pm = single_atom.calc_moments_of_inertia() @@ -264,6 +266,7 @@ def test_str_output(self, water): # ============================================================ +@pytest.mark.unit class TestCalcMomentsOfInertia: def test_returns_none_for_empty_structure(self): assert make_no_real_atoms_structure().calc_moments_of_inertia() is None @@ -320,26 +323,14 @@ def test_elem_masses_override(self, water): assert pm_default is not None and pm_elem is not None assert pm_default.Ic != pytest.approx(pm_elem.Ic, rel=1e-3) - def test_elem_masses_priority_over_atom_mass(self, water): - """`elem_masses` should take priority over `atom.mass` for the same atom.""" + def test_atom_mass_priority_over_elem_masses(self, water): + """`atom.mass` should take priority over `elem_masses` for the same atom.""" water.real_atoms[0].mass = 17.999 pm_atom_mass = water.calc_moments_of_inertia() pm_elem = water.calc_moments_of_inertia(elem_masses={"O": 16.0}) assert pm_atom_mass is not None and pm_elem is not None - assert pm_atom_mass.Ic != pytest.approx(pm_elem.Ic, rel=1e-3) - - def test_unknown_element_warns_and_excludes(self): - """Atoms with unknown elements should warn and be assigned mass 0.""" - structure = Structure( - atoms=[ - Atom(element=Element("C"), coordinates=Coordinates(coordinates=(0.0, 0.0, 0.0))), - Atom(element=Element("H"), coordinates=Coordinates(coordinates=(1.0, 0.0, 0.0))), - ] - ) - # Force C mass to 0 via atom.mass - structure.real_atoms[0].mass = 0.0 - pm = structure.calc_moments_of_inertia() - assert pm is not None + # atom.mass wins, so specifying elem_masses for O should NOT change the result + assert pytest.approx(pm_atom_mass.Ic, rel=1e-9) == pm_elem.Ic def test_zero_mass_atoms_filtered(self, water): """Atoms with zero mass are excluded; remaining atoms determine the rotor type.""" @@ -370,6 +361,7 @@ def test_water_moments_known_values(self, water): # ============================================================ +@pytest.mark.unit class TestCalcRotorType: def test_with_precomputed_moments(self, water): pm = water.calc_moments_of_inertia() @@ -396,6 +388,7 @@ def test_all_rotor_types(self, single_atom, hcl, methane, benzene, ammonia, ch3c # ============================================================ +@pytest.mark.unit class TestCalcRotationalConstants: def test_returns_none_for_empty_structure(self): assert make_no_real_atoms_structure().calc_rotational_constants() is None From 27755e78ce478a48747a27633b95abc836f4b776 Mon Sep 17 00:00:00 2001 From: haneug <38649381+haneug@users.noreply.github.com> Date: Mon, 13 Jul 2026 11:56:33 +0200 Subject: [PATCH 23/27] chore: update CHANGELOG.md --- CHANGELOG.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 61692e36..e8e4b015 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,9 +20,9 @@ - Added fallback keyword argument for `get_final_energy`, `get_gradient`, and `get_structure` that allows to parse these properties from the `.out` file if no JSON output is available (#237). - Added `get_frequencies`, `get_imaginary_frequencies`, `is_pes_minimum`, and `is_pes_transition_state` to the `Output` class (#247). - `Output.parse()`, `Output.parse_property()`, `Output.parse_gbw()`, and `Output.__init__()` now accept a `strict` parameter (default True). When set to False, output fields that fail Pydantic validation are silently set to None and a UserWarning is emitted instead of raising a ValidationError (#248). -- Added `rmsd` and `rmsd_kabsch` to calculate RMSD without and with rotational alignment, respectively (#230). -- Added `calc_rotational_constants` to calculate molecular rotational constants (#230). -- Added `calc_rotor_type` to classify a molecule's rotor type (#230). +- Added `rmsd` and `rmsd_kabsch` to calculate RMSD without and with rotational alignment to `Structure` class (#230). +- Added `calc_rotational_constants` to calculate molecular rotational constants to `Structure` class (#230). +- Added `calc_rotor_type` to classify a molecule's rotor type to `Structure` class (#230). ### Changed - Refactored methods from Runner into BaseRunner (#193) From 60f91e5fe530ef0ae201453d8fb0727341e1f07f Mon Sep 17 00:00:00 2001 From: haneug <38649381+haneug@users.noreply.github.com> Date: Mon, 13 Jul 2026 11:56:57 +0200 Subject: [PATCH 24/27] docs: small docstring improvement --- src/opi/input/structures/structure.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 08e6cfa8..4d3b52db 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -1505,7 +1505,7 @@ def calc_rotational_constants( Per-element mass overrides keyed by element symbol string or `Element` instance (e.g. `{"C": 13.003}`). Only applied to atoms that do not have a mass set directly via `atom.mass`. If `None`, - masses fall back to `atom.mass` or the default atomic masses. + masses fall back to default atomic masses. Returns ------- From c1428943fe874b58ccdac3feebf2cdb0a88b1749 Mon Sep 17 00:00:00 2001 From: haneug <38649381+haneug@users.noreply.github.com> Date: Mon, 13 Jul 2026 11:59:42 +0200 Subject: [PATCH 25/27] refactor: put cm conversion into numerator --- src/opi/utils/rotconst.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py index 62b49b0b..f2c9a634 100644 --- a/src/opi/utils/rotconst.py +++ b/src/opi/utils/rotconst.py @@ -187,4 +187,4 @@ def mhz_to_wavenumber(mhz: float | None) -> float | None: """ if mhz is None: return None - return mhz * 1e6 / constants.C / 100.0 + return mhz * 1e4 / constants.C From 6c1881211cd41e2ec7b3622c3e899d53adcd694c Mon Sep 17 00:00:00 2001 From: haneug <38649381+haneug@users.noreply.github.com> Date: Mon, 13 Jul 2026 12:01:36 +0200 Subject: [PATCH 26/27] test: rename test files to indicate structure class --- tests/unit/{test_rmsd.py => test_structure_rmsd.py} | 0 tests/unit/{test_rotconst.py => test_structure_rotconst.py} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename tests/unit/{test_rmsd.py => test_structure_rmsd.py} (100%) rename tests/unit/{test_rotconst.py => test_structure_rotconst.py} (100%) diff --git a/tests/unit/test_rmsd.py b/tests/unit/test_structure_rmsd.py similarity index 100% rename from tests/unit/test_rmsd.py rename to tests/unit/test_structure_rmsd.py diff --git a/tests/unit/test_rotconst.py b/tests/unit/test_structure_rotconst.py similarity index 100% rename from tests/unit/test_rotconst.py rename to tests/unit/test_structure_rotconst.py From 22bc8a662f731a2bdb2ccad85219c648bc2ff47e Mon Sep 17 00:00:00 2001 From: haneug <38649381+haneug@users.noreply.github.com> Date: Mon, 13 Jul 2026 12:14:04 +0200 Subject: [PATCH 27/27] refactor: change order of functions in Structure class --- src/opi/input/structures/structure.py | 492 +++++++++++++------------- 1 file changed, 246 insertions(+), 246 deletions(-) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 4d3b52db..6f012b84 100644 --- a/src/opi/input/structures/structure.py +++ b/src/opi/input/structures/structure.py @@ -106,6 +106,15 @@ def atoms( value = [value] self._atoms = list(value) + @property + def real_atoms(self) -> list[Atom]: + """ + Return only real `Atom` instances, excluding `GhostAtom`, `PointCharge`, + and `EmbeddingPotential`. Note that `GhostAtom` is a subclass of `Atom`, + so `type(a) is Atom` is used rather than `isinstance`. + """ + return [a for a in self.atoms if type(a) is Atom] + @property def charge(self) -> int: return self._charge @@ -222,6 +231,23 @@ def combine_molecules(cls, structure1: "Structure", structure2: "Structure") -> new_atoms = structure1.atoms + structure2.atoms return Structure(atoms=new_atoms) + def centered_structure(self) -> "Structure": + """ + Return a new `Structure` centered at its centroid. + Only `Atom` instances contribute to the centroid calculation; + `GhostAtom`, `PointCharge`, and `EmbeddingPotential` are ignored. + + Returns + ------- + Structure + New `Structure` with centered coordinates. + """ + new_structure = copy.deepcopy(self) + real_indices = [i for i, a in enumerate(new_structure.atoms) if type(a) is Atom] + centroid = new_structure.get_coordinates(only_atoms=real_indices).mean(axis=0) + new_structure.set_coordinates(new_structure.get_coordinates() - centroid) + return new_structure + def format_orca(self) -> str: """ Returns string representation of Molecule @@ -343,6 +369,73 @@ def extract_substructure(self, indexes: list[int]) -> "Structure": """ return Structure(atoms=[self.atoms[i] for i in indexes]) + def get_coordinates( + self, + only_atoms: Sequence[int] = (), + ) -> npt.NDArray[np.float64]: + """ + Return the coordinates of atoms as a numpy array. + + Parameters + ---------- + only_atoms : Sequence[int], default: () + Indices of atoms in `self.atoms` to include. If empty, all + atoms are included. + + Returns + ------- + npt.NDArray[np.float64], shape (N, 3) + Coordinates in the same order as the selected atoms. + + The `dtype` is explicitly set to `float64` to guarantee the return type + regardless of how coordinates are stored internally. + """ + atom_list = [self.atoms[i] for i in only_atoms] if only_atoms else self.atoms + return np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) + + def get_coordinates_at_centroid( + self, + only_atoms: Sequence[int] = (), + ) -> npt.NDArray[np.float64]: + """ + Return coordinates translated to the centroid. + + Parameters + ---------- + only_atoms : Sequence[int], default: () + Indices of atoms in `self.atoms` to include. If empty, all + atoms are included. + + Returns + ------- + npt.NDArray[np.float64], shape (N, 3) + Centroid-translated coordinates in the same order as the selected atoms. + """ + coords = self.get_coordinates(only_atoms=only_atoms) + return coords - np.mean(coords, axis=0, dtype=np.float64) + + def set_coordinates(self, coords: npt.NDArray[np.float64]) -> None: + """ + Update the coordinates of all atoms in-place. + + Parameters + ---------- + coords : npt.NDArray[np.float64], shape (N, 3) + New coordinates for all atoms in the same order as `self.atoms`. + + Raises + ------ + ValueError + If *coords* shape does not match the number of atoms. + """ + coords = np.asarray(coords, dtype=np.float64) + if coords.shape != (len(self.atoms), 3): + raise ValueError( + f"coords shape {coords.shape} does not match expected ({len(self.atoms)}, 3)" + ) + for atom, new_coord in zip(self.atoms, coords): + atom.coordinates = new_coord + def update_coordinates(self, array: npt.NDArray[np.float64]) -> None: """ Validates dimensions of array first @@ -1017,208 +1110,10 @@ def from_lists( return cls(atoms=atoms, charge=charge, multiplicity=multiplicity) - @classmethod - def _iter_xyz_structures( - cls, - tracked: TrackingTextIO, - charge: int = 0, - multiplicity: int = 1, - comment_symbols: str | Sequence[str] | None = None, - n_struc_limit: int | None = None, - ) -> Iterator["Structure"]: - """ - Yield properties from the buffer until exhausted or the limit is reached. - - Parameters - ---------- - tracked: TrackingTextIO - A buffer that contains XYZ file data. - charge : int, default: 0 - Optional charge for the structure - multiplicity : int, default: 1 - Optional multiplicity for the structure - comment_symbols: str | Sequence[str] | None, default: None - List of symbols that indicate user comments in the XYZ data. - User comments have to start with the given symbol, fill a whole line, and come before the actual XYZ data. - n_struc_limit: int | None, default: None - If >0, only read the first n structures. - - Returns - -------- - Iterator["Structure"] - Iterator of `Structure` object extracted from the buffer. - - Raises - -------- - ValueError - If `n_struc_limit` is negative or zero. - - """ - - if n_struc_limit is not None and n_struc_limit < 0: - raise ValueError("n_struc_limit must be None, 0, or a positive integer") - - n_struc = 0 - while True: - try: - struct = cls.from_xyz_buffer( - tracked, - charge=charge, - multiplicity=multiplicity, - comment_symbols=comment_symbols, - ) - except EOFError: - break - yield struct - n_struc += 1 - if n_struc_limit and n_struc >= n_struc_limit: - break - # ------------------------------------------------------------------ # # RMSD # # ------------------------------------------------------------------ # - @property - def real_atoms(self) -> list[Atom]: - """ - Return only real `Atom` instances, excluding `GhostAtom`, `PointCharge`, - and `EmbeddingPotential`. Note that `GhostAtom` is a subclass of `Atom`, - so `type(a) is Atom` is used rather than `isinstance`. - """ - return [a for a in self.atoms if type(a) is Atom] - - def get_coordinates( - self, - only_atoms: Sequence[int] = (), - ) -> npt.NDArray[np.float64]: - """ - Return the coordinates of atoms as a numpy array. - - Parameters - ---------- - only_atoms : Sequence[int], default: () - Indices of atoms in `self.atoms` to include. If empty, all - atoms are included. - - Returns - ------- - npt.NDArray[np.float64], shape (N, 3) - Coordinates in the same order as the selected atoms. - - The `dtype` is explicitly set to `float64` to guarantee the return type - regardless of how coordinates are stored internally. - """ - atom_list = [self.atoms[i] for i in only_atoms] if only_atoms else self.atoms - return np.array([a.coordinates.coordinates for a in atom_list], dtype=np.float64) - - def set_coordinates(self, coords: npt.NDArray[np.float64]) -> None: - """ - Update the coordinates of all atoms in-place. - - Parameters - ---------- - coords : npt.NDArray[np.float64], shape (N, 3) - New coordinates for all atoms in the same order as `self.atoms`. - - Raises - ------ - ValueError - If *coords* shape does not match the number of atoms. - """ - coords = np.asarray(coords, dtype=np.float64) - if coords.shape != (len(self.atoms), 3): - raise ValueError( - f"coords shape {coords.shape} does not match expected ({len(self.atoms)}, 3)" - ) - for atom, new_coord in zip(self.atoms, coords): - atom.coordinates = new_coord - - def centered_structure(self) -> "Structure": - """ - Return a new `Structure` centered at its centroid. - Only `Atom` instances contribute to the centroid calculation; - `GhostAtom`, `PointCharge`, and `EmbeddingPotential` are ignored. - - Returns - ------- - Structure - New `Structure` with centered coordinates. - """ - new_structure = copy.deepcopy(self) - real_indices = [i for i, a in enumerate(new_structure.atoms) if type(a) is Atom] - centroid = new_structure.get_coordinates(only_atoms=real_indices).mean(axis=0) - new_structure.set_coordinates(new_structure.get_coordinates() - centroid) - return new_structure - - def _filtered_atoms( - self, - only_atoms: Sequence[int], - ignore_hs: bool, - ) -> list[int]: - """ - Return indices of real `Atom` instances after applying `only_atoms` - and `ignore_hs` filters. - - The final `type(a) is Atom` check ensures that explicitly indexed atoms - (via `only_atoms`) cannot introduce `GhostAtom`, `PointCharge`, or - `EmbeddingPotential` instances into the returned list. Note that `GhostAtom` - is a subclass of `Atom`, so `isinstance` is intentionally avoided here. - """ - if only_atoms: - candidates = list(only_atoms) - elif ignore_hs: - candidates = [ - i for i, a in enumerate(self.atoms) if type(a) is Atom and a.element != Element.H - ] - else: - candidates = [i for i, a in enumerate(self.atoms) if type(a) is Atom] - return [i for i in candidates if type(self.atoms[i]) is Atom] - - @staticmethod - def _validate_rmsd_compatibility(atoms1: list[Atom], atoms2: list[Atom]) -> None: - """ - Raise `ValueError` if `atoms1` and `atoms2` differ in count or element order. - """ - if len(atoms1) != len(atoms2): - raise ValueError( - f"Structures have different number of atoms: {len(atoms1)} vs {len(atoms2)}" - ) - for i, (a, b) in enumerate(zip(atoms1, atoms2), start=1): - if a.element != b.element: - raise ValueError(f"Atom mismatch at position {i}: {a.element!r} != {b.element!r}") - - @staticmethod - def _rmsd_coords( - coords1: npt.NDArray[np.float64], - coords2: npt.NDArray[np.float64], - ) -> float: - """ - Compute RMSD between two aligned (N, 3) coordinate arrays. - """ - diff = coords1 - coords2 - return float(np.sqrt(np.sum(diff**2) / len(coords1))) - - def get_coordinates_at_centroid( - self, - only_atoms: Sequence[int] = (), - ) -> npt.NDArray[np.float64]: - """ - Return coordinates translated to the centroid. - - Parameters - ---------- - only_atoms : Sequence[int], default: () - Indices of atoms in `self.atoms` to include. If empty, all - atoms are included. - - Returns - ------- - npt.NDArray[np.float64], shape (N, 3) - Centroid-translated coordinates in the same order as the selected atoms. - """ - coords = self.get_coordinates(only_atoms=only_atoms) - return coords - np.mean(coords, axis=0, dtype=np.float64) - def rmsd( self, other: "Structure", @@ -1343,54 +1238,6 @@ def rmsd_kabsch( # Moment of inertia # # ------------------------------------------------------------------ # - def _resolve_masses( - self, - elem_masses: dict[str | Element, float] | None = None, - ) -> npt.NDArray[np.float64]: - """ - Resolve masses for all real `Atom` instances without mutating any atom. - - Parameters - ---------- - elem_masses : dict[str | Element, float] | None, default: None - Per-element mass overrides keyed by element symbol string or - `Element` instance (e.g. `{"C": 13.003}`). If `None`, default - atomic masses are used. - - Mass priority - ------------- - atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) - - Returns - ------- - npt.NDArray[np.float64], shape (N,) - Masses in amu in the same order as `self.real_atoms`. - """ - elem_masses = elem_masses or {} - # Normalise keys to Element for consistent lookup - normalised: dict[Element, float] = {} - for k, v in elem_masses.items(): - normalised[Element(k) if isinstance(k, str) else k] = v - - masses_list: list[float] = [] - for atom in self.real_atoms: - if atom.mass is not None: - # Mass set directly on the atom is the most specific → highest priority - m = atom.mass - elif atom.element in normalised: - # Per-element override is next - m = normalised[atom.element] - elif atom.element in ATOMIC_MASSES_FROM_ELEMENT: - # Fall back to default atomic masses - m = ATOMIC_MASSES_FROM_ELEMENT[atom.element] - else: - raise ValueError( - f"Unknown element '{atom.element.value}': no mass available. " - "Provide a mass via `elem_masses` or set `atom.mass` directly." - ) - masses_list.append(m) - return np.array(masses_list, dtype=np.float64) - def calc_moments_of_inertia( self, elem_masses: dict[str | Element, float] | None = None, @@ -1519,3 +1366,156 @@ def calc_rotational_constants( B = moment_to_mhz(pm.Ib) C = moment_to_mhz(pm.Ic) return RotationalConstants(A=A, B=B, C=C) + + @classmethod + def _iter_xyz_structures( + cls, + tracked: TrackingTextIO, + charge: int = 0, + multiplicity: int = 1, + comment_symbols: str | Sequence[str] | None = None, + n_struc_limit: int | None = None, + ) -> Iterator["Structure"]: + """ + Yield properties from the buffer until exhausted or the limit is reached. + + Parameters + ---------- + tracked: TrackingTextIO + A buffer that contains XYZ file data. + charge : int, default: 0 + Optional charge for the structure + multiplicity : int, default: 1 + Optional multiplicity for the structure + comment_symbols: str | Sequence[str] | None, default: None + List of symbols that indicate user comments in the XYZ data. + User comments have to start with the given symbol, fill a whole line, and come before the actual XYZ data. + n_struc_limit: int | None, default: None + If >0, only read the first n structures. + + Returns + -------- + Iterator["Structure"] + Iterator of `Structure` object extracted from the buffer. + + Raises + -------- + ValueError + If `n_struc_limit` is negative or zero. + + """ + + if n_struc_limit is not None and n_struc_limit < 0: + raise ValueError("n_struc_limit must be None, 0, or a positive integer") + + n_struc = 0 + while True: + try: + struct = cls.from_xyz_buffer( + tracked, + charge=charge, + multiplicity=multiplicity, + comment_symbols=comment_symbols, + ) + except EOFError: + break + yield struct + n_struc += 1 + if n_struc_limit and n_struc >= n_struc_limit: + break + + def _filtered_atoms( + self, + only_atoms: Sequence[int], + ignore_hs: bool, + ) -> list[int]: + """ + Return indices of real `Atom` instances after applying `only_atoms` + and `ignore_hs` filters. + + The final `type(a) is Atom` check ensures that explicitly indexed atoms + (via `only_atoms`) cannot introduce `GhostAtom`, `PointCharge`, or + `EmbeddingPotential` instances into the returned list. Note that `GhostAtom` + is a subclass of `Atom`, so `isinstance` is intentionally avoided here. + """ + if only_atoms: + candidates = list(only_atoms) + elif ignore_hs: + candidates = [ + i for i, a in enumerate(self.atoms) if type(a) is Atom and a.element != Element.H + ] + else: + candidates = [i for i, a in enumerate(self.atoms) if type(a) is Atom] + return [i for i in candidates if type(self.atoms[i]) is Atom] + + @staticmethod + def _validate_rmsd_compatibility(atoms1: list[Atom], atoms2: list[Atom]) -> None: + """ + Raise `ValueError` if `atoms1` and `atoms2` differ in count or element order. + """ + if len(atoms1) != len(atoms2): + raise ValueError( + f"Structures have different number of atoms: {len(atoms1)} vs {len(atoms2)}" + ) + for i, (a, b) in enumerate(zip(atoms1, atoms2), start=1): + if a.element != b.element: + raise ValueError(f"Atom mismatch at position {i}: {a.element!r} != {b.element!r}") + + @staticmethod + def _rmsd_coords( + coords1: npt.NDArray[np.float64], + coords2: npt.NDArray[np.float64], + ) -> float: + """ + Compute RMSD between two aligned (N, 3) coordinate arrays. + """ + diff = coords1 - coords2 + return float(np.sqrt(np.sum(diff**2) / len(coords1))) + + def _resolve_masses( + self, + elem_masses: dict[str | Element, float] | None = None, + ) -> npt.NDArray[np.float64]: + """ + Resolve masses for all real `Atom` instances without mutating any atom. + + Parameters + ---------- + elem_masses : dict[str | Element, float] | None, default: None + Per-element mass overrides keyed by element symbol string or + `Element` instance (e.g. `{"C": 13.003}`). If `None`, default + atomic masses are used. + + Mass priority + ------------- + atom.mass > elem_masses > default (ATOMIC_MASSES_FROM_ELEMENT) + + Returns + ------- + npt.NDArray[np.float64], shape (N,) + Masses in amu in the same order as `self.real_atoms`. + """ + elem_masses = elem_masses or {} + # Normalise keys to Element for consistent lookup + normalised: dict[Element, float] = {} + for k, v in elem_masses.items(): + normalised[Element(k) if isinstance(k, str) else k] = v + + masses_list: list[float] = [] + for atom in self.real_atoms: + if atom.mass is not None: + # Mass set directly on the atom is the most specific → highest priority + m = atom.mass + elif atom.element in normalised: + # Per-element override is next + m = normalised[atom.element] + elif atom.element in ATOMIC_MASSES_FROM_ELEMENT: + # Fall back to default atomic masses + m = ATOMIC_MASSES_FROM_ELEMENT[atom.element] + else: + raise ValueError( + f"Unknown element '{atom.element.value}': no mass available. " + "Provide a mass via `elem_masses` or set `atom.mass` directly." + ) + masses_list.append(m) + return np.array(masses_list, dtype=np.float64)