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/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 5844307a..5ef8d8e4 100644 --- a/cellconstructor/Phonons.py +++ b/cellconstructor/Phonons.py @@ -1618,9 +1618,8 @@ def ForcePositiveDefinite_2(self): matrix = np.einsum("i, ji, ki", fr**2, v, np.conj(v)) * np.sqrt(_m1_ * _m2_) self.dynmats[iq] = matrix - - - def GetRamanResponce(self, pol_in, pol_out, T = 0): + + def GetRamanResponce(self, pol_in, pol_out, T = 0, w_pols = None): r""" RAMAN RESPONSE ============== @@ -1657,8 +1656,11 @@ 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) for i in range(self.structure.N_atoms): @@ -2810,6 +2812,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): @@ -3111,9 +3115,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. """ @@ -3121,17 +3123,34 @@ 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__: - raise ImportError("Error, the SymmetrizeSupercell method of the Phonon class requires spglib") + if sym_mat is not None: + self.ForceSymmetries(sym_mat) + return - superdyn = self.GenerateSupercellDyn(supercell_size) + + if not __SPGLIB__ and sym_mat is None: + raise ImportError("Error, the SymmetrizeSupercell method of the Phonon class requires spglib") + + 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(ss_struct) - qe_sym = symmetries.QE_Symmetry(superdyn.structure) qe_sym.SetupFromSPGLIB() + qe_sym.ApplySymmetriesToV2(super_dynmat) + #qe_sym.SetupQPoint() qe_sym.ApplySymmetriesToV2(superdyn.dynmats[0]) @@ -3140,15 +3159,14 @@ def SymmetrizeSupercell(self, supercell_size = None): #superdyn.ForceSymmetries(syms) # Get the dynamical matrix back - fcq = GetDynQFromFCSupercell_parallel(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, :, :] + 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, :, :] # 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: @@ -3229,12 +3247,12 @@ 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) - - # 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) + #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) + def GetIRActive(self, use_spglib = False): """ @@ -3331,11 +3349,13 @@ 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: eq_atoms = irt + + #print(" IRT AP sym: ", eq_atoms) #print eq_atoms # Get the number of atoms @@ -3358,8 +3378,10 @@ 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())) new_m_sym = new_s_mat.transpose().dot(new_m.dot( new_s_mat)) @@ -3367,8 +3389,10 @@ 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))) # Write the matrix into the output @@ -3381,10 +3405,10 @@ def ApplySymmetry(self, symmat, irt = None): # Return the symmetrized result #print "Total distance:", np.sqrt(np.sum( (out_fc - np.real(in_fc))**2)) return out_fc - - - - 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 ========================================== @@ -3412,9 +3436,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 ) - - self.structure.fix_coords_in_unit_cell() - for i, sym in enumerate(symmetries): + print("EF CH HERE:", self.effective_charges) + #self.structure.fix_coords_in_unit_cell() + for i, sym in enumerate(syms): # Check if the structure satisfy the symmetry if not self.structure.check_symmetry(sym): print (sym) @@ -3427,10 +3451,8 @@ 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() #newP.dynmats[0] = current_fc @@ -3444,10 +3466,12 @@ def ForceSymmetries(self, symmetries, irt = None, apply_sum_rule = True): new_fc += current_fc # Average all the symmetrized structures - new_fc /= len(symmetries) - - - print ("DIST_SYM_FORC:", np.sqrt(np.sum( (new_fc - self.dynmats[0])**2))) + new_fc /= len(syms) + + + #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() @@ -3455,8 +3479,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""" @@ -3495,6 +3526,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 @@ -3516,7 +3550,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 @@ -4421,7 +4461,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/cellconstructor/Structure.py b/cellconstructor/Structure.py index ac4313c5..6ee8663a 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. @@ -682,38 +682,47 @@ def apply_symmetry(self, sym_mat, delete_original = False, thr = 1e-6, timer = T (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: - raise ValueError("The structure has no unit cell!") - + avoid_unit_cell = True + 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 +754,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,14 +763,15 @@ 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 + #print("DIST:", dist) if dist > thr: return False return True @@ -860,8 +870,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 +948,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 = [] @@ -1500,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 920352bf..946b1d9d 100644 --- a/cellconstructor/symmetries.py +++ b/cellconstructor/symmetries.py @@ -69,7 +69,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) @@ -115,12 +116,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() @@ -130,7 +137,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 @@ -385,6 +391,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 @@ -392,9 +399,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): @@ -412,6 +431,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 @@ -1217,63 +1239,132 @@ 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 + 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 + 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) + + # 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) + + + + # 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. - nsym = len(symmetries) + # TODO: add the q_point preparation by limitng the symmetries only to + # those that satisfies the specified q_point + # """ - self.QE_nsymq = np.intc(nsym) - self.QE_nsym = self.QE_nsymq + # 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]) + + # 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): """ @@ -1340,9 +1431,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): @@ -1535,6 +1641,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") @@ -1549,6 +1657,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 @@ -1775,10 +1891,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: diff --git a/setup.py b/setup.py index 6bab55d1..d0dc98a7 100644 --- a/setup.py +++ b/setup.py @@ -61,14 +61,13 @@ 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"], license = "MIT", include_package_data = True, scripts = ["scripts/symmetrize_dynmat.py", "scripts/cellconstructor_test.py", "scripts/view_scf_atoms.py"],