Inconsistency in SC->UC transformation operations

[materialsvirtuallab]

Inconsistency in SC->UC transformation operations

[jmmshn]

Summary of the problem

Based on the information in the notebook it looks like we took a UC structure and two endpoint positions:

base_struct_uc = Structure.from_file('VOPO4_uc_base.cif') # base unit cell with no Mg
sc_mat = [[2, 0, 0], [0, 2, -1], [0, 0, 2]] # matrix used to do UC->SC stored in DB
hop_isite = np.array([0.58798658, 0.14402022, 0.95959355]) # isite frac_coords in terms of UC, stored in DB
hop_esite = np.array([0.06753297, -0.33552617, 1.45959355]) # esite frac_coords in terms of UC, stored in DB

then created a supercell containing those two end points

sc_mat_inv = np.linalg.inv(sc_mat)
ipos_sc = np.dot(hop_isite, sc_mat_inv)
epos_sc = np.dot(hop_esite, sc_mat_inv)

inserted_sc = base_struct_uc.copy() * sc_mat
inserted_sc.insert(0, 'Mg', ipos_sc)
inserted_sc.insert(1, 'Mg', epos_sc)

at this point we calculate the supercell mapping sc_m and total_t using StructureMatcher

sm = StructureMatcher(ignored_species=["Mg"])
mapping = get_matched_structure_mapping(base=base_struct_uc, inserted=inserted_sc, sm=sm)
isite = inserted_sc[0]
esite = inserted_sc[1]
uc = base_struct_uc.copy()

if mapping is None:
    raise ValueError("Cannot obtain inverse mapping, consider lowering tolerances " "in StructureMatcher")
sc_m, total_t = mapping

which has a different (but equivalent) supercell transformation compared to sc_mat

ss1 = base_struct_uc * sc_m
print(ss1.lattice)
ss2 = base_struct_uc * sc_mat
print(ss2.lattice)

output:

# 
-0.592209 9.676808 3.136251
9.694910 0.000000 -4.250698
0.000000 0.000000 14.773895

9.694910 0.000000 3.136250
-0.592209 9.676808 -4.250697
0.000000 0.000000 14.773895

Then there appears to be some code in the notebook that tries to map positions in the inserted_sc structure back to the UC using only the trasformations derived from pymatgen.analysis.diffusion.utils.parse_entries.get_matched_structure_mapping.

sm = StructureMatcher(ignored_species=["Mg"])
mapping = get_matched_structure_mapping(base=base_struct_uc, inserted=inserted_sc, sm=sm)
isite = inserted_sc[0]
esite = inserted_sc[1]
uc = base_struct_uc.copy()

if mapping is None:
    raise ValueError("Cannot obtain inverse mapping, consider lowering tolerances " "in StructureMatcher")
sc_m, total_t = mapping
sc_ipos = isite.frac_coords
sc_ipos_t = sc_ipos - total_t
uc_ipos = sc_m.dot(sc_ipos_t)
image_trans = np.floor(uc_ipos)
uc_ipos = uc_ipos - image_trans
uc_ipos = _get_first_close_site(uc_ipos, uc)

sc_epos = esite.frac_coords
sc_epos_t = sc_epos - total_t
uc_epos = sc_epos_t.dot(sc_m)
uc_epos = uc_epos - image_trans
uc_epos = _get_first_close_site(uc_epos, uc)

This is not totally valid since the there are many different transformations between the same SC and UC so you are not always going to get back to the same position.

A better way to check if things are working correctly is to do:

from pymatgen.analysis.diffusion.utils.edge_data_from_sc import get_uc_pos
i,m,e = get_uc_pos(isite=inserted_sc.sites[0], esite=inserted_sc.sites[1], uc = base_struct_uc, sc = inserted_sc, sm = sm)

which will guarantee a consisten mapping.

uc = base_struct_uc.copy()
uc.sites.append(i)
orign_inserted = base_struct_uc.copy()
orign_inserted.insert(0, "Mg", hop_isite)
sm_w.fit(uc, orign_inserted)

which will return True.

