TABLE OF CONTENTS
ABINIT/cchi0 [ Functions ]
NAME
cchi0
FUNCTION
Main calculation of the independent-particle susceptibility chi0 for qpoint!=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. valence states are stored in Wfs_val and only resonant transitions are calculated (time reversal is assumed) Dtset <type(dataset_type)>=all input variables in this dataset Cryst<crystal_t>= data type gathering info on symmetries and unit cell %natom=number of atoms %nsym=number of symmetries %xred(3,natom)=reduced coordinated of atoms %typat(natom)=type of each atom %rprimd(3,3)=dimensional primitive translations in real space (bohr) %timrev= 2 if time reversal can be used, 1 otherwise qpoint(3)=reciprocal space coordinates of the q wavevector Ep<type(em1params_t_type)>= Parameters related to the calculation of the inverse dielectric matrix. %nbnds=number of bands summed over %npwe=number of planewaves for the irreducible polarizability X^0_GGp %npwvec=maximum number of G vectors used to define the dimension of some arrays e.g igfft %nsppol=1 for unpolarized, 2 for spin-polarized %nomega=total number of frequencies in X^0 (both real and imaginary) %nomegasf=number of real frequencies used to sample the imaginary part of X^0 (spectral method) %spmeth=1 if we use the spectral method, 0 for standard Adler-Wiser expression %spsmear=gaussian broadening used to approximate the delta distribution %zcut=small imaginary shift to avoid poles in X^0 Psps <type(pseudopotential_type)>=variables related to pseudopotentials Kmesh <kmesh_t>= datatype gathering parameters related to the k-point sampling %nibz=number of k-points in the IBZ %nbz=number of k-points in the BZ %bz(3,nbz)=reduced coordinates for k-points in the full Brillouin zone %ibz(3,nibz)=reduced coordinates for k-points in the irreducible wedge %tab(nbz)=mapping between a kpt in the BZ (array bz) and the irred point in the array ibz %tabi(nbz)= -1 if inversion is needed to obtain this particular kpt in the BZ, 1 means identity %tabo(nbz)= for each point in the BZ, the index of the symmetry operation S in reciprocal space which rotates k_IBZ onto \pm k_BZ (depending on tabi) %tabp(nbz)= For each k_BZ, it gives the phase factors associated to non-symmorphic operations, i.e e^{-i 2 \pi k_IBZ \cdot R{^-1}t} == e{-i 2\pi k_BZ cdot t} where : \transpose R{-1}=S and (S k_IBZ) = \pm k_BZ (depending on ktabi) %tabr(nfftot,nbz) For each point r on the real mesh and for each k-point in the BZ, tabr gives the index of (R^-1 (r-t)) in the FFT array where R=\transpose S^{-1} and k_BZ=S k_IBZ. t is the fractional translation associated to R Gsph_epsG0<gsphere_t data type> The G-sphere used to describe chi0/eps. (including umklapp G0 vectors) %ng=number of G vectors for chi0 %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G %phmGt(Ep%npwe,nsym)=phase factors e^{-iG \cdot t} needed to symmetrize oscillator matrix elements and epsilon %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1) %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$). nbvw=number of bands in the arrays wfrv ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths (see input variable) nfftot_gw=Total number of points in the GW FFT grid Ltg_q<Little group>=Data type gathering information on the little group of the q-points. Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data Pawang<pawang_type> angular mesh discretization and related data: 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 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 at wavevector qpoint and each frequeny defined by Ep%omega and Ep%nomega
PARENTS
screening
CHILDREN
accumulate_chi0sumrule,approxdelta,assemblychi0_sym,assemblychi0sf calc_wfwfg,chi0_bbp_mask,completechi0_deltapart,cwtime,flush_unit get_bz_diff,get_bz_item,get_gftt,gsph_fft_tabs,gsph_free,gsph_in_fftbox hilbert_transform,littlegroup_print,make_transitions,paw_cross_rho_tw_g paw_rho_tw_g,paw_symcprj,pawcprj_alloc,pawcprj_free,pawpwij_free pawpwij_init,read_plowannier,rho_tw_g,setup_spectral symmetrize_afm_chi0,timab,wfd_change_ngfft,wfd_distribute_kb_kpbp wfd_get_cprj,wfd_get_ur,wfd_paw_get_aeur,wrtout,xmpi_sum
SOURCE
93 #if defined HAVE_CONFIG_H 94 #include "config.h" 95 #endif 96 97 #include "abi_common.h" 98 99 subroutine cchi0(use_tr,Dtset,Cryst,qpoint,Ep,Psps,Kmesh,QP_BSt,Gsph_epsG0,& 100 & Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,nbvw,ngfft_gw,nfftot_gw,ngfftf,nfftf_tot,& 101 & chi0,ktabr,ktabrf,Ltg_q,chi0_sumrule,Wfd,Wfdf) 102 103 use defs_basis 104 use defs_datatypes 105 use defs_abitypes 106 use m_xmpi 107 use m_blas 108 use m_errors 109 use m_profiling_abi 110 use m_time 111 112 use m_gwdefs, only : GW_TOL_DOCC, GW_TOL_W0, czero_gw, em1params_t, g0g0w 113 use m_numeric_tools, only : imin_loc 114 use m_geometry, only : normv 115 use m_io_tools, only : flush_unit 116 use m_crystal, only : crystal_t 117 use m_bz_mesh, only : kmesh_t, get_BZ_item, get_BZ_diff, littlegroup_t, littlegroup_print 118 use m_gsphere, only : gsphere_t, gsph_fft_tabs, gsph_free, gsph_in_fftbox 119 use m_fft_mesh, only : get_gftt 120 use m_wfd, only : wfd_get_ur, wfd_t, wfd_distribute_kb_kpbp, wfd_get_cprj, wfd_barrier, wfd_change_ngfft,& 121 & wfd_paw_get_aeur, wfd_sym_ur 122 use m_oscillators, only : rho_tw_g, calc_wfwfg 123 use m_chi0, only : hilbert_transform, setup_spectral, assemblychi0_sym, assemblychi0sf, symmetrize_afm_chi0,& 124 & approxdelta, completechi0_deltapart, accumulate_chi0sumrule, make_transitions, chi0_bbp_mask 125 126 use m_pawang, only : pawang_type 127 use m_pawtab, only : pawtab_type 128 use m_pawfgrtab, only : pawfgrtab_type 129 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 130 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g 131 use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t 132 133 !This section has been created automatically by the script Abilint (TD). 134 !Do not modify the following lines by hand. 135 #undef ABI_FUNC 136 #define ABI_FUNC 'cchi0' 137 use interfaces_14_hidewrite 138 use interfaces_18_timing 139 use interfaces_65_paw 140 use interfaces_70_gw, except_this_one => cchi0 141 !End of the abilint section 142 143 implicit none 144 145 !Arguments ------------------------------------ 146 !scalars 147 integer,intent(in) :: nbvw,nfftot_gw,nfftf_tot 148 logical,intent(in) :: use_tr 149 type(ebands_t),target,intent(in) :: QP_BSt 150 type(kmesh_t),intent(in) :: Kmesh 151 type(crystal_t),intent(in) :: Cryst 152 type(Dataset_type),intent(in) :: Dtset 153 type(em1params_t),intent(in) :: Ep 154 type(gsphere_t),intent(in) :: Gsph_epsG0 155 type(littlegroup_t),intent(in) :: Ltg_q 156 type(Pawang_type),intent(in) :: Pawang 157 type(Pseudopotential_type),intent(in) :: Psps 158 type(wfd_t),target,intent(inout) :: Wfd,Wfdf 159 !arrays 160 integer,intent(in) :: ktabr(nfftot_gw,Kmesh%nbz),ktabrf(nfftf_tot*Dtset%pawcross,Kmesh%nbz) 161 integer,intent(in) :: ngfft_gw(18),ngfftf(18) 162 real(dp),intent(in) :: qpoint(3) 163 real(dp),intent(out) :: chi0_sumrule(Ep%npwe) 164 complex(gwpc),intent(out) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega) 165 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 166 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 167 type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw) 168 type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom) 169 170 !Local variables ------------------------------ 171 !scalars 172 integer,parameter :: tim_fourdp1=1,two_poles=2,one_pole=1,ndat1=1 173 integer :: bandinf,bandsup,dim_rtwg,band1,band2,ierr 174 integer :: ig1,ig2,iat,ik_bz,ik_ibz,ikmq_bz,ikmq_ibz 175 integer :: io,iomegal,iomegar,ispinor1,ispinor2,isym_k,itypatcor,nfft 176 integer :: isym_kmq,itim_k,itim_kmq,m1,m2,my_wl,my_wr,size_chi0 177 integer :: nfound,nkpt_summed,nspinor,nsppol,mband 178 integer :: comm,gw_mgfft,use_padfft,gw_fftalga,lcor,mgfftf,use_padfftf 179 integer :: my_nbbp,my_nbbpks,spin,nbmax,dummy 180 real(dp) :: cpu_time,wall_time,gflops 181 real(dp) :: deltaeGW_b1kmq_b2k,deltaeGW_enhigh_b2k,deltaf_b1kmq_b2k 182 real(dp) :: e_b1_kmq,en_high,fac,fac2,fac3,f_b1_kmq,factor,max_rest,min_rest,my_max_rest 183 real(dp) :: my_min_rest,numerator,spin_fact,weight,wl,wr 184 real(dp) :: gw_gsq,memreq 185 complex(dpc) :: ph_mkmqt,ph_mkt 186 complex(gwpc) :: local_czero_gw 187 logical :: qzero,isirred_k,isirred_kmq,luwindow 188 character(len=500) :: msg,allup 189 type(gsphere_t) :: Gsph_FFT 190 !arrays 191 integer :: G0(3),umklp_k(3),umklp_kmq(3) 192 integer :: ucrpa_bands(2) 193 integer :: wtk_ltg(Kmesh%nbz),got(Wfd%nproc) 194 integer,allocatable :: tabr_k(:),tabr_kmq(:),tabrf_k(:),tabrf_kmq(:) 195 integer,allocatable :: igfftepsG0(:),gspfft_igfft(:),igfftepsG0f(:) 196 integer,allocatable :: gw_gfft(:,:),gw_gbound(:,:),dummy_gbound(:,:),gboundf(:,:) 197 integer,allocatable :: bbp_ks_distrb(:,:,:,:) 198 real(dp) :: kbz(3),kmq_bz(3),spinrot_k(4),spinrot_kmq(4),q0(3),tsec(2) 199 real(dp),ABI_CONTIGUOUS pointer :: qp_energy(:,:,:),qp_occ(:,:,:) 200 real(dp),allocatable :: omegasf(:) 201 complex(dpc),allocatable :: green_enhigh_w(:),green_w(:),kkweight(:,:) 202 complex(gwpc),allocatable :: sf_chi0(:,:,:),rhotwg(:) 203 complex(gwpc),allocatable :: ur1_kmq_ibz(:),ur2_k_ibz(:),wfwfg(:) 204 complex(gwpc),allocatable :: usr1_kmq(:),ur2_k(:) 205 complex(gwpc),allocatable :: ur_ae1(:),ur_ae_onsite1(:),ur_ps_onsite1(:) 206 complex(gwpc),allocatable :: ur_ae2(:),ur_ae_onsite2(:),ur_ps_onsite2(:) 207 complex(dpc), allocatable :: coeffW_BZ(:,:,:,:,:,:) 208 logical,allocatable :: bbp_mask(:,:) 209 type(pawcprj_type),allocatable :: Cprj1_kmq(:,:),Cprj2_k(:,:) 210 type(pawpwij_t),allocatable :: Pwij(:),Pwij_fft(:) 211 !************************************************************************ 212 213 DBG_ENTER("COLL") 214 215 call timab(331,1,tsec) ! cchi0 216 call cwtime(cpu_time,wall_time,gflops,"start") 217 218 nsppol = Wfd%nsppol; nspinor = Wfd%nspinor 219 ucrpa_bands(1)=dtset%ucrpa_bands(1) 220 ucrpa_bands(2)=dtset%ucrpa_bands(2) 221 luwindow=.false. 222 if(abs(dtset%ucrpa_window(1)+1_dp)>tol8.and.(abs(dtset%ucrpa_window(2)+1_dp)>tol8)) then 223 luwindow=.true. 224 endif 225 ! write(6,*)"ucrpa_bands",ucrpa_bands 226 ! write(6,*)"ucrpa_window",dtset%ucrpa_window 227 ! write(6,*)"luwindow",luwindow 228 229 ! For cRPA calculation of U: read forlb.ovlp 230 if(dtset%ucrpa>=1) then 231 call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor,Kmesh,lcor,luwindow,& 232 & nspinor,nsppol,pawang,dtset%prtvol,ucrpa_bands) 233 endif 234 ! End of reading forlb.ovlp 235 236 if ( ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,ngfft_gw) 237 gw_mgfft = MAXVAL(ngfft_gw(1:3)) 238 gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10) 239 240 if (Dtset%pawcross==1) mgfftf = MAXVAL(ngfftf(1:3)) 241 242 ! == Copy some values === 243 comm = Wfd%comm 244 mband = Wfd%mband 245 nfft = Wfd%nfft 246 ABI_CHECK(Wfd%nfftot==nfftot_gw,"Wrong nfftot_gw") 247 248 dim_rtwg=1 !; if (nspinor==2) dim_rtwg=2 ! can reduce size depending on Ep%nI and Ep%nj 249 size_chi0 = Ep%npwe*Ep%nI*Ep%npwe*Ep%nJ*Ep%nomega 250 251 qp_energy => QP_BSt%eig; qp_occ => QP_BSt%occ 252 253 ! Initialize the completeness correction 254 if (Ep%gwcomp==1) then 255 en_high=MAXVAL(qp_energy(Ep%nbnds,:,:)) + Ep%gwencomp 256 write(msg,'(a,f8.2,a)')' Using completeness correction with the energy ',en_high*Ha_eV,' [eV]' 257 call wrtout(std_out,msg,'COLL') 258 259 ! Allocation of wfwfg and green_enhigh_w moved inside openmp loop 260 ! Init the largest G-sphere contained in the FFT box for the wavefunctions. 261 call gsph_in_fftbox(Gsph_FFT,Cryst,Wfd%ngfft) 262 263 !call print_gsphere(Gsph_FFT,unit=std_out,prtvol=10) 264 265 ABI_MALLOC(gspfft_igfft,(Gsph_FFT%ng)) 266 ABI_MALLOC(dummy_gbound,(2*gw_mgfft+8,2)) 267 268 ! Mapping between G-sphere and FFT box. 269 call gsph_fft_tabs(Gsph_FFT,(/0,0,0/),Wfd%mgfft,Wfd%ngfft,dummy,dummy_gbound,gspfft_igfft) 270 ABI_FREE(dummy_gbound) 271 272 if (Psps%usepaw==1) then ! * Prepare the onsite contributions on the GW FFT mesh. 273 ABI_MALLOC(gw_gfft,(3,nfft)) 274 q0=zero 275 call get_gftt(ngfft_gw,q0,Cryst%gmet,gw_gsq,gw_gfft) ! Get the set of plane waves in the FFT Box. 276 ABI_DT_MALLOC(Pwij_fft,(Psps%ntypat)) 277 call pawpwij_init(Pwij_fft,nfft,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 278 end if 279 end if 280 281 ! Setup weights (2 for spin unpolarized sistem, 1 for polarized). 282 ! spin_fact is used to normalize the occupation factors to one. Consider also the AFM case. 283 SELECT CASE (nsppol) 284 CASE (1) 285 weight=two/Kmesh%nbz; spin_fact=half 286 if (Wfd%nspden==2) then 287 weight=one/Kmesh%nbz; spin_fact=half 288 end if 289 if (nspinor==2) then 290 weight=one/Kmesh%nbz; spin_fact=one 291 end if 292 CASE (2) 293 weight=one/Kmesh%nbz; spin_fact=one 294 CASE DEFAULT 295 MSG_BUG("Wrong nsppol") 296 END SELECT 297 298 ! Weight for points in the IBZ_q. 299 wtk_ltg(:)=1 300 if (Ep%symchi==1) then 301 do ik_bz=1,Ltg_q%nbz 302 wtk_ltg(ik_bz)=0 303 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only k-points in the IBZ_q. 304 wtk_ltg(ik_bz)=SUM(Ltg_q%wtksym(:,:,ik_bz)) 305 end do 306 end if 307 308 write(msg,'(a,i2,2a,i2)')& 309 & ' Using spectral method for the imaginary part = ',Ep%spmeth,ch10,& 310 & ' Using symmetries to sum only over the IBZ_q = ',Ep%symchi 311 call wrtout(std_out,msg,'COLL') 312 313 if (use_tr) then 314 ! Special care has to be taken in metals and/or spin dependent systems 315 ! as Wfs_val might contain unoccupied states. 316 call wrtout(std_out,' Using faster algorithm based on time reversal symmetry.','COLL') 317 else 318 call wrtout(std_out,' Using slow algorithm without time reversal symmetry.') 319 end if 320 321 ! TODO this table can be calculated for each k-point 322 my_nbbpks=0; allup="All"; got=0 323 ABI_MALLOC(bbp_ks_distrb,(mband,mband,Kmesh%nbz,nsppol)) 324 ABI_MALLOC(bbp_mask,(mband,mband)) 325 326 do spin=1,nsppol 327 do ik_bz=1,Kmesh%nbz 328 if (Ep%symchi==1) then 329 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q 330 end if 331 332 ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz. 333 call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k) 334 335 ! Get index of k-q in the BZ, stop if not found as the weight=one/nkbz is not correct. 336 call get_BZ_diff(Kmesh,kbz,qpoint,ikmq_bz,g0,nfound) 337 ABI_CHECK(nfound==1,"Check kmesh") 338 339 ! Get ikmq_ibz, non-symmorphic phase, ph_mkmqt, and symmetries from ikmq_bz. 340 call get_BZ_item(Kmesh,ikmq_bz,kmq_bz,ikmq_ibz,isym_kmq,itim_kmq) 341 342 call chi0_bbp_mask(Ep,use_tr,QP_BSt,mband,ikmq_ibz,ik_ibz,spin,spin_fact,bbp_mask) 343 344 call wfd_distribute_kb_kpbp(Wfd,ikmq_ibz,ik_ibz,spin,allup,my_nbbp,bbp_ks_distrb(:,:,ik_bz,spin),got,bbp_mask) 345 my_nbbpks = my_nbbpks + my_nbbp 346 end do 347 end do 348 349 ABI_FREE(bbp_mask) 350 351 write(msg,'(a,i0,a)')" Will sum ",my_nbbpks," (b,b',k,s) states in chi0." 352 call wrtout(std_out,msg,'PERS') 353 354 if (Psps%usepaw==1) then 355 ABI_DT_MALLOC(Pwij,(Psps%ntypat)) 356 call pawpwij_init(Pwij,Ep%npwepG0,qpoint,Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 357 ! Allocate statements moved to inside openmp loop 358 end if 359 360 SELECT CASE (Ep%spmeth) 361 CASE (0) 362 call wrtout(std_out,' Calculating chi0(q,omega,G,G")','COLL') 363 ! Allocation of green_w moved inside openmp loop 364 365 CASE (1, 2) 366 call wrtout(std_out,' Calculating Im chi0(q,omega,G,G")','COLL') 367 368 ! Find Max and min resonant transitions for this q, report also treated by this proc. 369 call make_transitions(Wfd,1,Ep%nbnds,nbvw,nsppol,Ep%symchi,Cryst%timrev,GW_TOL_DOCC,& 370 & max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,qp_energy,qp_occ,qpoint,bbp_ks_distrb) 371 ! 372 ! Calculate frequency dependent weights for Hilbert transform. 373 ABI_MALLOC(omegasf,(Ep%nomegasf)) 374 ABI_MALLOC(kkweight,(Ep%nomegasf,Ep%nomega)) 375 !my_wl=1; my_wr=Ep%nomegasf 376 call setup_spectral(Ep%nomega,Ep%omega,Ep%nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,& 377 & 0,Ep%zcut,zero,my_wl,my_wr,kkweight) 378 379 if (.not.use_tr) then 380 MSG_BUG('spectral method requires time-reversal') 381 end if 382 383 memreq = two*gwpc*Ep%npwe**2*(my_wr-my_wl+1)*b2Gb 384 write(msg,'(a,f10.4,a)')' memory required per spectral point: ',two*gwpc*Ep%npwe**2*b2Mb,' [Mb]' 385 call wrtout(std_out,msg,'PERS') 386 write(msg,'(a,f10.4,a)')' memory required by sf_chi0: ',memreq,' [Gb]' 387 call wrtout(std_out,msg,'PERS') 388 if (memreq > two) then 389 MSG_WARNING(' Memory required for sf_chi0 is larger than 2.0 Gb!') 390 end if 391 ABI_STAT_MALLOC(sf_chi0,(Ep%npwe,Ep%npwe,my_wl:my_wr), ierr) 392 ABI_CHECK(ierr==0, 'out-of-memory in sf_chi0') 393 sf_chi0=czero_gw 394 395 CASE DEFAULT 396 MSG_BUG("Wrong spmeth") 397 END SELECT 398 399 nkpt_summed=Kmesh%nbz 400 if (Ep%symchi==1) then 401 nkpt_summed=Ltg_q%nibz_ltg 402 call littlegroup_print(Ltg_q,std_out,Dtset%prtvol,'COLL') 403 end if 404 405 write(msg,'(a,i6,a)')' Calculation status : ',nkpt_summed,' to be completed ' 406 call wrtout(std_out,msg,'COLL') 407 408 ! ============================================ 409 ! === Begin big fat loop over transitions === 410 ! ============================================ 411 chi0=czero_gw; chi0_sumrule=zero 412 413 ! === Loop on spin to calculate trace $\chi_{up,up}+\chi_{down,down}$ === 414 ! Only $\chi_{up,up} for AFM. 415 do spin=1,nsppol 416 if (ALL(bbp_ks_distrb(:,:,:,spin) /= Wfd%my_rank)) CYCLE 417 418 ! Allocation of arrays that are private to loop 419 if (Ep%gwcomp==1) then 420 ABI_MALLOC(wfwfg,(nfft*nspinor**2)) 421 end if 422 if (Ep%gwcomp==1) then 423 ABI_MALLOC(green_enhigh_w,(Ep%nomega)) 424 end if 425 if (Ep%spmeth==0) then 426 ABI_MALLOC(green_w,(Ep%nomega)) 427 end if 428 if (Psps%usepaw==1) then 429 ABI_DT_MALLOC(Cprj2_k ,(Cryst%natom,nspinor)) 430 call pawcprj_alloc(Cprj2_k, 0,Wfd%nlmn_atm) 431 ABI_DT_MALLOC(Cprj1_kmq,(Cryst%natom,nspinor)) 432 call pawcprj_alloc(Cprj1_kmq,0,Wfd%nlmn_atm) 433 if (Dtset%pawcross==1) then 434 ABI_MALLOC(ur_ae1,(nfftf_tot*nspinor)) 435 ABI_MALLOC(ur_ae_onsite1,(nfftf_tot*nspinor)) 436 ABI_MALLOC(ur_ps_onsite1,(nfftf_tot*nspinor)) 437 ABI_MALLOC(ur_ae2,(nfftf_tot*nspinor)) 438 ABI_MALLOC(ur_ae_onsite2,(nfftf_tot*nspinor)) 439 ABI_MALLOC(ur_ps_onsite2,(nfftf_tot*nspinor)) 440 ABI_MALLOC(igfftepsG0f,(Ep%npwepG0)) 441 ABI_MALLOC(tabrf_k,(nfftf_tot)) 442 ABI_MALLOC(tabrf_kmq,(nfftf_tot)) 443 end if 444 end if 445 446 ABI_MALLOC(rhotwg,(Ep%npwepG0*dim_rtwg)) 447 ABI_MALLOC(tabr_k,(nfft)) 448 ABI_MALLOC(tabr_kmq,(nfft)) 449 ABI_MALLOC(ur1_kmq_ibz,(nfft*nspinor)) 450 ABI_MALLOC(ur2_k_ibz,(nfft*nspinor)) 451 ABI_MALLOC(usr1_kmq,(nfft*nspinor)) 452 ABI_MALLOC(ur2_k, (nfft*nspinor)) 453 ABI_MALLOC(igfftepsG0,(Ep%npwepG0)) 454 455 ! Loop over k-points in the BZ. 456 do ik_bz=1,Kmesh%nbz 457 458 if (Ep%symchi==1) then 459 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q 460 end if 461 462 if (ALL(bbp_ks_distrb(:,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE 463 464 write(msg,'(2(a,i4),a,i2,a,i3)')' ik= ',ik_bz,'/',Kmesh%nbz,' spin= ',spin,' done by mpi rank:',Wfd%my_rank 465 call wrtout(std_out,msg,'PERS') 466 467 ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz. 468 call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt,umklp_k,isirred_k) 469 470 call get_BZ_diff(Kmesh,kbz,qpoint,ikmq_bz,G0,nfound) 471 if (nfound==0) then 472 MSG_ERROR("Cannot find kbz - qpoint in Kmesh") 473 end if 474 475 ! Get ikmq_ibz, non-symmorphic phase, ph_mkmqt, and symmetries from ikmq_bz. 476 call get_BZ_item(Kmesh,ikmq_bz,kmq_bz,ikmq_ibz,isym_kmq,itim_kmq,ph_mkmqt,umklp_kmq,isirred_kmq) 477 478 !BEGIN DEBUG 479 !if (ANY(umklp_k /=0)) then 480 ! write(msg,'(a,3i2)')" umklp_k /= 0 ",umklp_k 481 ! MSG_ERROR(msg) 482 !end if 483 !if (ANY( g0 /= -umklp_kmq + umklp_k) ) then 484 !if (ANY( g0 /= -umklp_kmq ) ) then 485 ! write(msg,'(a,3(1x,3i2))')" g0 /= -umklp_kmq + umklp_k ",g0, umklp_kmq, umklp_k 486 ! MSG_ERROR(msg) 487 !end if 488 !g0 = -umklp_k + umklp_kmq 489 !g0 = +umklp_k - umklp_kmq 490 !if (ANY (ABS(g0) > Ep%mg0) ) then 491 ! write(msg,'(a,6(1x,i0))')" ABS(g0) > Ep%mg0 ",g0,Ep%mg0 492 ! MSG_ERROR(msg) 493 !end if 494 !END DEBUG 495 496 ! Copy tables for rotated FFT points 497 tabr_k(:) =ktabr(:,ik_bz) 498 spinrot_k(:)=Cryst%spinrot(:,isym_k) 499 500 tabr_kmq(:)=ktabr(:,ikmq_bz) 501 spinrot_kmq(:)=Cryst%spinrot(:,isym_kmq) 502 503 if (Dtset%pawcross==1) then 504 tabrf_k(:) =ktabrf(:,ik_bz) 505 tabrf_kmq(:)=ktabrf(:,ikmq_bz) 506 end if 507 ! 508 ! Tables for the FFT of the oscillators. 509 ! a) FFT index of G-G0. 510 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 511 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2)) 512 call gsph_fft_tabs(Gsph_epsG0,g0,gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igfftepsG0) 513 if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 514 if (use_padfft==0) then 515 ABI_FREE(gw_gbound) 516 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft)) 517 end if 518 519 if (Dtset%pawcross==1) then 520 ABI_MALLOC(gboundf,(2*mgfftf+8,2)) 521 call gsph_fft_tabs(Gsph_epsG0,g0,mgfftf,ngfftf,use_padfftf,gboundf,igfftepsG0f) 522 if (ANY(gw_fftalga == [2, 4])) use_padfftf=0 523 if (use_padfftf==0) then 524 ABI_FREE(gboundf) 525 ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf)) 526 end if 527 end if 528 529 nbmax=Ep%nbnds 530 do band1=1,nbmax ! Loop over "conduction" states. 531 if (ALL(bbp_ks_distrb(band1,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE 532 533 call wfd_get_ur(Wfd,band1,ikmq_ibz,spin,ur1_kmq_ibz) 534 535 if (Psps%usepaw==1) then 536 call wfd_get_cprj(Wfd,band1,ikmq_ibz,spin,Cryst,Cprj1_kmq,sorted=.FALSE.) 537 call paw_symcprj(ikmq_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_kmq) 538 if (Dtset%pawcross==1) then 539 call wfd_paw_get_aeur(Wfdf,band1,ikmq_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae1,ur_ae_onsite1,ur_ps_onsite1) 540 end if 541 end if 542 543 e_b1_kmq=qp_energy(band1,ikmq_ibz,spin) 544 f_b1_kmq= qp_occ(band1,ikmq_ibz,spin) 545 546 do band2=1,nbmax ! Loop over "valence" states. 547 !debug if (.not.luwindow.AND.dtset%ucrpa==1.AND.band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)& 548 !debug& .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) CYClE 549 !write(6,*) "ik,band1,band2",ik_bz,band1,band2 550 if (luwindow.AND.dtset%ucrpa==1.AND.((QP_Bst%eig(band1,ik_ibz ,spin)-QP_Bst%fermie)<=dtset%ucrpa_window(2)) & 551 & .AND.((QP_Bst%eig(band1,ik_ibz ,spin)-QP_Bst%fermie)>=dtset%ucrpa_window(1)) & 552 & .AND.((QP_Bst%eig(band2,ikmq_ibz,spin)-QP_Bst%fermie)<=dtset%ucrpa_window(2)) & 553 & .AND.((QP_Bst%eig(band2,ikmq_ibz,spin)-QP_Bst%fermie)>=dtset%ucrpa_window(1))) CYCLE 554 555 if (bbp_ks_distrb(band1,band2,ik_bz,spin) /= Wfd%my_rank) CYCLE 556 557 deltaf_b1kmq_b2k=spin_fact*(f_b1_kmq-qp_occ(band2,ik_ibz,spin)) 558 559 if (Ep%gwcomp==0) then ! Skip negligible transitions. 560 if (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC) CYCLE 561 else 562 ! When the completeness correction is used, 563 ! we need to also consider transitions with vanishing deltaf 564 !if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC) CYCLE 565 ! 566 ! Rangel This is to compute chi correctly when using the extrapolar method 567 if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC .and. (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC .or. band1<band2)) CYCLE 568 end if 569 570 deltaeGW_b1kmq_b2k=e_b1_kmq-qp_energy(band2,ik_ibz,spin) 571 572 call wfd_get_ur(Wfd,band2,ik_ibz,spin,ur2_k_ibz) 573 574 if (Psps%usepaw==1) then 575 call wfd_get_cprj(Wfd,band2,ik_ibz,spin,Cryst,Cprj2_k,sorted=.FALSE.) 576 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj2_k) 577 if (Dtset%pawcross==1) then 578 call wfd_paw_get_aeur(Wfdf,band2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae2,ur_ae_onsite2,ur_ps_onsite2) 579 end if 580 end if 581 582 SELECT CASE (Ep%spmeth) 583 CASE (0) 584 ! Standard Adler-Wiser expression. 585 ! Add the small imaginary of the Time-Ordered RF only for non-zero real omega ! FIXME What about metals? 586 if (.not.use_tr) then 587 ! Have to sum over all possible resonant and anti-resonant transitions. 588 do io=1,Ep%nomega 589 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,one_pole) 590 end do 591 592 else 593 if (Ep%gwcomp==0) then ! cannot be completely skipped in case of completeness correction 594 if (band1<band2) CYCLE ! Here we GAIN a factor ~2 595 end if 596 597 do io=1,Ep%nomega 598 !Rangel: In metals, the intra-band transitions term does not contain the antiresonant part 599 !green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0) 600 if (band1==band2) then 601 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,one_pole) 602 else 603 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,two_poles) 604 end if 605 606 if (Ep%gwcomp==1) then ! Calculate the completeness correction 607 numerator= -spin_fact*qp_occ(band2,ik_ibz,spin) 608 deltaeGW_enhigh_b2k = en_high-qp_energy(band2,ik_ibz,spin) 609 610 if (REAL(Ep%omega(io))<GW_TOL_W0) then ! Completeness correction is NOT valid for real frequencies 611 green_enhigh_w(io) = g0g0w(Ep%omega(io),numerator,deltaeGW_enhigh_b2k,Ep%zcut,GW_TOL_W0,two_poles) 612 else 613 green_enhigh_w(io) = local_czero_gw 614 end if 615 ! 616 !Rangel Correction for metals 617 !if (deltaf_b1kmq_b2k<0.d0) then 618 if (band1>=band2 .and. ABS(deltaf_b1kmq_b2k) > GW_TOL_DOCC ) then 619 green_w(io)= green_w(io) - green_enhigh_w(io) 620 else ! Disregard green_w, since it is already accounted for through the time-reversal 621 green_w(io)= - green_enhigh_w(io) 622 end if 623 end if !gwcomp==1 624 end do !io 625 626 if (Ep%gwcomp==1.and.band1==band2) then 627 ! Add the "delta part" of the extrapolar method. TODO doesnt work for spinor 628 call calc_wfwfg(tabr_k,itim_k,spinrot_k,nfft,nspinor,ngfft_gw,ur2_k_ibz,ur2_k_ibz,wfwfg) 629 630 if (Psps%usepaw==1) then 631 call paw_rho_tw_g(nfft,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,gw_gfft,& 632 & Cprj2_k,Cprj2_k,Pwij_fft,wfwfg) 633 634 ! Add PAW cross term 635 if (Dtset%pawcross==1) then 636 call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,& 637 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_kmq,tabrf_kmq,ph_mkmqt,spinrot_kmq,& 638 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k ,tabrf_k ,ph_mkt ,spinrot_k,dim_rtwg,wfwfg) 639 end if 640 end if 641 642 qzero=.FALSE. 643 call completechi0_deltapart(ik_bz,qzero,Ep%symchi,Ep%npwe,Gsph_FFT%ng,Ep%nomega,nspinor,& 644 & nfft,ngfft_gw,gspfft_igfft,gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0) 645 646 end if 647 end if ! use_tr 648 649 CASE (1, 2) 650 ! Spectral method, WARNING time-reversal here is always assumed! 651 if (deltaeGW_b1kmq_b2k<0) CYCLE 652 call approxdelta(Ep%nomegasf,omegasf,deltaeGW_b1kmq_b2k,Ep%spsmear,iomegal,iomegar,wl,wr,Ep%spmeth) 653 END SELECT 654 655 ! Form rho-twiddle(r)=u^*_{b1,kmq_bz}(r) u_{b2,kbz}(r) and its FFT transform. 656 call rho_tw_g(nspinor,Ep%npwepG0,nfft,ndat1,ngfft_gw,1,use_padfft,igfftepsG0,gw_gbound,& 657 & ur1_kmq_ibz,itim_kmq,tabr_kmq,ph_mkmqt,spinrot_kmq,& 658 & ur2_k_ibz, itim_k ,tabr_k ,ph_mkt ,spinrot_k,dim_rtwg,rhotwg) 659 660 if (Psps%usepaw==1) then 661 ! Add PAW on-site contribution, projectors are already in the BZ. 662 call paw_rho_tw_g(Ep%npwepG0,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,& 663 & Cprj1_kmq,Cprj2_k,Pwij,rhotwg) 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_ae1,ur_ae_onsite1,ur_ps_onsite1,itim_kmq,tabrf_kmq,ph_mkmqt,spinrot_kmq,& 669 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k ,tabrf_k ,ph_mkt ,spinrot_k,dim_rtwg,rhotwg) 670 end if 671 end if 672 673 SELECT CASE (Ep%spmeth) 674 675 CASE (0) ! Adler-Wiser. 676 !debug if(dtset%ucrpa==2) then 677 if(dtset%ucrpa>=1.and..not.luwindow) then 678 fac=one 679 fac2=one 680 fac3=one 681 m1=-1 682 m2=-1 683 call flush_unit(std_out) 684 if(dtset%ucrpa<=2) then 685 if ( band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)& 686 & .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then 687 do iat=1, cryst%nattyp(itypatcor) 688 do ispinor1=1,nspinor 689 do ispinor2=1,nspinor 690 do m1=1,2*lcor+1 691 do m2=1,2*lcor+1 692 fac=fac - real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 693 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* & 694 & coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor2,m2)*& 695 & conjg(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor2,m2))) 696 enddo 697 enddo 698 enddo 699 enddo 700 enddo 701 if(dtset%ucrpa==1) fac=zero 702 endif 703 else if (dtset%ucrpa>=3) then 704 if (band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)) then 705 do iat=1, cryst%nattyp(itypatcor) 706 do ispinor1=1,nspinor 707 do m1=1,2*lcor+1 708 fac2=fac2-real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 709 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))) 710 enddo 711 enddo 712 enddo 713 if(dtset%ucrpa==4) fac2=zero 714 endif 715 if (band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then 716 do iat=1, cryst%nattyp(itypatcor) 717 do ispinor1=1,nspinor 718 do m1=1,2*lcor+1 719 fac3=fac3-real(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor1,m1)*& 720 & conjg(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor1,m1))) 721 enddo 722 enddo 723 enddo 724 if(dtset%ucrpa==4) fac3=zero 725 endif 726 fac=real(fac2*fac3) 727 endif 728 729 ! if(dtset%prtvol>=10) write(6,'(6i3,e15.5,a)') ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q=/0" 730 ! if(dtset%prtvol>=10.and.abs(fac-one)>0.00001) & 731 !& write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q/=0" 732 ! if(dtset%prtvol>=10) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q/=0" 733 green_w=green_w*fac 734 endif 735 736 call assemblychi0_sym(ik_bz,nspinor,Ep,Ltg_q,green_w,Ep%npwepG0,rhotwg,Gsph_epsG0,chi0) 737 738 CASE (1, 2) 739 ! Spectral method (not yet adapted for nspinor=2) 740 call assemblychi0sf(ik_bz,Ep%symchi,Ltg_q,Ep%npwepG0,Ep%npwe,rhotwg,Gsph_epsG0,& 741 & deltaf_b1kmq_b2k,my_wl,iomegal,wl,my_wr,iomegar,wr,Ep%nomegasf,sf_chi0) 742 743 CASE DEFAULT 744 MSG_BUG("Wrong spmeth") 745 END SELECT 746 747 ! Accumulating the sum rule on chi0. Eq. (5.284) in G.D. Mahan Many-Particle Physics 3rd edition. 748 ! TODO Does not work with spinor 749 factor=spin_fact*qp_occ(band2,ik_ibz,spin) 750 call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_b1kmq_b2k,& 751 & Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule) 752 753 ! Include also the completeness correction in the sum rule 754 if (Ep%gwcomp==1) then 755 factor=-spin_fact*qp_occ(band2,ik_ibz,spin) 756 call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_enhigh_b2k,& 757 & Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule) 758 if (band1==Ep%nbnds) then 759 chi0_sumrule(:)=chi0_sumrule(:) + wtk_ltg(ik_bz)*spin_fact*qp_occ(band2,ik_ibz,spin)*deltaeGW_enhigh_b2k 760 end if 761 end if 762 763 end do !band2 764 end do !band1 765 766 ABI_FREE(gw_gbound) 767 if (Dtset%pawcross==1) then 768 ABI_FREE(gboundf) 769 end if 770 end do !ik_bz 771 772 ! Deallocation of arrays private to the spin loop. 773 ABI_FREE(igfftepsG0) 774 ABI_FREE(ur1_kmq_ibz) 775 ABI_FREE(ur2_k_ibz) 776 ABI_FREE(usr1_kmq) 777 ABI_FREE(ur2_k) 778 ABI_FREE(rhotwg) 779 ABI_FREE(tabr_k) 780 ABI_FREE(tabr_kmq) 781 if (allocated(green_w)) then 782 ABI_FREE(green_w) 783 end if 784 if (allocated(wfwfg)) then 785 ABI_FREE(wfwfg) 786 end if 787 if (allocated(green_enhigh_w)) then 788 ABI_FREE(green_enhigh_w) 789 end if 790 if (Psps%usepaw==1) then 791 call pawcprj_free(Cprj2_k) 792 ABI_DT_FREE(Cprj2_k) 793 call pawcprj_free(Cprj1_kmq) 794 ABI_DT_FREE(Cprj1_kmq) 795 if (Dtset%pawcross==1) then 796 ABI_FREE(ur_ae1) 797 ABI_FREE(ur_ae_onsite1) 798 ABI_FREE(ur_ps_onsite1) 799 ABI_FREE(ur_ae2) 800 ABI_FREE(ur_ae_onsite2) 801 ABI_FREE(ur_ps_onsite2) 802 ABI_FREE(tabrf_k) 803 ABI_FREE(tabrf_kmq) 804 ABI_FREE(gboundf) 805 ABI_FREE(igfftepsG0f) 806 end if 807 end if 808 end do !spin 809 810 ! After big loop over transitions, now MPI 811 ! Master took care of the contribution in case of metallic|spin polarized systems. 812 SELECT CASE (Ep%spmeth) 813 CASE (0) 814 ! Adler-Wiser 815 ! Collective sum of the contributions of each node. 816 ! Looping on frequencies to avoid problems with the size of the MPI packet 817 do io=1,Ep%nomega 818 call xmpi_sum(chi0(:,:,io),comm,ierr) 819 end do 820 821 CASE (1, 2) 822 ! Spectral method. 823 call hilbert_transform(Ep%npwe,Ep%nomega,Ep%nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,Ep%spmeth) 824 825 ! Deallocate here before xmpi_sum 826 if (allocated(sf_chi0)) then 827 ABI_FREE(sf_chi0) 828 end if 829 830 ! Collective sum of the contributions. 831 ! Looping over frequencies to avoid problems with the size of the MPI packet 832 do io=1,Ep%nomega 833 call xmpi_sum(chi0(:,:,io),comm,ierr) 834 end do 835 836 CASE DEFAULT 837 MSG_BUG("Wrong spmeth") 838 END SELECT 839 840 ! Divide by the volume 841 !$OMP PARALLEL WORKSHARE 842 chi0=chi0*weight/Cryst%ucvol 843 !$OMP END PARALLEL WORKSHARE 844 845 ! === Collect the sum rule === 846 ! * The pi factor comes from Im[1/(x-ieta)] = pi delta(x) 847 call xmpi_sum(chi0_sumrule,comm,ierr) 848 chi0_sumrule=chi0_sumrule*pi*weight/Cryst%ucvol 849 ! 850 ! ************************************************* 851 ! **** Now each node has chi0(q,G,Gp,Ep%omega) **** 852 ! ************************************************* 853 854 ! Impose Hermiticity (valid only for zero or purely imaginary frequencies) 855 ! MG what about metals, where we have poles around zero? 856 do io=1,Ep%nomega 857 if (ABS(REAL(Ep%omega(io))) <0.00001) then 858 do ig2=1,Ep%npwe 859 do ig1=1,ig2-1 860 chi0(ig2,ig1,io) = GWPC_CONJG(chi0(ig1,ig2,io)) 861 end do 862 end do 863 end if 864 end do 865 866 ! === Symmetrize chi0 in case of AFM system === 867 ! Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$. 868 ! Works only in case of magnetic group Shubnikov type IV. 869 if (Cryst%use_antiferro) call symmetrize_afm_chi0(Cryst,Gsph_epsG0,Ltg_q,Ep%npwe,Ep%nomega,chi0) 870 871 ! ===================== 872 ! ==== Free memory ==== 873 ! ===================== 874 ABI_FREE(bbp_ks_distrb) 875 876 if (allocated(gw_gfft)) then 877 ABI_FREE(gw_gfft) 878 end if 879 if (allocated(kkweight)) then 880 ABI_FREE(kkweight) 881 end if 882 if (allocated(omegasf)) then 883 ABI_FREE(omegasf) 884 end if 885 if (allocated(gspfft_igfft)) then 886 ABI_FREE(gspfft_igfft) 887 end if 888 889 call gsph_free(Gsph_FFT) 890 891 ! deallocation for PAW. 892 if (Psps%usepaw==1) then 893 call pawpwij_free(Pwij) 894 ABI_DT_FREE(Pwij) 895 if (allocated(Pwij_fft)) then 896 call pawpwij_free(Pwij_fft) 897 ABI_DT_FREE(Pwij_fft) 898 end if 899 end if 900 901 if(dtset%ucrpa>=1) then 902 ABI_DEALLOCATE(coeffW_BZ) 903 endif 904 905 call timab(331,2,tsec) 906 call cwtime(cpu_time,wall_time,gflops,"stop") 907 write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time 908 909 DBG_EXIT("COLL") 910 911 end subroutine cchi0