diff --git a/ThermoScreening/thermo/system.py b/ThermoScreening/thermo/system.py index c395e6e..98ccc84 100644 --- a/ThermoScreening/thermo/system.py +++ b/ThermoScreening/thermo/system.py @@ -1,7 +1,9 @@ import logging +import warnings import numpy as np from beartype.typing import List +from numpy.exceptions import ComplexWarning from pymatgen.core import Molecule, Structure from pymatgen.symmetry.analyzer import SpacegroupAnalyzer @@ -15,6 +17,48 @@ from .cell import Cell +def _real_position(position: np.ndarray) -> np.ndarray: + """ + Return a real floating-point atom position for symmetry analysis. + """ + + real_position = np.real_if_close(np.asarray(position)) + if np.iscomplexobj(real_position): + if not np.allclose(np.imag(real_position), 0.0): + raise TSValueError("Atom positions must be real-valued.") + real_position = np.real(real_position) + + return np.asarray(real_position, dtype=float) + + +def _molecule_from_atoms(atoms: List[Atom]) -> Molecule: + """ + Build a pymatgen molecule with real-valued coordinates. + """ + + names = [] + coordinates = [] + for atom in atoms: + names.append(atom.symbol) + coordinates.append(_real_position(atom.position)) + + return Molecule(names, coordinates) + + +def _point_group_analyzer(molecule): + """ + Build a pymatgen point-group analyzer without surfacing benign cast warnings. + """ + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + category=ComplexWarning, + module=r"pymatgen\.core\.operations", + ) + return PointGroupAnalyzer(molecule) + + def linearity(atoms: List[Atom]) -> bool: """ Checks if the system is linear or non-linear. @@ -197,14 +241,8 @@ def rotational_symmetry_number(atoms: List[Atom]) -> int: int The symmetry number of the system. """ - name = [] - coord = [] - for atom in atoms: - name.append(atom.symbol) - coord.append(atom.position) - - mol = Molecule(name, coord) - symmetry_number = PointGroupAnalyzer(mol).get_rotational_symmetry_number + mol = _molecule_from_atoms(atoms) + symmetry_number = _point_group_analyzer(mol).get_rotational_symmetry_number if callable(symmetry_number): return symmetry_number() return symmetry_number @@ -285,14 +323,8 @@ def rotational_group_calc(atoms: List[Atom]) -> str: str The rotational group of the system. """ - name = [] - coord = [] - for atom in atoms: - name.append(atom.symbol) - coord.append(atom.position) - - mol = Molecule(name, coord) - symb = PointGroupAnalyzer(mol).sch_symbol + mol = _molecule_from_atoms(atoms) + symb = _point_group_analyzer(mol).sch_symbol return symb diff --git a/tests/thermo/test_system.py b/tests/thermo/test_system.py index d45a56f..49912f0 100644 --- a/tests/thermo/test_system.py +++ b/tests/thermo/test_system.py @@ -1,4 +1,5 @@ import unittest +import warnings import numpy as np from ase.atoms import Atoms import ThermoScreening.thermo.system as system_module @@ -44,6 +45,69 @@ def __init__(self, molecule): assert rotational_symmetry_number(atoms) == 7 +def test_symmetry_analysis_converts_zero_imaginary_positions(monkeypatch): + captured_coordinates = [] + + class FakePointGroupAnalyzer: + sch_symbol = "D*h" + + def __init__(self, molecule): + self.molecule = molecule + self.get_rotational_symmetry_number = 2 + captured_coordinates.append(np.array(molecule.cart_coords)) + + atoms = [ + Atom(symbol="H", position=np.array([0.0 + 0.0j, 0.0, 0.0])), + Atom(symbol="H", position=np.array([1.0 + 0.0j, 0.0, 0.0])), + ] + monkeypatch.setattr(system_module, "PointGroupAnalyzer", FakePointGroupAnalyzer) + + assert rotational_symmetry_number(atoms) == 2 + assert system_module.rotational_group_calc(atoms) == "D*h" + assert not np.iscomplexobj(captured_coordinates[0]) + assert not np.iscomplexobj(captured_coordinates[1]) + + +def test_real_position_converts_near_zero_imaginary_noise(): + position = system_module._real_position(np.array([1.0 + 1e-12j, 0.0, 0.0])) + + np.testing.assert_array_equal(position, np.array([1.0, 0.0, 0.0])) + assert not np.iscomplexobj(position) + + +def test_point_group_analyzer_suppresses_pymatgen_complex_warning(monkeypatch): + class FakePointGroupAnalyzer: + sch_symbol = "C1" + + def __init__(self, molecule): + self.molecule = molecule + warnings.warn_explicit( + "Casting complex values to real discards the imaginary part", + system_module.ComplexWarning, + filename="operations.py", + lineno=1, + module="pymatgen.core.operations", + ) + + monkeypatch.setattr(system_module, "PointGroupAnalyzer", FakePointGroupAnalyzer) + + with warnings.catch_warnings(): + warnings.simplefilter("error", system_module.ComplexWarning) + analyzer = system_module._point_group_analyzer(object()) + + assert analyzer.sch_symbol == "C1" + + +def test_symmetry_analysis_rejects_imaginary_positions(): + atoms = [ + Atom(symbol="H", position=np.array([0.0 + 1.0j, 0.0, 0.0])), + Atom(symbol="H", position=np.array([1.0, 0.0, 0.0])), + ] + + with pytest.raises(TSValueError, match="Atom positions must be real-valued"): + rotational_symmetry_number(atoms) + + class TestSystem(unittest.TestCase): atoms = [ Atom(symbol='H', position=np.array([0, 0, 0])),