change dEa2_dQ2

src_complex/cc_fssh.f90

@@ -87,7 +87,7 @@ subroutine nonadiabatic_transition_ccfssh(lfeedback,nfre,nfre_sh,nmodes,nq,&
     integer , intent(out)    :: surface_type
     real(kind=dp),intent(in) :: EE(nfre)
     complex(kind=dpc),intent(in) :: P(nfre,nfre),P0(nfre,nfre)
-    complex(kind=dpc),intent(in) :: DD(nfre_sh,nfre_sh,nmodes,nq)
+    complex(kind=dpc),intent(in) :: DD(2,nfre_sh,nmodes,nq)
     real(kind=dp),intent(out):: S_bi(nfre_sh)
     real(kind=dp),intent(in) :: GG(nfre_sh)
     complex(kind=dpc),intent(inout) :: VV(nmodes,nq)
@@ -138,8 +138,8 @@ subroutine nonadiabatic_transition_ccfssh(lfeedback,nfre,nfre_sh,nmodes,nq,&
           sumdd = 0.0
           do iq=1,nq
             do imode=1,nmodes
-              sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(isurface,ifre,imode,iq)) ! A
-              sumdd = sumdd + ABS(DD(isurface,ifre,imode,iq))**2  ! B      
+              sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(1,ifre,imode,iq)) ! A
+              sumdd = sumdd + ABS(DD(1,ifre,imode,iq))**2  ! B      
             enddo
           enddo
           detaE = EE(isurface_a)-EE(ifre) ! C
