[WIP] magnetic spacegroup analyzer

[mt-huebsch]

Summary

Added MagneticSpacegroupAnalyzer that given a magnetic structure

• finds all magnetic spacegroup operations
• determines the magnetic spacegroup
• can return the magnetic pointgroup
• can compute all magnetic domains

Added ApplyMagSymmOpTransformation to standard transformations

• apply MagSymmOp to a magnetic structure

Additional dependencies introduced

None

To Do

• add a test

I was not sure how/where to add a test. This snipped shows how to run the main functionality of the code, i.e. finding the magnetic spacegroup:
path_to_s = "0.13_Ca3Co2-xMnxO6.mcif"
s = mg.Structure.from_file( path_to_s)
ma = MagneticSpacegroupAnalyzer( s )
print( "OG info: ", ma.get_og() )

0.13_Ca3Co2-xMnxO6.mcif can be downloaded from MAGNDATA
http://webbdcrista1.ehu.es/magndata/index.php?index=0.13

The OG info for this material is ('R3c', '161.1.1300')

To obtain the magnetic domains, run
domain_dict = ma.get_domains()
print( "Magnetic domains: ", len(domain_dict), [mop.as_xyzt_string() for mop in domain_dict.keys()] )

This retruns
Magnetic domains: 4 ['x, y, z, +1', 'x, y, z, -1', '-x, -y, -z, +1', '-x, -y, -z, -1']

I'd be happy to provide more extensive tests, if you point me in the right direction.

• enable pip to get pymatgen/symmetry/mspg2mpg_map.txt

I noticed that this file does not get installed when installing pymatgen from my repository in a new enviroment. Maybe I should have chosen a different extension? Or does this need to be configured?

• code review

This is my first contribution, so I would appreciate any comment.

Comment

In the past I have submitted a contribution, where the search for the magnetic spacegroup was dependent on tools of the Bilbao Crystallographic Server. This has now been replaced by a search independent of any external sources. The code is robust enough to handle uncertainties of a typical ab initio run and can return all magnetic domains. Thus, if two ab initio results correspond to different magnetic domains, this feature will help identify the correspondence.

Looking forwards to your feedback!

Cheers,
Marie-Therese

[mkhorton]

Thank you @mt-huebsch! This looks like an incredibly useful contribution, I would be happy to see this merged.

I will take some time to review the PR carefully to provide comments and a response to your question on testing.

I noticed that this file does not get installed when installing pymatgen from my repository in a new enviroment. Maybe I should have chosen a different extension? Or does this need to be configured?

By default, only .json, .csv and .yaml files are packaged as specified by this line in the setup.py. We can either choose to include .txt files too (either all .txt files, or a specific file), or we can change the file to be one of these other formats instead.

[mt-huebsch]

Dear @mkhorton, thank you for your quick response!

Regarding the .txt file, I will change the format. No problem.
I will wait for some more comments from you and then make the changes.

Thank you for showing interest in my contribution.

[mkhorton]

I've made a partial code review with individual comments above. Most of my comments so far are not major, they're mainly regarding code conventions, rather than correctness. Please let me know if there are any questions!

[mkhorton]

Even for private methods, a quick docstring is helpful.

[mkhorton]

cartesian option is not documented in docstring

[mkhorton]

Capitalization should be avoided inside variable names, and should be spaced by underscore, e.g. this should be self._spacegroup_analyzer. This follows the official "PEP8" recommendations.

[mkhorton]

Perhaps worth discussing -- I think this functionality might already be partially present in StructureMatcher or perhaps CollinearMagneticStructureAnalyzer.match(), or maybe this method has a more narrow application than I understand? Where can struc1 and struc2 come from here?

[mkhorton]

An empty set evaluates as False in Python, so this statement could be if not res -- this is no more or less correct, but this latter notation is more commonly used. The variable res could be called diff to make its purpose more clear.

[mkhorton]

Following PEP8, and to be more explicit, this could be called _search_for_magnetic_spacegroup

[mkhorton]

exit() might make some users sad since it would quit their Python session, instead an Exception should be raised. For example, raise VallueError("Error in the definitions of the unit cell.").

For approximation comparisons, np.all_close() is a nice short-hand for the np.any(np.abs(...) > ...) construction.

[mkhorton]

