Add the create_magnetic_configuration function

[mbercx]

Fixes #639

In order to pass on the magnetic configuration between calculations, the
structure used for the follow-up calculation needs to have the right kinds
to be able to properly assign the magnetization to each site.

Here we add a calculation function that, based on the structure and
magnetic moments, returns a new StructureData with the required
magnetic kinds, as well as a Dict with the corresponding magnetic
moments for each kind.

To create the new list of kinds, the algorithm loops over all the elements
in the structure and makes a list of the sites with that element and their
corresponding magnetic moment. Next, it splits this list in three lists:

• Zero magnetic moments: Any site that has an absolute magnetic moment
  lower than ztol
• Positive magnetic moments
• Negative magnetic moments

The algorithm then sorts the positive and negative lists from large to small
absolute value, and loops over each list. New magnetic kinds will be created
when the absolute difference between the magnetic moment of the current
kind and the site exceeds atol.

The positive and negative magnetic moments are handled separately to
avoid assigning two sites with opposite signs in their magnetic moment to
the same kind and make sure that each kind has the correct magnetic
moment, i.e. the largest magnetic moment in absolute value of the sites
corresponding to that kind.

[greschd]

I have a high-level question that might be slightly off-topic, so apologies in advance:

The underlying issue here is that QE requires each magnetic moment to be specified as a separate kind. However, the AiiDA StructureData is a container for only the atomic types + positions, and doesn't contain magnetic moments. I feel this is partially tying the generic StructureData to QE implementation details.
So I wonder if we would get a more consistent interface by pushing this down to the PwCalculation by adding an input that roughly corresponds to magnetic_moment_per_site. As a result, the same StructureData could be used regardless of magnetization, and the splitting into different kind names with suffixes would happen there. The drawback would be that all magnetism-related inputs would need special treatment.

I have to admit I didn't think this through (so it's probably a bad idea), I'm mostly interested in your thoughts on how to split the interface between generic and QE-specific parts.

[mbercx]

I have to admit I didn't think this through (so it's probably a bad idea), I'm mostly interested in your thoughts on how to split the interface between generic and QE-specific parts.

Actually, I think your suggestion makes sense. While working on this, I was thinking it would be great if the StructureData contains the magnetic moments. In this case we could do exactly as you suggest: push down the creation of the new kinds using the algorithm of the create_magnetic_allotrope calcfunction down to the PwCalculation.

Originally the create_magnetic_allotrope calcfunction had an ArrayData input so the output of the calculation of the magnetic moments would be connected to the new StructureData in the provenance, which I think was one @sphuber's motivations for using a calcfunction. However, as we changed how the function works this got lost in favour of simplifying the inputs.

So, at this point I'm leaning in favour of your suggestion to move this logic into the PwCalculation, and provide an magnetic_moment_per_site input. I suppose we'd also have to provide inputs for the tolerances. Maybe there is something I'm missing though, let's see what @sphuber thinks. I'm also going to ping @yakutovicha and @espenfl for comments on how magnetism is dealt with in the aiida-cp2k and aiida-vasp plugins.

One thing I just thought of: what if the user provides the magnetic_moment_per_site input but doesn't set nspin to 2? Do we set it in the prepare_for_submission method? What if magnetic_moment_per_site is provided and nspin is also set to 1? Do we override the NSPIN value, or raise an error?

---

To comment some more on splitting QE-specific and generic parts: Isn't the reliance of StructureData on kinds also a consequence of historical reasons (i.e. @giovannipizzi using Quantum ESPRESSO? 🙃)? Should we consider removing kinds from StructureData at some point, in favor of a list of sites with occupations - possibly magnetic moments as well - and have the calculation job of the plugin deal with how to convert this information into the required input?

[greschd]

input but doesn't set nspin to 2

just a quick side note: non-collinear magnetism (noncolin=.TRUE., which automatically sets nspin=4) can also have initial magnetization.

[mbercx]

Another note: If we push the create_magnetic_allotrope logic down to the PwCalculation, it also becomes easier to pass the magnetic configuration between different calculations in a work chain. Currently this is not done for e.g. the PwRelaxWorkChain, the main reason being that this might require new magnetic kinds (I think). We could also call the create_magnetic_allotrope inside the work chain, but this would possibly create a lot of new StructureDatas that just have different kinds.

[greschd]

And yet another note: QE has a hard-coded limit on the number of kinds that is allowed (ntypx, defined in parameters.F90). The default is 10, to use more kinds QE needs to be re-compiled.

[espenfl]

Actually, I think your suggestion makes sense. While working on this, I was thinking it would be great if the StructureData contains the magnetic moments. In this case we could do exactly as you suggest: push down the creation of the new kinds using the algorithm of the create_magnetic_allotrope calcfunction down to the PwCalculation.