As far as I can see, get_uc_pos is still working properly, I think if want to figure out exactly why @acrutt 's code is breaking we should make multiple different calls to get_uc_pos and see where the i,m,e positions are carefully.

Let me know if this helps, and if we should discuss more.

[hmlli]

To summarize, in this particular case, get_uc_pos (an equivalent method to the one originally used in the notebook uploaded, which was copy-pasted from get_uc_pos) takes a pair of sites in SC, and returns a pair of sites in UC. However, this returned pair is not the original UC coords that were used to created the SC sites in the first, but a pair of symmetrically equivalent sites in UC. (This can be confirmed, as mentioned above, with the fit method from StructureMatcher. It was also visually inspected.)

With all above in mind, the issue becomes that inability to initiate MigrationHop with the returned pair from get_uc_pos. Specifically, the following:

mh_from_sc = MigrationHop(uc_isite, uc_esite, symm_structure=mg.symm_structure)

will return an error:

#
~/uw/repo/pymatgen-diffusion/pymatgen/analysis/diffusion/neb/pathfinder.py in __init__(self, isite, esite, symm_structure, symprec)
    344 
    345         if self.iindex is None:
--> 346             raise RuntimeError(f"No symmetrically equivalent site was found for {isite}")
    347         if self.eindex is None:
    348             raise RuntimeError(f"No symmetrically equivalent site was found for {esite}")

RuntimeError: No symmetrically equivalent site was found for [3.19210999 0.29513241 8.2409902 ] Mg

HOWEVER, as I was trying to put the symm_structure into a file to send here, I've found something interesting.

The symm_structure object itself has an as_dict function, but the returned dict object cannot be json dumped because something in there is not json serializable. Therefore, all the information needed to initiated the symm_structure object should be exacted. In particular, the spacegroup number of the structure returned international number 1, BUT the structure itself on MP is actually number 9.

For reference, one other compound that went through the same process has number 4 from the mg.symm_structure object, and it indeed has number 4 from MP.

[hmlli]

The error can be reproduced by the following:

MG_struct = Structure.from_file("MG_symm_struct.vasp")
a = SpacegroupAnalyzer(MG_struct, symprec=0.1)
MG_symm_struct = a.get_symmetrized_structure()
mh_from_sc = MigrationHop(i, e, symm_structure=MG_symm_struct)

it will give the error:

#
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-58-cc31c91cc9ce> in <module>
----> 1 mh_from_sc = MigrationHop(i, e, symm_structure=MG_symm_struct)

~/uw/repo/pymatgen-diffusion/pymatgen/analysis/diffusion/neb/pathfinder.py in __init__(self, isite, esite, symm_structure, symprec)
    344 
    345         if self.iindex is None:
--> 346             raise RuntimeError(f"No symmetrically equivalent site was found for {isite}")
    347         if self.eindex is None:
    348             raise RuntimeError(f"No symmetrically equivalent site was found for {esite}")

RuntimeError: No symmetrically equivalent site was found for [0.52402612 2.84491651 4.54287191] Mg

MG_symm_struct.vasp.zip

[jmmshn]

OK so summary of further investigation:

inserted1 = base_struct_uc.copy()  # using the data from the approxNEB database
inserted2 = base_struct_uc.copy() # using the data from the current database
inserted1.insert(0, 'Mg', hop_isite)
inserted1.insert(1, 'Mg', hop_esite)
inserted2.insert(0, 'Mg', i.frac_coords) # i, and e as defined above
inserted2.insert(1, 'Mg', e.frac_coords)
sm_w = StructureMatcher()
sm_w.fit(inserted1, inserted2)  # returns true

We see that these inserted structures at the same.
Additionally, we see can call this function to obtain the transformation between inserted1 and inserted2.
https://github.com/materialsvirtuallab/pymatgen-diffusion/blob/c3e80e90f157c7bf586d36b3e558163058503183/pymatgen/analysis/diffusion/utils/parse_entries.py#L132

