diff --git a/ThermoScreening/thermo/system.py b/ThermoScreening/thermo/system.py index 6cfee9f..0f16e64 100644 --- a/ThermoScreening/thermo/system.py +++ b/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) diff --git a/ThermoScreening/thermo/thermo.py b/ThermoScreening/thermo/thermo.py index 1f98158..670ec8f 100644 --- a/ThermoScreening/thermo/thermo.py +++ b/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,7 +281,7 @@ 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): @@ -268,7 +289,7 @@ 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): """ diff --git a/tests/thermo/test_thermo.py b/tests/thermo/test_thermo.py index 99a2896..aa5d7e9 100644 --- a/tests/thermo/test_thermo.py +++ b/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):