How to retain site_properties when converting from primitive to conventional cell?

Hello everybody!

I am trying to make some slabs in a FireWorks workflow. Initially I am getting the structures from the MaterialsProject, which has the nice feature that the magnetic moments are included as site_properties of the structure. E.g. if I search for the NiO structure with the lowest total energy I get the rocksalt prototype with the correct AFM-2 magnetic ordering and 4 atoms in the unit cell:

from pymatgen import MPRester
struct = MPRester().get_structure_by_material_id('mp-19009')
struct.sites
Out[3]: 
[PeriodicSite: Ni (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
 PeriodicSite: Ni (0.0022, -2.1073, -2.1072) [0.5000, 0.5000, 0.5000],
 PeriodicSite: O (0.0011, -2.1083, 0.0011) [0.2500, 0.7500, 0.2500],
 PeriodicSite: O (0.0032, -2.1062, -4.2156) [0.7500, 0.2500, 0.7500]]
struct.site_properties
Out[4]: {'magmom': [-1.742, 1.742, 0.0, 0.0]}

If I now want to generate a slab using the SlabGenerator of pymatgen.core.surface, I should use a conventional unit cell, so that the miller indices are what one would expect. I can do this by using the get_conventional_standard_structure method of the SpacegroupAnalyzer class, but the site_properties are not transfered!

from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
conv_struct = SpacegroupAnalyzer(struct).get_conventional_standard_structure()
conv_struct.sites
Out[6]: 
[PeriodicSite: Ni (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
 PeriodicSite: Ni (0.0000, 2.1083, 2.1083) [0.0000, 0.5000, 0.5000],
 PeriodicSite: Ni (2.1083, 0.0000, 2.1083) [0.5000, 0.0000, 0.5000],
 PeriodicSite: Ni (2.1083, 2.1083, 0.0000) [0.5000, 0.5000, 0.0000],
 PeriodicSite: O (2.1083, 0.0000, 0.0000) [0.5000, 0.0000, 0.0000],
 PeriodicSite: O (2.1083, 2.1083, 2.1083) [0.5000, 0.5000, 0.5000],
 PeriodicSite: O (0.0000, 0.0000, 2.1083) [0.0000, 0.0000, 0.5000],
 PeriodicSite: O (0.0000, 2.1083, 0.0000) [0.0000, 0.5000, 0.0000]]
conv_struct.site_properties
Out[7]: {}

Is there a way the magnetic moments stored in the magmom key of the site_properties could be transferred to the conventional cell? Is there another way to get the conventional cell, or to transfer the miller indices to the corresponding ones for the primitive structure? Since the SlabGenerator does transfer the site_properties correctly, it would be much easier to create slabs with complex magentic structures if this would work somehow.

BTW, If I use conventional_unit_cell=True to get the initial structure, the site_properties are empty as well:

struct = MPRester().get_structure_by_material_id('mp-19009', conventional_unit_cell=True)
struct.site_properties
Out[5]: {}

Hi @mwo, I’ve noticed this too. This behavior ultimately comes from the way we construct the conventional cell:

from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
sga = SpacegroupAnalyzer(struct)
print(sga.get_refined_structure())

The site properties are lost when we use spglib to do the symmetry refinement. We could add a patch there to add site properties back in, if you’d like to take a go at this yourself pull requests are very welcome. Without spending too much time on it, it looks like a call to add the properties back when the new Structure is constructed should be sufficient, perhaps using StructureMatcher to providing the mapping of sites between input and conventional structures if that information is lost.

Best,

Matt

Thinking about this more, it’s not trivial or guaranteed that such a mapping exists. Going from primitive to conventional should always be fine, but these functions can take any arbitrary input structure (including supercells), where the conventional crystallographic structure might not have enough sites to retain the magnetic moment information (e.g. imagine a magnetic cell with a k vector != (0, 0, 0)).

Thanks for the reply @mkhorton!

I agree that this will be difficult to implement as a general case. However, your idea of using the StructureMatcher to find a mapping between the sites of the conventional and the input structure is a good way to proceed. In the end, at least for my application, this could be done on the Fireworks level, where I can be sure that I go from a primitive to an conventional cell. Once I have a mapping it should be simple to transfer the site_properties.

Playing around with the structure matcher however got me nowhere. I cannot seem to get a mapping between the conventional and the primitive cell. I guess this has to do with point 3. in the explanation of the StructureMatcher algorithm?

If the number of sites do not match, return False

Staying with the example for NiO:

from pymatgen import MPRester
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.analysis.structure_matcher import *
struct = MPRester().get_structure_by_material_id('mp-19009')
conv_struct = SpacegroupAnalyzer(struct).get_conventional_standard_structure()
matcher = StructureMatcher(primitive_cell=False)
matcher.fit(struct1=conv_struct, struct2=struct)
Out[26]: False
print(matcher.get_mapping(conv_struct, struct))
None

If I turn on the conversion to primitive cells, the structure is found to be matched, but I cannot get a mapping, which would be useless anyhow:

matcher_prim = StructureMatcher(primitive_cell=True)
matcher_prim.fit(struct1=conv_struct, struct2=struct)
Out[31]: True
matcher_prim.get_mapping(conv_struct, struct)
Traceback (most recent call last):

  File "<ipython-input-33-751cb7f7d59d>", line 1, in <module>
    matcher_prim.get_mapping(conv_struct, struct)

  File "/home/mwo/FireWorks/atomate_env/lib/python3.7/site-packages/pymatgen/analysis/structure_matcher.py", line 1141, in get_mapping
    raise ValueError("cannot compute mapping with primitive cell "

ValueError: cannot compute mapping with primitive cell option

Using a different comparator and switching some possible flags on does not help either:

matcher = StructureMatcher(primitive_cell=False, allow_subset=True,
                           scale=False, comparator=ElementComparator())
print(matcher.get_mapping(conv_struct, struct))
None

Any further help with this would be greatly appreciated. Is there even a way to get a mapping between structures with different number of sites? Will I have to dive into spglib?
Thanks again,
Michael

The way I would do this, is label the primitive sites with site-properties and then generate the “supercell” that will get you to a structure with the same number of sites. Then you can use structure matcher to get your mapping and then follow back via the labels which are preserved when making a supercell. You’ll have to decide how you want to deal with site-properties then.

Another approach could be to query both the primitive and conventional structures and then feed the conventional structure to the function get_conventional_to_primitive_transformation_matrix. Then you could invert that matrix and apply it to the primitive structure (with SymmOp) that has the decorated magnetic properties. However, that’s just a rough idea, not sure how robust that would be.

Dear @shyamd, thanks for the idea! However, I am still not able to get this to work…
Even if I have the same number of sites in the conventional cell and the supercell (8, I tried 1x1x2, 1x2x1, and 2x1x1 supercells) the structure matcher returns a match only if I use the primitive_cell=True option, and than there is no mapping.

Here is the full example:

from pymatgen import MPRester
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer as sga
from pymatgen.analysis.structure_matcher import *

struct = MPRester().get_structure_by_material_id('mp-19009')
conv_struct = sga(struct).get_conventional_standard_structure()
doubled_struct = struct.copy()
doubled_struct.make_supercell([1, 1, 2])

matcher = StructureMatcher(primitive_cell=False, allow_subset=True,
                           scale=False, comparator=ElementComparator())
print(matcher.fit(struct1=conv_struct, struct2=doubled_struct))
print(matcher.get_mapping(conv_struct, doubled_struct))

prim_matcher = StructureMatcher(primitive_cell=True, allow_subset=True,
                           scale=False, comparator=ElementComparator())
print(prim_matcher.fit(struct1=conv_struct, struct2=doubled_struct))
print(prim_matcher.get_mapping(conv_struct, doubled_struct))

Dear @peterschindler, thanks for the input!

I am not sure if this should work even in theory, since the conventional cell has 8 sites, while the primitive cell has only 4 (in the AFM-2 magentic structure that is, in theory the crystal structure alone would also allow for a 2 atom cell). As I understand SymmOp, which is not very well) it should not add sites, but only rotate and translate, right?

I tried it out however, but there are only 4 atoms in the transformed primitive structure and the cell looks nothing like the conventional unit cell should. However, there is a good possibility that I screwed up with the affine transformation, so I would appreciate if you could check my attempt! Thanks!

Here is the full example:

import numpy as np
from pymatgen import MPRester
from pymatgen.io.vasp.inputs import Poscar
from pymatgen.core.operations import SymmOp
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer as sga

struct = MPRester().get_structure_by_material_id('mp-19009')
struct_conv = MPRester().get_structure_by_material_id('mp-19009',
                                                conventional_unit_cell=True)
TM = sga(struct_conv).get_conventional_to_primitive_transformation_matrix()
# invert the TM and add no translation to make it an affine 4x4 transformation matrix.
Affine_invM = np.hstack([np.vstack([np.linalg.inv(TM),[0,0,0]]),
                       [[0],[0],[0],[1]]])
prim_struct = struct.copy()
prim_struct.apply_operation(SymmOp(Affine_invM))
print(prim_struct)
print(struct_conv)

However, going back to my initial objective of retaining site_properties in the slab I am creating from a bulk structure:

I need the conventional unit cell for the slab generator, because otherwise the miller indices will not create the correct slab. However, maybe I can use a transition matrix from the conventional cell to the primitive cell to transform the vector of the miller indices to the basis of the primitive cell and use that to construct the slab? I am not sure if the slab_generator will like to work with non-conventional cells, but maybe this could be an option.

I think you might be skipping a step

1.) Figure out the transformation matrix from primitive to conventional and from primitive to your cell.
2.) Label primitive cell
3.) Generate conventional_cell_supercell and your_cell_supercell from the labeled primitive cell.
4.) Use structure matcher to get the mapping from your_cell to your_cell_supercell. This will let you match the labels to the magmoms.
5.) Transfer magmoms to conventional_cell_supercell or use mapping from conventional_cell to conventional_cell_supercell to transfer magmoms the two cells are not the same shape.

It’s not pretty but basically you’re labeling atoms in the primitive cell and then generating the conventional cell and the MP cell so that you can use these labels to map site-properties around.

True, I thought that the number of atoms could be an issue… My idea may also need to be combined with creating a supercell geometry (with equal or more atoms than the conventional cell) and then after the transformation, truncating the structure back to the original lattice lengths. But maybe there is a more elegant way to do this?
The latter idea you mentioned sounds interesting but I think the problem is that the transformation might lead to Miller indices that are not integers. And as you said, I think that would be problematic for the slab creation code.

Hi, I am not quite sure how to follow this description:

2.) Label primitive cell

How do I label a cell? I cannot find anything in the Structure documentation regarding labels? Also, the primitive cell of NiO only has 2 atoms and thus cannot accommodate the AFM-2 magnetic structure.

However, this got me thinking, and I realized that the conventional unit cell is also not able to contain the AFM-2 order where (111) planes are showing alternating up and down spins. Only a 2x2x2 supercell of the conventional structure is large enough.

Looking at the transformation matrix from the initial cell to the conventional cell I realized that it is easy to transform in a matrix consisting of near integers:

from pymatgen import MPRester
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer as sga

struct = MPRester().get_structure_by_material_id('mp-19009')
conv_struct = sga(struct).get_conventional_standard_structure()
my_lattice = struct.lattice
conv_lattice = conv_struct.lattice
#Get the transition matrix from the initial to the conventional lattice
TM_MyToConv = np.dot(my_lattice.inv_matrix, conv_lattice.matrix)
TM_MyToConv
Out[184]: 
array([[ 1.49948877,  0.4994643 , -0.50051017],
       [-0.49946216, -1.49948427, -0.50050899],
       [-0.49950552,  0.49950458, -0.50050116]])

