Skip to content

Imaginary vibrational modes silently produce NaN thermochemistry #60

Description

@galjos

An imaginary (negative-wavenumber) mode in the kept dof set makes the harmonic-oscillator formulas in thermo.py (_compute_vibrational_*, ~lines 293-398) return nan for entropy (np.log(1 - exp(-theta/T)) of a negative argument) and a downward-corrupted ZPE/internal energy. System.has_imaginary_frequencies is computed but Thermo never consults it, so the result is a silent nan total entropy / Gibbs rather than an error.

Repro:

import numpy as np
from ThermoScreening.thermo.atoms import Atom
from ThermoScreening.thermo.system import System
from ThermoScreening.thermo.thermo import Thermo

atoms = [Atom(symbol=s, position=np.array(p, float)) for s, p in
         zip(["O","H","H"], [[0,0,0.119],[0,0.763,-0.477],[0,-0.763,-0.477]])]
system = System(atoms, periodicity=False, cell=None, charge=0,
                electronic_energy=0.0, vibrational_frequencies=np.array([-500.0, 3657.0, 3756.0]))
t = Thermo(temperature=298.15, pressure=101325, system=system, engine="dftb+"); t.run()
print(system.has_imaginary_frequencies, t.total_entropy("cal/(mol*K)"))
# True nan

Imaginary modes are common in real DFTB+ output (transition states, unconverged minima), so silent NaN is a real hazard for screening.

Suggested fix: have _vibrational_contribution consult system.has_imaginary_frequencies and raise a TSValueError ("Imaginary vibrational frequencies are present.") when an imaginary mode is in the kept dof set (or explicitly exclude/floor them, with a documented policy). The HO formulas themselves are correct. Found alongside #59 (rotational geometry fix).

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions