diff --git a/src/quacc/atoms/phonons.py b/src/quacc/atoms/phonons.py index a634badfe1..eb90e8b292 100644 --- a/src/quacc/atoms/phonons.py +++ b/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,20 +60,38 @@ 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, @@ -80,7 +99,16 @@ def get_phonopy( **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: diff --git a/src/quacc/recipes/common/phonons.py b/src/quacc/recipes/common/phonons.py index 9da3b6920c..dd1d3668b9 100644 --- a/src/quacc/recipes/common/phonons.py +++ b/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,7 +124,7 @@ def _thermo_job( additional_fields=additional_fields, ) - phonopy = get_phonopy( + phonopy, fixed_atoms = get_phonopy( atoms, min_lengths=min_lengths, supercell_matrix=supercell_matrix, @@ -118,8 +132,12 @@ def _thermo_job( 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( diff --git a/src/quacc/runners/phonons.py b/src/quacc/runners/phonons.py index 3200a54c03..5128bb050d 100644 --- a/src/quacc/runners/phonons.py +++ b/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) diff --git a/tests/core/atoms/test_phonons.py b/tests/core/atoms/test_phonons.py index 1db5a3a7da..8bd544710c 100644 --- a/tests/core/atoms/test_phonons.py +++ b/tests/core/atoms/test_phonons.py @@ -7,6 +7,7 @@ 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 @@ -14,19 +15,27 @@ 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 diff --git a/tests/core/recipes/emt_recipes/test_emt_phonons.py b/tests/core/recipes/emt_recipes/test_emt_phonons.py index 46584eb629..7e097e753b 100644 --- a/tests/core/recipes/emt_recipes/test_emt_phonons.py +++ b/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 + )