Provide a mechanism to resolve imaginary modes

src/matcalc/_phonon.py

@@ -13,7 +13,7 @@
 from ._base import PropCalc
 from ._relaxation import RelaxCalc
 from .backend import run_pes_calc
-from .utils import to_pmg_structure
+from .utils import to_ase_atoms, to_pmg_structure
 
 if TYPE_CHECKING:
     from pathlib import Path
@@ -79,6 +79,9 @@ class PhononCalc(PropCalc):
         imaginary_freq_tol, then either raise a ValueError ("error") or log a
         warning ("warn").
     :type on_imaginary_modes: Literal["error", "warn"]
+    :ivar fix_imaginary_attempts: Number of attempts to resolve imaginary modes by rattling
+        the atoms and re-optimizing at fixed cell volume. 0 disables correction.
+    :type fix_imaginary_attempts: int
     :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.
@@ -114,6 +117,7 @@ def __init__(
         relax_calc_kwargs: dict | None = None,
         imaginary_freq_tol: float = -0.01,
         on_imaginary_modes: Literal["error", "warn"] = "warn",
+        fix_imaginary_attempts: int = 0,
         write_force_constants: bool | str | Path = False,
         write_band_structure: bool | str | Path = False,
         write_total_dos: bool | str | Path = False,
@@ -141,6 +145,9 @@ def __init__(
             a value below imaginary_freq_tol, it is considered imaginary.
         :param on_imaginary_modes: If there is a frequency with a value below imaginary_freq_tol, then
             raise a ValueError ("error") or log a warning ("warn"). Defaults to "warn".
+        :param fix_imaginary_attempts: Number of attempts to resolve imaginary modes. For each attempt, atoms
+            are rattled (ASE rattle, stdev=0.01 Å), the structure is re-optimized at fixed cell volume, and phonons
+            are recalculated. 0 disables correction.
         :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.
@@ -163,6 +170,7 @@ def __init__(
         self.relax_calc_kwargs = relax_calc_kwargs
         self.imaginary_freq_tol = imaginary_freq_tol
         self.on_imaginary_modes = on_imaginary_modes
+        self.fix_imaginary_attempts = fix_imaginary_attempts
         self.write_force_constants = write_force_constants
         self.write_band_structure = write_band_structure
         self.write_total_dos = write_total_dos
@@ -214,15 +222,18 @@ def calc(self, structure: Structure | Atoms | dict[str, Any]) -> dict:
 
         if self.relax_structure:
             logger.info("Relaxing input structure before phonon calculation")
-            relax_calc_kwargs = {"fmax": self.fmax, "optimizer": self.optimizer, "max_steps": self.max_steps} | (
-                self.relax_calc_kwargs or {}
-            )
-            relaxer = RelaxCalc(self.calculator, **relax_calc_kwargs)
-            result |= relaxer.calc(structure_in)
+            result |= self._relax_structure(structure_in)
             structure_in = result["final_structure"]
 
         phonon, frequencies, disp_supercells = self._run_phonopy(structure_in)
 
+        if self.fix_imaginary_attempts > 0 and np.any(frequencies < self.imaginary_freq_tol):
+            relax_result, phonon, frequencies, disp_supercells = self._resolve_imaginary_modes(
+                structure_in, phonon, frequencies, disp_supercells
+            )
+            result |= relax_result
+            structure_in = result["final_structure"]
+
         self._check_imaginary_modes(frequencies)
 
         phonon.run_thermal_properties(t_step=self.t_step, t_max=self.t_max, t_min=self.t_min)
@@ -251,7 +262,7 @@ def _run_phonopy(self, structure: Structure) -> tuple[phonopy.Phonopy, np.ndarra
             Tuple of (phonon, frequencies, disp_supercells).
         """
         cell = get_phonopy_structure(structure)
-        if self.supercell_matrix:
+        if self.supercell_matrix is not None:
             supercell_matrix = np.array(self.supercell_matrix, dtype=int)
         else:
             supercell_matrix = np.diag(np.ceil(self.min_length / np.array(structure.lattice.abc)).astype(int))
@@ -266,7 +277,7 @@ def _run_phonopy(self, structure: Structure) -> tuple[phonopy.Phonopy, np.ndarra
         logger.info("Computing forces on %d displaced supercells", len(disp_supercells))
         phonon.forces = [run_pes_calc(supercell, self.calculator).forces for supercell in disp_supercells]
         phonon.produce_force_constants()
-        phonon.run_mesh()
+        phonon.run_mesh(with_eigenvectors=True)
         frequencies = phonon.get_mesh_dict()["frequencies"]
         return phonon, frequencies, disp_supercells
 
@@ -282,13 +293,92 @@ def _check_imaginary_modes(self, frequencies: np.ndarray) -> None:
             n_imag = np.sum(imag_freq_mask)
             n_freqs = frequencies.size
             min_mode = np.min(frequencies)
-            pct = 100 * n_imag / n_freqs
+            pct_imag = 100 * n_imag / n_freqs
             msg = (
-                f"{n_imag}/{n_freqs} ({pct:.12}%) modes are imaginary (below {self.imaginary_freq_tol:.4f} THz). "
+                f"{n_imag}/{n_freqs} ({pct_imag:.12}%) modes are imaginary (below {self.imaginary_freq_tol:.4f} THz). "
                 f"Most negative: {min_mode:.2f} THz. This indicates a dynamically unstable structure. "
                 f"Thermal properties may not be reliable."
             )
             if self.on_imaginary_modes.lower() == "warn":
                 logger.warning(msg)
             else:
                 raise ValueError(msg)
+
+    def _resolve_imaginary_modes(
+        self,
+        structure_in: Structure,
+        phonon: phonopy.Phonopy,
+        frequencies: np.ndarray,
+        disp_supercells: list,
+    ) -> tuple[dict, phonopy.Phonopy, np.ndarray, list[Structure]]:
+        """Attempt to resolve imaginary modes by rattling and re-relaxing.
+
+        Args:
+            structure_in: Current primitive-cell structure.
+            phonon: Phonopy object from the most recent build.
+            frequencies: Frequencies from the most recent mesh run.
+            disp_supercells: Displaced supercells from the most recent build.
+
+        Returns:
+            Updated (RelaxCalc dictionary, phonon, frequencies, disp_supercells).
+        """
+        logger.warning(
+            "Imaginary modes detected (percent imaginary = %.2f%%, min freq: %.2f THz). "
+            "Attempting to resolve over %d attempt(s).",
+            100 * np.sum(frequencies < self.imaginary_freq_tol) / frequencies.size,
+            np.min(frequencies),
+            self.fix_imaginary_attempts,
+        )
+        relax_result: dict = {}
+        for attempt in range(self.fix_imaginary_attempts):
+            logger.info("Imaginary mode correction attempt %d/%d", attempt + 1, self.fix_imaginary_attempts)
+            structure_in = self._rattle_structure(structure_in)
+
+            logger.info("Re-relaxing structure at fixed cell volume following rattle.")
+            relax_result = self._relax_structure(structure_in)
+            structure_in = relax_result["final_structure"]
+
+            logger.info("Recomputing phonons after correction attempt %d", attempt + 1)
+            phonon, frequencies, disp_supercells = self._run_phonopy(structure_in)
+
+            if np.any(frequencies < self.imaginary_freq_tol):
+                logger.info(
+                    "Imaginary modes persist after attempt %d (percent imaginary = %.2f%%, min freq: %.2f THz).",
+                    attempt + 1,
+                    100 * np.sum(frequencies < self.imaginary_freq_tol) / frequencies.size,
+                    np.min(frequencies),
+                )
+            else:
+                logger.info("Imaginary modes resolved after %d attempt(s)", attempt + 1)
+                break
+        return relax_result, phonon, frequencies, disp_supercells
+
+    def _rattle_structure(self, structure_in: Structure, stdev: float = 0.01) -> Structure:
+        """Rattle the atoms to bump out of stationary point (stdev=0.01 Å).
+
+        Args:
+            structure_in: Pymatgen structure.
+            stdev: standard deviation in angstrom for the rattle.
+
+        Returns:
+            Pymatgen structure with rattled atomic positions.
+        """
+        atoms = to_ase_atoms(structure_in)
+        atoms.rattle(stdev=stdev)
+        atoms.wrap()
+        return to_pmg_structure(atoms)
+
+    def _relax_structure(self, structure_in: Structure) -> dict:
+        """Relax a structure.
+
+        Args:
+            structure_in: Pymatgen structure.
+
+        Returns:
+            Full result dict from RelaxCalc.
+        """
+        relax_calc_kwargs = {"fmax": self.fmax, "optimizer": self.optimizer, "max_steps": self.max_steps} | (
+            self.relax_calc_kwargs or {}
+        )
+        relaxer = RelaxCalc(self.calculator, **relax_calc_kwargs)
+        return relaxer.calc(structure_in)

src/matcalc/_qha.py

@@ -52,6 +52,9 @@ class QHACalc(PropCalc):
     :type optimizer: str
     :ivar eos: Equation of state used for fitting energy vs. volume data.
     :type eos: Literal["vinet", "birch_murnaghan", "murnaghan"]
+    :ivar allow_shape_change: Whether or not to allow the unit cell shape to
+        change at fixed cell volume during the EOS calculations. Default is True.
+    :type allow_shape_change: bool
     :ivar relax_structure: Whether to perform structure relaxation before phonon calculations.
     :type relax_structure: bool
     :ivar relax_calc_kwargs: Additional keyword arguments for structure relaxation calculations.
@@ -67,6 +70,9 @@ class QHACalc(PropCalc):
         imaginary_freq_tol, then either raise a ValueError ("error") or log a
         warning ("warn").
     :type on_imaginary_modes: Literal["error", "warn"]
+    :ivar fix_imaginary_attempts: Number of attempts passed to PhononCalc to resolve imaginary modes
+        at each scale factor. 0 disables correction.
+    :type fix_imaginary_attempts: int
     :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.
@@ -99,12 +105,14 @@ def __init__(
         max_steps: int = 5000,
         optimizer: str = "FIRE",
         eos: Literal["vinet", "birch_murnaghan", "murnaghan"] = "vinet",
+        allow_shape_change: bool = True,
         relax_structure: bool = True,
         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)),
+        scale_factors: Sequence[float] = (0.95, 0.96, 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03, 1.04, 1.05),
         imaginary_freq_tol: float = -0.01,
         on_imaginary_modes: Literal["error", "warn"] = "warn",
+        fix_imaginary_attempts: int = 0,
         write_helmholtz_volume: bool | str | os.PathLike = False,
         write_volume_temperature: bool | str | os.PathLike = False,
         write_thermal_expansion: bool | str | os.PathLike = False,
@@ -132,6 +140,8 @@ def __init__(
             "FIRE".
         :param eos: Equation of state to use for calculating energy vs. volume relationships.
             Default is "vinet".
+        :param allow_shape_change: Whether or not to allow the unit cell shape to
+            change at fixed cell volume during the EOS calculations. Default is True.
         :param relax_structure: A boolean flag indicating whether the initial atomic structure should be
             relaxed as part of the computation workflow. Note that subsequent relaxations on the
             volume-scaled structures will be carried out regardless.
@@ -145,6 +155,8 @@ def __init__(
             a value below imaginary_freq_tol, it is considered imaginary.
         :param on_imaginary_modes: If there is a frequency with a value below imaginary_freq_tol, then
             raise a ValueError ("error") or log a warning ("warn"). Defaults to "warn".
+        :param fix_imaginary_attempts: Number of attempts passed to PhononCalc to resolve imaginary modes
+            at each scale factor. 0 disables correction.
         :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
@@ -172,12 +184,26 @@ def __init__(
         self.max_steps = max_steps
         self.optimizer = optimizer
         self.eos = eos
+        self.allow_shape_change = allow_shape_change
         self.relax_structure = relax_structure
         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.on_imaginary_modes = on_imaginary_modes
+        self.fix_imaginary_attempts = fix_imaginary_attempts
+
+        # Needed to make sure the volume doesn't change during relaxations on scaled structures
+        self._fixed_cell_relax_calc_kwargs = {
+            "optimizer": self.optimizer,
+            "fmax": self.fmax,
+            "max_steps": self.max_steps,
+            **(self.relax_calc_kwargs or {}),
+        }
+        self._fixed_cell_relax_calc_kwargs["relax_cell"] = bool(self.allow_shape_change)
+        self._fixed_cell_relax_calc_kwargs["cell_filter_kwargs"] = (
+            {"constant_volume": True} if self.allow_shape_change else {}
+        )
 
         # Normalize write_* inputs to Optional[str | os.PathLike]:
         # - True  -> default filename (meaning "write to default file")
@@ -236,11 +262,13 @@ def calc(self, structure: Structure | Atoms | dict[str, Any]) -> dict:
 
         if self.relax_structure:
             logger.info("Relaxing input structure before QHA")
-            relax_calc_kwargs = {"fmax": self.fmax, "optimizer": self.optimizer, "max_steps": self.max_steps} | (
-                self.relax_calc_kwargs or {}
-            )
-            relaxer = RelaxCalc(self.calculator, **relax_calc_kwargs)
-            result |= relaxer.calc(structure_in)
+            result |= RelaxCalc(
+                self.calculator,
+                fmax=self.fmax,
+                optimizer=self.optimizer,
+                max_steps=self.max_steps,
+                **self.relax_calc_kwargs or {},
+            ).calc(structure_in)
             structure_in = result["final_structure"]
 
         logger.info(
@@ -299,30 +327,36 @@ def _collect_properties(self, structure: Structure) -> dict[str, list]:
         scaled_structures = []
         for scale_factor in self.scale_factors:
             # Apply linear strain
-            struct = self._scale_structure(structure, scale_factor)
-            volumes.append(struct.volume)
+            scaled_structure = self._scale_structure(structure, scale_factor)
 
             # Relax at fixed volume
             logger.info("Scale factor %.3f: relaxing at fixed volume", scale_factor)
-            relaxer = RelaxCalc(
-                self.calculator,
-                optimizer=self.optimizer,
-                fmax=self.fmax,
-                max_steps=self.max_steps,
-                relax_cell=False,
-                **(self.relax_calc_kwargs or {}),
-            )
-            relaxed_result = relaxer.calc(struct)
-            electronic_energies.append(relaxed_result["energy"])
-            scaled_structures.append(relaxed_result["final_structure"])
+            relaxer = RelaxCalc(self.calculator, **self._fixed_cell_relax_calc_kwargs)
+            relaxed_result = relaxer.calc(scaled_structure)
 
-            # Calculate thermal properties from phonon calculation
+            # Calculate thermal properties from phonon calculation and tabulate results.
+            # We must store the energy, structure, etc. from the phonon calculation if
+            # available in case there are any correction attempts made to resolve imaginary modes,
+            # which would update the relaxation output.
             logger.info(
                 "Scale factor %.3f: computing phonon thermal properties (V=%.1f ų)",
                 scale_factor,
                 relaxed_result["final_structure"].volume,
             )
             phonon_result = self._calculate_thermal_properties(relaxed_result["final_structure"])
+
+            # Collect properties
+            scaled_structures.append(phonon_result["final_structure"])
+            volume = phonon_result["final_structure"].volume
+            if (
+                np.abs(volume - scaled_structure.volume) / scaled_structure.volume > 1e-4  # noqa: PLR2004
+            ):
+                raise ValueError(
+                    f"Somehow the volume changed during relaxation. This is a bug! Before: {scaled_structure.volume}. "
+                    f"After: {volume}"
+                )
+            volumes.append(volume)
+            electronic_energies.append(phonon_result.get("energy", relaxed_result.get("energy")))
             thermal_properties = phonon_result["thermal_properties"]
             free_energies.append(thermal_properties["free_energy"])
             entropies.append(thermal_properties["entropy"])
@@ -365,11 +399,13 @@ def _calculate_thermal_properties(self, structure: Structure) -> dict:
             "t_step": self.t_step,
             "t_max": self.t_max,
             "t_min": self.t_min,
-            "relax_structure": False,
+            "relax_calc_kwargs": self._fixed_cell_relax_calc_kwargs,  # needed for _resolve_imaginary_modes
             "imaginary_freq_tol": self.imaginary_freq_tol,
             "on_imaginary_modes": self.on_imaginary_modes,
+            "fix_imaginary_attempts": self.fix_imaginary_attempts,
             "write_phonon": False,
         } | (self.phonon_calc_kwargs or {})
+        phonon_calc_kwargs["relax_structure"] = False  # already relaxed
         phonon_calc = PhononCalc(
             self.calculator,
             **phonon_calc_kwargs,

tests/test_phonon.py

@@ -193,3 +193,44 @@ def test_phonon_calc_imaginary_freq_tol(
     result = phonon_calc.calc(distorted_si_atoms)
     assert "frequencies" in result
     assert np.min(result["frequencies"]) < -0.1
+
+
+def test_phonon_calc_fix_imaginary_attempts(
+    Si_atoms: Atoms,
+    matpes_calculator: PESCalculator,
+    caplog: pytest.LogCaptureFixture,
+) -> None:
+    """Test fix_imaginary_attempts parameter for resolving imaginary modes."""
+    # Stable structure which needs no fixing — correction attempt should NOT be logged
+    phonon_calc = PhononCalc(
+        calculator=matpes_calculator,
+        supercell_matrix=((2, 0, 0), (0, 2, 0), (0, 0, 2)),
+        fmax=0.1,
+        t_step=50,
+        t_max=1000,
+        fix_imaginary_attempts=1,
+        imaginary_freq_tol=-0.01,
+        on_imaginary_modes="error",
+    )
+    with caplog.at_level(logging.INFO, logger="matcalc"):
+        result = phonon_calc.calc(Si_atoms)
+    assert np.min(result["frequencies"]) > -0.01
+    assert not any("Imaginary mode correction attempt" in r.message for r in caplog.records)
+
+    caplog.clear()
+
+    # 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,
+        on_imaginary_modes="error",
+        fix_imaginary_attempts=1,
+    )
+    with caplog.at_level(logging.INFO, logger="matcalc"), pytest.raises(ValueError, match="modes are imaginary"):
+        phonon_calc.calc(distorted_si_atoms)
+    assert any("Imaginary mode correction attempt" in r.message for r in caplog.records)