From 03cf64cd7a0f9774d474c50439d752a948d75bab Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Tue, 9 Aug 2022 15:55:53 +0200 Subject: [PATCH 1/6] Added new symmetries for the molecule --- cellconstructor/Phonons.py | 7 +- cellconstructor/Structure.py | 54 +++++++++---- cellconstructor/symmetries.py | 145 +++++++++++++++++++++++++--------- 3 files changed, 148 insertions(+), 58 deletions(-) diff --git a/cellconstructor/Phonons.py b/cellconstructor/Phonons.py index 665e972e..ccebd133 100644 --- a/cellconstructor/Phonons.py +++ b/cellconstructor/Phonons.py @@ -1618,7 +1618,7 @@ def ForcePositiveDefinite_2(self): - def GetRamanResponce(self, pol_in, pol_out, T = 0): + def GetRamanResponce(self, pol_in, pol_out, T = 0, w_pols = None): r""" RAMAN RESPONSE ============== @@ -1655,7 +1655,10 @@ def GetRamanResponce(self, pol_in, pol_out, T = 0): if self.raman_tensor is None: raise ValueError("Error, to get the raman responce the raman tensor must be defined") - w, pol_vects = self.DyagDinQ(0) + if w_pols is None: + w, pol_vects = self.DyagDinQ(0) + else: + w, pol_vects = w_pols # Get the mass array _m_ = np.zeros( 3*self.structure.N_atoms) diff --git a/cellconstructor/Structure.py b/cellconstructor/Structure.py index f059cb40..048362b6 100644 --- a/cellconstructor/Structure.py +++ b/cellconstructor/Structure.py @@ -664,7 +664,7 @@ def delete_copies(self, minimum_dist=1e-6, verbose=False): self.coords = np.delete(self.coords, list_pop, axis = 0) self.N_atoms -= N_rep - def apply_symmetry(self, sym_mat, delete_original = False, thr = 1e-6, timer = Timer.Timer()): + def apply_symmetry(self, sym_mat, delete_original = False, thr = 1e-6, avoid_unit_cell = False, timer = Timer.Timer()): """ This function apply the symmetry operation to the atoms of the current structure. @@ -683,37 +683,47 @@ def apply_symmetry(self, sym_mat, delete_original = False, thr = 1e-6, timer = T but bigger than the numerical error in the wyckoff positions of the structure). """ - if not self.has_unit_cell: + if not self.has_unit_cell and not avoid_unit_cell: raise ValueError("The structure has no unit cell!") if delete_original: #self.N_atoms *= 2 new_atoms = np.zeros( (self.N_atoms, 3)) - timer.execute_timed_function(self.fix_coords_in_unit_cell, delete_copies = False) - - old_coords = timer.execute_timed_function(Methods.covariant_coordinates, self.unit_cell, self.coords) + if not avoid_unit_cell: + timer.execute_timed_function(self.fix_coords_in_unit_cell, delete_copies = False) + old_coords = timer.execute_timed_function(Methods.covariant_coordinates, self.unit_cell, self.coords) + else: + old_coords = self.coords.copy() new_coords = sym_mat[:, :3].dot(old_coords.T).T new_coords += np.tile( sym_mat[:, 3], (self.N_atoms, 1)) - self.coords = new_coords.dot(self.unit_cell) - - timer.execute_timed_function(self.fix_coords_in_unit_cell, delete_copies = False) + if not avoid_unit_cell: + self.coords = new_coords.dot(self.unit_cell) + timer.execute_timed_function(self.fix_coords_in_unit_cell, delete_copies = False) + else: + self.coords = new_coords else: for i in range(self.N_atoms): # Convert the coordinates into covariant - old_coords = Methods.covariant_coordinates(self.unit_cell, self.coords[i, :]) + if not avoid_unit_cell: + old_coords = Methods.covariant_coordinates(self.unit_cell, self.coords[i, :]) + else: + old_coords = self.coords[i, :].copy() # Apply the symmetry new_coords = sym_mat[:, :3].dot(old_coords) new_coords += sym_mat[:, 3] # Return into the cartesian coordinates - coords = np.dot( np.transpose(self.unit_cell), new_coords) # Put the atoms into the unit cell - new_atoms[i, :] = Methods.put_into_cell(self.unit_cell, coords) + if not avoid_unit_cell: + coords = np.dot( np.transpose(self.unit_cell), new_coords) + new_atoms[i, :] = Methods.put_into_cell(self.unit_cell, coords) + else: + new_atoms[i, :] = new_coords # Add also the atom type if not delete_original: @@ -745,7 +755,7 @@ def check_symmetry(self, sym_mat, thr = 1e-6): new_struct = self.copy() # Apply the symmetry - new_struct.apply_symmetry(sym_mat, delete_original=True) + new_struct.apply_symmetry(sym_mat, delete_original=True, avoid_unit_cell = not self.has_unit_cell) # Get the equivalence eq_atoms = new_struct.get_equivalent_atoms(self) @@ -754,13 +764,14 @@ def check_symmetry(self, sym_mat, thr = 1e-6): new_struct.coords[eq_atoms, :] = new_struct.coords.copy() # Fix the structure in the unit cell - new_struct.fix_coords_in_unit_cell() + if self.has_unit_cell: + new_struct.fix_coords_in_unit_cell() # Get the displacements u_vect = self.get_displacement(new_struct) # Get the distance between the structures - dist = np.sqrt(np.sum(u_vect ** 2)) + dist = np.sqrt(np.sum(u_vect ** 2)) / self.N_atoms if dist > thr: return False @@ -860,8 +871,9 @@ def impose_symmetries(self, symmetries, threshold = 1.0e-6, verbose = True): for sym in symmetries: aux_struct = self.copy() - aux_struct.apply_symmetry(sym, delete_original = True) - aux_struct.fix_coords_in_unit_cell() + aux_struct.apply_symmetry(sym, delete_original = True, avoid_unit_cell= not self.has_unit_cell) + if self.has_unit_cell: + aux_struct.fix_coords_in_unit_cell() # Get the equivalent atoms eq_atoms = self.get_equivalent_atoms(aux_struct) @@ -937,7 +949,15 @@ def get_equivalent_atoms(self, target_structure, return_distances = False, debug if self.atoms.count(typ) != target_structure.atoms.count(typ): raise ValueError("Error, the target structure must be of the same type of the current one") - eq_atm = list(symph.get_equivalent_atoms(self.coords, target_structure.coords, self.unit_cell, self.get_ityp(), target_structure.get_ityp())) + if self.has_unit_cell: + eq_atm = list(symph.get_equivalent_atoms(self.coords, target_structure.coords, self.unit_cell, self.get_ityp(), target_structure.get_ityp())) + else: + uc_fake = np.zeros((3,3), dtype = np.double) + uc_fake[0,0] = 1000 + uc_fake[1,1] = 1000 + uc_fake[2,2] = 1000 + eq_atm = list(symph.get_equivalent_atoms(self.coords, target_structure.coords, uc_fake, self.get_ityp(), target_structure.get_ityp())) + if debug or return_distances: equiv_atoms = [] diff --git a/cellconstructor/symmetries.py b/cellconstructor/symmetries.py index 7f1139ea..4df3eb26 100644 --- a/cellconstructor/symmetries.py +++ b/cellconstructor/symmetries.py @@ -1216,63 +1216,130 @@ def ApplyTranslationsToVector(self, vector): sum_all /= self.QE_translation_nr vector[:,:] = sum_all + def InitFromSymmetries(self, symmetries): + """ + USE THE GIVEN SYMMETRIES TO SETUP THE SYMMETRIZATION + ====================================== + This function setup all the variables to perform the symmetrization inside the supercell. - - def InitFromSymmetries(self, symmetries, q_point = np.array([0,0,0])): """ - This function initialize the QE symmetries from the symmetries expressed in the - Cellconstructor format, i.e. a list of numpy array 3x4 where the last column is - the fractional translation. + + trans_irt = 0 + self.QE_s[:,:,:] = 0 + + + # Check how many point group symmetries do we have + n_syms = 0 + for i, sym in enumerate(symmetries): + # Extract the rotation and the fractional translation + rot = sym[:,:3] + + # Check if the rotation is equal to the first one + if np.sum( (rot - symmetries[0][:,:3])**2 ) < 0.1 and n_syms == 0 and i > 0: + # We got all the rotations + n_syms = i + break + + # Extract the point group + if n_syms == 0: + self.QE_s[:,:, i] = rot.T + + # Get the IRT (Atoms mapping using symmetries) + irt = GetIRT(self.structure, sym) + self.QE_irt[i, :] = irt + 1 #Py to Fort + - TODO: add the q_point preparation by limitng the symmetries only to - those that satisfies the specified q_point - """ + if n_syms == 0: + n_syms = len(symmetries) - nsym = len(symmetries) + # From the point group symmetries, get the supercell + n_supercell = len(symmetries) // n_syms + self.QE_translation_nr = n_supercell + self.QE_nsymq = n_syms + self.QE_nsym = n_syms + + self.QE_translations_irt = np.zeros( (self.structure.N_atoms, n_supercell), dtype = np.intc, order = "F") + self.QE_translations = np.zeros( (3, n_supercell), dtype = np.double, order = "F") + + # Now extract the translations + for i in range(n_supercell): + sym = symmetries[i * n_syms] + # Check if the symmetries are correctly setup + + I = np.eye(3) + ERROR_MSG=""" + Error, symmetries are not correctly ordered. + They must always start with the identity. + + N_syms = {}; N = {}; SYM = {} + """.format(n_syms,i*n_syms, sym) + assert np.sum( (I - sym[:,:3])**2) < 0.5, ERROR_MSG + + # Get the irt for the translation (and the translation) + irt = GetIRT(self.structure, sym) + self.QE_translations_irt[:, i] = irt + 1 + self.QE_translations[:, i] = sym[:,3] + + # For each symmetry operation, assign the inverse + self.QE_invs[:] = get_invs(self.QE_s, self.QE_nsym) - self.QE_nsymq = np.intc(nsym) - self.QE_nsym = self.QE_nsymq + + + # def InitFromSymmetries(self, symmetries, q_point = np.array([0,0,0])): + # """ + # This function initialize the QE symmetries from the symmetries expressed in the + # Cellconstructor format, i.e. a list of numpy array 3x4 where the last column is + # the fractional translation. + # TODO: add the q_point preparation by limitng the symmetries only to + # those that satisfies the specified q_point + # """ - for i, sym in enumerate(symmetries): - self.QE_s[:,:, i] = np.transpose(sym[:, :3]) + # nsym = len(symmetries) + + # self.QE_nsymq = np.intc(nsym) + # self.QE_nsym = self.QE_nsymq + + + # for i, sym in enumerate(symmetries): + # self.QE_s[:,:, i] = np.transpose(sym[:, :3]) - # Get the atoms correspondence - eq_atoms = GetIRT(self.structure, sym) + # # Get the atoms correspondence + # eq_atoms = GetIRT(self.structure, sym) - self.QE_irt[i, :] = eq_atoms + 1 + # self.QE_irt[i, :] = eq_atoms + 1 - # Get the inverse symmetry - inv_sym = np.linalg.inv(sym[:, :3]) - for k, other_sym in enumerate(symmetries): - if np.sum( (inv_sym - other_sym[:, :3])**2) < __EPSILON__: - break + # # Get the inverse symmetry + # inv_sym = np.linalg.inv(sym[:, :3]) + # for k, other_sym in enumerate(symmetries): + # if np.sum( (inv_sym - other_sym[:, :3])**2) < __EPSILON__: + # break - self.QE_invs[i] = k + 1 + # self.QE_invs[i] = k + 1 - # Setup the position after the symmetry application - for k in range(self.QE_nat): - self.QE_rtau[:, i, k] = self.structure.coords[eq_atoms[k], :].astype(np.float64) + # # Setup the position after the symmetry application + # for k in range(self.QE_nat): + # self.QE_rtau[:, i, k] = self.structure.coords[eq_atoms[k], :].astype(np.float64) - # Get the reciprocal lattice vectors - b_vectors = self.structure.get_reciprocal_vectors() + # # Get the reciprocal lattice vectors + # b_vectors = self.structure.get_reciprocal_vectors() - # Get the minus_q operation - self.QE_minusq = False + # # Get the minus_q operation + # self.QE_minusq = False - # NOTE: HERE THERE COULD BE A BUG + # # NOTE: HERE THERE COULD BE A BUG - # q != -q - # Get the q vectors in crystal coordinates - q = Methods.covariant_coordinates(b_vectors, q_point) - for k, sym in enumerate(self.QE_s): - new_q = self.QE_s[:,:, k].dot(q) - if np.sum( (Methods.put_into_cell(b_vectors, -q_point) - new_q)**2) < __EPSILON__: - self.QE_minus_q = True - self.QE_irotmq = k + 1 - break + # # q != -q + # # Get the q vectors in crystal coordinates + # q = Methods.covariant_coordinates(b_vectors, q_point) + # for k, sym in enumerate(self.QE_s): + # new_q = self.QE_s[:,:, k].dot(q) + # if np.sum( (Methods.put_into_cell(b_vectors, -q_point) - new_q)**2) < __EPSILON__: + # self.QE_minus_q = True + # self.QE_irotmq = k + 1 + # break def GetSymmetries(self, get_irt=False): """ From 010e40e53811eeec82d01d43b9696d984501be14 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Tue, 9 Aug 2022 17:46:57 +0200 Subject: [PATCH 2/6] Now symmetries work --- cellconstructor/Phonons.py | 54 +++++++++++++++++++++++------------ cellconstructor/Structure.py | 9 +++--- cellconstructor/symmetries.py | 9 +++--- 3 files changed, 44 insertions(+), 28 deletions(-) diff --git a/cellconstructor/Phonons.py b/cellconstructor/Phonons.py index ccebd133..ea05c455 100644 --- a/cellconstructor/Phonons.py +++ b/cellconstructor/Phonons.py @@ -2788,6 +2788,8 @@ def GetSupercell(self): supercell : list of 3 int The supercell in each direction. """ + if len(self.dynmats) == 1: + return [1,1,1] return symmetries.GetSupercellFromQlist(self.q_tot, self.structure.unit_cell) def InterpolateMesh(self, mesh_dim, lo_to_splitting = False): @@ -3088,10 +3090,7 @@ def SwapQPoints(self, other_dyn): - - - - def SymmetrizeSupercell(self, supercell_size = None): + def SymmetrizeSupercell(self, sym_mat = None, supercell_size = None): """ Testing function, it applies symmetries in the supercell. """ @@ -3099,34 +3098,47 @@ def SymmetrizeSupercell(self, supercell_size = None): if supercell_size == None: supercell_size = self.GetSupercell() + print("SUPERCELL:", supercell_size) + isgamma = np.prod(supercell_size) == 1 + - if not __SPGLIB__: + if not __SPGLIB__ and sym_mat is None: raise ImportError("Error, the SymmetrizeSupercell method of the Phonon class requires spglib") - superdyn = self.GenerateSupercellDyn(supercell_size) + if not isgamma: + superdyn = self.GenerateSupercellDyn(supercell_size) + super_dynmat = superdyn.dynmats[0] + ss_struct = superdyn.structure + else: + super_dynmat = self.dynmats[0] + ss_struct = self.structure # Apply the sum rule - symmetries.CustomASR(superdyn.dynmats[0]) + symmetries.CustomASR(super_dynmat) + - qe_sym = symmetries.QE_Symmetry(superdyn.structure) - qe_sym.SetupFromSPGLIB() + qe_sym = symmetries.QE_Symmetry(ss_struct) + + if sym_mat is None: + qe_sym.SetupFromSPGLIB() + else: + qe_sym.InitFromSymmetries(sym_mat) #qe_sym.SetupQPoint() - qe_sym.ApplySymmetriesToV2(superdyn.dynmats[0]) + qe_sym.ApplySymmetriesToV2(super_dynmat) #spgsym = spglib.get_symmetry(superdyn.structure.get_ase_atoms()) #syms = symmetries.GetSymmetriesFromSPGLIB(spgsym, False) #superdyn.ForceSymmetries(syms) # Get the dynamical matrix back - fcq = GetDynQFromFCSupercell(superdyn.dynmats[0], np.array(self.q_tot), self.structure, superdyn.structure) + if not isgamma: + fcq = GetDynQFromFCSupercell(superdyn.dynmats[0], np.array(self.q_tot), self.structure, superdyn.structure) - for iq, q in enumerate(self.q_tot): - self.dynmats[iq] = fcq[iq, :, :] + for iq, q in enumerate(self.q_tot): + self.dynmats[iq] = fcq[iq, :, :] # Symmetrize also the effective charges and the Raman Tensor if any # To do this, the symmetries must be initialized once again in the unit cell - qe_sym = symmetries.QE_Symmetry(self.structure) - qe_sym.SetupFromSPGLIB() if not self.effective_charges is None: qe_sym.ApplySymmetryToEffCharge(self.effective_charges) if not self.raman_tensor is None: @@ -3309,7 +3321,7 @@ def ApplySymmetry(self, symmat, irt = None): if irt is None: aux_struct = self.structure.copy() aux_struct.apply_symmetry(symmat, delete_original = True) - aux_struct.fix_coords_in_unit_cell() + #aux_struct.fix_coords_in_unit_cell() eq_atoms = self.structure.get_equivalent_atoms(aux_struct) else: @@ -3336,7 +3348,9 @@ def ApplySymmetry(self, symmat, irt = None): current_m = in_fc[3 * na : 3*na + 3, 3*nb : 3*nb + 3] # Conver the matrix in crystalline - new_m = Methods.convert_matrix_cart_cryst(current_m, self.structure.unit_cell * A_TO_BOHR) + new_m = current_m + if self.structure.has_unit_cell: + new_m = Methods.convert_matrix_cart_cryst(current_m, self.structure.unit_cell * A_TO_BOHR) # Apply the symmetry #new_m_sym = new_s_mat.dot(new_m.dot( new_s_mat.transpose())) @@ -3345,7 +3359,9 @@ def ApplySymmetry(self, symmat, irt = None): #new_m_sym =new_m.copy() # Convert back to cartesian coordinates - new_m = Methods.convert_matrix_cart_cryst(new_m_sym, self.structure.unit_cell * A_TO_BOHR, cryst_to_cart=True) + new_m = new_m_sym + if self.structure.has_unit_cell: + new_m = Methods.convert_matrix_cart_cryst(new_m_sym, self.structure.unit_cell * A_TO_BOHR, cryst_to_cart=True) #print "%d -> %d , %d -> %d)" % (na, s_na, nb, s_nb)#, "d = %.5f" % np.real(np.sqrt(np.sum( (new_m - current_m)**2))) @@ -3391,7 +3407,7 @@ def ForceSymmetries(self, symmetries, irt = None, apply_sum_rule = True): new_fc = np.zeros( np.shape(self.dynmats[0]), dtype = np.complex128 ) - self.structure.fix_coords_in_unit_cell() + #self.structure.fix_coords_in_unit_cell() for i, sym in enumerate(symmetries): # Check if the structure satisfy the symmetry if not self.structure.check_symmetry(sym): diff --git a/cellconstructor/Structure.py b/cellconstructor/Structure.py index 048362b6..f6d04fec 100644 --- a/cellconstructor/Structure.py +++ b/cellconstructor/Structure.py @@ -682,10 +682,9 @@ def apply_symmetry(self, sym_mat, delete_original = False, thr = 1e-6, avoid_uni (must be smaller than the minimum distance between two generic atoms in the struct, but bigger than the numerical error in the wyckoff positions of the structure). """ - - if not self.has_unit_cell and not avoid_unit_cell: - raise ValueError("The structure has no unit cell!") - + if not self.has_unit_cell: + avoid_unit_cell = True + if delete_original: #self.N_atoms *= 2 new_atoms = np.zeros( (self.N_atoms, 3)) @@ -772,7 +771,7 @@ def check_symmetry(self, sym_mat, thr = 1e-6): # Get the distance between the structures dist = np.sqrt(np.sum(u_vect ** 2)) / self.N_atoms - + print("DIST:", dist) if dist > thr: return False return True diff --git a/cellconstructor/symmetries.py b/cellconstructor/symmetries.py index 4df3eb26..822fe8ae 100644 --- a/cellconstructor/symmetries.py +++ b/cellconstructor/symmetries.py @@ -1841,10 +1841,11 @@ def GetIRT(structure, symmetry, timer = Timer.Timer(), debug = False): new_struct = structure.copy() - if timer is None: - new_struct.fix_coords_in_unit_cell(delete_copies = False, debug = debug) - else: - timer.execute_timed_function(new_struct.fix_coords_in_unit_cell, delete_copies = False, debug = debug) + if new_struct.has_unit_cell: + if timer is None: + new_struct.fix_coords_in_unit_cell(delete_copies = False, debug = debug) + else: + timer.execute_timed_function(new_struct.fix_coords_in_unit_cell, delete_copies = False, debug = debug) n_struct_2 = new_struct.copy() if timer is None: From cbbd581f77ec856620bc14d1c9645233375b3fe5 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Tue, 9 Aug 2022 19:41:19 +0200 Subject: [PATCH 3/6] Bugfix on the default value for bg --- cellconstructor/Phonons.py | 11 ++++++++++- cellconstructor/Structure.py | 8 ++++++++ cellconstructor/symmetries.py | 22 ++++++++++++++-------- 3 files changed, 32 insertions(+), 9 deletions(-) diff --git a/cellconstructor/Phonons.py b/cellconstructor/Phonons.py index ea05c455..20592d3b 100644 --- a/cellconstructor/Phonons.py +++ b/cellconstructor/Phonons.py @@ -3489,6 +3489,9 @@ def DiagonalizeSupercell(self, verbose = False, lo_to_split = None): supercell_size = len(self.q_tot) nat = self.structure.N_atoms + if not self.structure.has_unit_cell and supercell_size > 1: + raise ValueError("Error, the structure must have a defined supercell") + nmodes = 3*nat*supercell_size nat_sc = nat*supercell_size @@ -3510,7 +3513,13 @@ def DiagonalizeSupercell(self, verbose = False, lo_to_split = None): R_vec[3*i : 3*i+3, :] = np.tile(super_structure.coords[i, :] - self.structure.coords[itau[i], :], (3,1)) i_mu = 0 - bg = self.structure.get_reciprocal_vectors() / (2*np.pi) + + if self.structure.has_unit_cell: + bg = self.structure.get_reciprocal_vectors() / (2*np.pi) + else: + bg = np.eye(3) + + for iq, q in enumerate(self.q_tot): # Check if the current q point has been seen (we do not distinguish between q and -q) skip_this_q = False diff --git a/cellconstructor/Structure.py b/cellconstructor/Structure.py index f6d04fec..de5b0344 100644 --- a/cellconstructor/Structure.py +++ b/cellconstructor/Structure.py @@ -1519,6 +1519,14 @@ def generate_supercell(self, dim, itau = None, QE_convention = True, get_itau = if len(dim) != 3: raise ValueError("ERROR, dim must have 3 integers.") + if np.prod(dim) == 1: + if not get_itau: + return self.copy() + + itau = np.zeros(self.N_atoms, dtype = np.intc) + itau[:] = np.arange(self.N_atoms) + return self.copy(), itau + if not self.has_unit_cell: raise ValueError("ERROR, the specified system has not the unit cell.") diff --git a/cellconstructor/symmetries.py b/cellconstructor/symmetries.py index 822fe8ae..777e2ba9 100644 --- a/cellconstructor/symmetries.py +++ b/cellconstructor/symmetries.py @@ -68,7 +68,8 @@ def __init__(self, structure, threshold = 1e-5): """ if not structure.has_unit_cell: - raise ValueError("Error, symmetry operation can be initialize only if the structure has a unit cell") + warnings.warn("WARNING: the structure has no unit cell!!!!") + # raise ValueError("Error, symmetry operation can be initialize only if the structure has a unit cell") self.structure = structure self.threshold = np.float64(threshold) @@ -114,12 +115,18 @@ def __init__(self, structure, threshold = 1e-5): self.QE_at = np.zeros( (3,3), dtype = np.float64, order = "F") self.QE_bg = np.zeros( (3,3), dtype = np.float64, order = "F") - bg = structure.get_reciprocal_vectors() - for i in range(3): - for j in range(3): - self.QE_at[i,j] = structure.unit_cell[j,i] - self.QE_bg[i,j] = bg[j,i] / (2* np.pi) - + if structure.has_unit_cell: + bg = structure.get_reciprocal_vectors() + for i in range(3): + for j in range(3): + self.QE_at[i,j] = structure.unit_cell[j,i] + self.QE_bg[i,j] = bg[j,i] / (2* np.pi) + else: + bg = np.eye(3) + for i in range(3): + self.QE_at[i,i] = 1 + self.QE_bg[i,i] = 1 + # Here we define the quantities required to symmetrize the supercells self.QE_at_sc = self.QE_at.copy() self.QE_bg_sc = self.QE_bg.copy() @@ -129,7 +136,6 @@ def __init__(self, structure, threshold = 1e-5): # After the translation, which vector is transformed in which one? # This info is stored here as ndarray( size = (N_atoms, N_trans), dtype = np.intc, order = "F") self.QE_translations_irt = [] - def ForceSymmetry(self, structure): """ FORCE SYMMETRY From 1d708441c27b5a289e0aed24d351728a099efc82 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 10 Aug 2022 17:32:11 +0200 Subject: [PATCH 4/6] Now works --- FModules/symvector.f90 | 57 +++++++++++++++++++++++++++++++++++ cellconstructor/Phonons.py | 47 ++++++++++++++++++----------- cellconstructor/Structure.py | 2 +- cellconstructor/symmetries.py | 51 ++++++++++++++++++++++++++++--- 4 files changed, 135 insertions(+), 22 deletions(-) diff --git a/FModules/symvector.f90 b/FModules/symvector.f90 index 9fd45749..b6ecda81 100644 --- a/FModules/symvector.f90 +++ b/FModules/symvector.f90 @@ -51,3 +51,60 @@ SUBROUTINE symvector (nat, nsym, irt, s, at, bg, vect) DEALLOCATE (work) ! END SUBROUTINE symvector + + + + + SUBROUTINE symvector_double (nat, nsym, irt, s, at, bg, vect) + !----------------------------------------------------------------------- + ! Symmetrize a function f(i,na), i=cartesian component, na=atom index + ! e.g. : forces (in cartesian axis) + ! + IMPLICIT NONE + ! + INTEGER, INTENT(IN) :: nat, nsym + INTEGER, INTENT(IN) :: irt(48,nat) + double precision, INTENT(IN) :: s(3,3,48) + double precision, INTENT(IN) :: at(3,3), bg(3,3) + double precision, intent(INOUT) :: vect(3,nat) + ! + INTEGER :: na, isym, nar + double precision, ALLOCATABLE :: work (:,:) + ! + IF (nsym == 1) RETURN + ! + ALLOCATE (work(3,nat)) + ! + ! bring vector to crystal axis + ! + DO na = 1, nat + work(:,na) = vect(1,na)*at(1,:) + & + vect(2,na)*at(2,:) + & + vect(3,na)*at(3,:) + END DO + ! + ! symmetrize in crystal axis + ! + vect (:,:) = 0.0d0 + DO na = 1, nat + DO isym = 1, nsym + nar = irt (isym, na) + vect (:, na) = vect (:, na) + & + s (:, 1, isym) * work (1, nar) + & + s (:, 2, isym) * work (2, nar) + & + s (:, 3, isym) * work (3, nar) + END DO + END DO + work (:,:) = vect (:,:) / DBLE(nsym) + ! + ! bring vector back to cartesian axis + ! + DO na = 1, nat + vect(:,na) = work(1,na)*bg(:,1) + & + work(2,na)*bg(:,2) + & + work(3,na)*bg(:,3) + END DO + ! + DEALLOCATE (work) + ! + END SUBROUTINE symvector_double diff --git a/cellconstructor/Phonons.py b/cellconstructor/Phonons.py index 20592d3b..94e18c46 100644 --- a/cellconstructor/Phonons.py +++ b/cellconstructor/Phonons.py @@ -3101,6 +3101,10 @@ def SymmetrizeSupercell(self, sym_mat = None, supercell_size = None): print("SUPERCELL:", supercell_size) isgamma = np.prod(supercell_size) == 1 + if sym_mat is not None: + self.ForceSymmetries(sym_mat) + return + if not __SPGLIB__ and sym_mat is None: raise ImportError("Error, the SymmetrizeSupercell method of the Phonon class requires spglib") @@ -3119,12 +3123,10 @@ def SymmetrizeSupercell(self, sym_mat = None, supercell_size = None): qe_sym = symmetries.QE_Symmetry(ss_struct) - if sym_mat is None: - qe_sym.SetupFromSPGLIB() - else: - qe_sym.InitFromSymmetries(sym_mat) - #qe_sym.SetupQPoint() + qe_sym.SetupFromSPGLIB() qe_sym.ApplySymmetriesToV2(super_dynmat) + + #qe_sym.SetupQPoint() #spgsym = spglib.get_symmetry(superdyn.structure.get_ase_atoms()) #syms = symmetries.GetSymmetriesFromSPGLIB(spgsym, False) @@ -3219,11 +3221,11 @@ def ApplySumRule(self, kind = "custom"): # Apply the sum rule on the effective charge - if self.effective_charges is not None: - total_charge = np.sum(self.effective_charges, axis = 0) + #if self.effective_charges is not None: + # total_charge = np.sum(self.effective_charges, axis = 0) - # Subtract to each atom an average of the total charges - self.effective_charges = np.einsum("aij, ij -> aij", self.effective_charges, - total_charge / self.structure.N_atoms) + # # Subtract to each atom an average of the total charges + # self.effective_charges = np.einsum("aij, ij -> aij", self.effective_charges, - total_charge / self.structure.N_atoms) def GetIRActive(self, use_spglib = False): @@ -3326,6 +3328,8 @@ def ApplySymmetry(self, symmat, irt = None): eq_atoms = self.structure.get_equivalent_atoms(aux_struct) else: eq_atoms = irt + + #print(" IRT AP sym: ", eq_atoms) #print eq_atoms # Get the number of atoms @@ -3378,7 +3382,7 @@ def ApplySymmetry(self, symmat, irt = None): - def ForceSymmetries(self, symmetries, irt = None, apply_sum_rule = True): + def ForceSymmetries(self, syms, irt = None, apply_sum_rule = True): """ FORCE THE PHONON TO RESPECT THE SYMMETRIES ========================================== @@ -3406,9 +3410,9 @@ def ForceSymmetries(self, symmetries, irt = None, apply_sum_rule = True): # Apply the symmetries new_fc = np.zeros( np.shape(self.dynmats[0]), dtype = np.complex128 ) - + print("EF CH HERE:", self.effective_charges) #self.structure.fix_coords_in_unit_cell() - for i, sym in enumerate(symmetries): + for i, sym in enumerate(syms): # Check if the structure satisfy the symmetry if not self.structure.check_symmetry(sym): print (sym) @@ -3421,9 +3425,9 @@ def ForceSymmetries(self, symmetries, irt = None, apply_sum_rule = True): current_irt = None if not irt is None: current_irt = irt[i, :] + #print("I: ", i, end = "") current_fc = self.ApplySymmetry(sym, irt = current_irt) - print (i) # Try to add the sum rule here #newP = self.Copy() @@ -3438,10 +3442,12 @@ def ForceSymmetries(self, symmetries, irt = None, apply_sum_rule = True): new_fc += current_fc # Average all the symmetrized structures - new_fc /= len(symmetries) + new_fc /= len(syms) - print ("DIST_SYM_FORC:", np.sqrt(np.sum( (new_fc - self.dynmats[0])**2))) + #print ("DIST_SYM_FORC:", np.sqrt(np.sum( (new_fc - self.dynmats[0])**2))) + if apply_sum_rule: + symmetries.CustomASR(new_fc) self.dynmats[0] = new_fc.copy() @@ -3449,8 +3455,15 @@ def ForceSymmetries(self, symmetries, irt = None, apply_sum_rule = True): #print "\n".join( ["\t".join("%.4e" % (xval - freqs[0,j]) for xval in freqs[:, j]) for j in range(3 * self.structure.N_atoms)]) # Apply the acustic sum rule - if apply_sum_rule: - self.ApplySumRule() + + qe_sym = symmetries.QE_Symmetry(self.structure) + qe_sym.InitFromSymmetries(syms) + + print("EF CH HERE2:", self.effective_charges) + if not self.effective_charges is None: + qe_sym.ApplySymmetryToEffCharge(self.effective_charges) + if not self.raman_tensor is None: + qe_sym.ApplySymmetryToRamanTensor(self.raman_tensor) def DiagonalizeSupercell(self, verbose = False, lo_to_split = None): r""" diff --git a/cellconstructor/Structure.py b/cellconstructor/Structure.py index de5b0344..cfc3c16e 100644 --- a/cellconstructor/Structure.py +++ b/cellconstructor/Structure.py @@ -771,7 +771,7 @@ def check_symmetry(self, sym_mat, thr = 1e-6): # Get the distance between the structures dist = np.sqrt(np.sum(u_vect ** 2)) / self.N_atoms - print("DIST:", dist) + #print("DIST:", dist) if dist > thr: return False return True diff --git a/cellconstructor/symmetries.py b/cellconstructor/symmetries.py index 777e2ba9..479722dd 100644 --- a/cellconstructor/symmetries.py +++ b/cellconstructor/symmetries.py @@ -390,6 +390,7 @@ def ApplySymmetryToEffCharge(self, eff_charges): assert nat == self.QE_nat, "Error, the structure and effective charges are not compatible" + # Apply the sum rule tot_sum = np.sum(eff_charges, axis = 0) eff_charges -= np.tile(tot_sum, (nat, 1)).reshape((nat, 3,3 )) / nat @@ -397,9 +398,21 @@ def ApplySymmetryToEffCharge(self, eff_charges): new_eff_charges = np.zeros((nat, cart1, cart2), dtype = np.double) # Get the effective charges in crystal components + print("QE_at:", self.QE_at) + print("EC before:") + print(eff_charges[:, :, :]) for i in range(nat): eff_charges[i, :, :] = Methods.convert_matrix_cart_cryst(eff_charges[i, :, :], self.QE_at.T) + + print("SYmmetries:") + for i in range(self.QE_nsym): + print("I:", i) + print(self.QE_s[:, :, i].T) + print() + + print("TRANS:", self.QE_translation_nr) + # Apply translations if self.QE_translation_nr > 1: for i in range(self.QE_translation_nr): @@ -417,6 +430,9 @@ def ApplySymmetryToEffCharge(self, eff_charges): for j in range(nat): new_mat = self.QE_s[:,:, i].dot( eff_charges[irt[j], :, :].dot(self.QE_s[:,:,i].T)) + print("SYM {} AT {}, new effective charge:".format(i, j)) + print(new_mat) + print("OLD EFFECTIVE CHARGE:", eff_charges[irt[j], :, :]) new_eff_charges[j, :, :] += new_mat new_eff_charges /= self.QE_nsym @@ -1232,7 +1248,9 @@ def InitFromSymmetries(self, symmetries): """ trans_irt = 0 - self.QE_s[:,:,:] = 0 + warnings.warn("NOTE: use this initialization only for Effective Charges and Raman") + self.QE_s = np.zeros((3, 3, 48), dtype = np.double) + #self.QE_s[:,:,:] = 0 # Check how many point group symmetries do we have @@ -1412,9 +1430,24 @@ def SymmetrizeVector(self, vector): tmp_vector[0, i] = vector[i,0] tmp_vector[1, i] = vector[i,1] tmp_vector[2,i] = vector[i,2] - - symph.symvector(self.QE_nsymq, self.QE_irt, self.QE_s, self.QE_at, self.QE_bg, - tmp_vector, self.QE_nat) + + + print("Symmetries (total: {})".format(self.QE_nsym)) + for i in range(self.QE_nsym): + print("I:", i) + print(self.QE_s[:, :, i].T) + print() + + if isinstance(self.QE_s[0,0,0], np.intc): + symph.symvector(self.QE_nsymq, self.QE_irt, self.QE_s, self.QE_at, self.QE_bg, + tmp_vector, self.QE_nat) + elif isinstance(self.QE_s[0,0,0], np.double): + symph.symvector_double(self.QE_nsymq, self.QE_irt, self.QE_s, self.QE_at, self.QE_bg, + tmp_vector, self.QE_nat) + else: + print("TYPE:", type(self.QE_s[0,0,0])) + raise ValueError("Error, type of QE_s not recognized. Error while initializing the symmetries") + for i in range(self.QE_nat): @@ -1607,6 +1640,8 @@ def ApplySymmetriesToV2(self, v2, apply_translations = True): # Apply the Permutation symmetry v2[:,:] = 0.5 * (v2 + v2.T) + assert isinstance(self.QE_s[0,0,0], np.intc), "Error, initialize the symmetries correctly" + # First lets recall that the fortran subroutines # Takes the input as (3,3,nat,nat) new_v2 = np.zeros( (3,3, self.QE_nat, self.QE_nat), dtype = np.double, order ="F") @@ -1621,6 +1656,14 @@ def ApplySymmetriesToV2(self, v2, apply_translations = True): symph.trans_v2(new_v2, self.QE_translations_irt) # Apply the symmetrization + #print("QE IRT:", self.QE_irt[:self.QE_nsym, :]-1) + #print("QE BG:", self.QE_bg) + #print("QE AT:", self.QE_at) + #print("QE ALL SYMMETRIES:") + #for i in range(self.QE_nsym): + # print("{})".format(i)) + # print(self.QE_s[:, :, i].T) + # print() symph.sym_v2(new_v2, self.QE_at, self.QE_bg, self.QE_s, self.QE_irt, self.QE_nsym, self.QE_nat) # Return back From 1ed4a680119bd7edbdfd95ae743a7161501a6203 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 17 Nov 2022 10:06:33 +0100 Subject: [PATCH 5/6] Fix a problem with the unit cell --- cellconstructor/Phonons.py | 10 +++++++++- setup.py | 4 ++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/cellconstructor/Phonons.py b/cellconstructor/Phonons.py index 94e18c46..348fd44f 100644 --- a/cellconstructor/Phonons.py +++ b/cellconstructor/Phonons.py @@ -4357,7 +4357,15 @@ def compute_phonons_finite_displacements(structure, ase_calculator, epsilon = 0. """ - super_structure = structure.generate_supercell(supercell) + if np.prod(supercell) > 1: + super_structure = structure.generate_supercell(supercell) + else: + super_structure = structure.copy() + + if not super_structure.has_unit_cell: + super_structure.unit_cell = np.eye(3) * 200 + super_structure.has_unit_cell = True + final_dyn = Phonons(super_structure) nat3 = 3 * super_structure.N_atoms diff --git a/setup.py b/setup.py index 6bab55d1..0c0678c5 100644 --- a/setup.py +++ b/setup.py @@ -61,14 +61,14 @@ setup( name = "CellConstructor", - version = "1.1", + version = "1.2", description = "Python utilities that is interfaced with ASE for atomic crystal analysis", author = "Lorenzo Monacelli", url = "https://github.com/mesonepigreco/CellConstructor", packages = ["cellconstructor"], package_dir = {"cellconstructor": "cellconstructor"}, package_data = {"cellconstructor": ["SymData/*.dat"]}, - setup_requires = ["numpy", "ase", "scipy"], + install_requires = ["", "", "", ""], license = "MIT", include_package_data = True, scripts = ["scripts/symmetrize_dynmat.py", "scripts/cellconstructor_test.py", "scripts/view_scf_atoms.py"], From a1e57bdfa5df801dc4ef352ebd39e6988ac13311 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 10 Aug 2023 11:53:11 +0200 Subject: [PATCH 6/6] Fixed the encoding --- FModules/get_equivalent_atoms.f90 | 3 +-- setup.py | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/FModules/get_equivalent_atoms.f90 b/FModules/get_equivalent_atoms.f90 index 28befeb3..fe5796a3 100644 --- a/FModules/get_equivalent_atoms.f90 +++ b/FModules/get_equivalent_atoms.f90 @@ -199,7 +199,6 @@ end subroutine fix_coords_in_unit_cell subroutine matinv3(A, B) implicit none - !! Performs a direct calculation of the inverse of a 3×3 matrix. double precision, intent(in) :: A(3,3) !! Matrix double precision, intent(out) :: B(3,3) !! Inverse matrix double precision :: detinv @@ -219,4 +218,4 @@ subroutine matinv3(A, B) B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2)) B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1)) B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1)) -end subroutine matinv3 \ No newline at end of file +end subroutine matinv3 diff --git a/setup.py b/setup.py index 0c0678c5..d0dc98a7 100644 --- a/setup.py +++ b/setup.py @@ -68,7 +68,6 @@ packages = ["cellconstructor"], package_dir = {"cellconstructor": "cellconstructor"}, package_data = {"cellconstructor": ["SymData/*.dat"]}, - install_requires = ["", "", "", ""], license = "MIT", include_package_data = True, scripts = ["scripts/symmetrize_dynmat.py", "scripts/cellconstructor_test.py", "scripts/view_scf_atoms.py"],