change to the subspace to do sh

src_complex/dynamics.f90

@@ -2,7 +2,7 @@ module dynamics
   use kinds,only     : dp,dpc
   use constants,only : cone,czero
   use parameters,only : temp
-  use epwcom,only         : kqmap
+  use epwcom,only         : kqmap_sub
 
   implicit none
 
@@ -28,7 +28,7 @@ subroutine set_gamma(nmodes,nq,gamma,ld_fric,wf,ld_gamma)
   end subroutine
 
   subroutine get_dEa_dQ(nmodes,nq,nband,nk,P_nk,gmnvkq,lit_gmnvkq,isurface,dEa_dQ)
-    use elph2,only : iminusq
+    use elph2,only : iminusq_sub
     implicit none
     integer,intent(in) :: nmodes,nq,nband,nk
     integer,intent(in) :: isurface
@@ -47,18 +47,20 @@ subroutine get_dEa_dQ(nmodes,nq,nband,nk,P_nk,gmnvkq,lit_gmnvkq,isurface,dEa_dQ)
     do iq=1,nq   
       do imode=1,nmodes
         do ik=1,nk
-          ikq = kqmap(ik,iq)
+          ikq = kqmap_sub(ik,iq)
+          if(ikq /= 0) then
           do iband1=1,nband
             do iband2=1,nband
               epc = gmnvkq(iband1,iband2,imode,ik,iq)
               if(epc /= czero) then
-                dEa_dQ(imode,iminusq(iq)) = dEa_dQ(imode,iminusq(iq)) + &
+                dEa_dQ(imode,iminusq_sub(iq)) = dEa_dQ(imode,iminusq_sub(iq)) + &
                 epc*CONJG(P_nk(iband1,ik,isurface))*P_nk(iband2,ikq,isurface)                 
                 !dEa_dQ(imode,iq) = dEa_dQ(imode,iq) + &
                 !epc*CONJG(P_nk(iband1,ik,isurface))*P_nk(iband2,ikq,isurface)       
               endif
             enddo
           enddo
+          endif
         enddo
       enddo
     enddo
@@ -230,7 +232,7 @@ subroutine rk4_electron_diabatic(nfre,HH,cc,cc0,dc1,dc2,dc3,dc4,tt)
 
 
   subroutine get_dEa2_dQ2(nmodes,nq,nfre,nfre_sh,isurface,EE,dd,dEa2_dQ2)
-    use elph2 , only : iminusq
+    use elph2 , only : iminusq_sub
     implicit none
     integer,intent(in) :: nmodes,nq,nfre,nfre_sh
     integer,intent(in) :: isurface
@@ -251,8 +253,8 @@ subroutine get_dEa2_dQ2(nmodes,nq,nfre,nfre_sh,isurface,EE,dd,dEa2_dQ2)
             !(EE(isurface)-EE(ifre))*(DD(ifre,isurface,imode,iq)*DD(isurface,ifre,imode,iq)+&
             !                        (DD(isurface,ifre,imode,iq)*DD(ifre,isurface,imode,iq)))            
             dEa2_dQ2(imode,iq) = dEa2_dQ2(imode,iq) + &
-            (EE(isurface)-EE(ifre))*(DD(ifre,isurface,imode,iq)*DD(isurface,ifre,imode,iminusq(iq))+&
-                                    (DD(isurface,ifre,imode,iq)*DD(ifre,isurface,imode,iminusq(iq))))
+            (EE(isurface)-EE(ifre))*(DD(ifre,isurface,imode,iq)*DD(isurface,ifre,imode,iminusq_sub(iq))+&
+                                    (DD(isurface,ifre,imode,iq)*DD(ifre,isurface,imode,iminusq_sub(iq))))
           endif
         enddo
       enddo

src_complex/elph2.f90

@@ -39,6 +39,7 @@ MODULE elph2
     igk(:),                  &!  Index for k+G vector
     igkq(:),                 &!  Index for k+q+G vector
     iminusq(:),              &!  Index for -q vector
+    iminusq_sub(:),          &!
     iminusk(:),              &!  Index for -k vector
     igk_k_all(:, :),         &!  Global index (in case of parallel)
     ngk_all(:),              &!  Global number of plane wave for each global k-point