@@ -152,7 +152,7 @@ subroutine nonadiabatic_transition_ccfssh(lfeedback,nfre,nfre_sh,nmodes,nq,&
               flagd = sumvd/sumdd*(-1.0+dsqrt(flagd))
               do iq=1,nq
                 do imode=1,nmodes
-                  if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(isurface,ifre,imode,iq)
+                  if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(1,ifre,imode,iq)
                 enddo
               enddo
               isurface = isurface_k
@@ -176,7 +176,7 @@ subroutine nonadiabatic_transition_ccfssh(lfeedback,nfre,nfre_sh,nmodes,nq,&
                 flagd = sumvd/sumdd*(-1.0+dsqrt(flagd))
                 do iq=1,nq
                   do imode=1,nmodes
-                    if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(isurface,ifre,imode,iq)
+                    if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(1,ifre,imode,iq)
                   enddo
                 enddo
                 isurface = isurface_k

src_complex/dynamics.f90

@@ -27,40 +27,6 @@ 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
-    implicit none
-    integer,intent(in) :: nmodes,nq,nband,nk
-    integer,intent(in) :: isurface
-    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)
-
-    integer :: iq,imode,ik,ikq,iband1,iband2
-    complex(kind=dpc) :: epc
-    complex(kind=dpc) :: cmp_tmp1,cmp_tmp2
-    logical :: ltmp = .true.
-
-    ! F_qv = -dEa_dQ_-qv  (2.60)
-    dEa_dQ = czero
-    do iq=1,nq   
-      do imode=1,nmodes
-        do ik=1,nk
-          ikq = kqmap(ik,iq)
-          do iband1=1,nband ! m
-            do iband2=1,nband  ! n
-              epc = gmnvkq(iband1,iband2,imode,ik,iq)
-              if(epc /= czero) then
-                dEa_dQ(imode,iminusq(iq)) = dEa_dQ(imode,iminusq(iq)) + &
-                epc*CONJG(P_nk(iband1,ikq,isurface))*P_nk(iband2,ik,isurface)                      
-              endif
-            enddo
-          enddo
-        enddo
-      enddo
-    enddo
-
-  end subroutine get_dEa_dQ
 
   !=======================================================================!
   != rk4 method to obtain coordinate and velocitie after a time interval =!
@@ -135,14 +101,6 @@ subroutine derivs_nuclei(nmodes,nq,dEa_dQ,wf,ld_gamma,xx,vv,dx,dv)
     dv = -wf**2*xx - dEa_dQ - ld_gamma*vv
     !(2.60) PPT-94
 
-    ! F_qv = -dEa_dQ_-qv
-    !do iq=1,nq
-    !  do imode=1,nmodes
-    !    dx(imode,iq) = vv(imode,iq)
-    !    dv(imode,iq) = -wf(imode,iq)**2*xx(imode,iq)-dEa_dQ(imode,iq)-ld_gamma(imode,iq)*vv(imode,iq)
-    !  enddo
-    !enddo
-
   endsubroutine derivs_nuclei
 
 
@@ -219,20 +177,19 @@ subroutine get_dEa2_dQ2(nmodes,nq,nfre,nfre_sh,isurface,EE,dd,dEa2_dQ2)
     integer,intent(in) :: nmodes,nq,nfre,nfre_sh
     integer,intent(in) :: isurface
     real(kind=dp),intent(in)  :: EE(nfre)
-    complex(kind=dpc),intent(in)  :: dd(nfre_sh,nmodes,nq)
+    complex(kind=dpc),intent(in)  :: dd(2,nfre_sh,nmodes,nq)
     real(kind=dp),intent(out) :: dEa2_dQ2(nmodes,nq)
 
     integer :: iq,imode,ifre
 
-    ! dEa2_dQ2 = dEa2/d(Q_qv*)d(Q_qv)
     ! ref : (2.60) PPT-95
     dEa2_dQ2 = 0.0
     do iq=1,nq
       do imode=1,nmodes
         do ifre=1,nfre_sh
           if(ifre /= isurface) then         
             dEa2_dQ2(imode,iq) = dEa2_dQ2(imode,iq) + &
-            (EE(isurface)-EE(ifre))*(DD(ifre,imode,iq)*CONJG(DD(ifre,imode,iq)))*2.0
+            (EE(isurface)-EE(ifre))*REAL((DD(1,ifre,imode,iq)*DD(2,ifre,imode,iq)))*2.0
           endif
         enddo
       enddo

src_complex/fssh.f90

@@ -13,7 +13,7 @@ subroutine nonadiabatic_transition_fssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurfac
     integer , intent(inout)  :: isurface
     real(kind=dp),intent(in) :: EE(nfre)
     complex(kind=dpc),intent(in) :: P(nfre,nfre)
-    complex(kind=dpc),intent(in) :: DD(nfre_sh,nmodes,nq)
+    complex(kind=dpc),intent(in) :: DD(2,nfre_sh,nmodes,nq)
     real(kind=dp),intent(in) :: GG(nfre_sh)
     complex(kind=dp),intent(inout) :: VV(nmodes,nq)
 
@@ -37,8 +37,8 @@ subroutine nonadiabatic_transition_fssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurfac
           sumdd = 0.0
           do iq=1,nq
             do imode=1,nmodes
-              sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(ifre,imode,iq)) ! A
-              sumdd = sumdd + ABS(DD(ifre,imode,iq))**2  ! B                   
+              sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(1,ifre,imode,iq)) ! A
+              sumdd = sumdd + ABS(DD(1,ifre,imode,iq))**2  ! B                   
             enddo
           enddo
           detaE = EE(isurface)-EE(ifre) !C
@@ -51,7 +51,7 @@ subroutine nonadiabatic_transition_fssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurfac
             flagd = (sumvd/sumdd)*(-1.0+dsqrt(flagd))  ! effective scratting time
             do iq=1,nq
               do imode=1,nmodes
-                if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(ifre,imode,iq)
+                if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(1,ifre,imode,iq)
               enddo
             enddo
             isurface = ifre          

src_complex/lvcsh.f90

@@ -80,7 +80,7 @@ program lvcsh
   use cell_base,only      : bg
   use date_and_times,only : get_date_and_time
   use io      ,only       : stdout,io_time,time1,time2,time_,msg,io_error
-  use dynamics,only       : set_gamma,get_dEa_dQ,get_dEa2_dQ2,rk4_nuclei,&
+  use dynamics,only       : set_gamma,get_dEa2_dQ2,rk4_nuclei,&
                             rk4_electron_diabatic,ADD_BATH_EFFECT,pre_md
   use memory_report,only  : MB,GB,complex_size,real_size,int_size,ram,&
                             print_memory  
@@ -183,6 +183,7 @@ program lvcsh
         c_e_nk = reshape(c_e,(/ neband,nktotf /))
 
         call calculate_nonadiabatic_coupling(nmodes,nqtotf,neband,nktotf,E_e,P_e_nk,gmnvkq_e,nefre_sh,iesurface,d_e)
