Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Cannot extract the right mapping from spglib #228

Closed
XinYu73 opened this issue Aug 26, 2023 · 4 comments
Closed

Cannot extract the right mapping from spglib #228

XinYu73 opened this issue Aug 26, 2023 · 4 comments

Comments

@XinYu73
Copy link

XinYu73 commented Aug 26, 2023

the poscar

 C 
 1.0000000000000000
     8.5033854756483755    0.0000000000000000    0.0000000000000000
     0.0000000000000000    8.6029826114743475    0.0000000000000000
    -4.1863778586472247    0.0000000000000000    4.2881075098659860
 C  
  28
Cartesian
  0.0854925514737253  8.5942204736845600  4.2827473754786523
 -3.1297313316973261  7.5833872228015471  3.2125215031663004
  0.0555787782448595  6.4599538340082523  2.1613905735953804
  2.1786814807973376  0.0087621377897869  2.1386936205456597
  3.2151503354504740  5.3210866944099742  1.0684677482333074
  6.4004604453926603  6.4445200832032690  0.0173368186623873
  2.2202396384378775  4.4154076999870098  4.2873399386217201
  0.1114587507761341  6.4354482380394664  4.2066163147484925
  2.0854444328988992  6.4305015230378695  2.2167414653327300
  4.2360849769630384  0.1048832625077887  2.1462878589456320
  4.3134285677614903  4.1875749114873377  2.1432861836887263
  6.4563404179239345  6.4690256791720548  2.0625625598154986
  8.4303261000466989  6.4739723941736518  0.0726877103997371
  6.3292739062866508  8.4980993489665586  0.0022341040126392
 -4.1662001863504621  4.2927291679473871  4.2827473754786523
  1.1219614061268617  3.2818959170643729  3.2125215031663004
  4.3072715160690480  2.1584625282710785  2.1613905735953804
 -2.0730112570268497  4.3102534435269604  2.1386936205456597
 -1.0365424023737129  1.0195953886728006  1.0684677482333074
  2.1487677075684717  2.1430287774660952  0.0173368186623873
 -2.0314530993863102  0.1139163942498368  4.2873399386217201
  4.3631514886003222  2.1339569323022940  4.2066163147484925
  6.3371371707230857  2.1290102173006962  2.2167414653327300
 -0.0156077608611481  4.4063745682449627  2.1462878589456320
  0.0617358299373028  8.4890662172245115  2.1432861836887263
  2.2046476800997468  2.1675343734348798  2.0625625598154986
  4.1786333622225103  2.1724810884364776  0.0726877103997371
  2.0775811684624643  4.1966080432293849  0.0022341040126392

I run

crys = pyxtal()
crys.from_seed(crys_atom, tol=0.1, a_tol=5)

and I got following error

0
5
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[17], line 145
    131 # one_crystal(
    132 #     crys_atom,
    133 #     'mp-1182070',
   (...)
    141 #     5,
    142 # )
    144 crys = pyxtal()
--> 145 crys.from_seed(crys_atom, tol=0.1, a_tol=5)

File ~/workplace/PyXtal/pyxtal/__init__.py:409), in pyxtal.from_seed(self, seed, molecules, tol, a_tol, ignore_HH, add_H, backend, style, hn, standard)
    407     from pymatgen.io.ase import AseAtomsAdaptor
    408     pmg_struc = AseAtomsAdaptor.get_structure(seed)
--> 409     self._from_pymatgen(pmg_struc, tol, a_tol, style=style)
    410 elif isinstance(seed, Structure): #Pymatgen
    411     self._from_pymatgen(seed, tol, style=style)

File ~/workplace/PyXtal/pyxtal/__init__.py:485), in pyxtal._from_pymatgen(self, struc, tol, a_tol, style, hn)
    483     print(len(atom_sites))
    484     print(len(sym_struc.equivalent_sites))
--> 485     raise RuntimeError("Cannot extract the right mapping from spglib")
    486 else:
    487     self.atom_sites = atom_sites

RuntimeError: Cannot extract the right mapping from spglib

but findsym gives the result

data_findsym-output
_audit_creation_method FINDSYM

_cell_length_a    8.5033900000
_cell_length_b    8.6029800000
_cell_length_c    5.9928000000
_cell_angle_alpha 90.0000000000
_cell_angle_beta  134.3120000000
_cell_angle_gamma 90.0000000000
_cell_volume      313.6957407842

_symmetry_space_group_name_H-M "C 1 c 1"
_symmetry_Int_Tables_number 9
_space_group.reference_setting '009:C -2yc'
_space_group.transform_Pp_abc a,b,c;0,0,0

