SUBROUTINE green_linsys (ik0) 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, RYTOEV, eps8 USE cell_base, ONLY : tpiba2 USE ener, ONLY : ef USE klist, ONLY : xk, wk, nkstot USE gvect, ONLY : nrxx, g, nl, ngm, ecutwfc USE gsmooth, ONLY : doublegrid, nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, ngms USE lsda_mod, ONLY : lsda, nspin, current_spin, isk USE wvfct, ONLY : nbnd, npw, npwx, igk, g2kin, et 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, eta, tr2_green, w_green_start USE nlcc_gw, ONLY : nlcc_any USE units_gw, ONLY : iuwfc, lrwfc, iuwfcna, iungreen, lrgrn USE eqv, ONLY : evq, eprec USE qpoint, ONLY : xq, npwq, igkq, nksq, ikks, ikqs USE recover_mod, ONLY : read_rec, write_rec USE mp, ONLY : mp_sum USE disp, ONLY : nqs USE freq_gw, ONLY : fpol, fiu, nfs, nfsmax, nwgreen, wgreen USE gwsigma, ONLY : ngmgrn, ngmsco USE mp_global, ONLY : inter_pool_comm, intra_pool_comm, mp_global_end, mpime, & nproc_pool, nproc, me_pool, my_pool_id, npool USE mp, ONLY: mp_barrier, mp_bcast, mp_sum IMPLICIT NONE real(DP) :: thresh, anorm, averlt, dr2 logical :: conv_root COMPLEX(DP) :: rhs(npwx, 1) COMPLEX(DP) :: aux1(npwx) COMPLEX(DP) :: ci, cw, green(ngmgrn,ngmgrn) COMPLEX(DP), ALLOCATABLE :: etc(:,:) REAL(DP) :: eprecloc INTEGER :: iw, igp, iwi INTEGER :: iq, ik0 INTEGER :: ngvecs, ngsol INTEGER :: rec0, n1 REAL(DP) :: dirac, delta REAL(DP) :: x real(DP) :: k0mq(3) real(DP) :: w_ryd(nwgreen) external cg_psi, cch_psi_all_fix, cch_psi_all_green real(DP) , allocatable :: h_diag (:,:) COMPLEX(DP), ALLOCATABLE :: gr_A(:,:) 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 INTEGER, PARAMETER :: lmres = 1 !Arrays to handle case where nlsco does not contain all G vectors required for |k+G| < ecut INTEGER :: igkq_ig(npwx) INTEGER :: igkq_tmp(npwx) INTEGER :: counter !PARALLEL INTEGER :: igstart, igstop, ngpool, ngr, igs !LINALG COMPLEX(DP), EXTERNAL :: zdotc REAL(DP) :: ar, ai !HL NOTA BENE: !Green's function has dimensions npwx, 1 in current notation... !Since we store in G space there is no truncation error at this stage... !When we start going to real space we want to transform on to the full correlation density grid. ALLOCATE (h_diag (npwx, 1)) ALLOCATE (etc(nbnd, nkstot)) ci = (0.0d0, 1.0d0) !Convert freq array generated in freqbins into rydbergs. w_ryd(:) = wgreen(:)/RYTOEV CALL start_clock('greenlinsys') where_rec='no_recover' if (nksq.gt.1) rewind (unit = iunigk) ngvecs = igstart - igstop + 1 ALLOCATE(gr_A(npwx, 1)) !Loop over q in the IBZ_{k} DO iq = w_green_start, nksq if (lgamma) then ikq = iq else 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 !Need a loop to find all plane waves below ecutsco !when igkq takes us outside of this sphere. !igkq_tmp is gamma centered index up to ngmsco, !igkq_ig is the linear index for looping up to npwq. !need to loop over... !This is a hack version of gk_sort... is there something wrong with the logic here... !ngmsig is just the number of G vectors, really I should be packing an array sorted up to !ecutsig ... I have literally been stuck on this for years... counter = 0 igkq_tmp(:) = 0 igkq_ig(:) = 0 do ig = 1, npwx if((igkq(ig).le.ngmgrn).and.((igkq(ig)).gt.0)) then counter = counter + 1 !index in total G grid. igkq_tmp (counter) = igkq(ig) !index for loops igkq_ig (counter) = ig endif enddo !Difference in parallelization routine. Instead of parallelizing over the usual list of G-vectors as in the straight !forward pilot implementation I need to first generate the list of igkq's within my correlation cutoff !this gives the number of vectors that requires parallelizing over. Then I split the work (up to counter) between the !nodes as with the coulomb i.e. igstart and igstop. #ifdef __PARA npool = nproc / nproc_pool write(stdout,'("npool", i4, i5)') npool, counter if (npool.gt.1) then ! number of g-vec per pool and reminder ngpool = counter / npool ngr = counter - ngpool * npool ! the remainder goes to the first ngr pools if ( my_pool_id < ngr ) ngpool = ngpool + 1 igs = ngpool * my_pool_id + 1 if ( my_pool_id >= ngr ) igs = igs + ngr ! the index of the first and the last g vec in this pool igstart = igs igstop = igs - 1 + ngpool write (stdout,'(/4x,"Max n. of G vecs in Green_linsys per pool = ",i5)') igstop-igstart+1 else #endif igstart = 1 igstop = counter #ifdef __PARA endif #endif ! Now the G-vecs up to the correlation cutoff have been divided between pools. ! Calculates beta functions (Kleinman-Bylander projectors), with ! structure factor, for all atoms, in reciprocal space call init_us_2 (npwq, igkq, xk (1, ikq), vkb) ! psi_{k+q}(r) is every ikq entry call davcio (evq, lrwfc, iuwfc, ikq, - 1) 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 WRITE(6, '(4x,"k0-q = (",3f12.7," )",10(3x,f7.3))') xk(:,ikq), et(:,ikq)*RYTOEV WRITE(600+mpime, '(4x,"k0-q = (",3f12.7," )",10(3x,f7.3))') xk(:,ikq), et(:,ikq)*RYTOEV aux1=(0.d0,0.d0) DO ig = 1, npwq aux1 (ig) = g2kin (ig) * evq (ig,4) END DO eprecloc = zdotc(npwx*npol, evq(1,4), 1, aux1(1),1) gr_A(:,:) = (0.0d0, 0.0d0) DO iw = 1, nwgreen green = DCMPLX(0.0d0, 0.0d0) h_diag(:,:) = 0.d0 do ibnd = 1, 1 !For all elements up to use standard TPA: if (w_ryd(iw).le.0.0d0) then do ig = 1, npwq !The preconditioner needs to be real and symmetric to decompose as E^T*E x = (g2kin(ig)-w_ryd(iw)) h_diag(ig,ibnd) = (27.d0+18.d0*x+12.d0*x*x+8.d0*x**3.d0) & /(27.d0+18.d0*x+12.d0*x*x+8.d0*x**3.d0+16.d0*x**4.d0) enddo else !Really choosy preconditioner. do ig = 1, npwq if(g2kin(ig).gt.w_ryd(iw)) then !Trying without... x = (g2kin(ig)-w_ryd(iw)) !maybe eprecloc is no goood... ! x = (g2kin(ig)-w_ryd(iw)) h_diag(ig,ibnd) = (27.d0+18.d0*x+12.d0*x*x+8.d0*x**3.d0) & /(27.d0+18.d0*x+12.d0*x*x+8.d0*x**3.d0+16.d0*x**4.d0) else h_diag(ig,ibnd) = 1.0d0 endif enddo endif enddo do ig = igstart, igstop !allows us to use the solution for the previous frequency for the current frequency ngsol = ngvecs - (igstop-igstart) rhs(:,:) = DCMPLX(0.0d0, 0.0d0) rhs(igkq_ig(ig), 1) = DCMPLX(-1.0d0, 0.0d0) gr_A(:,:) = DCMPLX(0.0d0, 0.0d0) lter = 0 etc(:, :) = CMPLX(0.0d0, 0.0d0, kind=DP) cw = CMPLX(w_ryd(iw), eta, kind=DP) !TR !cw = CMPLX(w_ryd(iw), -eta, kind=DP) anorm = 0.0d0 !Doing Linear System with Wavefunction cutoff (full density). call cbicgstabl(cch_psi_all_green, cg_psi, etc(1,ikq), rhs, gr_A(:,1), h_diag, & npwx, ngmsco, tr2_green, ikq, lter, conv_root, anorm, 1, npol, cw, lmres, .true.) if(.not.conv_root) write(600+mpime,'(2f15.6)') wgreen(iw), anorm !In case of breakdown... ar = real(anorm) if (( ar .ne. ar )) then gr_A(:,:) = (0.0d0,0.0d0) write(600 + mpime,'("anorm imaginary")') endif do igp = 1, counter green (igkq_tmp(ig), igkq_tmp(igp)) = green (igkq_tmp(ig), igkq_tmp(igp)) + gr_A(igkq_ig(igp),1) !HLTR !green (igkq_tmp(ig), igkq_tmp(igp)) = green (igkq_tmp(ig), igkq_tmp(igp)) + gr_A(igkq_ig(igp),1) enddo enddo !ig !Green's Fxn Non-analytic Component: !PARA do ig = igstart, igstop do igp = 1, counter !should be nbnd_occ: do ibnd = 1, nbnd x = et(ibnd, ikq) - w_ryd(iw) dirac = eta / pi / (x**2.d0 + eta**2.d0) !Green should now be indexed (igkq_tmp(ig), igkq_tmp(igp)) according to the !large G-grid which extends out to 4*ecutwfc. Keep this in mind when doing !ffts and taking matrix elements (especially matrix elements in G space!). !Normal. green(igkq_tmp(ig), igkq_tmp(igp)) = green(igkq_tmp(ig), igkq_tmp(igp)) + & tpi*ci*conjg(evq(igkq_ig(ig), ibnd)) * & evq(igkq_ig(igp), ibnd) * dirac enddo enddo enddo !Collect G vectors across processors and then write the full green's function to file. #ifdef __PARA CALL mp_barrier(inter_pool_comm) !Collect all elements of green's matrix from different !processors. CALL mp_sum (green, inter_pool_comm ) CALL mp_barrier(inter_pool_comm) if(ionode) then #endif !HL Original: ! rec0 = (iw-1) * 1 * nqs + (ik0-1) * nqs + (iq-1) + 1 rec0 = (iw-1) * 1 * nksq + (iq-1) + 1 CALL davcio(green, lrgrn, iungreen, rec0, +1, ios) #ifdef __PARA endif CALL mp_barrier(inter_pool_comm) #endif ENDDO ! iw ENDDO ! iq DEALLOCATE(h_diag) DEALLOCATE(etc) CALL stop_clock('greenlinsys') RETURN END SUBROUTINE green_linsys