@@ -73,9 +74,11 @@ MODULE elph2
     xqf(:, :),               &!  fine q point grid
     wqf(:),                  &!  weights on the fine q grid
     etf(:, :),               &!  interpolated eigenvalues (nbnd, nkqf)
+    etf_sub(:,:),            &!
     etf_k(:, :),             &!  Saved interpolated KS eigenenergies for later used in q-parallelization (nbnd, nkqf)
     etf_ks(:, :),            &!  interpolated eigenvalues (nbnd, nkqf) KS eigenvalues in the case of eig_read
     wf(:, :),                &!  interpolated eigenfrequencies
+    wf_sub(:,:),             &!
     gmnvkq(:,:,:,:,:),          &!  !! g vectex
     gamma_all(:, :, :, :),      &!  Gamma
     gamma_nest(:, :),        &!  Nesting function in the case of q-parallelization

src_complex/epwcom.f90

@@ -420,7 +420,9 @@ module klist_epw
     use kinds, only :dp  
     implicit none
     integer, allocatable :: kqmap(:,:)
-    !!!! map of k+q grid into k grid
+    !! map of k+q grid into k grid
+    integer, allocatable :: kqmap_sub(:,:)
+    !! 
     integer, allocatable :: isk_all(:)
     !! Spin index of each k-point (used in LSDA calculations only)
     INTEGER, ALLOCATABLE :: isk_loc(:)

src_complex/getwcvk.f90

