Add optional check for imaginary modes in PhononCalc/QHACalc and return the frequencies

.gitignore

@@ -24,6 +24,7 @@ setuptools*
 .eggs/
 .coverage
 .*_cache
+.claude/*
 # VS Code
 .vscode/*
 _site

src/matcalc/_phonon.py

@@ -4,6 +4,7 @@
 
 from typing import TYPE_CHECKING
 
+import numpy as np
 import phonopy
 from phonopy.file_IO import write_FORCE_CONSTANTS as write_force_constants
 from pymatgen.io.phonopy import get_phonopy_structure, get_pmg_structure
@@ -63,6 +64,13 @@ class PhononCalc(PropCalc):
     :ivar relax_calc_kwargs: Optional dictionary containing additional
         arguments for the structural relaxation calculation.
     :type relax_calc_kwargs: dict | None
+    :ivar imaginary_freq_tol: Tolerance for imaginary frequency detection
+        in THz. If set to a float value, the calculator will raise a
+        ValueError when imaginary frequencies with magnitude exceeding this
+        threshold are found. Imaginary frequencies are represented as
+        negative values in phonopy. A value of None (default) disables
+        the check.
+    :type imaginary_freq_tol: float | None
     :ivar write_force_constants: Path, boolean, or string specifying whether
         to write the calculated force constants to an output file, and the
         path or name of the file if applicable.
@@ -94,6 +102,7 @@ def __init__(
         optimizer: str = "FIRE",
         relax_structure: bool = True,
         relax_calc_kwargs: dict | None = None,
+        imaginary_freq_tol: float | None = None,
         write_force_constants: bool | str | Path = False,
         write_band_structure: bool | str | Path = False,
         write_total_dos: bool | str | Path = False,
@@ -114,6 +123,9 @@ def __init__(
         :param optimizer: Name of the optimization algorithm for structural relaxation.
         :param relax_structure: Flag to indicate whether structure relaxation should be performed before calculations.
         :param relax_calc_kwargs: Additional keyword arguments for relaxation phase calculations.
+        :param imaginary_freq_tol: Tolerance for imaginary frequency detection in THz.
+            If a positive float, a ValueError is raised when any imaginary frequency with
+            magnitude exceeding this value is found. Defaults to None (no check).
         :param write_force_constants: File path or boolean flag to write force constants.
             Defaults to "force_constants".
         :param write_band_structure: File path or boolean flag to write band structure data.
@@ -132,6 +144,7 @@ def __init__(
         self.optimizer = optimizer
         self.relax_structure = relax_structure
         self.relax_calc_kwargs = relax_calc_kwargs
+        self.imaginary_freq_tol = imaginary_freq_tol
         self.write_force_constants = write_force_constants
         self.write_band_structure = write_band_structure
         self.write_total_dos = write_total_dos
@@ -195,6 +208,22 @@ def calc(self, structure: Structure | Atoms | dict[str, Any]) -> dict:
         ]
         phonon.produce_force_constants()
         phonon.run_mesh()
+        mesh_dict_results = phonon.get_mesh_dict()
+        frequencies = mesh_dict_results["frequencies"]
+
+        if self.imaginary_freq_tol is not None:
+            # In phonopy, imaginary frequencies are represented as negative values.
+            imag_freq_mask = frequencies < -self.imaginary_freq_tol
+            n_freqs = len(frequencies)
+            if np.any(imag_freq_mask):
+                n_imag = np.sum(imag_freq_mask)
+                min_mode = np.min(frequencies)
+                raise ValueError(
+                    f"{n_imag}/{n_freqs} are imaginary (below -{self.imaginary_freq_tol:.4f} THz). "
+                    f"Most negative: {min_mode:.4f} THz. This indicates a dynamically unstable structure. "
+                    f"Thermal properties may not be reliable."
+                )
+
         phonon.run_thermal_properties(t_step=self.t_step, t_max=self.t_max, t_min=self.t_min)
         if self.write_force_constants:
             write_force_constants(phonon.force_constants, filename=self.write_force_constants)  # type: ignore[arg-type]
@@ -204,4 +233,8 @@ def calc(self, structure: Structure | Atoms | dict[str, Any]) -> dict:
             phonon.auto_total_dos(write_dat=True, filename=self.write_total_dos)
         if self.write_phonon:
             phonon.save(filename=self.write_phonon)  # type: ignore[arg-type]
-        return result | {"phonon": phonon, "thermal_properties": phonon.get_thermal_properties_dict()}
+        return result | {
+            "phonon": phonon,
+            "thermal_properties": phonon.get_thermal_properties_dict(),
+            "frequencies": frequencies,
+        }

src/matcalc/_qha.py

@@ -56,6 +56,12 @@ class QHACalc(PropCalc):
     :type phonon_calc_kwargs: dict | None
     :ivar scale_factors: List of scale factors for lattice scaling.
     :type scale_factors: Sequence[float]
+    :ivar imaginary_freq_tol: Tolerance for imaginary frequency detection
+        in THz. Passed through to the internal PhononCalc instances. If set
+        to a float value, the calculator will raise a ValueError when
+        imaginary frequencies with magnitude exceeding this threshold are
+        found. A value of None (default) disables the check.
+    :type imaginary_freq_tol: float | None
     :ivar write_helmholtz_volume: Path or boolean to control saving Helmholtz free energy vs. volume data.
     :type write_helmholtz_volume: bool | str | Path
     :ivar write_volume_temperature: Path or boolean to control saving volume vs. temperature data.
@@ -91,6 +97,7 @@ def __init__(
         relax_calc_kwargs: dict | None = None,
         phonon_calc_kwargs: dict | None = None,
         scale_factors: Sequence[float] = tuple(np.arange(0.95, 1.05, 0.01)),
+        imaginary_freq_tol: float | None = None,
         write_helmholtz_volume: bool | str | Path = False,
         write_volume_temperature: bool | str | Path = False,
         write_thermal_expansion: bool | str | Path = False,
@@ -125,6 +132,10 @@ def __init__(
             phonon calculation routine.
         :param scale_factors: A sequence of scale factors for volume scaling during
             thermodynamic and phononic calculations.
+        :param imaginary_freq_tol: Tolerance for imaginary frequency detection in THz.
+            Passed through to internal PhononCalc instances. If a positive float, a ValueError
+            is raised when any imaginary frequency with magnitude exceeding this value is found.
+            Defaults to None, which means no check will be carried out.
         :param write_helmholtz_volume: Path, boolean, or string to indicate whether and where
             to save Helmholtz energy as a function of volume.
         :param write_volume_temperature: Path, boolean, or string to indicate whether and where
@@ -155,6 +166,7 @@ def __init__(
         self.relax_calc_kwargs = relax_calc_kwargs
         self.phonon_calc_kwargs = phonon_calc_kwargs
         self.scale_factors = scale_factors
+        self.imaginary_freq_tol = imaginary_freq_tol
         self.write_helmholtz_volume = write_helmholtz_volume
         self.write_volume_temperature = write_volume_temperature
         self.write_thermal_expansion = write_thermal_expansion
@@ -314,6 +326,7 @@ def _calculate_thermal_properties(self, structure: Structure) -> dict:
             relax_structure=True,
             relax_calc_kwargs={"relax_cell": False},
             write_phonon=False,
+            imaginary_freq_tol=self.imaginary_freq_tol,
             **(self.phonon_calc_kwargs or {}),
         )
         return phonon_calc.calc(structure)["thermal_properties"]

tests/test_lammps.py

@@ -93,7 +93,7 @@ def test_invalid_ensemble(Si: Structure) -> None:
         ValueError,
         match=(
             "The specified ensemble is not supported, choose from 'nve', 'nvt',"
-            " 'nvt_nose_hoover', 'npt', 'npt_nose_hoover'."
+            r" 'nvt_nose_hoover', 'npt', 'npt_nose_hoover'."
         ),
     ):
         LAMMPSMDCalc(

tests/test_phonon.py

@@ -6,6 +6,7 @@
 import os
 from typing import TYPE_CHECKING
 
+import numpy as np
 import pytest
 
 from matcalc import PhononCalc
@@ -95,3 +96,46 @@ def test_phonon_calc_atoms(
     thermal_props = result["thermal_properties"]
     ind = thermal_props["temperatures"].tolist().index(300)
     assert thermal_props["heat_capacity"][ind] == pytest.approx(43.3138042001517, rel=1e-1)
+
+
+def test_phonon_calc_imaginary_freq_tol(
+    Si_atoms: Atoms,
+    matpes_calculator: PESCalculator,
+) -> None:
+    """Test that imaginary frequency tolerance check works correctly."""
+    # Stable structure with tolerance should pass without error
+    phonon_calc = PhononCalc(
+        calculator=matpes_calculator,
+        supercell_matrix=((2, 0, 0), (0, 2, 0), (0, 0, 2)),
+        fmax=0.1,
+        imaginary_freq_tol=0.1,
+    )
+    result = phonon_calc.calc(Si_atoms)
+    assert "frequencies" in result
+    assert np.min(result["frequencies"]) >= -0.1
+
+    # Distort the structure to create imaginary modes, then check that
+    # the tolerance check raises ValueError
+    distorted_si_atoms = Si_atoms.copy()
+    distorted_si_atoms.cell += 0.5
+    phonon_calc = PhononCalc(
+        calculator=matpes_calculator,
+        supercell_matrix=((2, 0, 0), (0, 2, 0), (0, 0, 2)),
+        fmax=100.0,
+        imaginary_freq_tol=0.1,
+    )
+    with pytest.raises(ValueError, match=r"are imaginary"):
+        phonon_calc.calc(distorted_si_atoms)
+
+    # Distort the structure to create imaginary modes, but it's okay
+    distorted_si_atoms = Si_atoms.copy()
+    distorted_si_atoms.cell += 0.5
+    phonon_calc = PhononCalc(
+        calculator=matpes_calculator,
+        supercell_matrix=((2, 0, 0), (0, 2, 0), (0, 0, 2)),
+        fmax=100.0,
+        imaginary_freq_tol=None,
+    )
+    result = phonon_calc.calc(distorted_si_atoms)
+    assert "frequencies" in result
+    assert np.min(result["frequencies"]) < -0.1

tests/test_qha.py

@@ -191,3 +191,56 @@ def test_qha_calc_atoms(
     # Test values at 300 K
     ind = result["temperatures"].tolist().index(300)
     assert result["thermal_expansion_coefficients"][ind] == pytest.approx(5.191273165438463e-06, rel=1e-1)
+
+
+def test_phonon_calc_imaginary_freq_tol(
+    Si_atoms: Atoms,
+    matpes_calculator: PESCalculator,
+) -> None:
+
+    # Initialize QHACalc
+    qha_calc = QHACalc(
+        calculator=matpes_calculator,
+        t_step=50,
+        t_max=1000,
+        scale_factors=[0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03],
+        fmax=0.1,
+        imaginary_freq_tol=0.1,
+        phonon_calc_kwargs={"supercell_matrix": ((2, 0, 0), (0, 2, 0), (0, 0, 2))},
+    )
+
+    result = qha_calc.calc(Si_atoms)
+
+    ind = result["temperatures"].tolist().index(300)
+    assert result["thermal_expansion_coefficients"][ind] == pytest.approx(5.191273165438463e-06, rel=1e-1)
+
+    # Distorted
+    distorted_Si_atoms = Si_atoms.copy()
+    distorted_Si_atoms.cell += 0.5
+    qha_calc = QHACalc(
+        calculator=matpes_calculator,
+        t_step=50,
+        t_max=1000,
+        scale_factors=[0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03],
+        fmax=100,
+        imaginary_freq_tol=0.1,
+        phonon_calc_kwargs={"supercell_matrix": ((2, 0, 0), (0, 2, 0), (0, 0, 2))},
+    )
+
+    with pytest.raises(ValueError, match="are imaginary"):
+        qha_calc.calc(distorted_Si_atoms)
+
+    # Distorted no check
+    distorted_Si_atoms = Si_atoms.copy()
+    distorted_Si_atoms.cell += 0.5
+    qha_calc = QHACalc(
+        calculator=matpes_calculator,
+        t_step=50,
+        t_max=1000,
+        scale_factors=[0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03],
+        fmax=100,
+        imaginary_freq_tol=None,
+        phonon_calc_kwargs={"supercell_matrix": ((2, 0, 0), (0, 2, 0), (0, 0, 2))},
+    )
+
+    assert qha_calc.calc(distorted_Si_atoms)