!----------------------------------------------------------------------- ! Copyright (C) 2010-2015 Henry Lambert, Feliciano Giustino ! This file is distributed under the terms of the GNU General Public ! License. See the file `LICENSE' in the root directory of the ! present distribution, or http://www.gnu.org/copyleft.gpl.txt . !----------------------------------------------------------------------- SUBROUTINE coulomb(iq, igstart, igstop, scrcoul) !----------------------------------------------------------------------- ! This subroutine is the main driver of the COULOMB self consistent cycle ! which calculates the dielectric matrix by generating the density response ! to a charge dvbare(nl(ig)) = 1.00 + i*0.00 at a single fourier component (G). ! The dielectric matrix is given by: ! eps_{q}^{-1}(G,G',iw) = (\delta_{GG'} + drhoscfs^{scf}_{G,G',iw}) USE kinds, ONLY : DP USE ions_base, ONLY : nat USE constants, ONLY : e2, fpi, RYTOEV, pi, eps8 USE cell_base, ONLY : alat, tpiba2, omega USE lsda_mod, ONLY : nspin USE io_global, ONLY : stdout, ionode USE uspp, ONLY : okvan USE control_gw, ONLY : zue, convt, rec_code, modielec, eta, godbyneeds, padecont,& solve_direct, do_epsil, do_q0_only USE partial, ONLY : done_irr, comp_irr USE modes, ONLY : nirr, npert, npertx USE uspp_param, ONLY : nhm USE eqv, ONLY : drhoscfs, dvbare USE paw_variables, ONLY : okpaw USE noncollin_module, ONLY : noncolin, nspin_mag USE gwsigma, ONLY : sigma_c_st USE qpoint, ONLY : xq USE freq_gw, ONLY : fpol, fiu, nfs, nfsmax, nwcoul, wcoul USE units_gw, ONLY : iuncoul, lrcoul USE disp, ONLY : nqs, nq1, nq2, nq3 !Symmetry Stuff USE gwsymm, ONLY : ig_unique, ngmunique !FFTS USE gvect, ONLY : ngm, g, nl USE gvecs, ONLY : nls USE fft_base, ONLY : dfftp, dffts USE fft_interfaces, ONLY : invfft, fwfft USE klist, ONLY : lgauss USE mp_world, ONLY : mpime USE mp_pools, ONLY : me_pool, root_pool, inter_pool_comm USE mp, ONLY : mp_sum, mp_barrier USE mp_global, ONLY : inter_image_comm, intra_image_comm, & my_image_id, nimage, root_image IMPLICIT NONE REAL(DP) :: tcpu, get_clock ! timing variables REAL(DP) :: qg2, qg2coul INTEGER :: ig, igp, iw, npe, irr, icounter INTEGER :: igstart, igstop, igpert, isp COMPLEX(DP), allocatable :: drhoaux (:,:) COMPLEX(DP) :: padapp, w !HL temp variable for scrcoul to write to file. COMPLEX(DP) :: cw INTEGER :: unf_recl, recl, ios INTEGER :: iq, screening LOGICAL :: exst !again should decide if this should be allocated globally. COMPLEX(DP) :: scrcoul(sigma_c_st%ngmt, sigma_c_st%ngmt, nfs, 1) !modeps and spencer-alavi vars REAL(DP) :: wwp, eps0, q0, wwq, fac REAL(DP) :: qg, rcut, spal ! used to test the recover file EXTERNAL get_clock CALL start_clock ('coulomb') if(solve_direct) then ALLOCATE (drhoscfs(dffts%nnr, nfs)) else !for self-consistent solution we only consider one !frequency at a time. To save memory and time and lines of codes etc. !we use the frequency variable for multishift as the nspin_mag var. !to extend this to magnetic with multishift we need to add another !dimension to drhoscfrs WRITE(stdout, '(4x,4x,"nspinmag", i4)'), nspin_mag ALLOCATE (drhoscfs(dffts%nnr, nspin_mag)) endif irr=1 scrcoul(:,:,:,:) = (0.d0, 0.0d0) !LOOP OVER ig, unique g vectors only. !g is sorted in magnitude order. !WRITE(1000+mpime, '(2i4)') igstart, igstop DO ig = igstart, igstop ! if (do_q0_only.and.ig.gt.1) CYCLE qg2 = (g(1,ig_unique(ig))+xq(1))**2+(g(2,ig_unique(ig))+xq(2))**2+(g(3,ig_unique(ig))+xq(3))**2 IF(solve_direct) THEN !if(qg2.lt.eps8) CYCLE drhoscfs(:,:) = dcmplx(0.0d0, 0.0d0) dvbare(:) = dcmplx(0.0d0, 0.0d0) dvbare (nls(ig_unique(ig)) ) = dcmplx(1.d0, 0.d0) CALL invfft('Smooth', dvbare, dffts) CALL solve_lindir (dvbare, drhoscfs) CALL fwfft('Smooth', dvbare, dffts) do iw = 1, nfs CALL fwfft ('Dense', drhoscfs(:,iw), dffts) !if(ig.eq.1.or.mod(ig,10).eq.0) WRITE(stdout, '(4x,4x,"eps_{GG}(q,w) = ", 2f10.4)'),drhoscfs(nls(ig_unique(ig)),iw)+dvbare(nls(ig_unique(ig))) WRITE(stdout, '(4x,4x,"eps_{GG}(q,w) = ", 2f10.4)'), drhoscfs(nls(ig_unique(ig)),iw)+dvbare(nls(ig_unique(ig))) do igp = 1, sigma_c_st%ngmt if(igp.ne.ig_unique(ig)) then !diagonal elements drho(G,G'). scrcoul(ig_unique(ig), igp, iw, nspin_mag) = drhoscfs(nls(igp), iw) else !diagonal elements eps(\G,\G') = \delta(G,G') - drho(G,G'). scrcoul(ig_unique(ig), igp, iw, nspin_mag) = drhoscfs(nls(igp), iw) + dvbare(nls(ig_unique(ig))) endif enddo enddo !iw ELSE if(qg2.lt.0.001.AND.lgauss) then write(6,'("Not calculating static electric field applied to metal, cycling coulomb")') WRITE(stdout, '(4x,4x,"inveps_{GG}(q,w) = 0.000000 0.0000000")') CYCLE endif DO iw = 1, nfs drhoscfs(:,:) = dcmplx(0.0d0, 0.0d0) dvbare(:) = dcmplx(0.0d0, 0.0d0) dvbare (nls(ig_unique(ig)) ) = dcmplx(1.d0, 0.d0) CALL invfft('Smooth', dvbare, dffts) CALL solve_linter (dvbare, iw, drhoscfs) CALL fwfft('Smooth', dvbare, dffts) DO isp =1 , nspin_mag CALL fwfft('Dense', drhoscfs(:,isp), dffts) ENDDO IF(ionode) THEN WRITE(stdout, '(4x,4x,"inveps_{GG}(q,w) = ", 2f16.9)'), drhoscfs(nls(ig_unique(ig)), 1) + dvbare(nls(ig_unique(ig))) DO isp = 1, nspin_mag DO igp = 1, sigma_c_st%ngmt scrcoul(ig_unique(ig), igp, iw, isp) = drhoscfs(nl(igp), isp) ENDDO ENDDO ENDIF ENDDO ENDIF !solve_direct/solve_linter ENDDO 545 CONTINUE tcpu = get_clock ('GW') DEALLOCATE (drhoscfs) CALL stop_clock ('coulomb') RETURN END SUBROUTINE coulomb