+        dEa_dQ_e = d_e(1,iesurface,:,:)
         E0_e = E_e;P0_e=P_e;P0_e_nk=P_e_nk;d0_e=d_e;w0_e=w_e
       endif
 
@@ -199,6 +200,7 @@ program lvcsh
         c_h_nk = reshape(c_h,(/ nhband,nktotf /))
 
         call calculate_nonadiabatic_coupling(nmodes,nqtotf,nhband,nktotf,E_h,P_h_nk,gmnvkq_h,nhfre_sh,ihsurface,d_h)
+        dEa_dQ_h = d_h(1,ihsurface,:,:)
         E0_h = E_h;P0_h=P_h;P0_h_nk=P_h_nk;d0_h=d_h;w0_h=w_h
       endif
 
@@ -222,22 +224,15 @@ program lvcsh
           !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
           dEa_dQ = 0.0
           !dEa_dQ in time t0
-          if(l_dEa_dQ) then
-            if(lelecsh) then
-              call get_dEa_dQ(nmodes,nqtotf,neband,nktotf,P0_e_nk,gmnvkq_e,iesurface,dEa_dQ_e)
-              dEa_dQ = dEa_dQ + dEa_dQ_e
-            endif
-            if(lholesh) then
-              call get_dEa_dQ(nmodes,nqtotf,nhband,nktotf,P0_h_nk,gmnvkq_h,ihsurface,dEa_dQ_h)
-              dEa_dQ = dEa_dQ + dEa_dQ_h
-            endif
-          endif
+          if(lelecsh) dEa_dQ = dEa_dQ + dEa_dQ_e
+          if(lholesh) dEa_dQ = dEa_dQ + dEa_dQ_h
 
           !==============================!
           != update phQ,phP             =!
           !==============================!        
           !use rk4 to calculate the dynamical of phonon normal modes
           !update phQ,phP to time t0+dt
+          if(.not. l_dEa_dQ) dEa_dQ = 0.0
           call rk4_nuclei(nmodes,nqtotf,dEa_dQ,ld_gamma,wf,phQ,phP,dt)
 
 
@@ -266,7 +261,7 @@ program lvcsh
             ! Calculate non-adiabatic coupling vectors with the Hellmann-Feynman theorem.
             ! update d_e in time t0+dt
             call calculate_nonadiabatic_coupling(nmodes,nqtotf,neband,nktotf,E_e,p_e,gmnvkq_e,nefre_sh,iesurface,d_e)
-
+            dEa_dQ_e = d_e(1,iesurface,:,:)
 
             ! use p_e in time t0+dt, to convert c_e(t0+dt) to w_e(t0+dt) 
             call convert_diabatic_adiabatic(nefre,p_e,c_e,w_e)
@@ -326,7 +321,7 @@ program lvcsh
             ! Calculate non-adiabatic coupling vectors with the Hellmann-Feynman theorem.
             ! update d_h in time t0+dt
             call calculate_nonadiabatic_coupling(nmodes,nqtotf,nhband,nktotf,E_h,p_h,gmnvkq_h,nhfre_sh,ihsurface,d_h)          
-
+            dEa_dQ_h = d_h(1,ihsurface,:,:)
 
             ! use p_h in time t0+dt, to convert c_h(t0+dt) to w_h(t0+dt) 
             call convert_diabatic_adiabatic(nhfre,p_h,c_h,w_h)
@@ -368,20 +363,19 @@ program lvcsh
           !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
           dEa2_dQ2 = 0.0
           !dEa2_dQ2 in time t0
-          if(l_dEa2_dQ2) then
-            if(lelecsh) then
-              call get_dEa2_dQ2(nmodes,nqtotf,nefre,nefre_sh,iesurface,E0_e,d0_e,dEa2_dQ2_e)
-              dEa2_dQ2 = dEa2_dQ2 + dEa2_dQ2_e
-            endif
-            if(lholesh) then
-              call get_dEa2_dQ2(nmodes,nqtotf,nhfre,nhfre_sh,ihsurface,E0_h,d0_h,dEa2_dQ2_h)
-              dEa2_dQ2 = dEa2_dQ2 + dEa2_dQ2_h
-            endif
+          if(lelecsh) then
+            call get_dEa2_dQ2(nmodes,nqtotf,nefre,nefre_sh,iesurface,E0_e,d0_e,dEa2_dQ2_e)
+            dEa2_dQ2 = dEa2_dQ2 + dEa2_dQ2_e
+          endif
+          if(lholesh) then
+            call get_dEa2_dQ2(nmodes,nqtotf,nhfre,nhfre_sh,ihsurface,E0_h,d0_h,dEa2_dQ2_h)
+            dEa2_dQ2 = dEa2_dQ2 + dEa2_dQ2_h
           endif
 
           !===================!
           != add bath effect =!