loop_
_space_group_symop_id
_space_group_symop_operation_xyz
1 x,y,z
2 x,-y,z+1/2
3 x+1/2,y+1/2,z
4 x+1/2,-y+1/2,z+1/2

loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_symmetry_multiplicity
_atom_site_Wyckoff_label
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_occupancy
_atom_site_fract_symmform
C1 C   4 a 0.00176  0.49898 -0.00125 1.00000 Dx,Dy,Dz 
C2 C   4 a 0.50077  0.38148 0.74917  1.00000 Dx,Dy,Dz 
C3 C   4 a 0.75469  0.25090 0.50404  1.00000 Dx,Dy,Dz 
C4 C   4 a 0.25333  0.01324 -0.00018 1.00000 Dx,Dy,Dz 
C5 C   4 a -0.00393 0.24805 -0.01900 1.00000 Dx,Dy,Dz 
C6 C   4 a -0.00025 0.24747 0.51695  1.00000 Dx,Dy,Dz 
C7 C   4 a 0.74458  0.01219 0.50052  1.00000 Dx,Dy,Dz 
# end of cif
@XinYu73
Copy link
Author

XinYu73 commented Aug 26, 2023

I suppose that there maybe some error in _from_pymatgen

#if wp.index>0: print(wp)
                pos1 = wp.search_generator(pos, self.group[0], tol=tol)
                if pos1 is not None:
                    atom_sites.append(atom_site(wp, pos1, specie))
                else:
                    break

the way we choose representation location in an orbital may be wrong

@XinYu73
Copy link
Author

XinYu73 commented Aug 26, 2023

I have written a snippet of code to transform the parent crystal to a derivative crystal according to transformation matrix, I choose another way to select the representative site in one orbital

def transform_crystal_ugly_way(parent_crystal, Pmat, subgroup_num):
    Rns = get_R_in_P(Pmat)
    cry_ase = parent_crystal.to_ase()
    # get all scaled position in supercell, first 3 rows of site_mat are positions, the last row is element number
    site_mat = np.zeros((4, len(cry_ase.symbols.numbers) * Rns.shape[1]))
    for i in range(Rns.shape[1]):
        site_mat[
            :3,
            i * len(cry_ase.symbols.numbers) : (i + 1) * len(cry_ase.symbols.numbers),
        ] = (
            cry_ase.get_scaled_positions().T
            + np.tile(Rns[:, i], (len(cry_ase.symbols.numbers), 1)).T
        )
        site_mat[
            -1,
            i * len(cry_ase.symbols.numbers) : (i + 1) * len(cry_ase.symbols.numbers),
        ] = cry_ase.symbols.numbers
    site_mat[:3, :] = np.linalg.inv(Pmat[:3, :3]) @ (
        site_mat[:3, :] - np.tile(Pmat[:3, -1], (site_mat.shape[1], 1)).T
    )
    site_mat[:3, :] = rescale_zero_one(site_mat[:3, :])
    _, index = np.unique(site_mat.round(decimals=5), axis=1, return_index=True)
    site_mat = site_mat[:, index]
    # split all sites into several orbitals, through the orbital labels are not known at the moment
    equivalent = [-1] * site_mat.shape[1]
    for i in range(site_mat.shape[1]):
        if equivalent[i] != -1:
            continue
        else:
            other_site_on_orb = [
                rescale_zero_one(op.operate(site_mat[:3, i]))
                for op in Group(subgroup_num).wyckoffs[0]
            ]
            # any symm operation can only translate one site to places on the orbitals it belongs to
            # so we set the sites that could be related by symm ops with the mark(!=-1)
            for j in range(len(other_site_on_orb)):
                diff = (
                    site_mat[:3, :]
                    - np.tile(other_site_on_orb[j], (site_mat.shape[1], 1)).T
                )
                equivalent[np.argmin(np.sum(np.abs(diff), axis=0))] = i
    equivalent = np.array(equivalent)
    # set correct orbital label
    site_dict = {}
    for i in set(equivalent):
        site_dict[i] = [site_mat[:, equivalent == i]]
        for wyckoff in Group(subgroup_num).Wyckoff_positions:
            for site_index_on_orb in range(site_dict[i][0].shape[1]):
                # so, why there is a loop over sites on each orbital?
                # Group(49).Wyckoff_positions[1]
                """
                Wyckoff position 4q in space group 49 with site symmetry ..m
                x, y, 0
                -x, -y, 0
                -x, y, 1/2
                x, -y, 1/2
                """
                # Group(49).Wyckoff_positions[1].ops
                """
                [Rot:
                [[1. 0. 0.]
                [0. 1. 0.]
                [0. 0. 0.]]
                tau[0. 0. 0.],
                Rot:
                [[-1.  0.  0.]
                [ 0. -1.  0.]
                [ 0.  0.  0.]]
                tau
                [0. 0. 0.],
                Rot:
                [[-1.  0.  0.]
                [ 0.  1.  0.]
                [ 0.  0.  0.]]
                tau
                [0.  0.  0.5],
                Rot:
                [[ 1.  0.  0.]
                [ 0. -1.  0.]
                [ 0.  0.  0.]]
                tau
                [0.  0.  0.5]]
                """
                # say you have
                # (0.3,0.4,0),(0.7,0.6,0),(0.7,0.4,0.5),(0.3,0.6,0.5)
                # and you start with (0.7,0.4,0.5)
                # you will get
                # (0.7,0.4,0.0),(0.3,0.6,0.0),(0.7,0.4,0.5),(0.7,0.6,0.5)
                # which is not the same orb
                other_site_on_orb = [
                    rescale_zero_one(op.operate(site_dict[i][0][:3, site_index_on_orb]))
                    for op in wyckoff.ops
                ]
                other_site_on_orb = np.array(other_site_on_orb).T
                true_label, _ = test_if_2_array_column_wise_same(
                    other_site_on_orb, site_dict[i][0][:3, :]
                )
                # translate ith site to sites using ops of chosen wyckoff position
                if true_label:
                    site_dict[i].extend(
                        [str(wyckoff.multiplicity) + wyckoff.letter, site_index_on_orb]
                    )
                    break
    site_dict = dict(sorted(site_dict.items(), key=lambda x: -x[1][0][-1, 0]))
    # sort sites according to element number
    new_cry = pyxtal()
    # new_cry.build(
    #     group=subgroup_num,
    #     species=[
    #         Element(int(item[0][-1, 0])).short_name for item in site_dict.values()
    #     ],
    #     numIons=[int(item[1][:-1]) for item in site_dict.values()],
    #     lattice=Lattice.from_matrix(
    #         np.dot(Pmat[:3, :3].T, cry_ase.cell.array),
    #         ltype=Group(subgroup_num).lattice_type,
    #     ),
    #     sites=[{item[1]: item[0][:3, item[2]]} for item in site_dict.values()],
    # )
    new_cry.build(
        group=subgroup_num,
        species=[
            Element(int(item[0][-1, 0])).short_name for item in site_dict.values()
        ],
        numIons=[int(item[1][:-1]) for item in site_dict.values()],
        lattice=Lattice.from_matrix(
            np.dot(Pmat[:3, :3].T, cry_ase.cell.array),
            ltype=Group(subgroup_num).lattice_type,
        ),
        sites=[[{item[1]: item[0][:3, item[2]]}] for item in site_dict.values()],
    )
    return new_cry

@qzhu2017
Copy link
Owner

One more failed example

 C 
 1.0000000000000000
     4.1896699999999996    0.0000000000000003    0.0000000000000003
     0.0000000000000000    5.9243990000000002    0.0000000000000004
     0.0000000000000000    0.0000000000000000    5.9243949999999996
 C  
  16
Cartesian
  3.1422524999999997  0.7405202530050002  1.4810987500000001
 -3.1422524999999997 -0.7405202530050002  1.4810987499999997
 -3.1422524999999997  0.7405202530049999 -1.4810987500000001
  3.1422524999999997 -0.7405202530049999 -1.4810987499999997
  5.2370874999999995  3.7027197530050002  1.4810987500000004
 -1.0474174999999999  2.2216792469950004  1.4810987499999999
 -1.0474174999999999  3.7027197530049998 -1.4810987499999997
  5.2370874999999995  2.2216792469950004 -1.4810987499999995
  0.0000000000000000  0.0000000000000000  5.1838160030249991
  0.0000000000000000  0.0000000000000000 -5.1838160030249991
  2.0948349999999998  2.9621995000000001  5.1838160030250009
  2.0948349999999998  2.9621995000000001 -5.1838160030249991
  0.0000000000000000  2.9621995000000001  3.7027824213699998
  0.0000000000000000  2.9621995000000001 -3.7027824213699994
  2.0948349999999998  5.9243990000000002  3.7027824213699998
  2.0948349999999998  5.9243990000000002 -3.7027824213699989

@qzhu2017
Copy link
Owner

@XinYu73 This issue should have been fixed if you update the code. Let me know if you have any problem.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants