From 6c7933e3a661afe22fabcfcbe1e3c74b051932dc Mon Sep 17 00:00:00 2001 From: xh125 Date: Sun, 26 Sep 2021 19:59:26 +0800 Subject: [PATCH] fixed a bug in sub kpoint space sh --- src_complex/getwcvk.f90 | 35 ++++++++++++----------------------- src_complex/initialsh.f90 | 5 ++++- src_complex/lvcsh.f90 | 17 +++++++++-------- src_complex/saveinf.f90 | 20 ++++++++++---------- 4 files changed, 35 insertions(+), 42 deletions(-) diff --git a/src_complex/getwcvk.f90 b/src_complex/getwcvk.f90 index e36aadc..3c61fd0 100644 --- a/src_complex/getwcvk.f90 +++ b/src_complex/getwcvk.f90 @@ -6,18 +6,17 @@ module getwcvk contains !ref : 1 S. Butscher et al., Physical Review B 72 (2005) !ref : 2 <固体物理> (9-29)(9-31) - subroutine get_Wcvk(ihband_min,ieband_max,fwhm,w_center) + subroutine get_Wcvk(ihband_min,ieband_max,nk_sub,fwhm,w_center) !得到光激发下垂直跃迁的跃迁几率 - use elph2,only : vmef,nkf !vmef(3,nbndsub,nbndsub,nkf) - use readepw,only : etf,icbm - use surfacecom,only : nk_sub,indexk + use elph2,only : vmef,nkf,etf_sub !vmef(3,nbndsub,nbndsub,nkf) + use readepw,only : icbm + use surfacecom,only : indexk use io,only : stdout use constants,only : ryd2eV,ry_to_fs implicit none - integer , intent(in) :: ihband_min,ieband_max + integer , intent(in) :: ihband_min,ieband_max,nk_sub real(kind=dp),intent(in) :: fwhm real(kind=dp),intent(in) :: w_center - real(kind=dp),allocatable :: W_cvk_(:,:,:) real(kind=dp) :: fwhm_2T2 @@ -31,38 +30,28 @@ subroutine get_Wcvk(ihband_min,ieband_max,fwhm,w_center) integer :: ivbm ivbm = icbm-1 - - allocate(W_cvk_(icbm:ieband_max,ihband_min:ivbm,nkf),stat=ierr) + + allocate(W_cvk(icbm:ieband_max,ihband_min:ivbm,nk_sub),stat=ierr) if(ierr /=0) call errore('getmcvk','Error allocating W_cvk',1) !ref : 1 S. Butscher et al., Physical Review B 72 (2005) !ref : 1 S. Fernandez-Alberti et al., The Journal of Chemical Physics 137 (2012) fwhm_2T2 = fwhm**2.0/4.0*log(2.0) - W_cvk_ = 0.0 - do ik=1,nkf + W_cvk = 0.0 + do ik=1,nk_sub do ibnd=ihband_min,ivbm do jbnd=icbm,ieband_max Evmef = 0.0 - E_mnk = etf(jbnd,ik)-etf(ibnd,ik) + E_mnk = etf_sub(jbnd,ik)-etf_sub(ibnd,ik) do ipol=1,3 - Evmef =Evmef+ (efield_cart(ipol)*(vmef(ipol,jbnd,ibnd,ik))) + Evmef =Evmef+ (efield_cart(ipol)*(vmef(ipol,jbnd,ibnd,indexk(ik)))) enddo fcw = f_w(E_mnk,w_center,fwhm_2T2) - W_cvk_(jbnd,ibnd,ik) = REAL(Evmef*CONJG(Evmef))*fcw + W_cvk(jbnd,ibnd,ik) = REAL(Evmef*CONJG(Evmef))*fcw enddo enddo enddo - allocate(W_cvk(icbm:ieband_max,ihband_min:ivbm,nk_sub),stat=ierr) - W_cvk = 0.0 - do ik=1,nk_sub - W_cvk(icbm:ieband_max,ihband_min:ivbm,ik)=W_cvk_(icbm:ieband_max,ihband_min:ivbm,indexk(ik)) - enddo - - deallocate(W_cvk_) - - - !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! !% Write laser information %! diff --git a/src_complex/initialsh.f90 b/src_complex/initialsh.f90 index 693ae19..1905047 100644 --- a/src_complex/initialsh.f90 +++ b/src_complex/initialsh.f90 @@ -472,7 +472,8 @@ subroutine set_subkspace(lelecsh,lholesh,llaser,w_laser,nk_sub) do ik=1,nk_sub etf_sub(ibndmin:ibndmax,ik) = etf(ibndmin:ibndmax,indexk(ik)) enddo - + + deallocate(etf) end subroutine set_subkspace @@ -519,6 +520,8 @@ subroutine set_subqspace(nk_sub,indexk,nq,nq_sub) wf_sub(:,iq) = wf(:,indexq(iq)) enddo + deallocate(wf) + allocate(iminusq_sub(nq_sub)) iminusq_sub = 0 do iq=1,nq_sub diff --git a/src_complex/lvcsh.f90 b/src_complex/lvcsh.f90 index 69d3e54..cadf44f 100644 --- a/src_complex/lvcsh.f90 +++ b/src_complex/lvcsh.f90 @@ -77,7 +77,7 @@ program lvcsh calculate_nonadiabatic_coupling, & calculate_hopping_probability, & convert_adiabatic_diabatic,add_decoherence - use elph2,only : wf_sub,xkf,nqtotf,nktotf,nbndfst,gmnvkq,etf + use elph2,only : wf_sub,xkf,nqtotf,nktotf,nbndfst,gmnvkq,etf,epmatq use modes,only : nmodes use cell_base,only : bg use date_and_times,only : get_date_and_time @@ -139,6 +139,7 @@ program lvcsh endif deallocate(gmnvkq) + deallocate(epmatq) call allocatesh(methodsh,lelecsh,lholesh,nk_sub,nmodes,nq_sub) @@ -147,7 +148,7 @@ program lvcsh ! set the friction coefficient of Langevin dynamica of all phonon modes. call set_gamma(nmodes,nq_sub,gamma,ld_fric,wf_sub,ld_gamma) - if(llaser) call get_Wcvk(ihband_min,ieband_max,fwhm,w_laser) + if(llaser) call get_Wcvk(ihband_min,ieband_max,nk_sub,fwhm,w_laser) ! ref : https://journals.aps.org/prb/pdf/10.1103/PhysRevB.72.045314 !get W_cvk(icband,ivband,ik) @@ -527,7 +528,7 @@ program lvcsh call save_csit(nefre,nsnap,naver,csit_e,csit_e_file) call save_wsit(nefre,nsnap,naver,wsit_e,wsit_e_file) call save_psit(nefre,nsnap,naver,psit_e,psit_e_file) - call plot_band_occupatin_withtime(neband,nk_sub,Enk_e,xkf,nsnap,psit_e,csit_e,savedsnap,band_e_file) + call plot_band_occupatin_withtime(neband,nk_sub,Enk_e,nsnap,psit_e,csit_e,savedsnap,band_e_file) endif if(lholesh) then @@ -535,7 +536,7 @@ program lvcsh call save_csit(nhfre,nsnap,naver,csit_h,csit_h_file) call save_wsit(nhfre,nsnap,naver,wsit_h,wsit_h_file) call save_psit(nhfre,nsnap,naver,psit_h,psit_h_file) - call plot_band_occupatin_withtime(nhband,nk_sub,Enk_h,xkf,nsnap,psit_h,csit_h,savedsnap,band_h_file) + call plot_band_occupatin_withtime(nhband,nk_sub,Enk_h,nsnap,psit_h,csit_h,savedsnap,band_h_file) endif @@ -611,8 +612,8 @@ program lvcsh call plot_csit(nefre,nsnap,naver,csit_e,csit_e_file) call plot_wsit(nefre,nefre_sh,nsnap,naver,wsit_e,wsit_e_file) call plot_psit(nefre,nsnap,naver,psit_e,psit_e_file) - call plot_band_occupatin_withtime(neband,nk_sub,Enk_e,xkf,nsnap,psit_e,csit_e,savedsnap,band_e_file) - call plot_eh_temp(neband,nk_sub,Enk_e,xkf,nsnap,psit_e,csit_e,savedsnap,e_temp_file) + call plot_band_occupatin_withtime(neband,nk_sub,Enk_e,nsnap,psit_e,csit_e,savedsnap,band_e_file) + call plot_eh_temp(neband,nk_sub,Enk_e,nsnap,psit_e,csit_e,savedsnap,e_temp_file) deallocate(pes_e,pes_one_e,csit_e,wsit_e,psit_e) endif @@ -623,8 +624,8 @@ program lvcsh call plot_csit(nhfre,nsnap,naver,csit_h,csit_h_file) call plot_wsit(nhfre,nhfre_sh,nsnap,naver,wsit_h,wsit_h_file) call plot_psit(nhfre,nsnap,naver,psit_h,psit_h_file) - call plot_band_occupatin_withtime(nhband,nk_sub,Enk_h,xkf,nsnap,psit_h,csit_h,savedsnap,band_h_file) - call plot_eh_temp(nhband,nk_sub,Enk_h,xkf,nsnap,psit_h,csit_h,savedsnap,h_temp_file) + call plot_band_occupatin_withtime(nhband,nk_sub,Enk_h,nsnap,psit_h,csit_h,savedsnap,band_h_file) + call plot_eh_temp(nhband,nk_sub,Enk_h,nsnap,psit_h,csit_h,savedsnap,h_temp_file) deallocate(pes_h,pes_one_h,csit_h,wsit_h,psit_h) endif diff --git a/src_complex/saveinf.f90 b/src_complex/saveinf.f90 index e6fbaf5..44844ab 100644 --- a/src_complex/saveinf.f90 +++ b/src_complex/saveinf.f90 @@ -1050,11 +1050,11 @@ &different diabatic wave function to file:",trim(psit_file_) end subroutine plot_psit - subroutine plot_band_occupatin_withtime(nband,nk,Enk,xk,nsnap,psit,csit,dsnap,band_file) + subroutine plot_band_occupatin_withtime(nband,nk,Enk,nsnap,psit,csit,dsnap,band_file) + use elph2,only : xkf implicit none integer , intent(in) :: nband,nk,dsnap real(kind=dp),intent(in) :: Enk(nband,nk) - real(kind=dp),intent(in) :: xk(3,nk) integer , intent(in) :: nsnap real(kind=dp),intent(in) :: psit(1:nband*nk,0:nsnap),csit(1:nband*nk,0:nsnap) character(len=*),intent(in) :: band_file @@ -1075,9 +1075,9 @@ subroutine plot_band_occupatin_withtime(nband,nk,Enk,xk,nsnap,psit,csit,dsnap,ba (dt*nstep*isnap*ry_to_fs,dt*nstep*isnap*ry_to_fs,isnap=0,nsnap,dsnap) do ik=1,nk ifre = (ik-1)*nband+iband - kx = xk(1,indexk(ik)) - ky = xk(2,indexk(ik)) - kz = xk(3,indexk(ik)) + kx = xkf(1,indexk(ik)) + ky = xkf(2,indexk(ik)) + kz = xkf(3,indexk(ik)) if(kx > 0.5) kx=kx-1.0 if(ky > 0.5) ky=ky-1.0 if(kz > 0.5) kz=kz-1.0 @@ -1097,11 +1097,11 @@ subroutine plot_band_occupatin_withtime(nband,nk,Enk,xk,nsnap,psit,csit,dsnap,ba end subroutine plot_band_occupatin_withtime - subroutine plot_eh_temp(nband,nk,Enk,xk,nsnap,psit,csit,dsnap,temp_file) + subroutine plot_eh_temp(nband,nk,Enk,nsnap,psit,csit,dsnap,temp_file) + use elph2,only : xkf implicit none integer , intent(in) :: nband,nk,dsnap real(kind=dp),intent(in) :: Enk(nband,nk) - real(kind=dp),intent(in) :: xk(3,nk) integer , intent(in) :: nsnap real(kind=dp),intent(in) :: psit(1:nband*nk,0:nsnap),csit(1:nband*nk,0:nsnap) character(len=*),intent(in) :: temp_file @@ -1124,9 +1124,9 @@ subroutine plot_eh_temp(nband,nk,Enk,xk,nsnap,psit,csit,dsnap,temp_file) do isnap =0,nsnap do ik=1,nk ifre = (ik-1)*nband+iband - kx = xk(1,indexk(ik)) - ky = xk(2,indexk(ik)) - kz = xk(3,indexk(ik)) + kx = xkf(1,indexk(ik)) + ky = xkf(2,indexk(ik)) + kz = xkf(3,indexk(ik)) if(kx > 0.5) kx=kx-1.0 if(ky > 0.5) ky=ky-1.0 if(kz > 0.5) kz=kz-1.0