fix a typo for 220 k-subgroup #252 and some other issues in get_tab_rep

pyxtal/__init__.py

@@ -786,6 +786,14 @@ def subgroup_once(self, eps=0.1, H=None, perms=None, group_type='t', \
             a pyxtal structure with lower symmetries
         """
         idx, sites, t_types, k_types = self._get_subgroup_ids(H, group_type, None, max_cell, min_cell)
+        #print("Test subgroup_once", self.group.number, idx, sites)
+
+        if idx is None or len(idx) == 0:
+            msg = "Cannot find the valid splitter id\n"
+            msg += "This is most likely due to the cell size is not sufficient for k-subgroup\n"
+            msg += "You may increase the max_cell"
+            print(msg)
+            return None
 
         # Try 100 times to see if a valid split can be found
         count = 0
@@ -837,7 +845,9 @@ def subgroup_once(self, eps=0.1, H=None, perms=None, group_type='t', \
                             if ans.group.number > 1:
                                 return ans
             count += 1
-        raise RuntimeError("Cannot find the splitter")
+
+        print("Cannot find the splitter")
+        return None
 
     def _apply_substitution(self, splitter, perms):
         """
@@ -901,7 +911,7 @@ def _get_subgroup_ids(self, H, group_type, idx, max_cell, min_cell):
 
         if idx is None:
             idx = []
-            if not self.molecular or self.group.number>142:
+            if not self.molecular or self.group.number > 142:
                 for i, tran in enumerate(trans):
                     if min_cell<=np.linalg.det(tran[:3,:3])<=max_cell:
                         idx.append(i)
@@ -925,8 +935,8 @@ def _get_subgroup_ids(self, H, group_type, idx, max_cell, min_cell):
         if H is not None:
             idx = [id for id in idx if Hs[id] == H]
 
-        if len(idx) == 0:
-            raise RuntimeError("Cannot find the splitter")
+        #if len(idx) == 0:
+        #    raise RuntimeError("Cannot find the splitter")
 
         if self.molecular:
             struc_sites = self.mol_sites
@@ -3353,9 +3363,9 @@ def get_tabular_representation(self, ids=None, normalize=True, N_wp=6, max_abc=5
         """
 
         if ids is None:
-            ids = [0] * len(xtal.atom_sites)
+            ids = [0] * len(self.atom_sites)
         else:
-            assert(len(ids) == len(xtal.atom_sites))
+            assert(len(ids) == len(self.atom_sites))
 
         if len(ids) > N_wp:
             raise ValueError("Cannot handle too many atom sites", len(ids), N_wp)
@@ -3364,18 +3374,18 @@ def get_tabular_representation(self, ids=None, normalize=True, N_wp=6, max_abc=5
 
         # If the number of wp is less then 6, assign -1
         rep = -np.ones(N_cell + N_site) # 7 + 6 + 18 = 31
-        rep[0] = xtal.group.number
-        rep[1:7] = xtal.lattice.get_para()
+        rep[0] = self.group.number
+        rep[1:7] = self.lattice.get_para()
 
         if normalize:
             rep[0] /= 230
             rep[1:4] /= max_abc
             rep[4:7] /= max_angle
 
         count = N_cell
-        for id, site in zip(ids, xtal.atom_sites):
+        for id, site in zip(ids, self.atom_sites):
             rep[count] = site.wp.index
-            if normalize: rep[count] /= len(xtal.group)
+            if normalize: rep[count] /= len(self.group)
             xyz = site.coords[id] # position
             xyz -= np.floor(xyz)
             rep[count+1:count+4] = xyz

pyxtal/symmetry.py

@@ -777,7 +777,7 @@ def get_min_supergroup(self, group_type='t', G=None):
             for sg in sgs:
                 subgroups = None
                 if group_type == 't':
-                    if sg>self.number:
+                    if sg > self.number:
                         subgroups = Group(sg, quick=True).get_max_t_subgroup()
                 else:
                     g1 = Group(sg)

pyxtal/wyckoff_split.py

@@ -102,26 +102,26 @@ def parse_wp2(self, idx):
         """
         query the wp2 and transformation matrix from the given {G, H, wp1}
         """
-        #print("trans", idx)
         #print(self.wyc['transformation'])
         #subgroup_relations.reverse()
-        trans = self.wyc['transformation'][idx]
-        subgroup_relations = self.wyc['relations'][idx]
+        trans = self.wyc['transformation'][idx]#; print("trans", trans)
+        subgroup_relations = self.wyc['relations'][idx] #; print('subgroup_relations', subgroup_relations)
         subgroup_relations = list(reversed(subgroup_relations))
 
-        self.R = np.zeros([4,4])
-        self.R[:3,:3] += trans[:3,:3]
-        self.R[3,3] = 1
+        self.R = np.zeros([4, 4])
+        self.R[:3, :3] += trans[:3, :3]
+        self.R[3, 3] = 1
         self.inv_R = np.linalg.inv(self.R)
-        inv_t = np.dot(self.inv_R[:3,:3], trans[:,3].T)
-        self.inv_R[:3,3] = -1*inv_t.T
-        self.R[:3,3] = trans[:3,3]
+        inv_t = np.dot(self.inv_R[:3, :3], trans[:,3].T)
+        self.inv_R[:3, 3] = -1*inv_t.T
+        self.R[:3, 3] = trans[:3, 3]
         self.multi = np.linalg.det(self.R[:3, :3])
 
         wp2_lists = []
         for wp1_index in self.wp1_indices:
             wp2_list = []
             for letter in subgroup_relations[wp1_index]:
+                #if letter == '4a': letter = '8c'
                 id = sym.index_from_letter(letter[-1], self.H)
                 wp2_list.append(self.H[id])
             wp2_lists.append(wp2_list)
@@ -299,7 +299,7 @@ def split_t(self, wp1, wp2_lists, quadrant=None):
         return G1_orbits, G2_orbits
 
 
-    def split_k(self, wp1, wp2_lists):
+    def split_k(self, wp1, wp2_lists, tol=1e-5):
         """
         split the generators in w1 to different w2s for k-subgroup
         """
@@ -308,52 +308,57 @@ def split_k(self, wp1, wp2_lists):
 
         G1_orbits = []
         G2_orbits = []
-        quadrant=deepcopy(self.inv_R[:3,3])
-        quadrant[np.abs(quadrant)<1e-5] = 0 #finds the orientation of the subgroup_basis
+        quadrant = deepcopy(self.inv_R[:3, 3])
+        quadrant[np.abs(quadrant) < tol] = 0 #finds the orientation of the subgroup_basis
+
         for i in range(3):
-            if quadrant[i]>=0:
+            if quadrant[i] >= 0:
                 quadrant[i] = 1
             else:
                 quadrant[i] = -1
 
         all_g2_orbits = []
         translations = self.translation_generator()
+
         # the translation generator provides all the possible ways to translate
         # the starting positions, then they are shifted
         for translation in translations:
             for gen in wp1_generators:#into the proper orientation
                 orbit = np.matmul(self.inv_R,gen)
-                orbit[np.abs(orbit)<1e-5] = 0
-                orbit[np.abs(orbit-1)<1e-5] = 1
-                orbit[np.abs(orbit+1)<1e-5] = -1
+                orbit[np.abs(orbit) < tol] = 0
+                orbit[np.abs(orbit-1) < tol] = 1
+                orbit[np.abs(orbit+1) < tol] = -1
                 for i in range(3):
                     if quadrant[i] == 1:
                         orbit[i][3] += translation[i]
-                        orbit[i][3] = orbit[i][3]%1
-                        if np.abs(orbit[i][3]-1) < 1e-5:
+                        orbit[i][3] = orbit[i][3] % 1
+                        if np.abs(orbit[i][3]-1) < tol:
                             orbit[i][3] = 0
                     else:
                         orbit[i][3] += (translation[i])%-1
-                        orbit[i][3] = orbit[i][3]%-1
-                        if np.abs(orbit[i][3]) < 1e-5:
+                        orbit[i][3] = orbit[i][3] % -1
+                        if np.abs(orbit[i][3]) < tol:
                             orbit[i][3] = -1
                 all_g2_orbits.append(orbit)
 
         for wp2 in wp2_lists:
 
             #final_G2=[]
             temp = np.array(deepcopy(all_g2_orbits))
-            temp[np.abs(temp) <1e-5] = 0
+            temp[np.abs(temp) < tol] = 0
             temp = temp.tolist()
 
             for j, x in enumerate(temp):
                 temp[j] = SymmOp(x)
 
             for orbit in temp:
                 try_match = np.array([np.matmul(x.as_dict()['matrix'], orbit.as_dict()['matrix']) for x in wp2])
-                try_match[np.abs(try_match) <1e-5] = 0
-                try_match[np.abs(try_match-1) <1e-5] = 1
-                try_match[np.abs(try_match+1) <1e-5] = -1
+                try_match[np.abs(try_match) < tol] = 0
+                try_match[np.abs(try_match-1) < tol] = 1
+                try_match[np.abs(try_match+1) < tol] = -1
+
+                #print('\norbit\n', orbit)
+                #print('\ntry_match\n', try_match)
 
                 for j in range(len(try_match)):
                     for k in range(3):
@@ -385,7 +390,10 @@ def split_k(self, wp1, wp2_lists):
             G1_orbits.append(final_G1)
 
         if len(G1_orbits) != len(wp2_lists):
-            raise ValueError('inaccurate')
+            print('wp2_lists', wp2_lists)
+            print('G1_orbits', G1_orbits)
+            print('G2_orbits', G2_orbits)
+            raise ValueError('inconsistent G1_orbits and wp2_lists')
         else:
             return G1_orbits, G2_orbits
 
@@ -497,5 +505,10 @@ def in_lists(mat1, mat2, eps=1e-2, PBC=True):
 
 if __name__ == "__main__":
 
-    sp = wyckoff_split(G=14, idx=1, wp1=['2c', '4e'], group_type='t')
-    print(sp)
+    #sp = wyckoff_split(G=14, idx=1, wp1=['2c', '4e'], group_type='t')
+    #print(sp)
+    for idx in range(4):
+        #sp = wyckoff_split(G=210, idx=idx, wp1=['8b'], group_type='k')
+        sp = wyckoff_split(G=224, idx=idx, wp1=['24j'], group_type='k')
+        #print(sp)
+        print(sp.error)