SUBROUTINE sigma_extra(ik0) !Sigma Extra !now could add extra ananlytic piece: -\inf to wsigmamin-wcoul !and from wsigmamin+wcoul to +\inf !G(r,r';w'-w0) = FFT[\delta(\G,G')/(w-w0-w1)] !W(r,r';w') !\Sigma(r, r'; w0) = \Sigma(r,r'; w0) + \int_{\-inf}^{wcoul) G(r,r'; w0 - w')W(r,\r';w')dw' ! + \int_{wcoul}^{\inf) G(r,r'; w0 - w')W(r,\r';w')dw' ! USE kinds, ONLY : DP USE klist, ONLY : wk, xk, nkstot, nks USE constants, ONLY : e2, fpi, RYTOEV, tpi, eps8, pi USE disp, ONLY : nqs, nq1, nq2, nq3, wq, x_q, xk_kpoints, num_k_pts USE cell_base, ONLY : omega, tpiba2 USE units_gw, ONLY : iuncoul, iungreen, iunsigma, lrsigma, lrcoul, lrgrn, iuwfc, lrwfc,iunsigext,lrsigext USE gwsigma, ONLY : ngmsco, nrsco, nlsco, fft6_g2r, ecutsco, ngmsig,& nr1sco, nr2sco, nr3sco, ngmgrn, ngmpol USE eqv, ONLY : evq, eprec USE qpoint, ONLY : xq, npwq, igkq, nksq, ikks, ikqs USE wvfct, ONLY : nbnd, npw, npwx, igk, g2kin, et USE io_files, ONLY : prefix, iunigk USE fft_scalar, ONLY : cfft3d USE io_global, ONLY : stdout, ionode_id, ionode USE freq_gw, ONLY : fpol, fiu, nfs, nfsmax,& nwcoul, nwgreen, nwalloc, nwsigma, wtmp, wcoul,& wgreen, wsigma, wsigmamin, wsigmamax,& deltaw, wcoulmax, ind_w0mw, ind_w0pw, plasmon USE control_gw, ONLY : lgamma, eta USE gvect, ONLY : ngm, nrxx, g, nr1, nr2, nr3, nrx1, nrx2, nrx3, nl IMPLICIT NONE INTEGER :: rec0, ios INTEGER :: ikq !Integration Variable: COMPLEX(DP) :: cprefac !V arrays: COMPLEX(DP), ALLOCATABLE :: scrcoul(:,:) INTEGER, ALLOCATABLE :: gmapsym(:,:) COMPLEX(DP), ALLOCATABLE :: barcoul(:,:), barcoulr(:,:) REAL(DP) :: rcut, spal REAL(DP) :: qg2, xq0s(3), qg, xxq(3) !G arrays: COMPLEX(DP), ALLOCATABLE :: green_id (:,:), green_idr(:,:) COMPLEX(DP), ALLOCATABLE :: sigma_g(:,:), sigma(:,:) !Constants: COMPLEX(DP) :: cone COMPLEX(DP) :: ci, czero !Counters: INTEGER :: ig, igp, npe, irr, icounter, ir, irp INTEGER :: iq, ipol, ibnd, jbnd, counter INTEGER :: ikmq, ik0, ik INTEGER :: igkq_ig(npwx) INTEGER :: igkq_tmp(npwx) COMPLEX(DP) :: aux(nrsco) COMPLEX(DP) :: aux1(nrsco) REAL(DP) :: xq_coul(3) REAL(DP) :: dirac, delta, x ALLOCATE ( barcoul (ngmsco, ngmsco) ) ALLOCATE ( barcoulr (nrsco, nrsco) ) ALLOCATE ( green_id (ngmgrn, ngmgrn) ) ALLOCATE ( green_idr (nrsco, nrsco) ) if(.not.allocated(sigma_g)) ALLOCATE(sigma_g (ngmsco, ngmsco)) ALLOCATE ( sigma (nrsco, nrsco) ) cone = DCMPLX(1.0d0,0.0d0) ci = (0.0d0, 1.d0) czero = (0.0d0, 0.0d0) sigma(:,:) = (0.0d0, 0.0d0) !Need a loop to find all plane waves below ecutsco when igkq takes us outside of this sphere. counter = 0 igkq_tmp = 0 igkq_ig = 0 if (nksq.gt.1) rewind (unit = iunigk) wq(:) = 0.0d0 do iq = 1, nksq if(lgamma) then wq(iq) = 0.5d0*wk(iq) else wq(iq) = 0.5d0*wk(2*iq-1) endif enddo if (lgamma) then ikq = ik0 else ikq = 2*ik0 endif DO iq = 1, nksq if (lgamma) then ikq = iq else ! even k + q ikq = 2*iq endif if (nksq.gt.1) then read (iunigk, err = 100, iostat = ios) npw, igk 100 CALL errore ('green_linsys', 'reading igk', abs (ios) ) endif if (lgamma) npwq = npw if (.not.lgamma.and.nksq.gt.1) then read (iunigk, err = 200, iostat = ios) npwq, igkq 200 call errore ('green_linsys', 'reading igkq', abs (ios) ) endif !Should have a data_type which keeps track of these indices... counter = 0 igkq_tmp(:) = 0 igkq_ig(:) = 0 do ig = 1, npwq if((igkq(ig).le.ngmsco).and.((igkq(ig)).gt.0)) then counter = counter + 1 igkq_tmp (counter) = igkq(ig) igkq_ig (counter) = ig endif enddo !Green \delta(G,G')= green_id(:,:) = (0.0d0, 0.0d0) do ig = 1, counter green_id(igkq_tmp(ig),igkq_tmp(ig)) = (1.0d0, 0.0d0) enddo ! WRITE(6, '("FFT_GREEN")') green_idr(:,:) = czero do ig = 1, ngmgrn aux(:) = czero do igp = 1, ngmgrn aux(nlsco(igp)) = green_id(ig,igp) enddo call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1) do irp = 1, nrsco green_idr(ig, irp) = aux(irp) / omega enddo enddo ! the conjg/conjg is to calculate sum_G f(G) exp(-iGr) ! following the convention set in the paper ! [because the standard transform is sum_G f(G) exp(iGr) ] do irp = 1, nrsco aux = czero do ig = 1, ngmsco aux(nlsco(ig)) = conjg( green_idr(ig,irp) ) enddo call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1) green_idr(1:nrsco,irp) = conjg ( aux ) enddo !Now have G^A(\r,\r')... !Calculate v(r,r'): rcut = (float(3)/float(4)/pi*omega*float(nq1*nq2*nq3))**(float(1)/float(3)) ! -q xq0s = (/ -0.01 , 0.00, 0.00 /) ! this should be set from input barcoul(:,:) = (0.0d0,0.0d0) ! -q = k0 - (k0 + q) xq_coul(:) = xk_kpoints(:,ik0) - xk(:,ikq) do ig = 1, ngmgrn qg = sqrt((g(1,ig) + xq_coul(1))**2.d0 + (g(2,ig) + xq_coul(2))**2.d0 & + (g(3,ig ) + xq_coul(3))**2.d0) qg2 = (g(1,ig) + xq_coul(1))**2.d0 + (g(2,ig) + xq_coul(2))**2.d0 & + ((g(3,ig)) + xq_coul(3))**2.d0 ! These if conditions need to go... if (qg < eps8) qg = sqrt((g(1, ig) + xq0s(1))**2.d0 + (g(2, ig) + xq0s(2))**2.d0 & + (g(3, ig) + xq0s(3))**2.d0) if (qg2 < eps8) qg2 = ((g(1, ig) + xq0s(1))**2.d0 + (g(2, ig) + xq0s(2))**2.d0 & + (g(3, ig) + xq0s(3))**2.d0) spal = 1.0d0 - cos (rcut * sqrt(tpiba2) * qg) barcoul (ig, ig) = e2 * fpi / (tpiba2*qg2) * dcmplx(spal, 0.0d0) enddo ! WRITE(6, '("FFT_COULOMB")') barcoulr(:,:) = (0.0d0, 0.0d0) do ig = 1, ngmgrn aux(:) = czero do igp = 1, ngmgrn aux(nlsco(igp)) = barcoul(ig,igp) enddo call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1) do irp = 1, nrsco barcoulr(ig, irp) = aux(irp) / omega enddo enddo ! the conjg/conjg is to calculate sum_G f(G) exp(-iGr) ! following the convention set in the paper ! [because the standard transform is sum_G f(G) exp(iGr) ] do irp = 1, nrsco aux = czero do ig = 1, ngmgrn aux(nlsco(ig)) = conjg( barcoulr(ig,irp) ) enddo call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1) barcoulr(1:nrsco,irp) = conjg ( aux ) enddo !Calculate \Sigma_extra(\r,\r') = \frac{i}{2\pi}(\frac{\omega_p}{\omega_c})^2 G(\r,\r')v(\r,\r') !cprefac = (DCMPLX(plasmon, 0.0d0)/DCMPLX((wcoulmax), 0.0d0))**2 & ! * DCMPLX(0.0d0,1.0d0)/tpi cprefac = (DCMPLX(plasmon, 0.0d0)/DCMPLX((wcoulmax), 0.0d0))**2/tpi !sigma = sigma + cprefac * green_idr * barcoulr !sigma = sigma + cprefac * green_idr sigma(:,:) = sigma(:,:) + cprefac*wq(iq)*green_idr(:,:)*barcoulr(:,:) ENDDO ! on q sigma_g = (0.0d0,0.0d0) do ir = 1, nrsco aux = (0.0d0, 0.0d0) do irp = 1, nrsco aux(irp) = sigma(ir,irp) enddo call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, -1) do igp = 1, ngmsco sigma (ir, igp) = aux(nlsco(igp)) enddo enddo do igp = 1, ngmsco aux = czero do ir = 1, nrsco aux(ir) = conjg ( sigma(ir,igp) ) enddo call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, -1) do ig = 1, ngmsco sigma (ig,igp) = conjg ( aux( nlsco( ig )) ) * omega enddo enddo do ig = ngmsco + 1, nrsco do igp = ngmsco + 1, nrsco sigma (ig, igp) = (0.0d0, 0.0d0) enddo enddo !sigma_g = sigma(1:ngmsco,1:ngmsco,:) do ig = 1, ngmsco do igp = 1, ngmsco sigma_g(ig,igp) = sigma(ig,igp) !if(ig.eq.igp) sigma_g(ig,igp) = ((plasmon)/(wcoulmax))**2 enddo enddo CALL davcio (sigma_g, lrsigext, iunsigext, 1, 1) DEALLOCATE ( sigma ) DEALLOCATE ( sigma_g ) RETURN END SUBROUTINE sigma_extra