change to the full kpoint and qpoing space

src_complex/dynamics.f90

@@ -1,8 +1,8 @@
 module dynamics
   use kinds,only     : dp,dpc
   use constants,only : cone,czero
-  use parameters,only: temp
-  use epwcom,only    : kqmap_sub
+  use parameters,only : temp
+  use epwcom,only         : kqmap
 
   implicit none
 
@@ -27,11 +27,12 @@ subroutine set_gamma(nmodes,nq,gamma,ld_fric,wf,ld_gamma)
 
   end subroutine
 
-  subroutine get_dEa_dQ(nmodes,nq,nband,nk,P_nk,gmnvkq,isurface,dEa_dQ)
-    use elph2,only : iminusq_sub
+  subroutine get_dEa_dQ(nmodes,nq,nband,nk,P_nk,gmnvkq,lit_gmnvkq,isurface,dEa_dQ)
+    use elph2,only : iminusq
     implicit none
     integer,intent(in) :: nmodes,nq,nband,nk
     integer,intent(in) :: isurface
+    real(kind=dp),intent(in)      :: lit_gmnvkq
     complex(kind=dpc),intent(in)  :: P_nk(nband,nk,nband*nk)
     complex(kind=dpc),intent(in)  :: gmnvkq(nband,nband,nmodes,nk,nq)
     complex(kind=dpc),intent(out) :: dEa_dQ(nmodes,nq)
@@ -46,20 +47,18 @@ subroutine get_dEa_dQ(nmodes,nq,nband,nk,P_nk,gmnvkq,isurface,dEa_dQ)
     do iq=1,nq   
       do imode=1,nmodes
         do ik=1,nk
-          ikq = kqmap_sub(ik,iq)
-          if(ikq /= 0) then
+          ikq = kqmap(ik,iq)
           do iband1=1,nband
             do iband2=1,nband
               epc = gmnvkq(iband1,iband2,imode,ik,iq)
               if(epc /= czero) then
-                dEa_dQ(imode,iminusq_sub(iq)) = dEa_dQ(imode,iminusq_sub(iq)) + &
+                dEa_dQ(imode,iminusq(iq)) = dEa_dQ(imode,iminusq(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
@@ -137,7 +136,7 @@ subroutine rk4_nuclei(nmodes,nq,dEa_dQ,ld_gamma,wf,xx,vv,tt)
   ! ref : 1 J. Qiu, X. Bai, and L. Wang, The Journal of Physical Chemistry Letters 9 (2018) 4319.
   ! ref : "<Phonons Theory and Experiments I Lattice Dynamics and Models of Interatomic Forces by Dr. Peter Brüesch (auth.) (z-lib.org).pdf>."
   subroutine derivs_nuclei(nmodes,nq,dEa_dQ,wf,ld_gamma,xx,vv,dx,dv)
-    use elph2,only : iminusq_sub
+    use elph2,only : iminusq
     implicit none
     integer,intent(in) :: nmodes,nq
     complex(kind=dpc),intent(in) :: dEa_dQ(nmodes,nq)
@@ -231,7 +230,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_sub
+    use elph2 , only : iminusq
     implicit none
     integer,intent(in) :: nmodes,nq,nfre,nfre_sh
     integer,intent(in) :: isurface
@@ -252,8 +251,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_sub(iq))+&
-                                    (DD(isurface,ifre,imode,iq)*DD(ifre,isurface,imode,iminusq_sub(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))))
           endif
         enddo
       enddo
@@ -271,7 +270,7 @@ end subroutine get_dEa2_dQ2
   SUBROUTINE ADD_BATH_EFFECT(nmodes,nq,wf,ld_gamma,temp,dEa2_dQ2,TT,l_ph_quantum,XX,VV)
     use kinds,only : dp,dpc
     use parameters,only : lit_ephonon,eps_acustic
-    use elph2,only :  iminusq_sub
+    use elph2,only :  iminusq
     use randoms,only : GAUSSIAN_RANDOM_NUMBER_FAST
     use constants,only : KB=>K_B_Ryd,sqrt3,sqrt5,sqrt7,ryd2mev,cone,ci,tpi
     implicit none
