fixed a bug in sub kpoint space sh

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            %!

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

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,15 +528,15 @@ 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
       call save_pes(nhfre,nsnap,naver,pes_one_h,pes_h,pes_h_file)
       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
 

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