Phonopy: fixing atoms

[tomdemeyere]

Summary of Changes

I like the idea that you had, removing the parameter fixed_atoms and deal with constraints instead, some problems though: users might want to optimise these atoms before hand during the relax_job, or sometimes not, how to deal with the constraints then? Also, this might be seen as a "hidden" feature?

This MR also really removes the second get_phonopy which somehow evaded both of our attention during the last PR 😂

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]

Can one of the admins verify this patch?

[Andrew-S-Rosen]

I like the idea that you had, removing the parameter fixed_atoms and deal with constraints instead, some problems though: users might want to optimise these atoms before hand during the relax_job, or sometimes not, how to deal with the constraints then? Also, this might be seen as a "hidden" feature?

This is worth thinking about. Let me think out loud here.

In Scenario A, the user wants to relax all the atoms and then, afterwards, do a phonon calculation with a subset of atoms fixed. In this scenario, they should just run a relax_job() first and then pass the output to phonon_flow() with run_relax=False.

from quacc.recipes.emt.core import relax_job
from quacc.recipes.emt.phonons import phonon_flow

output = relax_job(atoms)
new_atoms = output["atoms"]
new_atoms.set_constraint(...)
output2 = phonon_flow(new_atoms, run_relax=False)

In Scenario B, the user wants to have a subset of atoms fixed during both the relax and the phonon flow. In this case, the only need to call phonon_flow() directly (with the run_relax=True default set) and their input Atoms objects with constraints.

So, in both scenarios, this is currently possible. While I agree that adding a dedicated keyword argument is more convenient for the user, the counterpoint to this is that it adds more lines to maintain even though it is not strictly necessary. Right now, virtually every recipe in quacc supports fixing atoms in any way you'd like simply by passing in an Atoms object with constraints pre-applied. We didn't have to do anything to achieve that; it comes "for free" by simply basing things around ASE. Admittedly, this could be communicated more clearly to end-users and should likely (at the very least) have a dedicated call out in the introductory section of the documentation: #2074.

This MR also really removes the second get_phonopy which somehow evaded both of our attention during the last PR 😂

Oops, yeah that was a mistake. I also like the addition you have made regarding symmetrization of the force constants.

[tomdemeyere]

the counterpoint to this is that it adds more lines to maintain even though it is not strictly necessary

To use ASE's constraints the only way that comes to my mind is to do:

fixed = []

for constr in atoms.constraints:
    if isinstance(constr, FixAtoms):
        fixed.extend(constr.get_indices())

Then remove the fixed_atoms parameter, but the rest of the code would not change.

[Andrew-S-Rosen]

@tomdemeyere: Thanks. My point was more about the need for these changes at all if one can pass in an Atoms object with constraints already applied. For instance, why is the following needed among all the other changes with the indexing?

    if len(fixed) > 0:
        fixed = get_phonopy_structure(AseAtomsAdaptor.get_structure(fixed))
        fixed = phonopy.structure.cells.get_supercell(fixed, supercell_matrix)
        fixed = AseAtomsAdaptor.get_atoms(get_pmg_structure(fixed))

I don't yet understand what the issues are (in main) with providing an ASE Atoms object with constraints (like mentioned in my prior comment). Maybe put another way: if we removed the fixed_atoms keyword argument proposed here, why would the other added lines still be needed? Please note that I have not used phonopy in practice so may need additional context here if there's motivation beyond what you put in your PR summary. I am going based on what you have written.

The other changes, specifically regarding getting rid of the second get_phonopy() call and adding symmetrization of the force constants, are quite logical.

[tomdemeyere]

Apologies for the lack of context. What we have to do to effectively fix some atoms when using phonopy is the following:

1. Construct the supercell but only sending the partial system (non-fixed) to phonopy, this way phonopy thinks the system being calculated is the non-fixed one. (lines 80 to 93 in quacc.atoms.phonons.get_phonopy)
2. Re-construct the full system manually by doing the same supercell transformation for the rest of the system (the fixed atoms) outside of phonopy (lines 95 to 99 in quacc.atoms.phonons.get_phonopy)
3. The generated displacements will only be for the non-fixed atoms due to point 1, but we still send the full systems to the calculator by adding all the fixed system to each displacement (lines 101 to 104 in quacc.recipes.common.phonons)
4. Phonopy expects forces with shape (len(disp), len(non_fixed_atoms), 3) so we manually extract these forces from the full calculations (lines 115 to 118 in quacc.recipes.common.phonons)
5. Post-processing as usual...
6. As a consequence, the symmetry of the force constant is broken, which we manually try to enforce later.

ASE's constraints unfortunately are not of help here. This whole trick has to be done manually.

This whole process drastically reduces the number of displacement required, for an adsorbate on a slab for example you might end up calculating the displacement of the adsorbate only. As you might guess this is a huge approximation because you assume that modes from the fixed system do not interact with the ones of the non-fixed system.

On the previous MR I did a test which showed that for water on Pt(111) the trick does not change frequencies much as the slab and waters are loosely coupled.

[Andrew-S-Rosen]

@tomdemeyere: Thank you very much for your helpful reply. I had (incorrectly) assumed based on your PR description that this was mainly centered around ease-of-use. However, I now understand that it is actually correcting the behavior for fixed atoms and that simply supplying the Atoms object with FixAtoms constraints is not sufficient. I am familiar with the typical assumption about decoupling of vibrational modes between surface and adsorbate but have normally used ASE's Vibrations module and ASE's thermo modules for this instead of phonon calculations via phonopy. If I recall correctly, ASE's Vibrations module is smart enough not to displace and vibrate atoms that are fixed.

