!
!----------------------------------------------------------------
subroutine dielec_mat ( xxq, nw )
!----------------------------------------------------------------
!
! the static dielectric function eps(G=0,G'=0,q,w=0)
!
!----------------------------------------------------------------
!
use parameters
use constants
use gspace
use kspace
implicit none
!
integer :: iq, count, i, j, k, ipol, ikk, ikq, nw
integer :: ig, iw, igp, ierr, ibnd, ios, recl, unf_recl, notconv
integer :: ideltag, ig1, ig2, ig3, ik, jbnd
real(dbl) :: xxq(3), ui, uj, uk, qg2, xxk(3), et(nbnd, nks), epsil, vqg, w
real(dbl) :: g2kin(ngm), arg, deltag(3), avg_iter, precondition (ngm)
real(DP), allocatable :: enval(:), u(:,:), fv1(:), fv2(:), hk(:,:), ss(:), vs(:)
complex(dbl) :: vr(nr), psi(ngm, nbnd), dvbare(nr), dvscf(nr), ZDOTC, chi0
complex(kind=DP) :: evc (ngm, nbnd), evq (ngm, nbnd), a(nksq,nbnd_occ,nbnd)
logical :: equiv
!
allocate ( enval(ngm0), u(ngm0,ngm0), fv1(ngm0), fv2(ngm0), &
ss(ngm), vs(ngm), hk (ngm0, ngm0) )
!
allocate ( xk (3,nks), wk(nks) )
!
recl = 2 * nbnd * ngm ! 2 stands for complex
unf_recl = DIRECT_IO_FACTOR * recl
open ( iunwfc, file = "./silicon.wfc", iostat = ios, form = 'unformatted', &
status = 'unknown', access = 'direct', recl = unf_recl)
!
! the empirical pseudopotential
!
vs = zero
! integer comparison - careful with other structures
do ig = 1, ngm
if ( int ( gl(igtongl(ig)) ) .eq. 3 ) then
vs (ig) = v3
elseif ( int ( gl(igtongl(ig)) ) .eq. 8 ) then
vs (ig) = v8
elseif ( int ( gl(igtongl(ig)) ) .eq. 11 ) then
vs (ig) = v11
endif
enddo
!
! the structure factor
!
do ig = 1, ngm
arg = twopi * ( g(1,ig) * tau( 1, 1) + g(2,ig) * tau( 2, 1) + g(3,ig) * tau( 3, 1) )
ss (ig) = cos ( arg )
enddo
!
! the empirical pseudopotential in real space
! for further use in h_psi
!
vr = czero
do ig = 1, ngm
vr ( nl ( ig ) ) = dcmplx ( ss (ig) * vs (ig), zero )
enddo
call cfft3 ( vr, nr1, nr2, nr3, 1)
!
!
! generate uniform {k} and {k+q} grids - no symmetry-reduction for now
! (MP method)
!
count = 0
do i = 1, nq1
ui = (q1 + 2.d0 * i - nq1 - 1.d0) / (2.d0 * nq1)
do j = 1, nq2
uj = (q2 + 2.d0 * j - nq2 - 1.d0) / (2.d0 * nq2)
do k = 1, nq3
uk = (q3 + 2.d0 * k - nq3 - 1.d0) / (2.d0 * nq3)
count = count + 1
xk (:, count) = ui * bg(:,1) + uj * bg(:,2) + uk * bg(:,3)
enddo
enddo
enddo
!
! include spin degeneracy
wk = 2.d0 / float ( count )
!
! double the grid and add k+q in the even positions
!
do ik = nksq, 1, -1
ikk = 2 * ik - 1
ikq = 2 * ik
xk(:,ikk) = xk(:,ik)
xk(:,ikq) = xk(:,ik)+xxq
wk(ikk) = wk(ik)
wk(ikq) = 0.d0
enddo
!
! write(6,'(/4x,a)') repeat('-',67)
! write(6,'(4x,"Uniform k-point grid for the screened coulomb interaction"/)')
! do ik = 1, nks
! write ( 6, '(4x,"k(",i4,") = (",3f12.7,"), wk =",f12.7)') &
! ik, (xk (ipol, ik) , ipol = 1, 3) , wk (ik)
! enddo
! write(6,'(4x,a/)') repeat('-',67)
!
do ik = 1, nks
!
if (ik/100*100-ik.eq.0) write(6,*) ik,nks
!
! find wavefunctions for this kpoint
!
xxk = xk(:, ik)
!
! the k-dependent kinetic energy in Ry
!
do ig = 1, ngm
g2kin(ig) = ( (xxk(1) + g(1,ig))**2.d0 + &
(xxk(2) + g(2,ig))**2.d0 + &
(xxk(3) + g(3,ig))**2.d0 ) * tpiba2
enddo
!
! starting guess from direct diagonalization - low cutoff
!
do ig = 1, ngm0
hk ( ig, ig) = g2kin ( ig )
enddo
!
do ig1 = 1, ngm0
do ig2 = ig1+1, ngm0
!
! define ideltag
!
deltag = g(:,ig1) - g(:,ig2)
ideltag = 1
do while ( .not. equiv ( deltag, g(:,ideltag) ) .and. ideltag.lt.ngm )
ideltag = ideltag + 1
enddo
!
if (ideltag.ne.ngm) then
hk ( ig1, ig2) = ss (ideltag) * vs (ideltag)
hk ( ig2, ig1) = hk ( ig1, ig2)
endif
!
enddo
enddo
!
call rs ( ngm0, ngm0, hk, enval, 1, u, fv1, fv2, ierr)
psi = czero
do ibnd = 1, nbnd
do ig = 1, ngm0
psi(ig,ibnd) = dcmplx ( u(ig,ibnd), zero)
enddo
enddo
!
! conjugate gradients diagonalization
!
precondition = max( 1.d0, g2kin )
!
call ccgdiagg (ngm, ngm, nbnd, psi, et(:,ik), precondition, eps, &
maxter, .true., notconv, avg_iter, g2kin, vr)
!
!
! write(6,'(4x,i3,8(2x,f6.3))') ik, et( 1: 8,ik)*ryd2ev
! write(6,'(4x,3x,8(2x,f6.3))') et( 9:16,ik)*ryd2ev
! write(6,'(4x,3x,8(2x,f6.3))') et(17:24,ik)*ryd2ev
! write(6,'(4x,a/)') repeat('-',67)
!
! direct write to file
!
write ( iunwfc, rec = ik, iostat = ios) psi
!
! to read: read ( iunwfc, rec = ik, iostat = ios) psi
!
enddo
!
do ik = 1, nksq
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! reads unperturbed wavefuctions u(k) and u(k+q)
!
read ( iunwfc, rec = ikk, iostat = ios) evc
read ( iunwfc, rec = ikq, iostat = ios) evq
!
! calculate the matrix elements <c,k+q|v,k> (G=0,G'=0)
!
do ibnd = 1, nbnd_occ
do jbnd = nbnd_occ+1, nbnd
a(ik,ibnd,jbnd) = ZDOTC(ngm, evq(:,jbnd), 1, evc(:,ibnd), 1)
enddo
enddo
!
enddo
!
! the dielectric matrix (G=0,G'=0,w=0)
!
! the factor 2 takes into accout the (c,v) and (v,c) orderings of the mat.
! elements. the crystal volume & spin factor is 2/(nksq*omega), with 2/nks
! being given by wk(ikk)
!
qg2 = xxq(1)**2.d0 + xxq(2)**2.d0 + xxq(3)**2.d0
vqg = fpi * e2 / omega / (tpiba2 * qg2)
!
do iw = 1, nw
!
w = fmin + ( wmax - wmin ) * float(iw-1) / float(nw-1) / ryd2ev
!
chi0 = czero
do ik = 1, nksq
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
do ibnd = 1, nbnd_occ
do jbnd = nbnd_occ+1, nbnd
chi0 = chi0 + 2.d0 * wk(ikk) * 0.5d0 &
* conjg ( a(ik,ibnd,jbnd) ) * a(ik,ibnd,jbnd) * &
( cone / ( et (ibnd,ikk) - et (jbnd,ikq) + ci * eta - w) &
+ cone / ( et (ibnd,ikk) - et (jbnd,ikq) + ci * eta + w) )
enddo
enddo
enddo
!
! epsil = dreal ( one - vqg * chi0 )
epsil = dreal ( one / ( one - vqg * chi0 ) )
!
! write(6,'(2f12.5)') xxq(1), epsil
!
write(6, *) w, epsil
write(16,*) w, epsil
!
enddo
!
deallocate ( xk, wk )
deallocate ( enval, hk, u, fv1, fv2, ss, vs )
!
return
end subroutine dielec_mat
!----------------------------------------------------------------
!