Just to put it into context. We decided that e.g. the degrees of freedom when doing relaxations should not go into the StructureData (or a more general container based on that) as it is a result of a calculation which imposes some restriction.

For the magnetics it goes deeper. As you know, the symmetry also changes when introducing magnetic moments. So then the question maybe becomes more, what is StructureData now and what should it be in the future? A following question is, should it be connected to symmetry? If so, it somehow needs to connect to the magnetics. If not, it does not have to be directly connected and one could think if the magnetics as part of the calculation, which results in a given StructureData (kind of like we do with the degrees of freedom for restrictive relaxation etc.). Right now I believe it is not really explicitly connected to symmetry, but implicitly due to the added k-points (which would obey symmetry). It is just a container for positions of some kind, in a container (cell) etc. Maybe I am wrong or have not inspected the code in too much detail.

[giovannipizzi]

Just trying to answer the historical part.
Admittedly, an inspiration from how things are done in QE is there. However, other codes that I had checked at the time also had the concept of "kinds" or "species". And, more generally, you can think of 'kinds' as a "compression" feature, where you don't need to repeat certain properties for each site. Imagine a huge supercell with 1000s of atoms, you would need repeat, for each atom, the mass, if it's a virtual alloy, etc.
Instead, by keeping this information on the 'species/kinds', sites are just defined by a list of strings (the kind_name) + a Nx3 array of floats (that is currently in the attributes, but could be moved to a numpy array on disk as it's done now for Trajectories.

An additional note: it's always possible to go in a well-defined and unique way from a description with kinds + sites, to a description with only sites (just assign all kind properties to each site with that kind name). The reverse, unfortunately, is not always possible (it's similar to the issue you have, of deciding a threshold to determine when to consider two sites being of the same species or not.

So, from a design point of view, and also for simplicity, I would keep the current design with kinds, unless there are obvious problems (that I don't see).

Now, whatever you do, you will still have some issues in some cases, so I don't think there is a perfect solution:

• in QE as @greschd mentioned, the number of kinds is limited at compile time, AND also the maximum length of their name (3 characters! so you cannot have e.g. Fe10). So even if possible to define one kind per site, then in the QE input one has to decide how to convert this to an actual input. HOWEVER, note that this is not a flaw with the StructureData design, but with QE's input design: if you need to have 10 different Fe kinds, it's tricky to do it in QE independent on whether you have kinds or not in StructureData
• if you remove Kinds from StructureData, you have the same issue of "merging" sites in few kinds for ANY QE calculation (imagine, e.g. that there is some numerical noise in definition of masses and you end up with multiple kinds for the same atom type)

A different question is whether to add information on magnetisation or not in StructureData. In the original design, I decided to keep StructureData simple and not do it, with the ideas that:

• most codes can enable/disable magnetism, so in a sense I see it more as a code input
• physically, the structure itself is enough to determine the material, and then "chemistry" will decide the magnetic behaviour, so I wanted to keep the StructureData definition closer to what you get, say, from XRD.

Now: when I made this decision, I wasn't 100% sure it was the perfect, and I am still not (even if I'm still 85-90% positive that's OK :-) ). So if you are all convinced it's better to move magnetic information in StructureData I won't oppose to it (but you need to balance the advantage - as mentioned above, it's going to require a change in all plugins that need to get the information from there and convert it to something the code understands, with the risk that some codes might not be able to convert in a well-definite way. One way, if you really want to do it (maybe the only reason I see now is the discussion above on the symmetry of the system when sites have different magnetisation, but also there one could think to pass an optional parameter with the magnetisation to any symmetry-detection function, together with the structure), is to have a subclass that extends StructureData in MagneticStructureData or something like this, so that it does not break existing code.
In this case, to answer @mbercx's question: if the user specifies incompatible inputs, you raise a ValidationError when creating the inputs (a bit like when a user defined a 'blocked' input like nsites). The only difference here is that you cannot just add nspin to the blocked keywords, as the user still needs to specify it, if it passes a standard StructureData. So you'll need to perform the validation in the prepare_for_submission.

To conclude: considering the amount of work that might be needed to do any change, involving all plugins, I would be personally hesitant to do changes to StructureData unless the benefits outweigh the issues. If it's really important to do a change, I would open an issue just to discuss that, and make sure the majority of active plugin developers can discuss if it's more of an improvement or a of burden for them.

[ltalirz]

Just my 2 cents: I've also recently struggled with the fact that StructureData isn't able to transport information on magnetic moments, which has forced me to use a non-AiiDA object (ASE Atoms) to pass this information around instead in aiida-lsmo.

In view of improving interoperability between ASE and AiiDA going forward, I would personally highly welcome if the information that can be stored in an ASE Atom could be stored in StructureData as well (i.e. besides the atomic mass that would be initial charge, magmom, and momentum).
Especially, if this could be achieved in a backwards-compatible way.

[sphuber]

In view of improving interoperability between ASE and AiiDA going forward, I would personally highly welcome if the information that can be stored in an ASE Atom could be stored in StructureData as well (i.e. besides the atomic mass that would be initial charge, magmom, and momentum).

Fair point and one that we hear more often from users, especially if they have worked with ase and/or pymatgen.

Not sure how much you have been involved in the discussions surrounding this topic, but we should of course not forget that the AiiDA data structures are fundamentally quite different from those in ASE, as in they are immutable. Not saying that they are therefore incompatible, but we always find out the hard way that the immutability of AiiDA's datastructures makes working with these objects way more complex than with the ASE analogue which can cause frustration. For example, if we were to incorporate the magmom for each site in the StructureData, if you want to change those moments, let's say for another calculation, you should create a new StructureData, and if you want to keep the provenance of that structure, you have to go through a calcfunction. People often ask why it has to be so complex in AiiDA, but that is simply the "cost" of provenance I think. Anyway, since we have recently discussed this a lot in the context of the common workflows (there also with regards to how to limit the movement of certain atoms in certain cartesian directions), justed wanted to highlight this here as well.

[espenfl]

* 

A different question is whether to add information on magnetisation or not in StructureData. In the original design, I decided to keep StructureData simple and not do it, with the ideas that:

* most codes can enable/disable magnetism, so in a sense I see it more as a code input

* physically, the structure itself is enough to determine the material, and then "chemistry" will decide the magnetic behaviour, so I wanted to keep the StructureData definition closer to what you get, say, from XRD.

Well this depends. It can be argued that symmetry is a much more fundamental concept than applied chemistry/physics. So if we base something on chemistry and physical principles, it should definitively obey symmetry. If not, there should be very good arguments for not doing this (and it should be clear to the users why that is). In my view we should have a representation of points in space, within a domain. This is somewhat what StructureData is now. Then, there is a more general object that attaches this to e.g. symmetry concepts, which will depend on external factors and the context of usage of the represenetation of points in space.

I think, as has been discussed before, it is probably a good time to start to relax the connection of the datastructure and AiiDA core and give the existing ones a bit more context. E.g. put the StructureData, KpointsData etc. into a plugin/domain etc. Doing so will make it much more easy to define what the respective data structures should represent and contain.

[giovannipizzi]

Hi Espen, thanks.

Let me try just to summarise here some discussions we had recently in the team about moving out these classes.

While I think most people agree that moving StructureData and similar in a different package is a good idea in principle, we've been discussing this multiple times in the various AiiDA meetings and it's not clear if the benefits outweigh the amount of work to make it happen. It's not just moving the code out, it's really making sure that everything works fine also in an intermediate situation, without cyclic dependencies between packages etc.

E.g. now in 1.6 one could move it out and make AiiDA depend on the new package for now, so as to make sure all code continue working when you do from aiida.orm import StructureData. But then if, say, in 2.0 we want to remove completely StructureData from AiiDA-core, you would essentially have to "invert" the dependencies and make the domain-specific package depend on AiiDA, and make sure that people's code (and plugins) do not break. It's not even clear to us if this is possible in a robust way. Another complication is that the entry point string and the node_type string in the DB will have to change and this will always be an "abrupt" change even if migrations are applied in either 1.6 or 2.0, and therefore will break queries of users.

In addition, moving out the plugins does not solve the real problem. These discussions will still have to take place before deciding possible extensions (and among the same people) with the additional disadvantage that if we change something, now we could still make a migration to "fix" things, but once it's in an external package, we can't touch the data anymore.

For all these reasons, it's still not clear to us if it's worth investing so much energy in moving them out at this stage, considering there is other work to do in the code that might have a much larger impact.

This is just to explain why this has not happened yet, and what are the technical details that make it complex (that might not be obvious when one first thinks about it).

[espenfl]

Thanks @giovannipizzi. Maybe it would make sense to continue this discussion on the discussion boards of AiiDA, but I will just quickly reply to your comment.

Yes, I fully understand this is a complex undertaking, but the longer we wait the more complex it is going to be. Also it only pertains to the domain specific data types. It is nice to see in the docs that the domain specific types are now covered in their own section. In my view there are three paths into the future: (i) leave it as is, (ii) leave as is, but start the development of other data types that are put into dedicated plugins and recommend users to use those (iii) try to approach the task of separating out the material science data types now.

With respect to this not really addressing the true problem, in some ways it will as it will give additional freedom to set more specific bounds what the data types are supposed to describe/represent. Now, it is not so easy to define this as it would maybe be difficult to define what StructureData is going to be going forward. I would assume if we want to keep it (and continue to use it) it should be very broad, which is suboptimal for portability and maintenance. Again, I am coming back to the symmetry question for this particular definition (there are other aspects as well). We would have to decide if this type of StructureData should contain symmetry. If so, it needs to respect magnetism and it should either connect closely to some magnetic data set, or fully embed it.

Also, with respect to the definitions of the data types, I personally think we should focus on this being as portable as possible given that the ontology is still pretty floating.

[yakutovicha]

I would like to also express my view on the matter.

• most codes can enable/disable magnetism, so in a sense I see it more as a code input

Personally, I don't think the ability to enable/disable a certain feature is a to-be-considered argument when we are talking about a class. I think one of the points of working with objects is to assign to it the properties, which by nature belong to it. To me, magnetization is the feature of a structure, so I would vote to add it there.

Additionally, a code can decide not to use magnetization that is coming from the structure.

• physically, the structure itself is enough to determine the material, and then "chemistry" will decide the magnetic behaviour, so I wanted to keep the StructureData definition closer to what you get, say, from XRD.

I think that magnetization is also a factor that influences the structure of a material (remember that the structures of ferro or antiferro iron are different).

To conclude: considering the amount of work that might be needed to do any change, involving all plugins, I would be personally hesitant to do changes to StructureData unless the benefits outweigh the issues.

I don't think we need to request plugin developers to change anything if we add the magnetic moment attribute and we don't remove kinds. I believe the atomic kinds can and should stay there to cover the corner cases when certain elements of the same type should be considered differently.

Finally, from my point of view it would be nice to add create_allotrope() method able to group atoms according to the selected properties (not only magnetization). Something like this:

s = StructureData(...)

allotrope = s.create_allotrope(properties = ['mass', 'magnetic_moment'])

# The content of the allotrope.
{
    kinds = [ 'Fe1', 'Fe2'],
    properties: {
        'Fe1': {'magmom': 1},
        'Fe2': {'magmom': -1},
    },
}

[mbercx]

@sphuber after seeing a copy of this code in the aiida-common-workflows project, I think it's time to merge this pull request. This feature is already being used in production in several projects, also some of the high-throughput projects (e.g. MC3D).

Although the discussion here is important in the context of aiidateam/team-compass#21, it's clear that this switch to a new StructureData is still not around the corner. And so until we have a StructureData type that is able to deal with magnetic moments this feature would be good to have in a proper release, and can be easily deprecated in the future.

If you agree, I'll add some documentation and we can finally merge this one.

[sphuber]

@mbercx Agreed. If you can update the branch making sure tests pass and add some minimal docs, I will approve. Thanks!

[mbercx]

Tests failing seem to be related to the same issue that caused the nightly build to fail. It seems to be a pydantic related error:

E   pydantic.errors.ConfigError: duplicate validator function "aiida_quantumespresso.common.hubbard.HubbardParameters.check_manifolds"; if this is intended, set `allow_reuse=True`

But that seems strange... The nightly build also succeeded when I ran it for Python 3.10, see https://github.com/aiidateam/aiida-quantumespresso/actions/runs/5052086649

[mbercx]

I think the issue was related to pydantic/pydantic#5821, i.e. a bug in pydantic version 1.10.X in connection with typing_extensions 4.6.0. Not 100% sure about the details, but should be fixed in 1.10.8 which was released yesterday. Rerunning the tests installs that version and so they pass.

Will open a PR to restrict our pydantic dependency.

[sphuber]

Thanks @mbercx . It has been a long time coming this one 😄 Let's get this things merged

[mbercx]

Hold it, still working on the documentation. 😅

[mbercx]

Blocked since built on #946

[mbercx]

@sphuber I've still made a few changes/additions as I was working on the documentation for this:

1. Changed the name of the calculation function to create_magnetic_configuration. Although perhaps allotrope is somewhat more precise, it just isn't used very often in this context I think.
2. Added a get_magnetic_configuration tool for the PwCalculation. This one felt quite natural as I was working on the documentation, and significantly simplifies getting the proper structure + magnetic configuration dict from a previous pw.x run.

Things (potentially) left to do:

• Add documentation for the create_magnetic_configuration calculation function to the How-to or Topics sections. Here I'm wondering where we should add this. I don't want the "calculations" how-to sections to be too numerous.
• Same for the PwCalculation tools. Here I'm thinking perhaps having a reference at the bottom of the "Topics" page for each CalcJob that has tools may be preferable.