Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
wants to merge 30 commits into
base: main
Choose a base branch
from

Conversation

tomdemeyere
Copy link
Contributor

Summary of Changes

Should be much better

Checklist

  • I have read the "Guidelines" section 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.

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. This is solely for the maintainers.
  • When your code is ready for review, ping one of the active maintainers.

@buildbot-princeton
Copy link
Collaborator

Can one of the admins verify this patch?

Copy link
Member

@Andrew-S-Rosen Andrew-S-Rosen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@tomdemeyere: Thank you for taking the time to revisit this. It is certainly a lot cleaner and is very easy to follow now! I have some minor comments below, none of which are anything crazy. Overall, very clever modification to the original approach!

I do just have one question: the additional_atoms is an interesting idea. I'm curious to hear more about how you intend to use this in practice (or see people using it in practice).

In surface catalysis, I am most familiar with the opposite scenario to some degree. To calculate thermochemical corrections, people will invoke weak coupling between adsorbate and the metal surface. They'll then vibrate the adsorbate atoms and, perhaps, a few of the adjacent atoms on the surface while fixing the rest. They'll then use the harmonic approximation (e.g. with ASE utilities) to calculate the corrections. This is implicitly assuming that when you compare the energies of two different states along the reaction coordinate, the vibrational modes of the slab are effectively unchanged.

Here, it looks like it's almost the opposite perspective. I guess that's simply because here you are trying to calculate the phonon spectra of the material? Is that right? If doing so, wouldn't fixing a few adsorbate atoms not really make a big difference compared to the size of the slab?

src/quacc/atoms/phonons.py Outdated Show resolved Hide resolved
supercell_matrix
The supercell matrix to use. If specified, it will override any
value specified by `min_lengths`.
Returns
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing a line break before Returns.

Comment on lines +113 to +114
The supercell matrix to use. If specified, it will override any
value specified by `min_lengths`.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is no min_lengths in this function. I think we can just ignore that part. Also, it must be specified since it is a positional argument.

Comment on lines +55 to +58
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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this is relevant anymore.

Comment on lines 28 to 29
if has_deps:
from phonopy import Phonopy
pass
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems to be unnecessary.

src/quacc/recipes/common/phonons.py Outdated Show resolved Hide resolved
Comment on lines +127 to +132
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,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I assume you just removed the type hints for brevity's sake since they are basically implied elsewhere? I was also thinking about this since it is an internal function. Although, if we at some point enabled static type checking, having the type hints would ensure there isn't an error. I doubt we will be able to get to the point where we enable static type checking though... certainly not anytime soon.

@@ -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,
additional_atoms: Atoms | None = None,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note to self: we'll want to add this parameter to the other recipes calling the phonon_subflow (tblite and MLP I think?).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this test still needs updating. get_phonopy() doesn't return a tuple anymore, and we aren't dealing with FixAtoms. That's the reason for the failure in the CI.

Comment on lines +99 to +104
def get_atoms_supercell_by_phonopy(
atoms: Atoms,
supercell_matrix: tuple[
tuple[int, int, int], tuple[int, int, int], tuple[int, int, int]
],
) -> Atoms:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be ideal to have a simple unit test for this specific function.

@@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Based on the test failure in CI, it looks like the original atoms object is getting mutated somewhere (since it has an EMT calculator attached)? That doesn't seem ideal... 😅

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this was related to #2230.

@tomdemeyere
Copy link
Contributor Author

Overall additional atoms can be seen as the same to "fixed_atoms". I renamed it to put an emphasis on the fact that these atoms are merely spectators of the Phonons calculation, only having an implicit influence through the forces. When you think about this has nothing to do with "fixing".

If you want to calculate an adsorbate on a slab that way you would just send the adsorbate part to regular atoms and the slab as "additional_atoms".

This could be renamed as "fixed_atoms" as well, I am trying to explore the semantic here

@Andrew-S-Rosen
Copy link
Member

Andrew-S-Rosen commented May 10, 2024

Overall additional atoms can be seen as the same to "fixed_atoms". I renamed it to put an emphasis on the fact that these atoms are merely spectators of the Phonons calculation, only having an implicit influence through the forces. When you think about this has nothing to do with "fixing".