@@ -10,12 +10,14 @@ subroutine get_Wcvk(ihband_min,ieband_max,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 io,only : stdout
     use constants,only : ryd2eV,ry_to_fs
     implicit none
     integer , intent(in) :: ihband_min,ieband_max
     real(kind=dp),intent(in) :: fwhm
     real(kind=dp),intent(in) :: w_center
+    real(kind=dp),allocatable :: W_cvk_(:,:,:)
 
     real(kind=dp) :: fwhm_2T2
 
@@ -30,13 +32,13 @@ subroutine get_Wcvk(ihband_min,ieband_max,fwhm,w_center)
 
     ivbm = icbm-1
 
-    allocate(W_cvk(icbm:ieband_max,ihband_min:ivbm,nkf),stat=ierr)
+    allocate(W_cvk_(icbm:ieband_max,ihband_min:ivbm,nkf),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
+    W_cvk_ = 0.0
     do ik=1,nkf
       do ibnd=ihband_min,ivbm
         do jbnd=icbm,ieband_max
@@ -46,11 +48,22 @@ subroutine get_Wcvk(ihband_min,ieband_max,fwhm,w_center)
             Evmef =Evmef+ (efield_cart(ipol)*(vmef(ipol,jbnd,ibnd,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/hamiltonian.f90

@@ -8,7 +8,7 @@ module hamiltonian
   use readepw,only : E_nk
 	use constants,only : ryd2mev,cone,czero
   use surfacecom,only: lelecsh,lholesh,ieband_min,ieband_max,&
-                       ihband_min,ihband_max
+                       ihband_min,ihband_max,nk_sub,nq_sub
 	use memory_report,only : MB,GB,complex_size, real_size,int_size,ram,print_memory  
 
   implicit none
@@ -40,16 +40,16 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
     integer,intent(in) :: ieband_min,ieband_max,ihband_min,ihband_max
     if(lelecsh) then
       neband = ieband_max - ieband_min + 1
-      nefre   = neband * nk
-			allocate(Enk_e(neband,nk),stat=ierr,errmsg=msg)
+      nefre   = neband * nk_sub
+			allocate(Enk_e(neband,nk_sub),stat=ierr,errmsg=msg)
 			if(ierr /= 0) call io_error(msg)
-      allocate(gmnvkq_e(neband,neband,nmodes,nk,nq),stat=ierr,errmsg=msg)
+      allocate(gmnvkq_e(neband,neband,nmodes,nk_sub,nq_sub),stat=ierr,errmsg=msg)
       if(ierr /=0) then
 				call errore('hamiltonian','Error allocating gmnvkq_e',1)
 				call io_error(msg)
 			endif
 			gmnvkq_e = 0.0
-			ram = complex_size*neband*neband*nmodes*nk*nq
+			ram = complex_size*neband*neband*nmodes*nk_sub*nq_sub
 			call print_memory("gmnvkq_e",ram)
       allocate(H0_e(nefre,nefre),stat=ierr,errmsg=msg)
       if(ierr /=0) then
@@ -59,7 +59,7 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
 			H0_e = 0.0
 			ram = complex_size*nefre*nefre
 			call print_memory("H0_e",ram)
-      allocate(H0_e_nk(neband,nk,neband,nk),stat=ierr,errmsg=msg)
+      allocate(H0_e_nk(neband,nk_sub,neband,nk_sub),stat=ierr,errmsg=msg)
       if(ierr /=0) then
 				call errore('hamiltonian','Error allocating H0_e_nk',1)
 				call io_error(msg)
@@ -73,7 +73,7 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
 			endif
 			H_e = 0.0
 			call print_memory("H_e",ram)
-      allocate(H_e_nk(neband,nk,neband,nk),stat=ierr,errmsg=msg)
+      allocate(H_e_nk(neband,nk_sub,neband,nk_sub),stat=ierr,errmsg=msg)
       if(ierr /=0) then
 				call errore('hamiltonian','Error allocating H_e_nk',1)
 				call io_error(msg)
@@ -85,16 +85,16 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
 
     if(lholesh) then
       nhband = ihband_max - ihband_min + 1
-      nhfre   = nhband * nk
-			allocate(Enk_h(nhband,nk),stat=ierr,errmsg=msg)
+      nhfre   = nhband * nk_sub
+			allocate(Enk_h(nhband,nk_sub),stat=ierr,errmsg=msg)
 			if(ierr /= 0) call io_error(msg)
-      allocate(gmnvkq_h(nhband,nhband,nmodes,nk,nq),stat=ierr,errmsg=msg)
+      allocate(gmnvkq_h(nhband,nhband,nmodes,nk_sub,nq_sub),stat=ierr,errmsg=msg)
       if(ierr /=0) then
 				call errore('hamiltonian','Error allocating gmnvkq_h',1)
 				call io_error(msg)
 		  endif
 			gmnvkq_h = 0.0
-			ram = complex_size*nhband*nhband*nmodes*nk*nq
+			ram = complex_size*nhband*nhband*nmodes*nk_sub*nq_sub
 			call print_memory("gmnvkq_h",ram)
       allocate(H0_h(nhfre,nhfre),stat=ierr,errmsg=msg)
       if(ierr /=0) then
@@ -104,7 +104,7 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
 			H0_h = 0.0
 			ram = complex_size*nhfre*nhfre
 			call print_memory("H0_h",ram)			
-      allocate(H0_h_nk(nhband,nk,nhband,nk),stat=ierr,errmsg=msg)
+      allocate(H0_h_nk(nhband,nk_sub,nhband,nk_sub),stat=ierr,errmsg=msg)
       if(ierr /=0) then
 				call errore('hamiltonian','Error allocating H0_h_nk',1)
 				call io_error(msg)
@@ -119,7 +119,7 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
 			H_h = 0.0
 			ram = complex_size*nhfre*nhfre
 			call print_memory("H_h",ram)
-      allocate(H_h_nk(nhband,nk,nhband,nk),stat=ierr,errmsg=msg)
+      allocate(H_h_nk(nhband,nk_sub,nhband,nk_sub),stat=ierr,errmsg=msg)
       if(ierr /=0) then
 				call errore('hamiltonian','Error allocating H_h_nk',1)
 				call io_error(msg)
@@ -131,14 +131,14 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
 
   end subroutine allocate_hamiltonian
 
-  subroutine set_H0_nk(nk,ibandmin,ibandmax,Enk,H0_nk,gmnvkq_eh)
-    use elph2, only : etf,gmnvkq
+  subroutine set_H0_nk(nk,nk_sub,indexk,nq_sub,indexq,ibandmin,ibandmax,Enk,H0_nk,gmnvkq_eh)
+    use elph2, only : etf_sub,gmnvkq
     implicit none
-    integer , intent(in) :: nk
-    integer , intent(in) :: ibandmin,ibandmax
-    real(kind=dp), intent(out) :: Enk(ibandmax-ibandmin+1,nk)
-		complex(kind=dp), intent(out) :: H0_nk(ibandmax-ibandmin+1,nk,ibandmax-ibandmin+1,nk)
-    complex(kind=dp), intent(out) :: gmnvkq_eh(ibandmax-ibandmin+1,ibandmax-ibandmin+1,nmodes,nk,nq)
+    integer , intent(in) :: nk,nk_sub,nq_sub
+    integer , intent(in) :: indexk(nk_sub),indexq(nq_sub),ibandmin,ibandmax
+    real(kind=dp), intent(out) :: Enk(ibandmax-ibandmin+1,nk_sub)
+		complex(kind=dp), intent(out) :: H0_nk(ibandmax-ibandmin+1,nk_sub,ibandmax-ibandmin+1,nk_sub)
+    complex(kind=dp), intent(out) :: gmnvkq_eh(ibandmax-ibandmin+1,ibandmax-ibandmin+1,nmodes,nk_sub,nq_sub)
 
 		integer :: nband
     integer :: ik,iband
@@ -147,15 +147,20 @@ subroutine set_H0_nk(nk,ibandmin,ibandmax,Enk,H0_nk,gmnvkq_eh)
 		nband = ibandmax-ibandmin+1
 
     H0_nk = 0.0
-    do ik=1,nk
+    do ik=1,nk_sub
       do iband=1,nband
-				Enk(iband,ik) = etf(iband+ibandmin-1,ik)
-        H0_nk(iband,ik,iband,ik)=etf(iband+ibandmin-1,ik)*cone
+				Enk(iband,ik) = etf_sub(iband+ibandmin-1,ik)
+        H0_nk(iband,ik,iband,ik)=etf_sub(iband+ibandmin-1,ik)*cone
       enddo
     enddo  
-
-    gmnvkq_eh = epmatq(ibandmin:ibandmax,ibandmin:ibandmax,:,:,:)    
-
+
+    do iq=1,nq_sub
+      do ik=1,nk_sub
+        gmnvkq_eh(:,:,:,ik,iq) = epmatq(ibandmin:ibandmax,ibandmin:ibandmax,:,indexk(ik),indexq(iq))    
+      enddo
+    enddo
+
+
   end subroutine set_H0_nk
 
 
@@ -164,7 +169,7 @@ end subroutine set_H0_nk
   !ref: 1 F. Giustino, Reviews of Modern Physics 89 (2017) 015003.
   !     (1)  
   subroutine set_H_nk(nband,nk,nmodes,nq,ph_Q,gmnvkq_eh,H0_nk,H_nk)
-    use epwcom,only  : kqmap
+    use epwcom,only  : kqmap_sub
     implicit none
     integer , intent(in) :: nband,nk,nmodes,nq
     complex(kind=dpc),intent(in) :: ph_Q(nmodes,nq)
@@ -179,15 +184,17 @@ subroutine set_H_nk(nband,nk,nmodes,nq,ph_Q,gmnvkq_eh,H0_nk,H_nk)
 
     do iq=1,nq
       do ik=1,nk
-        ikq=kqmap(ik,iq)
-				do nu=1,nmodes					
-          do iband=1,nband !|iband,ik>
-            do jband=1,nband !|jband,ikq>
-              H_nk(jband,ikq,iband,ik) = H_nk(jband,ikq,iband,ik)+&
-              gmnvkq_eh(iband,jband,nu,ik,iq)*ph_Q(nu,iq)
+        ikq=kqmap_sub(ik,iq)
+        if(ikq/=0) then
+          do nu=1,nmodes					
+            do iband=1,nband !|iband,ik>
+              do jband=1,nband !|jband,ikq>
+                H_nk(jband,ikq,iband,ik) = H_nk(jband,ikq,iband,ik)+&
+                gmnvkq_eh(iband,jband,nu,ik,iq)*ph_Q(nu,iq)
+              enddo
             enddo
           enddo
-        enddo
+        endif
       enddo
     enddo