SUBROUTINE green_linsys(ik0, iq) ! SGW: call green_linsys ( vr, g2kin, k0mq, nwgreen, wgreen, igstart, igstop, ik0, iq ) ! "In order to calculate the analytic component G^{A} we consider G^{A}(r,r',w) as a ! parametric function of the first space variable and of the frequency: G^{A}_[r,w](r') ! G^{A} = \sum_{n}\psi_n(r)\psi_n(r')/(w - ev + idelta) (16) ! (H - w^{+})G^{A}_[r,w] = -\delta_[r](r') (18) ! Now we will extend our analysis to an explicitly plane wave basis. ! let n -> nk so we are dealing with bloch states and expand the functions in terms of ! the plane waves: ! G_[r,w]^{A}(r') = \frac{1}{(N_{k}\Omega)}\sum_{kg}g^{A}_[k,\omega,G](\r')e^{-i(k + G)\r} e^{ik\r'} ! (H_{k} - \omega^{+})g_{[k,G,\omega]}(G') = -\delta_{GG'} USE kinds, ONLY : DP USE ions_base, ONLY : nat, ntyp => nsp, ityp USE io_global, ONLY : stdout, ionode USE io_files, ONLY : prefix, iunigk USE check_stop, ONLY : check_stop_now USE wavefunctions_module, ONLY : evc USE constants, ONLY : degspin, pi, tpi USE cell_base, ONLY : tpiba2 USE ener, ONLY : ef USE klist, ONLY : lgauss, degauss, ngauss, xk, wk, nkstot USE gvect, ONLY : nrxx, g, nl, ngm USE gsmooth, ONLY : doublegrid, nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, ngms USE lsda_mod, ONLY : lsda, nspin, current_spin, isk USE spin_orb, ONLY : domag USE wvfct, ONLY : nbnd, npw, npwx, igk, g2kin, et USE scf, ONLY : rho USE uspp, ONLY : okvan, vkb USE uspp_param, ONLY : upf, nhm, nh USE noncollin_module, ONLY : noncolin, npol, nspin_mag USE paw_variables, ONLY : okpaw USE paw_onecenter, ONLY : paw_dpotential, paw_dusymmetrize, & paw_dumqsymmetrize USE control_gw, ONLY : rec_code, niter_gw, nmix_gw, tr2_gw, & alpha_pv, lgamma, lgamma_gamma, convt, & nbnd_occ, alpha_mix, ldisp, rec_code_read, & where_rec, flmixdpot, current_iq, & ext_recover USE nlcc_gw, ONLY : nlcc_any USE units_gw, ONLY : iuwfc, lrwfc USE eqv, ONLY : evq, eprec USE qpoint, ONLY : xq, npwq, igkq, nksq, ikks, ikqs USE recover_mod, ONLY : read_rec, write_rec ! used oly to write the restart file USE mp_global, ONLY : inter_pool_comm, intra_pool_comm USE mp, ONLY : mp_sum ! USE qpoint, ONLY : xq, npwq, igkq, nksq, ikks, ikqs USE disp, ONLY : nqs USE freq_gw, ONLY : fpol, fiu, nfs, nfsmax USE gwsigma, ONLY : green USE units_gw, ONLY : iungreen implicit none real(DP) :: thresh, anorm, averlt, dr2 logical :: conv_root COMPLEX(DP) :: gr_A(ngms), gr_N(ngms), rhs(ngms), gr(ngms), ci INTEGER :: iw, igp INTEGER :: iq, ik0 INTEGER :: rec0, n1 REAL(DP) :: dirac, x, delta real(DP) :: k0mq(3) real(DP) :: eta external cch_psi_all, ch_psi_all, cg_psi real(DP) , allocatable :: h_diag (:,:) integer :: kter, & ! counter on iterations iter0, & ! starting iteration ipert, & ! counter on perturbations ibnd, & ! counter on bands iter, & ! counter on iterations lter, & ! counter on iterations of linear system ltaver, & ! average counter lintercall, & ! average number of calls to cgsolve_all ik, ikk, & ! counter on k points ikq, & ! counter on k+q points ig, & ! counter on G vectors ndim, & ! integer actual row dimension of dpsi is, & ! counter on spin polarizations nt, & ! counter on types na, & ! counter on atoms nrec, nrec1,& ! the record number for dvpsi and dpsi ios, & ! integer variable for I/O control mode ! mode index !Parameter for the dirac delta function ! real(DP), parameter :: eta = 0.04 allocate (h_diag ( npwx*npol, nbnd)) ci = (0.0d0, 1.0d0) eta = 0.04 ! write(6,*)ik0 ikq = ikqs(ik0) ! if (.not.lgamma.and.nksq.gt.1) then ! read (iunigk, err = 200, iostat = ios) npwq, igkq !200 call errore ('solve_linter', 'reading igkq', abs (ios) ) ! endif ! Should technically be looking at the wave vector k0mq = k0 - q ! for now i'm just reading in k+q as the wave vector. Some thinking ! needs to be done here re: commensurate meshes. For test cases shouldn't ! be to hard just generate the wavefunctions at xk - xq for the small test systems k0mq(:) = xk(:,ik0) - xq(:) call davcio (evq, lrwfc, iuwfc, ikq, -1) !do ibnd = 1,4 ! do ig=1,npwx ! write(6,*)evq(ig, ibnd) ! enddo !enddo !stop do ig = 1, npwq g2kin (ig) = ( (xk (1,ikq) + g (1, igkq(ig)) ) **2 + & (xk (2,ikq) + g (2, igkq(ig)) ) **2 + & (xk (3,ikq) + g (3, igkq(ig)) ) **2 ) * tpiba2 enddo h_diag = 0.d0 do ibnd = 1, nbnd_occ (ikq) do ig = 1, npwq h_diag(ig,ibnd)=1.d0/max(1.0d0,g2kin(ig)/eprec(ibnd,ikq)) enddo IF (noncolin) THEN do ig = 1, npwq h_diag(ig+npwx,ibnd)=1.d0/max(1.0d0,g2kin(ig)/eprec(ibnd,ikq)) enddo END IF enddo !write(6,*)npwq, npwx, ngm, ngms !stop do iw = 1, nfs ! do ig = 1, ngm do ig = 1, npwq !zero after this point ! Solve linear system to obtain: G^{A}_[r,w](r') rhs(:) = (0.0d0, 0.0d0) rhs(ig) = -(1.0d0, 0.0d0) gr_A = (0.0d0, 0.0d0) conv_root = .true. call cbcg_solve(cch_psi_all, cg_psi, et(:,ikq), rhs, gr_A, h_diag, & npwx, npwq, thresh, ikq, lter, conv_root, anorm, nbnd_occ(ikq), npol, .true.) write(6,*)gr_A ! Evaluate the G_{N} part of the Green's function. gr_N(:) = (0.0d0, 0.0d0) !Do ibnd = 1, nbnd_occ(k0mq) do ibnd = 1, 4 ! x = w_ryd(iw) - eval(ibnd) x = fiu(iw) - et(ibnd, ikq) dirac = eta / pi / (x**2.d0 + eta**2.d0) ! no spin factor (see above) ! gr_N = gr_N + 2.d0 * twopi * ci * conjg ( psi(ig,ibnd) ) * psi(1:ngm,ibnd) * dirac ! gr_N = gr_N + tpi * ci * conjg( evq(ig,ibnd) ) * evq(1:ngm,ibnd) * dirac do igp = 1, npwq gr_N(ig) = gr_N(ig) + tpi * ci * conjg( evq(ig,ibnd) ) * evq(igp,ibnd) * dirac enddo ! write(6,*) gr_N(ig) ! write(6,*) et(ibnd,ikq) ! write(6,*) evq(ig,ibnd) ! write(6,*) dirac enddo gr = gr_A + gr_N ! write(6,*)gr ! write(stdout,*) gr do igp = 1, ngms green (ig,igp,iw,nspin_mag) = gr (igp) enddo enddo ! Collect G vectors across processors and then write the full green's function to file. rec0 = (iw-1) * nksq * nqs + (ik0-1) * nqs + (iq-1) + 1 write ( iungreen, rec = rec0, iostat = ios) green enddo END SUBROUTINE green_linsys