SUBROUTINE green_linsys(ik0, iq) 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 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 USE nlcc_gw, ONLY : nlcc_any USE units_gw, ONLY : iuwfc, lrwfc, iuwfcna, iungreen USE eqv, ONLY : evq, eprec USE qpoint, ONLY : xq, npwq, igkq, nksq, ikks, ikqs USE recover_mod, ONLY : read_rec, write_rec USE mp_global, ONLY : inter_pool_comm, intra_pool_comm USE mp, ONLY : mp_sum USE disp, ONLY : nqs USE freq_gw, ONLY : fpol, fiu, nfs, nfsmax, nwgreen, wgreen USE gwsigma, ONLY : ngmsig, ecutsco, ngmsco IMPLICIT NONE real(DP) :: thresh, anorm, averlt, dr2 logical :: conv_root COMPLEX(DP) :: gr_A(ngm, 1), rhs(ngm , 1) COMPLEX(DP) :: gr_N(ngm, 1), gr(ngm, 1), ci, cw, green(ngmsig,ngmsig) COMPLEX(DP), ALLOCATABLE :: etc(:,:) INTEGER :: iw, igp, iwi INTEGER :: iq, ik0 INTEGER :: rec0, n1 REAL(DP) :: dirac, x, delta real(DP) :: k0mq(3) real(DP) :: w_ryd(nwgreen) external cch_psi_all, ch_psi_all, cg_psi, cch_psi_all_fix, cch_psi_all_green 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 !HL need a threshold here for the linear system solver. This could also go in the punch card !with some default at a later date. REAL(DP) :: tr_cgsolve = 1.0d-5 !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 allocate (h_diag ( ngm, nbnd)) allocate (etc(nbnd, nkstot)) !Convert freq array generated in freqbins into rydbergs. ci = (0.0d0, 1.0d0) w_ryd(:) = wgreen(:)/RYTOEV CALL start_clock('greenlinsys') where_rec='no_recover' !HL for q = Gamma the k+q list coincides with the klist if (lgamma) write(6, '("lgamma=.true.")') if (nksq.gt.1) rewind (unit = iunigk) if (nksq.gt.1) then read (iunigk, err = 100, iostat = ios) npw, igk 100 call errore ('green_linsys', 'reading igk', abs (ios) ) endif 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. counter = 0 igkq_tmp = 0 igkq_ig = 0 do ig = 1, npwx 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 if (lgamma) then ikq = 1 else ikq = 2 endif ! 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) 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 h_diag = 1.d0 ibnd = 1 do ig = 1, npwq x = g2kin(ig)/(eprec(2,1)) 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 WRITE( 6, '(4x,"k0-q = (",3f12.7," )",10(3x,f7.3))') xq, et(:,ikq)*RYTOEV DO iw = 1, nwgreen green =(0.0d0, 0.0d0) !PARA case !do ig = igstart, igstop do ig = 1, counter !zero after this point rhs(:,:) = (0.0d0, 0.0d0) !rhs(ig, 1) = -(1.0d0, 0.0d0) !should ensure we calculate all G-vectors within sco cutoff. rhs(igkq_ig(ig), 1) = -(1.0d0, 0.0d0) gr_A(:,:) = (0.0d0, 0.0d0) lter = 0 etc(:, :) = CMPLX( 0.0d0, 0.0d0, kind=DP) cw = CMPLX( w_ryd(iw), eta, kind=DP) !Doing Linear System with Wavefunction cutoff (and full density). call cbcg_solve_fix(cch_psi_all_green, cg_psi, etc(1,ikq), rhs, gr_A, h_diag, & npwx, npwq, tr_cgsolve, ikq, lter, conv_root, anorm, 1, npol, cw, .true.) ! IF (.not.conv_root) then ! WRITE(6,'(" Root Not Converged ")') ! ENDIF gr = gr_A do igp = 1, counter green (igkq_tmp(ig), igkq_tmp(igp)) = green (igkq_tmp(ig), igkq_tmp(igp)) + gr(igkq_ig(igp),1) enddo enddo !ig !NON-ANALYTIC do ig = 1, counter gr_N(:,1) = (0.0d0, 0.0d0) do igp = 1, counter do ibnd = 1, 4 x = w_ryd(iw) - et(ibnd, ikq) dirac = eta / pi / (x**2.d0 + eta**2.d0) green(igkq_tmp(ig), igkq_tmp(igp)) = green(igkq_tmp(ig), igkq_tmp(igp)) + & tpi*ci*(evq(igkq_ig(ig), ibnd)) * & conjg(evq(igkq_ig(igp), ibnd)) * dirac enddo ! if((igkq_tmp(ig).eq.1).and.(igkq_tmp(igp).eq.1)) write(100, '(3f15.10)') wgreen(iw), green(igkq_tmp(ig), igkq_tmp(igp)) ! if((igkq_tmp(ig).eq.7).and.(igkq_tmp(igp).eq.7)) write(101,'(3f15.10)') wgreen(iw), green(igkq_tmp(ig), igkq_tmp(igp)) ! if((igkq_tmp(ig).eq.7).and.(igkq_tmp(igp).eq.11)) write(102,'(3f15.10)') wgreen(iw), green(igkq_tmp(ig), igkq_tmp(igp)) enddo enddo ! Collect G vectors across processors and then write the full green's function to file. ! #ifdef __PARA ! use poolreduce to bring together the results from each pool ! call poolreduce ( 2 * ngms * ngms, green) ! if (me.eq.1.and.mypool.eq.1) then ! #endif rec0 = (iw-1) * 1 * nqs + (ik0-1) * nqs + (iq-1) + 1 write ( iungreen, rec = rec0, iostat = ios) green !#ifdef __PARA ! endif !#endif ENDDO !iw CALL stop_clock('greenlinsys') RETURN END SUBROUTINE green_linsys