Add support for having non-displaced atoms in Phonopy routines

src/quacc/atoms/phonons.py

@@ -9,12 +9,12 @@
 from monty.dev import requires
 from pymatgen.io.ase import AseAtomsAdaptor
 from pymatgen.io.phonopy import get_phonopy_structure, get_pmg_structure
-from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 
 has_phonopy = find_spec("phonopy")
 
 if has_phonopy:
     from phonopy import Phonopy
+    from phonopy.structure.cells import get_supercell
 
 if TYPE_CHECKING:
     from ase.atoms import Atoms
@@ -26,7 +26,7 @@
 @requires(has_phonopy, "Phonopy not installed.")
 def get_phonopy(
     atoms: Atoms,
-    min_lengths: float | tuple[float, float, float] | None = None,
+    min_lengths: float | tuple[float, float, float] | None = 20.0,
     supercell_matrix: (
         tuple[tuple[int, int, int], tuple[int, int, int], tuple[int, int, int]] | None
     ) = None,
@@ -60,22 +60,21 @@ def get_phonopy(
     """
     phonopy_kwargs = phonopy_kwargs or {}
 
-    symmetrized_structure = SpacegroupAnalyzer(
-        AseAtomsAdaptor().get_structure(atoms), symprec=symprec
-    ).get_symmetrized_structure()
-
     if supercell_matrix is None and min_lengths is not None:
         supercell_matrix = np.diag(
-            np.round(np.ceil(min_lengths / np.array(symmetrized_structure.lattice.abc)))
+            np.round(np.ceil(min_lengths / atoms.cell.lengths()))
         )
 
+    structure = AseAtomsAdaptor.get_structure(atoms)
+
     phonon = Phonopy(
-        get_phonopy_structure(symmetrized_structure),
+        get_phonopy_structure(structure),
         symprec=symprec,
         supercell_matrix=supercell_matrix,
         **phonopy_kwargs,
     )
     phonon.generate_displacements(distance=displacement)
+
     return phonon
 
 
@@ -95,3 +94,33 @@ def phonopy_atoms_to_ase_atoms(phonpy_atoms: PhonopyAtoms) -> Atoms:
     """
     pmg_structure = get_pmg_structure(phonpy_atoms)
     return pmg_structure.to_ase_atoms()
+
+
+def get_atoms_supercell_by_phonopy(
+    atoms: Atoms,
+    supercell_matrix: tuple[
+        tuple[int, int, int], tuple[int, int, int], tuple[int, int, int]
+    ],
+) -> Atoms:
+    """
+    Get the supercell of an ASE atoms object using a supercell matrix.
+
+    Parameters
+    ----------
+    atoms
+        ASE atoms object.
+    supercell_matrix
+        The supercell matrix to use. If specified, it will override any
+        value specified by `min_lengths`.
+    Returns
+    -------
+    Atoms
+        ASE atoms object of the supercell.
+    """
+
+    return phonopy_atoms_to_ase_atoms(
+        get_supercell(
+            get_phonopy_structure(AseAtomsAdaptor.get_structure(atoms)),
+            supercell_matrix,
+        )
+    )

src/quacc/recipes/common/phonons.py

@@ -5,10 +5,15 @@
 from importlib.util import find_spec
 from typing import TYPE_CHECKING
 
+from ase.atoms import Atoms
 from monty.dev import requires
 
 from quacc import job, subflow
-from quacc.atoms.phonons import get_phonopy, phonopy_atoms_to_ase_atoms
+from quacc.atoms.phonons import (
+    get_atoms_supercell_by_phonopy,
+    get_phonopy,
+    phonopy_atoms_to_ase_atoms,
+)
 from quacc.runners.phonons import run_phonopy
 from quacc.schemas.phonons import summarize_phonopy
 
@@ -17,22 +22,21 @@
 if TYPE_CHECKING:
     from typing import Any
 
-    from ase.atoms import Atoms
-
     from quacc import Job
     from quacc.schemas._aliases.phonons import PhononSchema
 
     if has_deps:
-        from phonopy import Phonopy
+        pass
 
 
 @subflow
 @requires(
     has_deps, "Phonopy and seekpath must be installed. Run `pip install quacc[phonons]`"
 )
 def phonon_subflow(
-    atoms: Atoms,
+    displaced_atoms: Atoms,
     force_job: Job,
+    non_displaced_atoms: Atoms | None = None,
     symprec: float = 1e-4,
     min_lengths: float | tuple[float, float, float] | None = 20.0,
     supercell_matrix: (
@@ -48,12 +52,22 @@ def phonon_subflow(
     """
     Calculate phonon properties.
 
+    In Quacc the ASE constraints can be used to fix atoms. These atoms will
+    not be displaced during the phonon calculation. This will greatly reduce
+    the computational cost of the calculation. However, this is an important
+    approximation and should be used with caution.
+
     Parameters
     ----------
-    atoms
+    displaced_atoms
         Atoms object with calculator attached.
     force_job
         The static job to calculate the forces.
+    non_displaced_atoms
+        Additional atoms to add to the supercells i.e. fixed atoms.
+        These atoms will not be displaced during the phonon calculation.
+        Useful for adsorbates on surfaces with weak coupling etc.
+        Important approximation, use with caution.
     symprec
         Precision for symmetry detection.
     min_lengths
@@ -80,6 +94,27 @@ def phonon_subflow(
         Dictionary of results from [quacc.schemas.phonons.summarize_phonopy][]
     """
 
+    non_displaced_atoms = non_displaced_atoms or Atoms()
+
+    phonopy = get_phonopy(
+        displaced_atoms,
+        min_lengths=min_lengths,
+        supercell_matrix=supercell_matrix,
+        symprec=symprec,
+        displacement=displacement,
+        phonopy_kwargs=phonopy_kwargs,
+    )
+
+    if non_displaced_atoms:
+        non_displaced_atoms = get_atoms_supercell_by_phonopy(
+            non_displaced_atoms, phonopy.supercell_matrix
+        )
+
+    supercells = [
+        phonopy_atoms_to_ase_atoms(s) + non_displaced_atoms
+        for s in phonopy.supercells_with_displacements
+    ]
+
     @subflow
     def _get_forces_subflow(supercells: list[Atoms]) -> list[dict]:
         return [
@@ -89,17 +124,25 @@ def _get_forces_subflow(supercells: list[Atoms]) -> list[dict]:
     @job
     def _thermo_job(
         atoms: Atoms,
-        phonopy: Phonopy,
+        phonopy,
         force_job_results: list[dict],
-        t_step: float,
-        t_min: float,
-        t_max: float,
-        additional_fields: dict[str, Any] | None,
+        t_step,
+        t_min,
+        t_max,
+        additional_fields,
     ) -> PhononSchema:
         parameters = force_job_results[-1].get("parameters")
-        forces = [output["results"]["forces"] for output in force_job_results]
+        forces = [
+            output["results"]["forces"][: len(phonopy.supercell)]
+            for output in force_job_results
+        ]
         phonopy_results = run_phonopy(
-            phonopy, forces, t_step=t_step, t_min=t_min, t_max=t_max
+            phonopy,
+            forces,
+            symmetrize=bool(non_displaced_atoms),
+            t_step=t_step,
+            t_min=t_min,
+            t_max=t_max,
         )
 
         return summarize_phonopy(
@@ -110,18 +153,13 @@ def _thermo_job(
             additional_fields=additional_fields,
         )
 
-    phonopy = get_phonopy(
-        atoms,
-        min_lengths=min_lengths,
-        supercell_matrix=supercell_matrix,
-        symprec=symprec,
-        displacement=displacement,
-        phonopy_kwargs=phonopy_kwargs,
-    )
-    supercells = [
-        phonopy_atoms_to_ase_atoms(s) for s in phonopy.supercells_with_displacements
-    ]
     force_job_results = _get_forces_subflow(supercells)
     return _thermo_job(
-        atoms, phonopy, force_job_results, t_step, t_min, t_max, additional_fields
+        displaced_atoms,
+        phonopy,
+        force_job_results,
+        t_step,
+        t_min,
+        t_max,
+        additional_fields,
     )

src/quacc/recipes/emt/phonons.py

@@ -20,13 +20,14 @@
 
 @flow
 def phonon_flow(
-    atoms: Atoms,
+    displaced_atoms: Atoms,
     symprec: float = 1e-4,
     min_lengths: float | tuple[float, float, float] | None = 20.0,
     supercell_matrix: (
         tuple[tuple[int, int, int], tuple[int, int, int], tuple[int, int, int]] | None
     ) = None,
     displacement: float = 0.01,
+    non_displaced_atoms: Atoms | None = None,
     t_step: float = 10,
     t_min: float = 0,
     t_max: float = 1000,
@@ -51,7 +52,7 @@ def phonon_flow(
 
     Parameters
     ----------
-    atoms
+    displaced_atoms
         Atoms object
     symprec
         Precision for symmetry detection.
@@ -62,6 +63,11 @@ def phonon_flow(
         value specified by `min_lengths`.
     displacement
         Atomic displacement (A).
+    non_displaced_atoms
+        Additional atoms to add to the supercells i.e. fixed atoms.
+        These atoms will not be displaced during the phonon calculation.
+        Useful for adsorbates on surfaces with weak coupling etc.
+        Important approximation, use with caution.
     t_step
         Temperature step (K).
     t_min
@@ -92,15 +98,16 @@ def phonon_flow(
         parameters=job_params,
         decorators=job_decorators,
     )
-    if run_relax:
+    if run_relax and not non_displaced_atoms:
         atoms = relax_job_(atoms)["atoms"]
 
     return phonon_subflow(
-        atoms,
+        displaced_atoms,
         static_job_,
         symprec=symprec,
         min_lengths=min_lengths,
         supercell_matrix=supercell_matrix,
+        non_displaced_atoms=non_displaced_atoms,
         displacement=displacement,
         t_step=t_step,
         t_min=t_min,

src/quacc/runners/phonons.py

@@ -26,6 +26,7 @@
 def run_phonopy(
     phonon: Phonopy,
     forces: NDArray,
+    symmetrize: bool = False,
     t_step: float = 10,
     t_min: float = 0,
     t_max: float = 1000,
@@ -40,6 +41,8 @@ def run_phonopy(
         Phonopy object
     forces
         Forces on the atoms
+    symmetrize
+        Whether to symmetrize the force constants
     t_step
         Temperature step
     t_min
@@ -60,6 +63,11 @@ def run_phonopy(
     # Run phonopy
     phonon.forces = forces
     phonon.produce_force_constants()
+
+    if symmetrize:
+        phonon.symmetrize_force_constants()
+        phonon.symmetrize_force_constants_by_space_group()
+
     phonon.run_mesh(with_eigenvectors=True)
     phonon.run_total_dos()
     phonon.run_thermal_properties(t_step=t_step, t_max=t_max, t_min=t_min)

tests/core/atoms/test_phonons.py

@@ -7,26 +7,35 @@
 
 import numpy as np
 from ase.build import bulk
+from ase.constraints import FixAtoms
 from numpy.testing import assert_almost_equal, assert_array_equal
 
 from quacc.atoms.phonons import get_phonopy
 
 
 def test_get_phonopy():
     atoms = bulk("Cu")
-    phonopy = get_phonopy(atoms)
+    phonopy, _ = get_phonopy(atoms)
     assert_array_equal(phonopy.supercell_matrix, [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
 
-    phonopy = get_phonopy(atoms, min_lengths=5)
+    phonopy, _ = get_phonopy(atoms, min_lengths=5)
     assert_array_equal(phonopy.supercell_matrix, [[2, 0, 0], [0, 2, 0], [0, 0, 2]])
 
-    phonopy = get_phonopy(atoms, min_lengths=[5, 10, 5])
+    phonopy, _ = get_phonopy(atoms, min_lengths=[5, 10, 5])
     assert_array_equal(phonopy.supercell_matrix, [[2, 0, 0], [0, 4, 0], [0, 0, 2]])
 
-    phonopy = get_phonopy(atoms, displacement=1)
+    phonopy, _ = get_phonopy(atoms, displacement=1)
     assert_almost_equal(
         phonopy.displacements, [[0, 0.0, np.sqrt(2) / 2, np.sqrt(2) / 2]]
     )
 
-    phonopy = get_phonopy(atoms, symprec=1e-8)
+    phonopy, _ = get_phonopy(atoms, symprec=1e-8)
     assert phonopy.symmetry.tolerance == 1e-8
+
+    atoms = bulk("Cu") * (2, 2, 2)
+
+    atoms.set_constraint(FixAtoms(indices=[0, 1, 2, 3]))
+
+    phonopy, fixed_atoms = get_phonopy(atoms, min_lengths=5)
+
+    assert len(fixed_atoms) == 4

tests/core/recipes/emt_recipes/test_emt_phonons.py

@@ -7,7 +7,7 @@
 pytest.importorskip("phonopy")
 pytest.importorskip("seekpath")
 
-from ase.build import bulk
+from ase.build import bulk, molecule
 
 from quacc.recipes.emt.phonons import phonon_flow
 
@@ -70,4 +70,22 @@ def test_phonon_flow_v4(tmp_path, monkeypatch):
     assert output["results"]["thermal_properties"]["temperatures"][-1] == 1000
     assert output["results"]["force_constants"].shape == (8, 8, 3, 3)
     assert "mesh_properties" in output["results"]
-    assert output["atoms"] != atoms
+    assert output["atoms"] == atoms
+
+
+def test_phonon_flow_fixed(tmp_path, monkeypatch):
+    monkeypatch.chdir(tmp_path)
+    atoms = molecule("H2", vacuum=20.0)
+    output = phonon_flow(atoms, run_relax=True, min_lengths=5.0)
+
+    atoms1 = molecule("H2", vacuum=20.0)
+    atoms2 = molecule("H2", vacuum=20.0)
+
+    atoms2.positions += [10, 10, 10]
+
+    output_fixed = phonon_flow(atoms1, non_displaced_atoms=atoms2, min_lengths=5.0)
+
+    # Should be very close but not exactly the same
+    assert output["results"]["mesh_properties"]["frequencies"] == pytest.approx(
+        output_fixed["results"]["mesh_properties"]["frequencies"], rel=0.0, abs=1e-5
+    )