-          !===================!  
+          !===================!
+          if(.not. l_dEa2_dQ2) dEa2_dQ2 = 0.0
           call add_bath_effect(nmodes,nqtotf,wf,ld_gamma,temp,dEa2_dQ2,dt,l_ph_quantum,phQ,phP)
 
 

src_complex/sc_fssh.f90

@@ -70,7 +70,7 @@ subroutine nonadiabatic_transition_scfssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurf
     integer , intent(inout)  :: isurface
     real(kind=dp),intent(in) :: EE(nfre)
     complex(kind=dpc),intent(in) :: P(nfre,nfre)
-    complex(kind=dpc),intent(in) :: DD(nfre_sh,nfre_sh,nmodes,nq)
+    complex(kind=dpc),intent(in) :: DD(2,nfre_sh,nmodes,nq)
     real(kind=dp),intent(in) :: GG(nfre_sh)
     complex(kind=dpc),intent(inout) :: VV(nmodes,nq)
 
@@ -94,8 +94,8 @@ subroutine nonadiabatic_transition_scfssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurf
           sumdd = 0.0
           do iq=1,nq
             do imode=1,nmodes
-              sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(isurface,ifre,imode,iq)) ! A
-              sumdd = sumdd + ABS(DD(isurface,ifre,imode,iq))**2  ! B       
+              sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(1,ifre,imode,iq)) ! A
+              sumdd = sumdd + ABS(DD(1,ifre,imode,iq))**2  ! B       
             enddo
           enddo
           detaE = EE(isurface)-EE(ifre) ! C
@@ -107,7 +107,7 @@ subroutine nonadiabatic_transition_scfssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurf
             flagd = sumvd/sumdd*(-1.0+dsqrt(flagd))
             do iq=1,nq
               do imode=1,nmodes
-                if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(isurface,ifre,imode,iq)
+                if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(1,ifre,imode,iq)
               enddo
             enddo
             isurface = ifre          

src_complex/surfacecom.f90

@@ -79,7 +79,7 @@ module surfacecom
   real(kind=dp) :: E_ph_CA_sum,E_ph_QA_sum
   ! averager crystal energy and temperature T,in classical and quantum . 
 
-  complex(kind=dpc),allocatable :: d_e(:,:,:),d_h(:,:,:),d0_e(:,:,:),d0_h(:,:,:)
+  complex(kind=dpc),allocatable :: d_e(:,:,:,:),d_h(:,:,:,:),d0_e(:,:,:,:),d0_h(:,:,:,:)
   real(kind=dp),allocatable :: g_e(:) ,g_h(:)
   real(kind=dp),allocatable :: g1_e(:),g1_h(:)  
 

src_complex/surfacehopping.f90

@@ -91,9 +91,9 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)
       if(ierr /=0) call errore('surfacehopping','Error allocating P_e_nk',1)
       allocate(P0_e_nk(neband,nktotf,nefre),stat=ierr,errmsg=msg)
       if(ierr /=0) call errore('surfacehopping','Error allocating P0_e_nk',1)			
