scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Simple_bse code 

P.Umari G.Prandini N.Marzari
This commit is contained in:
paoloumari 2018-06-30 10:30:41 +02:00
parent 49778475df
commit 5c33e6c21f
13 changed files with 2506 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
GWW/simple_bse/build_eemat.f90 Normal file
View File

@ -0,0 +1,76 @@
subroutine build_eemat(data_input,bd,eemat)
USE input_simple_exc
USE simple_objects
USE derived_objects
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode
USE constants, ONLY : rytoev
USE mp, ONLY : mp_barrier, mp_sum
USE mp_world, ONLY : world_comm,mpime
use io_files, only : prefix, tmp_dir
USE constants, ONLY : rytoev
implicit none
TYPE(input_options) :: data_input
TYPE(epe) :: element
TYPE(bands) :: bd
COMPLEX(KIND=DP),ALLOCATABLE::Etempc(:,:,:),Etempv(:,:,:),mate(:,:,:,:)
COMPLEX(KIND=DP) :: eemat(bd%numv,bd%numc,bd%nk_loc,3)
INTEGER :: a, i, j, v, c, k
allocate(Etempv(bd%ntot_e,bd%numv,bd%nk_loc))!vettori E_{iv}^k
allocate(Etempc(bd%ntot_e,bd%numc,bd%nk_loc))!vettori E_{ic}^k
allocate(mate(bd%ntot_e,bd%numc,bd%nk_loc,3))!E_{jc}^k <e_i|p_a|e_j>
call initialize_epe(element)
write(stdout,*)'reading epe'
call read_epe(data_input,element)
write(stdout ,*)'EPE MATRIX DIMENSION', element%ntot_e
!costruisco E_{jc}^k
do k=1,bd%nk_loc
do c=bd%numv+1,bd%num
do j=1,bd%ntot_e
Etempc(j,c-bd%numv,k) = bd%omat(j,c,k)
end do
end do
end do
write(stdout,*) 'E_{jc}^k'
!prodotto <e_i|p_a|e_j>E_{jc}^k
do a=1,3
do k=1,bd%nk_loc
call zgemm('N','N',bd%ntot_e,bd%numc,bd%ntot_e,(1.d0,0.d0),element%mat(:,:,a),bd%ntot_e,Etempc(:,:,k), &
&bd%ntot_e,(0.d0,0.d0),mate(:,:,k,a),bd%ntot_e)
end do
end do
write(stdout,*)'E_{jc}^k <e_i|p_a|e_j>'
!costruisco E_{iv}^k
do k=1,bd%nk_loc
do v=1,bd%numv
do i=1,bd%ntot_e
Etempv(i,v,k) = bd%omat(i,v,k)
end do
end do
end do
write(stdout,*) 'E_{iv}^k'
!prodotto E_{iv}^k^* E_{jc}^k <e_i|p_a|e_j>
do a=1,3
do k=1,bd%nk_loc
call zgemm('C','N',bd%numv,bd%numc,bd%ntot_e,(1.d0,0.d0),Etempv(1,1,k),bd%ntot_e,mate(1,1,k,a),bd%ntot_e,(0.d0,0.d0),&
&eemat(1,1,k,a),bd%numv)
end do
end do
write(stdout,*) ' E_{iv}^k^* E_{jc}^k <e_i|p_a|e_j>'
deallocate(Etempv)
deallocate(Etempc)
deallocate(mate)
call deallocate_epe(element)
return
end subroutine build_eemat

330
GWW/simple_bse/derived_object.f90 Normal file
View File

