Skip to content

Molecular symmetries #57

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.

Sign up for GitHub

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

Open
mesonepigreco wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Molecular symmetries #57

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
03cf64c
Added new symmetries for the molecule
mesonepigreco Aug 9, 2022
010e40e
Now symmetries work
mesonepigreco Aug 9, 2022
cbbd581
Bugfix on the default value for bg
mesonepigreco Aug 9, 2022
1d70844
Now works
mesonepigreco Aug 10, 2022
1ed4a68
Fix a problem with the unit cell
mesonepigreco Nov 17, 2022
fe4ee26
Merge branch 'master' into molecular_symmetries
mesonepigreco Nov 17, 2022
a1e57bd
Fixed the encoding
mesonepigreco Aug 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions FModules/get_equivalent_atoms.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
end subroutine matinv3
57 changes: 57 additions & 0 deletions FModules/symvector.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
142 changes: 95 additions & 47 deletions cellconstructor/Phonons.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
==============
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -3111,27 +3115,42 @@ 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.
"""

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])

Expand All @@ -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:
Expand Down Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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
Expand All @@ -3358,17 +3378,21 @@ 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))

#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
Expand All @@ -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
==========================================
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -3444,19 +3466,28 @@ 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()


# Print the phonons all toghether
#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"""
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading