Phonopy: fixing atoms

src/quacc/atoms/phonons.py

@@ -5,6 +5,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+from ase.constraints import FixAtoms
 from monty.dev import requires
 from pymatgen.io.ase import AseAtomsAdaptor
 from pymatgen.io.phonopy import get_phonopy_structure, get_pmg_structure
@@ -35,7 +36,7 @@ def get_phonopy(
     symprec: float = 1e-5,
     displacement: float = 0.01,
     phonopy_kwargs: dict | None = None,
-) -> Phonopy:
+) -> tuple[Phonopy, Atoms]:
     """
     Convert an ASE atoms object to a phonopy object with displacements generated.
 
@@ -59,28 +60,55 @@ def get_phonopy(
     -------
     Phonopy
         Phonopy object
+    Atoms
+        The fixed atoms as an ASE atoms object.
     """
     phonopy_kwargs = phonopy_kwargs or {}
 
-    structure = AseAtomsAdaptor().get_structure(atoms)
+    fixed_indices = np.array(
+        [
+            constr.get_indices()
+            for constr in atoms.constraints
+            if isinstance(constr, FixAtoms)
+        ]
+    )
+
+    fixed_indices.flatten()
+    is_fixed_atoms = np.array([i in fixed_indices for i in range(len(atoms))])
+
+    structure = AseAtomsAdaptor.get_structure(atoms)
     structure = SpacegroupAnalyzer(
         structure, symprec=symprec
     ).get_symmetrized_structure()
+    atoms = structure.to_ase_atoms()
+
+    fixed_atoms, non_fixed_atoms = atoms[is_fixed_atoms], atoms[~is_fixed_atoms]
+
+    non_fixed_atoms = AseAtomsAdaptor.get_structure(non_fixed_atoms)
 
     if supercell_matrix is None and min_lengths is not None:
         supercell_matrix = np.diag(
             np.round(np.ceil(min_lengths / atoms.cell.lengths()))
         )
 
-    phonopy_atoms = get_phonopy_structure(structure)
+    phonopy_atoms = get_phonopy_structure(non_fixed_atoms)
     phonon = phonopy.Phonopy(
         phonopy_atoms,
         symprec=symprec,
         supercell_matrix=supercell_matrix,
         **phonopy_kwargs,
     )
     phonon.generate_displacements(distance=displacement)
-    return phonon
+
+    if fixed_atoms:
+        fixed_atoms = phonopy_atoms_to_ase_atoms(
+            phonopy.structure.cells.get_supercell(
+                get_phonopy_structure(AseAtomsAdaptor.get_structure(fixed_atoms)),
+                supercell_matrix,
+            )
+        )
+
+    return phonon, fixed_atoms
 
 
 def phonopy_atoms_to_ase_atoms(phonpy_atoms: PhonopyAtoms) -> Atoms:

src/quacc/recipes/common/phonons.py

@@ -5,6 +5,7 @@
 from importlib.util import find_spec
 from typing import TYPE_CHECKING
 
+import numpy as np
 from monty.dev import requires
 
 from quacc import job, subflow
@@ -48,6 +49,11 @@ 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
@@ -97,9 +103,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"][~fixed_indices, :]
+            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=fixed_indices.any(),
+            t_step=t_step,
+            t_min=t_min,
+            t_max=t_max,
         )
 
         return summarize_phonopy(
@@ -110,16 +124,20 @@ def _thermo_job(
             additional_fields=additional_fields,
         )
 
-    phonopy = get_phonopy(
+    phonopy, fixed_atoms = get_phonopy(
         atoms,
         min_lengths=min_lengths,
         supercell_matrix=supercell_matrix,
         symprec=symprec,
         displacement=displacement,
         phonopy_kwargs=phonopy_kwargs,
     )
+    fixed_indices = np.array(
+        [False] * len(phonopy.supercell) + [True] * len(fixed_atoms)
+    )
     supercells = [
-        phonopy_atoms_to_ase_atoms(s) for s in phonopy.supercells_with_displacements
+        phonopy_atoms_to_ase_atoms(s) + fixed_atoms
+        for s in phonopy.supercells_with_displacements
     ]
     force_job_results = _get_forces_subflow(supercells)
     return _thermo_job(

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,8 @@
 pytest.importorskip("phonopy")
 pytest.importorskip("seekpath")
 
-from ase.build import bulk
+from ase.build import bulk, molecule
+from ase.constraints import FixAtoms
 
 from quacc.recipes.emt.phonons import phonon_flow
 
@@ -71,3 +72,33 @@ def test_phonon_flow_v4(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, run_relax=True, min_lengths=5.0)
+
+    atoms1 = molecule("H2", vacuum=20.0)
+    atoms2 = molecule("H2", vacuum=20.0)
+
+    atoms2.positions += [10, 10, 10]
+
+    full_system = atoms1 + atoms2
+    full_system.set_constraint(FixAtoms(indices=[0, 1]))
+
+    output_fixed = phonon_flow(full_system, run_relax=True, 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
+    )
+
+    full_system = atoms1 + atoms2
+    full_system.set_constraint(FixAtoms(indices=[0, 2]))
+
+    output_wrong = phonon_flow(full_system, run_relax=True, min_lengths=5.0)
+
+    assert output["results"]["mesh_properties"]["frequencies"] != pytest.approx(
+        output_wrong["results"]["mesh_properties"]["frequencies"], rel=0.0, abs=1e-2
+    )