! !---------------------------------------------------------------- subroutine dielec_mat ( xxq ) !---------------------------------------------------------------- ! ! 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 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 (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 !---------------------------------------------------------------- !