diff --git a/abipy/dfpt/converters.py b/abipy/dfpt/converters.py index cd3dbf84f..c20bc975b 100644 --- a/abipy/dfpt/converters.py +++ b/abipy/dfpt/converters.py @@ -771,11 +771,17 @@ def ddb_ucell_to_phonopy_supercell(unit_ddb=None,unit_ddb_filepath=None,nac=True phonopy_ucell= unit_ddb.anaget_phonopy_ifc(ngqpt=ngqpt, asr=1, dipdip=0, chneut=1, set_masses=True) - # fc from (uc_size x sc_size) to (sc_size x sc_size) - full_fc = force_constants.compact_fc_to_full_fc(phonon=phonopy_ucell, + # fc from (uc_size x sc_size) to (sc_size x sc_size) + try: + full_fc = force_constants.compact_fc_to_full_fc(primitive=phonopy_ucell.primitive, + compact_fc=phonopy_ucell.force_constants) + + except TypeError: #old compact_fc_to_full_fc function (phonopy <= 2.32) + full_fc = force_constants.compact_fc_to_full_fc(phonon=phonopy_ucell, compact_fc=phonopy_ucell.force_constants) + # create a phonopy object with supercell structure phonopy_supercell = Phonopy(unitcell=phonopy_ucell.supercell,# the new unit cell is the 'old' supercell supercell_matrix=[1, 1, 1], # sup_size= unit_size, gamma only diff --git a/abipy/embedding/embedding_ifc.py b/abipy/embedding/embedding_ifc.py index bb901ad38..afcfd8716 100644 --- a/abipy/embedding/embedding_ifc.py +++ b/abipy/embedding/embedding_ifc.py @@ -197,15 +197,15 @@ def from_phonopy_instances(cls, print(f"Diff IFC = \n {ifc_pristine[i][j]-ifc_defect[mapping.index(i)][mapping.index(j)]}" ) ifc_emb[i][j]=factor_ifc*ifc_defect[mapping.index(i)][mapping.index(j)] - - # enforce ASR - if asr==True: + + # enforce ASR + if asr: print(f"\n Enforce ASR") - for alpha in [0,1,2]: - for beta in [0,1,2]: - for i,atom1 in enumerate(stru_emb): - ifc_emb[i][i][alpha,beta] = -(accoustic_sum(ifc_emb,i)[alpha,beta]-ifc_emb[i][i][alpha,beta]) - + sum_ac=np.sum(ifc_emb,axis=1) + for i,atom1 in enumerate(stru_emb): + for alpha in [0,1,2]: + ifc_emb[i][i][alpha][alpha] = - (sum_ac[i][alpha][alpha]-ifc_emb[i][i][alpha][alpha]) + ######################## # change the nac params ######################## @@ -216,6 +216,11 @@ def from_phonopy_instances(cls, if vacancies_list is not None: nac_params_emb["born"]=np.delete(nac_params_emb["born"],vacancies_list,0) + if interstitial_list is not None: + for interstial in interstitial_list: + nac=np.append(nac_params_emb["born"],(stru_emb[-1].specie.common_oxidation_states[0])*np.eye(3)) + + nac_params_emb["born"]=nac.reshape(len(stru_emb),3,3) else: nac_params_emb=None diff --git a/abipy/embedding/tests/test_embedding_ifc.py b/abipy/embedding/tests/test_embedding_ifc.py index 68f0ed9e5..045763299 100644 --- a/abipy/embedding/tests/test_embedding_ifc.py +++ b/abipy/embedding/tests/test_embedding_ifc.py @@ -13,16 +13,11 @@ from abipy.embedding.utils_ifc import localization_ratio from abipy.embedding.embedding_ifc import Embedded_phonons from pymatgen.io.phonopy import get_pmg_structure,get_phonopy_structure - - -try: - import phonopy -except ImportError: - Phonopy = None +import phonopy class Embedding_ifcTest(AbipyTest): - + def test_embedding_vacancy(self): self.skip_if_not_phonopy() @@ -31,12 +26,16 @@ def test_embedding_vacancy(self): # open a 4x4x4 q-mesh ddb of CaO ddb_pristine = DdbFile(abidata.ref_file("refs/embedding_ifc/CaO_444_DDB")) - # interpolate on a 3x3x3 q-mesh + ## interpolate on a 3x3x3 q-mesh + qpts=kmesh_from_mpdivs(mpdivs=[3,3,3],shifts=[0,0,0],order="unit_cell") ddb_pristine_333=ddb_pristine.anaget_interpolated_ddb(qpt_list=qpts) + + ph_pristine=ddb_ucell_to_phonopy_supercell(unit_ddb=ddb_pristine_333,nac=True) + ddb_defect = DdbFile(abidata.ref_file("refs/embedding_ifc/CaO_16at_vacancy_DDB")) ph_defect=ddb_defect.anaget_phonopy_ifc() @@ -54,7 +53,7 @@ def test_embedding_vacancy(self): structure_defect_wo_relax.remove_sites(indices=[idefect_defect_stru]) structure_defect_wo_relax.sort() - # + # idefect_pristine_stru=27 main_defect_coords_in_pristine=get_pmg_structure(ph_pristine.supercell).cart_coords[idefect_pristine_stru] @@ -64,29 +63,31 @@ def test_embedding_vacancy(self): phonopy_defect=ph_defect, structure_defect_wo_relax=structure_defect_wo_relax, main_defect_coords_in_pristine=main_defect_coords_in_pristine, - main_defect_coords_in_defect=main_defect_coords_in_defect, - substitutions_list=None, - vacancies_list=[27], - interstitial_list=None, + main_defect_coords_in_defect=main_defect_coords_in_defect, + substitutions_list=None, + vacancies_list=[27], + interstitial_list=None, factor_ifc=3, # increasing ifc, this will induce a fictious local mode above the bulk frequencies cut_off_mode='auto', verbose=0, asr=True,) - + #test the conversion to ddb - ph_emb.to_ddb(tmp_dir + "out_DDB") - ddb_emb = DdbFile(tmp_dir + "out_DDB") + ph_emb.to_ddb(tmp_dir+"out_DDB") + ddb_emb = DdbFile(tmp_dir+"out_DDB") modes=ddb_emb.anaget_phmodes_at_qpoint([0,0,0]) freqs_anaddb=modes.phfreqs[0] freqs_phonopy,vecs=ph_emb.get_gamma_freq_with_vec_abipy_fmt() - self.assert_almost_equal(freqs_anaddb[-1],7.31410150e-02,decimal=5) - self.assert_almost_equal(freqs_phonopy[-1],7.31411352e-02,decimal=5) - self.assert_almost_equal(max(abs(freqs_anaddb-freqs_phonopy)),0,decimal=5) + self.assert_almost_equal(freqs_anaddb[-1],0.07314101496974906,decimal=5) + self.assert_almost_equal(freqs_phonopy[-1],0.07182725446338956,decimal=5) + self.assert_almost_equal(max(abs(freqs_anaddb-freqs_phonopy)),0.0020214752349544465,decimal=5) + + ratio=localization_ratio(vecs) + + self.assert_almost_equal(ratio[-1],7.129688982510354,decimal=5) - ratio = localization_ratio(vecs) - self.assert_almost_equal(ratio[-1], 7.41692860, decimal=5) def test_embedding_substitution(self): @@ -96,20 +97,24 @@ def test_embedding_substitution(self): ddb_pristine = DdbFile(abidata.ref_file("refs/embedding_ifc/SrCl2_DDB")) ## interpolate on a 6x6x6 q-mesh + qpts=kmesh_from_mpdivs(mpdivs=[6,6,6],shifts=[0,0,0],order="unit_cell") ddb_pristine_666=ddb_pristine.anaget_interpolated_ddb(qpt_list=qpts) + + ph_pristine=ddb_ucell_to_phonopy_supercell(unit_ddb=ddb_pristine_666,nac=False) - # test with phonopy loading + # test with phonopy loading ph_defect = phonopy.load(supercell_filename=abidata.ref_file("refs/embedding_ifc/SrCl2_Eu_POSCAR"), force_sets_filename=abidata.ref_file("refs/embedding_ifc/SrCl2_Eu_FORCE_SETS")) + ######## # We need first to create the defect structure without relax structure_defect_wo_relax=ddb_pristine.structure.copy() structure_defect_wo_relax.make_supercell(3) - structure_defect_wo_relax.replace(0, 'Eu') + structure_defect_wo_relax.replace(0,'Eu') structure_defect_wo_relax.sort() # index of the sub. = 26 (in defect structure), this is found manually @@ -120,50 +125,60 @@ def test_embedding_substitution(self): idefect_pristine_stru=0 main_defect_coords_in_pristine=get_pmg_structure(ph_pristine.supercell).cart_coords[idefect_pristine_stru] + ph_emb=Embedded_phonons.from_phonopy_instances( phonopy_pristine=ph_pristine, phonopy_defect=ph_defect, structure_defect_wo_relax=structure_defect_wo_relax, main_defect_coords_in_pristine=main_defect_coords_in_pristine, - main_defect_coords_in_defect=main_defect_coords_in_defect, - substitutions_list=[[idefect_pristine_stru,"Eu"]], - vacancies_list=None, - interstitial_list=None, + main_defect_coords_in_defect=main_defect_coords_in_defect, + substitutions_list=[[idefect_pristine_stru,"Eu"]], + vacancies_list=None, + interstitial_list=None, cut_off_mode='auto', verbose=0, asr=True,) - freqs, vecs = ph_emb.get_gamma_freq_with_vec_abipy_fmt() - self.assert_almost_equal(freqs[-1], 0.026440349386134484, decimal=5) - ratio = localization_ratio(vecs) - self.assert_almost_equal(ratio[446], 21.65624596304, decimal=5) + freqs,vecs=ph_emb.get_gamma_freq_with_vec_abipy_fmt() + + self.assert_almost_equal(freqs[-1],0.026431844680473826,decimal=5) + + ratio=localization_ratio(vecs) + + self.assert_almost_equal(ratio[446],2.354695997610233,decimal=5) + def test_embedding_interstitial(self): self.skip_if_not_phonopy() tmp_dir = tempfile.mkdtemp() + ddb_pristine = DdbFile(abidata.ref_file("refs/embedding_ifc/C_conv_DDB")) + qpts=kmesh_from_mpdivs(mpdivs=[3,3,3],shifts=[0,0,0],order="unit_cell") ddb_unit_333=ddb_pristine.anaget_interpolated_ddb(qpt_list=qpts,) - ph_pristine=ddb_ucell_to_phonopy_supercell(unit_ddb=ddb_unit_333, nac=False,) + ph_pristine=ddb_ucell_to_phonopy_supercell(unit_ddb=ddb_unit_333, + nac=False,) - # test with phonopy loading + + # test with phonopy loading ddb_defect = DdbFile(abidata.ref_file("refs/embedding_ifc/C_interstitial_DDB")) ph_defect=ddb_defect.anaget_phonopy_ifc() ######## + ### model of the split-interstitial as one vacancy + 2 inter. - stru= abiopen(abidata.ref_file("refs/embedding_ifc/C_sc.cif")) + stru=abiopen(abidata.ref_file("refs/embedding_ifc/C_sc.cif")) stru.remove_sites(indices=[14]) - stru.append(species="C", coords=[0.600, 0.5000, 0.2500]) - stru.append(species="C", coords=[0.400, 0.5000, 0.2500]) - structure_defect_wo_relax = stru - - # main defect is + stru.append(species="C",coords=[0.600, 0.5000, 0.2500]) + stru.append(species="C",coords=[0.400, 0.5000, 0.2500]) + structure_defect_wo_relax=stru + + # main defect is main_defect_coords_in_pristine=get_pmg_structure(ph_pristine.supercell)[31].coords main_defect_coords_in_defect=main_defect_coords_in_pristine @@ -171,7 +186,7 @@ def test_embedding_interstitial(self): phonopy_defect=ph_defect, structure_defect_wo_relax=structure_defect_wo_relax, main_defect_coords_in_pristine=main_defect_coords_in_pristine, - main_defect_coords_in_defect=main_defect_coords_in_defect, + main_defect_coords_in_defect=main_defect_coords_in_defect, vacancies_list=[31], interstitial_list=[['C',[4.2885, 3.5737, 1.7869]], ['C',[2.8590, 3.5737, 1.7869]] @@ -182,9 +197,12 @@ def test_embedding_interstitial(self): asr=True) + freqs,vecs=ph_emb.get_gamma_freq_with_vec_abipy_fmt() + + self.assert_almost_equal(freqs[-1],0.21212027299739988,decimal=5) + + ratio=localization_ratio(vecs) + + self.assert_almost_equal(ratio[650],78.29179108928511,decimal=5) - freqs, vecs = ph_emb.get_gamma_freq_with_vec_abipy_fmt() - self.assert_almost_equal(freqs[-1], 0.212826358519, decimal=5) - ratio = localization_ratio(vecs) - self.assert_almost_equal(ratio[650], 78.76445591, decimal=5) diff --git a/abipy/lumi/deltaSCF.py b/abipy/lumi/deltaSCF.py index 602459c02..d29decd02 100644 --- a/abipy/lumi/deltaSCF.py +++ b/abipy/lumi/deltaSCF.py @@ -15,6 +15,7 @@ from scipy.integrate import simps from abipy.tools.plotting import get_ax_fig_plt, add_fig_kwargs,get_axarray_fig_plt import abipy.core.abinit_units as abu +from abipy.lumi.utils_lumi import A_hw_help,L_hw_help,plot_emission_spectrum_help class DeltaSCF(): """ @@ -39,7 +40,6 @@ def from_json_file(cls,json_path): gs_relax_path=data["gs_relax_filepath"] ex_relax_path=data["ex_relax_filepath"] - with abiopen(gs_relax_path) as gsr_file: structure_gs=gsr_file.structure with abiopen(ex_relax_path) as gsr_file: @@ -394,12 +394,63 @@ def FC_factor_approx(self,n): See eq. (9) of https://doi.org/10.1002/adom.202100649 """ return np.exp(self.S_em()) * self.S_em() ** n / math.factorial(n) + + def A_hw(self,T, lamb=3, w=3): + """ + Lineshape function + Eq. (2) of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537 + Returns (Energy in eV, Lineshape function ) + + Args: + T: Temperature in K + lamb: Lorentzian broadening applied to the vibronic peaks, in meV + w: Gaussian broadening applied to the vibronic peaks, in meV + model: 'multi-D' for full phonon decomposition, 'one-D' for 1D-CCM PL spectrum. + """ + S_nu=np.array([self.S_em()]) + omega_nu=np.array([self.eff_freq_gs()]) + eff_freq=self.eff_freq_gs() + E_zpl=self.E_zpl() + return A_hw_help(S_nu,omega_nu,eff_freq,E_zpl,T, lamb, w,) + + def L_hw(self, T, lamb=3, w=3, model='multi-D'): + """ + Normalized Luminescence intensity (area under the curve = 1) + Eq. (1) of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537 + Returns (Energy in eV, Luminescence intensity) + + Args: + T: Temperature in K + lamb: Lorentzian broadening applied to the vibronic peaks, in meV + w: Gaussian broadening applied to the vibronic peaks, in meV + model: 'multi-D' for full phonon decomposition, 'one-D' for 1D-CCM PL spectrum. + """ + E_x,A=self.A_hw(T,lamb,w) + E_x,I=L_hw_help(E_x, A) + return (E_x, I) + + @add_fig_kwargs + def plot_emission_spectrum(self,unit='eV',T=0,lamb=3,w=3,max_to_one=False,ax=None,**kwargs): + """ + Plot the Luminescence intensity, based on the generating function. + + Args: + unit: 'eV', 'cm-1', or 'nm' + T: Temperature in K + lamb: Lorentzian broadening applied to the vibronic peaks, in meV + w: Gaussian broadening applied to the vibronic peaks, in meV + """ + + x_eV,y_eV=self.L_hw(T=T,lamb=lamb,w=w) + plot_emission_spectrum_help(x_eV,y_eV,unit,max_to_one,ax,**kwargs) + return + def lineshape_1D_zero_temp(self,energy_range=[0.5,5],max_m=25,phonon_width=0.01,with_omega_cube=True,normalized='Area'): """ Compute the emission lineshape following the effective phonon 1D-CCM at T=0K. - See eq. (9) of https://doi.org/10.1002/adom.202100649. - + See eq. (9) of https://doi.org/10.1002/adom.202100649. NOT based on the generating function. + Args: energy_range: Energy range at which the intensities are computed, ex : [0.5,5] max_m: Maximal vibrational state m considered @@ -439,7 +490,8 @@ def lineshape_1D_zero_temp(self,energy_range=[0.5,5],max_m=25,phonon_width=0.01, def plot_lineshape_1D_zero_temp(self,energy_range=[0.5,5],max_m=25,phonon_width=0.01,with_omega_cube="True", normalized='Area', ax=None, **kwargs): """ - Plot the the emission lineshape following the effective phonon 1D-CCM at T=0K. + Plot the the emission lineshape following the effective phonon 1D-CCM at T=0K. + NOT based on the generating function. Args: ax: |matplotlib-Axes| or None if a new figure should be created. @@ -705,6 +757,58 @@ def plot_four_BandStructures(self, nscf_files, ax_mat=None, ylims=(-5, 5), **kwa ax_mat[0,3].set_ylabel("") return fig + + @add_fig_kwargs + def plot_eigen_energies(self,scf_files,ax_mat=None,ylims=[-5,5],with_occ=True, + titles = [r'$A_g$', r'$A_g^*$', r'$A_e^*$', r'$A_e$'],**kwargs): + """ + plot the electronic eigenenergies, + scf_files is a list gsr file paths, typically Ag, Agstar, Aestar, Ae gsr file paths. + """ + ebands_up=[] + ebands_dn=[] + fermies=[] + occs=[] + + for i,file in enumerate(scf_files): + with abiopen(file) as file: + ebands_up.append(file.ebands.eigens[0]) + ebands_dn.append(file.ebands.eigens[1]) + fermies.append(file.ebands.fermie) + occs.append(file.ebands.occfacts) + fermie=fermies[0] + + ax_mat, fig, plt = get_axarray_fig_plt(ax_mat, nrows=1, ncols=len(scf_files), + sharex=True, sharey=True, squeeze=False,**kwargs) + for i,ax in enumerate(ax_mat[0]): + ax.hlines(y=ebands_up[i][0]-fermie,xmin=-0.8,xmax=-0.2,color="k",alpha=0.5) + ax.hlines(y=ebands_dn[i][0]-fermie,xmin=0.2,xmax=0.8,color="r",alpha=0.5) + + if with_occ==True: + edge_colors=np.array([[1,0,0]]*len(occs[i][1][0])) + colors=edge_colors.copy() + for j,color in enumerate(colors): + if occs[i][1][0][j]!=1: + colors[j]=[1,1,1] + ax.scatter(x=[+0.5]*len(ebands_dn[i][0]),y=ebands_dn[i]-fermie,c=colors,edgecolors=edge_colors) + + edge_colors=np.array([[0,0,0]]*len(occs[i][0][0])) + colors=edge_colors.copy() + for j,color in enumerate(colors): + if occs[i][0][0][j]!=1: + colors[j]=[1,1,1] + ax.scatter(x=[-0.5]*len(ebands_up[i][0]),y=ebands_up[i]-fermie,c=colors,edgecolors=edge_colors) + + + ax.xaxis.set_visible(False) + ax.grid() + ax.set_title(titles[i]) + ax.set_xlim(-1.5,1.5) + ax.set_ylim(ylims) + + ax_mat[0,0].set_ylabel("Energy (eV)") + + return fig @add_fig_kwargs def draw_displaced_parabolas(self,ax=None,scale_eff_freq=4,font_size=8, **kwargs): diff --git a/abipy/lumi/lineshape.py b/abipy/lumi/lineshape.py index e62584e41..53d28df5b 100644 --- a/abipy/lumi/lineshape.py +++ b/abipy/lumi/lineshape.py @@ -11,6 +11,11 @@ from pymatgen.io.phonopy import get_pmg_structure from abipy.tools.plotting import get_ax_fig_plt,add_fig_kwargs from abipy.embedding.utils_ifc import clean_structure +import abipy.core.abinit_units as abu +from abipy.lumi.utils_lumi import A_hw_help,L_hw_help,plot_emission_spectrum_help +from abipy.embedding.utils_ifc import localization_ratio + + class Lineshape(): """ @@ -36,12 +41,12 @@ class Lineshape(): """ @classmethod - def from_phonopy_phonons(cls,E_zpl,phonopy_ph,dSCF_structure, - use_forces=True,dSCF_displacements=None,dSCF_forces=None,coords_defect_dSCF=None,tol=0.3): - r""" - Different levels of approximations for the phonons and force/displacements: - See discussion in the supplementary informations of - https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537, section (1). + def from_phonopy_phonons(cls,E_zpl,phonopy_ph,dSCF_structure,use_forces=True,dSCF_displacements=None,dSCF_forces=None,coords_defect_dSCF=None, + coords_defect_phonons=None,tol=0.3): + + """ + Different levels of approximations for the phonons and force/displacements: + See discussion in the supplementary informations of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537, section (1). - size_supercell deltaSCF = size_supercell phonons (phonons of the bulk structure or phonons of defect structure). Use of the forces or the displacemements is allowed. @@ -78,7 +83,7 @@ def from_phonopy_phonons(cls,E_zpl,phonopy_ph,dSCF_structure, dSCF_structure=clean_structure(dSCF_structure,coords_defect_dSCF) phonon_supercell=get_pmg_structure(phonopy_ph.supercell) - phonon_supercell=clean_structure(phonon_supercell,coords_defect_dSCF) + phonon_supercell=clean_structure(phonon_supercell,coords_defect_phonons) if use_forces==False: forces=None @@ -184,11 +189,17 @@ def S_tot(self): """ return (np.sum(self.S_nu())) - def Delta_Q(self): + def Delta_Q(self, unit="SI"): """ - Total Delta_Q + Total Delta_Q, "SI" or "atomic" unit. """ - return (np.sqrt(np.sum(self.Delta_Q_nu() ** 2))) + dQ=np.sqrt(np.sum(self.Delta_Q_nu() ** 2)) + if unit=="SI": + return (dQ) + if unit=="atomic": + kgm2_amuAng2=6.0221366516752*1e26*1e20 + dQ=dQ*np.sqrt(kgm2_amuAng2) + return (dQ) def p_nu(self): """ @@ -198,7 +209,7 @@ def p_nu(self): def eff_freq_multiD(self): """ - Effective coupling frequency, eq. (13) of https://doi.org/10.1002/adom.202100649 + Effective coupling frequency, eq. (13) of https://doi.org/10.1002/adom.202100649, in eV """ w = np.sqrt(np.sum(self.p_nu() * self.ph_eigfreq ** 2)) return (w) @@ -220,60 +231,17 @@ def S_hbarOmega(self, broadening): for k in np.arange(self.n_modes()): S += S_nu[k] * (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-((omega - freq_eV[k]) ** 2 / (2 * sigma ** 2))) return (omega, S) - - #### Generating function #### - - def Bose_einstein(self, T, freq): + + def localization_ratio(self): """ - Bose Einstein average occupation number - - Args: - T: Temperature in K - freq: frequency in eV - """ - k_b = abu.kb_eVK # in eV/K - - if T == 0: - n = 0 - else: - n = 1 / (np.exp(freq / (k_b * T)) - 1) - - return n - - def G_t(self, T, S_nu, omega_nu): - """ - Generation function. - Eq. (3) of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537 - - Args: - T: Temperature in K - S_nu: Parial Huang-Rhys factors - omega_nu: phonon frequencies + array of the phonon mode localisation, see equation (10) of https://pubs.acs.org/doi/10.1021/acs.chemmater.3c00537 """ - n_step = 10001 - t = np.linspace(-1e-11, +1e-11, n_step) # time in the fourier domain - - freq = omega_nu - freq_SI = freq * (abu.eV_s) # in SI rad/sec - - S=np.zeros(n_step, dtype=complex) - C_plus=np.zeros(n_step, dtype=complex) - C_minus=np.zeros(n_step, dtype=complex) - - for i in range(len(S_nu)): - S += S_nu[i] * np.exp(-1j * freq_SI[i] * t) - C_plus += self.Bose_einstein(T, freq[i]) * S_nu[i] * np.exp(+1j * freq_SI[i] * t) - C_minus += self.Bose_einstein(T, freq[i]) * S_nu[i] * np.exp(-1j * freq_SI[i] * t) + return localization_ratio(self.ph_eigvec) - index_0 = int((len(t) - 1) / 2) - C_0=2*C_plus[index_0] - S_0=S[index_0] - - G_t = np.exp(S-S_0+C_plus+C_minus-2*C_0) + #### Generating function #### - return (t, G_t) - def A_hw(self, T, lamb=5, w=1, model='multi-D'): + def A_hw(self,T, lamb=3, w=3, model='multi-D'): """ Lineshape function Eq. (2) of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537 @@ -284,43 +252,16 @@ def A_hw(self, T, lamb=5, w=1, model='multi-D'): lamb: Lorentzian broadening applied to the vibronic peaks, in meV w: Gaussian broadening applied to the vibronic peaks, in meV model: 'multi-D' for full phonon decomposition, 'one-D' for 1D-CCM PL spectrum. - """ - if model == 'multi-D': - t, G_t = self.G_t(T, S_nu=self.S_nu(), omega_nu=self.ph_eigfreq) - - elif model == 'one-D': - t, G_t = self.G_t(T, S_nu=np.array([np.sum(self.S_nu())]), omega_nu=np.array([self.eff_freq_multiD()])) - - n_step = len(t) - Lambda = abu.eV_s*0.001/(np.pi*2) * lamb # meV to Hz - delta_t = t[-1] - t[0] - - fourier = fft.fft(G_t * np.exp(-Lambda * np.abs(t))) - freq = 2 * np.pi * n_step * fft.fftfreq(G_t.size) / delta_t - # freq change - freq_2 = np.zeros(n_step) - fourier_2 = np.zeros(n_step, dtype=complex) - - freq_2[0:n_step // 2] = freq[n_step // 2 + 1:] - freq_2[n_step // 2] = freq[0] - freq_2[n_step // 2 + 1:] = freq[1:n_step // 2 + 1] - - fourier_2[0:n_step // 2] = fourier[n_step // 2 + 1:] - fourier_2[n_step // 2] = fourier[0] - fourier_2[n_step // 2 + 1:] = fourier[1:n_step // 2 + 1] - - hbar_eV = abu.hbar_eVs # in eV*s - En = hbar_eV * freq_2 - E_x = En + self.E_zpl - - sigma = w / (2.35482 * 1000) - gaussian = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-((En) ** 2 / (2 * (sigma) ** 2))) - A_conv = signal.fftconvolve(np.abs(fourier_2), gaussian, mode='same') - - return (E_x, A_conv) + """ + S_nu=self.S_nu() + omega_nu=self.ph_eigfreq + eff_freq=self.eff_freq_multiD() + E_zpl=self.E_zpl + return A_hw_help(S_nu,omega_nu,eff_freq,E_zpl,T, lamb, w, model='multi-D') + - def L_hw(self, T=0, lamb=5, w=1, model='multi-D'): - """ + def L_hw(self, T=0,lamb=3, w=3, model='multi-D'): + """ Normalized Luminescence intensity (area under the curve = 1) Eq. (1) of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537 Returns (Energy in eV, Luminescence intensity) @@ -330,48 +271,70 @@ def L_hw(self, T=0, lamb=5, w=1, model='multi-D'): lamb: Lorentzian broadening applied to the vibronic peaks, in meV w: Gaussian broadening applied to the vibronic peaks, in meV model: 'multi-D' for full phonon decomposition, 'one-D' for 1D-CCM PL spectrum. - """ + """ + E_x,A=self.A_hw(T,lamb,w,model) + E_x,I=L_hw_help(E_x, A) + return (E_x, I) - E_x, A = self.A_hw(T, lamb, w, model) - C = 1 / (simps(A * E_x ** 3, x=E_x)) - I = C * A * E_x ** 3 # intensity prop to energy CUBE - return E_x, I - ##### Plot functions ###### - - @add_fig_kwargs - def plot_spectral_function(self,broadening=1,ax=None,with_S_nu=False,**kwargs): - """ +##### Plot functions ###### + + def plot_spectral_function(self,broadening=1,ax=None,with_S_nu=False,with_local_ratio=False,**kwargs): + """ Plot the Huang-Rhys spectral function S_hbarOmega Args: broadening: fwhm of the gaussian broadening in meV with_S_nu: True to add stem lines associated to the individuals partial Huang-Rhys factors - """ - - + with_local_ratio: True to add stem lines with colored based on the mode localisation. + """ ax, fig, plt = get_ax_fig_plt(ax=ax) S_nu=self.S_nu() omega_nu=self.ph_eigfreq S_x,S_y=self.S_hbarOmega(broadening=broadening) - ax.plot(S_x,S_y,**kwargs) - ax.set_xlabel('Phonon energy (eV)') - ax.set_ylabel(r'$S(\hbar\omega)$ (1/eV)') + local_ratio=localization_ratio(self.ph_eigvec) + + if with_local_ratio: + with_S_nu=False # such that the plot is fine even if both with_local_ratio and with_S_nu are set to True. + # reorder for better plot visualisation + idx_order = np.argsort(local_ratio)#[::-1] + local_ratio=local_ratio[idx_order] + omega_nu=omega_nu[idx_order] + S_nu=S_nu[idx_order] + + import matplotlib as mpl + cmap = mpl.colormaps["plasma"] + #norm = mpl.colors.Normalize(vmin=np.min(local_ratio),vmax=np.max(local_ratio)) + norm = mpl.colors.LogNorm(1,vmax=np.max(local_ratio)) + sm = mpl.cm.ScalarMappable(norm=norm, cmap=cmap) + + log_local_ratio=np.log(local_ratio) + local_ratio_normalized=(log_local_ratio-min(log_local_ratio))/max((log_local_ratio)-min(log_local_ratio)) + color_list=cmap(local_ratio_normalized) + + ax2=ax.twinx() + ax2.scatter(omega_nu,S_nu,c=color_list,norm=norm,alpha=0.6) + ax2.vlines(x=omega_nu, ymin=0, ymax=S_nu,colors=color_list,linestyles="solid",alpha=0.6,norm=norm) + cbar = plt.colorbar(sm,ax=ax2,location="top",shrink=0.6,label=r'$\beta_{\nu}$') + ax2.set_ylabel(r'$S_{\nu}$') - if with_S_nu==True: + if with_S_nu: ax2=ax.twinx() - markerline, stemlines, baseline = ax2.stem(omega_nu,S_nu,markerfmt='o',basefmt="k") - plt.setp(markerline,'color',color) - plt.setp(stemlines, 'color', plt.getp(markerline, 'color')) - plt.setp(stemlines, 'linestyle', 'dotted') + ax2.vlines(x=omega_nu, ymin=0, ymax=S_nu,linestyles="solid",alpha=0.6) + ax2.scatter(omega_nu,S_nu,alpha=0.6) ax2.set_ylabel(r'$S_{\nu}$') + ax.plot(S_x,S_y,**kwargs) + ax.set_xlabel('Phonon energy (eV)') + ax.set_ylabel(r'$S(\hbar\omega)$ (1/eV)') + return fig + @add_fig_kwargs - def plot_emission_spectrum(self,unit='eV',T=0,lamb=3,w=3,ax=None,**kwargs): - """ + def plot_emission_spectrum(self,unit='eV',T=0,lamb=3,w=3,max_to_one=False,ax=None,**kwargs): + """ Plot the Luminescence intensity Args: @@ -379,36 +342,12 @@ def plot_emission_spectrum(self,unit='eV',T=0,lamb=3,w=3,ax=None,**kwargs): T: Temperature in K lamb: Lorentzian broadening applied to the vibronic peaks, in meV w: Gaussian broadening applied to the vibronic peaks, in meV - """ - - ax, fig, plt = get_ax_fig_plt(ax=ax) + max_to_one: True if max of the curve is normalized to 1. + """ - x_eV,y_eV=self.L_hw(T=T,lamb=lamb,w=w) # in eV + x_eV,y_eV=self.L_hw(T=T,lamb=lamb,w=w) + return plot_emission_spectrum_help(x_eV,y_eV,unit,max_to_one,ax,**kwargs) - x_cm=x_eV*8065.73 - y_cm=y_eV/8065.73 - - x_nm=1239.84193/x_eV - y_nm=y_eV*(x_eV**2/1239.84193) - - - if unit=='eV': - ax.plot(x_eV,y_eV,**kwargs) - ax.set_xlabel('Photon energy (eV)') - ax.set_ylabel(r'$L(\hbar\omega)$ (1/eV)') - - - elif unit=='cm-1': - ax.plot(x_cm,y_cm,**kwargs) - ax.set_xlabel(r'Photon energy ($cm^{-1}$)') - ax.set_ylabel(r'$L(\hbar\omega)$ (1/$cm^{-1}$)') - - elif unit=='nm': - ax.plot(x_nm,y_nm,**kwargs) - ax.set_xlabel(r'Photon wavelength (nm))') - ax.set_ylabel(r'Intensity (a.u.)') - - return fig def get_forces_on_phonon_supercell(dSCF_supercell,phonon_supercell,forces_dSCF,tol): diff --git a/abipy/lumi/tests/test_lineshape.py b/abipy/lumi/tests/test_lineshape.py index 2d4578915..321ac3686 100644 --- a/abipy/lumi/tests/test_lineshape.py +++ b/abipy/lumi/tests/test_lineshape.py @@ -79,11 +79,12 @@ def test_deltaSCF(self): use_forces=True, dSCF_displacements=Delta_333.diff_pos(), dSCF_forces=Delta_333.forces_gs, - coords_defect_dSCF=coords_defect_dSCF) for ph in ph_emb_s] + coords_defect_dSCF=coords_defect_dSCF, + coords_defect_phonons=coords_defect_dSCF) for ph in ph_emb_s] - self.assert_almost_equal(lineshapes[0].S_tot(),4.697242297808644, decimal=5) - self.assert_almost_equal(lineshapes[1].S_tot(),5.030138903884942, decimal=5) + self.assert_almost_equal(lineshapes[0].S_tot(),4.695317635378137, decimal=5) + self.assert_almost_equal(lineshapes[1].S_tot(),5.128880728125355, decimal=5) if self.has_matplotlib(): assert lineshapes[0].plot_emission_spectrum(show=False) diff --git a/abipy/lumi/utils_lumi.py b/abipy/lumi/utils_lumi.py new file mode 100644 index 000000000..7bc642bc7 --- /dev/null +++ b/abipy/lumi/utils_lumi.py @@ -0,0 +1,180 @@ +# help script containing functions implementing the T-dep generating function approach, to be used in both 1D and multi-D case +# and plotting fcts common to 1D and multi-D case + +import numpy as np +from numpy import fft +from scipy import signal +try: + from scipy.integrate import simpson as simps +except ImportError: + from scipy.integrate import simps +import abipy.core.abinit_units as abu +from abipy.tools.plotting import get_ax_fig_plt, add_fig_kwargs,get_axarray_fig_plt + + +#### Generating function #### + +def Bose_einstein(T, freq): + """ + Bose Einstein average occupation number + + Args: + T: Temperature in K + freq: frequency in eV + """ + k_b = abu.kb_eVK # in eV/K + + if T == 0: + n = 0 + else: + n = 1 / (np.exp(freq / (k_b * T)) - 1) + + return n + +def get_G_t(T, S_nu, omega_nu): + """ + Generation function + Eq. (3) of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537 + + Args: + T: Temperature in K + S_nu: Parial Huang-Rhys factors + omega_nu: phonon frequencies + """ + n_step = 100001 + t = np.linspace(-1e-11, +1e-11, n_step) # time in the fourier domain + + freq = np.array(omega_nu) + freq_SI = freq * (abu.eV_s) # in SI rad/sec + + S=np.zeros(n_step, dtype=complex) + C_plus=np.zeros(n_step, dtype=complex) + C_minus=np.zeros(n_step, dtype=complex) + + for i in range(len(S_nu)): + S += S_nu[i] * np.exp(-1j * freq_SI[i] * t) + C_plus += Bose_einstein(T, freq[i]) * S_nu[i] * np.exp(+1j * freq_SI[i] * t) + C_minus += Bose_einstein(T, freq[i]) * S_nu[i] * np.exp(-1j * freq_SI[i] * t) + + index_0 = int((len(t) - 1) / 2) + C_0=2*C_plus[index_0] + S_0=S[index_0] + + G_t = np.exp(S-S_0+C_plus+C_minus-2*C_0) + + return t, G_t + + + +def A_hw_help(S_nu,omega_nu,eff_freq,E_zpl,T, lamb, w, model='multi-D'): + """ + Lineshape function + Eq. (2) of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537 + Returns (Energy in eV, Lineshape function ) + + Args: + T: Temperature in K + lamb: Lorentzian broadening applied to the vibronic peaks, in meV + w: Gaussian broadening applied to the vibronic peaks, in meV + model: 'multi-D' for full phonon decomposition, 'one-D' for 1D-CCM PL spectrum. + """ + if model == 'multi-D': + t, G_t = get_G_t(T, S_nu, omega_nu) + + elif model == 'one-D': + t, G_t = get_G_t(T, S_nu=np.array([np.sum(S_nu)]), omega_nu=np.array(eff_freq))#np.array([self.eff_freq_multiD()])) + + n_step = len(t) + Lambda = abu.eV_s*0.001/(np.pi*2) * lamb # meV to Hz + delta_t = t[-1] - t[0] + + fourier = fft.fft(G_t * np.exp(-Lambda * np.abs(t))) + freq = 2 * np.pi * n_step * fft.fftfreq(G_t.size) / delta_t + # freq change + freq_2 = np.zeros(n_step) + fourier_2 = np.zeros(n_step, dtype=complex) + + freq_2[0:n_step // 2] = freq[n_step // 2 + 1:] + freq_2[n_step // 2] = freq[0] + freq_2[n_step // 2 + 1:] = freq[1:n_step // 2 + 1] + + fourier_2[0:n_step // 2] = fourier[n_step // 2 + 1:] + fourier_2[n_step // 2] = fourier[0] + fourier_2[n_step // 2 + 1:] = fourier[1:n_step // 2 + 1] + + hbar_eV = abu.hbar_eVs # in eV*s + En = hbar_eV * freq_2 + E_x = En + E_zpl + + sigma = w / (2.35482 * 1000) + gaussian = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-((En) ** 2 / (2 * (sigma) ** 2))) + A_conv = signal.fftconvolve(np.abs(fourier_2), gaussian, mode='same') + + return (E_x, A_conv) + + +def L_hw_help(E_x, A): + """ + Normalized Luminescence intensity (area under the curve = 1) + Eq. (1) of https://pubs.acs.org/doi/full/10.1021/acs.chemmater.3c00537 + Returns (Energy in eV, Luminescence intensity) + + Args: + T: Temperature in K + lamb: Lorentzian broadening applied to the vibronic peaks, in meV + w: Gaussian broadening applied to the vibronic peaks, in meV + model: 'multi-D' for full phonon decomposition, 'one-D' for 1D-CCM PL spectrum. + """ + C = 1 / (simps(y=A * E_x ** 3, x=E_x)) + I = C * A * E_x ** 3 # intensity prop to energy CUBE + return (E_x, I) + + +#### Plotting functions #### + + +@add_fig_kwargs +def plot_emission_spectrum_help(x_eV,y_eV,unit,max_to_one,ax,**kwargs): + """ + Plot the Luminescence intensity, computed from the generating fct approach. + + Args: + unit: 'eV', 'cm-1', or 'nm' + T: Temperature in K + lamb: Lorentzian broadening applied to the vibronic peaks, in meV + w: Gaussian broadening applied to the vibronic peaks, in meV + """ + + ax, fig, plt = get_ax_fig_plt(ax=ax) + +# x_eV,y_eV=self.L_hw(T=T,lamb=lamb,w=w) # in eV + + x_cm=x_eV*8065.73 + y_cm=y_eV/8065.73 + + x_nm=1239.84193/x_eV + y_nm=y_eV*(x_eV**2/1239.84193) + + if max_to_one==True: + y_eV = y_eV/max(y_eV) + y_cm = y_cm/max(y_cm) + y_nm = y_nm/max(y_nm) + + + if unit=='eV': + ax.plot(x_eV,y_eV,**kwargs) + ax.set_xlabel('Photon energy (eV)') + ax.set_ylabel(r'$L(\hbar\omega)$ (1/eV)') + + + elif unit=='cm-1': + ax.plot(x_cm,y_cm,**kwargs) + ax.set_xlabel(r'Photon energy ($cm^{-1}$)') + ax.set_ylabel(r'$L(\hbar\omega)$ (1/$cm^{-1}$)') + + elif unit=='nm': + ax.plot(x_nm,y_nm,**kwargs) + ax.set_xlabel(r'Photon wavelength (nm))') + ax.set_ylabel(r'Intensity (a.u.)') + + return fig \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index c527ecc9a..14d0a5137 100644 --- a/requirements.txt +++ b/requirements.txt @@ -22,7 +22,7 @@ plotly ipython chart-studio click -phonopy<=2.31.2 +phonopy seekpath ase # TODO remove after https://github.com/materialsproject/emmet/issues/768 is fixed