diff --git a/src_complex/hamiltonian.f90 b/src_complex/hamiltonian.f90 index f376a16..90520be 100644 --- a/src_complex/hamiltonian.f90 +++ b/src_complex/hamiltonian.f90 @@ -3,7 +3,7 @@ module hamiltonian use io, only : msg,io_error,stdout use types use parameters - use elph2, only : nbndfst,epmatq,wf + use elph2, only : nbndfst,epmatq,wf,wf_sub,nqtotf use modes, only : nmodes use readepw,only : E_nk use constants,only : ryd2mev,cone,czero @@ -24,6 +24,10 @@ module hamiltonian complex(kind=dp),allocatable :: gmnvkq_e(:,:,:,:,:) complex(kind=dp),allocatable :: gmnvkq_h(:,:,:,:,:) + type(gmnvkq_n0),allocatable :: gijk_e(:),gijk_h(:) + + integer :: gnzfree_e,gnzfree_h + integer :: ngfre_e,ngfre_h !平衡位置哈密顿量，电声耦合项，总的H和临时H @@ -131,18 +135,21 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihb end subroutine allocate_hamiltonian - subroutine set_H0_nk(nk_sub,indexk,nq_sub,indexq,ibandmin,ibandmax,Enk,H0_nk,gmnvkq_eh) + subroutine set_H0_nk(nk_sub,indexk,nq_sub,indexq,ibandmin,ibandmax,Enk,H0_nk,gmnvkq_eh,lit_g,gnzfree) use elph2, only : etf_sub,gmnvkq + use epwcom,only : kqmap_sub implicit none integer , intent(in) :: 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 , intent(out) :: gnzfree + real(kind=dp) , intent(in) :: lit_g + integer :: nband - integer :: ik,iband - integer :: m,n,nu,iq + integer :: ik,iband,jband + integer :: m,n,nu,iq,ikq nband = ibandmax-ibandmin+1 @@ -160,9 +167,61 @@ subroutine set_H0_nk(nk_sub,indexk,nq_sub,indexq,ibandmin,ibandmax,Enk,H0_nk,gmn enddo enddo + gnzfree = 0 + do iq=1,nq_sub + do ik=1,nk_sub + ikq = kqmap_sub(ik,iq) + if(ikq /= 0) then + do nu=1,nmodes + do iband=1,nband + do jband=1,nband + if(ABS(gmnvkq_eh(jband,iband,nu,ik,iq)) > sqrt(2.0*wf_sub(nu,iq)/nqtotf)*lit_g) gnzfree = gnzfree +1 + enddo + enddo + enddo + endif + enddo + enddo + end subroutine set_H0_nk + subroutine set_gijk(nband,nmodes,nk_sub,nq_sub,gmnvkq,lit_g,gnzfree,gijk) + use epwcom,only : kqmap_sub + implicit none + integer, intent(in) :: nband,nmodes,nk_sub,nq_sub,gnzfree + complex(kind=dpc),intent(in) :: gmnvkq(nband,nband,nmodes,nk_sub,nq_sub) + real(kind=dp) , intent(in) :: lit_g + type(gmnvkq_n0),intent(out) :: gijk(gnzfree) + + integer :: iq,ik,ikq,nu,iband,jband,ink,iqv,imkq + integer :: ig + + ig = 0 + do iq=1,nq_sub + do ik=1,nk_sub + ikq = kqmap_sub(ik,iq) + if(ikq /= 0) then + do nu=1,nmodes + do iband=1,nband + do jband=1,nband + if(ABS(gmnvkq(jband,iband,nu,ik,iq)) > sqrt(2.0*wf_sub(nu,iq)/nqtotf)*lit_g) then + ig = ig +1 + gijk(ig)% iqv = (iq-1)*nmodes + nu + gijk(ig)% ink = (ik-1)*nband + iband + gijk(ig)% imkq= (ikq-1)*nband + jband + gijk(ig)% g = gmnvkq(jband,iband,nu,ik,iq) + endif + enddo + enddo + enddo + endif + enddo + enddo + + + + end subroutine !ref: 1 G. GRIMvall, 1981), ! (3.20) (6.4) diff --git a/src_complex/lvcsh.f90 b/src_complex/lvcsh.f90 index 879b44e..af3c099 100644 --- a/src_complex/lvcsh.f90 +++ b/src_complex/lvcsh.f90 @@ -39,10 +39,12 @@ program lvcsh 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 + gmnvkq_e,gnzfree_e,gijk_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,gnzfree_h,gijk_h,Enk_h,H0_h_nk,E0_h,P0_h, & + allocate_hamiltonian,set_H_nk, set_gijk, & + 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 @@ -74,7 +76,8 @@ program lvcsh dc3_e,dc4_e,& w_h,w0_h,g_h,g1_h,hsurface_type,cc0_h,dc1_h,dc2_h, & dc3_h,dc4_h,allocatesh,convert_diabatic_adiabatic, & - calculate_nonadiabatic_coupling, & + calculate_nonadiabatic_coupling, & + calculate_nonadiabatic_coupling_,& calculate_hopping_probability, & convert_adiabatic_diabatic,add_decoherence use elph2,only : wf_sub,xkf,nqtotf,nktotf,nbndfst,gmnvkq,etf,epmatq @@ -126,16 +129,19 @@ program lvcsh call allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihband_min,ihband_max) if(lelecsh) then - call set_H0_nk(nk_sub,indexk,nq_sub,indexq,ieband_min,ieband_max,Enk_e,H0_e_nk,gmnvkq_e) + call set_H0_nk(nk_sub,indexk,nq_sub,indexq,ieband_min,ieband_max,Enk_e,H0_e_nk,gmnvkq_e,lit_gmnvkq,gnzfree_e) !set Enk_e(nband,nk_sub) H0_e_nk and gmnvkq_e + allocate(gijk_e(gnzfree_e)) + call set_gijk(neband,nmodes,nk_sub,nq_sub,gmnvkq_e,lit_gmnvkq,gnzfree_e,gijk_e) if(nefre_sh == 0 .or. nefre_sh > nefre) nefre_sh = nefre endif if(lholesh) then - call set_H0_nk(nk_sub,indexk,nq_sub,indexq,ihband_min,ihband_max,Enk_h,H0_h_nk,gmnvkq_h) + call set_H0_nk(nk_sub,indexk,nq_sub,indexq,ihband_min,ihband_max,Enk_h,H0_h_nk,gmnvkq_h,lit_gmnvkq,gnzfree_h) H0_h_nk = -1.0 * H0_h_nk gmnvkq_h = -1.0 * gmnvkq_h - !set Enk_h(nband,nk_sub) H0_h_nk and gmnvkq_h(:::,nk_sub,nq_sub) + !set Enk_h(nband,nk_sub) H0_h_nk and gmnvkq_h(:::,nk_sub,nq_sub) + allocate(gijk_h(gnzfree_h)) if(nhfre_sh == 0 .or. nhfre_sh > nhfre) nhfre_sh = nhfre endif @@ -191,7 +197,8 @@ program lvcsh call convert_adiabatic_diabatic(nefre,P_e,w_e,c_e) c_e_nk = reshape(c_e,(/ neband,nk_sub /)) - call calculate_nonadiabatic_coupling(nmodes,nq_sub,neband,nk_sub,E_e,P_e_nk,gmnvkq_e,lit_gmnvkq,nefre_sh,d_e) + call calculate_nonadiabatic_coupling(nmodes,nq_sub,neband,nk_sub,E_e,P_e_nk,gmnvkq_e,nefre_sh,d_e) + !call calculate_nonadiabatic_coupling_(nmodes,nq_sub,neband,nk_sub,E_e,P_e_nk,gnzfree_e,gijk_e,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 @@ -207,7 +214,7 @@ program lvcsh call convert_adiabatic_diabatic(nefre,P_h,w_h,c_h) c_h_nk = reshape(c_h,(/ nhband,nk_sub /)) - call calculate_nonadiabatic_coupling(nmodes,nq_sub,nhband,nk_sub,E_h,P_h_nk,gmnvkq_h,lit_gmnvkq,nhfre_sh,d_h) + call calculate_nonadiabatic_coupling(nmodes,nq_sub,nhband,nk_sub,E_h,P_h_nk,gmnvkq_h,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 @@ -274,7 +281,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,nq_sub,neband,nk_sub,E_e,p_e,gmnvkq_e,lit_gmnvkq,nefre_sh,d_e) + call calculate_nonadiabatic_coupling(nmodes,nq_sub,neband,nk_sub,E_e,p_e,gmnvkq_e,nefre_sh,d_e) ! use p_e in time t0+dt, to convert c_e(t0+dt) to w_e(t0+dt) @@ -334,7 +341,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,nq_sub,nhband,nk_sub,E_h,p_h,gmnvkq_h,lit_gmnvkq,nhfre_sh,d_h) + call calculate_nonadiabatic_coupling(nmodes,nq_sub,nhband,nk_sub,E_h,p_h,gmnvkq_h,nhfre_sh,d_h) ! use p_h in time t0+dt, to convert c_h(t0+dt) to w_h(t0+dt) diff --git a/src_complex/surfacehopping.f90 b/src_complex/surfacehopping.f90 index ba92bd7..bf3c96a 100644 --- a/src_complex/surfacehopping.f90 +++ b/src_complex/surfacehopping.f90 @@ -349,19 +349,18 @@ end subroutine add_decoherence !==================================================================! ! ref : PPT-91 ! The most time-consuming part of the program - subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit_gmnvkq,nfre_sh,dd) + subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfre_sh,dd) use kinds,only : dp implicit none integer, intent(in) :: nmodes,nq,nband,nk,nfre_sh real(kind=dp),intent(in) :: ee(nband*nk) complex(kind=dpc),intent(in) :: p_nk(nband,nk,nband*nk) - real(kind=dp),intent(in) :: lit_gmnvkq complex(kind=dpc),intent(in) :: gmnvkq(nband,nband,nmodes,nk,nq) complex(kind=dpc),intent(out):: dd(nfre_sh,nfre_sh,nmodes,nq) integer :: nfre,ifre,jfre,iq,imode ,igfre integer :: ik,ikq,iband1,iband2 - complex(kind=dpc) :: epc + complex(kind=dpc) :: epc,pnk12 !nfre = nband*nk nfre = nfre_sh @@ -379,7 +378,7 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit do ifre=1,nfre do jfre=1,nfre dd(ifre,jfre,imode,iq) = dd(ifre,jfre,imode,iq)+& - &epc*CONJG(p_nk(iband1,ik,ifre))*p_nk(iband2,ikq,jfre) + epc*CONJG(p_nk(iband1,ik,ifre))*p_nk(iband2,ikq,jfre) enddo enddo endif @@ -389,7 +388,7 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit endif enddo enddo - + do ifre=1,nfre do jfre=1,nfre if(jfre/= ifre) then @@ -402,6 +401,51 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit end subroutine calculate_nonadiabatic_coupling + subroutine calculate_nonadiabatic_coupling_(nmodes,nq,nband,nk,ee,p_nk,gnzfree,gijk,nfre_sh,dd) + use types + use kinds,only : dp + implicit none + integer, intent(in) :: nmodes,nq,nband,nk,nfre_sh,gnzfree + real(kind=dp),intent(in) :: ee(nband*nk) + complex(kind=dpc),intent(in) :: p_nk(nband*nk,nband*nk) + type(gmnvkq_n0),intent(in) :: gijk(gnzfree) + complex(kind=dpc),intent(out):: dd(nfre_sh,nfre_sh,nmodes*nq) + + integer :: nfre,ifre,jfre,iq,imode ,ig + integer :: ik,ikq,iband1,iband2 ,iqv,ink,imkq + complex(kind=dpc) :: epc,pnk12,g + + nfre = nfre_sh + + dd=czero + do ig=1,gnzfree + iqv = gijk(ig) % iqv + ink = gijk(ig) % ink + imkq= gijk(ig) % imkq + epc = gijk(ig) % g + do ifre=1,nfre + do jfre=1,nfre + dd(ifre,jfre,iqv) = dd(ifre,jfre,iqv)+& + epc*CONJG(p_nk(ink,ifre))*p_nk(imkq,jfre) + enddo + enddo + + enddo + + + + do ifre=1,nfre + do jfre=1,nfre + if(jfre/= ifre) then + dd(ifre,jfre,:) = dd(ifre,jfre,:)/(ee(jfre)-ee(ifre)) + else + !dd(ifre,jfre,:,:) = czero + endif + enddo + enddo + + end subroutine calculate_nonadiabatic_coupling_ + function bolziman(womiga,temp) use kinds ,only : dp diff --git a/src_complex/types.f90 b/src_complex/types.f90 index 7c7f0bf..56a0b1d 100644 --- a/src_complex/types.f90 +++ b/src_complex/types.f90 @@ -1,15 +1,12 @@ module types - use kinds,only :dp + use kinds,only :dp,dpc implicit none !declare type of gmnvkq_n0 with gmnvkq>0.0 type :: gmnvkq_n0 - integer :: m - integer :: n - integer :: v - integer :: ik - integer :: iq - integer :: ikq - real(kind=dp) :: g + integer :: iqv + integer :: ink + integer :: imkq + complex(kind=dpc) :: g end type gmnvkq_n0 end module types \ No newline at end of file