change define initial states

README.md

@@ -146,8 +146,8 @@ lehpairsh     = .true.
 !ieband_max    = 4
 !ihband_min    = 1
 !ihband_max    = 2
-lsortpes      = .false.
-mix_thr       = 0.8
+!lsortpes      = .false.
+!mix_thr       = 0.8
 !nefre_sh      = 40
 !nhfre_sh      = 40
 epwoutname    = "./QEfiles/epw40.out"

src_complex/initialsh.f90

@@ -80,7 +80,7 @@ end subroutine set_subband
 
   ! ref: 1 S. Fernandez-Alberti et al., The Journal of Chemical Physics 137 (2012) 
   ! ref: 2. HuangKun <固体物理> (9-29) (9-31)
-  subroutine init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband)
+  subroutine init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband,e_en,h_en)
     use parameters,only : init_ikx,init_iky,init_ikz,init_kx,init_ky,init_kz
     use epwcom,only : nkf1,nkf2,nkf3
     use readepw,only : etf,icbm,xkf
@@ -96,6 +96,7 @@ subroutine init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband)
     logical,intent(in)    :: llaser
     integer,intent(out)   :: init_ik
     integer,intent(inout) :: init_eband,init_hband
+    real(kind=dp),intent(out) :: e_en,h_en
     integer :: ik,ihband,ieband
     real(kind=dp) :: flagr,flagd,W_cvk_all
     real(kind=dp) :: obsorbEn
@@ -121,19 +122,12 @@ subroutine init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband)
                  enddo
                enddo
       enddo outter
-      obsorbEn = (etf(init_eband,init_ik)-etf(init_hband,init_ik))*ryd2eV
-      write(stdout,"(/5X,A)") "Initial eh_KSstate:  "
-      write(stdout,"(5X,A3,I5,1X,A10,3(F12.6,1X),A2)") "ik=",init_ik, "coord.: ( ",(xkf(ipol,init_ik),ipol=1,3)," )"
-      write(stdout,"(5X,A11,I5,A21,F12.5,A3)") "init_hband=",init_hband," Initial hole Energy:",&
-                                                etf(init_hband,init_ik)*ryd2eV," eV"
-      write(stdout,"(5X,A11,I5,A21,F12.5,A3)") "init_eband=",init_eband," Initial elec Energy:",&
-                                                etf(init_eband,init_ik)*ryd2eV," eV"      
-      write(stdout,"(5X,A18,F12.5,A3)")"elec-hole energy= ",obsorbEn," eV"
+
     else
       init_ikx = get_ik(init_kx,nkf1)
       init_iky = get_ik(init_ky,nkf2)
       init_ikz = get_ik(init_kz,nkf3)
-      init_ik  =  (init_ikx - 1) * nkf2 * nkf3 + (init_iky - 1) * nkf3 + init_ikz
+      init_ik  =  (init_ikx - 1) * nkf2 * nkf3 + (init_iky - 1) * nkf3 + init_ikz     
       if(init_eband < icbm .or. init_eband > ieband_max) then
         write(stdout,"(/5X,A29,I5)") "Wrong! The init_eband set as:",init_eband
         write(stdout,"(5X,A29,I5,A16,I5)") "The init_eband need to set : ",icbm,"<= init_eband <=",ieband_max
@@ -143,6 +137,17 @@ subroutine init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband)
         write(stdout,"(5X,A29,I5,A16,I5)") "The init_hband need to set : ",ihband_min,"<= init_hband <=",ivbm
       endif      
     endif
+
+    obsorbEn = (etf(init_eband,init_ik)-etf(init_hband,init_ik))*ryd2eV
+    e_en = etf(init_eband,init_ik)
+    h_en = -1.0*etf(init_hband,init_ik)
+    write(stdout,"(/5X,A)") "Initial eh_KSstate:  "
+    write(stdout,"(5X,A3,I5,1X,A10,3(F12.6,1X),A2)") "ik=",init_ik, "coord.: ( ",(xkf(ipol,init_ik),ipol=1,3)," )"
+    write(stdout,"(5X,A11,I5,A21,F12.5,A3)") "init_hband=",init_hband," Initial hole Energy:",&
+                                              etf(init_hband,init_ik)*ryd2eV," eV"
+    write(stdout,"(5X,A11,I5,A21,F12.5,A3)") "init_eband=",init_eband," Initial elec Energy:",&
+                                              etf(init_eband,init_ik)*ryd2eV," eV"      
+    write(stdout,"(5X,A18,F12.5,A3)")"elec-hole energy= ",obsorbEn," eV"
 
   end subroutine init_eh_KSstat
 
@@ -155,11 +160,35 @@ subroutine init_stat_diabatic(init_ik,init_band,iband_min,nband,nk,c_nk)
     complex(kind=dpc),intent(out) :: c_nk(nband,nk) 
 
     init_band = init_band - iband_min + 1
-    c_nk = 0.0d0
-    c_nk(init_band,init_ik) = 1.0d0
+    c_nk = czero
+    c_nk(init_band,init_ik) = cone
 
   end subroutine init_stat_diabatic
-
+
+  subroutine init_stat_adiabatic(nfre,EE,nfre_sh,init_En,isurface)
+    implicit none
+    integer, intent(in) :: nfre,nfre_sh
+    real(kind=dp),intent(in) :: EE(nfre)
+    real(kind=dp),intent(in) :: init_En
+    integer, intent(out)     :: isurface
+
+    integer :: ifre
+    real(kind=dp),allocatable :: Ptsur(:)
+    integer :: localp(1)
+
+    allocate(Ptsur(1:nfre_sh))
+    Ptsur = 0.0
+    do ifre = 1,nfre_sh
+      Ptsur(ifre) = (EE(ifre)-init_En)**2
+    enddo
+
+    localp = MinLOC(Ptsur)
+    isurface = localp(1)
+
+    deallocate(Ptsur)
+
+  end subroutine
+
   function get_ik(kx,nkx)
     use kinds,only : dp
     implicit none

src_complex/lvcsh.f90

@@ -27,29 +27,28 @@ program lvcsh
   use omp_lib
   use f95_precision
   use blas95  
-  use constants,only      : maxlen,ryd2eV,ry_to_fs,ryd2meV
+  use constants,only      : maxlen,ryd2eV,ry_to_fs,ryd2meV,czero,cone
   use environments,only   : environment_start,mkl_threads,&
                             set_mkl_threads,lsetthreads
   use readinput,only      : get_inputfile
   use readscf,only        : readpwscf_out
   use readphout,only      : readph_out
   use readepw,only        : readepwout
-  use parameters, only    : lreadscfout,scfoutname,lreadphout,phoutname,       &
+  use parameters, only    : lreadscfout,scfoutname,lreadphout,phoutname,nnode, &
                             epwoutname,inputfilename,llaser,init_ik,init_eband,&
-                            init_hband,mix_thr,lsortpes,calculation,verbosity, &
-                            nnode,ncore,naver_sum,savedsnap
+                            init_hband,init_e_en,init_h_en,mix_thr,lsortpes,   &
+                            calculation,verbosity,naver_sum,savedsnap,ncore
   use hamiltonian,only    : nefre,neband,H_e,H_e_nk,E_e,P_e,P_e_nk,P0_e_nk,    &
                             gmnvkq_e,Enk_e,H0_e_nk,E0_e,P0_e,nhfre,nhband,H_h, &
                             H_h_nk,E_h,P_h,P_h_nk,P0_h_nk,gmnvkq_h,Enk_h,      &
                             H0_h_nk,E0_h,P0_h,allocate_hamiltonian,set_H_nk,   &
                             set_H0_nk,calculate_eigen_energy_state    
-
   use sortting,only       : resort_eigen_energy_stat
   use randoms,only        : init_random_seed
   use lasercom,only       : fwhm,w_laser
   use getwcvk,only        : get_Wcvk
   use initialsh,only      : set_subband,init_normalmode_coordinate_velocity,   &
-                            init_eh_KSstat,init_stat_diabatic,init_surface
+                            init_eh_KSstat,init_stat_adiabatic,init_surface
   use surfacecom,only     : methodsh,lfeedback,lit_gmnvkq,naver,nsnap,nstep,dt,&
                             pre_nstep,pre_dt,l_ph_quantum,gamma,temp,iaver,    &
                             isnap,istep,ldecoherence,Cdecoherence,             &
@@ -123,8 +122,8 @@ program lvcsh
 
   if(lholesh) then
     call set_H0_nk(nktotf,ihband_min,ihband_max,Enk_h,H0_h_nk,gmnvkq_h)
-    H0_h_nk = -1.0 * H0_h_nk
-    gmnvkq_h  = -1.0 * gmnvkq_h    
+    H0_h_nk  = -1.0 * H0_h_nk
+    gmnvkq_h = -1.0 * gmnvkq_h    
     !set H0_h_nk and gmnvkq_h  
     if(nhfre_sh == 0 .or. nhfre_sh > nhfre) nhfre_sh = nhfre
   endif
@@ -188,32 +187,37 @@ program lvcsh
 
 
     !!得到初始电子和空穴的初始的KS状态 init_ik,init_eband,init_hband(in the diabatic states)
-    call init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband)
+    call init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband,init_e_en,init_h_en)
+
 
     if(lelecsh) then
-      call init_stat_diabatic(init_ik,init_eband,ieband_min,neband,nktotf,c_e_nk)
-      c_e = reshape(c_e_nk,(/nefre/))
       call set_H_nk(neband,nktotf,nmodes,nqtotf,phQ,gmnvkq_e,H0_e_nk,H_e_nk)
       H_e = reshape(H_e_nk,(/ nefre,nefre /))
       call calculate_eigen_energy_state(nefre,H_e,E_e,P_e)
       P_e_nk = reshape(P_e,(/ neband,nktotf,nefre /))
