From 36d7ebaa168c1a616a712b6a28a3860fab624a05 Mon Sep 17 00:00:00 2001 From: DjordjeDangic Date: Mon, 24 Jul 2023 16:30:41 +0200 Subject: [PATCH 1/8] Small bug in group velocities --- cellconstructor/ThermalConductivity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cellconstructor/ThermalConductivity.py b/cellconstructor/ThermalConductivity.py index a20b97a4..a897d42b 100644 --- a/cellconstructor/ThermalConductivity.py +++ b/cellconstructor/ThermalConductivity.py @@ -2511,7 +2511,7 @@ def get_group_velocity(self, q, freqs, eigvecs): for iuc in range(len(self.force_constants)): for iat in range(len(uc_positions)): for jat in range(len(uc_positions)): - ruc = -self.ruc[iuc] #+ uc_positions[iat] - uc_positions[jat] + ruc = -self.ruc[iuc] + uc_positions[iat] - uc_positions[jat] phase = np.dot(ruc, q)*2.0*np.pi auxfc[3*iat:3*(iat+1),3*jat:3*(jat+1)] += complex(0.0,1.0)*ruc[icart]*self.force_constants[iuc,3*iat:3*(iat+1),3*jat:3*(jat+1)]*np.exp(1j*phase) ddynmat.append(auxfc * mm_inv_mat) From 78afc81c3118560eb5f55fceff31167a4cf849de Mon Sep 17 00:00:00 2001 From: DjordjeDangic Date: Fri, 4 Aug 2023 14:36:48 +0200 Subject: [PATCH 2/8] Cleaning spectral_kappa --- cellconstructor/ThermalConductivity.py | 75 +++++++++++++++++--------- 1 file changed, 50 insertions(+), 25 deletions(-) diff --git a/cellconstructor/ThermalConductivity.py b/cellconstructor/ThermalConductivity.py index a897d42b..a5d95d87 100644 --- a/cellconstructor/ThermalConductivity.py +++ b/cellconstructor/ThermalConductivity.py @@ -27,12 +27,12 @@ import time -__SEEKPATH=False +__SEEKPATH__ = False try: import seekpath - __SEEKPATH=True + __SEEKPATH__ = True except: - __SEEKPATH=False + __SEEKPATH__ = False try: from mpi4py import MPI @@ -298,7 +298,7 @@ def get_kpoints_in_path(path, nkpts, kprim): else: distance.append(distance[-1]) start_dist = distance[-1] - nkpts1 = np.int(np.floor(length/dl)) + nkpts1 = int(np.floor(length/dl)) for j in range(nkpts1): newqpt = qpt2 + (qpt1 - qpt2)/float(nkpts1)*float(j + 1) kpoints.append(newqpt) @@ -853,11 +853,12 @@ def calc_spectral_kappa_gk_diag(self, temperature): ls_key = format(temperature, '.1f') spec_kappa = np.zeros((3,3,self.lineshapes[ls_key].shape[-1])) + spec_kappa_off = np.zeros((3,3,self.lineshapes[ls_key].shape[-1])) energies = np.arange(spec_kappa.shape[-1], dtype=float)*self.delta_omega + self.delta_omega - exponents = np.exp(energies*SSCHA_TO_THZ*1.0e12*HPLANCK/KB/temperature) - integrands_plus = self.lineshapes[ls_key]**2*exponents/(exponents - 1.0)**2 - exponents = np.exp(-1.0*energies*SSCHA_TO_THZ*1.0e12*HPLANCK/KB/temperature) - integrands_minus = self.lineshapes[ls_key]**2*exponents/(exponents - 1.0)**2 + exponents_plus = np.exp(energies*SSCHA_TO_THZ*1.0e12*HPLANCK/KB/temperature) + integrands_plus = self.lineshapes[ls_key]**2*exponents_plus/(exponents_plus - 1.0)**2 + exponents_minus = np.exp(-1.0*energies*SSCHA_TO_THZ*1.0e12*HPLANCK/KB/temperature) + integrands_minus = self.lineshapes[ls_key]**2*exponents_minus/(exponents_minus - 1.0)**2 integrands = (integrands_plus + integrands_minus) for istar in self.qstar: @@ -865,14 +866,32 @@ def calc_spectral_kappa_gk_diag(self, temperature): for iband in range(self.nband): if(self.freqs[iqpt, iband] != 0.0): if(self.off_diag): - gvel = np.zeros_like(self.gvels[iqpt, iband, iband]) - gvel_sum = np.zeros((3,3), dtype=complex) - for r in self.rotations: - rot_q = np.dot(self.reciprocal_lattice.T, np.dot(r.T, np.linalg.inv(self.reciprocal_lattice.T))) - gvel = np.dot(rot_q, self.gvels[iqpt, iband, iband]) - gvel_sum += np.outer(gvel.conj(), gvel) - gvel_sum = gvel_sum.real*SSCHA_TO_MS**2/float(len(self.rotations)) - spec_kappa += np.einsum('ij,k->ijk',gvel_sum,integrands[iqpt, iband])*self.freqs[iqpt, iband]**2 + for jband in range(self.nband): + if(self.freqs[iqpt, jband] != 0.0): + if(self.group_velocity_mode != 'wigner'): + vel_fact = 1.0 + else: + vel_fact = 2.0*np.sqrt(self.freqs[iqpt, jband]*self.freqs[iqpt, iband])/(self.freqs[iqpt, jband] + self.freqs[iqpt, iband]) # as per Eq.34 in Caldarelli et al + if(iband == jband): + gvel = np.zeros_like(self.gvels[iqpt, iband, iband]) + gvel_sum = np.zeros((3,3), dtype=complex) + for r in self.rotations: + rot_q = np.dot(self.reciprocal_lattice.T, np.dot(r.T, np.linalg.inv(self.reciprocal_lattice.T))) + gvel = np.dot(rot_q, self.gvels[iqpt, iband, iband]) + gvel_sum += np.outer(gvel.conj(), gvel) + gvel_sum = gvel_sum.real*SSCHA_TO_MS**2/float(len(self.rotations)) + spec_kappa += np.einsum('ij,k->ijk',gvel_sum,integrands[iqpt, iband])*self.freqs[iqpt, iband]**2 + else: + integrands_plus1 = self.lineshapes[ls_key][iqpt, iband]*self.lineshapes[ls_key][iqpt, jband]*exponents_plus/(exponents_plus - 1.0)**2 + integrands_minus1 = self.lineshapes[ls_key][iqpt, iband]*self.lineshapes[ls_key][iqpt, jband]*exponents_minus/(exponents_minus - 1.0)**2 + gvel = np.zeros_like(self.gvels[iqpt, iband, jband]) + gvel_sum = np.zeros((3,3), dtype=complex) + for r in self.rotations: + rot_q = np.dot(self.reciprocal_lattice.T, np.dot(r.T, np.linalg.inv(self.reciprocal_lattice.T))) + gvel = np.dot(rot_q, self.gvels[iqpt, iband, jband]) + gvel_sum += np.outer(gvel.conj(), gvel) + gvel_sum = gvel_sum.real/vel_fact**2/float(len(self.rotations)) + spec_kappa_off += np.einsum('ij,k->ijk',gvel_sum,integrands_plus1 + integrands_minus1)*self.freqs[iqpt,iband]*self.freqs[iqpt,jband]*SSCHA_TO_MS**2*(SSCHA_TO_THZ)*1.0e12*2.0*np.pi/2.0 else: gvel = np.zeros_like(self.gvels[iqpt, iband]) gvel_sum = np.zeros((3,3), dtype=complex) @@ -882,11 +901,18 @@ def calc_spectral_kappa_gk_diag(self, temperature): gvel_sum += np.outer(gvel.conj(), gvel) gvel_sum = gvel_sum.real*SSCHA_TO_MS**2/float(len(self.rotations)) spec_kappa += np.einsum('ij,k->ijk',gvel_sum,integrands[iqpt, iband])*self.freqs[iqpt, iband]**2 + for ie in range(np.shape(spec_kappa)[-1]): + spec_kappa[:,:,ie] += spec_kappa[:,:,ie].T + spec_kappa[:,:,ie] /= 2.0 + spec_kappa_off[:,:,ie] += spec_kappa_off[:,:,ie].T + spec_kappa_off[:,:,ie] /= 2.0 + spec_kappa = spec_kappa*HBAR_JS**2/KB/temperature**2/self.volume/float(self.nkpt)*1.0e30*np.pi*(SSCHA_TO_THZ*2.0*np.pi)*1.0e12/2.0 - tot_kappa = np.sum(spec_kappa, axis = len(spec_kappa) - 1)*self.delta_omega + spec_kappa_off = spec_kappa_off*HBAR_JS**2/KB/temperature**2/self.volume/float(self.nkpt)*1.0e30*np.pi + tot_kappa = np.sum(spec_kappa + spec_kappa_off, axis = len(spec_kappa) - 1)*self.delta_omega print('Total kappa is: ', np.diag(tot_kappa)) - return energies*SSCHA_TO_THZ, spec_kappa/SSCHA_TO_THZ + return energies*SSCHA_TO_THZ, spec_kappa/SSCHA_TO_THZ, spec_kappa_off/SSCHA_TO_THZ ################################################################################################################################## @@ -1191,7 +1217,7 @@ def calculate_kappa_gk_offdiag_mode_mixing(self, temperature, write_lineshapes, if(self.freqs[iqpt, kband] != 0.0): for lband in range(self.nband): if(self.freqs[iqpt, lband] != 0.0): - if(self.group_velocity_mode == 'wigner'): + if(self.group_velocity_mode != 'wigner'): vel_fact = 1.0 else: vel_fact = 2.0*np.sqrt(self.freqs[iqpt, jband]*self.freqs[iqpt, iband])/(self.freqs[iqpt, jband] + self.freqs[iqpt, iband]) # as per Eq.34 in Caldarelli et al @@ -1213,7 +1239,7 @@ def calculate_kappa_gk_offdiag_mode_mixing(self, temperature, write_lineshapes, gvel1 = np.dot(rot_q, self.gvels[iqpt, iband, jband]) gvel2 = np.dot(rot_q, self.gvels[iqpt, kband, lband]) gvel_sum += np.outer(gvel1.conj(), gvel2) - gvel_sum = gvel_sum.real*vel_fact**2/float(len(self.rotations)) + gvel_sum = gvel_sum.real/vel_fact**2/float(len(self.rotations)) if(iband == jband and iband == kband and kband == lband): kappa_diag += integrals*np.sqrt(self.freqs[iqpt,iband]*self.freqs[iqpt, kband]*self.freqs[iqpt,jband]*self.freqs[iqpt, lband])*\ gvel_sum*SSCHA_TO_MS**2*(SSCHA_TO_THZ)*1.0e12*2.0*np.pi/4.0 @@ -1282,7 +1308,7 @@ def calculate_kappa_gk_offdiag(self, temperature, write_lineshapes, energies, ga for iband in range(self.nband): if(self.freqs[iqpt, iband] != 0.0): for jband in range(self.nband): - if(self.group_velocity_mode == 'wigner'): + if(self.group_velocity_mode != 'wigner'): vel_fact = 1.0 else: vel_fact = 2.0*np.sqrt(self.freqs[iqpt, jband]*self.freqs[iqpt, iband])/(self.freqs[iqpt, jband] + self.freqs[iqpt, iband]) # as per Eq.34 in Caldarelli et al @@ -1296,7 +1322,7 @@ def calculate_kappa_gk_offdiag(self, temperature, write_lineshapes, energies, ga rot_q = np.dot(self.reciprocal_lattice.T, np.dot(r.T, np.linalg.inv(self.reciprocal_lattice.T))) gvel = np.dot(rot_q, self.gvels[iqpt, iband, jband]) gvel_sum += np.outer(gvel.conj(), gvel) - gvel_sum = gvel_sum.real*vel_fact**2/float(len(self.rotations)) + gvel_sum = gvel_sum.real/vel_fact**2/float(len(self.rotations)) #kappa_nondiag += integrals*self.freqs[iqpt,iband]**2*(self.freqs[iqpt, iband]**2 + self.freqs[iqpt, jband]**2)/self.freqs[iqpt][jband]/self.freqs[iqpt][iband]*\ # gvel_sum*SSCHA_TO_MS**2*(SSCHA_TO_THZ)*1.0e12*2.0*np.pi/4.0 kappa_nondiag += integrals*self.freqs[iqpt,iband]*self.freqs[iqpt,jband]*gvel_sum*SSCHA_TO_MS**2*(SSCHA_TO_THZ)*1.0e12*2.0*np.pi/2.0 @@ -1312,7 +1338,6 @@ def calculate_kappa_gk_offdiag(self, temperature, write_lineshapes, energies, ga integrands_minus = self.lineshapes[ls_key][iqpt, iband]**2*exponents_minus/(exponents_minus - 1.0)**2 integrals = (np.sum(integrands_plus, axis = len(integrands_plus.shape) - 1) + np.sum(integrands_minus, axis = len(integrands_plus.shape) - 1))*self.delta_omega kappa_diag += gvel_sum*integrals*self.freqs[iqpt][iband]**2*SSCHA_TO_MS**2*(SSCHA_TO_THZ)*1.0e12*2.0*np.pi/2.0 - # Here freqs[iqpt,iband] is squared instead of power of 4 because imag self-energy in SSCHA is defined as 2*freqs[iqpt,iband]*Gamma[iqpt, iband] # Factor of 1/2 comes from the fact that we multiplied the lineshapes with 2.0 after we calculated them! kappa_diag += kappa_diag.T @@ -1354,7 +1379,7 @@ def calculate_kappa_gk_AC(self, temperature, write_lineshapes, energies, gauss_s for iband in range(self.nband): if(self.freqs[iqpt, iband] != 0.0): for jband in range(self.nband): - if(self.group_velocity_mode == 'wigner'): + if(self.group_velocity_mode != 'wigner'): vel_fact = 1.0 else: vel_fact = 2.0*np.sqrt(self.freqs[iqpt, jband]*self.freqs[iqpt, iband])/(self.freqs[iqpt, jband] + self.freqs[iqpt, iband]) # as per Eq.34 in Caldarelli et al @@ -1650,7 +1675,7 @@ def get_lineshapes_along_the_line(self, temperature, ne = 1000, filename = 'spec """ - if(not __SEEKPATH and kpoints is None): + if(not __SEEKPATH__ and kpoints is None): raise RuntimeError('To automatically generated a line in reciprocal space one need seekpath. First do "pip install seekpath"!') start_time = time.time() From 24c8546b531b02f98d5ba9dc8f88b2796689a982 Mon Sep 17 00:00:00 2001 From: DjordjeDangic Date: Thu, 24 Aug 2023 08:13:06 +0200 Subject: [PATCH 3/8] Adding symmetrizetion of group velocities --- cellconstructor/ThermalConductivity.py | 78 ++++++++++++++++++++++++-- 1 file changed, 72 insertions(+), 6 deletions(-) diff --git a/cellconstructor/ThermalConductivity.py b/cellconstructor/ThermalConductivity.py index a5d95d87..721b13fe 100644 --- a/cellconstructor/ThermalConductivity.py +++ b/cellconstructor/ThermalConductivity.py @@ -1438,16 +1438,68 @@ def calculate_kappa_gk_AC(self, temperature, write_lineshapes, energies, gauss_s for ie in range(np.shape(kappa_diag)[-1]): kappa_diag[:,:,ie] += kappa_diag[:,:,ie].T kappa_diag = kappa_diag/2.0*HBAR_JS**2/KB/temperature**2/self.volume/float(self.nkpt)*1.0e30*np.pi - + kappa_diag *= (np.exp(energies*HBAR_RY*2.0*np.pi/KB/temperature) - 1.0)/(energies*HBAR_RY*2.0*np.pi/KB/temperature) for ie in range(np.shape(kappa_nondiag)[-1]): kappa_nondiag[:,:,ie] += kappa_nondiag[:,:,ie].T kappa_nondiag = kappa_nondiag/2.0*HBAR_JS**2/KB/temperature**2/self.volume/float(self.nkpt)*1.0e30*np.pi + kappa_nondiag *= (np.exp(energies*HBAR_RY*2.0*np.pi/KB/temperature) - 1.0)/(energies*HBAR_RY*2.0*np.pi/KB/temperature) return kappa_diag, kappa_nondiag ################################################################################################################################# + def get_self_energy_at_q(self, iqpt, temperature, energies, mode_mixing = 'no', gauss_smearing = False, write_self_energy = False): + + if(self.delta_omega == 0.0 and energies is not None): + self.delta_omega = energies[1] - energies[0] + + if(not self.set_up_scattering_grids): + self.set_scattering_grids_simple() + + is_q_gamma = CC.Methods.is_gamma(self.fc2.unitcell_structure.unit_cell, self.k_points[iqpt]) + + if(mode_mixing == 'mode_mixing'): + self_energy = thermal_conductivity.get_lf.calculate_self_energy_full(self.freqs[iqpt], self.k_points[iqpt], self.eigvecs[iqpt], is_q_gamma, \ + self.scattering_grids[iqpt].T, self.scattering_weights[iqpt], self.fc2.tensor, self.fc3.tensor, self.fc2.r_vector2, self.fc3.r_vector2, \ + self.fc3.r_vector3, self.dyn.structure.coords.T, self.reciprocal_lattice, self.dyn.structure.get_masses_array(), self.sigmas[iqpt], \ + temperature, energies, True, gauss_smearing, False, len(self.scattering_grids[iqpt]), self.dyn.structure.N_atoms, len(self.fc2.tensor), len(self.fc3.tensor), len(energies)) + elif(mode_mixing == 'no'): + self_energy = thermal_conductivity.get_lf.calculate_self_energy_p(self.freqs[iqpt], self.k_points[iqpt], self.eigvecs[iqpt], is_q_gamma, \ + self.scattering_grids[iqpt].T, self.scattering_weights[iqpt], self.fc2.tensor, self.fc3.tensor, self.fc2.r_vector2, self.fc3.r_vector2, \ + self.fc3.r_vector3, self.dyn.structure.coords.T, self.reciprocal_lattice, self.dyn.structure.get_masses_array(), self.sigmas[iqpt], \ + temperature, energies, True, gauss_smearing, False, len(self.scattering_grids[iqpt]), self.dyn.structure.N_atoms, len(self.fc2.tensor), len(self.fc3.tensor), len(energies)) + else: + raise RuntimeError('The chosen option for mode_mixing(' + mode_mixing + ') does not exist!') + if(gauss_smearing): + from scipy.signal import hilbert + real_part = hilbert(self_energy.imag) + self_energy.real = -1.0*real_part.imag + + if(write_self_energy): + with open('Self_energy_' + str(iqpt + 1), 'w+') as outfile: + outfile.write('# ' + format('Energy (THz)', STR_FMT)) + for iband in range(self.nband): + if(mode_mixing == 'mode_mixing'): + for jband in range(self.nband): + outfile.write(' ' + format('Self energy ' + str(iband) + ' - ' + str(jband) + ' (THz)', STR_FMT)) + else: + outfile.write(' ' + format('Self energy ' + str(iband) +' (THz)', STR_FMT)) + outfile.write('\n') + for ie in range(len(energies)): + outfile.write(3*' ' + format(energies[ie]*SSCHA_TO_THZ, '.12e')) + for iband in range(self.nband): + if(mode_mixing == 'mode_mixing'): + for jband in range(self.nband): + outfile.write(3*' ' + format(self_energy[ie, iband, jband].real*SSCHA_TO_THZ, '.12e') + ' ' + format(self_energy[ie, iband, jband].imag*SSCHA_TO_THZ, '.12e')) + else: + outfile.write(3*' ' + format(self_energy[ie, iband].real*SSCHA_TO_THZ, '.12e') + ' ' + format(self_energy[ie, iband].imag*SSCHA_TO_THZ, '.12e')) + outfile.write('\n') + + return self_energy + + ################################################################################################################################# + def get_lineshapes(self, temperature, write_lineshapes, energies, method = 'fortran', mode_mixing = 'no', gauss_smearing = False): """ @@ -2406,7 +2458,7 @@ def set_lifetimes(self, temperatures, lifetimes, freqs_shifts): #################################################################################################################################### - def setup_harmonic_properties(self, smearing_value = 0.00005, symmetrize = False): + def setup_harmonic_properties(self, smearing_value = 0.00005, symmetrize = True): """ @@ -2428,8 +2480,8 @@ def setup_harmonic_properties(self, smearing_value = 0.00005, symmetrize = False if(symmetrize): self.gvels[ikpt] = self.symmetrize_group_velocity_by_index(self.gvels[ikpt], ikpt) - if(symmetrize): - self.symmetrize_group_velocities_over_star() + #if(symmetrize): + # self.symmetrize_group_velocities_over_star() #self.check_group_velocities() #self.check_frequencies() self.setup_smearings(smearing_value) @@ -2540,14 +2592,22 @@ def get_group_velocity(self, q, freqs, eigvecs): phase = np.dot(ruc, q)*2.0*np.pi auxfc[3*iat:3*(iat+1),3*jat:3*(jat+1)] += complex(0.0,1.0)*ruc[icart]*self.force_constants[iuc,3*iat:3*(iat+1),3*jat:3*(jat+1)]*np.exp(1j*phase) ddynmat.append(auxfc * mm_inv_mat) + ddynmat[-1] += ddynmat[-1].conj().T + ddynmat[-1] /= 2.0 if(icart == 0): dirdynmat = ddynmat.copy() + + new_eigvecs = np.zeros_like(eigvecs) + for deg in degs: + _, eigvecs1 = np.linalg.eigh(np.dot(eigvecs[:,deg].T.conj(), np.dot(np.sum(ddynmat,axis=0)/3.0, eigvecs[:,deg]))) + new_eigvecs[:,deg] = np.dot(eigvecs[:,deg], eigvecs1) + ddynmat = np.array(ddynmat) tmp_gvel = [] freqs_matrix = np.einsum('i,j->ij', freqs, freqs) freqs_matrix = np.divide(np.ones_like(freqs_matrix), freqs_matrix, out=np.zeros_like(freqs_matrix), where=freqs_matrix!=0.0) for icart in range(3): - tmp_gvel.append(np.dot(eigvecs.T.conj(), np.dot(ddynmat[icart], eigvecs))/2.0*np.sqrt(freqs_matrix)) + tmp_gvel.append(np.dot(new_eigvecs.T.conj(), np.dot(ddynmat[icart], new_eigvecs))/2.0*np.sqrt(freqs_matrix)) tmp_gvel = np.array(tmp_gvel).transpose((1,2,0)) if(not self.off_diag): @@ -2675,11 +2735,17 @@ def get_group_velocity_finite_difference(self, q, freqs, eigvecs): if(icart == 0): dirdynmat = ddynmat[0].copy() ddynmat = np.array(ddynmat) + + new_eigvecs = np.zeros_like(eigvecs) + for deg in degs: + _, eigvecs1 = np.linalg.eigh(np.dot(eigvecs[:,deg].T.conj(), np.dot(np.sum(ddynmat,axis=0)/3.0, eigvecs[:,deg]))) + new_eigvecs[:,deg] = np.dot(eigvecs[:,deg], eigvecs1) + tmp_gvel = [] freqs_matrix = np.einsum('i,j->ij', freqs, freqs) freqs_matrix = np.divide(np.ones_like(freqs_matrix), freqs_matrix, out=np.zeros_like(freqs_matrix), where=freqs_matrix!=0.0) for icart in range(3): - tmp_gvel.append(np.dot(eigvecs.T.conj(), np.dot(ddynmat[icart], eigvecs))/2.0*np.sqrt(freqs_matrix)) + tmp_gvel.append(np.dot(new_eigvecs.T.conj(), np.dot(ddynmat[icart], new_eigvecs))/2.0*np.sqrt(freqs_matrix)) tmp_gvel = np.array(tmp_gvel).transpose((1,2,0)) if(not self.off_diag): From 65e01acb98b7900c72b0f7523d7b27fbd88dd02e Mon Sep 17 00:00:00 2001 From: DjordjeDangic Date: Tue, 5 Sep 2023 08:48:31 +0200 Subject: [PATCH 4/8] Added optional phase factor --- FModules/get_lf.f90 | 2 +- cellconstructor/ThermalConductivity.py | 55 ++++++++++++++++++++------ 2 files changed, 45 insertions(+), 12 deletions(-) diff --git a/FModules/get_lf.f90 b/FModules/get_lf.f90 index 35364227..d51c179b 100644 --- a/FModules/get_lf.f90 +++ b/FModules/get_lf.f90 @@ -345,7 +345,7 @@ subroutine calculate_lineshapes_cartesian(irrqgrid, qgrid, weights, scatt_events print*, 'Negative eigenvalue of dynamical matrix!' w_neg_freqs = .True. endif - print*, 'Interpolate frequency' +! print*, 'Interpolate frequency' if(.not. w_neg_freqs) then w_q = sqrt(w2_q) self_energy(:,:,:) = complex(0.0_DP, 0.0_DP) diff --git a/cellconstructor/ThermalConductivity.py b/cellconstructor/ThermalConductivity.py index 721b13fe..d80c7b73 100644 --- a/cellconstructor/ThermalConductivity.py +++ b/cellconstructor/ThermalConductivity.py @@ -452,7 +452,7 @@ def load_thermal_conductivity(filename = 'tc.pkl'): class ThermalConductivity: - def __init__(self, dyn, tensor3, kpoint_grid = 2, scattering_grid = None, smearing_scale = 1.0, smearing_type = 'adaptive', cp_mode = 'quantum', group_velocity_mode = 'analytical', off_diag = False): + def __init__(self, dyn, tensor3, kpoint_grid = 2, scattering_grid = None, smearing_scale = 1.0, smearing_type = 'adaptive', cp_mode = 'quantum', group_velocity_mode = 'analytical', off_diag = False, phase_conv = 'smooth'): """ @@ -472,7 +472,8 @@ def __init__(self, dyn, tensor3, kpoint_grid = 2, scattering_grid = None, smeari smearing_type : Type of smearing used. Could be constant (same for all phonon modes) or adaptive (scaled by the phonon group velocity and the q point density). cp_mode : Flag determining how phonon occupation factors are calculated (quantum/classical), default is quantum group_velocity_mode : How to calculate group velocities. 'analytical', 'finite_difference', 'wigner' - off_diag : Boolean parameter for the calculation of the off-diagonal elements of group velocity. + off_diag : Boolean parameter for the calculation of the off-diagonal elements of group velocity. + phase_conv : Phase convention for Fourier interpolation. Smooth (wrt atomic positions) or step (wrt to lattice vectors) """ @@ -498,6 +499,7 @@ def __init__(self, dyn, tensor3, kpoint_grid = 2, scattering_grid = None, smeari self.supercell = self.dyn.structure.generate_supercell(dyn.GetSupercell()).unit_cell #print('Supercell: ') #print(self.supercell) + self.phase_conv = phase_conv self.smearing_type = smearing_type self.cp_mode = cp_mode self.off_diag = off_diag @@ -1094,7 +1096,7 @@ def calculate_kappa(self, temperatures = [300.0], write_lifetimes = True, mode = self.delta_omega = np.amax(self.freqs)*2.0/float(ne) energies = np.arange(ne, dtype=float)*self.delta_omega + self.delta_omega if(self.off_diag): - kappa, kappa_nondiag = self.calculate_kappa_gk_AC(temperatures[itemp], write_lineshapes, energies, gauss_smearing = gauss_smearing) + kappa, kappa_nondiag, im_kappa, im_kappa_nondiag = self.calculate_kappa_gk_AC(temperatures[itemp], write_lineshapes, energies, gauss_smearing = gauss_smearing) kappa_file.write(3*' ' + format(temperatures[itemp], '.12e')) for icart in range(3): kappa_file.write(3*' ' + format(kappa[icart][icart][0], '.12e')) @@ -1125,6 +1127,16 @@ def calculate_kappa(self, temperatures = [300.0], write_lifetimes = True, mode = outfile.write(3*' ' + format(kappa_nondiag[0][1][ien], '.12e')) outfile.write(3*' ' + format(kappa_nondiag[1][2][ien], '.12e')) outfile.write(3*' ' + format(kappa_nondiag[2][0][ien], '.12e')) + for icart in range(3): + outfile.write(3*' ' + format(im_kappa[icart][icart][ien], '.12e')) + outfile.write(3*' ' + format(im_kappa[0][1][ien], '.12e')) + outfile.write(3*' ' + format(im_kappa[1][2][ien], '.12e')) + outfile.write(3*' ' + format(im_kappa[2][0][ien], '.12e')) + for icart in range(3): + outfile.write(3*' ' + format(im_kappa_nondiag[icart][icart][ien], '.12e')) + outfile.write(3*' ' + format(im_kappa_nondiag[0][1][ien], '.12e')) + outfile.write(3*' ' + format(im_kappa_nondiag[1][2][ien], '.12e')) + outfile.write(3*' ' + format(im_kappa_nondiag[2][0][ien], '.12e')) outfile.write('\n') else: print('Can not recognize this method of calculating kappa! ') @@ -1370,6 +1382,9 @@ def calculate_kappa_gk_AC(self, temperature, write_lineshapes, energies, gauss_s ne = len(energies) exponents_plus = np.exp(energies*SSCHA_TO_THZ*1.0e12*HPLANCK/KB/temperature) exponents_minus = np.exp(-1.0*energies*SSCHA_TO_THZ*1.0e12*HPLANCK/KB/temperature) + x = energies*SSCHA_TO_THZ*1.0e12*HPLANCK/KB/temperature + scale = (np.exp(x) - 1.0)/x + scale = np.insert(scale, 0, 1.0) kappa_diag = np.zeros((3,3, ne + 1)) kappa_diag1 = np.zeros((3,3)) @@ -1438,14 +1453,26 @@ def calculate_kappa_gk_AC(self, temperature, write_lineshapes, energies, gauss_s for ie in range(np.shape(kappa_diag)[-1]): kappa_diag[:,:,ie] += kappa_diag[:,:,ie].T kappa_diag = kappa_diag/2.0*HBAR_JS**2/KB/temperature**2/self.volume/float(self.nkpt)*1.0e30*np.pi - kappa_diag *= (np.exp(energies*HBAR_RY*2.0*np.pi/KB/temperature) - 1.0)/(energies*HBAR_RY*2.0*np.pi/KB/temperature) for ie in range(np.shape(kappa_nondiag)[-1]): kappa_nondiag[:,:,ie] += kappa_nondiag[:,:,ie].T kappa_nondiag = kappa_nondiag/2.0*HBAR_JS**2/KB/temperature**2/self.volume/float(self.nkpt)*1.0e30*np.pi - kappa_nondiag *= (np.exp(energies*HBAR_RY*2.0*np.pi/KB/temperature) - 1.0)/(energies*HBAR_RY*2.0*np.pi/KB/temperature) - return kappa_diag, kappa_nondiag + for i in range(3): + for j in range(3): + kappa_diag[i,j] *= scale + kappa_nondiag[i,j] *= scale + + from scipy.signal import hilbert + im_kappa_diag = np.zeros_like(kappa_diag) + im_kappa_nondiag = np.zeros_like(kappa_nondiag) + for i in range(3): + for j in range(3): + im_kappa_diag[i,j] = hilbert(kappa_diag[i,j]).imag + im_kappa_nondiag[i,j] = hilbert(kappa_nondiag[i,j]).imag + + + return kappa_diag, kappa_nondiag, im_kappa_diag, im_kappa_nondiag ################################################################################################################################# @@ -2933,11 +2960,17 @@ def get_dynamical_matrix(self, q, q_direct = None, lo_to_splitting = True): mm_inv_mat = 1.0 / mm_mat dynmat = np.zeros_like(self.force_constants[0], dtype = complex) for ir in range(len(self.ruc)): - for iat in range(len(uc_positions)): - for jat in range(len(uc_positions)): - r = -1.0*self.ruc[ir] + uc_positions[iat] - uc_positions[jat] - phase = np.dot(r, q)*2.0*np.pi - dynmat[3*iat:3*(iat+1),3*jat:3*(jat+1)] += self.force_constants[ir,3*iat:3*(iat+1),3*jat:3*(jat+1)]*np.exp(1j*phase) + if(self.phase_conv == 'smooth'): + for iat in range(len(uc_positions)): + for jat in range(len(uc_positions)): + r = -1.0*self.ruc[ir] + uc_positions[iat] - uc_positions[jat] + phase = np.dot(r, q)*2.0*np.pi + dynmat[3*iat:3*(iat+1),3*jat:3*(jat+1)] += self.force_constants[ir,3*iat:3*(iat+1),3*jat:3*(jat+1)]*np.exp(1j*phase) + else: + r = -1.0*self.ruc[ir] + phase = np.dot(r, q)*2.0*np.pi + dynmat += self.force_constants[ir]*np.exp(1j*phase) + dynmat = dynmat#*mm_inv_mat #dynmat = (dynmat + dynmat.conj().T)/2.0 From 905f9053071f33373e5d656961377ac54ea91e00 Mon Sep 17 00:00:00 2001 From: DjordjeDangic Date: Mon, 25 Sep 2023 11:56:34 +0200 Subject: [PATCH 5/8] Added functions to symmetrize dynamical matrices --- FModules/get_lf.f90 | 2 +- cellconstructor/ThermalConductivity.py | 244 +++++++++++++++++++++++-- 2 files changed, 233 insertions(+), 13 deletions(-) diff --git a/FModules/get_lf.f90 b/FModules/get_lf.f90 index d51c179b..4d8facdb 100644 --- a/FModules/get_lf.f90 +++ b/FModules/get_lf.f90 @@ -361,7 +361,7 @@ subroutine calculate_lineshapes_cartesian(irrqgrid, qgrid, weights, scatt_events curr_grid(:,jqpt) = qgrid(:,prev_events + jqpt) curr_w(jqpt) = weights(prev_events + jqpt) enddo - print*, 'Got grids' +! print*, 'Got grids' call calculate_self_energy_full(w_q, qpt, pols_q, is_q_gamma, scatt_events(iqpt), nat, nfc2, & nfc3, ne, curr_grid, curr_w, fc2, fc3, r2_2, r3_2, r3_3, pos, kprim, masses, smear(:,iqpt), T, & energies, .not. parallelize, gaussian, classical, self_energy) diff --git a/cellconstructor/ThermalConductivity.py b/cellconstructor/ThermalConductivity.py index d80c7b73..2a4ceb85 100644 --- a/cellconstructor/ThermalConductivity.py +++ b/cellconstructor/ThermalConductivity.py @@ -131,7 +131,7 @@ def get_spglib_cell(dyn): for iat in range(len(dyn.structure.atoms)): for jat in range(len(symbols_unique)): if(dyn.structure.atoms[iat] == symbols_unique[jat]): - numbers[iat] == jat + numbers[iat] = jat return (dyn.structure.unit_cell, pos, numbers) ######################################################################################################################################## @@ -449,6 +449,71 @@ def load_thermal_conductivity(filename = 'tc.pkl'): infile.close() return tc +def get_mapping_of_q_points(star, points, rotations): + + mapping = [] + for istar in range(len(star)): + star_mapping = [] + iqpt1 = star[istar][0] + qpt1 = points[iqpt1] + for iqpt in range(len(star[istar])): + found = False + qpt_mapping = [] + iqpt2 = star[istar][iqpt] + qpt2 = points[iqpt2] + if(np.linalg.norm(qpt2 + qpt1 - np.rint(qpt2 + qpt1)) < 1.0e-6): + qpt_mapping.append((-1, True)) + found = True + else: + for irot in range(len(rotations)): + qpt21 = np.dot(rotations[irot].T, qpt1) + diffq = qpt21 - qpt2 + addq = qpt21 + qpt2 + if(np.linalg.norm(diffq - np.rint(diffq)) < 1.0e-6): + qpt_mapping.append((irot, False)) + found = True + elif(np.linalg.norm(addq - np.rint(addq)) < 1.0e-6): + qpt_mapping.append((irot, True)) + found = True + if(found): + star_mapping.append(qpt_mapping) + else: + print('QPT1: ', qpt1) + print('QPT2: ', qpt2) + raise RuntimeError('Could not find any mapping between qpoints!') + mapping.append(star_mapping) + + return mapping + +def construct_symmetry_matrix(rotation, translation, qvec, pos, atom_map, cell): + matrix = np.zeros((len(pos)*3, len(pos)*3), dtype=complex) + rx1 = [] + for iat in range(len(pos)): + for jat in range(len(pos)): + if(iat == atom_map[jat]): + r = np.zeros_like(pos[iat]) + phase = 0.0 + r = pos[iat] - np.dot(rotation, pos[jat]) - translation + rx = np.dot(r, np.linalg.inv(cell)) + rx1.append(rx) + if(np.linalg.norm(rx - np.rint(rx))>1.0e-6): + raise RuntimeError('The atom is translated different than the translation vector!') + #else: + # print(rx) + #phase = np.dot(qvec, r)*2.0*np.pi + phase = -1.0*np.dot(qvec, r)*2.0*np.pi + #matrix[3*jat:3*(jat+1), 3*iat:3*(iat+1)] = rotation*np.exp(complex(0.0, phase)) + matrix[3*iat:3*(iat+1), 3*jat:3*(jat+1)] = rotation*np.exp(complex(0.0, phase)) + try: + imatrix = np.linalg.inv(matrix) + #print('Gamma is invertible!') + if(not np.linalg.norm(imatrix - matrix.conj().T)/np.linalg.norm(matrix) < 1.0e-5): + raise RuntimeError('The transformation matrix is not unitary!') + except: + print(np.matmul(matrix, matrix.conj().T)) + raise RuntimeError('Gamma is not invertible!') + return matrix#, rx1, np.exp(complex(0.0, phase)), phase, qvec, np.dot(qvec, r) + class ThermalConductivity: @@ -524,6 +589,7 @@ def __init__(self, dyn, tensor3, kpoint_grid = 2, scattering_grid = None, smeari self.freqs = np.zeros((self.nkpt, self.nband)) self.gruneisen = np.zeros((self.nkpt, self.nband)) self.eigvecs = np.zeros((self.nkpt, self.nband, self.nband), dtype=complex) + self.dynmats = np.zeros((self.nkpt, self.nband, self.nband), dtype=complex) if(self.off_diag): self.gvels = np.zeros((self.nkpt, self.nband, self.nband, 3), dtype = complex) else: @@ -554,13 +620,34 @@ def set_kpoints_spglib(self): cell = get_spglib_cell(self.dyn) mapping, grid = spglib.get_ir_reciprocal_mesh(self.kpoint_grid, cell, is_shift=[0, 0, 0]) - rotations = spglib.get_symmetry_dataset(cell)['rotations'] + symmetry_dataset = spglib.get_symmetry_dataset(cell, symprec=1e-5, angle_tolerance=-1.0, hall_number=0) + rotations = symmetry_dataset['rotations'] + translations = symmetry_dataset['translations'] + print(cell) + print('Spacegroup: ' + symmetry_dataset['international'] + str(symmetry_dataset['number'])) + self.atom_map = np.zeros((len(rotations), len(cell[1])), dtype=int) + nrot = len(rotations) + for irot in range(nrot): + found = [False for jat in range(len(cell[1]))] + new_pos = np.einsum('ij,jk->ik', cell[1], rotations[irot].T) + translations[irot] + for iat in range(len(cell[1])): + for jat in range(len(cell[1])): + diff = cell[1][jat] - new_pos[iat] + if(np.linalg.norm(diff - np.rint(diff)) < 1.0e-5): + self.atom_map[irot, iat] = jat + if(not found[iat]): + found[iat] = True + else: + print('Again found mapping to this atom!') + if(not all(found)): + raise RuntimeError('Could not find atom mapping for spacegroup symmetry: ' + str(irot + 1)) self.irr_k_points = np.array(grid[np.unique(mapping)] / np.array(self.kpoint_grid, dtype=float)) self.irr_k_points = np.dot(self.irr_k_points, self.reciprocal_lattice) self.qpoints = grid / np.array(self.kpoint_grid, dtype=float) self.k_points = np.dot(self.qpoints, self.reciprocal_lattice) self.rotations = np.array(rotations).copy() + self.translations = np.array(translations).copy() self.qstar = [] for i in np.unique(mapping): @@ -581,8 +668,8 @@ def set_kpoints_spglib(self): if(np.linalg.norm(diffq) < 1.0e-6): curr_little_group.append(irot) self.little_group[iqpt].extend(curr_little_group) - if(len(istar) * len(self.little_group[iqpt]) != len(rotations)): - raise RuntimeError('Number of symmetry operation wrong!', len(istar), len(self.little_group[iqpt]), len(rotations)) + #if(len(istar) * len(self.little_group[iqpt]) != len(rotations)): + # raise RuntimeError('Number of symmetry operation wrong!', len(istar), len(self.little_group[iqpt]), len(rotations)) mapping1, grid1 = spglib.get_ir_reciprocal_mesh(self.scattering_grid, cell, is_shift=[0, 0, 0]) self.scattering_qpoints = np.array(grid1 / np.array(self.scattering_grid, dtype=float)) @@ -1629,7 +1716,7 @@ def get_lineshapes(self, temperature, write_lineshapes, energies, method = 'fort else: print('Selected mode_mixing approach: ', mode_mixing) raise RuntimeError('Do not recognize the selected mode_mixing approach!') - + scaled_positions = np.dot(self.dyn.structure.coords, np.linalg.inv(self.unitcell)) for ikpt in range(self.nirrkpt): jkpt = self.qstar[ikpt][0] if(mode_mixing == 'no'): @@ -1640,11 +1727,49 @@ def get_lineshapes(self, temperature, write_lineshapes, energies, method = 'fort for iqpt in range(len(self.qstar[ikpt])): jqpt = self.qstar[ikpt][iqpt] + found = False if(mode_mixing != 'no'): #for iband in range(self.nband): # for jband in range(self.nband): # curr_ls[ikpt, iband,jband,:] = curr_ls[ikpt, iband,jband,:]/np.sum(curr_ls[ikpt,iband,jband,:])/(energies[1]-energies[0]) # Forcing the normalization. Not sure if the best option! - lineshapes[jqpt,:,:,:] = curr_ls[ikpt,:,:,:] + if(iqpt == 0): + lineshapes[jqpt,:,:,:] = curr_ls[ikpt,:,:,:] + found = True + else: + qpt1 = self.qpoints[self.qstar[ikpt][0]] + qpt2 = self.qpoints[jqpt] + if(np.linalg.norm(qpt2 + qpt1 - np.rint(qpt2 + qpt1)) < 1.0e-6): + lineshapes[jqpt,:,:,:] = 2.0*curr_ls[ikpt,:,:,:].conj() + found = True + else: + found = False + for irot in range(len(self.rotations)): + qpt21 = np.dot(self.rotations[irot].T, qpt1) + kpt21 = np.dot(qpt21, self.reciprocal_lattice) + diffq = qpt21 - qpt2 + addq = qpt21 + qpt2 + if(np.linalg.norm(diffq - np.rint(diffq)) < 1.0e-6): + rotation = self.rotations[irot] + translation = self.translations[irot] + atom_map = self.atom_map[irot] + gamma = construct_symmetry_matrix_for_dyn(rotation, translation, kpt21, scaled_positions, atom_map, self.unitcell, self.reciprocal_lattice) + #rotation = np.matmul(self.reciprocal_lattice.T, np.matmul(self.rotations[irot].T, np.linalg.inv(self.reciprocal_lattice.T))) + #translation = np.dot(self.translations[irot], self.unitcell) + lineshapes[jqpt,:,:,:] = 2.0*np.einsum('ij,jkl,km->iml', gamma, curr_ls[ikpt,:,:,:], gamma.conj().T) #np.matmul(gamma, np.matmul(dyn1, gamma.conj().T)) + found = True + break + elif(np.linalg.norm(addq - np.rint(addq)) < 1.0e-6): + rotation = self.rotations[irot] + translation = self.translations[irot] + atom_map = self.atom_map[irot] + gamma = construct_symmetry_matrix_for_dyn(rotation, translation, kpt21, scaled_positions, atom_map, self.unitcell, self.reciprocal_lattice) + #rotation = np.matmul(self.reciprocal_lattice.T, np.matmul(self.rotations[irot].T, np.linalg.inv(self.reciprocal_lattice.T))) + #translation = np.dot(self.translations[irot], self.unitcell) + lineshapes[jqpt,:,:,:] = 2.0*np.einsum('ij,jkl,km->iml', gamma, curr_ls[ikpt,:,:,:], gamma.conj().T).conj() #np.matmul(gamma, np.matmul(dyn1, gamma.conj().T)) + found = True + break + if(not found): + raise RuntimeError('Could not find mapping between q points in the star!') #tot_const_diag = 0.0 #tot_const_nondiag = 0.0 #for iband in range(len(lineshapes[jqpt])): @@ -1873,7 +1998,7 @@ def get_lineshapes_along_the_line(self, temperature, ne = 1000, filename = 'spec if(freqs[ikpt, iband] < np.amax(freqs[ikpt])*1.0e-6): curr_ls[ikpt, iband] = 0.0 if(mode_mixing != 'no'): - lineshapes[ikpt,:,:,:] = curr_ls[ikpt,:,:,:] + lineshapes[ikpt,:,:,:] = curr_ls[ikpt,:,:,:]*2.0 tot_const_diag = 0.0 tot_const_nondiag = 0.0 for iband in range(len(lineshapes[ikpt])): @@ -1952,6 +2077,7 @@ def get_lineshapes_along_the_line(self, temperature, ne = 1000, filename = 'spec outfile.write('\n') print('Calculated SSCHA lineshapes in: ', time.time() - start_time) + return energies, lineshapes ################################################################################################################################## def write_lineshape(self, filename, curr_ls, jkpt, energies, mode_mixing): @@ -2495,7 +2621,7 @@ def setup_harmonic_properties(self, smearing_value = 0.00005, symmetrize = True) """ for ikpt, kpt in enumerate(self.k_points): - self.freqs[ikpt], self.eigvecs[ikpt] = self.get_frequency_at_q(kpt) + self.freqs[ikpt], self.eigvecs[ikpt], self.dynmats[ikpt] = self.get_frequency_at_q(kpt) if(self.group_velocity_mode == 'wigner'): self.gvels[ikpt], self.ddms[ikpt] = self.get_group_velocity_wigner(kpt) elif(self.group_velocity_mode == 'analytical'): @@ -2511,6 +2637,7 @@ def setup_harmonic_properties(self, smearing_value = 0.00005, symmetrize = True) # self.symmetrize_group_velocities_over_star() #self.check_group_velocities() #self.check_frequencies() + #self.check_dynamical_matrices() self.setup_smearings(smearing_value) print('Harmonic properties are set up!') @@ -2526,6 +2653,97 @@ def symmetrize_eigenvectors(self): self.eigvecs[pairs[0]] = (self.eigvecs[pairs[0]] + self.eigvecs[pairs[1]].conj())/2.0 self.eigvecs[pairs[1]] = self.eigvecs[pairs[0]].conj() + ################################################################################################################################# + + def get_sg_in_cartesian(self): + + rotations = np.zeros_like(self.rotations, dtype=float) + translations = np.zeros_like(self.translations, dtype=float) + for irot in range(len(rotations)): + rotations[irot] = np.matmul(self.reciprocal_lattice.T, np.matmul(self.rotations[irot].T, np.linalg.inv(self.reciprocal_lattice.T))) + translations = np.dot(self.translations, self.unitcell) + + return rotations, translations + + ################################################################################################################################# + + def check_dynamical_matrices(self): + + rotations, translations = self.get_sg_in_cartesian() + #for irot in range(len(rotations)): + # print('SG' + str(irot + 1)) + # print(rotations[irot]) + # print(translations[irot]) + # print(tc.atom_map[irot]) + mapping = get_mapping_of_q_points(self.qstar, self.qpoints, self.rotations) + for istar in range(len(self.qstar)): + iqpt1 = self.qstar[istar][0] + qpt1 = self.qpoints[iqpt1] + kpt1 = self.k_points[iqpt1] + dyn1 = self.dynmats[iqpt1].copy() + for iqpt in range(1, len(self.qstar[istar])): + iqpt2 = self.qstar[istar][iqpt] + qpt2 = self.qpoints[iqpt2] + kpt2 = self.k_points[iqpt2] + dyn2 = self.dynmats[iqpt2].copy() + for imap in range(len(mapping[istar][iqpt])): + irot = mapping[istar][iqpt][imap][0] + atom_map = self.atom_map[irot] + conjugate = mapping[istar][iqpt][imap][1] + qpt21 = np.dot(self.rotations[irot].T, qpt1) + kpt22 = np.dot(rotations[irot], kpt1) + kpt21 = np.dot(self.reciprocal_lattice.T, qpt21) + if(np.linalg.norm(kpt21 - kpt22) > 1.0e-6): + print(np.linalg.norm(kpt21 - kpt22)) + print('Rotation in cartesian and reduced coordinates gives different results!') + diffq = qpt21 - qpt2 + addq = qpt21 + qpt2 + if(irot == -1): + dyn21 = dyn1.conj() + if(np.any(np.abs(dyn21 - dyn2)/np.amax(np.abs(dyn2)) > 1.0e-2)): + print('Some differences between rotated and original dynamical matrices!') + print(np.abs(dyn21 - dyn2)/np.amax(np.abs(dyn2)) > 1.0e-2) + print(qpt1, qpt2) + for iband in range(len(dyn1)): + print(dyn2[iband]) + print(dyn21[iband]) + print(dyn2[iband] - dyn21[iband]) + print('') + #raise RuntimeError('Mapping dynamical matrices from q to -q did not work!') + else: + if(np.linalg.norm(diffq - np.rint(diffq)) < 1.0e-6 and not conjugate): + gamma = construct_symmetry_matrix(rotations[irot], translations[irot], kpt21, self.dyn.structure.coords, atom_map, self.unitcell) + dyn21 = np.matmul(gamma, np.matmul(dyn1, gamma.conj().T)) + if(np.any(np.abs(dyn21 - dyn2)/np.amax(np.abs(dyn2)) > 1.0e-2)): + print('Some differences between rotated and original dynamical matrices!') + print(rotations[irot]) + print(np.abs(dyn21 - dyn2)/np.amax(np.abs(dyn2)) > 1.0e-2) + print(qpt1, qpt2) + for iband in range(len(dyn1)): + print(dyn2[iband]) + print(dyn21[iband]) + print(dyn2[iband] - dyn21[iband]) + print('') + #raise RuntimeError('Mapping dynamical matrices from q to q star did not work!') + elif(np.linalg.norm(addq - np.rint(addq)) < 1.0e-6 and conjugate): + gamma = construct_symmetry_matrix(rotations[irot], translations[irot], kpt21, self.dyn.structure.coords, atom_map, self.unitcell) + dyn21 = np.matmul(gamma, np.matmul(dyn1, gamma.conj().T)) + dyn21 = dyn21.conj() + if(np.any(np.abs(dyn21 - dyn2)/np.amax(np.abs(dyn2)) > 1.0e-2)): + print('Some differences between rotated and original dynamical matrices!') + print(rotations[irot]) + print(np.abs(dyn21 - dyn2)/np.amax(np.abs(dyn2)) > 1.0e-2) + print(qpt1, qpt2) + for iband in range(len(dyn1)): + print(dyn2[iband]) + print(dyn21[iband]) + print(dyn2[iband] - dyn21[iband]) + print('') + #raise RuntimeError('Mapping dynamical matrices from q to q star through -q did not work!') + else: + raise RuntimeError('The mapping was wrong! This rotation does not give expected q point!') + + print('Dynamical matrices satisfy symmetries!') ################################################################################################################################# @@ -2938,7 +3156,7 @@ def get_frequency_at_q(self, q): #else: w_q = np.sqrt(w2_q) - return w_q, pols_q + return w_q, pols_q, dynmat ################################################################################################################################### @@ -2948,7 +3166,7 @@ def get_dynamical_matrix(self, q, q_direct = None, lo_to_splitting = True): Get dynamical matrix at wave vector. - q : wave vector without 2pi factor + q : wave vector in cartesian coordinates without 2pi factor """ @@ -2966,12 +3184,14 @@ def get_dynamical_matrix(self, q, q_direct = None, lo_to_splitting = True): r = -1.0*self.ruc[ir] + uc_positions[iat] - uc_positions[jat] phase = np.dot(r, q)*2.0*np.pi dynmat[3*iat:3*(iat+1),3*jat:3*(jat+1)] += self.force_constants[ir,3*iat:3*(iat+1),3*jat:3*(jat+1)]*np.exp(1j*phase) - else: + elif(self.phase_conv == 'step'): r = -1.0*self.ruc[ir] phase = np.dot(r, q)*2.0*np.pi dynmat += self.force_constants[ir]*np.exp(1j*phase) + else: + raise RuntimeError('Can not recognize phase convention!') - dynmat = dynmat#*mm_inv_mat + #dynmat = dynmat#*mm_inv_mat #dynmat = (dynmat + dynmat.conj().T)/2.0 if self.fc2.effective_charges is not None: From 9ed0e7160d166cd8c9358b3dc9ed53d18467090c Mon Sep 17 00:00:00 2001 From: DjordjeDangic Date: Mon, 25 Sep 2023 15:31:24 +0200 Subject: [PATCH 6/8] Transform of lineshapes included! --- cellconstructor/ThermalConductivity.py | 41 +++++++------------------- 1 file changed, 11 insertions(+), 30 deletions(-) diff --git a/cellconstructor/ThermalConductivity.py b/cellconstructor/ThermalConductivity.py index 2a4ceb85..5f2cd245 100644 --- a/cellconstructor/ThermalConductivity.py +++ b/cellconstructor/ThermalConductivity.py @@ -623,7 +623,6 @@ def set_kpoints_spglib(self): symmetry_dataset = spglib.get_symmetry_dataset(cell, symprec=1e-5, angle_tolerance=-1.0, hall_number=0) rotations = symmetry_dataset['rotations'] translations = symmetry_dataset['translations'] - print(cell) print('Spacegroup: ' + symmetry_dataset['international'] + str(symmetry_dataset['number'])) self.atom_map = np.zeros((len(rotations), len(cell[1])), dtype=int) nrot = len(rotations) @@ -1717,6 +1716,8 @@ def get_lineshapes(self, temperature, write_lineshapes, energies, method = 'fort print('Selected mode_mixing approach: ', mode_mixing) raise RuntimeError('Do not recognize the selected mode_mixing approach!') scaled_positions = np.dot(self.dyn.structure.coords, np.linalg.inv(self.unitcell)) + rotations, translations = self.get_sg_in_cartesian() + mapping = get_mapping_of_q_points(self.qstar, self.qpoints, self.rotations) for ikpt in range(self.nirrkpt): jkpt = self.qstar[ikpt][0] if(mode_mixing == 'no'): @@ -1742,34 +1743,14 @@ def get_lineshapes(self, temperature, write_lineshapes, energies, method = 'fort lineshapes[jqpt,:,:,:] = 2.0*curr_ls[ikpt,:,:,:].conj() found = True else: - found = False - for irot in range(len(self.rotations)): - qpt21 = np.dot(self.rotations[irot].T, qpt1) - kpt21 = np.dot(qpt21, self.reciprocal_lattice) - diffq = qpt21 - qpt2 - addq = qpt21 + qpt2 - if(np.linalg.norm(diffq - np.rint(diffq)) < 1.0e-6): - rotation = self.rotations[irot] - translation = self.translations[irot] - atom_map = self.atom_map[irot] - gamma = construct_symmetry_matrix_for_dyn(rotation, translation, kpt21, scaled_positions, atom_map, self.unitcell, self.reciprocal_lattice) - #rotation = np.matmul(self.reciprocal_lattice.T, np.matmul(self.rotations[irot].T, np.linalg.inv(self.reciprocal_lattice.T))) - #translation = np.dot(self.translations[irot], self.unitcell) - lineshapes[jqpt,:,:,:] = 2.0*np.einsum('ij,jkl,km->iml', gamma, curr_ls[ikpt,:,:,:], gamma.conj().T) #np.matmul(gamma, np.matmul(dyn1, gamma.conj().T)) - found = True - break - elif(np.linalg.norm(addq - np.rint(addq)) < 1.0e-6): - rotation = self.rotations[irot] - translation = self.translations[irot] - atom_map = self.atom_map[irot] - gamma = construct_symmetry_matrix_for_dyn(rotation, translation, kpt21, scaled_positions, atom_map, self.unitcell, self.reciprocal_lattice) - #rotation = np.matmul(self.reciprocal_lattice.T, np.matmul(self.rotations[irot].T, np.linalg.inv(self.reciprocal_lattice.T))) - #translation = np.dot(self.translations[irot], self.unitcell) - lineshapes[jqpt,:,:,:] = 2.0*np.einsum('ij,jkl,km->iml', gamma, curr_ls[ikpt,:,:,:], gamma.conj().T).conj() #np.matmul(gamma, np.matmul(dyn1, gamma.conj().T)) - found = True - break - if(not found): - raise RuntimeError('Could not find mapping between q points in the star!') + irot = mapping[ikpt][iqpt][0][0] + atom_map = self.atom_map[irot] + qpt21 = np.dot(self.rotations[irot].T, qpt1) + kpt21 = np.dot(qpt21, self.reciprocal_lattice) + gamma = construct_symmetry_matrix(rotations[irot], translations[irot], kpt21, self.dyn.structure.coords, atom_map, self.unitcell) + lineshapes[jqpt,:,:,:] = 2.0*np.einsum('ij,jkl,km->iml', gamma, curr_ls[ikpt,:,:,:], gamma.conj().T) + if(mapping[ikpt][iqpt][0][1]): + lineshapes[jqpt,:,:,:] = lineshapes[jqpt,:,:,:].conj() #tot_const_diag = 0.0 #tot_const_nondiag = 0.0 #for iband in range(len(lineshapes[jqpt])): @@ -1918,7 +1899,7 @@ def get_lineshapes_along_the_line(self, temperature, ne = 1000, filename = 'spec nkpts = len(kpoints) freqs = np.zeros((nkpts, self.nband)) for ikpt in range(nkpts): - freqs[ikpt], _ = self.get_frequency_at_q(kpoints[ikpt]) + freqs[ikpt], _, _ = self.get_frequency_at_q(kpoints[ikpt]) maxfreq = np.amax(freqs)*2.1 energies = (np.arange(ne, dtype=float) + 1.0)/float(ne)*maxfreq From 088f6bc3c5ef4ad44f9a99c6dd6984346a9aed0e Mon Sep 17 00:00:00 2001 From: DjordjeDangic Date: Tue, 26 Sep 2023 11:32:05 +0200 Subject: [PATCH 7/8] Wrong scaling factor in cartesian lineshapes! --- cellconstructor/ThermalConductivity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cellconstructor/ThermalConductivity.py b/cellconstructor/ThermalConductivity.py index 5f2cd245..98090917 100644 --- a/cellconstructor/ThermalConductivity.py +++ b/cellconstructor/ThermalConductivity.py @@ -1740,7 +1740,7 @@ def get_lineshapes(self, temperature, write_lineshapes, energies, method = 'fort qpt1 = self.qpoints[self.qstar[ikpt][0]] qpt2 = self.qpoints[jqpt] if(np.linalg.norm(qpt2 + qpt1 - np.rint(qpt2 + qpt1)) < 1.0e-6): - lineshapes[jqpt,:,:,:] = 2.0*curr_ls[ikpt,:,:,:].conj() + lineshapes[jqpt,:,:,:] = curr_ls[ikpt,:,:,:].conj() found = True else: irot = mapping[ikpt][iqpt][0][0] @@ -1748,7 +1748,7 @@ def get_lineshapes(self, temperature, write_lineshapes, energies, method = 'fort qpt21 = np.dot(self.rotations[irot].T, qpt1) kpt21 = np.dot(qpt21, self.reciprocal_lattice) gamma = construct_symmetry_matrix(rotations[irot], translations[irot], kpt21, self.dyn.structure.coords, atom_map, self.unitcell) - lineshapes[jqpt,:,:,:] = 2.0*np.einsum('ij,jkl,km->iml', gamma, curr_ls[ikpt,:,:,:], gamma.conj().T) + lineshapes[jqpt,:,:,:] = np.einsum('ij,jkl,km->iml', gamma, curr_ls[ikpt,:,:,:], gamma.conj().T) if(mapping[ikpt][iqpt][0][1]): lineshapes[jqpt,:,:,:] = lineshapes[jqpt,:,:,:].conj() #tot_const_diag = 0.0 From 3d55795513028d7f21aa74c10c80b2b190a7471d Mon Sep 17 00:00:00 2001 From: DjordjeDangic Date: Tue, 26 Sep 2023 15:46:28 +0200 Subject: [PATCH 8/8] Some changes --- FModules/get_lf.f90 | 3 ++- cellconstructor/ThermalConductivity.py | 8 +++++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/FModules/get_lf.f90 b/FModules/get_lf.f90 index 4d8facdb..9a36fcf7 100644 --- a/FModules/get_lf.f90 +++ b/FModules/get_lf.f90 @@ -1315,8 +1315,9 @@ subroutine calculate_self_energy_full(w_q, qpt, pols_q, is_q_gamma, nqpt, nat, n if(any(self_energy .ne. self_energy)) then print*, 'NaN for jqpt', jqpt endif + deallocate(intermediate) deallocate(ifc3, d3, d3_pols, selfnrg) - deallocate(pols_k, pols_mk_mq, intermediate) + deallocate(pols_k, pols_mk_mq) deallocate(pols_k2, pols_mk_mq2) deallocate(kpt, mkpt) deallocate(w2_k, w2_mk_mq, w_k, w_mk_mq) diff --git a/cellconstructor/ThermalConductivity.py b/cellconstructor/ThermalConductivity.py index 98090917..0e9343c2 100644 --- a/cellconstructor/ThermalConductivity.py +++ b/cellconstructor/ThermalConductivity.py @@ -490,13 +490,18 @@ def construct_symmetry_matrix(rotation, translation, qvec, pos, atom_map, cell): rx1 = [] for iat in range(len(pos)): for jat in range(len(pos)): - if(iat == atom_map[jat]): + if(jat == atom_map[iat]): r = np.zeros_like(pos[iat]) phase = 0.0 r = pos[iat] - np.dot(rotation, pos[jat]) - translation rx = np.dot(r, np.linalg.inv(cell)) rx1.append(rx) if(np.linalg.norm(rx - np.rint(rx))>1.0e-6): + print(rotation) + print(translation) + print(pos[iat], pos[jat]) + print(atom_map) + print(rx) raise RuntimeError('The atom is translated different than the translation vector!') #else: # print(rx) @@ -625,6 +630,7 @@ def set_kpoints_spglib(self): translations = symmetry_dataset['translations'] print('Spacegroup: ' + symmetry_dataset['international'] + str(symmetry_dataset['number'])) self.atom_map = np.zeros((len(rotations), len(cell[1])), dtype=int) + nrot = len(rotations) for irot in range(nrot): found = [False for jat in range(len(cell[1]))]