Fixed bug: Classification of q-points for zone-border points out from the 2x2x2 grid

Modules/Classify.py

@@ -436,10 +436,10 @@ def recognize_triplet(dyn, orbit1a, map_uc, symprec=1e-5, verbose=False):
             Maps an arbitrary triplet to its reference triplet, permutaion, and symmtry operation after
             the classification. Dimension [Natom,Natom_sc,Natom_sc,3]
     """
-    start_time = time.time()
+    # start_time = time.time()
 
-    if verbose:
-        print("===== STARTING TRIPLET CLASSIFICATION =====")
+    # if verbose:
+    #     print("===== STARTING TRIPLET CLASSIFICATION =====")
 
     # Get the symmetries of the unit cell
     spg_syms = spglib.get_symmetry(dyn.structure.get_spglib_cell(), symprec)
@@ -536,21 +536,21 @@ def recognize_triplet(dyn, orbit1a, map_uc, symprec=1e-5, verbose=False):
     nref3 = SCHAModules.module_hess.get_nref3(nat,ntot,tot3,nsym,orbit1a,map_uc)
     orbit3a, orbit3s, norbit, indep_fc, n_indep_fc, kernel, mapping_triplet = SCHAModules.module_hess.recognize_triplet(nat,ntot,tot3,nref3,nsym,orbit1a,map_uc,nontrivial,M,verbose)
 
-    end_time = time.time()
-    execution_time = end_time - start_time
+    # end_time = time.time()
+    # execution_time = end_time - start_time
 
-    # Print detailed info
-    if verbose:
-        print(" ")
-        print("Total triplets:", tot3)
-        print("Number of Orbits:", nref3)
-        print("Total number of independent elements:", sum(n_indep_fc[:nref3]))
-        print(" ")
+    # # Print detailed info
+    # if verbose:
+    #     print(" ")
+    #     print("Total triplets:", tot3)
+    #     print("Number of Orbits:", nref3)
+    #     print("Total number of independent elements:", sum(n_indep_fc[:nref3]))
+    #     print(" ")
 
 
-        print("execution_time in triplet recognition:", execution_time, " s")
-        print("===== TRIPLET CLASSIFICATION and GJ ELIMINATION FINISHED ======")
-        print(" ")
+    #     print("execution_time in triplet recognition:", execution_time, " s")
+    #     print("===== TRIPLET CLASSIFICATION and GJ ELIMINATION FINISHED ======")
+    #     print(" ")
 
     sys.stdout.flush()
 
@@ -603,10 +603,10 @@ def recognize_quadruplet(dyn, orbit1a, map_uc, verbose=False, symprec=1e-5):
             the classification. Dimension [Natom,Natom_sc,Natom_sc,3]
     """
 
-    start_time = time.time()
+    # start_time = time.time()
 
-    if verbose:
-        print("===== STARTING QUADRUPLET CLASSIFICATION =====")
+    # if verbose:
+    #     print("===== STARTING QUADRUPLET CLASSIFICATION =====")
 
     # Get the symmetries of the unit cell
     spg_syms = spglib.get_symmetry(dyn.structure.get_spglib_cell(), symprec)
@@ -671,19 +671,19 @@ def recognize_quadruplet(dyn, orbit1a, map_uc, verbose=False, symprec=1e-5):
     nref4 = SCHAModules.module_hess.get_nref4(nat,ntot,tot4,nsym,orbit1a,map_uc)
     orbit4a, orbit4s, norbit, indep_fc, n_indep_fc, kernel, mapping_quadruplet = SCHAModules.module_hess.recognize_quadruplet(nat,ntot,tot4,nref4,nsym,orbit1a,map_uc,nontrivial,M,verbose)
 
-    end_time = time.time()
-    execution_time = end_time - start_time
+    # end_time = time.time()
+    # execution_time = end_time - start_time
 
-    if verbose:
-        print(" ")
-        print("Total quadruplets:", tot4)
-        print("Number of Orbits:",nref4)
-        print("Total number of independent elements:", sum(n_indep_fc[:nref4]))
-        print(" ")
+    # if verbose:
+    #     print(" ")
+    #     print("Total quadruplets:", tot4)
+    #     print("Number of Orbits:",nref4)
+    #     print("Total number of independent elements:", sum(n_indep_fc[:nref4]))
+    #     print(" ")
 
-        print("execution_time in quadruplet recognition:", execution_time, " s")
-        print("===== QUADRUPLET CLASSIFICATION and GJ ELIMINATION FINISHED ======")
-        print(" ")
+    #     print("execution_time in quadruplet recognition:", execution_time, " s")
+    #     print("===== QUADRUPLET CLASSIFICATION and GJ ELIMINATION FINISHED ======")
+    #     print(" ")
     return(orbit4a[:nref4], orbit4s[:nref4], norbit[:nref4], indep_fc[:nref4], n_indep_fc[:nref4], kernel[:nref4], Rot, mapping_quadruplet)
 
 # Small functions required in main subroutines.

Modules/Ensemble.py

@@ -3869,7 +3869,7 @@ def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, eps = 1
         nat_sc = dyn_supercell.structure.N_atoms
         n_modes = nat_sc*3
 
-        mapping, rot_cart, map_uc, map_tr, T_list, T_list_frac = Classify.map_singlet(self.current_dyn, verbose = verbose)
+        mapping, rot_cart, map_uc, map_tr, T_list, T_list_frac = Classify.map_singlet(self.current_dyn, verbose = False)
 
         orbit2a, orbit2s, norbit, indep_elem, n_indep_elem, tensor = Classify.recognize_doublet(self.current_dyn, mapping, map_uc, verbose = verbose)
 
@@ -4712,5 +4712,4 @@ def _wrapper_julia_vector_vector_fourier(*args, **kwargs):
     """
 
     return julia.Main.multiply_vector_vector_fourier(*args,
-            **kwargs)
-
+            **kwargs)

Modules/qClassify.py

@@ -51,25 +51,31 @@ def map_singlet(q_list, q_list_frac, rcell, rot_cart):
 
     q_list_frac_fixed = np.empty(q_list_frac.shape)
     its_zb = np.empty(q_list_frac.shape[0], dtype=np.int8)
+    count1, count2 = 0, 0 
     for qi,q in enumerate(q_list):
-        q_list_frac_fixed[qi] = np.round(CC.Methods.cart_to_cryst(rcell, _map_q_to_1st_bz(rcell, q)),6)
+        q_list_frac_fixed[qi] = np.round(CC.Methods.cart_to_cryst(rcell, _map_q_to_frac_bz(rcell, q)[0]),6)
         if np.all((np.abs(np.round(2*q_list_frac_fixed[qi], 6))%1)<1e-6):
             its_zb[qi] = 0
         else:
             for alpha in range(3):
-                if np.abs(q_list[qi,alpha])>1e-6:
-                    if q_list[qi,alpha] > 1e-6:
+                if np.abs(q_list_frac_fixed[qi,alpha])>1e-6 and np.abs(q_list_frac_fixed[qi,alpha]+0.5)>1e-6:
+                    if q_list_frac_fixed[qi,alpha] > 1e-6:
                         its_zb[qi] = 1 # Non zone border. Positive class.
+                        count1 += 1
                         break
                     else:
                         its_zb[qi] = 2 # Non zone border. Negative class.
+                        count2 += 1
                         break
+    if count1 != count2:
+        print("q-point classification failing...")
+        raise ValueError
 
     mapping = np.zeros([len(q_list), rot_cart.shape[0]], dtype=np.int32)
     for qi, q in enumerate(q_list):
         for isym in range(rot_cart.shape[0]):
             q_sym_cart = rot_cart[isym] @ q
-            q_sym_cart_1bz = _map_q_to_1st_bz(rcell, q_sym_cart)
+            q_sym_cart_1bz = _map_q_to_frac_bz(rcell, q_sym_cart)[0]
             q_sym = np.round(CC.Methods.cart_to_cryst(rcell, q_sym_cart_1bz),6)
             match = np.all(np.abs(q_list_frac_fixed-q_sym)<1e-3, axis=1)
             qii = np.where(match)
@@ -336,6 +342,45 @@ def _map_q_to_1st_bz(rcell, q_cart, atol=1e-5):
 
     return tied_candidates[idx[0]]
 
+def _map_q_to_frac_bz(rcell, q_cart, atol=1e-5):
+    """
+    Maps a Cartesian q-vector into the [-0.5, 0.5) fractional reciprocal unit cell.
+    Forces +0.5 to -0.5 so opposite BZ faces map consistently.
+
+    Parameters
+    ----------
+    rcell : np.ndarray
+        Reciprocal unit cell vectors [3, 3].
+    q_cart : np.ndarray
+        Cartesian coordinates of wave-vector [3].
+    atol : float
+        Tolerance for floating point noise near boundary.
+
+    Returns
+    -------
+    q_cart_mapped : np.ndarray
+        Mapped Cartesian coordinates [3].
+    q_frac_mapped : np.ndarray
+        Mapped fractional coordinates in [-0.5, 0.5) [3].
+    """
+    q_cart = np.array(q_cart).flatten()
+    rcell = np.array(rcell)
+
+    inv_rcell = np.linalg.inv(rcell)
+    q_frac = q_cart @ inv_rcell
+
+    q_frac_round = np.round(q_frac, 6)
+
+    # Shift coordinates into the [-0.5, 0.5) range
+    q_frac_mapped = q_frac_round - np.floor(q_frac_round + 0.5)
+
+    # Force values near +0.5 to strictly -0.5
+    q_frac_mapped = np.where(np.abs(q_frac_mapped - 0.5) < atol, -0.5, q_frac_mapped)
+
+    q_cart_mapped = q_frac_mapped @ rcell
+
+    return q_cart_mapped, q_frac_mapped
+
 def _doublet_in_list(doublet, llist, nlist):
     """
     Return True if doublet is found in llist[:,:nlist].