-      call convert_diabatic_adiabatic(nefre,P_e,c_e,w_e)
-      call init_surface(nefre,nefre_sh,w_e,iesurface)
+
+      call init_stat_adiabatic(nefre,E_e,nefre_sh,init_e_en,iesurface)
+      w_e = czero
+      w_e(iesurface) = cone
+      call convert_adiabatic_diabatic(nefre,P_e,w_e,c_e)
+      c_e_nk = reshape(c_e,(/ neband,nktotf /))
 
       call calculate_nonadiabatic_coupling(nmodes,nqtotf,neband,nktotf,E_e,P_e_nk,gmnvkq_e,lit_gmnvkq,nefre_sh,d_e)
       E0_e = E_e;P0_e=P_e;P0_e_nk=P_e_nk;d0_e=d_e;w0_e=w_e
     endif
 
     if(lholesh) then
-      call init_stat_diabatic(init_ik,init_hband,ihband_min,nhband,nktotf,c_h_nk)
-      c_h = reshape(c_h_nk,(/nhfre/))
       call set_H_nk(nhband,nktotf,nmodes,nqtotf,phQ,gmnvkq_h,H0_h_nk,H_h_nk)
       H_h = reshape(H_h_nk,(/ nhfre,nhfre /))    
       call calculate_eigen_energy_state(nhfre,H_h,E_h,P_h)
       P_h_nk = reshape(P_h,(/ nhband,nktotf,nhfre /))
-      call convert_diabatic_adiabatic(nhfre,P_h,c_h,w_h)
-      call init_surface(nhfre,nhfre_sh,w_h,ihsurface)
-
+
+      call init_stat_adiabatic(nhfre,E_h,nhfre_sh,init_h_en,ihsurface)
+      w_h = czero
+      w_h(ihsurface) = cone
+      call convert_adiabatic_diabatic(nefre,P_h,w_h,c_h)      
+      c_h_nk = reshape(c_h,(/ nhband,nktotf /))
+
       call calculate_nonadiabatic_coupling(nmodes,nqtotf,nhband,nktotf,E_h,P_h_nk,gmnvkq_h,lit_gmnvkq,nhfre_sh,d_h)
       E0_h = E_h;P0_h=P_h;P0_h_nk=P_h_nk;d0_h=d_h;w0_h=w_h
     endif

src_complex/parameters.f90

@@ -18,6 +18,7 @@ module parameters
 
   real(kind=dp)   :: init_kx,init_ky,init_kz  !激发后初始的电子(空穴)的k坐标
   integer         :: init_hband,init_eband    !激发后初始的电子(空穴)所处的能带
+  real(kind=dp)   :: init_e_en,init_h_en      !激发后初始电子(空穴)的能量
   integer         :: init_ikx,init_iky,init_ikz,init_ik 
   real(kind=dp)   :: mix_thr
 	! threshold to find the mixxing states

src_complex/surfacehopping.f90

@@ -11,8 +11,8 @@ module surfacehopping
                          phQ,phP,phQ0,phP0,phK,phU,SUM_phU,SUM_phK,SUM_phK0,SUM_phE,&
                          phQsit,phPsit,phKsit,phUsit,ld_gamma,&
                          dEa_dQ,dEa_dQ_e,dEa_dQ_h,dEa2_dQ2,dEa2_dQ2_e,dEa2_dQ2_h,&
-                         d_e,g_e,g1_e,c_e_nk,w_e,w0_e,d0_e,nefre_sh,&
-                         d_h,g_h,g1_h,c_h_nk,w_h,w0_h,d0_h,nhfre_sh,&
+                         d_e,g_e,g1_e,c_e,c_e_nk,w_e,w0_e,d0_e,nefre_sh,&
+                         d_h,g_h,g1_h,c_h,c_h_nk,w_h,w0_h,d0_h,nhfre_sh,&
                          pes_one_e,pes_e,csit_e,wsit_e,psit_e,&
                          pes_one_h,pes_h,csit_h,wsit_h,psit_h 
   use cc_fssh,only : S_ai_e,S_ai_h,S_bi_e,S_bi_h
@@ -98,6 +98,8 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)
 
       allocate(c_e_nk(neband,nktotf),stat=ierr,errmsg=msg)
       if(ierr /=0) call errore('surfacehopping','Error allocating c_e_nk',1)
+      allocate(c_e(neband*nktotf),stat=ierr,errmsg=msg)
+      if(ierr /=0) call errore('surfacehopping','Error allocating c_e',1)      
 
       allocate(pes_one_e(0:nefre,0:nsnap),pes_e(0:nefre,0:nsnap))
 			ram = real_size*(1+nefre)*(1+nsnap)
@@ -173,6 +175,8 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)
 
       allocate(c_h_nk(nhband,nktotf),stat=ierr,errmsg=msg)
       if(ierr /=0) call errore('surfacehopping','Error allocating c_h_nk',1)  
+      allocate(c_h(nhband*nktotf),stat=ierr,errmsg=msg)
+      if(ierr /=0) call errore('surfacehopping','Error allocating c_h',1)  
 
       allocate(pes_one_h(0:nhfre,0:nsnap),pes_h(0:nhfre,0:nsnap))
 			ram = real_size*(1+nhfre)*(1+nsnap)