diff --git a/FModules/fc_supercell_from_dyn.f90 b/FModules/fc_supercell_from_dyn.f90 index d94cb126..752cf39b 100644 --- a/FModules/fc_supercell_from_dyn.f90 +++ b/FModules/fc_supercell_from_dyn.f90 @@ -81,7 +81,7 @@ subroutine fast_ft_real_space_from_dynq(unit_cell_coords, super_cell_coords, ita double complex :: im, phase im = (0.0d0,1.0d0) - twopi = 6.283185307179586d0 + twopi = 6.283185307179586d0 fc_supercell(:,:) = 0.0d0 diff --git a/cellconstructor/Methods.py b/cellconstructor/Methods.py index 635f3330..adc1c13c 100644 --- a/cellconstructor/Methods.py +++ b/cellconstructor/Methods.py @@ -1422,6 +1422,38 @@ def convert_3tensor_to_cryst(tensor, unit_cell, cryst_to_cart = False): comp_matrix = np.linalg.inv(comp_matrix) return np.einsum("ia, jb, kc, ijk -> abc", comp_matrix, comp_matrix, comp_matrix, tensor) + + +def convert_4tensor_to_cryst(tensor, unit_cell, cryst_to_cart = False): + """ + Convert to crystal coordinates + ============================== + + This subroutine converts the 4 rank tensor to crystal coordinates and vice versa. + + Paramters + --------- + tensor : ndarray( size = (3,3,3,3)) + The 4rank tensor to be converted + unit_cell : ndarray + The unit cell of the structure + cryst_to_cart : bool + If true, reverse convert crystal to cartesian. + """ + + # Get the metric tensor from the unit_cell + metric_tensor = np.zeros((3,3)) + for i in range(0, 3): + for j in range(i, 3): + metric_tensor[i, j] = metric_tensor[j,i] = unit_cell[i,:].dot(unit_cell[j, :]) + + comp_matrix = np.einsum("ij, jk", np.linalg.inv(metric_tensor), unit_cell) + if not cryst_to_cart: + comp_matrix = np.linalg.inv(comp_matrix) + + return np.einsum("ia, jb, kc, ld, ijkl -> abcd", comp_matrix, comp_matrix, comp_matrix, comp_matrix, tensor) + + def convert_fc(fc_matrix, unit_cell, cryst_to_cart = False): """ diff --git a/cellconstructor/Phonons.py b/cellconstructor/Phonons.py index 33f6921e..92644b8c 100644 --- a/cellconstructor/Phonons.py +++ b/cellconstructor/Phonons.py @@ -4102,6 +4102,7 @@ def GetSupercellFCFromDyn(dynmat, q_tot, unit_cell_structure, supercell_structur itau : Ndarray(nat_sc) , optional the correspondance between the supercell atoms and the unit cell one. If None is recomputed + img_thr: 1e-5 def Returns ------- fc_supercell : ndarray 3nat_sc x 3nat_sc @@ -4164,9 +4165,9 @@ def GetSupercellFCFromDyn(dynmat, q_tot, unit_cell_structure, supercell_structur imag = np.max(np.abs(np.imag(fc))) ASSERT_ERROR = """ Error, the imaginary part of the real space force constant - is not zero. IMAG={} + is not zero. IMAG={} imthr={} """ - assert imag < img_thr, ASSERT_ERROR.format(imag) + assert imag < img_thr, ASSERT_ERROR.format(imag, img_thr) # Remove anyway the imaginary part return fc - 1j*np.imag(fc) @@ -4377,8 +4378,10 @@ def InterpolateDynFC(starting_fc, coarse_grid, unit_cell_structure, super_cell_s natsc = np.shape(starting_fc)[0] // 3 nat = natsc // supercell_size - #print "nat:", nat - #print "natsc:", natsc + print('shape',starting_fc.shape) + print('sc',natsc) + print('uc',nat) + # Get the force constant in an appropriate supercell diff --git a/cellconstructor/symmetries.py b/cellconstructor/symmetries.py index fda3801e..d6fe10d8 100644 --- a/cellconstructor/symmetries.py +++ b/cellconstructor/symmetries.py @@ -477,7 +477,7 @@ def ApplySymmetryToRamanTensor(self, raman_tensor): Parameters ---------- - - raman_tensor : ndarray (size = (3, 3, 3*nat)) + - raman_tensor : ndarray (size = (3, 3, 3 * nat)) The raman tensor. The first two indices indicate the polarization of the incoming/outcoming field, while the last one is the atomic/cartesian coordinate @@ -509,6 +509,8 @@ def ApplySymmetryToRamanTensor(self, raman_tensor): # Now we apply the sum rule rt_reshaped -= shift / self.QE_nat + + # Auxiliary variable new_tensor = np.zeros(np.shape(rt_reshaped), dtype = np.double) # Get the raman tensor in crystal components @@ -521,7 +523,7 @@ def ApplySymmetryToRamanTensor(self, raman_tensor): irt = self.QE_translations_irt[:, i] - 1 for j in range(self.QE_nat): new_mat = rt_reshaped[:,:, irt[j], :] - new_tensor += new_mat + new_tensor[:,:,j,:] += new_mat[:,:,:] rt_reshaped = new_tensor / self.QE_translation_nr new_tensor[:,:,:,:] = 0. @@ -532,7 +534,7 @@ def ApplySymmetryToRamanTensor(self, raman_tensor): for j in range(self.QE_nat): # Apply the symmetry to the 3 order tensor - new_mat = np.einsum("ai, bj, ck, ijk", self.QE_s[:,:,i], self.QE_s[:,:,i], self.QE_s[:,:,i], rt_reshaped[:,:, irt[j], :]) + new_mat = np.einsum("ai, bj, ck, ijk -> abc", self.QE_s[:,:,i], self.QE_s[:,:,i], self.QE_s[:,:,i], rt_reshaped[:,:, irt[j], :]) #new_mat = self.QE_s[:,:, i].dot( eff_charges[irt[j], :, :].dot(self.QE_s[:,:,i].T)) new_tensor[:,:,j,:] += new_mat @@ -559,8 +561,9 @@ def ApplySymmetryToSecondOrderEffCharge(self, dM_drdr, apply_asr = True): Parameters ---------- - dM_drdr : ndarray (size = (3 nat_sc, 3nat_sc, 3)) + dM_drdr : ndarray (size = (3 nat_sc, 3 nat_sc, 3)) The derivative of effective charges. + The last index refers to electric field apply_asr : bool If True the sum rule is applied. The sum rule is the 'custom' one where translations are projected @@ -574,12 +577,8 @@ def ApplySymmetryToSecondOrderEffCharge(self, dM_drdr, apply_asr = True): assert cart == 3 nat = int(nat3 / 3) - - # Apply hermitianity - #print("Original:") - #print(dM_drdr[:,:,0]) - dM_drdr += np.einsum("abc->bac", dM_drdr) + dM_drdr += np.einsum("abc -> bac", dM_drdr) dM_drdr /= 2 # Apply the Sum Rule @@ -587,19 +586,12 @@ def ApplySymmetryToSecondOrderEffCharge(self, dM_drdr, apply_asr = True): for pol in range(3): CustomASR(dM_drdr[:,:,pol]) - #print("After the sum rule:") - #print(dM_drdr[:,:,0]) - - # Convert in crystal coordinates + # Convert in crystal coordinates MAKE SURE THAT THIS IS CONVERTED BACK IN CARTESIAN COORDINATES for i in range(nat): for j in range(nat): dM_drdr[3*i : 3*i + 3, 3*j: 3*j+3, :] = Methods.convert_3tensor_to_cryst(dM_drdr[3*i:3*i+3, 3*j:3*j+3,:], self.QE_at.T) - #print("Crystal:") - #print(dM_drdr[:,:,0]) - - # Apply translations new_dM = np.zeros(np.shape(dM_drdr), dtype = np.double) if self.QE_translation_nr > 1: @@ -612,47 +604,157 @@ def ApplySymmetryToSecondOrderEffCharge(self, dM_drdr, apply_asr = True): dM_drdr[:,:,:] = new_dM / self.QE_translation_nr new_dM[:,:,:] = 0 - - #print("After transl:") - #print(dM_drdr[:,:,0]) - - #self.PrintSymmetries() + # # DEBUG VARIABLE + # debug = np.zeros(np.shape(new_dM)) # Apply rotations for i in range(self.QE_nsym): irt = self.QE_irt[i, :] - 1 - - #print("") - #print("--------------------") - #print("symmetry: {:d}, irt: {}".format(i+1, irt +1)) - - #prova = np.zeros(np.shape(new_dM)) - for jat in range(nat): for kat in range(nat): new_mat = dM_drdr[3*irt[jat]: 3*irt[jat]+3, 3*irt[kat]:3*irt[kat] + 3,:] # Apply the symmetries - new_mat = np.einsum("ck, ijk->ijc", self.QE_s[:,:,i], new_mat) - new_mat = np.einsum("bj, ijc->ibc", self.QE_s[:,:,i], new_mat) - new_mat = np.einsum("ai, ibc->abc", self.QE_s[:,:,i], new_mat) - #prova[3*jat:3*jat+3, 3*kat:3*kat+3,:] = new_mat + new_mat = np.einsum("ck, ijk -> ijc", self.QE_s[:,:,i], new_mat) + new_mat = np.einsum("bj, ijc -> ibc", self.QE_s[:,:,i], new_mat) + new_mat = np.einsum("ai, ibc -> abc", self.QE_s[:,:,i], new_mat) new_dM[3*jat:3*jat+3, 3*kat:3*kat+3,:] += new_mat - #print(np.einsum("abc->cab", prova)) - #print("--------------------") dM_drdr[:,:,:] = new_dM / self.QE_nsym + + # # CONVERT IN CARTESIAN COORDINATES TO DEBUG + # for _i_ in range(nat): + # for _j_ in range(nat): + # debug[3*_i_ : 3*_i_ + 3, 3*_j_ : 3*_j_ + 3, :] = Methods.convert_3tensor_to_cryst(new_dM[3*_i_ : 3*_i_ + 3, 3*_j_ : 3*_j_ + 3, :],\ + # self.QE_at.T, cryst_to_cart = True) + # np.save('CC_new_{}'.format(i), debug) - - # Convert in crystal coordinates + # Convert in cartesian coordinates for i in range(nat): for j in range(nat): dM_drdr[3*i : 3*i + 3, 3*j: 3*j+3, :] = Methods.convert_3tensor_to_cryst(dM_drdr[3*i:3*i+3, 3*j:3*j+3,:], self.QE_at.T, True) + def ApplySymmetryToSecondOrderRamanTensor(self, dalpha_drdr, apply_asr = True): + """ + SYMMETRIZE TWO PHONON EFFECTIVE CHARGES + ======================================= + + This subroutine applies simmetries to the second order Raman tensor. + + Note, to symmetrize this tensor, symmetries must be imposed on the supercell. + + Parameters + ---------- + dalpha_drdr : ndarray (size = (3, 3, 3 nat_sc, 3 nat_sc)) + The second derivative of polarizability. + apply_asr : bool + If True the sum rule is applied. + The sum rule is the 'custom' one where translations are projected + out from the space for each polarization components. + """ + # Check the shape of the tensor + E1, E2, nat3, nat3_ = np.shape(dalpha_drdr) + + assert nat3 == nat3_, "Error on the shape of the argument for the atomic indices" + assert nat3 == 3 * self.QE_nat, "Wrong number of atoms (Symmetries must be setup in the supercell)" + assert E1 == E2, "Error on the shape of the argument electric field" + assert E1 == 3, "The first two entries are assosciated with electric field" + + # Get the number of atoms in the supercell + nat = int(nat3 /3) + + # Apply hermitianity on the atomic indices + dalpha_drdr += np.einsum("abcd->abdc", dalpha_drdr) + dalpha_drdr /= 2 + + # Apply hermitianity on the electric field + dalpha_drdr += np.einsum("abcd->bacd", dalpha_drdr) + dalpha_drdr /= 2 + + # SUM RULE ranging on the electric field components + if apply_asr: + for pol1 in range(3): + for pol2 in range(3): + CustomASR(dalpha_drdr[pol1, pol2, :, :]) + + # CONVERT TO CRYSTAL COORDINATES + for i in range(nat): + for j in range(nat): + dalpha_drdr[:, :, 3*i : 3*i + 3, 3*j: 3*j + 3] = Methods.convert_4tensor_to_cryst(dalpha_drdr[:, :, 3*i : 3*i + 3, 3*j: 3*j + 3], self.QE_at.T) + + + # Get a ZERO second order Raman tensor (AUXILIARY) + new_dalpha_drdr = np.zeros(np.shape(dalpha_drdr), dtype = np.double) + + # TRANSLATIONS + if self.QE_translation_nr > 1: + for i in range(self.QE_translation_nr): + # irt[at1] is the atom on which the translation i maps at1 + irt = self.QE_translations_irt[:, i] - 1 + for at1 in range(nat): + for at2 in range(nat): + # Get the part of the tensor that is equivalent by translations of atom at1 at2 + new_mat = dalpha_drdr[:, :, 3*irt[at1]: 3*irt[at1] + 3, 3 * irt[at2]: 3*irt[at2] + 3] + # Fill with the symmetric counterparts + new_dalpha_drdr[:, :, 3*at1: 3*at1+3, 3*at2:3*at2+3] += new_mat + + # OVERWRITE the second order Raman tensor + dalpha_drdr[:,:,:,:] = new_dalpha_drdr / self.QE_translation_nr + # SET TO ZERO THE AUXILIARY VARIABLE + new_dalpha_drdr[:,:,:,:] = 0 + + # # DEBUG VARIABLE + # debug = np.zeros(np.shape(new_dalpha_drdr)) + + # ROTATIONS + for i in range(self.QE_nsym): + # the symmetry applied on irt[at] gives the atom at + irt = self.QE_irt[i, :] - 1 + + for at1 in range(nat): + for at2 in range(nat): + # Get the part of the tensor that is equivalent by rotations of atom at1 at2 + # This has shape = (3, 3, 3, 3) + new_mat = dalpha_drdr[:, :, 3*irt[at1] : 3*irt[at1] + 3, 3*irt[at2] : 3*irt[at2] + 3] + + # Apply the symmetries + new_mat = np.einsum("dl, ijkl -> ijkd", self.QE_s[:,:,i], new_mat) + new_mat = np.einsum("ck, ijkd -> ijcd", self.QE_s[:,:,i], new_mat) + new_mat = np.einsum("bj, ijcd -> ibcd", self.QE_s[:,:,i], new_mat) + new_mat = np.einsum("ai, ibcd -> abcd", self.QE_s[:,:,i], new_mat) + + new_dalpha_drdr[:, :, 3*at1: 3*at1 + 3, 3*at2: 3*at2 + 3] += new_mat + + # # CONVERT IN CARTESIAN COORDINATES TO DEBUG + # for _i_ in range(nat): + # for _j_ in range(nat): + # debug[:, :, 3*_i_ : 3*_i_ + 3, 3*_j_ : 3*_j_ + 3] = Methods.convert_4tensor_to_cryst(new_dalpha_drdr[:, :, 3*_i_ : 3*_i_ + 3, 3*_j_ : 3*_j_ + 3],\ + # self.QE_at.T, cryst_to_cart = True) + # np.save('CC_new_{}'.format(i), debug) + + # OVERWRITE THE second order Raman tensor + dalpha_drdr[:,:,:,:] = new_dalpha_drdr /self.QE_nsym + + + # Convert BACK in crystal coordinates + for i in range(nat): + for j in range(nat): + dalpha_drdr[:, :, 3*i : 3*i + 3, 3*j: 3*j+3] = Methods.convert_4tensor_to_cryst(dalpha_drdr[:, :, 3*i : 3*i + 3, 3*j : 3*j + 3],\ + self.QE_at.T, cryst_to_cart = True) + + # # TO DEBUG + # np.save('CC_raman', dalpha_drdr) + + return + + + + + def ApplySymmetryToTensor4(self, v4, initialize_symmetries = True): """ SYMMETRIZE A RANK-4 TENSOR diff --git a/tests/TestASR/asr_scha.py b/tests/TestASR/asr_scha.py new file mode 100644 index 00000000..6298fc3e --- /dev/null +++ b/tests/TestASR/asr_scha.py @@ -0,0 +1,126 @@ +import cellconstructor as CC +import cellconstructor.Phonons +import cellconstructor.Structure +import cellconstructor.Units +import os +import numpy as np +import sys + +if __name__ == "__main__": + """ + THE MEANING OF THE ASR IN THE SCHA + """ + struct = CC.Structure.Structure(2) + struct.atoms = ['H', 'C'] + struct.coords = np.array([[-1, 0., 0.], + [+1, 0., 0.]]) + struct.masses = {'H' : 1837, 'C' : 1837 * 12} + struct.unit_cell = np.eye(3) * 100 + dyn = CC.Phonons.Phonons(struct) + dyn.dynmats[0] = np.random.uniform(size = (6,6)) + dyn.dynmats[0] += dyn.dynmats[0] + dyn.ForcePositiveDefinite() + dyn.ApplySumRule() + Temp = 100 # Kelvin + masses = np.tile(struct.get_masses_array(), (3, 1)).T.ravel() + + T = np.zeros((3, 6)) + for i in range(3): + v = np.zeros(3) + v[i] = 1 + T[i, :] = np.tile(v, 2) + + Y = dyn.GetUpsilonMatrix(Temp) + + # Multiply by a translation + FC_T = np.einsum('ab, ib -> ia', dyn.dynmats[0], T) + + Y_T = np.einsum('ab, ib -> ia', Y, T) + + print('\n### TEST ###') + print('Phi dot T') + for i in range(3): + print('trasl {}'.format(i)) + print(FC_T[i,:]) + if np.any(np.abs(FC_T[i, :]) > 1e-3): + raise ValueError('The traslation of FC does not work') + + print('\n### TEST ###') + print('Y dot T') + for i in range(3): + print('trasl {}'.format(i)) + print(Y_T[i,:]) + if np.any(np.abs(Y_T[i, :]) > 1e-3): + raise ValueError('The traslation of Y does not work') + + + # The matrix A of SCHA + def get_A(T, w, pols, masses): + A = np.zeros((6,6)) + if T < 1e-3: + return A + + a_mu = np.zeros(6) + mask = np.abs(w) < 1e-4 + a_mu[mask] = T * (CC.Units.RY_TO_KELVIN**-1) # beta^-1 + # The BE occupations + n = (np.exp(w[~mask] * CC.Units.RY_TO_KELVIN /T) - 1)**-1 + # The other eigenvalues + a_mu[~mask] = 2 * w[~mask] * n * (1 + n) /(1 + 2 * n) + + A = np.einsum('m, am, bm -> ab', a_mu, pols, pols) + # Multiply by the masses + return np.einsum('a, ab, b -> ab', np.sqrt(masses), A, np.sqrt(masses)) + + # The matrix Y of SCHA + def get_Y(T, w, pols, masses): + Y = np.zeros((6,6)) + + if T < 1e-3: + return A + + y_mu = np.zeros(6) + mask = np.abs(w) < 1e-4 + y_mu[mask] = 0. + n = (np.exp(w[~mask] * CC.Units.RY_TO_KELVIN /T) - 1)**-1 + y_mu[~mask] = 2 * w[~mask] /(1 + 2 * n) + + Y = np.einsum('m, am, bm -> ab', y_mu, pols, pols) + + return np.einsum('a, ab, b -> ab', np.sqrt(masses), Y, np.sqrt(masses)) + + + + w, pols = dyn.DyagDinQ(0) + A = get_A(Temp, w, pols, masses) + print('\nSCHA frequencies cm-1') + print(w[np.abs(w) > 1e-3] * CC.Units.RY_TO_CM) + + print('\n\nTest of the COM kinetic energy @ {} K'.format(Temp)) + # Check the sum in the COM kinetic energy + totA = 0. + _A_ = A.reshape((2, 3, 2, 3)) + for i in range(2): + for j in range(2): + for alpha in range(3): + totA += _A_[i, alpha, j, alpha] + + print('>The A contribution') + print( totA) + + + # Check the sum in the COM kinetic energy + totY = 0. + _Y_ = Y.reshape((2, 3, 2, 3)) + for i in range(2): + for j in range(2): + for alpha in range(3): + totY += _Y_[i, alpha, j, alpha] + + print(' >The Y contribution') + print( totY) + + print(np.einsum('am, bm -> ab', pols, pols)) + + # print(np.abs(Y - get_Y(Temp, w, pols, masses)).max()) + # print(np.abs(Y - get_Y(Temp, w, pols, masses)).min()) \ No newline at end of file