Then I can create a supercell using the scaled (~doubled) transition matrix directly:

scale_factor = 1/np.amin(np.absolute(TM_MyToConv))
SC_matrix = TM_MyToConv*scale_factor
my_supercell = struct.copy()
my_supercell.make_supercell(SC_matrix)

(Note that I get an LinAlgError: Singular matrix when I do not rescale the matrix first.)
Then I create a supercell for the conventional structure, with the same scaling factor, so that I end up with a 64 atom supercell also.
Now the matching works and I can add the correct site properties:

conv_supercell = conv_struct.copy()
conv_supercell.make_supercell(scale_factor)

matcher = StructureMatcher(primitive_cell=False)
matched_list = matcher.get_mapping(my_supercell, conv_supercell)

for i, site in enumerate(conv_supercell.sites):
    mag_moment = my_supercell.site_properties['magmom'][matched_list[i]]
    conv_supercell.replace(i, conv_supercell.species[i], 
                           properties={'magmom': mag_moment})

This results in a conventional cell with the correct magnetic structure, but I am pretty certain that this is not a general way to go. I Tried the same with a simple bcc Fe cell and the whole thing failed because the scaling factor gets huge and the resulting supercells do not fit into memory :upside_down_face:
Just looking at the code it is clear that it must fail if a single entry in the transition matrix is 0 or close to 0!

I give the whole process that works for NiO (‘mp-19009’) as a function below, so you can play around with this if you want. However, I think I will abandon this issue for now. First I am going on holidays for two weeks, and second, for ferromagnetic crystals the process should be much easier. Just taking the averages of the magnetic moments of each species as magmom flag for all the atoms of the same species in the conventional standard structure will be enough to converge to the correct results.

I guess for now I will leave AFM structures alone…

Anyhow, here is the function that works at least for NiO and fails for Fe maybe someone can find the correct way to scale the transition matrix!

def MakeMagneticStandardStructure(struct):    
    conv_struct = sga(struct).get_conventional_standard_structure()
    my_lattice = struct.lattice
    conv_lattice = conv_struct.lattice
    #Get the tranisition matrix from the initial to the conventional lattice
    TM_MyToConv = np.dot(my_lattice.inv_matrix, conv_lattice.matrix)
    #This scaling factor is not working in general! It will approach inf if one
    #or more entries in the TM are close to 0!
    scale_factor = 1/np.amin(np.absolute(TM_MyToConv))
    SC_matrix = TM_MyToConv*scale_factor
  
    my_supercell = struct.copy()
    my_supercell.make_supercell(SC_matrix)
    
    conv_supercell = conv_struct.copy()
    conv_supercell.make_supercell(scale_factor)

    matcher = StructureMatcher(primitive_cell=False)
    matched_list = matcher.get_mapping(my_supercell, conv_supercell)

    for i, site in enumerate(conv_supercell.sites):
        mag_moment = my_supercell.site_properties['magmom'][matched_list[i]]
        conv_supercell.replace(i, conv_supercell.species[i], 
                               properties={'magmom': mag_moment})
    return conv_supercell