-      allocate(d_e(nefre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ajk
+      allocate(d_e(2,nefre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ajk,d_jak
       if(ierr /=0) call errore('surfacehopping','Error allocating d_e',1)
-			ram=real_size*nmodes*nq*nefre_sh
+			ram=real_size*nmodes*nq*nefre_sh*2
 			call print_memory("d_e",ram)
       allocate(dEa_dQ_e(nmodes,nq))
       allocate(dEa2_dQ2_e(nmodes,nq))
@@ -149,7 +149,7 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)
       if(ierr /=0) call errore('surfacehopping','Error allocating E0_e',1)
       allocate(P0_e(nefre,nefre),stat=ierr,errmsg=msg)
       if(ierr /=0) call errore('surfacehopping','Error allocating P0_e',1)
-      allocate(d0_e(nefre_sh,nmodes,nq),stat=ierr,errmsg=msg)
+      allocate(d0_e(2,nefre_sh,nmodes,nq),stat=ierr,errmsg=msg)
       if(ierr /=0) call errore('surfacehopping','Error allocating d0_e',1)
 
       if(methodsh == "CC-FSSH") then
@@ -168,7 +168,7 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)
       if(ierr /=0) call errore('surfacehopping','Error allocating P_h_nk',1)
       allocate(P0_h_nk(nhband,nktotf,nhfre),stat=ierr,errmsg=msg)
       if(ierr /=0) call errore('surfacehopping','Error allocating P0_h_nk',1)			
-      allocate(d_h(nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ajk
+      allocate(d_h(2,nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ajk
       if(ierr /=0) call errore('surfacehopping','Error allocating d_h',1)
 			ram=real_size*nmodes*nq*nhfre_sh
 			call print_memory("d_h",ram)
@@ -225,7 +225,7 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)
 
       allocate(P0_h(nhfre,nhfre),stat=ierr,errmsg=msg)
       if(ierr /=0) call errore('surfacehopping','Error allocating P0_h',1)
-      allocate(d0_h(nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg)
+      allocate(d0_h(2,nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg)
       if(ierr /=0) call errore('surfacehopping','Error allocating d0_h',1)
 
       if(methodsh == "CC-FSSH") then
@@ -356,7 +356,7 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr
     real(kind=dp),intent(in)     :: ee(nband*nk)
     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):: dd(nfre_sh,nmodes,nq)
+    complex(kind=dpc),intent(out):: dd(2,nfre_sh,nmodes,nq)
 
 		integer :: nfre,ifre,iq,imode
     integer :: ik,ikq,m,n
@@ -374,8 +374,11 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr
             do n=1,nband ! n
               epc = gmnvkq(m,n,imode,ik,iq)
               if(epc /= czero) then
-                do ifre=1,nfre                   
-                  dd(ifre,imode,iq) = dd(ifre,imode,iq)+epc*CONJG(p_nk(m,ikq,isurface))*p_nk(n,ik,ifre)
+                do ifre=1,nfre
+                  ! d_ajk
+                  dd(1,ifre,imode,iq) = dd(1,ifre,imode,iq)+epc*CONJG(p_nk(m,ikq,isurface))*p_nk(n,ik,ifre)
+                  ! d_jak
+                  dd(2,ifre,imode,iq) = dd(2,ifre,imode,iq)+epc*CONJG(p_nk(m,ikq,ifre))*p_nk(n,ik,isurface)
                 enddo
               endif
             enddo
@@ -386,7 +389,8 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr
 
 		do ifre=1,nfre
       if(ifre/= isurface) then
-        dd(ifre,:,:) = dd(ifre,:,:)/(ee(ifre)-ee(isurface))
+        dd(1,ifre,:,:) = dd(1,ifre,:,:)/(ee(ifre)-ee(isurface))
+        dd(2,ifre,:,:) = dd(2,ifre,:,:)/(ee(isurface)-ee(ifre))
       endif
 		enddo
 
@@ -416,7 +420,7 @@ subroutine calculate_hopping_probability(isurface,nfre,nfre_sh,nmodes,nq,WW,VV,d
     integer,intent(in)           :: isurface,nfre,nfre_sh,nmodes,nq
     complex(kind=dpc),intent(in) :: WW(nfre)
     complex(kind=dpc),intent(in) :: VV(nmodes,nq)
-    complex(kind=dpc),intent(in) :: dd(nfre_sh,nmodes,nq)
+    complex(kind=dpc),intent(in) :: dd(2,nfre_sh,nmodes,nq)
     real(kind=dp),intent(in)     :: tt
     real(kind=dp),intent(out)    :: gg(nfre_sh)
     real(kind=dp),intent(out)    :: gg1(nfre_sh)
@@ -433,7 +437,7 @@ subroutine calculate_hopping_probability(isurface,nfre,nfre_sh,nmodes,nq,WW,VV,d
         sumvd = czero
         do iq=1,nq
           do imode=1,nmodes
-            sumvd = sumvd+VV(imode,iq)*dd(ifre,imode,iq)
+            sumvd = sumvd+VV(imode,iq)*dd(1,ifre,imode,iq)
           enddo
         enddo
         ! in adiabatic representation：the switching probabilities from the active surface isurface to another surface iefre