Resolve inertia/COM issues (#35, #36, #37) (#73)

Bundles three small, related inertia/centre-of-mass issues (separate PRs
would conflict on `thermo.py`/`test_thermo.py`).

- **#37** — rename the two non-snake_case inertia attributes
`_eigenvalues_I` → `_inertia_eigenvalues` and `_eigenvalues_I_SI` →
`_inertia_eigenvalues_si` (clears pylint C0103; pure rename, no
behaviour change; the one test reference updated).
- **#36** — validate the principal moments of inertia against ASE's
`get_moments_of_inertia` for a linear (CO₂) and nonlinear (H₂O) molecule
(masses pinned to the tool's), plus a translation-invariance test
proving the COM relocation makes the moments position-independent.
- **#35** — document that `System.coord()` returns the *original* input
coordinates and that the COM relocation is an internal detail of the
thermo pipeline (`Thermo._reloc_coord`).

Closes #35, closes #36, closes #37.

Full suite: 255 passed / 11 skipped.

Author: galjos

ThermoScreening/thermo/system.py

@@ -739,12 +739,18 @@ def atomic_masses(self) -> List[float]:
 
     def coord(self) -> np.ndarray:
         """
-        The atom positions of the system.
+        The atom positions of the system, as originally provided.
+
+        These are the input coordinates; they are not centred on the centre of
+        mass. The centre-of-mass relocation used for the moment-of-inertia /
+        rotational analysis is an internal detail of the thermo pipeline
+        (``Thermo._relocate_to_cm`` stores it in ``Thermo._reloc_coord``) and
+        does not modify the coordinates returned here.
 
         Returns
         -------
         np.ndarray
-            The atom positions of the system.
+            The (original, non-relocated) atom positions of the system.
         """
         return np.array([atom.position for atom in self._atoms])
 

ThermoScreening/thermo/thermo.py

@@ -227,7 +227,7 @@ def _compute_rotational_partition_function(self):
 
         if self._n_rot == 0:
             # monatomic: no rotational degrees of freedom
-            self._eigenvalues_I_SI = np.array([])
+            self._inertia_eigenvalues_si = np.array([])
             self._rotational_temperature = np.array([])
             self._rotational_constant = np.array([])
             self._rotational_temperature_xyz = np.nan
@@ -238,17 +238,17 @@ def _compute_rotational_partition_function(self):
         # rotor, the two equal perpendicular moments for a linear molecule
         # (eigenvalues are sorted ascending, so the smallest, ~zero for a linear
         # molecule, is dropped)
-        self._eigenvalues_I_SI = (
-            self._eigenvalues_I[-self._n_rot:]
+        self._inertia_eigenvalues_si = (
+            self._inertia_eigenvalues[-self._n_rot:]
             * PhysicalConstants["u"]
             * PhysicalConstants["A"] ** 2
         )
         self._rotational_temperature = (
             PhysicalConstants["h"] ** 2 / (8 * np.pi**2 * PhysicalConstants["kB"])
-        ) / self._eigenvalues_I_SI
+        ) / self._inertia_eigenvalues_si
 
         self._rotational_constant = (
-            (((PhysicalConstants["hbar"] ** 2) / 2) / self._eigenvalues_I_SI)
+            (((PhysicalConstants["hbar"] ** 2) / 2) / self._inertia_eigenvalues_si)
             / PhysicalConstants["h"]
             * PhysicalConstants["HztoGHz"]
         )
@@ -322,7 +322,7 @@ def _rotational_contribution(self):
         self._compute_inertia_tensor()
         # the inertia tensor is symmetric by construction; eigvalsh returns real,
         # ascending eigenvalues (eig may emit spurious imaginary parts)
-        self._eigenvalues_I = np.linalg.eigvalsh(self._inertia_tensor)
+        self._inertia_eigenvalues = np.linalg.eigvalsh(self._inertia_tensor)
         self._n_rot = self._rotational_dof()
         self._compute_rotational_partition_function()
         self._compute_rotational_entropy()

tests/thermo/test_thermo.py

@@ -180,6 +180,47 @@ def test_total_entropy_matches_ase(symbols, positions, freqs, dof, geometry, sig
     assert ts_total == pytest.approx(ase_total, abs=0.05)
 
 
+@pytest.mark.parametrize(
+    "symbols,positions,freqs,dof",
+    [
+        # nonlinear (H2O) and linear (CO2) -> one ~zero moment for the linear case
+        (["O", "H", "H"], [[0, 0, 0.12], [0, 0.76, -0.48], [0, -0.76, -0.48]],
+         [1595.0, 3657.0, 3756.0], 3),
+        (["C", "O", "O"], [[0, 0, 0], [0, 0, 1.16], [0, 0, -1.16]],
+         [667.0, 667.0, 1333.0, 2349.0], 4),
+    ],
+)
+def test_inertia_eigenvalues_match_ase(symbols, positions, freqs, dof):
+    # the principal moments of inertia (amu*A^2, before the SI conversion) must
+    # match ASE's get_moments_of_inertia for the same geometry and masses
+    thermo = _ts_thermo(symbols, positions, freqs, dof)
+
+    ase_atoms = Atoms(symbols=symbols, positions=positions)
+    ase_atoms.set_masses([Atom(symbol=s, position=np.zeros(3)).mass for s in symbols])
+    ase_moments = np.sort(ase_atoms.get_moments_of_inertia())
+
+    np.testing.assert_allclose(
+        np.sort(thermo._inertia_eigenvalues), ase_moments, atol=1e-6
+    )
+
+
+def test_inertia_eigenvalues_are_translation_invariant():
+    # relocating to the centre of mass must make the moments independent of where
+    # the molecule sits in space
+    symbols = ["O", "H", "H"]
+    positions = np.array([[0, 0, 0.12], [0, 0.76, -0.48], [0, -0.76, -0.48]])
+    freqs, dof = [1595.0, 3657.0, 3756.0], 3
+
+    centred = _ts_thermo(symbols, positions, freqs, dof)
+    shifted = _ts_thermo(symbols, positions + np.array([10.0, -5.0, 3.0]), freqs, dof)
+
+    np.testing.assert_allclose(
+        np.sort(centred._inertia_eigenvalues),
+        np.sort(shifted._inertia_eigenvalues),
+        atol=1e-8,
+    )
+
+
 from ase.build import molecule  # noqa: E402
 
 
@@ -465,7 +506,7 @@ def test_thermo_aq_2(self):
         np.testing.assert_allclose(thermo._inertia_tensor, np.array([[477.579476, -0.002293, 0.023808], [-0.002293, 1612.635725, 0.000317],[0.023808, 0.000317, 1135.056249]]), atol=1e-4)
 
         # somewhere in the array the values should be the same
-        np.testing.assert_allclose(np.sort(thermo._eigenvalues_I), np.sort(
+        np.testing.assert_allclose(np.sort(thermo._inertia_eigenvalues), np.sort(
             np.array([1612.635724886585, 477.579474789699, 1135.056250097402])), atol=1e-10)
 
         np.testing.assert_allclose(np.sort(thermo._rotational_temperature), np.sort(