Anyway, I have left several comments below. Apologies in advance about some of them being a bit pedantic, but it will help increase the clarity for future me.

[Andrew-S-Rosen]

I think I need to return to this later to fully digest what's going on, but I'll trust you on it for the time being.

[tomdemeyere]

I don't really like the fact that get_phonopy now returns a tuple but I don't see any other non-convoluted way to do that :/

[Andrew-S-Rosen]

I don't really like the fact that get_phonopy now returns a tuple but I don't see any other non-convoluted way to do that :/

I agree that it feels a bit wrong since the Atoms object isn't needed in most cases and the function is primarily about getting a Phonopy object. If you end up thinking of a better alternative, by all means go ahead. If not though, I won't let it block the PR.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.06%. Comparing base (a095785) to head (5bcd7e4).
Report is 1 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #2073   +/-   ##
=======================================
  Coverage   99.06%   99.06%           
=======================================
  Files          81       81           
  Lines        3311     3325   +14     
=======================================
+ Hits         3280     3294   +14     
  Misses         31       31           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[tomdemeyere]

Last minute thought:

"Fixing" atoms during phonons calculations is not a real thing and many people would simply call it blasphemy (to be fair
it seems that many people just use it so that they can say "we did the free energy calculations" in their paper). I think it might be wise to add an extra warning with the logger, especially considering that many Phonopy users will not expect this feature to exist, what do you think?

[Andrew-S-Rosen]

@tomdemeyere: Thanks for the suggestion. I have been gravitating back towards having the fixed_indices keyword argument for this reason; the user really has to directly request it and know what they are getting themselves into. I have this taken care of in a custom branch I am working on right now (not public yet).

My main hesitations right now are that quacc.atoms.phonons.get_phonopy() is just super messy/complicated. There is a ton of back-and-forth converting happening, and it is not trivial to follow along. I am trying to make some headway in cleaning that up, although I don't know if I'll be successful... in any case, I'll post here later this afternoon whether I'm successful or not.

[tomdemeyere]

"I have been gravitating back towards having the fix_indices keyword argument for this reason; the user really has to directly request it and know what they are getting themselves into."

I think this might be better. When you suggested to use ASE's constraints I got hooked as it is definitely more consistent and "appealing", but after thinking about it a little bit more I would vote for keeping fix_indices for the reason you invoked.

About Phonopy in Quacc: It might be more interesting to make simpler Phonopy functions instead of whole flows, aiming for increased modularity, even at the cost of increased complexity. Phonopy is a very complex package and in the current state available flows can only scratch the package's surface due to the lack of possible parameters customization. This might also feel better in a workflow sense(from what I understand): if you would be better off resolving in the middle of a flow something is off (you know better, just shooting an idea here)

Users seeking something simpler might be better off being redirected to ASE's phonon/vibrations workflows (which are still powerfuls but arguably less complex).

[Andrew-S-Rosen]

@tomdemeyere: Apologies again that this hasn't been as trivial to merge as I initially believed.

I have opened up PR #2108 to demonstrate some areas for improvement. Namely, I want to draw your attention to quacc.atoms.phonons. Note that by sticking with the pymatgen Structure object, we can avoid back-and-forths between different atoms representations. Additionally, throughout this PR there were several areas where clearer variable names were needed; it got very hectic.

There are still some problems with PR #2108, however:

• The tests are failing, so I clearly made a mistake somewhere. But I don't have time to dig into this further.
• The pre-existing tests from Phonopy: fixing atoms #2073 are pretty easy to have pass even if we made a mistake. I am not 100% confident I haven't made an error. Additional confirmation that everything is running appropriately is needed, and, ideally, more rigorous tests should be made.
• The quacc.recipes.common.phonons.phonon_subflow is getting increasingly complicated with this fix atoms behavior. While not a dealbreaker, it is a bit unfortunate.
• The fact that quacc.atoms.phonons.get_phonopy returns a tuple even when fixed_indices=None is a bit annoying, although it is less severe of a problem in Rosen phonopy fixed atoms (attempt 2) #2108 than here in Phonopy: fixing atoms #2073.

Feel free to copy or pick things up from there. Unfortunately, I don't feel comfortable merging #2073 until this gets streamlined, especially for what amounts to a crude approximation in many cases (although one that is quite common).

---

About Phonopy in Quacc: It might be more interesting to make simpler Phonopy functions instead of whole flows, aiming for increased modularity, even at the cost of increased complexity. Phonopy is a very complex package and in the current state available flows can only scratch the package's surface due to the lack of possible parameters customization. This might also feel better in a workflow sense(from what I understand): if you would be better off resolving in the middle of a flow something is off (you know better, just shooting an idea here)

Agreed. There is a constant, delicate balance between recipes doing too much and too little. At the end of the day, quacc has been intentionally designed to give advanced users all the tools they need to build their own custom workflows. Quacc isn't necessarily meant to codify every option available to users, and my vision is that end users can make their own scripts as they see fit. At the same time, when flexibility is available and makes sense to support, it is good to have. I don't think we need to try to overload the workflow to access every phonopy option available.

[tomdemeyere]

I think I just got an idea that might avoid quite a lot of the current mess to fix atoms, Also avoiding get_phonopy returning a tuple.

Let me have a look tomorrow.

[Andrew-S-Rosen]

Closing in favor of #2110.