Add support for having non-displaced atoms in Phonopy routines

src/quacc/atoms/phonons.py

@@ -14,6 +14,7 @@
 
 if has_phonopy:
     from phonopy import Phonopy
+    from phonopy.structure.cells import get_supercell
 
 if TYPE_CHECKING:
     from ase.atoms import Atoms
@@ -25,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,
@@ -91,3 +92,32 @@ 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(Structure.from_ase_atoms(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 PhonopyRunner
 from quacc.schemas.phonons import summarize_phonopy
 
@@ -18,21 +23,17 @@
 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_phonopy:
-        from phonopy import Phonopy
-
 
 @subflow
 @requires(has_phonopy, "Phonopy must be installed. Run `pip install quacc[phonons]`")
 @requires(has_seekpath, "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 +49,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 +91,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 +121,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 = PhonopyRunner().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 +150,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

@@ -19,13 +19,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,
@@ -50,7 +51,7 @@ def phonon_flow(
 
     Parameters
     ----------
-    atoms
+    displaced_atoms
         Atoms object
     symprec
         Precision for symmetry detection.
@@ -61,6 +62,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
@@ -90,15 +96,16 @@ def phonon_flow(
         param_swaps=job_params,
         decorators=job_decorators,
     )
-    if run_relax:
-        atoms = relax_job_(atoms)["atoms"]
+    if run_relax and not non_displaced_atoms:
+        displaced_atoms = relax_job_(displaced_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

@@ -36,6 +36,7 @@ def run_phonopy(
         self,
         phonon: Phonopy,
         forces: NDArray,
+        symmetrize: bool = False,
         t_step: float = 10,
         t_min: float = 0,
         t_max: float = 1000,
@@ -50,6 +51,8 @@ def run_phonopy(
             Phonopy object
         forces
             Forces on the atoms
+        symmetrize
+            Whether to symmetrize the force constants
         t_step
             Temperature step
         t_min
@@ -66,6 +69,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
+    )