Fix rotational thermochemistry for linear and monatomic species

ThermoScreening/thermo/system.py

@@ -222,7 +222,8 @@ def dof(atoms: List[Atom]) -> int:
     """
     number_of_atoms = len(atoms)
     if number_of_atoms == 1:
-        return 3
+        # a monatomic species has no rotational or vibrational degrees of freedom
+        return 0
     if number_of_atoms > 1:
         return (
             (3 * number_of_atoms - 5) if linearity(atoms) else (3 * number_of_atoms - 6)

ThermoScreening/thermo/thermo.py

@@ -9,7 +9,7 @@
 from ThermoScreening.utils.physicalConstants import PhysicalConstants
 from ThermoScreening import __package_name__
 
-from .system import System
+from .system import System, linearity
 
 
 def _real_scalar(value) -> float:
@@ -209,8 +209,23 @@ def _compute_rotational_partition_function(self):
         None
         """
 
+        if self._n_rot == 0:
+            # monatomic: no rotational degrees of freedom
+            self._eigenvalues_I_SI = np.array([])
+            self._rotational_temperature = np.array([])
+            self._rotational_constant = np.array([])
+            self._rotational_temperature_xyz = np.nan
+            self._rotational_partition_function = 1.0
+            return
+
+        # the non-zero principal moments of inertia: all three for a nonlinear
+        # 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 * PhysicalConstants["u"] * PhysicalConstants["A"] ** 2
+            self._eigenvalues_I[-self._n_rot:]
+            * PhysicalConstants["u"]
+            * PhysicalConstants["A"] ** 2
         )
         self._rotational_temperature = (
             PhysicalConstants["h"] ** 2 / (8 * np.pi**2 * PhysicalConstants["kB"])
@@ -222,18 +237,20 @@ def _compute_rotational_partition_function(self):
             * PhysicalConstants["HztoGHz"]
         )
 
-        self._rotational_temperature_xyz = (
-            self._rotational_temperature[0]
-            * self._rotational_temperature[1]
-            * self._rotational_temperature[2]
-        )
+        sigma = self._system.rotational_symmetry_number
 
-        self._rotational_partition_function = (
-            np.pi ** (1 / 2) / self._system.rotational_symmetry_number
-        ) * (
-            self._temperature ** (3 / 2)
-            / (np.power(self._rotational_temperature_xyz, 1 / 2))
-        )
+        if self._n_rot == 2:
+            # linear molecule: q_rot = T / (sigma * theta)
+            theta = self._rotational_temperature[-1]
+            self._rotational_temperature_xyz = theta
+            self._rotational_partition_function = self._temperature / (sigma * theta)
+        else:
+            # nonlinear molecule: q_rot = (sqrt(pi)/sigma) * T^(3/2) / sqrt(theta_xyz)
+            self._rotational_temperature_xyz = np.prod(self._rotational_temperature)
+            self._rotational_partition_function = (np.pi ** (1 / 2) / sigma) * (
+                self._temperature ** (3 / 2)
+                / (np.power(self._rotational_temperature_xyz, 1 / 2))
+            )
 
 
     def _compute_rotational_entropy(self):
@@ -244,9 +261,13 @@ def _compute_rotational_entropy(self):
         -------
         None
         """
+        if self._n_rot == 0:
+            self._rotational_entropy = 0.0
+            return
+
         self._rotational_entropy = (
             PhysicalConstants["R"]
-            * (np.log(self._rotational_partition_function) + 3 / 2)
+            * (np.log(self._rotational_partition_function) + self._n_rot / 2)
             / PhysicalConstants["cal"]
         )
 
@@ -260,15 +281,15 @@ def _compute_rotational_energy(self):
         None
         """
         self._rotational_energy = (
-            (3 / 2) * PhysicalConstants["R"] * self._temperature
+            (self._n_rot / 2) * PhysicalConstants["R"] * self._temperature
         ) / PhysicalConstants["cal"]
 
     def _compute_rotational_heat_capacity(self):
         """
         Computes the rotational heat capacity of the system.
         """
         self._rotational_heat_capacity = (
-            (3 / 2) * PhysicalConstants["R"] / PhysicalConstants["cal"]
+            (self._n_rot / 2) * PhysicalConstants["R"] / PhysicalConstants["cal"]
         )
 
 
@@ -283,12 +304,25 @@ def _rotational_contribution(self):
 
         self._relocate_to_cm()
         self._compute_inertia_tensor()