We can make this a module-level variable so that the file doesn't have to be re-loaded and parsed every time. I can make a further practical suggestion here for how to structure the file, so that it's not necessary to iterate over lines individually.

[mkhorton]

I've attached a possible alternative as a yaml file (will have to be unzipped, since GitHub doesn't allow me to upload .yaml directly).

Then, at the top of this file, we could write:

from monty.serialization import loadfn

MAGSPACEGROUP_TO_MAGPOINTGROUP = loadfn(os.path.join(os.path.dirname(__file__), "magspacegroup_to_magpointgroup.yaml"))

This will be loaded as a dictionary keyed by the OG no., so in this get_magnetic_pointgroup function it can just return MAGSPACEGROUP_TO_MAGPOINTGROUP[glob_og_no]["pointgroup_label"], ... for example.

(monty.serialization.loadfn and monty.serialization.dumpfn are convenience functions for loading and dumping Python objects, pymatgen objects etc. to JSON or YAML)

magspacegroup_to_magpointgroup.yaml.zip

[guycmoore]

This is fantastic. Thank you @mt-huebsch!

I have a couple of questions. The major roadblock that I encountered while trying to develop a magnetic symmetry analyzer was for cases in which the "setting" of the msg operations is non-standard. From doing a search on the MAGNDATA server, you can see that the majority of magnetic structures (822 out of the 1217 structures), the msg operations are defined in a non-standard setting (defined by a similarity transform from the standard lattice or vector space).

