TABLE OF CONTENTS
ABINIT/cchi0q0 [ Functions ]
NAME
cchi0q0
FUNCTION
Calculate chi0 in the limit q-->0
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (GMR, VO, LR, RWG, MG, RShaltaf) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
use_tr=If .TRUE. Wfs_val are allocate and only resonant transitions are evaluated (assumes time reversal symmetry) Dtset <type(dataset_type)>=all input variables in this dataset Ep= datatype gathering differening parameters related to the calculation of the inverse dielectric matrix Gsph_epsG0<gvectors_data_type>: Info on the G-sphere used to describe chi0/espilon (including umklapp) %ng=number of G vectors %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G in the array gvec %phmGt(ng,nsym)=phase factor e^{-iG.\tau} needed to symmetrize oscillator matrix elements and chi0 %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$). %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1) Ep%inclvkb=flag to include (or not) the grad of Vkb Ltg_q= little group datatype nbvw=number of bands in the arrays wfrv,wfgv Kmesh<kmesh_t> The k-point mesh %kbz(3,nbz)=k-point coordinates, full Brillouin zone %tab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding irreducible point (kIBZ), where kBZ= (IS) kIBZ and I is either the inversion or the identity %tabi(nbzx)= for each point in the BZ defines whether inversion has to be considered in the relation kBZ=(IS) kIBZ (1 => only S; -1 => -S) %tabo(nbzx)= the symmetry operation S that takes kIBZ to each kBZ %tabp(nbzx)= phase factor associated to tnons e^{-i 2 \pi k\cdot R{^-1}t} ktabr(nfftot_gw,Kmesh%nbz) index of R^-(r-t) in the FFT array, where k_BZ = (IS) k_IBZ and S = \transpose R^{-1} Ep%nbnds=number of bands ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths. Ep%nomega=number of frequencies Cryst<crystal_t>= data type gathering info on symmetries and unit cell %natom=number of atoms %nsym=number of symmetry operations %symrec(3,3,nsym)=symmetry operations in reciprocal space %typat(natom)=type of each atom %xred(3,natom)=reduced coordinated of atoms %rprimd(3,3)=dimensional primitive translations in real space (bohr) %timrev=2 if time-reversal symmetry can be used, 1 otherwise Ep%npwe=number of planewaves for sigma exchange (input variable) nfftot_gw=Total number of points in the GW FFT grid Ep%omega(Ep%nomega)=frequencies Psps <type(pseudopotential_type)>=variables related to pseudopotentials %mpsang=1+maximum angular momentum for nonlocal pseudopotential Pawang<pawang_type> angular mesh discretization and related data: Pawrad(ntypat*usepaw)<Pawrad_type>=paw radial mesh and related data Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels QP_BSt<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities) %mband=MAX number of bands over k-points and spin (==Ep%nbnds) %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes KS_BSt<ebands_t>=KS energies and occupations. %eig(mband,nkpt,nsppol)=KS energies Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave. Wfd<wfd_t>=Object used to access the wavefunctions
OUTPUT
chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix for wavevector qq, and frequencies defined by Ep%omega chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)= Lower wings chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)= Upper wings chi0_head(3,3,Ep%nomega)=Head of chi0.
NOTES
*) The terms "head", "wings" and "body" of chi(G,Gp) refer to
G=Gp=0, either G or Gp=0, and neither=0 respectively
*) Symmetry conventions:
1) symmetry in real space is defined as: R_t f(r) = f(R^-1(r-t))
2) S=\transpose R^-1
3) kbz=S kibz
The wavefunctions for the k-point in the BZ are (assuming nondegenerate states):
u(G,b, Sk) = u ( S^-1G,b,k)* e^{-i(Sk+G)*t)
u(G,b,-Sk) = u*(-S^-1G,b,k)* e^{ i(Sk-G)*t)
u(r,b, Sk) = u (R^-1(r-t),b,k) e^{-iSk*t}
u(r,b,-Sk) = u*(R^-1(r-t),b,k) e^{ iSK*t}
The gradient of Vnl(K,Kp) for the k-point in the BZ should be:
gradvnl(SG,SGp,Sk)=S gradvnl(G,Gp,kibz)
/***********************************************************************/
TODO
Check npwepG0 before Switching on umklapp
PARENTS
screening
CHILDREN
accumulate_chi0_q0,accumulate_chi0sumrule,accumulate_sfchi0_q0 approxdelta,calc_wfwfg,chi0_bbp_mask,completechi0_deltapart,cwtime flush_unit,get_bz_item,get_gftt,gsph_fft_tabs,gsph_free,gsph_in_fftbox hilbert_transform,hilbert_transform_headwings,littlegroup_print make_transitions,paw_cross_ihr_comm,paw_cross_rho_tw_g,paw_rho_tw_g paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawhur_free pawhur_init,pawpwij_free,pawpwij_init,print_gsphere,read_plowannier rho_tw_g,setup_spectral,symmetrize_afm_chi0,vkbr_free,vkbr_init wfd_change_ngfft,wfd_distribute_bbp,wfd_get_cprj,wfd_get_ur wfd_paw_get_aeur,wrtout,xmpi_sum
SOURCE
114 #if defined HAVE_CONFIG_H 115 #include "config.h" 116 #endif 117 118 #include "abi_common.h" 119 120 subroutine cchi0q0(use_tr,Dtset,Cryst,Ep,Psps,Kmesh,QP_BSt,KS_BSt,Gsph_epsG0,& 121 & Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,Pawfgrtab,Paw_onsite,ktabr,ktabrf,nbvw,ngfft_gw,& 122 & nfftot_gw,ngfftf,nfftf_tot,chi0,chi0_head,chi0_lwing,chi0_uwing,Ltg_q,chi0_sumrule,Wfd,Wfdf) 123 124 use defs_basis 125 use defs_datatypes 126 use defs_abitypes 127 use m_profiling_abi 128 use m_xmpi 129 use m_errors 130 use m_time 131 132 use m_gwdefs, only : GW_TOL_DOCC, GW_TOL_W0, czero_gw, em1params_t, g0g0w 133 use m_numeric_tools, only : imin_loc 134 use m_geometry, only : normv, vdotw 135 use m_crystal, only : crystal_t 136 use m_fft_mesh, only : get_gftt 137 use m_bz_mesh, only : kmesh_t, get_BZ_item, littlegroup_t, littlegroup_print 138 use m_gsphere, only : gsphere_t, gsph_fft_tabs, gsph_in_fftbox, gsph_free, print_gsphere 139 use m_io_tools, only : flush_unit 140 use m_wfd, only : wfd_get_ur, wfd_t, wfd_distribute_bbp, wfd_get_cprj, & 141 & wfd_barrier, wfd_change_ngfft,wfd_paw_get_aeur, wfd_sym_ur 142 use m_oscillators, only : rho_tw_g, calc_wfwfg 143 use m_vkbr, only : vkbr_t, vkbr_free, vkbr_init, nc_ihr_comm 144 use m_chi0, only : hilbert_transform, setup_spectral, symmetrize_afm_chi0, approxdelta, & 145 accumulate_chi0_q0, accumulate_sfchi0_q0, hilbert_transform_headwings, & 146 completechi0_deltapart, accumulate_chi0sumrule, make_transitions, chi0_bbp_mask 147 use m_pawang, only : pawang_type 148 use m_pawrad, only : pawrad_type 149 use m_pawtab, only : pawtab_type 150 use m_paw_ij, only : paw_ij_type 151 use m_pawfgrtab, only : pawfgrtab_type 152 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy 153 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g 154 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t 155 use m_pawhr, only : pawhur_t, pawhur_free, pawhur_init, paw_ihr, paw_cross_ihr_comm 156 157 !This section has been created automatically by the script Abilint (TD). 158 !Do not modify the following lines by hand. 159 #undef ABI_FUNC 160 #define ABI_FUNC 'cchi0q0' 161 use interfaces_14_hidewrite 162 use interfaces_65_paw 163 use interfaces_70_gw, except_this_one => cchi0q0 164 !End of the abilint section 165 166 implicit none 167 168 !Arguments ------------------------------------ 169 !scalars 170 integer,intent(in) :: nbvw,nfftot_gw,nfftf_tot 171 logical,intent(in) :: use_tr 172 type(ebands_t),target,intent(in) :: QP_BSt,KS_BSt 173 type(crystal_t),intent(in) :: Cryst 174 type(Dataset_type),intent(in) :: Dtset 175 type(littlegroup_t),intent(in) :: Ltg_q 176 type(em1params_t),intent(in) :: Ep 177 type(kmesh_t),intent(in) :: Kmesh 178 type(gsphere_t),intent(in) :: Gsph_epsG0 179 type(Pseudopotential_type),intent(in) :: Psps 180 type(Pawang_type),intent(in) :: Pawang 181 type(wfd_t),target,intent(inout) :: Wfd,Wfdf 182 !arrays 183 integer,intent(in) :: ktabr(nfftot_gw,Kmesh%nbz),ktabrf(nfftf_tot*Dtset%pawcross,Kmesh%nbz) 184 integer,intent(in) :: ngfft_gw(18),ngfftf(18) 185 real(dp),intent(out) :: chi0_sumrule(Ep%npwe) 186 complex(gwpc),intent(out) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega) 187 complex(dpc),intent(out) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3) 188 complex(dpc),intent(out) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3) 189 complex(dpc),intent(out) :: chi0_head(3,3,Ep%nomega) 190 type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw) 191 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 192 type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw) 193 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 194 type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw) 195 type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom) 196 197 !Local variables ------------------------------ 198 !scalars 199 integer,parameter :: tim_fourdp=1,enough=10,two_poles=2,one_pole=1,ndat1=1 200 integer :: bandinf,bandsup,lcor,nspinor,npw_k,istwf_k,mband,nfft 201 integer :: band1,band2,iat,ig,itim_k,ik_bz,ik_ibz,io,iqlwl,ispinor1,ispinor2,isym_k 202 integer :: itypatcor,m1,m2,nkpt_summed,dim_rtwg,use_padfft,gw_fftalga,use_padfftf,mgfftf 203 integer :: my_nbbp,my_nbbpks,spin,nsppol !ig1,ig2, 204 integer :: comm,ierr,my_wl,my_wr,iomegal,iomegar,gw_mgfft,dummy 205 real(dp) :: cpu_time,wall_time,gflops 206 real(dp) :: fac,fac1,fac2,fac3,fac4,spin_fact,deltaf_b1b2,weight,factor 207 real(dp) :: max_rest,min_rest,my_max_rest,my_min_rest 208 real(dp) :: en_high,deltaeGW_enhigh_b2 209 real(dp) :: wl,wr,numerator,deltaeGW_b1b2 210 real(dp) :: gw_gsq,memreq 211 complex(dpc) :: deltaeKS_b1b2 212 logical :: qzero,luwindow 213 character(len=500) :: msg_tmp,msg,allup 214 type(gsphere_t) :: Gsph_FFT 215 !arrays 216 integer,ABI_CONTIGUOUS pointer :: kg_k(:,:) 217 integer :: ucrpa_bands(2) 218 integer :: wtk_ltg(Kmesh%nbz) 219 integer :: got(Wfd%nproc) 220 integer,allocatable :: tabr_k(:),tabrf_k(:) 221 integer,allocatable :: igffteps0(:),gspfft_igfft(:),igfftepsG0f(:) 222 integer,allocatable :: gw_gfft(:,:),gw_gbound(:,:),dummy_gbound(:,:),gboundf(:,:) 223 integer,allocatable :: bbp_ks_distrb(:,:,:,:) 224 real(dp) :: kbz(3),spinrot_kbz(4),q0(3) 225 real(dp),ABI_CONTIGUOUS pointer :: ks_energy(:,:,:),qp_energy(:,:,:),qp_occ(:,:,:) 226 real(dp),allocatable :: omegasf(:) 227 complex(gwpc) :: rhotwx(3,Wfd%nspinor**2) 228 complex(gwpc),allocatable :: rhotwg(:) 229 complex(dpc),allocatable :: green_w(:),green_enhigh_w(:) 230 complex(dpc),allocatable :: sf_lwing(:,:,:),sf_uwing(:,:,:),sf_head(:,:,:) 231 complex(dpc) :: wng(3),chq(3) 232 complex(dpc) :: ph_mkt 233 complex(gwpc),allocatable :: sf_chi0(:,:,:) 234 complex(dpc),allocatable :: kkweight(:,:) 235 complex(gwpc),allocatable :: ur1_kibz(:),ur2_kibz(:) 236 complex(gwpc),allocatable :: usr1_k(:),ur2_k(:) 237 complex(gwpc),allocatable :: wfwfg(:) 238 complex(gwpc),allocatable :: ur_ae1(:),ur_ae_onsite1(:),ur_ps_onsite1(:) 239 complex(gwpc),allocatable :: ur_ae2(:),ur_ae_onsite2(:),ur_ps_onsite2(:) 240 complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:),ug2(:) 241 complex(dpc), allocatable :: coeffW_BZ(:,:,:,:,:,:) 242 logical :: gradk_not_done(Kmesh%nibz) 243 logical,allocatable :: bbp_mask(:,:) 244 type(pawcprj_type),allocatable :: Cprj1_bz(:,:),Cprj2_bz(:,:) 245 type(pawcprj_type),allocatable :: Cprj1_ibz(:,:),Cprj2_ibz(:,:) 246 type(pawpwij_t),allocatable :: Pwij(:),Pwij_fft(:) 247 type(pawhur_t),allocatable :: Hur(:) 248 type(vkbr_t),allocatable :: vkbr(:) 249 !************************************************************************ 250 251 DBG_ENTER("COLL") 252 call cwtime(cpu_time,wall_time,gflops,"start") 253 254 ! Change FFT mesh if needed 255 if (ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3))) call wfd_change_ngfft(Wfd,Cryst,Psps,ngfft_gw) 256 257 gw_mgfft = MAXVAL(ngfft_gw(1:3)); gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10) 258 if (Dtset%pawcross==1) mgfftf = MAXVAL(ngfftf(1:3)) 259 260 ! Copy important variables. 261 comm = Wfd%comm; nsppol = Wfd%nsppol; nspinor = Wfd%nspinor; mband = Wfd%mband 262 nfft = Wfd%nfft 263 ABI_CHECK(Wfd%nfftot==nfftot_gw,"Wrong nfftot_gw") 264 dim_rtwg=1 !; if (nspinor==2) dim_rtwg=2 ! Can reduce size depending on Ep%nI and Ep%nj 265 266 ucrpa_bands(1)=dtset%ucrpa_bands(1) 267 ucrpa_bands(2)=dtset%ucrpa_bands(2) 268 luwindow=.false. 269 if(abs(dtset%ucrpa_window(1)+1_dp)>tol8.and.(abs(dtset%ucrpa_window(2)+1_dp)>tol8)) then 270 luwindow=.true. 271 endif 272 273 ! For cRPA calculation of U: read forlb.ovlp 274 if(dtset%ucrpa>=1) then 275 call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor,Kmesh,lcor,luwindow,& 276 & nspinor,nsppol,pawang,dtset%prtvol,ucrpa_bands) 277 endif 278 279 ks_energy => KS_BSt%eig 280 qp_energy => QP_BSt%eig; qp_occ => QP_BSt%occ 281 282 chi0_lwing = czero; chi0_uwing= czero; chi0_head = czero 283 284 if (Psps%usepaw==0) then 285 if (Ep%inclvkb/=0) then 286 ! Include the term <n,k|[Vnl,iqr]|n"k>' for q->0. 287 ABI_CHECK(nspinor==1,"nspinor+inclvkb not coded") 288 else 289 MSG_WARNING('Neglecting <n,k|[Vnl,iqr]|m,k>') 290 end if 291 else 292 ! For PAW+LDA+U, precalculate <\phi_i|[Hu,r]|phi_j\> 293 ABI_DT_MALLOC(HUr,(Cryst%natom)) 294 if (Dtset%usepawu/=0) then 295 call pawhur_init(hur,nsppol,Dtset%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij) 296 end if 297 end if 298 299 ! Initialize the completeness correction. 300 ABI_MALLOC(green_enhigh_w,(Ep%nomega)) 301 green_enhigh_w=czero 302 303 if (Ep%gwcomp==1) then 304 en_high=MAXVAL(qp_energy(Ep%nbnds,:,:))+Ep%gwencomp 305 write(msg,'(a,f8.2,a)')' Using completeness correction with energy ',en_high*Ha_eV,' [eV] ' 306 call wrtout(std_out,msg,'COLL') 307 ABI_MALLOC(wfwfg,(nfft*nspinor**2)) 308 309 ! Init the largest G-sphere contained in the FFT box for the wavefunctions. 310 call gsph_in_fftbox(Gsph_FFT,Cryst,Wfd%ngfft) 311 call print_gsphere(Gsph_FFT,unit=std_out,prtvol=10) 312 313 ABI_MALLOC(gspfft_igfft,(Gsph_FFT%ng)) 314 ABI_MALLOC(dummy_gbound,(2*gw_mgfft+8,2)) 315 316 ! Mapping between G-sphere and FFT box. 317 call gsph_fft_tabs(Gsph_FFT, [0, 0, 0],Wfd%mgfft,Wfd%ngfft,dummy,dummy_gbound,gspfft_igfft) 318 ABI_FREE(dummy_gbound) 319 320 if (Psps%usepaw==1) then 321 ! Prepare the onsite contributions on the GW FFT mesh. 322 ABI_MALLOC(gw_gfft,(3,nfft)) 323 q0=zero 324 call get_gftt(ngfft_gw,q0,Cryst%gmet,gw_gsq,gw_gfft) ! The set of plane waves in the FFT Box. 325 ABI_DT_MALLOC(Pwij_fft,(Psps%ntypat)) 326 call pawpwij_init(Pwij_fft,nfft,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 327 end if 328 end if 329 330 ! Setup weight (2 for spin unpolarized systems, 1 for polarized). 331 ! spin_fact is used to normalize the occupation factors to one. 332 ! Consider also the AFM case. 333 SELECT CASE (nsppol) 334 CASE (1) 335 weight=two/Kmesh%nbz; spin_fact=half 336 if (Wfd%nspden==2) then 337 weight=one/Kmesh%nbz; spin_fact=half 338 end if 339 if (nspinor==2) then 340 weight=one/Kmesh%nbz; spin_fact=one 341 end if 342 343 CASE (2) 344 weight=one/Kmesh%nbz; spin_fact=one 345 346 CASE DEFAULT 347 MSG_BUG("Wrong nsppol") 348 END SELECT 349 350 ! Weight for points in the IBZ_q 351 wtk_ltg(:)=1 352 if (Ep%symchi==1) then 353 do ik_bz=1,Ltg_q%nbz 354 wtk_ltg(ik_bz)=0 355 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only k in IBZ_q 356 wtk_ltg(ik_bz)=SUM(Ltg_q%wtksym(:,:,ik_bz)) 357 end do 358 end if 359 360 write(msg,'(a,i3,a)')' Q-points for long wave-length limit. # ',Ep%nqlwl,ch10 361 do iqlwl=1,Ep%nqlwl 362 write(msg_tmp,'(1x,i5,a,2x,3f12.6,a)') iqlwl,')',Ep%qlwl(:,iqlwl),ch10 363 msg=TRIM(msg)//msg_tmp 364 end do 365 call wrtout(std_out,msg,'COLL') 366 367 write(msg,'(a,i2,2a,i2)')& 368 & ' Using spectral method for the imaginary part = ',Ep%spmeth,ch10,& 369 & ' Using symmetries to sum only over the IBZ_q = ',Ep%symchi 370 call wrtout(std_out,msg,'COLL') 371 372 if (use_tr) then 373 ! Special care has to be taken in metals and/or spin dependent systems 374 ! as Wfs_val might contain unoccupied states. 375 call wrtout(std_out,' Using faster algorithm based on time reversal symmetry.') 376 else 377 call wrtout(std_out,' Using slow algorithm without time reversal symmetry.') 378 end if 379 380 ! Evaluate oscillator matrix elements btw partial waves. Note q=Gamma 381 if (Psps%usepaw==1) then 382 ABI_DT_MALLOC(Pwij,(Psps%ntypat)) 383 call pawpwij_init(Pwij,Ep%npwepG0,(/zero,zero,zero/),Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 384 385 ABI_DT_MALLOC(Cprj1_bz,(Cryst%natom,nspinor)) 386 call pawcprj_alloc(Cprj1_bz,0,Wfd%nlmn_atm) 387 ABI_DT_MALLOC(Cprj2_bz,(Cryst%natom,nspinor)) 388 call pawcprj_alloc(Cprj2_bz,0,Wfd%nlmn_atm) 389 390 ABI_DT_MALLOC(Cprj1_ibz,(Cryst%natom,nspinor)) 391 call pawcprj_alloc(Cprj1_ibz,0,Wfd%nlmn_atm) 392 ABI_DT_MALLOC(Cprj2_ibz,(Cryst%natom,nspinor)) 393 call pawcprj_alloc(Cprj2_ibz,0,Wfd%nlmn_atm) 394 if (Dtset%pawcross==1) then 395 ABI_MALLOC(ur_ae1,(nfftf_tot*nspinor)) 396 ABI_MALLOC(ur_ae_onsite1,(nfftf_tot*nspinor)) 397 ABI_MALLOC(ur_ps_onsite1,(nfftf_tot*nspinor)) 398 ABI_MALLOC(ur_ae2,(nfftf_tot*nspinor)) 399 ABI_MALLOC(ur_ae_onsite2,(nfftf_tot*nspinor)) 400 ABI_MALLOC(ur_ps_onsite2,(nfftf_tot*nspinor)) 401 ABI_MALLOC(igfftepsG0f,(Ep%npwepG0)) 402 ABI_MALLOC(tabrf_k,(nfftf_tot)) 403 end if 404 end if 405 406 ABI_MALLOC(rhotwg,(Ep%npwe*dim_rtwg)) 407 ABI_MALLOC(tabr_k,(nfft)) 408 ABI_MALLOC(ur1_kibz,(nfft*nspinor)) 409 ABI_MALLOC(ur2_kibz,(nfft*nspinor)) 410 411 ABI_MALLOC(usr1_k,(nfft*nspinor)) 412 ABI_MALLOC(ur2_k,(nfft*nspinor)) 413 ! 414 ! Tables for the FFT of the oscillators. 415 ! a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere). 416 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 417 ABI_MALLOC(igffteps0,(Gsph_epsG0%ng)) 418 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2)) 419 call gsph_fft_tabs(Gsph_epsG0, [0, 0, 0], gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igffteps0) 420 if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 421 #ifdef FC_IBM 422 ! XLF does not deserve this optimization (problem with [v67mbpt][t03]) 423 use_padfft = 0 424 #endif 425 if (use_padfft==0) then 426 ABI_FREE(gw_gbound) 427 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft)) 428 end if 429 if (Dtset%pawcross==1) then 430 ABI_MALLOC(gboundf,(2*mgfftf+8,2)) 431 call gsph_fft_tabs(Gsph_epsG0,(/0,0,0/),mgfftf,ngfftf,use_padfftf,gboundf,igfftepsG0f) 432 if ( ANY(gw_fftalga == (/2,4/)) ) use_padfftf=0 433 if (use_padfftf==0) then 434 ABI_FREE(gboundf) 435 ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf)) 436 end if 437 end if 438 439 ! TODO this table can be calculated for each k-point 440 my_nbbpks=0; allup="All"; got=0 441 ABI_MALLOC(bbp_ks_distrb,(mband,mband,Kmesh%nbz,nsppol)) 442 ABI_MALLOC(bbp_mask,(mband,mband)) 443 444 do spin=1,nsppol 445 do ik_bz=1,Kmesh%nbz 446 447 if (Ep%symchi==1) then 448 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q 449 end if 450 ik_ibz=Kmesh%tab(ik_bz) 451 452 call chi0_bbp_mask(Ep,use_tr,QP_BSt,mband,ik_ibz,ik_ibz,spin,spin_fact,bbp_mask) 453 454 call wfd_distribute_bbp(Wfd,ik_ibz,spin,allup,my_nbbp,bbp_ks_distrb(:,:,ik_bz,spin),got,bbp_mask) 455 my_nbbpks = my_nbbpks + my_nbbp 456 end do 457 end do 458 459 ABI_FREE(bbp_mask) 460 461 write(msg,'(a,i0,a)')" Will sum ",my_nbbpks," (b,b',k,s) states in chi0q0." 462 call wrtout(std_out,msg,'PERS') 463 464 SELECT CASE (Ep%spmeth) 465 CASE (0) 466 call wrtout(std_out,' Calculating chi0(q=(0,0,0),omega,G,G")',"COLL") 467 ABI_MALLOC(green_w,(Ep%nomega)) 468 469 CASE (1, 2) 470 call wrtout(std_out,' Calculating Im chi0(q=(0,0,0),omega,G,G")','COLL') 471 ! 472 ! === Find max and min resonant transitions for this q, report values for this processor === 473 call make_transitions(Wfd,1,Ep%nbnds,nbvw,nsppol,Ep%symchi,Cryst%timrev,GW_TOL_DOCC,& 474 & max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,qp_energy,qp_occ,(/zero,zero,zero/),bbp_ks_distrb) 475 476 !FIXME there is a problem in make_transitions due to MPI_enreg 477 !ltest=(MPI_enreg%gwpara==0.or.MPI_enreg%gwpara==2) 478 !ABI_CHECK(ltest,"spectral method with gwpara==1 not coded") 479 ! 480 ! === Calculate frequency dependent weights for Kramers Kronig transform === 481 ABI_MALLOC(omegasf,(Ep%nomegasf)) 482 ABI_MALLOC(kkweight,(Ep%nomegasf,Ep%nomega)) 483 !my_wl=1; my_wr=Ep%nomegasf 484 call setup_spectral(Ep%nomega,Ep%omega,Ep%nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,& 485 & 0,Ep%zcut,zero,my_wl,my_wr,kkweight) 486 487 if (.not.use_tr) then 488 MSG_BUG('Hilbert transform requires time-reversal') 489 end if 490 491 ! allocate heads and wings of the spectral function. 492 ABI_MALLOC(sf_head,(3,3,my_wl:my_wr)) 493 ABI_MALLOC(sf_lwing,(Ep%npwe,my_wl:my_wr,3)) 494 ABI_MALLOC(sf_uwing,(Ep%npwe,my_wl:my_wr,3)) 495 sf_head=czero; sf_lwing=czero; sf_uwing=czero 496 497 memreq = two*gwpc*Ep%npwe**2*(my_wr-my_wl+1)*b2Gb 498 write(msg,'(a,f10.4,a)')' memory required per spectral point: ',two*gwpc*Ep%npwe**2*b2Mb,' [Mb]' 499 call wrtout(std_out,msg,'PERS') 500 write(msg,'(a,f10.4,a)')' memory required by sf_chi0q0: ',memreq,' [Gb]' 501 call wrtout(std_out,msg,'PERS') 502 if (memreq > two) then 503 MSG_WARNING(' Memory required for sf_chi0q0 is larger than 2.0 Gb!') 504 end if 505 ABI_STAT_MALLOC(sf_chi0,(Ep%npwe,Ep%npwe,my_wl:my_wr), ierr) 506 ABI_CHECK(ierr==0, 'out-of-memory in sf_chi0q0') 507 sf_chi0=czero_gw 508 509 CASE DEFAULT 510 MSG_BUG("Wrong spmeth") 511 END SELECT 512 513 nkpt_summed=Kmesh%nbz 514 if (Ep%symchi/=0) then 515 nkpt_summed=Ltg_q%nibz_ltg 516 call littlegroup_print(Ltg_q,std_out,Dtset%prtvol,'COLL') 517 end if 518 519 ABI_DT_MALLOC(vkbr,(Kmesh%nibz)) 520 gradk_not_done=.TRUE. 521 522 write(msg,'(a,i6,a)')' Calculation status ( ',nkpt_summed,' to be completed):' 523 call wrtout(std_out,msg,'COLL') 524 ! 525 ! ============================================ 526 ! === Begin big fat loop over transitions ==== 527 ! ============================================ 528 chi0=czero_gw; chi0_sumrule =zero 529 530 ! Loop on spin to calculate $\chi_{\up,\up} + \chi_{\down,\down}$ 531 do spin=1,nsppol 532 if (ALL(bbp_ks_distrb(:,:,:,spin) /= Wfd%my_rank)) CYCLE 533 534 ! Loop over k-points in the BZ. 535 do ik_bz=1,Kmesh%nbz 536 if (Ep%symchi==1) then 537 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q 538 end if 539 540 if (ALL(bbp_ks_distrb(:,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE 541 542 write(msg,'(2(a,i4),a,i2,a,i3)')' ik= ',ik_bz,'/',Kmesh%nbz,' spin=',spin,' done by mpi rank:',Wfd%my_rank 543 call wrtout(std_out,msg,'PERS') 544 545 ! Get ik_ibz, non-symmorphic phase and symmetries from ik_bz. 546 call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt) 547 tabr_k=ktabr(:,ik_bz) ! Table for rotated FFT points 548 spinrot_kbz(:)=Cryst%spinrot(:,isym_k) 549 if (Dtset%pawcross==1) tabrf_k(:) = ktabrf(:,ik_bz) 550 551 istwf_k = Wfd%istwfk(ik_ibz) 552 npw_k = Wfd%npwarr(ik_ibz) 553 kg_k => Wfd%Kdata(ik_ibz)%kg_k 554 555 if (Psps%usepaw==0.and.Ep%inclvkb/=0.and.gradk_not_done(ik_ibz)) then 556 ! Include term <n,k|[Vnl,iqr]|n"k>' for q->0. 557 call vkbr_init(vkbr(ik_ibz),Cryst,Psps,Ep%inclvkb,istwf_k,npw_k,Kmesh%ibz(:,ik_ibz),kg_k) 558 gradk_not_done(ik_ibz)=.FALSE. 559 end if 560 561 ! Loop over "conduction" states. 562 do band1=1,Ep%nbnds 563 if (ALL(bbp_ks_distrb(band1,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE 564 565 ug1 => Wfd%Wave(band1,ik_ibz,spin)%ug 566 call wfd_get_ur(Wfd,band1,ik_ibz,spin,ur1_kibz) 567 568 if (Psps%usepaw==1) then 569 call wfd_get_cprj(Wfd,band1,ik_ibz,spin,Cryst,Cprj1_ibz,sorted=.FALSE.) 570 call pawcprj_copy(Cprj1_ibz,Cprj1_bz) 571 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_bz) 572 if (Dtset%pawcross==1) then 573 call wfd_paw_get_aeur(Wfdf,band1,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae1,ur_ae_onsite1,ur_ps_onsite1) 574 end if 575 end if 576 577 do band2=1,Ep%nbnds ! Loop over "valence" states. 578 ! write(6,*) "ik,band1,band2",ik_bz,band1,band2 579 580 ! ----------------- cRPA for U 581 !debug if (.not.luwindow.AND.dtset%ucrpa==1.AND.band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)& 582 !debug& .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) CYCLE 583 if (luwindow.AND.dtset%ucrpa==1.AND.((KS_Bst%eig(band1,ik_ibz,spin)-KS_Bst%fermie)<=dtset%ucrpa_window(2)) & 584 & .AND.((KS_Bst%eig(band1,ik_ibz,spin)-KS_Bst%fermie)>=dtset%ucrpa_window(1)) & 585 & .AND.((KS_Bst%eig(band2,ik_ibz,spin)-KS_Bst%fermie)<=dtset%ucrpa_window(2)) & 586 & .AND.((KS_Bst%eig(band2,ik_ibz,spin)-KS_Bst%fermie)>=dtset%ucrpa_window(1))) CYCLE 587 ! ----------------- cRPA for U 588 589 if (bbp_ks_distrb(band1,band2,ik_bz,spin) /= Wfd%my_rank) CYCLE 590 591 deltaf_b1b2 =spin_fact*(qp_occ(band1,ik_ibz,spin)-qp_occ(band2,ik_ibz,spin)) 592 deltaeKS_b1b2= ks_energy(band1,ik_ibz,spin) - ks_energy(band2,ik_ibz,spin) 593 deltaeGW_b1b2= qp_energy(band1,ik_ibz,spin) - qp_energy(band2,ik_ibz,spin) 594 595 if (Ep%gwcomp==0) then ! Skip negligible transitions. 596 if (ABS(deltaf_b1b2) < GW_TOL_DOCC) CYCLE 597 else 598 ! when the completeness trick is used, we need to also consider transitions with vanishing deltaf 599 !Rangel Correction for metals 600 !if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC) CYCLE 601 if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC .and. ( ABS(deltaf_b1b2)< GW_TOL_DOCC .or. band1<band2)) CYCLE 602 end if 603 604 ug2 => Wfd%Wave(band2,ik_ibz,spin)%ug 605 call wfd_get_ur(Wfd,band2,ik_ibz,spin,ur2_kibz) 606 607 if (Psps%usepaw==1) then 608 call wfd_get_cprj(Wfd,band2,ik_ibz,spin,Cryst,Cprj2_ibz,sorted=.FALSE.) 609 call pawcprj_copy(Cprj2_ibz,Cprj2_bz) 610 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj2_bz) 611 if (Dtset%pawcross==1) then 612 call wfd_paw_get_aeur(Wfdf,band2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae2,ur_ae_onsite2,ur_ps_onsite2) 613 end if 614 end if 615 616 SELECT CASE (Ep%spmeth) 617 CASE (0) 618 ! Adler-Wiser expression. 619 ! Add small imaginary of the Time-Ordered resp function but only for non-zero real omega FIXME What about metals? 620 621 if (.not.use_tr) then 622 ! Adler-Wiser without time-reversal. 623 do io=1,Ep%nomega 624 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0,one_pole) 625 end do 626 else 627 if (Ep%gwcomp==0) then ! cannot be completely skipped in case of completeness correction 628 if (band1<band2) CYCLE ! Here we GAIN a factor ~2 629 end if 630 631 do io=1,Ep%nomega 632 !Rangel: In metals, the intra-band transitions term does not contain the antiresonant part 633 !if(abs(deltaeGW_b1b2)>GW_TOL_W0) green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0) 634 if (band1==band2) green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0,one_pole) 635 if (band1/=band2) green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0,two_poles) 636 637 if (Ep%gwcomp==1) then ! Calculate the completeness correction 638 numerator= -spin_fact*qp_occ(band2,ik_ibz,spin) 639 deltaeGW_enhigh_b2=en_high-qp_energy(band2,ik_ibz,spin) 640 ! Completeness correction is NOT valid for real frequencies 641 if (REAL(Ep%omega(io))<GW_TOL_W0) then 642 green_enhigh_w(io) = g0g0w(Ep%omega(io),numerator,deltaeGW_enhigh_b2,Ep%zcut,GW_TOL_W0,two_poles) 643 else 644 green_enhigh_w(io) = czero_gw 645 endif 646 ! 647 !Rangel Correction for metals 648 !if (deltaf_b1b2<0.d0) then 649 if (band1>=band2 .and. abs(deltaf_b1b2) > GW_TOL_DOCC) then 650 green_w(io)= green_w(io) - green_enhigh_w(io) 651 else ! Disregard green_w, since it is already accounted for through the time-reversal 652 green_w(io)= - green_enhigh_w(io) 653 end if 654 end if !gwcomp==1 655 end do !io 656 657 if (Ep%gwcomp==1.and.band1==band2) then 658 ! Add the "delta part", symmetrization is done inside the routine. 659 call calc_wfwfg(tabr_k,itim_k,spinrot_kbz,nfft,nspinor,ngfft_gw,ur2_kibz,ur2_kibz,wfwfg) 660 661 if (Psps%usepaw==1) then 662 call paw_rho_tw_g(nfft,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,gw_gfft,& 663 & Cprj2_bz,Cprj2_bz,Pwij_fft,wfwfg) 664 665 ! Add PAW cross term 666 if (Dtset%pawcross==1) then 667 call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,& 668 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,& 669 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,& 670 & dim_rtwg,wfwfg) 671 end if 672 end if 673 674 qzero=.TRUE. 675 call completechi0_deltapart(ik_bz,qzero,Ep%symchi,Ep%npwe,Gsph_FFT%ng,Ep%nomega,nspinor,& 676 & nfft,ngfft_gw,gspfft_igfft,Gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0) 677 end if 678 end if ! use_tr 679 680 CASE (1, 2) 681 ! Spectral method, here time-reversal is always assumed. 682 if (deltaeGW_b1b2<0) CYCLE 683 call approxdelta(Ep%nomegasf,omegasf,deltaeGW_b1b2,Ep%spsmear,iomegal,iomegar,wl,wr,Ep%spmeth) 684 END SELECT 685 686 ! FFT of u^*_{b1,k}(r) u_{b2,k}(r) and (q,G=0) limit using small q and k.p perturbation theory 687 call rho_tw_g(nspinor,Ep%npwe,nfft,ndat1,ngfft_gw,1,use_padfft,igffteps0,gw_gbound,& 688 & ur1_kibz,itim_k,tabr_k,ph_mkt,spinrot_kbz,& 689 & ur2_kibz,itim_k,tabr_k,ph_mkt,spinrot_kbz,& 690 & dim_rtwg,rhotwg) 691 692 if (Psps%usepaw==0) then 693 ! Matrix elements of i[H,r] for NC pseudopotentials. 694 rhotwx = nc_ihr_comm(vkbr(ik_ibz),cryst,psps,npw_k,nspinor,istwf_k,Ep%inclvkb,Kmesh%ibz(:,ik_ibz),ug1,ug2,kg_k) 695 696 else 697 ! 1) Add PAW onsite contribution, projectors are already in the BZ. 698 call paw_rho_tw_g(Ep%npwe,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,& 699 & Cprj1_bz,Cprj2_bz,Pwij,rhotwg) 700 701 ! 2) Matrix elements of i[H,r] for PAW. 702 rhotwx = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug1,ug2,kg_k,Cprj1_ibz,Cprj2_ibz,HUr) 703 704 ! Add PAW cross term 705 if (Dtset%pawcross==1) then 706 call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,& 707 & ur_ae1,ur_ae_onsite1,ur_ps_onsite1,itim_k,tabrf_k,ph_mkt,spinrot_kbz,& 708 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,& 709 & dim_rtwg,rhotwg) 710 ! Add cross-term contribution to the commutator 711 if (Dtset%userib/=111) then 712 call paw_cross_ihr_comm(rhotwx,nspinor,nfftf_tot,Cryst,Pawfgrtab,Paw_onsite,& 713 & ur_ae1,ur_ae2,ur_ae_onsite1,ur_ae_onsite2,Cprj1_ibz,Cprj2_ibz) 714 end if 715 end if 716 end if 717 718 ! Treat a possible degeneracy between v and c. 719 if (ABS(deltaeKS_b1b2)>GW_TOL_W0) then 720 rhotwx=-rhotwx/deltaeKS_b1b2 721 else 722 rhotwx=czero_gw 723 end if 724 725 SELECT CASE (Ep%spmeth) 726 CASE (0) 727 728 ! ---------------- Ucrpa (begin) 729 if(dtset%ucrpa>=1.and..not.luwindow) then 730 fac=one 731 fac1=zero 732 fac2=zero 733 fac3=zero 734 fac4=one 735 m1=-1 736 m2=-1 737 if(dtset%ucrpa<=2) then 738 call flush_unit(std_out) 739 if ( band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)& 740 & .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then 741 ! if(dtset%prtvol>=10)write(6,*)"calculation is in progress",band1,band2,ucrpa_bands(1),ucrpa_bands(2) 742 do iat=1, cryst%nattyp(itypatcor) 743 do ispinor1=1,nspinor 744 do m1=1,2*lcor+1 745 fac1=fac1+ real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 746 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))) 747 fac2=fac2+ real(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)*& 748 & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1))) 749 do ispinor2=1,nspinor 750 do m2=1,2*lcor+1 751 fac=fac - real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 752 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* & 753 & coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)*& 754 & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2))) 755 fac3=fac3 + real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 756 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* & 757 & coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)*& 758 & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2))) 759 ! if(dtset%prtvol>=10)write(6,*) fac,fac3 760 enddo 761 enddo 762 ! if(dtset%prtvol>=10)write(6,*) fac,fac3,fac1,fac2,fac1*fac2 763 enddo 764 enddo 765 enddo 766 fac4=fac 767 ! fac=zero 768 if(dtset%ucrpa==1) fac=zero 769 ! write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0" 770 endif 771 else if (dtset%ucrpa==3) then 772 if (band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)) then 773 do iat=1, cryst%nattyp(itypatcor) 774 do ispinor1=1,nspinor 775 do m1=1,2*lcor+1 776 fac2=fac2-real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 777 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))) 778 enddo 779 enddo 780 enddo 781 if(dtset%ucrpa==4) fac2=zero 782 endif 783 if (band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then 784 do iat=1, cryst%nattyp(itypatcor) 785 do ispinor1=1,nspinor 786 do m1=1,2*lcor+1 787 fac3=fac3-real(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)*& 788 & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1))) 789 enddo 790 enddo 791 enddo 792 if(dtset%ucrpa==4) fac3=zero 793 endif 794 fac=real(fac2*fac3) 795 endif 796 ! if(dtset%prtvol>=10) write(6,'(6i4,e15.5,a)') ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0" 797 ! if(abs(fac-one)>0.00001) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0" 798 ! if(dtset%prtvol>=10) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac4," q==0" 799 green_w=green_w*fac 800 endif 801 ! ---------------- Ucrpa (end) 802 803 ! Adler-Wiser expression, to be consistent here we use the KS eigenvalues (?) 804 call accumulate_chi0_q0(ik_bz,isym_k,itim_k,Ep%gwcomp,nspinor,Ep%npwepG0,Ep,& 805 & Cryst,Ltg_q,Gsph_epsG0,chi0,rhotwx,rhotwg,green_w,green_enhigh_w,deltaf_b1b2,chi0_head,chi0_lwing,chi0_uwing) 806 807 CASE (1, 2) 808 ! Spectral method, to be consistent here we use the KS eigenvalues. 809 call accumulate_sfchi0_q0(ik_bz,isym_k,itim_k,nspinor,Ep%symchi,Ep%npwepG0,Ep%npwe,Cryst,Ltg_q,& 810 & Gsph_epsG0,deltaf_b1b2,my_wl,iomegal,wl,my_wr,iomegar,wr,rhotwx,rhotwg,Ep%nomegasf,sf_chi0,sf_head,sf_lwing,sf_uwing) 811 812 CASE DEFAULT 813 MSG_BUG("Wrong spmeth") 814 END SELECT 815 816 ! Accumulating the sum rule on chi0. Eq. (5.284) in G.D. Mahan Many-Particle Physics 3rd edition. 817 factor=spin_fact*qp_occ(band2,ik_ibz,spin) 818 819 call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_b1b2,& 820 & Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule) 821 822 if (Ep%gwcomp==1) then 823 ! Include also the completeness correction in the sum rule. 824 factor=-spin_fact*qp_occ(band2,ik_ibz,spin) 825 call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_enhigh_b2,& 826 & Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule) 827 if (band1==Ep%nbnds) then 828 chi0_sumrule(:)=chi0_sumrule(:) + wtk_ltg(ik_bz)*spin_fact*qp_occ(band2,ik_ibz,spin)*deltaeGW_enhigh_b2 829 end if 830 end if 831 832 end do !band2 833 end do !band1 834 835 if (Psps%usepaw==0.and.Ep%inclvkb/=0.and.Ep%symchi==1) then 836 call vkbr_free(vkbr(ik_ibz)) ! Not need anymore as we loop only over IBZ. 837 end if 838 end do !ik_bz 839 end do !spin 840 841 ABI_FREE(igffteps0) 842 843 call vkbr_free(vkbr) 844 ABI_DT_FREE(vkbr) 845 846 ! === After big fat loop over transitions, now MPI === 847 ! * Master took care of the contribution in case of (metallic|spin) polarized systems. 848 SELECT CASE (Ep%spmeth) 849 CASE (0) 850 ! Adler-Wiser expression 851 ! Sum contributions from each proc. 852 ! Looping on frequencies to avoid problems with the size of the MPI packet. 853 do io=1,Ep%nomega 854 call xmpi_sum(chi0(:,:,io),comm,ierr) 855 end do 856 857 CASE (1, 2) 858 ! Spectral method. 859 call hilbert_transform(Ep%npwe,Ep%nomega,Ep%nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,Ep%spmeth) 860 861 if (allocated(sf_chi0)) then 862 ABI_FREE(sf_chi0) 863 end if 864 865 ! Sum contributions from each proc === 866 ! Looping on frequencies to avoid problems with the size of the MPI packet 867 do io=1,Ep%nomega 868 call xmpi_sum(chi0(:,:,io),comm,ierr) 869 end do 870 871 call hilbert_transform_headwings(Ep%npwe,Ep%nomega,Ep%nomegasf,& 872 & my_wl,my_wr,kkweight,sf_lwing,sf_uwing,sf_head,chi0_lwing,& 873 & chi0_uwing,chi0_head,Ep%spmeth) 874 875 CASE DEFAULT 876 MSG_BUG("Wrong spmeth") 877 END SELECT 878 879 ! Divide by the volume 880 !$OMP PARALLEL WORKSHARE 881 chi0=chi0*weight/Cryst%ucvol 882 !$OMP END PARALLEL WORKSHARE 883 884 ! Collect sum rule. pi comes from Im[1/(x-ieta)] = pi delta(x) 885 call xmpi_sum(chi0_sumrule,comm,ierr) 886 chi0_sumrule=chi0_sumrule * pi * weight / Cryst%ucvol 887 888 ! Collect heads and wings. 889 call xmpi_sum(chi0_head,comm,ierr) 890 call xmpi_sum(chi0_lwing,comm,ierr) 891 call xmpi_sum(chi0_uwing,comm,ierr) 892 893 chi0_head = chi0_head * weight/Cryst%ucvol 894 do io=1,Ep%nomega ! Tensor in the basis of the reciprocal lattice vectors. 895 chi0_head(:,:,io) = MATMUL(chi0_head(:,:,io),Cryst%gmet) * (two_pi**2) 896 end do 897 chi0_lwing = chi0_lwing * weight/Cryst%ucvol 898 chi0_uwing = chi0_uwing * weight/Cryst%ucvol 899 900 ! =============================================== 901 ! ==== Symmetrize chi0 in case of AFM system ==== 902 ! =============================================== 903 ! Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$. 904 ! Works only in the case of magnetic group Shubnikov type IV. 905 if (Cryst%use_antiferro) then 906 call symmetrize_afm_chi0(Cryst,Gsph_epsG0,Ltg_q,Ep%npwe,Ep%nomega,chi0,chi0_head,chi0_lwing,chi0_uwing) 907 end if 908 909 ! =================================================== 910 ! ==== Construct heads and wings from the tensor ==== 911 ! =================================================== 912 do io=1,Ep%nomega 913 do ig=2,Ep%npwe 914 wng = chi0_uwing(ig,io,:) 915 chi0(1,ig,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G") 916 wng = chi0_lwing(ig,io,:) 917 chi0(ig,1,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G") 918 end do 919 chq = MATMUL(chi0_head(:,:,io), Ep%qlwl(:,1)) 920 chi0(1,1,io) = vdotw(Ep%qlwl(:,1),chq,Cryst%gmet,"G") ! Use user-defined small q 921 end do 922 923 ! Impose Hermiticity (valid only for zero or purely imaginary frequencies) 924 ! MG what about metals, where we have poles around zero? 925 ! do io=1,Ep%nomega 926 ! if (ABS(REAL(Ep%omega(io)))<0.00001) then 927 ! do ig2=1,Ep%npwe 928 ! do ig1=1,ig2-1 929 ! chi0(ig2,ig1,io)=GWPC_CONJG(chi0(ig1,ig2,io)) 930 ! end do 931 ! end do 932 ! end if 933 ! end do 934 ! 935 ! ===================== 936 ! ==== Free memory ==== 937 ! ===================== 938 ABI_FREE(bbp_ks_distrb) 939 ABI_FREE(rhotwg) 940 ABI_FREE(tabr_k) 941 ABI_FREE(ur1_kibz) 942 ABI_FREE(ur2_kibz) 943 ABI_FREE(usr1_k) 944 ABI_FREE(ur2_k) 945 ABI_FREE(gw_gbound) 946 947 if (Dtset%pawcross==1) then 948 ABI_FREE(gboundf) 949 end if 950 951 if (allocated(green_enhigh_w)) then 952 ABI_FREE(green_enhigh_w) 953 end if 954 if (allocated(gw_gfft)) then 955 ABI_FREE(gw_gfft) 956 end if 957 if (allocated(wfwfg)) then 958 ABI_FREE(wfwfg) 959 end if 960 if (allocated(kkweight)) then 961 ABI_FREE(kkweight) 962 end if 963 if (allocated(omegasf)) then 964 ABI_FREE(omegasf) 965 end if 966 if (allocated(green_w)) then 967 ABI_FREE(green_w) 968 end if 969 970 if (allocated(sf_head)) then 971 ABI_FREE(sf_head) 972 end if 973 if (allocated(sf_lwing)) then 974 ABI_FREE(sf_lwing) 975 end if 976 if (allocated(sf_uwing)) then 977 ABI_FREE(sf_uwing) 978 end if 979 if (allocated(gspfft_igfft)) then 980 ABI_FREE(gspfft_igfft) 981 end if 982 983 call gsph_free(Gsph_FFT) 984 985 if (Psps%usepaw==1) then ! deallocation for PAW. 986 call pawcprj_free(Cprj1_bz ) 987 ABI_DT_FREE(Cprj1_bz) 988 call pawcprj_free(Cprj2_bz ) 989 ABI_DT_FREE(Cprj2_bz) 990 call pawcprj_free(Cprj1_ibz ) 991 ABI_DT_FREE(Cprj1_ibz) 992 call pawcprj_free(Cprj2_ibz ) 993 ABI_DT_FREE(Cprj2_ibz) 994 call pawpwij_free(Pwij) 995 ABI_DT_FREE(Pwij) 996 if (allocated(Pwij_fft)) then 997 call pawpwij_free(Pwij_fft) 998 ABI_DT_FREE(Pwij_fft) 999 end if 1000 call pawhur_free(Hur) 1001 ABI_DT_FREE(Hur) 1002 if (Dtset%pawcross==1) then 1003 ABI_FREE(ur_ae1) 1004 ABI_FREE(ur_ae_onsite1) 1005 ABI_FREE(ur_ps_onsite1) 1006 ABI_FREE(ur_ae2) 1007 ABI_FREE(ur_ae_onsite2) 1008 ABI_FREE(ur_ps_onsite2) 1009 ABI_FREE(tabrf_k) 1010 ABI_FREE(gboundf) 1011 ABI_FREE(igfftepsG0f) 1012 end if 1013 end if 1014 1015 if(dtset%ucrpa>=1) then 1016 ABI_DEALLOCATE(coeffW_BZ) 1017 endif 1018 1019 call cwtime(cpu_time,wall_time,gflops,"stop") 1020 write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time 1021 1022 DBG_EXIT("COLL") 1023 1024 end subroutine cchi0q0