Add support for having non-displaced atoms in Phonopy routines (#2110)

## Summary of Changes

Should be much better

### Checklist

- [ ] I have read the ["Guidelines"
section](https://quantum-accelerators.github.io/quacc/dev/contributing.html#guidelines)
of the contributing guide. Don't lie! 😉
- [ ] My PR is on a custom branch and is _not_ named `main`.
- [ ] I have added relevant, comprehensive [unit
tests](https://quantum-accelerators.github.io/quacc/dev/contributing.html#unit-tests).

### Notes

- Your PR will likely not be merged without proper and thorough tests.
- If you are an external contributor, you will see a comment from
[@buildbot-princeton](https://github.com/buildbot-princeton). This is
solely for the maintainers.
- When your code is ready for review, ping one of the [active
maintainers](https://quantum-accelerators.github.io/quacc/about/contributors.html#active-maintainers).

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Andrew S. Rosen <[email protected]>

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
@@ -94,3 +95,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.
+
+    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,21 +5,26 @@
 from importlib.util import find_spec
 from typing import TYPE_CHECKING
 
+import numpy as np
+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
+from quacc.utils.dicts import recursive_dict_merge
 
 has_phonopy = bool(find_spec("phonopy"))
 has_seekpath = bool(find_spec("seekpath"))
 
 if TYPE_CHECKING:
     from typing import Any
 
-    from ase.atoms import Atoms
-
     from quacc import Job
     from quacc.types import PhononSchema
 
@@ -33,6 +38,7 @@
 def phonon_subflow(
     atoms: Atoms,
     force_job: Job,
+    fixed_atom_indices: list[int] | None = None,
     symprec: float = 1e-4,
     min_lengths: float | tuple[float, float, float] | None = 20.0,
     supercell_matrix: (
@@ -46,14 +52,19 @@ def phonon_subflow(
     additional_fields: dict[str, Any] | None = None,
 ) -> PhononSchema:
     """
-    Calculate phonon properties.
+    Calculate phonon properties using the Phonopy package.
 
     Parameters
     ----------
     atoms
         Atoms object with calculator attached.
     force_job
         The static job to calculate the forces.
+    fixed_atom_indices
+        Indices of 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
@@ -79,6 +90,33 @@ def phonon_subflow(
     PhononSchema
         Dictionary of results from [quacc.schemas.phonons.summarize_phonopy][]
     """
+    mask_to_fix = np.zeros(len(atoms), dtype=bool)
+
+    if fixed_atom_indices:
+        mask_to_fix[fixed_atom_indices] = True
+
+    displaced_atoms, non_displaced_atoms = atoms[~mask_to_fix], atoms[mask_to_fix]
+
+    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_supercell = get_atoms_supercell_by_phonopy(
+            non_displaced_atoms, phonopy.supercell_matrix
+        )
+    else:
+        non_displaced_atoms_supercell = Atoms()
+
+    supercells = [
+        phonopy_atoms_to_ase_atoms(s) + non_displaced_atoms_supercell
+        for s in phonopy.supercells_with_displacements
+    ]
 
     @subflow
     def _get_forces_subflow(supercells: list[Atoms]) -> list[dict]:
@@ -97,9 +135,17 @@ def _thermo_job(
         additional_fields: dict[str, Any] | None,
     ) -> 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,17 +156,15 @@ 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
-    ]
+    if non_displaced_atoms:
+        additional_fields = recursive_dict_merge(
+            additional_fields,
+            {
+                "displaced_atoms": displaced_atoms,
+                "non_displaced_atoms": non_displaced_atoms,
+            },
+        )
+
     force_job_results = _get_forces_subflow(supercells)
     return _thermo_job(
         atoms, phonopy, force_job_results, t_step, t_min, t_max, additional_fields

src/quacc/recipes/emt/phonons.py

@@ -27,6 +27,7 @@ def phonon_flow(
         tuple[tuple[int, int, int], tuple[int, int, int], tuple[int, int, int]] | None
     ) = None,
     displacement: float = 0.01,
+    fixed_atom_indices: list[int] | None = None,
     t_step: float = 10,
     t_min: float = 0,
     t_max: float = 1000,
@@ -63,6 +64,11 @@ def phonon_flow(
         value specified by `min_lengths`.
     displacement
         Atomic displacement (A).
+    fixed_atom_indices
+        Indices of 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
@@ -96,6 +102,7 @@ def phonon_flow(
         symprec=symprec,
         min_lengths=min_lengths,
         supercell_matrix=supercell_matrix,
+        fixed_atom_indices=fixed_atom_indices,
         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

@@ -9,7 +9,7 @@
 from ase.build import bulk
 from numpy.testing import assert_almost_equal, assert_array_equal
 
-from quacc.atoms.phonons import get_phonopy
+from quacc.atoms.phonons import get_atoms_supercell_by_phonopy, get_phonopy
 
 
 def test_get_phonopy():
@@ -30,3 +30,13 @@ def test_get_phonopy():
 
     phonopy = get_phonopy(atoms, symprec=1e-8)
     assert phonopy.symmetry.tolerance == 1e-8
+
+
+def test_get_supercell_by_phonopy():
+    atoms = bulk("Cu")
+    get_atoms_supercell_by_phonopy(atoms, np.eye(3) * 2)
+
+    cell = [[-1, 1, 1], [1, -1, 1], [1, 1, -1]]
+
+    supercell = get_atoms_supercell_by_phonopy(atoms, cell)
+    assert_almost_equal(np.diag(np.diag(supercell.cell)), supercell.cell)

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
 
@@ -58,3 +58,33 @@ def test_phonon_flow_v3(tmp_path, monkeypatch):
     assert output["results"]["force_constants"].shape == (8, 8, 3, 3)
     assert "mesh_properties" in output["results"]
     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, min_lengths=5.0)
+
+    atoms1 = molecule("H2", vacuum=20.0)
+    atoms2 = molecule("H2", vacuum=20.0)
+
+    atoms2.positions += [10, 10, 10]
+
+    atoms = atoms1 + atoms2
+
+    output_fixed = phonon_flow(atoms, fixed_atom_indices=[2, 3], min_lengths=1.0)
+    # Should be very close but not exactly the same, also check the size is correct
+    assert output["results"]["mesh_properties"]["frequencies"] == pytest.approx(
+        output_fixed["results"]["mesh_properties"]["frequencies"], rel=0.0, abs=1e-5
+    )
+
+    assert len(output_fixed["displaced_atoms"]) == 2
+    assert len(output_fixed["non_displaced_atoms"]) == 2
+
+    output_fixed_wrong = phonon_flow(atoms, fixed_atom_indices=[0, 2], min_lengths=1.0)
+
+    assert output_fixed_wrong["results"]["mesh_properties"][
+        "frequencies"
+    ] != pytest.approx(
+        output_fixed["results"]["mesh_properties"]["frequencies"], rel=0.0, abs=1e-5
+    )