diff --git a/CHANGELOG.md b/CHANGELOG.md index e5e21f08..e8e4b015 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +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 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) diff --git a/src/opi/input/structures/structure.py b/src/opi/input/structures/structure.py index 8dccf329..6f012b84 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 @@ -19,7 +20,13 @@ PointCharge, ) from opi.input.structures.coordinates import Coordinates -from opi.utils.element import Element +from opi.utils.element import ATOMIC_MASSES_FROM_ELEMENT, Element +from opi.utils.rotconst import ( + PrincipalMoments, + RotationalConstants, + RotorType, + moment_to_mhz, +) from opi.utils.tracking_text_io import TrackingTextIO __all__ = ("Structure",) @@ -99,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 @@ -215,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 @@ -336,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 @@ -1010,6 +1110,263 @@ def from_lists( return cls(atoms=atoms, charge=charge, multiplicity=multiplicity) + # ------------------------------------------------------------------ # + # RMSD # + # ------------------------------------------------------------------ # + + 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. + Neither `self` nor *other* is modified; all coordinate operations are + performed on temporary arrays. + + 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. + """ + 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) + + 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. + Neither `self` nor *other* is modified; all coordinate operations are + performed on temporary arrays. + + 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. + """ + 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(only_atoms=indices1) + coords2 = other.get_coordinates_at_centroid(only_atoms=indices2) + + # Kabsch algorithm: find optimal rotation matrix via SVD + # + + # ------------------------------------------------------------------ + # Build covariance matrix + # ------------------------------------------------------------------ + # H = A^T B + # This matrix captures how coordinates from B (coords2) map onto A (coords1) + H = coords1.T @ coords2 + + # ------------------------------------------------------------------ + # 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 + + return self._rmsd_coords(coords1, coords2 @ R) + + # ------------------------------------------------------------------ # + # Moment of inertia # + # ------------------------------------------------------------------ # + + def calc_moments_of_inertia( + self, + elem_masses: dict[str | Element, 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. Note that + `GhostAtom` is a subclass of `Atom`, so an exact type check + (`type(a) is Atom`) is used rather than `isinstance`. + + 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 + ------- + 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. + """ + 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) + + masses = self._resolve_masses(elem_masses) + + # --- Filter out zero-mass atoms before computing the inertia tensor --- + # This handles atoms explicitly assigned a mass of 0.0 via elem_masses + mask = masses > 0.0 + if not np.any(mask): + return None + masses = masses[mask] + coords = coords[mask] + + total_mass = float(masses.sum()) + com = (masses[:, None] * coords).sum(axis=0) / total_mass + coords -= com + + 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)) + + 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( + Ia=float(moments_raw[0]), + Ib=float(moments_raw[1]), + Ic=float(moments_raw[2]), + axes=axes, + ) + + def calc_rotor_type( + self, + moments: PrincipalMoments | None = None, + **mass_kwargs: Any, + ) -> RotorType | None: + """ + Classify the molecular rotor type. + + Parameters + ---------- + 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*. + **mass_kwargs + Forwarded to :meth:`calc_moments_of_inertia()` when `moments` is + `None` (i.e. `elem_masses`). + + Returns + ------- + RotorType | None + `None` if the structure has no real atoms or all masses vanish. + """ + if moments is None: + moments = self.calc_moments_of_inertia(**mass_kwargs) + if moments is None: + return None + return moments.rotor_type() + + def calc_rotational_constants( + self, + elem_masses: dict[str | Element, 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 + ------------- + 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}`). Only applied to atoms + that do not have a mass set directly via `atom.mass`. If `None`, + masses fall back to default atomic masses. + + Returns + ------- + RotationalConstants | None + `None` if no `Atom` instances are present or all masses are zero. + """ + pm = self.calc_moments_of_inertia(elem_masses=elem_masses) + 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) + @classmethod def _iter_xyz_structures( cls, @@ -1066,3 +1423,99 @@ def _iter_xyz_structures( 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) diff --git a/src/opi/utils/constants.py b/src/opi/utils/constants.py new file mode 100644 index 00000000..20b0980c --- /dev/null +++ b/src/opi/utils/constants.py @@ -0,0 +1,8 @@ +""" +Contains universal constants. +""" + +# Planck's Constant [m^2*kg/s] +H_PLANCK = 6.62607015e-34 +# Speed of light in vacuum [m/s] +C = 2.99792458e8 diff --git a/src/opi/utils/element.py b/src/opi/utils/element.py index a8d00c6d..5983c513 100644 --- a/src/opi/utils/element.py +++ b/src/opi/utils/element.py @@ -900,3 +900,246 @@ 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, +} diff --git a/src/opi/utils/rotconst.py b/src/opi/utils/rotconst.py new file mode 100644 index 00000000..f2c9a634 --- /dev/null +++ b/src/opi/utils/rotconst.py @@ -0,0 +1,190 @@ +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", + "moment_to_mhz", + "mhz_to_wavenumber", +) + + +# ============================================================ +# 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" + + 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 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)) + + 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 + + # --- 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(kappa - 1.0) < kappa_tol: + return RotorType.OBLATE_TOP + return RotorType.ASYMMETRIC_TOP + + 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 +# ============================================================ + + +class RotationalConstants: + """ + Stores rotational constants in MHz. + + Attributes + ---------- + A, B, C : float | None + Rotational constants in MHz (`None` for a degenerate axis). + """ + + def __init__( + self, + A: float | None, + B: float | None, + C: float | None, + ) -> None: + self.A = A + self.B = B + self.C = C + + 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(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⁻¹')})" + ) + + +# ============================================================ +# Helper functions (used by Structure methods) +# ============================================================ + + +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 + 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⁻¹. + Returns `None` if `mhz` is `None`. + """ + if mhz is None: + return None + return mhz * 1e4 / constants.C 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 diff --git a/tests/unit/test_structure_rmsd.py b/tests/unit/test_structure_rmsd.py new file mode 100644 index 00000000..c5b8df7e --- /dev/null +++ b/tests/unit/test_structure_rmsd.py @@ -0,0 +1,489 @@ +""" +Unit tests for RMSD and coordinate utilities. + +Covers: +- real_atoms property +- 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 copy + +import numpy as np +import pytest + +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 + +# ============================================================ +# 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 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 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" + ) + + +@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], + ] + ) + rotated = copy.deepcopy(centered) + rotated.set_coordinates(coords @ R.T) + return rotated + + +@pytest.fixture() +def ethanol() -> Structure: + """ + Ethanol (C2H5OH) for ignore_hs tests. + Coordinates are approximate — correctness of geometry is not critical here. + """ + 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" + ) + + +@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), + ] + ) + + +# ============================================================ +# real_atoms +# ============================================================ + + +@pytest.mark.unit +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))), + 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 structure.real_atoms == [] + + +# ============================================================ +# get_coordinates +# ============================================================ + + +@pytest.mark.unit +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]) + + 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): + """Default call includes all atom types.""" + 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): + """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]) + + +# ============================================================ +# get_coordinates_at_centroid +# ============================================================ + + +@pytest.mark.unit +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_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) + + 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 +# ============================================================ + + +@pytest.mark.unit +class TestSetCoordinates: + def test_mutates_in_place(self, water): + """set_coordinates should modify the structure in place.""" + new_coords = np.zeros((3, 3)) + 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_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) + + +# ============================================================ +# centered_structure +# ============================================================ + + +@pytest.mark.unit +class TestCenteredStructure: + def test_centroid_at_origin(self, water): + centered = water.centered_structure() + 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): + 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_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): + c1 = water.centered_structure().get_coordinates() + c2 = water_translated.centered_structure().get_coordinates() + np.testing.assert_allclose(c1, c2, atol=1e-6) + + +# ============================================================ +# _filtered_atoms +# ============================================================ + + +@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.""" + 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): + 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): + 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 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.""" + indices_with = ethanol._filtered_atoms([3, 4], True) + indices_without = ethanol._filtered_atoms([3, 4], False) + assert indices_with == indices_without + + +# ============================================================ +# _validate_rmsd_compatibility +# ============================================================ + + +@pytest.mark.unit +class TestValidateRmsdCompatibility: + def test_compatible_structures(self, water): + atoms = water.real_atoms + Structure._validate_rmsd_compatibility(atoms, atoms) # should not raise + + def test_raises_on_different_count(self, water, ethanol): + with pytest.raises(ValueError, match="different number of atoms"): + 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)] + 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 +# ============================================================ + + +@pytest.mark.unit +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 one heavy atom — RMSD with and without H should differ. + """ + 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.""" + 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 = 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 + + def test_result_is_float(self, water): + assert isinstance(water.rmsd(water), float) + + +# ============================================================ +# rmsd_kabsch +# ============================================================ + + +@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 + + 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 = 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): + 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): + 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): + 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) + + 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_structure_rotconst.py b/tests/unit/test_structure_rotconst.py new file mode 100644 index 00000000..efb371cb --- /dev/null +++ b/tests/unit/test_structure_rotconst.py @@ -0,0 +1,476 @@ +""" +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 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, mass: float | None = None) -> Atom: + return Atom( + element=Element(element), + coordinates=Coordinates(coordinates=(x, y, z)), + mass=mass, + ) + + +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)] + ) + + +# ============================================================ +# Fixtures +# ============================================================ + + +@pytest.fixture() +def water() -> Structure: + """H2O — asymmetric top.""" + 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 co2() -> Structure: + """ + CO2 — linear molecule aligned exactly on Z. + After COM shift, both O atoms lie on the Z axis so Ia = 0 exactly. + """ + 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 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 + 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. Avogadro optimized geometry.""" + return Structure.from_xyz_block( + "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" + "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 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 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 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" + ) + + +# ============================================================ +# moment_to_mhz and mhz_to_wavenumber +# ============================================================ + + +@pytest.mark.unit +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 +# ============================================================ + + +@pytest.mark.unit +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): + """str(PrincipalMoments) should match the expected formatted output exactly.""" + pm = water.calc_moments_of_inertia() + assert pm is not None + assert str(pm) == ( + "Moments of inertia (amu·Å²):\n Ia = 0.641840\n Ib = 1.155114\n Ic = 1.796955" + ) + + +# ============================================================ +# calc_moments_of_inertia +# ============================================================ + + +@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 + + def test_returns_none_for_all_zero_masses(self, water): + """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.""" + 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=[ + 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_atom_mass_override(self, water): + """Setting `atom.mass` directly should change the moments.""" + pm_default = water.calc_moments_of_inertia() + 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.Ic != pytest.approx(pm_custom.Ic, rel=1e-3) + + 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_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 + 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 + # 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.""" + # 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 + + def test_water_moments_known_values(self, water): + """ + Check moments are self-consistent: Ia <= Ib <= Ic and all positive. + 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 + assert 0.1 < pm.Ia < 5.0 + assert 0.1 < pm.Ib < 5.0 + assert 0.1 < pm.Ic < 5.0 + + +# ============================================================ +# calc_rotor_type +# ============================================================ + + +@pytest.mark.unit +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 +# ============================================================ + + +@pytest.mark.unit +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): + """str(RotationalConstants) should match the expected formatted output exactly.""" + rc = water.calc_rotational_constants() + assert rc is not None + 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): + """ + 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 + + 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)