TABLE OF CONTENTS
ABINIT/cohsex_me [ Functions ]
NAME
cohsex_me
FUNCTION
Calculate diagonal and off-diagonal matrix elements of the SEX or COHSEX self-energy operator.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (FB, 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
sigmak_ibz=Index of the k-point in the IBZ. minbnd, maxbnd= min and Max band index for GW correction (for this k-point) iomode=Option defining the file format of the SCR file (Fortran, NETCDF) Er <Epsilonm1_results> (see the definition of this structured datatype) %mqmem=if 0 use out-of-core method in which a single q-slice of espilon is read inside the loop over k %nomega_i=Number of imaginary frequencies. %nomega_r=Number of real frequencies. %nomega=Total number of frequencies. Gsph_c<gsphere_t>= info on the G-sphere for Sigma_x %nsym=number of symmetry operations %rottb(ng,timrev,nsym)=index of (IS) G where I is the identity or the inversion operation and G is one of the ng vectors in reciprocal space %timrev=2 if time-reversal symmetry is used, 1 otherwise %gvec(3,Sigp%npwc)=integer coordinates of each plane wave in reciprocal space ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated Ltg_k datatype containing information on the little group Kmesh <kmesh_t> %nbz=Number of points in the BZ %nibz=Number of points in IBZ %kibz(3,nibz)=k-point coordinates, irreducible Brillouin zone %kbz(3,nbz)=k-point coordinates, full Brillouin zone %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered %ktabp(nbz)= phase factor associated to tnons gwc_ngfft(18)=Information about 3D FFT for the oscillator strengths used for the correlation part, Vcp <vcoul_t datatype> containing information on the cutoff technique %vc_sqrt(npwc,nqibz)= square-root of the coulombian potential for q-points in the IBZ Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data Pawang <type(pawang_type)>=paw angular mesh and related data Psps <type(pseudopotential_type)>=variables related to pseudopotentials %usepaw=1 for PAW, 0 for NC pseudopotentials. Qmesh <kmesh_t> : datatype gathering information of the q-mesh used %ibz=q points where $\tilde\epsilon^{-1}$ has been computed %bz(3,nqbz)=coordinates of all q-points in BZ Sigp <sigparams_t> (see the definition of this structured datatype) Cryst<crystal_t>=Info on unit cell and symmetries %natom=number of atoms in unit cell %ucvol=unit cell volume %nsym=number of symmetry operations %typat(natom)=type of each atom much slower but it requires less memory QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,Kmesh%nibz,%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,Kmesh%nibz,%nsppol)=occupation numbers, for each k point in IBZ, each band and spin Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave. allQP_sym(%nkibz,%nsppol)<esymm_t>=Datatype collecting data on the irreducible representaions of the little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states. Sr=sigma_t (see the definition of this structured datatype)
OUTPUT
NOTES
1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) is included. 2) The calculation of energy derivative is based on finite elements. 3) On the symmetrization of Sigma matrix elements ***/ If Sk = k+G0 then M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G} (k,q) If -Sk = k+G0 then M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q) Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ. 4) The unitary transformation relating wavefunctions at symmetric k-points should be taken into account during the symmetrization of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however, it is possible to make an invariant by just including all the degenerate states and averaging the final results over the degenerate subset.
PARENTS
sigma
CHILDREN
calc_coh,calc_wfwfg,epsm1_symmetrizer,esymm_symmetrize_mels,findqg0 get_bz_item,get_epsm1,get_gftt,gsph_fft_tabs,hermitianize littlegroup_print,paw_rho_tw_g,paw_symcprj,pawcprj_alloc,pawcprj_copy pawcprj_free,pawpwij_free,pawpwij_init,rho_tw_g,rotate_fft_mesh sigma_distribute_bks,timab,wfd_change_ngfft,wfd_get_cprj wfd_get_many_ur,wfd_get_ur,wrtout,xgemv,xmpi_sum
SOURCE
99 #if defined HAVE_CONFIG_H 100 #include "config.h" 101 #endif 102 103 #include "abi_common.h" 104 105 subroutine cohsex_me(sigmak_ibz,ikcalc,nomega_sigc,minbnd,maxbnd,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_c,Vcp,& 106 & Kmesh,Qmesh,Ltg_k,Pawtab,Pawang,Paw_pwff,Psps,Wfd,allQP_sym,gwc_ngfft,iomode,prtvol,sigcme_tmp) 107 108 use defs_basis 109 use defs_datatypes 110 use m_defs_ptgroups 111 use m_gwdefs !, only : czero_gw, cone_gw, j_gw, sigparams_t, sigma_type_from_key, sigma_is_herm 112 use m_xmpi 113 use m_errors 114 use m_profiling_abi 115 116 use m_blas, only : xdotc, xgemv 117 use m_numeric_tools, only : hermitianize, imin_loc 118 use m_geometry, only : normv 119 use m_crystal, only : crystal_t 120 use m_bz_mesh, only : kmesh_t, get_BZ_item, findqg0, littlegroup_t, littlegroup_print 121 use m_gsphere, only : gsphere_t, gsph_fft_tabs 122 use m_fft_mesh, only : get_gftt, rotate_fft_mesh, cigfft 123 use m_vcoul, only : vcoul_t 124 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g 125 use m_wfd, only : wfd_get_ur, wfd_t, wfd_get_cprj, wfd_change_ngfft, wfd_get_many_ur, wfd_sym_ur 126 use m_oscillators, only : rho_tw_g, calc_wfwfg 127 use m_screening, only : epsm1_symmetrizer, get_epsm1, epsilonm1_results 128 use m_esymm, only : esymm_t, esymm_symmetrize_mels, esymm_failed 129 use m_sigma, only : sigma_t, sigma_distribute_bks 130 use m_pawang, only : pawang_type 131 use m_pawtab, only : pawtab_type 132 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap 133 134 !This section has been created automatically by the script Abilint (TD). 135 !Do not modify the following lines by hand. 136 #undef ABI_FUNC 137 #define ABI_FUNC 'cohsex_me' 138 use interfaces_14_hidewrite 139 use interfaces_18_timing 140 use interfaces_65_paw 141 use interfaces_70_gw, except_this_one => cohsex_me 142 !End of the abilint section 143 144 implicit none 145 146 !Arguments ------------------------------------ 147 !scalars 148 integer,intent(in) :: sigmak_ibz,ikcalc,prtvol,iomode,nomega_sigc,minbnd,maxbnd 149 type(crystal_t),intent(in) :: Cryst 150 type(ebands_t),target,intent(in) :: QP_BSt 151 type(kmesh_t),intent(in) :: Kmesh,Qmesh 152 type(vcoul_t),intent(in) :: Vcp 153 type(Epsilonm1_results),intent(inout) :: Er 154 type(gsphere_t),intent(in) :: Gsph_c 155 type(littlegroup_t),intent(in) :: Ltg_k 156 type(Pseudopotential_type),intent(in) :: Psps 157 type(pawang_type),intent(in) :: pawang 158 type(sigparams_t),target,intent(in) :: Sigp 159 type(sigma_t),intent(in) :: Sr 160 type(wfd_t),target,intent(inout) :: Wfd 161 !arrays 162 integer,intent(in) :: gwc_ngfft(18) 163 complex(dpc),intent(out) :: sigcme_tmp(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Wfd%nsppol*Sigp%nsig_ab) 164 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat) 165 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 166 type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol) 167 168 !Local variables ------------------------------ 169 !scalars 170 integer,parameter :: tim_fourdp=2,ndat1=1 171 integer :: iab,ib,ib1,ib2,ierr,ig,ii,iik,itim_q,i1,i2,npwc 172 integer :: ik_bz,ik_ibz,io,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx 173 integer :: band,band1,band2,idle,rank 174 integer :: jik,jk_bz,jk_ibz,kb,nspinor,nsppol 175 integer :: nomega_tot,nq_summed,ispinor,ibsp,dimcprj_gw 176 integer :: spad,spadc,spadc1,spadc2,irow,my_nbks 177 integer :: ndegs,wtqm,wtqp,mod10 178 integer :: isym_kgw,isym_ki,gwc_mgfft,use_padfft,gwc_fftalga,gwc_nfftot,ifft,npw_k 179 real(dp) :: fact_sp,theta_mu_minus_e0i,tol_empty,gw_gsq 180 complex(dpc) :: ctmp,ph_mkgwt,ph_mkt 181 logical :: iscompatibleFFT,q_is_gamma 182 character(len=500) :: msg 183 !arrays 184 integer :: g0(3),spinor_padc(2,4),nbv_ks(Kmesh%nibz,Wfd%nsppol) 185 integer,allocatable :: proc_distrb(:,:,:),coh_distrb(:,:,:,:),degtab(:,:,:) 186 integer,allocatable :: igfftcg0(:),gw_gfft(:,:),gw_gbound(:,:),irottb(:,:),ktabr(:,:) 187 integer :: got(Wfd%nproc) 188 real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),q0(3),tsec(2),qbz(3),spinrot_kbz(4),spinrot_kgw(4) 189 real(dp),ABI_CONTIGUOUS pointer :: qp_ene(:,:,:),qp_occ(:,:,:) 190 complex(gwpc) :: sigcohme(Sigp%nsig_ab) 191 complex(gwpc),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:),sigsex(:) 192 complex(gwpc),allocatable :: epsm1_qbz(:,:,:) 193 complex(gwpc),allocatable :: sigc_ket(:,:) 194 complex(gwpc),allocatable :: rhotwg_ki(:,:) 195 complex(gwpc),allocatable :: sigctmp(:,:) 196 complex(gwpc),allocatable :: wfr_bdgw(:,:),ur_sum(:),wf1swf2_g(:) 197 complex(gwpc),ABI_CONTIGUOUS pointer :: cg_jb(:),cg_sum(:) 198 complex(dpc),allocatable :: sym_cme(:,:,:,:),sigc(:,:,:,:,:) 199 logical :: rank_mask(Wfd%nproc),can_symmetrize(Wfd%nsppol) 200 logical,allocatable :: bks_mask(:,:,:) 201 type(sigijtab_t),pointer :: Sigcij_tab(:) 202 type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:) 203 type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:) 204 type(esymm_t),pointer :: QP_sym(:) 205 206 !************************************************************************ 207 208 DBG_ENTER("COLL") 209 210 call timab(423,1,tsec) ! cohsex_me 211 212 ! Initial check 213 ABI_CHECK(Sr%nomega_r == Sigp%nomegasr,"") 214 ABI_CHECK(Sr%nomega4sd == Sigp%nomegasrd,"") 215 !ABI_CHECK(Sigp%npwc==Gsph_c%ng,"") 216 217 ! Initialize some values 218 nspinor = Wfd%nspinor; nsppol = Wfd%nsppol 219 npwc = sigp%npwc 220 spinor_padc = RESHAPE([0, 0, npwc, npwc, 0, npwc, npwc,0], [2, 4]) 221 222 qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ 223 224 ! Extract the symmetries of the bands for this k-point 225 QP_sym => allQP_sym(sigmak_ibz,1:nsppol) 226 227 ! Index of the GW point in the BZ array, its image in IBZ and time-reversal === 228 jk_bz=Sigp%kptgw2bz(ikcalc) 229 call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt) 230 !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk) 231 spinrot_kgw=Cryst%spinrot(:,isym_kgw) 232 ib1 = minbnd; ib2 = maxbnd 233 234 write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,& 235 & ' Calculating <nk|Sigma_c(omega)|nk> at k = ',kgw(:),ch10,& 236 & ' bands n = from ',ib1,' to ',ib2,ch10 237 call wrtout(std_out,msg,'COLL') 238 239 if (ANY(gwc_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,gwc_ngfft) 240 gwc_mgfft = MAXVAL(gwc_ngfft(1:3)) 241 gwc_fftalga = gwc_ngfft(7)/100 !; gwc_fftalgc=MOD(gwc_ngfft(7),10) 242 243 can_symmetrize = .FALSE. 244 if (Sigp%symsigma>0) then 245 can_symmetrize = .TRUE. 246 if (Sigp%gwcalctyp >= 20) then 247 do spin=1,Wfd%nsppol 248 can_symmetrize(spin) = .not. esymm_failed(QP_sym(spin)) 249 if (.not.can_symmetrize(spin)) then 250 write(msg,'(a,i0,4a)')" Symmetrization cannot be performed for spin: ",spin,ch10,& 251 & " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg) 252 MSG_WARNING(msg) 253 end if 254 end do 255 end if 256 if (nspinor == 2) MSG_WARNING('Symmetrization with nspinor=2 not implemented') 257 end if 258 259 mod10=MOD(Sigp%gwcalctyp, 10) 260 261 call timab(491,1,tsec) ! csigme(tot) Overall clock. TODO check this 262 call timab(495,1,tsec) ! csigme (SigC) 263 264 ! Normalization of theta_mu_minus_e0i 265 ! If nsppol==2, qp_occ $\in [0,1]$ 266 SELECT CASE (nsppol) 267 CASE (1) 268 fact_sp=half; tol_empty=0.01 ! below this value the state is assumed empty 269 if (nspinor==2) then 270 fact_sp=one; tol_empty=0.005 ! below this value the state is assumed empty 271 end if 272 CASE (2) 273 fact_sp=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic 274 CASE DEFAULT ! spin unpolarized system is treated using nsppol==2 275 MSG_BUG('Wrong nsppol') 276 END SELECT 277 278 call timab(442,1,tsec) ! csigme(init0) 279 280 ! Precalculate the FFT index of $(R^{-1}(r-\tau))$ === 281 ! S=\transpose R^{-1} and k_BZ = S k_IBZ 282 ! irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk. 283 gwc_nfftot = PRODUCT(gwc_ngfft(1:3)) 284 ABI_MALLOC(irottb,(gwc_nfftot,Cryst%nsym)) 285 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwc_ngfft,irottb,iscompatibleFFT) 286 if (.not.iscompatibleFFT) then 287 MSG_WARNING("FFT mesh is not compatible with symmetries. Results might be affected by large errors!") 288 end if 289 290 ABI_MALLOC(ktabr,(gwc_nfftot, Kmesh%nbz)) 291 do ik_bz=1,Kmesh%nbz 292 isym=Kmesh%tabo(ik_bz) 293 do ifft=1,gwc_nfftot 294 ktabr(ifft,ik_bz)=irottb(ifft,isym) 295 end do 296 end do 297 ABI_FREE(irottb) 298 299 ! The number of occupied states for each point in the IBZ and spin. 300 ! nbv_ks(:,:) = COUNT(qp_occ>=tol_empty,DIM=1) MG: g95 returns random numbers, likely a bug in the compiler 301 do spin=1,nsppol 302 do ik_ibz=1,Kmesh%nibz 303 nbv_ks(ik_ibz,spin) = COUNT(qp_occ(:,ik_ibz,spin)>=tol_empty) 304 end do 305 end do 306 307 ! (b,k,s) mask for MPI distribution of the sum over occupied states in the BZ. 308 ABI_MALLOC(bks_mask,(Wfd%mband,Kmesh%nbz,nsppol)) 309 bks_mask=.FALSE. 310 do spin=1,nsppol 311 do ik_bz=1,Kmesh%nbz 312 ik_ibz = Kmesh%tab(ik_bz) 313 bks_mask(1:nbv_ks(ik_ibz,spin),ik_bz,spin) = .TRUE. 314 end do 315 end do 316 317 ! Distribute the individual terms of the sum over the BZ taking into account symmetries and MPI memory distribution. 318 ! got is used to optimize the distribution if more than one node can calculate the same (b,k,s) element. 319 got=0 320 ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,nsppol)) 321 call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,& 322 & proc_distrb,got,bks_mask,global=.TRUE.) 323 324 ABI_FREE(bks_mask) 325 326 write(msg,'(a,i0,a)')" Will sum ",my_nbks," (b,k,s) occupied states in (COHSEX|SEX)." 327 call wrtout(std_out,msg,'PERS') 328 329 Sigcij_tab => Sigp%Sigcij_tab(ikcalc,1:nsppol) 330 331 if (mod10==SIG_COHSEX) then 332 ! Distribute the COHSEX terms, taking into account the symmetries of the Sigma_ij matrix. 333 ABI_MALLOC(coh_distrb,(ib1:ib2,ib1:ib2,Kmesh%nbz,nsppol)) 334 335 coh_distrb = xmpi_undefined_rank 336 do spin=1,nsppol 337 do ik_bz=1,Kmesh%nbz 338 if (ANY(proc_distrb(:,ik_bz,spin) /= xmpi_undefined_rank) ) then ! This BZ point will be calculated. 339 rank_mask = .FALSE. ! To select only those nodes that will treat (k,s). 340 do band=1,Wfd%mband 341 rank = proc_distrb(band,ik_bz,spin) 342 if (rank /= xmpi_undefined_rank) rank_mask(rank+1)=.TRUE. 343 end do 344 do band2=ib1,ib2 345 do irow=1,Sigcij_tab(spin)%col(band2)%size1 ! Looping over the upper triangle of sigma_ij with non-zero elements. 346 band1 = Sigcij_tab(spin)%col(band2)%bidx(irow) 347 idle = imin_loc(got,mask=rank_mask) 348 got(idle) = got(idle)+1 349 coh_distrb(band1,band2,ik_bz,spin) = idle-1 350 end do 351 end do 352 end if 353 end do 354 end do 355 356 write(msg,'(a,i0,a)')" will treat ",COUNT(coh_distrb==Wfd%my_rank)," COH terms." 357 call wrtout(std_out,msg,'PERS') 358 end if 359 360 ABI_MALLOC(rhotwg_ki, (npwc * nspinor, minbnd:maxbnd)) 361 rhotwg_ki=czero_gw 362 ABI_MALLOC(rhotwg, (npwc * nspinor)) 363 ABI_MALLOC(rhotwgp ,(npwc * nspinor)) 364 ABI_MALLOC(vc_sqrt_qbz, (npwc)) 365 366 ! Additional allocations for PAW 367 if (Psps%usepaw==1) then 368 ABI_DT_MALLOC(Cprj_ksum,(Cryst%natom,nspinor)) 369 call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm) 370 371 ! For COHSEX we need the onsite terms of the PW on the FFT mesh. 372 ! gw_gfft is the set of plane waves in the FFT Box for the oscillators. 373 if (mod10==SIG_COHSEX) then 374 ABI_MALLOC(gw_gfft,(3,gwc_nfftot)) 375 q0=zero 376 call get_gftt(gwc_ngfft,q0,Cryst%gmet,gw_gsq,gw_gfft) 377 ABI_DT_MALLOC(Pwij_fft,(Psps%ntypat)) 378 call pawpwij_init(Pwij_fft,gwc_nfftot,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 379 end if 380 end if ! usepaw==1 381 382 ! === Calculate total number of frequencies and allocate related arrays === 383 ! sigcme2 is used to accumulate the diagonal matrix elements over k-points and 384 ! GW bands, used only in case of ppmodel 3 and 4 (TODO save memory) 385 nomega_tot=Sr%nomega_r+Sr%nomega4sd 386 387 ABI_MALLOC(sigctmp, (nomega_sigc,Sigp%nsig_ab)) 388 sigctmp = czero_gw 389 ABI_MALLOC(sigc_ket, (npwc*nspinor, nomega_sigc)) 390 391 if (mod10==SIG_COHSEX) then 392 ABI_MALLOC(wf1swf2_g,(gwc_nfftot*nspinor)) 393 end if 394 395 ! Arrays storing the contribution given by the Hermitian/anti-Hermitian part of \Sigma_c 396 !allocate(aherm_sigc_ket(npwc*nspinor,nomega_sigc)) 397 !allocate( herm_sigc_ket(npwc*nspinor,nomega_sigc)) 398 ABI_MALLOC(sigsex,(npwc)) 399 sigcme_tmp=czero 400 401 ABI_MALLOC(sigc,(2,nomega_sigc,ib1:ib2,ib1:ib2,nsppol*Sigp%nsig_ab)) 402 sigc=czero 403 404 ! Here we divide the states where the QP energies are required into complexes. Note however that this approach is not 405 ! based on group theory, and it might lead to spurious results in case of accidental degeneracies. 406 nq_summed=Kmesh%nbz 407 if (Sigp%symsigma>0) then 408 call littlegroup_print(Ltg_k,std_out,prtvol,'COLL') 409 nq_summed=SUM(Ltg_k%ibzq(:)) 410 ! 411 ! Find number of degenerate states and number of bands in each subspace 412 ! The tolerance is a little bit arbitrary (0.001 eV) 413 ! It could be reduced, in particular in case of nearly accidental degeneracies 414 ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,nsppol)) 415 degtab=0 416 do spin=1,nsppol 417 do ib=ib1,ib2 418 do jb=ib1,ib2 419 if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin))<0.001/Ha_ev) then 420 degtab(ib,jb,spin)=1 421 end if 422 end do 423 end do 424 end do 425 end if !symsigma 426 427 write(msg,'(2a,i6,a)')ch10,' calculation status ( ',nq_summed,' to be completed):' 428 call wrtout(std_out,msg,'COLL') 429 430 ! TODO if single q (ex molecule) dont allocate epsm1q, avoid waste of memory 431 ABI_STAT_MALLOC(epsm1_qbz, (npwc, npwc, 1), ierr) 432 ABI_CHECK(ierr==0, "out-of-memory in epsm1_qbz") 433 ABI_MALLOC(igfftcg0,(Gsph_c%ng)) 434 435 ! Out-of-core solution for epsilon. 436 if (Er%mqmem==0) then 437 MSG_COMMENT('Reading q-slices from file. Slower but less memory.') 438 end if 439 440 call timab(442,2,tsec) 441 442 ! ========================================== 443 ! ==== Fat loop over k_i in the full BZ ==== 444 ! ========================================== 445 ABI_MALLOC(ur_sum,(gwc_nfftot*nspinor)) 446 447 do spin=1,nsppol 448 if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE 449 450 ABI_MALLOC(wfr_bdgw,(gwc_nfftot*nspinor,ib1:ib2)) 451 call wfd_get_many_ur(Wfd, [(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw) 452 453 if (Wfd%usepaw==1) then 454 ! Load cprj for GW states, note the indexing. 455 dimcprj_gw=nspinor*(ib2-ib1+1) 456 ABI_DT_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1)) 457 call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm) 458 ibsp=ib1 459 do jb=ib1,ib2 460 call wfd_get_cprj(Wfd,jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 461 call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 462 call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1))) 463 ibsp=ibsp+nspinor 464 end do 465 end if 466 467 do ik_bz=1,Kmesh%nbz 468 ! Parallelization over k-points and spin 469 ! For the spin there is another check in the inner loop. 470 if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE 471 472 call timab(443,1,tsec) ! csigme (initq) 473 474 ! Find the corresponding irreducible k-point 475 call get_BZ_item(Kmesh,ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt) 476 spinrot_kbz(:)=Cryst%spinrot(:,isym_ki) 477 478 ! Identify q and G0 where q+G0=k_GW-k_i 479 kgw_m_ksum=kgw-ksum 480 call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0) 481 482 ! Symmetrize the matrix elements. 483 ! Sum only q"s in IBZ_k. In this case elements are weighted 484 ! according to wtqp and wtqm. wtqm is for time-reversal. 485 wtqp=1; wtqm=0 486 if (can_symmetrize(spin)) then 487 if (Ltg_k%ibzq(iq_bz)/=1) CYCLE 488 wtqp=0; wtqm=0 489 do isym=1,Ltg_k%nsym_sg 490 wtqp=wtqp+Ltg_k%wtksym(1,isym,iq_bz) 491 wtqm=wtqm+Ltg_k%wtksym(2,isym,iq_bz) 492 end do 493 end if 494 495 !%write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank 496 !%call wrtout(std_out,msg,'PERS') 497 498 ! Find the corresponding irred q-point. 499 call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q) 500 q_is_gamma = (normv(qbz,Cryst%gmet,"G") < GW_TOL_W0) 501 502 ! Tables for the FFT of the oscillators. 503 ! a) FFT index of the G-G0. 504 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 505 ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2)) 506 call gsph_fft_tabs(Gsph_c,g0,gwc_mgfft,gwc_ngfft,use_padfft,gw_gbound,igfftcg0) 507 if ( ANY(gwc_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 508 if (use_padfft==0) then 509 ABI_FREE(gw_gbound) 510 ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2*use_padfft)) 511 end if 512 513 if (Psps%usepaw==1) then 514 ! Get PAW oscillator matrix elements $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form. 515 ABI_DT_MALLOC(Pwij_qg,(Psps%ntypat)) 516 q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT 517 call pawpwij_init(Pwij_qg, npwc, q0, Gsph_c%gvec, Cryst%rprimd, Psps, Pawtab, Paw_pwff) 518 end if 519 520 if (Er%mqmem==0) then 521 ! Read q-slice of epsilon^{-1}|chi0 in Er%epsm1(:,:,:,1) (much slower but less memory). 522 call get_epsm1(Er,Vcp,0,0,iomode,xmpi_comm_self,iqibzA=iq_ibz) 523 end if 524 525 ! Only omega==0 for SEX or COHSEX 526 call Epsm1_symmetrizer(iq_bz, 1, npwc, Er, Gsph_c, Qmesh, .True., epsm1_qbz) 527 528 ! Get Fourier components of the Coulomb interaction in the BZ. 529 ! In 3D systems, neglecting umklapp, vc(Sq,sG)=vc(q,G)=4pi/|q+G| 530 ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S. 531 do ig=1,npwc 532 vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q)) = Vcp%vc_sqrt(ig,iq_ibz) 533 end do 534 535 call timab(443,2,tsec) ! csigme (initq) 536 537 ! Sum over bands. 538 do ib=1,Sigp%nbnds 539 ! Parallelism over spin 540 ! This processor has this k-point but what about spin? 541 if (proc_distrb(ib,ik_bz,spin)/=Wfd%my_rank) CYCLE 542 543 ! Skip empty state ib for HF, SEX, and COHSEX. 544 if (qp_occ(ib,ik_ibz,spin)<tol_empty) CYCLE 545 546 theta_mu_minus_e0i=fact_sp*qp_occ(ib,ik_ibz,spin) 547 548 call wfd_get_ur(Wfd,ib,ik_ibz,spin,ur_sum) 549 550 if (Psps%usepaw==1) then 551 ! Load cprj for point ksum, this spin or spinor and *THIS* band. 552 ! TODO MG I could avoid doing this but I have to exchange spin and bands ??? 553 ! For sure there is a better way to do this! 554 call wfd_get_cprj(Wfd,ib,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 555 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 556 end if 557 558 do jb=ib1,ib2 559 ! Get all <k-q,ib,s|e^{-i(q+G).r}|s,jb,k>, at once. 560 call rho_tw_g(nspinor,npwc,gwc_nfftot,ndat1,gwc_ngfft,1,use_padfft,igfftcg0,gw_gbound,& 561 & ur_sum ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, & 562 & wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,& 563 & nspinor,rhotwg_ki(:,jb)) 564 565 if (Psps%usepaw==1) then 566 ! Add on-site contribution, projectors are already in BZ !TODO Recheck this! 567 i2=jb; if (nspinor==2) i2=(2*jb-1) 568 spad=(nspinor-1) 569 call paw_rho_tw_g(npwc,nspinor,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_c%gvec,& 570 & Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),Pwij_qg,rhotwg_ki(:,jb)) 571 end if 572 573 ! Multiply by the square root of the Coulomb term. 574 ! In 3-D systems, the factor sqrt(4pi) is included) 575 do ii=1,nspinor 576 spad = (ii-1) * npwc 577 rhotwg_ki(spad+1:spad+npwc,jb) = rhotwg_ki(spad+1:spad+npwc,jb)*vc_sqrt_qbz(1:npwc) 578 end do 579 580 ! === Treat analytically the case q --> 0 === 581 ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma 582 ! while the Colulomb term is integrated out 583 ! * In the scalar case we have nonzero contribution only if ib==jb 584 ! * For nspinor==2 evalute <ib,up|jb,up> and <ib,dwn|jb,dwn>, 585 ! impose orthonormalization since npwwfn might be < npwvec. 586 if (ik_bz == jk_bz) then 587 if (nspinor == 1) then 588 rhotwg_ki(1, jb) = czero_gw 589 if (ib==jb) rhotwg_ki(1, jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) 590 else 591 npw_k = Wfd%npwarr(ik_ibz) 592 rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero 593 if (ib == jb) then 594 cg_sum => Wfd%Wave(ib, ik_ibz, spin)%ug 595 cg_jb => Wfd%Wave(jb, jk_ibz, spin)%ug 596 ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1) 597 rhotwg_ki(1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp) 598 ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1) 599 rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp) 600 ! PAW is missing 601 602 !rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half) 603 !rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half) 604 end if 605 end if 606 end if 607 end do !jb Got all matrix elements from minbnd up to maxbnd. 608 609 do kb=ib1,ib2 610 ! Get the ket \Sigma|\phi_{k,kb}> according to the method. 611 rhotwgp(:) = rhotwg_ki(:,kb) 612 sigc_ket = czero_gw 613 614 ! SEX part. TODO add check on theta_mu_minus_e0i 615 do ispinor=1,nspinor 616 spadc = (ispinor-1) * npwc 617 call XGEMV('N',npwc,npwc,cone_gw,epsm1_qbz(:,:,1),npwc,rhotwgp(1+spadc:),1,czero_gw,sigsex,1) 618 619 sigsex(:)= -theta_mu_minus_e0i*sigsex(:) 620 621 do io=1,nomega_tot ! nomega==1 as SEX is energy independent. 622 sigc_ket(spadc+1:spadc+npwc,io) = sigsex(:) 623 end do 624 end do 625 626 ! Loop over the non-zero row elements of this column. 627 ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS. 628 ! 2) If gwcalctyp>=20: 629 ! * Only off-diagonal elements connecting states with same character. 630 ! * Only the upper triangle if HF, SEX, or COHSEX. 631 do irow=1,Sigcij_tab(spin)%col(kb)%size1 632 jb = Sigcij_tab(spin)%col(kb)%bidx(irow) 633 rhotwg = rhotwg_ki(:,jb) 634 635 ! Calculate <\phi_j|\Sigma_c|\phi_k> 636 ! Different freqs according to method (AC or Perturbative), see nomega_sigc. 637 do iab=1,Sigp%nsig_ab 638 spadc1=spinor_padc(1,iab); spadc2=spinor_padc(2,iab) 639 do io=1,nomega_sigc 640 sigctmp(io,iab) = XDOTC(npwc,rhotwg(spadc1+1:),1,sigc_ket(spadc2+1:,io),1) 641 end do 642 end do 643 644 ! TODO: save wf1swf2_g to avoid having to recalculate it at each q-point. 645 if (mod10==SIG_COHSEX) then 646 ! Evaluate Static COH. TODO add spinor. 647 if (coh_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank) then 648 ! COH term is done only once for each k-point. 649 ! It does not depend on the index ib summed over. 650 coh_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank 651 652 #if 1 653 call calc_wfwfg(ktabr(:,jk_ibz),jik, spinrot_kgw, & ! TODO why jk_ibz? 654 & gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g) 655 #else 656 ABI_CHECK(jik==1,"jik") 657 call calc_wfwfg(ktabr(:,jk_bz),jik, spinrot_kgw, & 658 gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g) 659 #endif 660 661 if (Psps%usepaw==1) then 662 i1=jb; i2=kb 663 if (nspinor==2) then 664 i1=(2*jb-1); i2=(2*kb-1) 665 end if 666 spad=(nspinor-1) 667 call paw_rho_tw_g(gwc_nfftot,Sigp%nsig_ab,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,& 668 & gw_gfft,Cprj_kgw(:,i1:i1+spad),Cprj_kgw(:,i2:i2+spad),Pwij_fft,wf1swf2_g) 669 end if 670 671 call calc_coh(nspinor,Sigp%nsig_ab,gwc_nfftot,gwc_ngfft,npwc,Gsph_c%gvec,wf1swf2_g,epsm1_qbz(:,:,1),& 672 & vc_sqrt_qbz,Vcp%i_sz,iq_ibz,(jb==kb),sigcohme) 673 674 do io=1,nomega_sigc ! Should be 1 675 sigctmp(io,:) = sigctmp(io,:)+sigcohme(:) 676 end do 677 678 end if 679 end if ! COHSEX 680 681 ! Accumulate and, in case, symmetrize matrix elements of Sigma_c. 682 do iab=1,Sigp%nsig_ab 683 is_idx = spin; if (nspinor==2) is_idx=iab 684 685 sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + & 686 & (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab)) 687 688 sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp* sigctmp(:,iab) 689 sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab)) 690 ! TODO this should be the contribution coming from the anti-hermitian part. 691 end do 692 end do !jb used to calculate matrix elements of $\Sigma$ 693 694 end do !kb to calculate matrix elements of $\Sigma$ 695 end do !ib 696 697 ! Deallocate k-dependent quantities. 698 ABI_FREE(gw_gbound) 699 if (Psps%usepaw==1) then 700 call pawpwij_free(Pwij_qg) 701 ABI_DT_FREE(Pwij_qg) 702 end if 703 end do !ik_bz 704 705 ABI_FREE(wfr_bdgw) 706 if (Wfd%usepaw==1) then 707 call pawcprj_free(Cprj_kgw) 708 ABI_DT_FREE(Cprj_kgw) 709 end if 710 end do !spin 711 712 ABI_FREE(igfftcg0) 713 714 ! Gather contributions from all the CPUs. 715 call xmpi_sum(sigcme_tmp, wfd%comm, ierr) 716 call xmpi_sum(sigc, wfd%comm, ierr) 717 718 ! Multiply by constants 719 ! For 3D systems sqrt(4pi) is included in vc_sqrt_qbz. 720 sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz) 721 sigc = sigc /(Cryst%ucvol*Kmesh%nbz) 722 723 ! If we have summed over the IBZ_q now we have to average over degenerate states. 724 ! Presently only diagonal terms are considered 725 ! TODO it does not work if nspinor==2. 726 do spin=1,nsppol 727 if (can_symmetrize(spin)) then 728 ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, sigp%nsig_ab)) 729 sym_cme=czero 730 731 ! Average over degenerate diagonal elements. 732 ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results. 733 ! another good reason to use a strict criterion for the tolerance on eigenvalues. 734 do ib=ib1,ib2 735 ndegs=0 736 do jb=ib1,ib2 737 if (degtab(ib,jb,spin)==1) then 738 if (nspinor == 1) then 739 sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:, :, jb, jb, spin), dim=1) 740 else 741 do ii=1,sigp%nsig_ab 742 sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:, :, jb, jb, ii), dim=1) 743 end do 744 end if 745 746 end if 747 ndegs = ndegs + degtab(ib,jb,spin) 748 end do 749 sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs 750 end do 751 752 if (Sigp%gwcalctyp >= 20) then 753 call esymm_symmetrize_mels(QP_sym(spin),ib1,ib2,sigc(:,1,:,:,spin),sym_cme(1,:,:,1)) 754 end if 755 756 ! Copy symmetrized values. 757 do ib=ib1,ib2 758 do jb=ib1,ib2 759 if (nspinor == 1) then 760 sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1) 761 else 762 sigcme_tmp(:,ib,jb,:) = sym_cme(:,ib,jb,:) 763 end if 764 end do 765 end do 766 ABI_FREE(sym_cme) 767 end if 768 end do 769 770 ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX) 771 if (Sigp%gwcalctyp>=20 .and. sigma_is_herm(Sigp) ) then 772 ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!") 773 do spin=1,nsppol 774 do io=1,nomega_sigc 775 call hermitianize(sigcme_tmp(io,:,:,spin),"Upper") 776 end do 777 end do 778 end if 779 780 ! =========================== 781 ! ==== Deallocate memory ==== 782 ! =========================== 783 if (Psps%usepaw==1) then 784 if (allocated(gw_gfft)) then 785 ABI_FREE(gw_gfft) 786 end if 787 call pawcprj_free(Cprj_ksum) 788 ABI_DT_FREE(Cprj_ksum) 789 if (allocated(Pwij_fft)) then 790 call pawpwij_free(Pwij_fft) 791 ABI_DT_FREE(Pwij_fft) 792 end if 793 end if 794 795 ABI_FREE(ktabr) 796 ABI_FREE(ur_sum) 797 ABI_FREE(rhotwg_ki) 798 ABI_FREE(rhotwg) 799 ABI_FREE(rhotwgp) 800 ABI_FREE(vc_sqrt_qbz) 801 ABI_FREE(sigc_ket) 802 ABI_FREE(epsm1_qbz) 803 ABI_FREE(sigctmp) 804 ABI_FREE(sigc) 805 ABI_FREE(sigsex) 806 ABI_FREE(proc_distrb) 807 if (mod10==SIG_COHSEX) then 808 ABI_FREE(wf1swf2_g) 809 ABI_FREE(coh_distrb) 810 end if 811 812 if (allocated(degtab)) then 813 ABI_FREE(degtab) 814 end if 815 816 call timab(495,2,tsec) ! csigme(SigC) 817 call timab(491,2,tsec) 818 call timab(423,2,tsec) ! cohsex_me 819 820 DBG_EXIT("COLL") 821 822 end subroutine cohsex_me