From e99d41aa0ebbff8e3892bea140357b0cc97ca6b7 Mon Sep 17 00:00:00 2001 From: Tom Demeyere <115232841+tomdemeyere@users.noreply.github.com> Date: Fri, 11 Oct 2024 16:40:27 +0100 Subject: [PATCH] Add support for having non-displaced atoms in Phonopy routines (#2110) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ## 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 --- src/quacc/atoms/phonons.py | 30 +++++++ src/quacc/recipes/common/phonons.py | 78 +++++++++++++++---- src/quacc/recipes/emt/phonons.py | 7 ++ src/quacc/runners/phonons.py | 8 ++ tests/core/atoms/test_phonons.py | 12 ++- .../recipes/emt_recipes/test_emt_phonons.py | 32 +++++++- 6 files changed, 148 insertions(+), 19 deletions(-) diff --git a/src/quacc/atoms/phonons.py b/src/quacc/atoms/phonons.py index a81e19b7b9..438a84909d 100644 --- a/src/quacc/atoms/phonons.py +++ b/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 + ) + ) diff --git a/src/quacc/recipes/common/phonons.py b/src/quacc/recipes/common/phonons.py index 14d272a9d5..c8899d2ee9 100644 --- a/src/quacc/recipes/common/phonons.py +++ b/src/quacc/recipes/common/phonons.py @@ -5,12 +5,19 @@ 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")) @@ -18,8 +25,6 @@ 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,7 +52,7 @@ def phonon_subflow( additional_fields: dict[str, Any] | None = None, ) -> PhononSchema: """ - Calculate phonon properties. + Calculate phonon properties using the Phonopy package. Parameters ---------- @@ -54,6 +60,11 @@ def phonon_subflow( 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 diff --git a/src/quacc/recipes/emt/phonons.py b/src/quacc/recipes/emt/phonons.py index 2a111bb1f9..35de90d1f9 100644 --- a/src/quacc/recipes/emt/phonons.py +++ b/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, diff --git a/src/quacc/runners/phonons.py b/src/quacc/runners/phonons.py index 552416f5be..da76808251 100644 --- a/src/quacc/runners/phonons.py +++ b/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) diff --git a/tests/core/atoms/test_phonons.py b/tests/core/atoms/test_phonons.py index 1db5a3a7da..c7a8076cbc 100644 --- a/tests/core/atoms/test_phonons.py +++ b/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) diff --git a/tests/core/recipes/emt_recipes/test_emt_phonons.py b/tests/core/recipes/emt_recipes/test_emt_phonons.py index 5e4f5e9a8c..4c26e8f341 100644 --- a/tests/core/recipes/emt_recipes/test_emt_phonons.py +++ b/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 + )