If you want to calculate an adsorbate on a slab that way you would just send the adsorbate part to regular atoms and the slab as "additional_atoms".

This could be renamed as "fixed_atoms" as well, I am trying to explore the semantic here

To clarify, I understand that part. My question was about understanding why and when this is relevant, given that it is the opposite from what is usually invoked in the literature.

Edit: Ah, you are saying the slab would go to additional atoms. Sure, that makes sense. I think that conflicts with what you wrote in one of the docstrings. I think additional_atoms makes sense.

additional_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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was my point of confusion. I had thought you were implying that additional atoms would be for the adsorbates (I now see that's not the case). In reality, it would be for the surface of course.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I completely see where this confusion can get from, this is why I am not entirely sure about the current name.

Although this forces the users to stop and think before using this important approximation

EDIT: yes, the doc does not help as well

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One downside of fixed_atoms is that it potentially implies that they are sites within atoms (the first positional argument) that are fixed, whereas in reality this is a separate Atoms object. That's one reason why I like additional_atoms since it makes it clear that this is a separate Atoms object, although one would understand that if they read the type hint. Of course, additional_atoms doesn't inherently specify anything about them being fixed until you read the docstring.

I think either could be okay, personally.

@Andrew-S-Rosen
Copy link
Member

Historically, in such scenarios, I have used ASE's Vibrations module (specifying the atoms to vibrate with indices) and then passed that information into ASE's HarmonicThermo class. However, using phonopy for this seems like it might be even better because it would be possible to run the displacement calculations in parallel rather than sequentially.

At some point, I was going to add HarmonicThermo to quacc.runners.thermochemistry to support this, but the phonopy approach is great.

@Andrew-S-Rosen
Copy link
Member

What if we rename the first positional argument to unfixed_atoms and then additional_atoms as fixed_atoms? Not sure if that helps or hurts...

@tomdemeyere
Copy link
Contributor Author

tomdemeyere commented May 10, 2024

Claude 3 Opus is suggesting (😅):

  1. primary_atoms (instead of atoms) and auxiliary_atoms
  2. atoms and fixed_atoms
  3. displaced_atoms (instead of atoms) and nondisplaced_atoms

In this order of "preference", I might actually lean towards 2, (let's just be pragmatic?)

EDIT: To be fair, 3 looks like your last suggestion but without the double negative of "unfixed", this might be the best option to clearly distinguish between what is being moved and the rest, "nondisplaced_atoms" might just be "fixed_atoms" for simplicity?

@Andrew-S-Rosen
Copy link
Member

Ooh I like 3! It's better than my suggestion because it doesn't accidentally imply a FixAtoms constraint. But I'm okay with 2 if you have a strong opinion about it.

@tomdemeyere
Copy link
Contributor Author

That's a good argument let's go with that, I just added an extra underscore

@@ -92,15 +98,16 @@ def phonon_flow(
parameters=job_params,
decorators=job_decorators,
)
if run_relax:
if run_relax and not non_displaced_atoms:
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a little bit problematic, another option is to send both displaced and non_displaced and optimise them and separate them again...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. I mean, I guess we can ask if we really even need a relaxation step beforehand? Presumably one could just call relax_job() before calling the flow and it would be exactly the same (although without a tight force tolerance).

If there's a desire to keep this though, then yeah I think the only route would be to relax displaced_atoms+non_displaced_atoms as a single combined_atoms and then pass in combined_atoms[:len(displaced_atoms)] etc. to phonon_subflow. It does, admittedly, complicate things a bit because it seems weird to relax something that is called "non-displaced", but... 🤷

@Andrew-S-Rosen Andrew-S-Rosen mentioned this pull request May 10, 2024
3 tasks
@Andrew-S-Rosen Andrew-S-Rosen changed the title Demeyere phonopy fixed atoms v2 Add support for having non-displaced atoms in Phonopy routines May 10, 2024
@Andrew-S-Rosen Andrew-S-Rosen added the enhancement New feature or request label May 13, 2024
Andrew-S-Rosen added a commit that referenced this pull request Jun 12, 2024
Incorporate some of the changes from #2110.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Development

Successfully merging this pull request may close these issues.

None yet

3 participants