Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 48 additions & 16 deletions ThermoScreening/thermo/system.py
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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


Expand Down
64 changes: 64 additions & 0 deletions tests/thermo/test_system.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
import unittest
import warnings
import numpy as np
from ase.atoms import Atoms
import ThermoScreening.thermo.system as system_module
Expand Down Expand Up @@ -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])),
Expand Down
Loading