TABLE OF CONTENTS
ABINIT/chi0q0_intraband [ Functions ]
NAME
chi0q0_intraband
FUNCTION
Calculate chi0 in the limit q-->0
COPYRIGHT
Copyright (C) 2010-2018 ABINIT group (MG) 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) 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%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) Ep%nsppol=1 for unpolarized, 2 for spin-polarized 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 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 Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
OUTPUT
chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix for wavevector qq, and frequencies defined by Ep%omega
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
TODO
Check npwepG0 before Switching on umklapp
PARENTS
screening
CHILDREN
assemblychi0_sym,get_bz_item,getnel,gsph_fft_tabs,kmesh_free,kmesh_init littlegroup_free,littlegroup_init,littlegroup_print,pack_eneocc paw_rho_tw_g,paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free pawhur_free,pawhur_init,pawpwij_free,pawpwij_init,print_arr,rho_tw_g rotate_fft_mesh,symmetrize_afm_chi0,unpack_eneocc,vkbr_free,vkbr_init wfd_change_ngfft,wfd_distribute_bands,wfd_get_cprj,wfd_get_ur,wrtout xmpi_sum
SOURCE
75 #if defined HAVE_CONFIG_H 76 #include "config.h" 77 #endif 78 79 #include "abi_common.h" 80 81 subroutine chi0q0_intraband(Wfd,Cryst,Ep,Psps,BSt,Gsph_epsG0,Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,use_tr,usepawu,& 82 & ngfft_gw,chi0,chi0_head,chi0_lwing,chi0_uwing) 83 84 use defs_basis 85 use defs_datatypes 86 use m_xmpi 87 use m_errors 88 use m_profiling_abi 89 use m_wfd 90 91 use m_gwdefs, only : GW_TOL_DOCC, GW_TOL_W0, czero_gw, em1params_t, g0g0w 92 use m_geometry, only : vdotw 93 use m_numeric_tools, only : print_arr 94 use m_crystal, only : crystal_t 95 use m_fft_mesh, only : rotate_FFT_mesh 96 use m_occ, only : getnel 97 use m_ebands, only : pack_eneocc, unpack_eneocc 98 use m_bz_mesh, only : kmesh_t, kmesh_init, kmesh_free, get_BZ_item, & 99 & littlegroup_t, littlegroup_print, littlegroup_free, littlegroup_init 100 use m_gsphere, only : gsphere_t, gsph_fft_tabs 101 use m_oscillators, only : rho_tw_g 102 use m_pawhr, only : pawhur_t, pawhur_free, pawhur_init, paw_ihr 103 use m_vkbr, only : vkbr_t, vkbr_free, vkbr_init, nc_ihr_comm 104 use m_chi0, only : assemblychi0_sym, symmetrize_afm_chi0 105 use m_pawang, only : pawang_type 106 use m_pawrad, only : pawrad_type 107 use m_pawtab, only : pawtab_type 108 use m_paw_ij, only : paw_ij_type 109 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy 110 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g 111 112 !This section has been created automatically by the script Abilint (TD). 113 !Do not modify the following lines by hand. 114 #undef ABI_FUNC 115 #define ABI_FUNC 'chi0q0_intraband' 116 use interfaces_14_hidewrite 117 use interfaces_65_paw 118 !End of the abilint section 119 120 implicit none 121 122 !Arguments ------------------------------------ 123 !scalars 124 integer,intent(in) :: usepawu 125 logical,intent(in) :: use_tr 126 type(ebands_t),intent(in) :: BSt 127 type(crystal_t),intent(in) :: Cryst 128 type(em1params_t),intent(in) :: Ep 129 type(gsphere_t),intent(in) :: Gsph_epsG0 130 type(Pseudopotential_type),intent(in) :: Psps 131 type(Pawang_type),intent(in) :: Pawang 132 type(wfd_t),target,intent(inout) :: Wfd 133 !arrays 134 integer,intent(in) :: ngfft_gw(18) 135 complex(gwpc),intent(out) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega) 136 complex(dpc),intent(out) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3) 137 complex(dpc),intent(out) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3) 138 complex(dpc),intent(out) :: chi0_head(3,3,Ep%nomega) 139 type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw) 140 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 141 type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw) 142 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 143 144 !Local variables ------------------------------ 145 !scalars 146 integer,parameter :: tim_fourdp1=1,two_poles=2,one_pole=1,ndat1=1 147 integer,parameter :: unitdos0=0,option1=1,NOMEGA_PRINTED=15 148 integer :: nqlwl,nband_k,iomega,istwf_k,npw_k,my_nband,lbidx 149 integer :: band,itim_k,ik_bz,ik_ibz,io,isym_k,spin,iqlwl!,gw_eet !ig,ig1,ig2,my_nbbp,my_nbbpks 150 integer :: nkpt_summed,dim_rtwg,use_padfft,gw_fftalga,ifft 151 integer :: kptopt,isym,nsppol,nspinor 152 integer :: comm,ierr,gw_mgfft,use_umklp,inclvkb 153 real(dp) :: spin_fact,deltaf_b1b2,weight 154 real(dp) :: deltaeGW_b1b2,zcut 155 real(dp),parameter :: dummy_dosdeltae=HUGE(zero) 156 real(dp) :: o_entropy,o_nelect,maxocc 157 complex(dpc) :: ph_mkt 158 logical :: iscompatibleFFT !,ltest 159 character(len=500) :: msg,msg_tmp !,allup 160 type(kmesh_t) :: Kmesh 161 type(littlegroup_t) :: Ltg_q 162 type(vkbr_t) :: vkbr 163 !arrays 164 integer :: my_band_list(Wfd%mband) 165 integer,ABI_CONTIGUOUS pointer :: kg_k(:,:) 166 integer,allocatable :: ktabr(:,:),irottb(:,:) 167 !integer :: got(Wfd%nproc) 168 integer,allocatable :: tabr_k(:),igffteps0(:),gw_gbound(:,:) 169 real(dp),parameter :: q0(3)=(/zero,zero,zero/) 170 real(dp) :: kpt(3),dedk(3),kbz(3),spinrot_kbz(4) 171 !real(dp),ABI_CONTIGUOUS pointer :: ks_energy(:,:,:),qp_energy(:,:,:),qp_occ(:,:,:) 172 real(dp) :: shift_ene(BSt%mband,BSt%nkpt,BSt%nsppol) 173 real(dp) :: delta_occ(BSt%mband,BSt%nkpt,BSt%nsppol) 174 !real(dp) :: eigen_vec(BSt%bantot) 175 real(dp) :: o_doccde(BSt%bantot) 176 real(dp) :: eigen_pdelta_vec(BSt%bantot),eigen_mdelta_vec(BSt%bantot) 177 real(dp) :: o_occ_pdelta(BSt%bantot),o_occ_mdelta(BSt%bantot) 178 real(dp) :: delta_ene(BSt%mband,BSt%nkpt,BSt%nsppol) 179 real(dp) :: test_docc(BSt%mband,BSt%nkpt,BSt%nsppol) 180 real(dp),allocatable :: qlwl(:,:) 181 complex(gwpc) :: comm_kbbs(3,Wfd%nspinor**2) 182 complex(dpc),allocatable :: ihr_comm(:,:,:,:,:) 183 complex(gwpc),allocatable :: rhotwg(:) 184 complex(dpc) :: green_w(Ep%nomega) 185 complex(gwpc),allocatable :: ur1(:) 186 complex(gwpc),ABI_CONTIGUOUS pointer :: ug(:) 187 logical :: bmask(Wfd%mband) 188 type(pawcprj_type),allocatable :: Cprj1_bz(:,:),Cprj1_ibz(:,:),Cp_bks(:,:) 189 type(pawpwij_t),allocatable :: Pwij(:) 190 type(pawhur_t),allocatable :: Hur(:) 191 192 !************************************************************************ 193 194 DBG_ENTER("COLL") 195 196 nsppol = Wfd%nsppol 197 nspinor = Wfd%nspinor 198 199 gw_mgfft = MAXVAL(ngfft_gw(1:3)) 200 gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10) 201 202 ! Calculate <k,b1|i[H,r]|k',b2>. 203 inclvkb=2; if (Wfd%usepaw==1) inclvkb=0 204 ABI_MALLOC(ihr_comm,(3,nspinor**2,Wfd%mband,Wfd%nkibz,nsppol)) 205 ihr_comm = czero 206 207 if (Wfd%usepaw==1) then 208 ABI_DT_MALLOC(Cp_bks,(Cryst%natom,nspinor)) 209 call pawcprj_alloc(Cp_bks,0,Wfd%nlmn_atm) 210 ABI_DT_MALLOC(HUr,(Cryst%natom)) 211 if (usepawu/=0) then ! For PAW+LDA+U, precalculate <\phi_i|[Hu,r]|phi_j\>. 212 call pawhur_init(hur,nsppol,Wfd%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij) 213 end if 214 end if 215 216 do spin=1,nsppol 217 do ik_ibz=1,Wfd%nkibz 218 npw_k = Wfd%npwarr(ik_ibz) 219 nband_k= Wfd%nband(ik_ibz,spin) 220 kpt = Wfd%kibz(:,ik_ibz) 221 kg_k => Wfd%Kdata(ik_ibz)%kg_k 222 istwf_k = Wfd%istwfk(ik_ibz) 223 ABI_CHECK(istwf_k==1,"istwf_k/=1 not coded") 224 225 ! Distribute bands. 226 bmask=.FALSE.; bmask(1:nband_k)=.TRUE. ! TODO only bands around EF should be included. 227 call wfd_distribute_bands(Wfd,ik_ibz,spin,my_nband,my_band_list,bmask=bmask) 228 if (my_nband==0) CYCLE 229 230 if (Wfd%usepaw==0.and.inclvkb/=0) then ! Include term <n,k|[Vnl,iqr]|n"k>' for q->0. 231 call vkbr_init(vkbr,Cryst,Psps,inclvkb,istwf_k,npw_k,kpt,kg_k) 232 end if 233 234 do lbidx=1,my_nband 235 band=my_band_list(lbidx) 236 ug => Wfd%Wave(band,ik_ibz,spin)%ug 237 238 if (Wfd%usepaw==0) then 239 ! Matrix elements of i[H,r] for NC pseudopotentials. 240 comm_kbbs = nc_ihr_comm(vkbr,cryst,psps,npw_k,nspinor,istwf_k,inclvkb,Kmesh%ibz(:,ik_ibz),ug,ug,kg_k) 241 else 242 ! Matrix elements of i[H,r] for PAW. 243 call wfd_get_cprj(Wfd,band,ik_ibz,spin,Cryst,Cp_bks,sorted=.FALSE.) 244 comm_kbbs = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug,ug,kg_k,Cp_bks,Cp_bks,HUr) 245 end if 246 247 ihr_comm(:,:,band,ik_ibz,spin) = comm_kbbs 248 end do 249 250 call vkbr_free(vkbr) ! Not need anymore as we loop only over IBZ. 251 end do 252 end do 253 ! 254 ! Gather the commutator on each node. 255 call xmpi_sum(ihr_comm,Wfd%comm,ierr) 256 257 if (Wfd%usepaw==1) then 258 call pawcprj_free(Cp_bks) 259 ABI_DT_FREE(Cp_bks) 260 call pawhur_free(Hur) 261 ABI_DT_FREE(Hur) 262 end if 263 264 nqlwl=1 265 ABI_MALLOC(qlwl,(3,nqlwl)) 266 !qlwl = GW_Q0_DEFAULT(3) 267 qlwl(:,1) = (/0.00001_dp, 0.00002_dp, 0.00003_dp/) 268 ! 269 write(msg,'(a,i3,a)')' Q-points for long wave-length limit in chi0q_intraband. # ',nqlwl,ch10 270 do iqlwl=1,nqlwl 271 write(msg_tmp,'(1x,i5,a,2x,3f12.6,a)') iqlwl,')',qlwl(:,iqlwl),ch10 272 msg=TRIM(msg)//msg_tmp 273 end do 274 call wrtout(std_out,msg,'COLL') 275 ! 276 ! delta_ene = e_{b,k-q} - e_{b,k} = -q. <b,k| i[H,r] |b,k> + O(q^2). 277 delta_ene = zero 278 do spin=1,nsppol 279 do ik_ibz=1,Wfd%nkibz 280 do band=1,Wfd%nband(ik_ibz,spin) 281 dedk = REAL(ihr_comm(:,1,band,ik_ibz,spin)) 282 delta_ene(band,ik_ibz,spin) = -vdotw(qlwl(:,1),dedk,Cryst%gmet,"G") 283 end do 284 end do 285 end do 286 287 maxocc=two/(nsppol*nspinor) 288 289 ! Calculate the occupations at f(e+delta/2). 290 shift_ene = BSt%eig + half*delta_ene 291 292 call pack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,BSt%bantot,shift_ene,eigen_pdelta_vec) 293 294 call getnel(o_doccde,dummy_dosdeltae,eigen_pdelta_vec,o_entropy,BSt%fermie,maxocc,BSt%mband,BSt%nband,& 295 & o_nelect,BSt%nkpt,BSt%nsppol,o_occ_pdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk) 296 write(std_out,*)"nelect1: ",o_nelect 297 ! 298 ! Calculate the occupations at f(e-delta/2). 299 shift_ene = BSt%eig - half*delta_ene 300 301 call pack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,BSt%bantot,shift_ene,eigen_mdelta_vec) 302 303 call getnel(o_doccde,dummy_dosdeltae,eigen_mdelta_vec,o_entropy,BSt%fermie,maxocc,BSt%mband,BSt%nband,& 304 & o_nelect,BSt%nkpt,BSt%nsppol,o_occ_mdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk) 305 write(std_out,*)"nelect2: ",o_nelect 306 ! 307 ! f(e-delta/2) - f(e+delta/2). 308 o_occ_pdelta = o_occ_mdelta - o_occ_pdelta 309 310 call unpack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,o_occ_pdelta,delta_occ) 311 ! 312 ! Expand f(e-delta/2) - f(e+delta/2) up to the first order in the small q. 313 do spin=1,nsppol 314 do ik_ibz=1,Wfd%nkibz 315 do band=1,Wfd%nband(ik_ibz,spin) 316 dedk = REAL(ihr_comm(:,1,band,ik_ibz,spin)) 317 test_docc(band,ik_ibz,spin) = +vdotw(qlwl(:,1),dedk,Cryst%gmet,"G") * BSt%doccde(band,ik_ibz,spin) 318 write(std_out,'(a,3(i0,1x),1x,3es16.8)')" spin,ik_ibz,band, delta_occ: ",& 319 & spin,ik_ibz,band,delta_occ(band,ik_ibz,spin),& 320 & test_docc(band,ik_ibz,spin),delta_occ(band,ik_ibz,spin)-test_docc(band,ik_ibz,spin) 321 end do 322 end do 323 end do 324 325 ! MSG_ERROR("DONE") 326 ! do spin=1,nsppol 327 ! do ik_ibz=1,Wfd%nkibz 328 ! nband_k = Wfd%nband(ik_ibz,spin) 329 ! do band=1,nband_k 330 ! write(std_out,'(a,3i3,2es14.6)')" spin, band, ik_ibz, delta_ene, delta_occ ",& 331 !& spin,band,ik_ibz,delta_ene(band,ik_ibz,spin),delta_occ(band,ik_ibz,spin) 332 ! end do 333 ! end do 334 ! end do 335 336 ABI_FREE(ihr_comm) 337 ABI_FREE(qlwl) 338 339 if ( ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,ngfft_gw) 340 341 ! TODO take into account the case of random k-meshes. 342 kptopt=3 343 call kmesh_init(Kmesh,Cryst,Wfd%nkibz,Wfd%kibz,kptopt) 344 ! 345 !=== Get the FFT index of $ (R^{-1}(r-\tau)) $ === 346 !* S= $\transpose R^{-1}$ and k_BZ = S k_IBZ 347 !* irottb is the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Sk. 348 ABI_MALLOC(irottb,(Wfd%nfftot,Cryst%nsym)) 349 350 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,irottb,iscompatibleFFT) 351 ABI_CHECK(iscompatibleFFT,"FFT mesh not compatible with symmetries") 352 353 ABI_MALLOC(ktabr,(Wfd%nfftot,Kmesh%nbz)) 354 do ik_bz=1,Kmesh%nbz 355 isym=Kmesh%tabo(ik_bz) 356 do ifft=1,Wfd%nfftot 357 ktabr(ifft,ik_bz)=irottb(ifft,isym) 358 end do 359 end do 360 ABI_FREE(irottb) 361 ! 362 ! === Setup weight (2 for spin unpolarized systems, 1 for polarized) === 363 ! spin_fact is used to normalize the occupation factors to one. 364 ! Consider also the AFM case. 365 SELECT CASE (nsppol) 366 CASE (1) 367 weight=two/Kmesh%nbz; spin_fact=half 368 if (Wfd%nspden==2) then 369 weight=one/Kmesh%nbz; spin_fact=half 370 end if 371 if (nspinor==2) then 372 weight=one/Kmesh%nbz; spin_fact=one 373 end if 374 CASE (2) 375 weight=one/Kmesh%nbz; spin_fact=one 376 CASE DEFAULT 377 MSG_BUG("Wrong nsppol") 378 END SELECT 379 380 use_umklp=0 381 call littlegroup_init(q0,Kmesh,Cryst,use_umklp,Ltg_q,Ep%npwepG0,gvec=Gsph_epsG0%gvec) 382 383 write(msg,'(a,i2)')' Using symmetries to sum only over the IBZ_q = ',Ep%symchi 384 call wrtout(std_out,msg,'COLL') 385 ! 386 ! === Evaluate oscillator matrix elements btw partial waves. Note that q=Gamma is used. 387 if (Psps%usepaw==1) then 388 ABI_DT_MALLOC(Pwij,(Psps%ntypat)) 389 call pawpwij_init(Pwij,Ep%npwepG0, [zero,zero,zero], Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 390 391 ABI_DT_MALLOC(Cprj1_bz ,(Cryst%natom,nspinor)) 392 call pawcprj_alloc(Cprj1_bz, 0,Wfd%nlmn_atm) 393 ABI_DT_MALLOC(Cprj1_ibz,(Cryst%natom,nspinor)) 394 call pawcprj_alloc(Cprj1_ibz,0,Wfd%nlmn_atm) 395 end if 396 397 ABI_MALLOC(rhotwg,(Ep%npwe*nspinor**2)) 398 ABI_MALLOC(tabr_k,(Wfd%nfftot)) 399 ABI_MALLOC(ur1,(Wfd%nfft*nspinor)) 400 ! 401 ! Tables for the FFT of the oscillators. 402 ! a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere). 403 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 404 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2)) 405 ABI_MALLOC(igffteps0,(Gsph_epsG0%ng)) 406 407 call gsph_fft_tabs(Gsph_epsG0, [0, 0, 0], gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igffteps0) 408 if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 409 if (use_padfft==0) then 410 ABI_FREE(gw_gbound) 411 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft)) 412 end if 413 414 nkpt_summed=Kmesh%nbz 415 if (Ep%symchi/=0) then 416 nkpt_summed=Ltg_q%nibz_ltg 417 call littlegroup_print(Ltg_q,std_out,Wfd%prtvol,'COLL') 418 end if 419 ! 420 ! ============================================ 421 ! === Begin big fat loop over transitions ==== 422 ! ============================================ 423 chi0 = czero_gw 424 chi0_head = czero_gw; chi0_lwing = czero_gw; chi0_uwing = czero_gw 425 dim_rtwg=1; if (nspinor==2) dim_rtwg=2 !can reduce size depending on Ep%nI and Ep%nj 426 427 zcut = Ep%zcut 428 zcut = 0.1/Ha_eV 429 write(std_out,*)" using zcut ",zcut*Ha_eV," [eV]" 430 431 ! Loop on spin to calculate $ \chi_{\up,\up} + \chi_{\down,\down} 432 do spin=1,nsppol 433 ! Loop over k-points in the BZ. 434 do ik_bz=1,Kmesh%nbz 435 if (Ep%symchi==1) then 436 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q 437 end if 438 439 ! Get ik_ibz, non-symmorphic phase and symmetries from ik_bz. 440 call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt) 441 tabr_k=ktabr(:,ik_bz) ! Table for rotated FFT points 442 spinrot_kbz(:)=Cryst%spinrot(:,isym_k) 443 nband_k=Wfd%nband(ik_ibz,spin) 444 445 ! Distribute bands. 446 bmask=.FALSE.; bmask(1:nband_k)=.TRUE. ! TODO only bands around EF should be included. 447 call wfd_distribute_bands(Wfd,ik_ibz,spin,my_nband,my_band_list,bmask=bmask) 448 if (my_nband==0) CYCLE 449 450 write(msg,'(2(a,i4),a,i2,a,i3)')' ik = ',ik_bz,' / ',Kmesh%nbz,' spin = ',spin,' done by processor ',Wfd%my_rank 451 call wrtout(std_out,msg,'PERS') 452 453 do lbidx=1,my_nband 454 ! Loop over bands treated by this node. 455 band=my_band_list(lbidx) 456 call wfd_get_ur(Wfd,band,ik_ibz,spin,ur1) 457 458 if (Psps%usepaw==1) then 459 call wfd_get_cprj(Wfd,band,ik_ibz,spin,Cryst,Cprj1_ibz,sorted=.FALSE.) 460 call pawcprj_copy(Cprj1_ibz,Cprj1_bz) 461 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_bz) 462 end if 463 464 deltaf_b1b2 = spin_fact*delta_occ(band,ik_ibz,spin) 465 deltaeGW_b1b2= delta_ene(band,ik_ibz,spin) 466 467 ! Add small imaginary of the Time-Ordered resp function but only for non-zero real omega FIXME What about metals? 468 if (.not.use_tr) then 469 do io=1,Ep%nomega 470 !green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,-one,one_pole) 471 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,GW_TOL_W0,one_pole) 472 end do 473 else 474 do io=1,Ep%nomega ! This expression implements time-reversal even when the input k-mesh breaks it. 475 !green_w(io) = half * g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,-one,two_poles) 476 green_w(io) = half * g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,GW_TOL_W0,two_poles) 477 end do !io 478 end if ! use_tr 479 480 ! FFT of u^*_{b1,k}(r) u_{b2,k}(r). 481 call rho_tw_g(nspinor,Ep%npwe,Wfd%nfft,ndat1,ngfft_gw,1,use_padfft,igffteps0,gw_gbound,& 482 & ur1,itim_k,tabr_k,ph_mkt,spinrot_kbz,& 483 & ur1,itim_k,tabr_k,ph_mkt,spinrot_kbz,& 484 & dim_rtwg,rhotwg) 485 486 if (Psps%usepaw==1) then 487 ! Add PAW onsite contribution, projectors are already in the BZ. 488 call paw_rho_tw_g(Ep%npwe,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,& 489 & Cprj1_bz,Cprj1_bz,Pwij,rhotwg) 490 end if 491 492 ! ==== Adler-Wiser expression, to be consistent here we use the KS eigenvalues (?) ==== 493 call assemblychi0_sym(ik_bz,nspinor,Ep,Ltg_q,green_w,Ep%npwepG0,rhotwg,Gsph_epsG0,chi0) 494 end do !band 495 end do !ik_bz 496 end do !spin 497 498 ! Collect body, heads and wings within comm 499 comm=Wfd%comm 500 do io=1,Ep%nomega 501 call xmpi_sum(chi0(:,:,io),comm,ierr) 502 end do 503 call xmpi_sum(chi0_head,comm,ierr) 504 call xmpi_sum(chi0_lwing,comm,ierr) 505 call xmpi_sum(chi0_uwing,comm,ierr) 506 507 ! Divide by the volume 508 chi0 = chi0 * weight/Cryst%ucvol 509 chi0_head = chi0_head * weight/Cryst%ucvol 510 do io=1,Ep%nomega ! Tensor in the basis of the reciprocal lattice vectors. 511 chi0_head(:,:,io) = MATMUL(chi0_head(:,:,io),Cryst%gmet) * (two_pi**2) 512 end do 513 chi0_lwing = chi0_lwing * weight/Cryst%ucvol 514 chi0_uwing = chi0_uwing * weight/Cryst%ucvol 515 ! 516 ! =============================================== 517 ! ==== Symmetrize chi0 in case of AFM system ==== 518 ! =============================================== 519 ! * Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$. 520 ! * Works only in the case of magnetic group Shubnikov type IV. 521 if (Cryst%use_antiferro) then 522 call symmetrize_afm_chi0(Cryst,Gsph_epsG0,Ltg_q,Ep%npwe,Ep%nomega,chi0,chi0_head,chi0_lwing,chi0_uwing) 523 end if 524 ! 525 ! =================================================== 526 ! ==== Construct heads and wings from the tensor ==== 527 ! =================================================== 528 !do io=1,Ep%nomega 529 ! do ig=2,Ep%npwe 530 ! wng = chi0_uwing(ig,io,:) 531 ! chi0(1,ig,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G") 532 ! wng = chi0_lwing(ig,io,:) 533 ! chi0(ig,1,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G") 534 ! end do 535 ! chq = MATMUL(chi0_head(:,:,io), Ep%qlwl(:,1)) 536 ! chi0(1,1,io) = vdotw(Ep%qlwl(:,1),chq,Cryst%gmet,"G") ! Use user-defined small q 537 !end do 538 !call wfd_barrier(Wfd) 539 540 ! Impose Hermiticity (valid only for zero or purely imaginary frequencies) 541 ! MG what about metals, where we have poles around zero? 542 !if (dtset%gw_eet/=-1) then 543 ! do io=1,Ep%nomega 544 ! if (ABS(REAL(Ep%omega(io)))<0.00001) then 545 ! do ig2=1,Ep%npwe 546 ! do ig1=1,ig2-1 547 ! chi0(ig2,ig1,io)=CONJG(chi0(ig1,ig2,io)) 548 ! end do 549 ! end do 550 ! end if 551 ! end do 552 !end if 553 554 do iomega=1,MIN(Ep%nomega,NOMEGA_PRINTED) 555 write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0_intra(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]' 556 call wrtout(std_out,msg,'COLL') 557 call print_arr(chi0(:,:,iomega),unit=std_out) 558 end do 559 560 ! ===================== 561 ! ==== Free memory ==== 562 ! ===================== 563 ABI_FREE(rhotwg) 564 ABI_FREE(tabr_k) 565 ABI_FREE(ur1) 566 ABI_FREE(gw_gbound) 567 ABI_FREE(ktabr) 568 ABI_FREE(igffteps0) 569 570 ! deallocation for PAW. 571 if (Psps%usepaw==1) then 572 call pawcprj_free(Cprj1_bz) 573 ABI_DT_FREE(Cprj1_bz) 574 call pawcprj_free(Cprj1_ibz) 575 ABI_DT_FREE(Cprj1_ibz) 576 call pawpwij_free(Pwij) 577 ABI_DT_FREE(Pwij) 578 end if 579 580 call littlegroup_free(Ltg_q) 581 call kmesh_free(Kmesh) 582 583 DBG_EXIT("COLL") 584 585 end subroutine chi0q0_intraband