-        self._eigenvalues_I, self._eigenvectors_I = np.linalg.eig(self._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._n_rot = self._rotational_dof()
         self._compute_rotational_partition_function()
         self._compute_rotational_entropy()
         self._compute_rotational_energy()
         self._compute_rotational_heat_capacity()
 
+    def _rotational_dof(self) -> int:
+        """
+        Number of rotational degrees of freedom: 0 (monatomic), 2 (linear) or
+        3 (nonlinear).
+        """
+        atoms = self._system.atoms
+        if len(atoms) == 1:
+            return 0
+        return 2 if linearity(atoms) else 3
+
 
     def _compute_vibrational_partition_function(self):
         """

tests/thermo/test_thermo.py

@@ -115,6 +115,73 @@ def test_thermo_rejects_negative_pressure():
         Thermo(temperature=298.15, pressure=-1.0, system=_valid_system(), engine="dftb+")
 
 
+# --- geometry-dependent thermochemistry, validated against ASE IdealGasThermo --- #
+
+from ase import Atoms  # noqa: E402
+from ase.thermochemistry import IdealGasThermo  # noqa: E402
+
+_CM_TO_EV = 1.23984198e-4
+_EVK_TO_CALMOLK = 1.602176634e-19 * 6.02214076e23 / 4.184
+_T, _P = 298.15, 101325.0
+_R_CALMOLK = 8.314462618 / 4.184
+
+
+def _ts_thermo(symbols, positions, real_freqs, dof):
+    atoms = [Atom(symbol=s, position=np.array(p, float)) for s, p in zip(symbols, positions)]
+    pad = np.concatenate([np.zeros(3 * len(atoms) - dof), np.asarray(real_freqs, float)])
+    system = System(
+        atoms, periodicity=False, cell=None, charge=0,
+        electronic_energy=0.0, vibrational_frequencies=pad,
+    )
+    thermo = Thermo(temperature=_T, pressure=_P, system=system, engine="dftb+")
+    thermo.run()
+    return thermo
+
+
+def _ase_total_entropy(symbols, positions, real_freqs, geometry, sigma):
+    ase_thermo = IdealGasThermo(
+        vib_energies=[f * _CM_TO_EV for f in real_freqs],
+        geometry=geometry,
+        atoms=Atoms(symbols=symbols, positions=positions),
+        symmetrynumber=sigma, spin=0, potentialenergy=0.0,
+    )
+    return ase_thermo.get_entropy(_T, _P, verbose=False) * _EVK_TO_CALMOLK
+
+
+@pytest.mark.parametrize(
+    "symbols,positions,freqs,dof,geometry,sigma",
+    [
+        (["O", "H", "H"], [[0, 0, 0.119], [0, 0.763, -0.477], [0, -0.763, -0.477]],
+         [1595.0, 3657.0, 3756.0], 3, "nonlinear", 2),
+        (["C", "O", "O"], [[0, 0, 0], [0, 0, 1.16], [0, 0, -1.16]],
+         [667.0, 667.0, 1333.0, 2349.0], 4, "linear", 2),
+        (["N", "N"], [[0, 0, 0], [0, 0, 1.10]], [2359.0], 1, "linear", 2),
+        (["Ar"], [[0, 0, 0]], [], 0, "monatomic", 1),
+    ],
+)
+def test_total_entropy_matches_ase(symbols, positions, freqs, dof, geometry, sigma):
+    thermo = _ts_thermo(symbols, positions, freqs, dof)
+    ts_total = thermo.total_entropy("cal/(mol*K)")
+    ase_total = _ase_total_entropy(symbols, positions, freqs, geometry, sigma)
+
+    assert np.isfinite(ts_total)
+    assert ts_total == pytest.approx(ase_total, abs=0.05)
+
+
+def test_rotational_contribution_handles_linear_and_monatomic():
+    # linear and monatomic species used to give -inf rotational entropy
+    co2 = _ts_thermo(["C", "O", "O"], [[0, 0, 0], [0, 0, 1.16], [0, 0, -1.16]],
+                     [667.0, 667.0, 1333.0, 2349.0], 4)
+    assert co2._n_rot == 2
+    assert np.isfinite(co2._rotational_entropy)
+    assert co2._rotational_heat_capacity == pytest.approx(_R_CALMOLK)  # Cv_rot = R (linear)
+
+    argon = _ts_thermo(["Ar"], [[0, 0, 0]], [], 0)
+    assert argon._n_rot == 0
+    assert argon._rotational_entropy == 0.0
+    assert argon._rotational_heat_capacity == 0.0
+
+
 class TestThermo(unittest.TestCase):
 
     def test_thermo_aq_2(self):