@@ -332,16 +331,16 @@ SUBROUTINE ADD_BATH_EFFECT(nmodes,nq,wf,ld_gamma,temp,dEa2_dQ2,TT,l_ph_quantum,X
   ENDSUBROUTINE  
 
 
-  subroutine pre_md(nmodes,nq,wf,ld_gamma,temp,phQ,phP,l_ph_quantum,pre_dt)
+  subroutine pre_md(nmodes,nqtotf,wf,ld_gamma,temp,phQ,phP,l_ph_quantum,pre_dt)
     use surfacecom,only : dEa_dQ,dEa2_dQ2,pre_nstep
     use constants,only  : ry_to_fs,ryd2eV
     use io,only : stdout
     implicit none
-    integer ,intent(in) :: nmodes,nq
+    integer ,intent(in) :: nmodes,nqtotf
     real(kind=dp),intent(in) :: pre_dt,temp
     logical,intent(in) :: l_ph_quantum
-    real(kind=dp),intent(in) :: wf(nmodes,nq),ld_gamma(nmodes,nq)
-    complex(kind=dpc),intent(inout) :: phQ(nmodes,nq),phP(nmodes,nq)
+    real(kind=dp),intent(in) :: wf(nmodes,nqtotf),ld_gamma(nmodes,nqtotf)
+    complex(kind=dpc),intent(inout) :: phQ(nmodes,nqtotf),phP(nmodes,nqtotf)
 
 
     integer :: istep
@@ -351,8 +350,8 @@ subroutine pre_md(nmodes,nq,wf,ld_gamma,temp,phQ,phP,l_ph_quantum,pre_dt)
     dEa_dQ = 0.0
     dEa2_dQ2 = 0.0
     do istep=1,pre_nstep
-      call rk4_nuclei(nmodes,nq,dEa_dQ,ld_gamma,wf,phQ,phP,pre_dt)
-      call add_bath_effect(nmodes,nq,wf,ld_gamma,temp,dEa2_dQ2,pre_dt,l_ph_quantum,phQ,phP)
+      call rk4_nuclei(nmodes,nqtotf,dEa_dQ,ld_gamma,wf,phQ,phP,pre_dt)
+      call add_bath_effect(nmodes,nqtotf,wf,ld_gamma,temp,dEa2_dQ2,pre_dt,l_ph_quantum,phQ,phP)
     enddo
 
     time = pre_nstep*pre_dt*ry_to_fs

src_complex/elph2.f90

@@ -39,7 +39,6 @@ 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
@@ -74,11 +73,9 @@ 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,9 +420,7 @@ module klist_epw
     use kinds, only :dp  
     implicit none
     integer, allocatable :: kqmap(:,:)
-    !! map of k+q grid into k grid
-    integer, allocatable :: kqmap_sub(:,:)
-    !! 
+    !!!! map of k+q grid into k grid
     integer, allocatable :: isk_all(:)
     !! Spin index of each k-point (used in LSDA calculations only)
     INTEGER, ALLOCATABLE :: isk_loc(:)

src_complex/getwcvk.f90

@@ -6,15 +6,14 @@ 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,nk_sub,fwhm,w_center)
+  subroutine get_Wcvk(ihband_min,ieband_max,fwhm,w_center)
     !得到光激发下垂直跃迁的跃迁几率
-    use elph2,only  : vmef,nkf,etf_sub !vmef(3,nbndsub,nbndsub,nkf)
-    use readepw,only : icbm
-    use surfacecom,only : indexk
+    use elph2,only  : vmef,nkf  !vmef(3,nbndsub,nbndsub,nkf)
+    use readepw,only : etf,icbm
     use io,only : stdout
     use constants,only : ryd2eV,ry_to_fs
     implicit none
-    integer , intent(in) :: ihband_min,ieband_max,nk_sub
+    integer , intent(in) :: ihband_min,ieband_max
     real(kind=dp),intent(in) :: fwhm
     real(kind=dp),intent(in) :: w_center
 
@@ -30,29 +29,28 @@ subroutine get_Wcvk(ihband_min,ieband_max,nk_sub,fwhm,w_center)
     integer :: ivbm
 
     ivbm = icbm-1
-
-    allocate(W_cvk(icbm:ieband_max,ihband_min:ivbm,nk_sub),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
-    do ik=1,nk_sub
+    do ik=1,nkf
       do ibnd=ihband_min,ivbm
         do jbnd=icbm,ieband_max
           Evmef = 0.0
-          E_mnk = etf_sub(jbnd,ik)-etf_sub(ibnd,ik)
+          E_mnk = etf(jbnd,ik)-etf(ibnd,ik)
           do ipol=1,3
-            Evmef =Evmef+ (efield_cart(ipol)*(vmef(ipol,jbnd,ibnd,indexk(ik))))
+            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
         enddo
       enddo
     enddo  
 
-
     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
     !% Write laser information            %!
     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!