@ -0,0 +1,330 @@
!this module contains objected derived and contracted
MODULE derived_objects
USE kinds, ONLY : DP
TYPE prod_proj
!terms <\mathcal{E}_\alpha|u_{k,c}u_{k,v}^*>
!k ponits distributed on MPI tasks as bands object
INTEGER :: numv !number of valence states (those considered for excitons only)
INTEGER :: numc !number of conduction states
INTEGER :: nk!total number of k, points
INTEGER :: nk_loc!local number of k points
INTEGER :: ik_first!first local k point
INTEGER :: ik_last!last local k point
INTEGER :: ntot_e!dimension of global to all k, basis for KS states
INTEGER :: nprod_e!number of product terms
COMPLEX(kind=DP), DIMENSION(:,:,:,:), POINTER :: javc! (nprod_e,numv,numc,nk_loc)
END TYPE prod_proj
TYPE prod_mix
!terms <\mathcal{E}_\alpha|(u_{k,v}u_{k',v'}^*)>
!and terms <\mathcal{E}_\alpha|(u_{k,c}u_{k',c'}^*)>
!k' distributed over MPI tasks
!k NOT distributed
INTEGER :: numv !number of valence states (those considered for excitons only)
INTEGER :: numc !number of conduction states
INTEGER :: nk!total number of k, points
INTEGER :: nk_loc!local number of k points
INTEGER :: ik_first!first local k point
INTEGER :: ik_last!last local k point
INTEGER :: ntot_e!dimension of global to all k, basis for KS states
INTEGER :: nprod_e!number of product terms
COMPLEX(kind=DP), DIMENSION(:,:,:,:,:), POINTER :: gvv! (nprod_e,numv,nk,numv',nk_loc)! ' means relative to nk_loc
COMPLEX(kind=DP), DIMENSION(:,:,:,:,:), POINTER :: gcc! (nprod_e,numc,nk,numc',nk_loc)
END TYPE prod_mix
CONTAINS
SUBROUTINE initialize_prod_proj(pp)
implicit none
TYPE(prod_proj) :: pp
nullify(pp%javc)
return
END SUBROUTINE initialize_prod_proj
SUBROUTINE deallocate_prod_proj(pp)
implicit none
TYPE(prod_proj) :: pp
if(associated(pp%javc)) deallocate(pp%javc)
nullify(pp%javc)
return
END SUBROUTINE deallocate_prod_proj
SUBROUTINE initialize_prod_mix(pm)
implicit none
TYPE(prod_mix) :: pm
nullify(pm%gvv)
nullify(pm%gcc)
return
END SUBROUTINE initialize_prod_mix
SUBROUTINE deallocate_prod_mix(pm)
implicit none
TYPE(prod_mix) :: pm
if(associated(pm%gvv)) deallocate(pm%gvv)
nullify(pm%gvv)
if(associated(pm%gcc)) deallocate(pm%gcc)
nullify(pm%gcc)
return
END SUBROUTINE deallocate_prod_mix
SUBROUTINE build_prod_proj(bd,pd,pp)
!this subroutine constructs the prod_proj object
USE simple_objects, ONLY : bands,product
implicit none
TYPE(bands), INTENT(in) :: bd
TYPE(product), INTENT(in) :: pd
TYPE(prod_proj), INTENT(out) :: pp
INTEGER :: ik,iv,ic
COMPLEX(kind=DP), ALLOCATABLE :: tmp_mat(:,:,:),tmp_mat2(:,:),tmp_mat3(:,:)
COMPLEX(kind=DP), ALLOCATABLE :: zmat(:,:,:,:),zmat0(:,:),zmat1(:,:,:,:)
LOGICAL, parameter :: debug = .false.
pp%numv=bd%numv
pp%numc=bd%numc
pp%ntot_e=bd%ntot_e
pp%nk=bd%nk
pp%nk_loc=bd%nk_loc
pp%ik_first=bd%ik_first
pp%ik_last=bd%ik_last
pp%nprod_e=pd%nprod_e
if(pp%nk_loc>0) then
allocate(pp%javc(pp%nprod_e,pp%numv, pp%numc,pp%nk_loc))
allocate(tmp_mat(pp%nprod_e,pp%ntot_e,pp%numc))
allocate(tmp_mat2(pp%ntot_e,pp%numv))
allocate(tmp_mat3(pp%ntot_e,pp%numc))
do ik=1,pp%nk_loc
tmp_mat2(1:pp%ntot_e,1:pp%numv)=conjg(bd%omat(1:pp%ntot_e,1:pp%numv,ik))
tmp_mat3(1:pp%ntot_e,1:pp%numc)=bd%omat(1:pp%ntot_e,pp%numv+1:pp%numv+pp%numc,ik)
call ZGEMM('N','N',pp%nprod_e*pp%ntot_e,pp%numc,pp%ntot_e,(1.d0,0.d0),pd%fij,&
&pp%nprod_e*pp%ntot_e,tmp_mat3,pp%ntot_e,(0.d0,0.d0),tmp_mat,pp%nprod_e*pp%ntot_e)
do ic=1,pp%numc
call ZGEMM('N','N',pp%nprod_e,pp%numv,pp%ntot_e,(1.d0,0.d0),tmp_mat(1,1,ic),pp%nprod_e,tmp_mat2,pp%ntot_e,&
&(0.d0,0.d0),pp%javc(1,1,ic,ik),pp%nprod_e)
enddo
enddo
deallocate(tmp_mat,tmp_mat2,tmp_mat3)
else
nullify(pp%javc)
endif
if(debug) then
!test for consistency
allocate(zmat(pp%numc,pp%numv,pp%numv,pp%numc))
allocate(zmat0(bd%num,bd%num))
allocate(zmat1(pd%ntot_e,pd%ntot_e,pd%ntot_e,pd%ntot_e))
do ik=1,bd%nk_loc
call ZGEMM('C','N',bd%num,bd%num,bd%ntot_e,(1.d0,0.d0),bd%omat(1,1,ik),bd%ntot_e,bd%omat(1,1,ik),bd%ntot_e,&
&(0.d0,0.d0),zmat0,bd%num)
do iv=1,bd%num
do ic=1,bd%num
write(*,*) 'CHECK OMAT ik, ic,iv', ik, ic, iv, zmat0(ic,iv)
enddo
enddo
enddo
call ZGEMM('C','N',pd%ntot_e*pd%ntot_e,pd%ntot_e*pd%ntot_e,pd%nprod_e,(1.d0,0.d0),pd%fij,&
pd%nprod_e,pd%fij,pd%nprod_e,(0.d0,0.d0),zmat1,pd%ntot_e*pd%ntot_e)
do iv=1,pd%ntot_e
do ic=1,pd%ntot_e
write(*,*) 'CHECK FIJ ic,iv', ic, iv, zmat1(iv,ic,iv,ic)
enddo
enddo
do ik=1,pp%nk_loc
call ZGEMM('C','N',pp%numc*pp%numv,pp%numc*pp%numv,pp%nprod_e,(1.d0,0.d0),pp%javc(1,1,1,ik),pp%nprod_e,&
& pp%javc(1,1,1,ik),pp%nprod_e,(0.d0,0.0),zmat, pp%numc*pp%numv)
do iv=1,pp%numv
do ic=1,pp%numc
write(*,*) 'CHECK JAVC ik, ic,iv', ik, ic, iv, zmat(iv,ic,iv,ic)
enddo
enddo
enddo
deallocate(zmat,zmat0,zmat1)
endif
return
END SUBROUTINE build_prod_proj
SUBROUTINE build_prod_mix(sin,bd,pd,pm,pt)
!this subroutine constructs the prod_mix object
!COMPLEX(kind=DP), DIMENSION(:,:,:,:,:), POINTER :: gvv! (nprod_e,numv,nk,numv,nk_loc)
!COMPLEX(kind=DP), DIMENSION(:,:,:,:,:), POINTER :: gcc! (nprod_e,numc,nk,numc,nk_loc)
USE input_simple_exc
USE simple_objects, ONLY : bands,product,potential
USE mp_world, ONLY : mpime, world_comm
USE mp, ONLY : mp_sum, mp_bcast
USE io_global, ONLY : stdout
implicit none
TYPE(input_options) :: sin
TYPE(bands), INTENT(in) :: bd
TYPE(product), INTENT(in) :: pd
TYPE(prod_mix), INTENT(out) :: pm
TYPE(potential) :: pt
INTEGER :: ik,iv,ic
COMPLEX(kind=DP), ALLOCATABLE :: tmp_mat(:,:,:),tmp_mat2(:,:),tmp_mat3(:,:,:)
LOGICAL, parameter :: debug = .true.
COMPLEX(kind=DP), ALLOCATABLE :: emat(:,:)
INTEGER :: is_mine
INTEGER :: jk
COMPLEX(kind=DP), ALLOCATABLE :: tmp_pot(:,:),tmp_fij(:,:,:)
INTEGER :: ii,jj,kk
pm%numv=bd%numv
pm%numc=bd%numc
pm%ntot_e=bd%ntot_e
pm%nk=bd%nk
pm%nk_loc=bd%nk_loc
pm%ik_first=bd%ik_first
pm%ik_last=bd%ik_last
pm%nprod_e=pd%nprod_e
if(pm%nk_loc>0) then
allocate( pm%gvv(pm%nprod_e,pm%numv,pm%nk,pm%numv,pm%nk_loc))
allocate( pm%gcc(pm%nprod_e,pm%numc,pm%nk,pm%numc,pm%nk_loc))
else
nullify(pm%gcc)
nullify(pm%gvv)
endif
!now do gvv
!loop on nk
!if ik is my distribute to others
!loop on nk_loc
!calculate terms
allocate(emat(pm%ntot_e,pm%numv))
allocate(tmp_mat2(pm%ntot_e,pm%numv))
allocate(tmp_mat(pm%nprod_e,pm%ntot_e,pm%numv))
allocate(tmp_mat3(pm%nprod_e,pm%numv,pm%numv))
allocate(tmp_pot(pm%nprod_e,pm%nprod_e))
allocate(tmp_fij(pm%nprod_e,pm%ntot_e,pm%ntot_e))
do ik=1,pm%nk
if(ik>=pm%ik_first.and.ik<=pm%ik_last) then
emat(1:pm%ntot_e,1:pm%numv)=bd%omat(1:pm%ntot_e,1:pm%numv,ik-pm%ik_first+1)
is_mine=mpime+1
else
is_mine=0
endif
call mp_sum(is_mine,world_comm)
is_mine=is_mine-1
call mp_bcast( emat,is_mine, world_comm )
do jk=1,pm%nk_loc!on k' local
!find out q=k_i-k_j
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
tmp_pot(1:pm%nprod_e,1:pm%nprod_e)= pt%vpotq(1:pm%nprod_e,1:pm%nprod_e,ii,jj,kk)
!tmp_pot(1:pm%nprod_e,1:pm%nprod_e)=1.d0!DEBUG
if(sin%h_level >= 3) then
tmp_pot(1:pm%nprod_e,1:pm%nprod_e)= tmp_pot(1:pm%nprod_e,1:pm%nprod_e) +&
&pt%wpotq(1:pm%nprod_e,1:pm%nprod_e,ii,jj,kk)
endif
call ZGEMM('N','N',pm%nprod_e,pm%ntot_e*pm%ntot_e,pm%nprod_e,(-1.d0,0.d0),&
&tmp_pot,pm%nprod_e,pd%fij,pm%nprod_e,(0.d0,0.d0),tmp_fij,pm%nprod_e)
tmp_mat2(1:pm%ntot_e,1:pm%numv)=conjg(bd%omat(1:pm%ntot_e,1:pm%numv,jk))
!call ZGEMM('N','N',pm%nprod_e*pm%ntot_e,pm%numv,pm%ntot_e,(1.d0,0.d0),pd%fij,&
! &pm%nprod_e*pm%ntot_e,emat,pm%ntot_e,(0.d0,0.d0),tmp_mat,pm%nprod_e*pm%ntot_e)
call ZGEMM('N','N',pm%nprod_e*pm%ntot_e,pm%numv,pm%ntot_e,(1.d0,0.d0),tmp_fij,&
&pm%nprod_e*pm%ntot_e,emat,pm%ntot_e,(0.d0,0.d0),tmp_mat,pm%nprod_e*pm%ntot_e)
do iv=1,pm%numv
call ZGEMM('N','N',pm%nprod_e,pm%numv,pm%ntot_e,(1.d0,0.d0),tmp_mat(1,1,iv),pm%nprod_e,tmp_mat2,pm%ntot_e,&
&(0.d0,0.d0),tmp_mat3(1,1,iv),pm%nprod_e)
enddo
do iv=1,pm%numv
pm%gvv(1:pm%nprod_e,1:pm%numv,ik,iv,jk)=tmp_mat3(1:pm%nprod_e,iv,1:pm%numv)
enddo
enddo
enddo
deallocate(emat)
deallocate(tmp_mat2)
deallocate(tmp_mat)
deallocate(tmp_mat3)
deallocate(tmp_pot)
deallocate(tmp_fij)
allocate(emat(pm%ntot_e,pm%numc))
allocate(tmp_mat2(pm%ntot_e,pm%numc))
allocate(tmp_mat(pm%nprod_e,pm%ntot_e,pm%numc))
allocate(tmp_mat3(pm%nprod_e,pm%numc,pm%numc))
do ik=1,pm%nk
if(ik>=pm%ik_first.and.ik<=pm%ik_last) then
emat(1:pm%ntot_e,1:pm%numc)=bd%omat(1:pm%ntot_e,pm%numv+1:pm%numv+pm%numc,ik-pm%ik_first+1)
is_mine=mpime+1
else
is_mine=0
endif
call mp_sum(is_mine,world_comm)
is_mine=is_mine-1
call mp_bcast( emat,is_mine, world_comm )
do jk=1,pm%nk_loc!on k' local
tmp_mat2(1:pm%ntot_e,1:pm%numc)=conjg(bd%omat(1:pm%ntot_e,pm%numv+1:pm%numv+pm%numc,jk))
call ZGEMM('N','N',pm%nprod_e*pm%ntot_e,pm%numc,pm%ntot_e,(1.d0,0.d0),pd%fij,&
&pm%nprod_e*pm%ntot_e,emat,pm%ntot_e,(0.d0,0.d0),tmp_mat,pm%nprod_e*pm%ntot_e)
do ic=1,pm%numc
call ZGEMM('N','N',pm%nprod_e,pm%numc,pm%ntot_e,(1.d0,0.d0),tmp_mat(1,1,ic),pm%nprod_e,tmp_mat2,pm%ntot_e,&
&(0.d0,0.d0),tmp_mat3(1,1,ic),pm%nprod_e)
enddo
do ic=1,pm%numc
pm%gcc(1:pm%nprod_e,1:pm%numc,ik,ic,jk)=tmp_mat3(1:pm%nprod_e,ic,1:pm%numc)
enddo
enddo
enddo
deallocate(emat)
deallocate(tmp_mat2)
deallocate(tmp_mat)
deallocate(tmp_mat3)
END SUBROUTINE build_prod_mix
END MODULE derived_objects

319
GWW/simple_bse/diago_exc_cg.f90 Normal file
View File

@ -0,0 +1,319 @@
subroutine diago_exc_cg(data_input, bd, pp, pt, pm, x)
USE input_simple_exc
USE simple_objects
USE derived_objects
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode
USE constants, ONLY : rytoev
USE mp, ONLY : mp_barrier
USE mp_world, ONLY : world_comm
implicit none
TYPE(input_options), INTENT(in) :: data_input
TYPE(bands), INTENT(in) :: bd
TYPE(prod_proj), INTENT(in) :: pp
TYPE(potential), INTENT(in) :: pt
TYPE(prod_mix), INTENT(in) :: pm
TYPE(exc) :: xd, x2d, Hxd, Hx2d, h, Hx, Hx1, Hh, xtest
TYPE(exc) :: x(data_input%nvec)
INTEGER :: i, j, k, count=1, beginning, rate, end,kf
COMPLEX(kind=DP) :: d=(0.1d0,0.0d0), a, b, c, fd, f2d, alpha,stim, minus=(-1.0d0,0.0d0), echeck
COMPLEX(kind=DP), allocatable :: energy(:)
LOGICAL :: test
TYPE(exc) :: etmp1,etmp2,etmp3,etmp4
REAL(kind=DP) :: ene0,dene0,ene1,ene2,passop,passo,stima,enew,gamma,esse,essenew,ene
COMPLEX(kind=DP) :: passot=(1.d0,0.d0),passos,sca, passot2=(0.4d0,0.d0)
LOGICAL :: l_ok
REAL(kind=DP) :: norm
INTEGER :: icon
call setup_exc(bd,h)
call setup_exc(bd,Hx)
call setup_exc(bd,Hx1)
call setup_exc(bd,Hh)
call setup_exc(bd,xtest)
call setup_exc(bd,etmp1)
call setup_exc(bd,etmp2)
call setup_exc(bd,etmp3)
call setup_exc(bd,etmp4)
call mp_barrier(world_comm)
allocate(energy(data_input%nvec))
write(stdout,*)data_input%nvec
vectors: do j=1,data_input%nvec
icon=0
call randomize_exc(x(j))
call normalize_exc(x(j))
!x(j)%avc=0.d0
!x(j)%avc(1,1,1)=(1.d0,0.d0)
call hamiltonian(data_input, 1, bd, pp, pt, pm, x(j), Hx,0)
call mp_barrier(world_comm)
energy(j)=x(j)*Hx*rytoev!energy initialization
h=(-1.d0,0.d0)*Hx!search direction initialization
!call normalize_exc(h)
esse=dble(h*h)
! h=minus*h
l_ok=.true.
iterate:do i=1,data_input%max_nstep
echeck=energy(j)
call mp_barrier(world_comm)
call hamiltonian(data_input, 1, bd, pp, pt, pm, h, Hh,0)
ene0=dble(x(j)*Hx)
dene0=2.d0*dble(x(j)*Hh)-2.d0*ene0*dble(x(j)*h)
norm=sqrt(dble(h*h))
if(.true.) then
if(dene0>0.d0) then
!dene0=-dene0
passot=-(1.d0,0.d0)
passot2=-(2.d0,0.d0)
else
passot=(1.d0,0.d0)
passot2=(2.,0.d0)
endif
else
if(dene0>0.d0) then
!dene0=-dene0
passot=-(0.5d0,0.d0)/norm
passot2=-(1.d0,0.d0)/norm
else
passot=(0.5d0,0.d0)/norm
passot2=(1.,0.d0)/norm
endif
endif
etmp1=passot*h
call sum_exc_sub(etmp2,x(j),etmp1)
call normalize_exc(etmp2)
call hamiltonian(data_input, 1, bd, pp, pt, pm, etmp2, etmp3,0)
ene1=dble(etmp2*etmp3)
etmp1=passot2*h
call sum_exc_sub(etmp2,x(j),etmp1)
call normalize_exc(etmp2)
call hamiltonian(data_input, 1, bd, pp, pt, pm, etmp2, etmp3,0)
ene2=dble(etmp2*etmp3)
!call minparabola(ene0,dene0,ene1,dble(passot),passo,stima,l_ok)
call minparabola2(ene0,ene1,ene2,dble(passot),dble(passot2),passo,stima,l_ok)
! do k=-1000,1000
! passos=cmplx(dble(k)/20.d0,0.d0)
! etmp1=passos*h
! call sum_exc_sub(etmp2,x(j),etmp1)
! call normalize_exc(etmp2)
! call hamiltonian(data_input, 1, bd, pp, pt, pm, etmp2, etmp3,0)
! ene=dble(etmp2*etmp3)
! write(stdout,*) dble(k)/20.d0, ene*rytoev
! enddo
if(l_ok) then
passos=cmplx(passo,0.d0)
etmp1=passos*h
call sum_exc_sub(xtest,x(j),etmp1)
else
!write(stdout,*) l_ok,norm
h=minus*Hx
!ene2=ene0
!kf=0
!do k=1,40
! sca=cmplx(dble(k)/10.d0,0.d0)
! etmp1=sca*h
! call sum_exc_sub(xtest,x(j),etmp1)
! call normalize_exc(xtest)
! call hamiltonian(data_input, 1, bd, pp, pt, pm, xtest, Hx1,0)
! enew=dble(xtest*Hx1)
!
! if(enew<ene2) then
! ene2=enew
! kf=k
! endif
!enddo
!write(stdout,*) 'KF:', kf
!sca=cmplx(dble(kf)/10.d0,0.d0)
sca=(2.d0,0.d0)
etmp1=sca*h
call sum_exc_sub(xtest,x(j),etmp1)
endif
call normalize_exc(xtest)
call hamiltonian(data_input, 1, bd, pp, pt, pm, xtest, Hx1,0)
call mp_barrier(world_comm)
! write(stdout,*)i,'check:',real(echeck),real(xtest*Hx1*rytoev)
enew=dble(xtest*Hx1)
! write(stdout,*) 'LIN MIN', dble(echeck),stima*rytoev,enew*rytoev,passo
write(stdout,*) i,enew*rytoev
if (dble(echeck)-enew*rytoev<0.d0 .or. mod(i,20)==0 .or. .not.l_ok) then
!write(stdout,*) 'NO GOOD'
h=minus*Hx
sca=cmplx(2.d0,0.d0)
etmp1=sca*h
call sum_exc_sub(xtest,x(j),etmp1)
x(j) = xtest
call normalize_exc(x(j))
call mp_barrier(world_comm)
call hamiltonian(data_input, 1, bd, pp, pt, pm, x(j), Hx,0)
call mp_barrier(world_comm)
h=Hx
esse=dble(h*h)
! write(stdout,*)j,i,'sd',real(echeck),real(x(j)*Hx*rytoev)
!write(stdout,*)j,i,real(x(j)*Hx*rytoev)
energy(j)=x(j)*Hx*rytoev
else
x(j)=xtest
Hx=Hx1
gamma=dble(Hx*Hx)
essenew=gamma
gamma=gamma/esse
esse=essenew
sca=cmplx(gamma,0.d0)
etmp1=sca*h
etmp2=(1.d0,0.d0)*Hx
call sum_exc_sub(h,etmp2,etmp1)
!write(stdout,*)j,i,'cg', real(echeck),real(x(j)*Hx*rytoev)
! write(stdout,*)j,i,real(x(j)*Hx*rytoev)
energy(j)=x(j)*Hx*rytoev
if (i==data_input%max_nstep .or. (real(-energy(j)+echeck)<data_input%thr_evc .and. i/=1)) then
write(stdout,*)j,'converged'
exit
end if
end if
if(abs(-energy(j)+echeck)<data_input%thr_evc) then
icon=icon+1
else
icon=0
endif
if(icon==10) exit
end do iterate
x(j)%ene(1)=dble(energy(j))
call hamiltonian(data_input, 1, bd, pp, pt, pm, x(j), Hx,1)
x(j)%ene(2)=dble(x(j)*Hx*rytoev)
call hamiltonian(data_input, 1, bd, pp, pt, pm, x(j), Hx,2)
x(j)%ene(3)=dble(x(j)*Hx*rytoev)
call hamiltonian(data_input, 1, bd, pp, pt, pm, x(j), Hx,3)
x(j)%ene(4)=dble(x(j)*Hx*rytoev)
energy(j)=energy(j)/rytoev!Ry
end do vectors
write(stdout,*) 'calling routine spectrum'
!spectrum
!call spectrum(data_input,x,bd,energy)
!lanczos
deallocate(energy)
call deallocate_exc(Hx)
call deallocate_exc(h)
call deallocate_exc(Hx1)
call deallocate_exc(Hh)
call deallocate_exc(xtest)
call deallocate_exc(etmp1)
call deallocate_exc(etmp2)
call deallocate_exc(etmp3)
call deallocate_exc(etmp4)
return
end subroutine diago_exc_cg
subroutine minparabola(ene0,dene0,ene1,passop,passo,stima,l_ok)
!this subroutines finds the minimum of a quadratic real function \
use kinds, only : dp
use io_global, only :stdout
implicit none
real(dp) ene0,dene0,ene1,passop,passo,stima
real(dp) a,b,c!a*x^2+b*x+c
logical :: l_ok
c=ene0
b=dene0
a=(ene1-b*passop-c)/(passop**2.d0)
passo = -b/(2.d0*a)
if( a.lt.0.d0) then
l_ok=.false.
! write(stdout,*) 'CAZZI'
if(ene1.lt.ene0) then
passo=passop
else
passo=0.5d0*passop
endif
else
l_ok=.true.
endif
stima=a*passo**2.d0+b*passo+c
return
end subroutine minparabola
subroutine minparabola2(ene0,ene1,ene2,x1,x2,x,stima,l_ok)
use kinds, only : DP
use io_global, only : stdout
implicit none
real(kind=dp) :: ene0,ene1,ene2,x1,x2,x,stima
logical :: l_ok
real(kind=dp) :: a,b,c
c=ene0
a=(ene2-ene1*x2/x1+c*x2/x1-c)/(x2**2.d0-x1*x2)
b=(ene1-c-a*x1**2.d0)/x1
x = -b/(2.d0*a)
if( a.lt.0.d0) then
! write(stdout,*) 'CAZZI'
l_ok=.false.
if(ene1<=ene0) then
x=x1
else
x=0.5d0*x1
endif
else
l_ok=.true.
endif
stima=a*x**2.d0+b*x+c
return
end subroutine minparabola2

83
GWW/simple_bse/diago_exc_sd.f90 Normal file
View File

@ -0,0 +1,83 @@
SUBROUTINE diago_exc_sd(sin,bd,pp,pt,pm,a)
!this subroutine find the lowest (AT THE MOMENT)
!eigen values /vector of the excitonic Hamiltonia
!with the simple steepest descent scheme
USE input_simple_exc
USE simple_objects
USE derived_objects
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode
USE constants, ONLY : rytoev
USE mp, ONLY : mp_barrier
USE mp_world, ONLY : world_comm
implicit none
TYPE(input_options), INTENT(in) :: sin
TYPE(bands), INTENT(in) :: bd
TYPE(prod_proj), INTENT(in) :: pp
TYPE(potential), INTENT(in) :: pt
TYPE(prod_mix), INTENT(in) :: pm
TYPE(exc) :: a(sin%nvec)
INTEGER :: ii,jj,kk,is
TYPE(exc) :: ha,b
COMPLEX(kind=DP) ::passop,ene
call setup_exc(bd,ha)
call setup_exc(bd,b)
if(ionode) then
write(stdout,*) 'Routine diago_exc_sd'
endif
!now just one vector
!randoimize vectors
call mp_barrier(world_comm)
write(stdout,*) 'ATT-1'
do ii=1,sin%nvec
call randomize_exc(a(ii))
end do
!orthonormalize them
call mp_barrier(world_comm)
write(stdout,*) 'ATT-2'
do ii=1,sin%nvec
call normalize_exc(a(ii))
enddo
passop=cmplx(sin%l_step,0)
!lopp
do is=1,sin%max_nstep
!!find gradient
call mp_barrier(world_comm)
write(stdout,*) 'ATT6'
call hamiltonian(sin,1,bd,pp,pt,pm,a(1),ha,0)
call mp_barrier(world_comm)
write(stdout,*) 'ATT7'
ene=a(1)*ha
if(ionode) then
write(stdout,*) 'SD step energy :',is, ene*rytoev
endif
write(stdout,*) 'ATT1'
b=(-passop)*ha
write(stdout,*) 'ATT2'
ha=a(1)+b
write(stdout,*) 'ATT3'
call normalize_exc(ha)
write(stdout,*) 'ATT4'
a(1)=ha
write(stdout,*) 'ATT5'
!!find gradient
!!find new vector
!!normalize
!!calculate energy
!!print it
end do
call deallocate_exc(ha)
call deallocate_exc(b)
return
END SUBROUTINE diago_exc_sd

164
GWW/simple_bse/hamiltonian.f90 Normal file
View File

@ -0,0 +1,164 @@
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

73
GWW/simple_bse/input_simple_exc.f90 Normal file
View File

@ -0,0 +1,73 @@
MODULE input_simple_exc
!this module provides input file routines
USE kinds, ONLY: DP
TYPE input_options
CHARACTER(len=256) :: prefix = 'prefix'!prefix to designate the files same as in PW
CHARACTER(len=256) :: outdir = './'!outdir to designate the files same as in PW
INTEGER :: task !0 find eigenvectors/value , 1 calculate spectrum
INTEGER :: diago!for task==0, diagonalization scheme: 0=steepest descent, 1=conjugate gradient
INTEGER :: nvec!for task=0 number of eigen vectors to be found
REAL(kind=DP) :: thr_evc!threshold for the eigenvectors/value
INTEGER :: max_nstep!maximum number of steps during minimization
REAL(kind=DP) :: l_step!length of trial steps for minimization
INTEGER :: h_level!level of excitonic hamiltonian: 0 RPA (diagonal), 1 TD-H (+exchange), 2 TD-HF(+direct-bare), 3-BSE (+direct-correlation)
INTEGER :: spin_state!only for non spin polarised systems: 0:triplet, 1:singlet, 2:full (spin-orbit)
INTEGER :: spectrum_points!number of points in the spectrum
REAL(KIND=DP) :: omega_min!max omega in the spectrum
REAL(KIND=DP) :: omega_max!min omega in the spectrum
REAL(KIND=DP) :: eta!infinitesimal
INTEGER :: lanczos_step!!number of steps for haydoch recursive method
REAL(kind=DP) :: scissor!Gap scissor (eV)
END TYPE input_options
CONTAINS
SUBROUTINE read_input_simple_exc( simple_in )
USE io_global, ONLY : stdout, ionode, ionode_id
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE io_files, ONLY : tmp_dir,prefix
implicit none
CHARACTER(LEN=256), EXTERNAL :: trimcheck
TYPE(input_options) :: simple_in !in output the input parameters
NAMELIST/inputsimple/simple_in
CHARACTER(LEN=256) :: outdir
if(ionode) then
read(*,NML=inputsimple)
outdir = trimcheck(simple_in%outdir)
tmp_dir = outdir
prefix = trim(simple_in%prefix)
endif
CALL mp_bcast( outdir,ionode_id, world_comm )
CALL mp_bcast( tmp_dir,ionode_id, world_comm )
CALL mp_bcast( prefix,ionode_id, world_comm )
CALL mp_bcast( simple_in%task,ionode_id, world_comm )
CALL mp_bcast( simple_in%diago,ionode_id, world_comm )
CALL mp_bcast( simple_in%nvec,ionode_id, world_comm )
CALL mp_bcast( simple_in%thr_evc,ionode_id, world_comm )
CALL mp_bcast( simple_in%max_nstep,ionode_id, world_comm )
CALL mp_bcast( simple_in%l_step,ionode_id, world_comm )
CALL mp_bcast( simple_in%h_level,ionode_id, world_comm )
CALL mp_bcast( simple_in%spin_state,ionode_id, world_comm )
CALL mp_bcast( simple_in%spectrum_points,ionode_id, world_comm )
CALL mp_bcast( simple_in%omega_max,ionode_id, world_comm )
CALL mp_bcast( simple_in%omega_min,ionode_id, world_comm )
CALL mp_bcast( simple_in%eta,ionode_id, world_comm )
CALL mp_bcast( simple_in%lanczos_step,ionode_id, world_comm )
CALL mp_bcast( simple_in%scissor,ionode_id, world_comm )
return
END SUBROUTINE read_input_simple_exc
END MODULE input_simple_exc

186
GWW/simple_bse/lanczos.f90 Normal file
View File

@ -0,0 +1,186 @@
subroutine lanczos(data_input)
USE input_simple_exc
USE simple_objects
USE derived_objects
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode
USE constants, ONLY : rytoev
USE mp, ONLY : mp_barrier, mp_sum
USE mp_world, ONLY : world_comm
use io_files, only : prefix, tmp_dir
USE constants, ONLY : rytoev
implicit none
TYPE(input_options):: data_input
TYPE(exc):: x,Hx,xm1,x1,etmp1,etmp2,etmp3,etmp4
TYPE(epe):: element
TYPE(bands):: bd
TYPE(product) :: pd
TYPE(prod_proj) :: pp
TYPE(potential) :: pt
TYPE(prod_mix) :: pm
COMPLEX(KIND=DP) :: fraz,zeta
REAL(KIND=DP) :: omin,omax,step, cost, omega
COMPLEX(KIND=DP),ALLOCATABLE:: cabs(:,:),eemat(:,:,:,:),a(:),b(:)
INTEGER :: i,dir,v,c,k,io,iun,inizio,fine,count_rate
INTEGER, EXTERNAL :: find_free_unit
REAL(KIND=DP),ALLOCATABLE::diel(:,:)
COMPLEX(kind=DP) :: csca
REAL(kind=DP) :: norm(3)
INTEGER :: idir
call initialize_product(pd)
call initialize_potential(pt)
call initialize_prod_mix(pm)
call initialize_prod_proj(pp)
call system_clock(inizio,count_rate)
call read_bands(data_input,bd)
call read_product(data_input,pd)
call read_potential(data_input,pt)
write(stdout,*)'building derived objects mix'
call build_prod_mix(data_input,bd,pd,pm,pt)
write(stdout,*)'building derived objects proj'
call build_prod_proj(bd,pd,pp)
write(stdout,*)'setup exc'
call setup_exc(bd,x)
call setup_exc(bd,xm1)
call setup_exc(bd,Hx)
call setup_exc(bd,x1)
call setup_exc(bd,etmp1)
call setup_exc(bd,etmp2)
call setup_exc(bd,etmp3)
call setup_exc(bd,etmp4)
allocate(eemat(bd%numv,bd%numc,bd%nk_loc,3))
allocate(diel(data_input%spectrum_points,3))
allocate(cabs(data_input%spectrum_points,3))
allocate(a(data_input%lanczos_step))
allocate(b(data_input%lanczos_step-1))
diel=0.d0
!haydoch recursive method
! call system_clock(inizio,count_rate)
do dir=1,3
! dir=1
write(stdout,*)'Dir',dir,'constructing starting vector'
!constructing starting vector
call build_eemat(data_input,bd,eemat)
write(stdout,*)'eemat done'
if ( x%nk_loc > 0 ) then
do k=1,bd%nk_loc
do c=1,bd%numc
do v=1,bd%numv
x%avc(v,c,k) = conjg(eemat(v,c,k,dir))/(bd%en_v(v,k)-bd%en_c(c,k))
end do
end do
end do
end if
norm(dir)=x*x
write(stdout,*) 'Normalizzazione', norm(dir)
call normalize_exc(x)
!iterations
write(stdout,*)'Dir',dir,'matrix construction'
do i =1,data_input%lanczos_step
write(stdout,*) 'Step ', i
call hamiltonian(data_input, 1, bd, pp, pt, pm, x, Hx,0)
csca=x*Hx
a(i) = cmplx(dble(csca),0.d0)
if (i==1) then
etmp1=a(i)*(-1.d0,0.d0)*x
call sum_exc_sub(etmp2,Hx,etmp1)
!etmp2=Hx+etmp1
csca=etmp2*etmp2
b(i)=cmplx(sqrt(dble(csca)),0.d0)
! b(i) =cmplx(dble( sqrt((Hx+(-1.d0,0.d0)*(a(i)*x))*(Hx+(-1.d0,0.d0)*(a(i)*x)))),0.d0)
xm1 = x
etmp1=((-1.d0,0.d0)*a(i))*x
call sum_exc_sub(etmp2,Hx,etmp1)
x=(1/b(i))*etmp2
! x =(1/b(i))*(Hx+(-1.d0,0.d0)*(a(i)*x) )
else if (i ==data_input%lanczos_step) then
exit
else
etmp1=(-1.d0,0.d0)*a(i)*x
etmp2=(-1.d0,0.d0)*b(i-1)*xm1
call sum_exc_sub(etmp3,etmp1,etmp2)
! etmp3=etmp1+etmp2
call sum_exc_sub(etmp4,Hx,etmp3)
! etmp4=Hx+etmp3
!b(i)=cmplx(dble(sqrt(etmp4*etmp4)),0.d0)
b(i)=cmplx(sqrt(dble(etmp4*etmp4)),0.d0)
! b(i) =cmplx(dble( sqrt((Hx+(-1.d0,0.d0)*(a(i)*x+b(i-1)*xm1))*(Hx+(-1.d0,0.d0)*(a(i)*x+b(i-1)*xm1)))),0.d0)
! write(*,*)dir,i,'b',b(i)
etmp1=(-1.d0,0.d0)*a(i)*x
etmp2=(-1.d0,0.d0)*b(i-1)*xm1
call sum_exc_sub(etmp3,etmp1,etmp2)
call sum_exc_sub(etmp4,Hx,etmp3)
x1= (1/b(i))*etmp4
! x1 = (1/b(i))*(Hx+(-1.d0,0.d0)*(a(i)*x+b(i-1)*xm1))! (|i+1>= H|i>-a_i|i>-b_{i-1}|i-1>)/b_i
xm1 = x!|i-1> = |i>
x = x1!|i> = |i+1>
end if
write(stdout,*)dir,i
end do
write(stdout,*)'writing matrix',dir
end do
call system_clock(fine)
write(stdout,*)'elapsed time',real(fine-inizio)/real(count_rate)
if(ionode) then
iun=find_free_unit()
open( unit= iun, file='ab_coeff.dat',status='unknown')
write(iun,*) norm(1:3)
do idir=1,3
do i=1,data_input%lanczos_step-1
write(iun,*) idir, i, real(a(i))
enddo
enddo
do idir=1,3
do i=1,data_input%lanczos_step-1
write(iun,*) idir, i, real(b(i))
enddo
end do
close(iun)
endif
call deallocate_bands(bd)
call deallocate_product(pd)
call deallocate_potential(pt)
call deallocate_prod_mix(pm)
call deallocate_prod_proj(pp)
call deallocate_exc(x1)
call deallocate_exc(x)
call deallocate_exc(xm1)
call deallocate_exc(Hx)
call deallocate_exc(etmp1)
call deallocate_exc(etmp2)
call deallocate_exc(etmp3)
call deallocate_exc(etmp4)
deallocate(a)
deallocate(b)
deallocate(cabs)
deallocate(diel)
deallocate(eemat)
return
end subroutine lanczos

168
GWW/simple_bse/lanczos.ultimo.f90 Normal file
View File

@ -0,0 +1,168 @@
subroutine lanczos(data_input)
USE input_simple_exc
USE simple_objects
USE derived_objects
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode
USE constants, ONLY : rytoev
USE mp, ONLY : mp_barrier, mp_sum
USE mp_world, ONLY : world_comm
use io_files, only : prefix, tmp_dir
USE constants, ONLY : rytoev
implicit none
TYPE(input_options):: data_input
TYPE(exc):: x,Hx,xm1,x1
TYPE(epe):: element
TYPE(bands):: bd
TYPE(product) :: pd
TYPE(prod_proj) :: pp
TYPE(potential) :: pt
TYPE(prod_mix) :: pm
COMPLEX(KIND=DP) :: fraz,zeta
REAL(KIND=DP) :: omin,omax,step, cost, omega
COMPLEX(KIND=DP),ALLOCATABLE:: cabs(:,:),eemat(:,:,:,:),a(:),b(:)
INTEGER :: i,dir,v,c,k,io,iun,inizio,fine,count_rate,idir
INTEGER, EXTERNAL :: find_free_unit
REAL(KIND=DP),ALLOCATABLE::diel(:,:)
REAL(kind=DP) :: norm(3)
call initialize_product(pd)
call initialize_potential(pt)
call initialize_prod_mix(pm)
call initialize_prod_proj(pp)
call system_clock(inizio,count_rate)
call read_bands(data_input,bd)
call read_product(data_input,pd)
call read_potential(data_input,pt)
write(stdout ,*)'building derived objects'
call build_prod_mix(data_input,bd,pd,pm,pt)
call build_prod_proj(bd,pd,pp)
call setup_exc(bd,x)
call setup_exc(bd,xm1)
call setup_exc(bd,Hx)
call setup_exc(bd,x1)
data_input%lanczos_step=data_input%lanczos_step+1
allocate(eemat(bd%numv,bd%numc,bd%nk_loc,3))
allocate(diel(data_input%spectrum_points,3))
allocate(cabs(data_input%spectrum_points,3))
allocate(a(data_input%lanczos_step))
allocate(b(data_input%lanczos_step-1))
!haydoch recursive method
! call system_clock(inizio,count_rate)
! do dir=1,3
dir=1
write(stdout,*)'Dir',dir,'constructing starting vector'
!constructing starting vector
call build_eemat(data_input,bd,eemat)
write(stdout,*)'eemat done'
if ( x%nk_loc > 0 ) then
do k=1,bd%nk_loc
do c=1,bd%numc
do v=1,bd%numv
x%avc(v,c,k) = eemat(v,c,k,dir)/(bd%en_v(v,k)-bd%en_c(c,k))
!x%avc(v,c,k) = conjg(eemat(v,c,k,dir))/(bd%en_v(v,k)-bd%en_c(c,k))
end do
end do
end do
end if
norm(dir)=x*x
write(stdout,*) 'Normalizzazione', norm(dir)
call normalize_exc(x)
!iterations
write(*,*)'Dir',dir,'matrix construction'
do i =1,data_input%lanczos_step
call mp_barrier(world_comm)
call hamiltonian(data_input, 1, bd, pp, pt, pm, x, Hx)
call mp_barrier(world_comm)
a(i) = cmplx(dble(x*Hx),0.d0)
! write(*,*)dir,i,'a',a(i)
if (i==1) then
b(i) =cmplx(dble( sqrt((Hx+(-1.d0,0.d0)*(a(i)*x))*(Hx+(-1.d0,0.d0)*(a(i)*x)))),0.d0)
! write(*,*)dir,i,'b',b(i)
xm1 = x
x =(1/b(i))*(Hx+(-1.d0,0.d0)*(a(i)*x))
else if (i ==data_input%lanczos_step) then
exit
else
b(i) =cmplx(dble( sqrt((Hx+(-1.d0,0.d0)*(a(i)*x+b(i-1)*xm1))*(Hx+(-1.d0,0.d0)*(a(i)*x+b(i-1)*xm1)))),0.d0)
! write(*,*)dir,i,'b',b(i)
x1 = (1/b(i))*(Hx+(-1.d0,0.d0)*(a(i)*x+b(i-1)*xm1))! (|i+1>= H|i>-a_i|i>-b_{i-1}|i-1>)/b_i
xm1 = x!|i-1> = |i>
x = x1!|i> = |i+1>
end if
write(*,*)dir,i
end do
!construction of the continued fraction
omin = data_input%omega_min/rytoev!Ry
omax = data_input%omega_max/rytoev!Ry
step = (omax-omin)/(data_input%spectrum_points-1)
cost = - 8/3.14159/10.26!general case for V?
write(*,*)'Dir',dir,'continued fraction'
do io=1,data_input%spectrum_points
fraz=(0.d0,0.d0)
zeta =cmplx( omin + (io-1)*step, data_input%eta ) !Ry
fraz = (b(data_input%lanczos_step-1)*b(data_input%lanczos_step-1))/(zeta-a(data_input%lanczos_step))
do i=2,data_input%lanczos_step-1
fraz = (b(data_input%lanczos_step-i)*b(data_input%lanczos_step-i))/&
&(zeta - a(data_input%lanczos_step-i+1)-fraz)
diel(io,dir)=cost*aimag(1/(zeta-a(1)-fraz))!imaginary part
end do
end do
call system_clock(fine)
write(stdout,*)'elapsed time',real(fine-inizio)/real(count_rate)
!relation for the total dielectric function
!saving spectrum
if(ionode) then
iun=find_free_unit()
open( unit= iun, file='ab_coeff.dat',status='unknown')
write(iun,*) norm(1:3)
do idir=1,3
do i=1,data_input%lanczos_step-1
write(iun,*) idir, i, real(a(i))
enddo
enddo
do idir=1,3
do i=1,data_input%lanczos_step-1
write(iun,*) idir, i, real(b(i))
enddo
end do
close(iun)
endif
if(ionode) then
iun=find_free_unit()
open( unit= iun, file='spectrum_lanczos.dat',status='unknown')
do io=1,data_input%spectrum_points
omega = omin + (io-1)*step
! cabs(io,dir)=omega*aimag(diel(io))/(sqrt(4*3.14159*diel(io)))/(137)
! write(iun,*) omega*rytoev, cabs(io)
write(iun,*) omega*rytoev, (diel(io,1)+diel(io,2)+diel(io,3))/3
end do
close(iun)
end if
call deallocate_bands(bd)
call deallocate_product(pd)
call deallocate_potential(pt)
call deallocate_prod_mix(pm)
call deallocate_prod_proj(pp)
call deallocate_exc(x1)
call deallocate_exc(x)
call deallocate_exc(xm1)
call deallocate_exc(Hx)
deallocate(a)
deallocate(b)
deallocate(cabs)
deallocate(diel)
deallocate(eemat)
return
end subroutine lanczos

25
GWW/simple_bse/simple_bse.f90 Normal file
View File

@ -0,0 +1,25 @@
program simple_bse
USE start_end
USE input_simple_exc
implicit none
TYPE(input_options) :: sin
!setup MPI environment
call startup
call read_input_simple_exc( sin )
select case(sin%task)
case(0)!solve eigen-problem
call simple_eigen(sin)
case(1)!find spectrum lanczos
call lanczos(sin)
end select
call stop_run
stop
end program simple_bse

86
GWW/simple_bse/simple_eigen.f90 Normal file
View File

@ -0,0 +1,86 @@
!this subroutines calls solvers for eigen-value problem of the excitions
SUBROUTINE simple_eigen(sin)
USE input_simple_exc
USE simple_objects
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE derived_objects
implicit none
TYPE(input_options) :: sin
TYPE(exc), POINTER :: a(:)!excitons to be found
TYPE(bands) :: bd
TYPE(product) :: pd
TYPE(potential) :: pt
TYPE(prod_proj) :: pp
TYPE(prod_mix) :: pm
COMPLEX(kind=DP), ALLOCATABLE :: ene(:)!their energies
INTEGER :: i
!read in excitonic Hamiltonian stuff
! write(stdout,*) 'DEBUG 1'
call read_bands(sin,bd)
! write(stdout,*) 'DEBUG 2'
!read in product stuff
call initialize_product(pd)
call read_product(sin,pd)
! write(stdout,*) 'DEBUG 3'
!read in potential stuff
call initialize_potential(pt)
call read_potential(sin,pt)
! write(stdout,*) 'DEBUG 4'
!set up product contractions
call initialize_prod_proj(pp)
call build_prod_proj(bd,pd,pp)
! write(stdout,*) 'DEBUG 5'
!set up mixed contractions
call initialize_prod_mix(pm)
call build_prod_mix(sin,bd,pd,pm,pt)
! write(stdout,*) 'DEBUG 6'
allocate(a(sin%nvec))
do i=1,sin%nvec
call setup_exc(bd,a(i))
enddo
allocate(ene(sin%nvec))
!call solver
select case(sin%diago)
case(0)!steepest descent
call diago_exc_sd(sin,bd,pp,pt,pm,a)
case(1)!CG
call diago_exc_cg(sin,bd,pp,pt,pm,a)
end select
do i=1,sin%nvec
call nice_write_exc(bd,sin,a(i),i)
enddo
do i=1,sin%nvec
call deallocate_exc(a(i))
enddo
deallocate(a)
deallocate(ene)
call deallocate_bands(bd)
call deallocate_product(pd)
call deallocate_potential(pt)
call deallocate_prod_proj(pp)
call deallocate_prod_mix(pm)
return
END SUBROUTINE simple_eigen

794
GWW/simple_bse/simple_object.f90 Normal file
View File

@ -0,0 +1,794 @@
MODULE simple_objects
!this module describes the most important objects
USE kinds, ONLY : DP
TYPE bands
!data regarding bands and wavefunctions
INTEGER :: numv !number of valence states (those considered for excitons only)
INTEGER :: numc !number of conduction states
INTEGER :: num!numv+numc
INTEGER :: ntot_e!dimension of global to all k, basis for KS states
INTEGER :: nk!total number of k, points
INTEGER :: nk_loc!local number of k points
INTEGER :: ik_first!first local k point
INTEGER :: ik_last!last local k point
REAL(kind=DP) :: scissor!scissor in Ry a.u.
REAL(kind=DP) :: ene(4)!tot,diagonal,exchange,direct
REAL(kind=DP), DIMENSION(:,:), POINTER :: k!coordinates of local k points in atomic units
COMPLEX(kind=DP), DIMENSION(:,:,:), POINTER :: omat!overlap matrix (ntot_e,num,nk_loc)
REAL(kind=DP), DIMENSION(:,:), POINTER :: en_v!KS or GW energies of valence states (numv,nk_loc)
REAL(kind=DP), DIMENSION(:,:), POINTER :: en_c!KS or GW energies of valence states (numc,nk_loc)
END TYPE bands
TYPE exc
!this describes an electron-hole excitons
INTEGER :: numv !number of valence states (those considered for excitons only)
INTEGER :: numc !number of conduction states
INTEGER :: num!numv+numc
INTEGER :: nk!total number of k, points
INTEGER :: nk_loc!local number of k points
INTEGER :: ik_first!first local k point
INTEGER :: ik_last!last local k point
REAL(kind=DP) :: ene(4)!tot,diagonal,exchange,direct
COMPLEX(kind=DP), DIMENSION(:,:,:), POINTER :: avc!A_vck terms (numv,numc,nk_loc)
END TYPE exc
TYPE product
!this type describes the basis for the products of global wave-function basis vectors
!this is not distributed with k-points parallelization
INTEGER :: nprod_e!number of product terms
INTEGER :: ntot_e!dimension of global to all k, basis for KS states
COMPLEX(kind=DP), DIMENSION(:,:,:), POINTER :: fij!<\mathcal{E_\alpha}|(e_{i}^*e_{j}> terms
END TYPE product
TYPE potential
!this object described a potential matrix in the space of the common basis for products
INTEGER :: nprod_e!number of product terms
COMPLEX(kind=DP), DIMENSION(:,:), POINTER :: vpot
INTEGER :: nkp(3)!k-points grid equally spaced
INTEGER :: nk!total number of k, points
INTEGER, DIMENSION(:,:,:), POINTER :: ijk!correspondence k'-k-->i,j,k (3,k,k')
COMPLEX(kind=DP), DIMENSION(:,:,:,:,:), POINTER :: vpotq! V_{k'-k} (nprod_e,nprod_e,i,j,k), ijk from 1 to nkp() IMPORTANT
COMPLEX(kind=DP), DIMENSION(:,:,:,:,:), POINTER :: wpotq! Wc_{k'-k} (nprod_e,nprod_e,i,j,k), ijk from 1 to nkp() IMPORTANT
END TYPE potential
TYPE epe
INTEGER :: ntot_e
COMPLEX(KIND=DP), DIMENSION(:,:,:), POINTER :: mat
END TYPE epe
INTERFACE OPERATOR(+)
MODULE PROCEDURE sum_exc
END INTERFACE
INTERFACE OPERATOR(*)
MODULE PROCEDURE prod_exc
MODULE PROCEDURE prod_c_exc
END INTERFACE
INTERFACE ASSIGNMENT(=)
MODULE PROCEDURE assign_exc
END INTERFACE
INTERFACE OPERATOR(.hd.)
MODULE PROCEDURE h_diag
END INTERFACE
CONTAINS
!!!!!!!
! INTEGER :: nk_loc
! INTEGER :: ik_first
! INTEGER :: ik_last
! COMPLEX(kind=DP), DIMENSION(:,:,:), POINTER :: avc
! END TYPE exc
SUBROUTINE nice_write_exc(bd,simple_in,a,label)
!write in characters an exciton, with label
USE input_simple_exc, ONLY : input_options
USE mp, ONLY : mp_sum
USE mp_world, ONLY : world_comm
USE io_files, ONLY : tmp_dir, prefix
USE io_global, ONLY : ionode_id, ionode, stdout
USE mp_world, ONLY : mpime, nproc
implicit none
TYPE(bands) :: bd
TYPE(input_options) :: simple_in
TYPE(exc) :: a
INTEGER :: label
INTEGER, EXTERNAL :: find_free_unit
INTEGER :: iun
CHARACTER(4) :: nfile
INTEGER :: ik, ikl,iv,ic
REAL(kind=DP) :: kappa(3)
COMPLEX(kind=DP), ALLOCATABLE :: avc(:,:)
allocate(avc(a%numv,a%numc))
write(nfile,'(4i1)') &
& label/1000,mod(label,1000)/100,mod(label,100)/10,mod(label,10)
if(ionode) then
iun = find_free_unit()
open( unit=iun, file=trim(tmp_dir)//trim(simple_in%prefix)//'.exciton.'//nfile, status='unknown')
write(iun,*) a%numv,a%numc,a%nk
write(iun,*) a%ene(1:4)
end if
do ik= 1,a%nk
kappa=0.d0
do ikl=a%ik_first,a%ik_last
if(ikl==ik) then
kappa(1:3)=bd%k(1:3,ik-a%ik_first+1)
endif
enddo
call mp_sum(kappa,world_comm)
if(ionode) write(iun,*) kappa(1:3)
enddo
do ik= 1,a%nk
avc=(0.d0,0.d0)
do ikl=a%ik_first,a%ik_last
if(ikl==ik) then
avc(1:a%numv,1:a%numc)=a%avc(1:a%numv,1:a%numc,ik-a%ik_first+1)
endif
enddo
call mp_sum(avc,world_comm)
if(ionode) then
do ic=1,a%numc
do iv=1,a%numv
write(iun,*) avc(iv,ic)
enddo
enddo
endif
enddo
if(ionode) close(iun)
deallocate(avc)
return
END SUBROUTINE nice_write_exc
!!!!!!
subroutine initialize_epe(element)
implicit none
TYPE(epe) :: element
nullify(element%mat)
return
end subroutine initialize_epe
subroutine deallocate_epe(element)
implicit none
TYPE(epe) :: element
if(associated(element%mat)) deallocate(element%mat)
nullify(element%mat)
return
end subroutine deallocate_epe
subroutine read_epe(simple_in,element)
USE input_simple_exc, ONLY : input_options
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE io_files, ONLY : tmp_dir, prefix
USE io_global, ONLY : ionode_id, ionode, stdout
USE mp_world, ONLY : mpime, nproc
implicit none
TYPE(input_options) :: simple_in
TYPE(epe) :: element
INTEGER, EXTERNAL :: find_free_unit
INTEGER :: iun,i,a
write(*,*)'epe:opening file'
if(ionode) then
iun = find_free_unit()
open( unit=iun, file=trim(tmp_dir)//trim(simple_in%prefix)//'.epe', status='old',form='unformatted')
write(*,*)'file opened'
read(iun) element%ntot_e
end if
call mp_bcast(element%ntot_e,ionode_id,world_comm)
write(*,*)element%ntot_e
allocate(element%mat(element%ntot_e,element%ntot_e,3))
if(ionode) then
do a=1,3
do i=1,element%ntot_e
read(iun) element%mat(1:element%ntot_e,i,a)
end do
end do
close(iun)
end if
call mp_bcast(element%mat,ionode_id,world_comm)
write(*,*)'epe:read all'
return
end subroutine read_epe
subroutine deallocate_bands(bd)
implicit none
TYPE(bands) :: bd
if(associated(bd%k)) deallocate(bd%k)
if(associated(bd%omat)) deallocate(bd%omat)
if(associated(bd%en_v)) deallocate(bd%en_v)
if(associated(bd%en_c)) deallocate(bd%en_c)
nullify(bd%k)
nullify(bd%omat)
nullify(bd%en_v)
nullify(bd%en_c)
return
end subroutine deallocate_bands
subroutine deallocate_exc(a)
implicit none
TYPE(exc) :: a
if(associated(a%avc)) deallocate(a%avc)
nullify(a%avc)
return
end subroutine deallocate_exc
subroutine deallocate_product(pd)
implicit none
TYPE(product) :: pd
if(associated(pd%fij)) deallocate(pd%fij)
nullify(pd%fij)
return
end subroutine deallocate_product
subroutine initialize_product(pd)
implicit none
TYPE(product) :: pd
nullify(pd%fij)
end subroutine initialize_product
subroutine initialize_potential(pt)
implicit none
TYPE(potential) :: pt
nullify(pt%vpot)
nullify(pt%vpotq)
nullify(pt%ijk)
return
end subroutine initialize_potential
subroutine deallocate_potential(pt)
implicit none
TYPE(potential) :: pt
if(associated(pt%vpot)) deallocate(pt%vpot)
nullify(pt%vpot)
if(associated(pt%vpotq)) deallocate(pt%vpotq)
nullify(pt%vpotq)
if(associated(pt%wpotq)) deallocate(pt%wpotq)
nullify(pt%wpotq)
if(associated(pt%ijk)) deallocate(pt%ijk)
nullify(pt%ijk)
return
end subroutine deallocate_potential
subroutine read_bands(simple_in,bd)
!this subroutine read in the bands structure from disk
!and distribute it among processors
USE input_simple_exc, ONLY : input_options
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE io_files, ONLY : tmp_dir
USE io_global, ONLY : ionode_id, ionode, stdout
USE mp_world, ONLY : mpime, nproc
USE constants, ONLY : rytoev
implicit none
INTEGER, EXTERNAL :: find_free_unit
TYPE(input_options) :: simple_in
TYPE(bands) :: bd
INTEGER :: iun,ik,i
INTEGER :: l_blk
REAL(kind=DP) :: xk(3)
REAL(kind=DP), allocatable :: et(:)
COMPLEX(kind=DP), allocatable :: omat(:,:)
if(ionode) then
iun = find_free_unit()
open( unit=iun, file=trim(tmp_dir)//trim(simple_in%prefix)//'.wfc_basis', status='old',form='unformatted')
read(iun) bd%nk
read(iun) bd%numv
read(iun) bd%numc
read(iun) bd%ntot_e
end if
call mp_bcast(bd%nk, ionode_id, world_comm)
call mp_bcast(bd%numv, ionode_id, world_comm)
call mp_bcast(bd%numc, ionode_id, world_comm)
call mp_bcast(bd%ntot_e, ionode_id, world_comm)
write(stdout,*) 'NUMBER OF K POINTS : ', bd%nk
write(stdout,*) 'NUMBER OF VALENCE STATES : ', bd%numv
write(stdout,*) 'NUMBER OF CONDUCTION STATES : ', bd%numc
write(stdout,*) 'NUMBER OF GLOBAL STATES : ', bd%ntot_e
bd%num=bd%numv+bd%numc
l_blk=bd%nk/nproc
if(l_blk*nproc<bd%nk) l_blk=l_blk+1
if(l_blk*mpime+1 <= bd%nk) then
bd%ik_first=l_blk*mpime+1
bd%ik_last=bd%ik_first+l_blk-1
if(bd%ik_last>bd%nk) bd%ik_last=bd%nk
bd%nk_loc=bd%ik_last-bd%ik_first+1
else
bd%nk_loc=0
bd%ik_first=0
bd%ik_last=-1
endif
! write(stdout,*) 'DEBUG nk_loc :', bd%nk_loc
!nk_loc ik_first ik_last
if(bd%nk_loc>0) then
allocate(bd%k(3,bd%nk_loc))
allocate(bd%omat(bd%ntot_e,bd%num,bd%nk_loc))
allocate(bd%en_v(bd%numv,bd%nk_loc))
allocate(bd%en_c(bd%numc,bd%nk_loc))
else
nullify(bd%k)
nullify(bd%omat)
nullify(bd%en_v)
nullify(bd%en_c)
endif
allocate(et(bd%num))
allocate(omat(bd%ntot_e,bd%num))
do ik=1,bd%nk
if(ionode) then
read(iun) xk(1:3)
read(iun) et(1:bd%num)
do i=1,bd%num
read(iun) omat(1:bd%ntot_e,i)
enddo
endif
call mp_bcast(xk, ionode_id, world_comm )
call mp_bcast(et, ionode_id, world_comm )
call mp_bcast(omat, ionode_id, world_comm )
if(ik>=bd%ik_first .and. ik <= bd%ik_last) then
bd%k(1:3,ik-bd%ik_first+1)=xk(1:3)
bd%en_v(1:bd%numv,ik-bd%ik_first+1)=et(1:bd%numv)
bd%en_c(1:bd%numc,ik-bd%ik_first+1)=et(bd%numv+1:bd%num)
bd%omat(1:bd%ntot_e,1:bd%num,ik-bd%ik_first+1)=omat(1:bd%ntot_e,1:bd%num)
endif
enddo
if(ionode) close(iun)
deallocate(et,omat)
bd%scissor=simple_in%scissor/rytoev
return
end subroutine read_bands
subroutine read_product(simple_in,pd)
!read product terms from disk
USE input_simple_exc, ONLY : input_options
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE io_files, ONLY : tmp_dir
USE io_global, ONLY : ionode_id, ionode, stdout
USE mp_world, ONLY : mpime, nproc
implicit none
TYPE(input_options) :: simple_in
TYPE(product) :: pd
INTEGER, EXTERNAL :: find_free_unit
INTEGER :: iun
INTEGER :: ii,jj
if(ionode) then
iun = find_free_unit()
open( unit=iun, file=trim(tmp_dir)//trim(simple_in%prefix)//'.product_basis', status='old',form='unformatted')
read(iun) pd%nprod_e
read(iun) pd%ntot_e
end if
call mp_bcast(pd%nprod_e, ionode_id, world_comm)
call mp_bcast(pd%ntot_e, ionode_id, world_comm)
write(stdout,*) 'NUMBER OF PRODUCTS : ', pd%nprod_e
write(stdout,*) 'NUMBER OF GLOBAL STATES : ', pd%ntot_e
allocate(pd%fij(pd%nprod_e,pd%ntot_e,pd%ntot_e))
!note the order important for index corresponding to complex conjugate
if(ionode) then
do ii=1,pd%ntot_e
do jj=1,pd%ntot_e
read(iun) pd%fij(1:pd%nprod_e,ii,jj)
enddo
enddo
close(iun)
endif
call mp_bcast(pd%fij, ionode_id, world_comm)
return
end subroutine read_product
subroutine read_potential(simple_in,pt)
USE input_simple_exc, ONLY : input_options
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE io_files, ONLY : tmp_dir
USE io_global, ONLY : ionode_id, ionode, stdout
USE mp_world, ONLY : mpime, nproc
implicit none
INTEGER, EXTERNAL :: find_free_unit
TYPE(input_options) :: simple_in
TYPE(potential) :: pt
INTEGER :: iun,ii,ik,jk,kk
INTEGER :: nktot
if(ionode) then
iun = find_free_unit()
open( unit=iun, file=trim(tmp_dir)//trim(simple_in%prefix)//'.v_mat0', status='old',form='unformatted')
read(iun) pt%nprod_e
end if
call mp_bcast(pt%nprod_e, ionode_id, world_comm)
allocate(pt%vpot(pt%nprod_e,pt%nprod_e))
if(ionode) then
do ii=1,pt%nprod_e
read(iun) pt%vpot(1:pt%nprod_e,ii)
enddo
endif
call mp_bcast(pt%vpot, ionode_id, world_comm)
if(ionode) read(iun) pt%nkp(1:3)
call mp_bcast(pt%nkp, ionode_id, world_comm)
pt%nk=pt%nkp(1)*pt%nkp(2)*pt%nkp(3)
allocate(pt%ijk(3,pt%nk,pt%nk),pt%vpotq(pt%nprod_e,pt%nprod_e,pt%nkp(1),pt%nkp(2),pt%nkp(3)))
allocate(pt%wpotq(pt%nprod_e,pt%nprod_e,pt%nkp(1),pt%nkp(2),pt%nkp(3)))
if(ionode) then
do ik=1,pt%nk!on k'
do jk=1,pt%nk! on k
read(iun) pt%ijk(1:3,jk,ik)
enddo
enddo
endif
call mp_bcast(pt%ijk, ionode_id, world_comm)
do ik=1,pt%nkp(1)
do jk=1,pt%nkp(2)
do kk=1,pt%nkp(3)
if(ionode) then
do ii=1,pt%nprod_e
read(iun) pt%vpotq(1:pt%nprod_e,ii,ik,jk,kk)
enddo
endif
call mp_bcast(pt%vpotq(:,:,ik,jk,kk), ionode_id, world_comm)
enddo
enddo
enddo
do ik=1,pt%nkp(1)
do jk=1,pt%nkp(2)
do kk=1,pt%nkp(3)
if(ionode) then
do ii=1,pt%nprod_e
read(iun) pt%wpotq(1:pt%nprod_e,ii,ik,jk,kk)
enddo
endif
call mp_bcast(pt%wpotq(:,:,ik,jk,kk), ionode_id, world_comm)
enddo
enddo
enddo
if(ionode) close(iun)
return
end subroutine read_potential
subroutine setup_exc(bd,a)
!from bands object sets up exciton a
implicit none
TYPE(bands), INTENT(in) :: bd
TYPE(exc), INTENT(out) :: a
a%numv=bd%numv
a%numc=bd%numc
a%num=bd%num
a%nk=bd%nk
a%nk_loc=bd%nk_loc
a%ik_first=bd%ik_first
a%ik_last=bd%ik_last
if(a%nk_loc>0) then
allocate(a%avc(a%numv,a%numc,a%nk_loc))
else
nullify(a%avc)
endif
return
end subroutine setup_exc
FUNCTION sum_exc(a,b)
USE io_global, ONLY : stdout
! a+b
implicit none
TYPE(exc), INTENT(in) :: a,b
TYPE(exc) :: sum_exc
!check for consistency
if(a%numv/=b%numv .or. a%numc/=b%numc.or. a%num/=b%num .or. a%nk/=b%nk &
&.or.a%nk_loc/=b%nk_loc .or. a%ik_first/=b%ik_first .or. a%ik_last/=b%ik_last) then
write(stdout,*) 'Problem with sum_exc: inconsistency'
stop
endif
sum_exc%numv=b%numv
sum_exc%numc=b%numc
sum_exc%num=b%num
sum_exc%nk=b%nk
sum_exc%nk_loc=b%nk_loc
sum_exc%ik_first=b%ik_first
sum_exc%ik_last=b%ik_last
if(associated(sum_exc%avc)) deallocate(sum_exc%avc)
if(sum_exc%nk_loc>0) then
allocate(sum_exc%avc(sum_exc%numv,sum_exc%numc,sum_exc%nk_loc))
else
nullify(sum_exc%avc)
endif
if(a%nk_loc>0) then
sum_exc%avc(1:a%numv,1:a%numc,1:a%nk_loc)=a%avc(1:a%numv,1:a%numc,1:a%nk_loc)+&
&b%avc(1:a%numv,1:a%numc,1:a%nk_loc)
endif
return
END FUNCTION sum_exc
SUBROUTINE sum_exc_sub(sum_exc, a,b)
USE io_global, ONLY : stdout
! a+b
implicit none
TYPE(exc), INTENT(in) :: a,b
TYPE(exc) :: sum_exc
!check for consistency
if(a%numv/=b%numv .or. a%numc/=b%numc.or. a%num/=b%num .or. a%nk/=b%nk &
&.or.a%nk_loc/=b%nk_loc .or. a%ik_first/=b%ik_first .or. a%ik_last/=b%ik_last) then
write(stdout,*) 'Problem with sum_exc: inconsistency'
stop
endif
sum_exc%numv=b%numv
sum_exc%numc=b%numc
sum_exc%num=b%num
sum_exc%nk=b%nk
sum_exc%nk_loc=b%nk_loc
sum_exc%ik_first=b%ik_first
sum_exc%ik_last=b%ik_last
if(associated(sum_exc%avc)) deallocate(sum_exc%avc)
if(sum_exc%nk_loc>0) then
allocate(sum_exc%avc(sum_exc%numv,sum_exc%numc,sum_exc%nk_loc))
else
nullify(sum_exc%avc)
endif
if(a%nk_loc>0) then
sum_exc%avc(1:a%numv,1:a%numc,1:a%nk_loc)=a%avc(1:a%numv,1:a%numc,1:a%nk_loc)+&
&b%avc(1:a%numv,1:a%numc,1:a%nk_loc)
endif
return
END SUBROUTINE sum_exc_sub
FUNCTION prod_exc(a,b)
USE io_global, ONLY : stdout
USE mp, ONLY : mp_sum
USE mp_world, ONLY : world_comm
! scalar produc
implicit none
TYPE(exc), INTENT(in) :: a,b
COMPLEX(kind=DP) :: prod_exc
COMPLEX(kind=DP), EXTERNAL :: zdotc
!check for consistency
if(a%numv/=b%numv .or. a%numc/=b%numc.or. a%num/=b%num .or. a%nk/=b%nk &
&.or.a%nk_loc/=b%nk_loc .or. a%ik_first/=b%ik_first .or. a%ik_last/=b%ik_last) then
write(*,*) 'Problem with prod_exc: inconsistency'
stop
endif
if(a%nk_loc>0) then
prod_exc= zdotc(a%numv*a%numc*a%nk_loc,a%avc,1,b%avc,1)
else
prod_exc=(0.d0,0.d0)
endif
call mp_sum(prod_exc,world_comm)
return
END FUNCTION prod_exc
FUNCTION prod_c_exc(z,a)
!complex * exciton
USE io_global, ONLY : stdout
USE mp, ONLY : mp_sum
USE mp_world, ONLY : world_comm
implicit none
TYPE(exc), INTENT(in) :: a
COMPLEX(kind=DP), INTENT(in) :: z
TYPE(exc) :: prod_c_exc
prod_c_exc%numv=a%numv
prod_c_exc%numc=a%numc
prod_c_exc%num=a%num
prod_c_exc%nk=a%nk
prod_c_exc%nk_loc=a%nk_loc
prod_c_exc%ik_first=a%ik_first
prod_c_exc%ik_last=a%ik_last
!if(associated(prod_c_exc%avc)) deallocate(prod_c_exc%avc)
if(prod_c_exc%nk_loc>0) then
allocate(prod_c_exc%avc(prod_c_exc%numv,prod_c_exc%numc,prod_c_exc%nk_loc))
else
nullify(prod_c_exc%avc)
endif
if(prod_c_exc%nk_loc>0) then
prod_c_exc%avc(1:a%numv,1:a%numc,1:a%nk_loc)=z*a%avc(1:a%numv,1:a%numc,1:a%nk_loc)
endif
return
END FUNCTION prod_c_exc
SUBROUTINE assign_exc(a,b)
!x=a
implicit none
TYPE(exc), INTENT(out) :: a
TYPE(exc),INTENT(in) ::b
a%numv=b%numv
a%numc=b%numc
a%num=b%num
a%nk=b%nk
a%nk_loc=b%nk_loc
a%ik_first=b%ik_first
a%ik_last=b%ik_last
if(associated(a%avc)) deallocate(a%avc)
if(a%nk_loc>0) then
allocate(a%avc(a%numv,a%numc,a%nk_loc))
else
nullify(a%avc)
endif
if(a%nk_loc>0) then
a%avc(1:a%numv,1:a%numc,1:a%nk_loc)=b%avc(1:a%numv,1:a%numc,1:a%nk_loc)
endif
return
END SUBROUTINE assign_exc
SUBROUTINE randomize_exc(a)
!this subroutine set an exc object randomly
USE random_numbers, ONLY : randy
implicit none
TYPE(exc), INTENT(inout) :: a
INTEGER :: i,j,k
if(a%nk_loc > 0) then
if(a%nk_loc > 0 ) then
do i=1,a%numv
do j=1,a%numc
do k=1,a%nk_loc
a%avc(i,j,k)=dcmplx(randy(),randy())
enddo
enddo
enddo
endif
endif
return
END SUBROUTINE randomize_exc
SUBROUTINE normalize_exc(a)
!normalize the avc vector
implicit none
TYPE(exc), INTENT(inout) :: a
COMPLEX(kind=DP) :: csca
REAL(kind=DP) :: sca
csca=a*a
sca=dble(csca)
csca=cmplx(1.d0/sqrt(sca),0.d0)
if(a%nk_loc>0) a%avc=csca*a%avc
return
END SUBROUTINE normalize_exc
FUNCTION h_diag(bd,a)
!this function applies the diagonal part of the excitonic Hamiltonian
!TO DO SCISSOR OR SIMILA
USE io_global, ONLY : stdout,ionode
implicit none
TYPE(exc) :: h_diag
TYPE(exc), INTENT(in) :: a
TYPE(bands), INTENT(in) :: bd
INTEGER :: iv,ic,ik
h_diag%numv=a%numv
h_diag%numc=a%numc
h_diag%num=a%num
h_diag%nk=a%nk
h_diag%nk_loc=a%nk_loc
h_diag%ik_first=a%ik_first
h_diag%ik_last=a%ik_last
if(h_diag%nk_loc>0) then
!if(associated(h_diag%avc)) deallocate(h_diag%avc)
allocate(h_diag%avc(h_diag%numv,h_diag%numc,h_diag%nk_loc))
else
nullify(h_diag%avc)
endif
if(ionode) then
! write(stdout,*) 'DEBUG:',h_diag%nk_loc,h_diag%numc,h_diag%numv,bd%en_c(1,1)
! write(stdout,*) 'DEBUG:',bd%nk_loc,bd%numc,bd%numv,bd%en_v(1,1)
endif
if(h_diag%nk_loc>0) then
do ik=1,h_diag%nk_loc
do ic=1,h_diag%numc
do iv=1,h_diag%numv
h_diag%avc(iv,ic,ik)=(bd%en_c(ic,ik)-bd%en_v(iv,ik)+bd%scissor)*a%avc(iv,ic,ik)
enddo
enddo
enddo
endif
return
END FUNCTION h_diag
END MODULE simple_objects

134
GWW/simple_bse/spectrum.f90 Normal file
View File

@ -0,0 +1,134 @@
subroutine spectrum(data_input,x,bd,energy)
USE input_simple_exc
USE simple_objects
USE derived_objects
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode
USE constants, ONLY : rytoev
USE mp, ONLY : mp_barrier, mp_sum
USE mp_world, ONLY : world_comm,mpime
use io_files, only : prefix, tmp_dir
USE constants, ONLY : rytoev
implicit none
TYPE(input_options):: data_input
TYPE(exc):: x(data_input%nvec)
TYPE(bands):: bd
COMPLEX(KIND=DP) :: energy(data_input%nvec)
COMPLEX(KIND=DP), ALLOCATABLE :: diel(:),D(:,:),N(:),Etempc(:,:,:),Etempv(:,:,:),eemat(:,:,:,:),mate(:,:,:,:),coeff(:,:), cabs(:)
REAL(KIND=DP)::step,omega,cost,omin,omax, norm
INTEGER :: i,j,v,c,k,S,a,kk,io,iun
INTEGER, EXTERNAL :: find_free_unit
allocate(diel(data_input%spectrum_points))!costante dielettrica(omega)
allocate(D(data_input%nvec,data_input%spectrum_points))!denominatore(S,omega)
allocate(N(data_input%nvec))!numeratore(S)
allocate(eemat(bd%numv,bd%numc,bd%nk_loc,3))!E_{iv}^k^* E_{jc}^k <e_i|p_a|e_j>
allocate(coeff(3,data_input%nvec))!!elemento con coeff fissati dir,S
allocate(cabs(data_input%spectrum_points))
!comodità
!energy(1) =2.49911d-2
!energy(2) =2.57741d-2
!energy(3) =2.59408d-2
!energy(4) =2.64996d-2
!energy(5) =2.66287d-2
!energy(6) =2.68860d-2
!energy(7) =3.38354d-2
!energy(8) =3.52200d-2
!energy(9) =3.60030d-2
!energy(10) =3.62000d-2
!energy(11) =3.73634d-2
!energy(12) =3.88258d-2
!energy(13) =3.74831d-2
!energy(14) =3.93831d-2
!energy(15) =3.98621d-2
!energy(16) =4.02376d-2
write(stdout,*)'Energies eV'
do j=1,data_input%nvec
write(stdout,*) energy(j)*rytoev
end do
write(stdout,*) 'Calculating spectrum'
call build_eemat(data_input,bd,eemat)
!costruzione coefficienti*elemento
do a=1,3
do S=1,data_input%nvec
coeff(a,S)=(0.d0,0.d0)
do k=1,bd%nk_loc
do v=1,bd%numv
do c=1,bd%numc
coeff(a,S) = coeff(a,S) + x(S)%avc(v,c,k)*eemat(v,c,k,a)
! norm = norm + 1
! write(*,*)coeff(a,S)
end do
end do
end do
end do
end do
write(*,*)'coeff'
call mp_sum(coeff,world_comm)
do S=1,data_input%nvec
N(S)=(0.d0,0.d0)
do a=1,3
N(S) = N(S) + conjg(coeff(a,S))*coeff(a,S)
! write(*,*)N(S)
end do
N(S)=N(S)/3
write(*,*)N(S)
end do
write(*,*)'numerator'
!denominator
omin = data_input%omega_min/rytoev!Ry
omax = data_input%omega_max/rytoev!Ry
! write(*,*)'omega_max',omax,data_input%omega_max
step = (omax-omin)/(data_input%spectrum_points-1)
! write(*,*)'step',step
do io=1,data_input%spectrum_points
omega = omin + (io-1)*step!Ry
do S=1,data_input%nvec
! write(*,*)energy(S),omega
D(S,io) = energy(S)*( energy(S)*energy(S) - ( cmplx(omega,data_input%eta))&
&*(cmplx(omega,data_input%eta)))/data_input%eta
! write(*,*)omega,energy(S),D(S,io)
end do
end do
write(*,*)'denominator'
!constant
!
cost = 16/3.14159/10.26!general case for V?
!epsilon
do io=1,data_input%spectrum_points
diel(io) = (0.d0,0.d0)
do S = 1,data_input%nvec
diel(io) = diel(io) + N(S)/D(S,io)
end do
diel(io) = 1 + diel(io)*cost!all quantities in Rydberg atomic units
! write(*,*)io, diel(io)
end do
write(*,*)'epsilon'
!saving spectrum
if(ionode) then
iun=find_free_unit()
open( unit= iun, file='spectrum.txt',status='unknown')
do io=1,data_input%spectrum_points
omega = omin + (io-1)*step
cabs(io)=omega*aimag(diel(io))/(sqrt(4*3.14159*diel(io)))/(137)
! write(iun,*) omega*rytoev, cabs(io)
write(iun,*)omega*rytoev, aimag(diel(io))
end do
close(iun)
end if
write(*,*)'saved epsilon'
deallocate(diel)
deallocate(N)
deallocate(D)
deallocate(eemat)
deallocate(coeff)
deallocate(cabs)
return
end subroutine spectrum

68
GWW/simple_bse/start_end.f90 Normal file
View File

@ -0,0 +1,68 @@
!
! Copyright (C) 2001-2013 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
MODULE start_end
!this module contains routines to initialize the MPI environment
IMPLICIT NONE
CHARACTER (len=10), PARAMETER :: code = 'SIMPLE_BSE'
#if defined(_OPENMP)
INTEGER, SAVE :: ntids
#endif
CONTAINS
SUBROUTINE startup
!
USE io_global, ONLY : stdout, ionode
USE mp_world, ONLY : nproc
USE mp_global, ONLY : mp_startup
USE environment, ONLY: environment_start
IMPLICIT NONE
#if defined(__MPI)
CALL mp_startup()
#endif
CALL environment_start ( code )
#if defined(__MPI)
if(ionode) then
write(stdout,*) 'MPI PARALLEL VERSION'
write(stdout,*) 'Number of procs: ', nproc
write(stdout,*) 'SIMPLE_BSE: Version 1.00'
endif
#else
write(stdout,*) 'SIMPLE_BSE: Version 1.00'
#endif
return
END SUBROUTINE startup
SUBROUTINE stop_run
!this subroutine kills the MPI environment
USE io_global, ONLY : stdout, ionode
USE mp_global, ONLY : mp_global_end
IMPLICIT NONE
#if defined(__MPI)
if(ionode) write(stdout,*) 'Stopping MPI environment'
call mp_global_end( )
#endif
return
END SUBROUTINE stop_run
END MODULE start_end