mapping = get_matched_structure_mapping(base=inserted1, inserted=inserted2, sm=sm_w)
print(mapping)
# (array([[0, 1, 0],
#       [1, 0, 0],
#       [0, 0, 1]]), array([-8.63725124e-18,  1.06332807e-17,  5.00000000e-01]))

Since is operation is in the set of operations we find for the base:

a_base = SpacegroupAnalyzer(base_struct_uc, symprec=0.1)
a_base.get_space_group_operations()

which gives:

[Rot:
 [[1. 0. 0.]
  [0. 1. 0.]
  [0. 0. 1.]]
 tau
 [0. 0. 0.],
 Rot:
 [[0. 1. 0.]
  [1. 0. 0.]
  [0. 0. 1.]]
 tau
 [0.  0.  0.5]]

It looks like the symmetry part of the code is doing what it's supposed to.
However, the sites that correspond to the endpoints of hop_isite are missing from the MG_symm_struct.vasp file.
Can you send a JSON of the entries used to generate that?

So my current hypothesis is that this function is causing problem: https://github.com/materialsvirtuallab/pymatgen-diffusion/blob/c3e80e90f157c7bf586[…]163058503183/pymatgen/analysis/diffusion/utils/parse_entries.py

The SITE_MERGE_R radius is 1.0 Angstroms which might be causing sites to shift after merging.
Since there are only 2 symmetry operations, I would suggest you look at the plausible endpoints you find from the inserted positions and track how they change with the symmetry operations.
by manually applying the operations:
https://github.com/materialsvirtuallab/pymatgen-diffusion/blob/c3e80e90f157c7bf586d36b3e558163058503183/pymatgen/analysis/diffusion/utils/parse_entries.py#L237

Either way we need to look at how the sites from the entries get mapped more carefully to understand exactly what is going on.

[hmlli]

hop_isite and hop_esite are store in DB, so they are presented as coords in np.array directly, like above:

hop_isite = np.array([0.58798658, 0.14402022, 0.95959355]) # isite frac_coords in terms of UC, stored in DB
hop_esite = np.array([0.06753297, -0.33552617, 1.45959355]) # esite frac_coords in terms of UC, stored in DB

Attached is a JSON file of the base UC structure with hop_isite and hop_eiste inserted.

In the structure generated from MG_symm_struct.vasp file, hop_isite corresponds to site 0 (exact), while hop_eisite corresponds to site 3 (off by image of (0, -1, 1)).

The symmetry operations will take some more time to go through.
inserted_hop_uc.json.zip

[hmlli]

If we focus on the structure inserted2, which has the sites that were results of the UC->SC->UC transformation, we can get the following:

inserted2[0]
# PeriodicSite: Mg (0.5240, 2.8449, 4.5429) [0.1440, 0.5880, 0.4596]

this site is the initial site of the hop, which is also the site in the error message (see comment above).

Now we can manually apply transformation to inserted2:

a_base = SpacegroupAnalyzer(base_struct_uc, symprec=0.1)
[op0, op1] = a_base.get_space_group_operations() # note that op0 is simply the identity transformation, it's trivial here
symmtest_inserted2 = inserted2.copy()
symmtest_inserted2.apply_operation(symmop=op1, fractional=True)
symmtest_inserted2[0]
# PeriodicSite: Mg (0.5240, 2.8449, 8.2363) [0.5880, 0.1440, 0.9596]

We can see that this site after the symmetry operation has the same fractional coordinates as site0 in inserted1, which is from the DB directly (and also is included in the MigrationGraph object:

inserted1[0]
# PeriodicSite: Mg (2.8076, 0.6968, 8.2363) [0.5880, 0.1440, 0.9596]

So it seems like that applying the symmetry operation

Rot:
[[0. 1. 0.]
 [1. 0. 0.]
 [0. 0. 1.]]
tau
[0.  0.  0.5]

to the transformed UC coords will give back the coords in MigrationGraph. That being said, the initiation of MigrationHop somehow does not think that these sites are symmetrically equivalent.

PS: I've tried to change the merge radius in the source code and the MigrationHop initiation still throws an error.

[hmlli]

resolved with #228 .