While I get the same very nice result for Ca3Co2-xMnxO6, I get an error:
Error in the list of symmetry operations.
if I try a structure in which the msg is in a non-standard setting, such as FePO4, which has an msg P212121 (#19.25) within the non-standard setting (a,b,c;0,1/2,3/4).

I'm curious if I'm doing something wrong? Here is how I got the error:

>>> filename_a = "0.13_Ca3Co2-xMnxO6.mcif"
>>> struct_a = Structure.from_file( filename_a )
>>> ma_a = MagneticSpacegroupAnalyzer( struct_a )
>>> print( "OG info: ", ma_a.get_og() )
OG info:  ('R3c', '161.1.1300')
>>> 
>>> filename_b = "0.17_FePO4.mcif"
>>> struct_b = Structure.from_file( filename_b )
>>> ma_b = MagneticSpacegroupAnalyzer( struct_b )
Error in the list of symmetry operations.

It wasn't trivial to me how to find the similarity transform for a msg. I'm curious if that is what's happening here. I'll dive deeper into your code to see if you account for this. Did you encounter a similar issue? It's possible to check this setting transformation with the symmetry finder in the ISOTROPY suite: https://stokes.byu.edu/iso/findsym.php.

I believe that there are ways of identifying if two sg representations are "similar" using matrix normal forms, which appears to be the basis for many of the symmetry finding algorithms developed by Stokes et. al. and in spglib. I would be happy to talk about this in more depth if you also explored this.

If you found a work-around for finding the transformation to another setting, I am very curious to learn what your approach was. I don't see where you account for this in the code but I'll look more closely. If this is an issue, maybe this feature can be implemented in the future.

Thank you again!

[mt-huebsch]

Dear @guymoore13,

Thank you for running tests on the code!!

Indeed, you have discovered a bug. I need to transform to the standardized conventional cell. As far as I understand the same ambiguity is already present in the spacegroup operations, so spglib can handle it. I am already working on the solution and can perform all necessary transformations, but at the moment I haven't figured out how "standardized" is defined exactly. Maybe I am wrong but spglib and Stokes et al seem to have a different definition of standardized.

Please let me know in case you have any thoughts on that.

[guycmoore]

Dear @mt-huebsch, thank you for your response!

You are correct that spglib accounts for this, although I'm also still trying to identify how they go about addressing the setting/basis issue. Due to the ambiguity in the "standard setting" in which to compare operations, I was trying a different approach - because I encountered this question as well.

For this reason, I was trying to find a way of comparing "reduced representations" of the msg operations - which should be independent from the settings of each representation (and therefore the similarity transform to the standard setting). I haven't made much progress on this front, but I will let you know if I do - it may be nice to compare notes at some point. I was having issues with the "uniqueness" of the diagonalization scheme that I was using.

Your approach to find the standard setting may be more robust. I look forward to hearing how it goes!

[maxmarkov]

Great work and a very important contribution! I want to report just a couple of minor things I have noticed so far

1. In _paramag_symmetry_operations of MagneticSpacegroupAnalyzer class:
   rotation = self._spacegroupAnalyzer._space_group_data['rotations']
translation = self._spacegroupAnalyzer._space_group_data['translations']

Your rotations and translations come directly from spglib. However, it often returns small translation vectors that are duplicated i,e. [0, 0, 0], [1e-17,1e-17,1e-17] etc You can reproduce it on 1.0.1_Ag2CrO2.mcif from MAGNDATA. As a result, the number of symmetry operations is larger than it should be (40 symmetries instead of 4 in case of Ag2CrO2). A similar issue is described in _get_symmetry of SpacegroupAnalyzer. One can avoid it by using

rotation = self._spacegroupAnalyzer._get_symmetry()[0]
translation = self._spacegroupAnalyzer._get_symmetry()[1]

2. You got the following line twice:
   if np.any(np.abs(self._structure._lattice.matrix - self._spacegroupAnalyzer._cell[0]) > 1e-12): exit("Error in the definitions of the unit cell.")

First, in __init__ and then in _paramag_symmetry_operations.

[mkhorton]

@maxmarkov thanks for the testing. Could you clarify how:

rotation = self._spacegroupAnalyzer._get_symmetry()[0]
translation = self._spacegroupAnalyzer._get_symmetry()[1]

avoids the duplicated SymmOp issue? Related to this, we use np.allclose() to define whether two SymmOps are equal but the tolerance is currently very tight, do you think this should be relaxed to make it easier to filter out duplicates?

[mt-huebsch]

Dear @maxmarkov! Thank you for testing. Honestly, I haven't had a look at q>0 magnetic structures, but that is just a side remark. I will add it to my to do list.

Regarding _spacegroupAnalyzer._get_symmetry(), it is always good to make the code more robust. But rather than hardcoding a rounding, maybe it should be implemented using pos_tolerance with a good fall back default value. What is pymatgen's policy on that, @mkhorton?

Bdw, previously I used _spacegroupAnalyzer._get_symmetry_operations() at that point, but this function always returns the conventional symmetries. This is unexpected behavior to me, if the input structure is primitive. That is why I wrote _paramag_symmetry_operations() in the first place. Although it is probably not good if the SpacegroupAnalyzer and the MagneticSpacegroupAnalyzer behave differently. Harmonizing the behavior should maybe be added to the to do list, too.

@guymoore13, regarding the setting issue, I have a small update:

Every MSPG symbol is derived from a non-magnetic SPG symbol. That is why fortunately the same alternative settings apply. But the standard setting for MSPGs is not the one closest to the standard setting of the parent SPG [1].
Furthermore, Hall symbols [2] represent an explicit-origin spacegroup notation. There is a method to decode Hall symbols into matrix representations described in ITC volume B. The resulting matrices are precomputed and stored in spglib [3].
Spglib determines matrix representations for SPG operations and then compares to that dataset. Both tasks, determining the symmetries and comparing representations, are quite involved and currently not available individually.
It should be possible to overload the initializer of spglib to accept a set of symmetry operations as an input instead of a structure. Then, the transformation matrix to a chosen setting (by selecting the Hall number) can be returned. I will contact Togo-sensei and look into spglib's repo to check how difficult this will be.

[1] https://www.annualreviews.org/doi/10.1146/annurev-matsci-070214-021008
[2] https://doi.org/10.1107/S0567739481001228
[3] https://arxiv.org/pdf/1808.01590.pdf

I will try to carve out more time to get this functionality ready for distribution. Thank you all for your feedback and patience.

Best regards!

[guycmoore]

Thank you for the detailed response @mt-huebsch. Yes, I think we are saying more or less the same thing. You mention the Hall symbol, and my approach has been to store the matrix normal forms required to identify the Hall symbols (on page 7 of the spglib paper that you shared). They use the Smith canonical form. I saw this a bit later - I was using a different matrix diagonalization. However, I recently discovered that the Smith canonical form has some very nice properties, especially for integer matrices. They use SageMath to compute the Smith normal forms. I have been experimenting with other Python packages to compute normal forms.

Take care, Guy