quantum-espresso/GWW/simple_bse/hamiltonian.f90

165 lines
5.2 KiB
Fortran

SUBROUTINE hamiltonian(sin,n,bd,pp,pt,pm,a,ha,ilevel)
!this subroutine applies the excitonic hamiltonian
USE simple_objects
USE derived_objects
USE input_simple_exc
USE mp, ONLY : mp_sum,mp_barrier
USE mp_world, ONLY : world_comm
USE io_global, ONLY : stdout
implicit none
TYPE(input_options), INTENT(in) :: sin
INTEGER :: n!number of vectors
TYPE(bands), INTENT(in) :: bd
TYPE(prod_proj), INTENT(in) :: pp
TYPE(potential), INTENT(in) :: pt
TYPE(prod_mix), INTENT(in) :: pm
TYPE(exc), INTENT(in) :: a(n)
TYPE(exc), INTENT(inout) :: ha(n)
INTEGER :: ilevel!for analysis:0, All, 1 Diagonal, 2 Exchange ,3 Direct
INTEGER :: i,ik,jk,ii,jj,kk
COMPLEX(kind=DP), ALLOCATABLE :: emat(:,:),vemat(:,:)
COMPLEX(kind=DP), ALLOCATABLE :: hmat(:,:,:), imat(:,:,:)
COMPLEX(kind=DP), ALLOCATABLE :: bvc(:,:),bvc_t(:,:)
COMPLEX(kind=DP), ALLOCATABLE :: gcc_t(:,:,:),imat_t(:,:,:), gvv_t(:,:,:)
COMPLEX(kind=DP) :: fact
INTEGER :: ia,ikk,iv,ic,icp,ivp
!Diagonal part
do i=1,n
ha(i)= bd .hd. a(i)
end do
if(ilevel>1) ha(1)%avc(1:ha(1)%numv,1:ha(1)%numc,1:ha(1)%nk_loc)=(0.d0,0.d0)
if(sin%h_level >= 1 .and. sin%spin_state>=1 .and. (ilevel/=1)) then
!Exchange term
allocate(emat(pp%nprod_e,n),vemat(pp%nprod_e,n))
emat=(0.d0,0.d0)
if(sin%spin_state==1) then
fact=(2.d0,0.d0)
else
fact=(1.d0,0.d0)
endif
do i=1,n
if(pp%nk_loc>0) then
call ZGEMM('N','N',pp%nprod_e,1,pp%numc*pp%numv*pp%nk_loc,(1.d0,0.d0),pp%javc,pp%nprod_e,&
&a(i)%avc,pp%numc*pp%numv*pp%nk_loc,(0.d0,0.d0),emat(1,i),pp%nprod_e)
else
emat(1:pp%nprod_e,i)=(0.d0,0.d0)
endif
end do
call mp_sum(emat,world_comm)
call ZGEMM('N','N',pp%nprod_e,n,pp%nprod_e,fact,pt%vpot,pt%nprod_e,emat,pp%nprod_e,(0.d0,0.d0),&
&vemat, pp%nprod_e)
do i=1,n
if(pp%nk_loc>0) then
call ZGEMM('C','N',pp%numc*pp%numv*pp%nk_loc,1,pp%nprod_e,(1.d0,0.d0),pp%javc,pp%nprod_e,vemat(1,i),pp%nprod_e,&
(1.d0,0.d0),ha(i)%avc,pp%numc*pp%numv*pp%nk_loc)
endif
enddo
deallocate(emat,vemat)
endif
if(ilevel==3) ha(1)%avc=0.d0
if(sin%h_level >= 2 .and. ( ilevel==0 .or. ilevel==3)) then
!TD-HF term
allocate(hmat(pm%nprod_e,pm%numv,pm%numv))
allocate(imat(pm%nprod_e,pm%numv,pm%numc))
allocate(bvc(pm%numv,pm%numc))
allocate(imat_t(pm%nprod_e,pm%numc, pm%numv))
allocate(gcc_t(pm%nprod_e,pm%numc, pm%numc))
allocate(bvc_t(pm%numc,pm%numv))
allocate(gvv_t(pm%nprod_e,pm%numv,pm%numv))
!write(stdout,*) 'DEBUG',pm%nk,pm%nk_loc
do i=1,n
!loop on k
do ik=1,pm%nk
!loop on k'
bvc=(0.d0,0.d0)
do jk=1,pm%nk_loc
!multiply V*Gvv'
ii=pt%ijk(1,ik,jk+pm%ik_first-1)+1
jj=pt%ijk(2,ik,jk+pm%ik_first-1)+1
kk=pt%ijk(3,ik,jk+pm%ik_first-1)+1
!change indices
do ia=1,pm%nprod_e
do iv=1,pm%numv
do ivp=1,pm%numv
! gvv_t(ia,iv,ivp)=pm%gvv(ia,iv,ik,ivp,jk)
hmat(ia,iv,ivp)=pm%gvv(ia,iv,ik,ivp,jk)
enddo
enddo
enddo
! hmat=1.d0!DEBUG
! call ZGEMM('N','N',pm%nprod_e,pm%numv*pm%numv, pm%nprod_e,&
! &(-1.d0,0.d0),pt%vpotq(1,1,ii,jj,kk),pm%nprod_e,gvv_t, pm%nprod_e,(0.d0,0.d0),hmat,pm%nprod_e)
!in case add W_c term
! if(sin%h_level >= 3) then
! call ZGEMM('N','N',pm%nprod_e,pm%numv*pm%numv, pm%nprod_e,&
! &(-1.d0,0.d0),pt%wpotq(1,1,ii,jj,kk),pm%nprod_e,gvv_t, pm%nprod_e,(1.d0,0.d0),hmat,pm%nprod_e)
! endif
!multiply *Akv
!multiply Gcc'*
call ZGEMM('N','N',pm%nprod_e*pm%numv,pm%numc,pm%numv,(1.d0,0.d0)&
&,hmat,pm%nprod_e*pm%numv,a(i)%avc(1,1,jk),pm%numv,(0.d0,0.d0),imat,pm%nprod_e*pm%numv)
!now invert indices
do ia=1,pm%nprod_e
do icp=1,pm%numc
do iv=1,pm%numv
imat_t(ia,icp,iv)=imat(ia,iv,icp)
enddo
enddo
enddo
do ia=1,pm%nprod_e
do icp=1,pm%numc
do ic=1,pm%numc
gcc_t(ia,icp,ic)=pm%gcc(ia,ic,ik,icp,jk)
enddo
enddo
enddo
call ZGEMM('C','N',pm%numc,pm%numv,pm%nprod_e*pm%numc,(1.d0,0.d0),&
&gcc_t,pm%nprod_e*pm%numc,imat_t,pm%nprod_e*pm%numc,(0.d0,0.d0),&
&bvc_t,pm%numc)
!put indices in correct order
do iv=1,pm%numv
do ic=1,pm%numc
bvc(iv,ic)=bvc(iv,ic)+bvc_t(ic,iv)
enddo
enddo
enddo!on jk
call mp_sum(bvc,world_comm)
!DEBUG START
!bvc=bvc/dble(pm%nk)
!DEBUG END
if(ik>=pm%ik_first .and. ik<=pm%ik_last) then
ha(i)%avc(1:pm%numv,1:pm%numc,ik-pm%ik_first+1)=ha(i)%avc(1:pm%numv,1:pm%numc,ik-pm%ik_first+1)+&
&bvc(1:pm%numv,1:pm%numc)
endif
enddo!ok ik
enddo!on i
deallocate(hmat,imat,bvc)
deallocate(gcc_t,imat_t,bvc_t,gvv_t)
endif
return
END SUBROUTINE hamiltonian