TABLE OF CONTENTS
ABINIT/calc_sigx_me [ Functions ]
NAME
calc_sigx_me
FUNCTION
Calculate diagonal and off-diagonal matrix elements of the exchange part of the 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) Gsph_x<gsphere_t>= Info on the G-sphere used 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,ng)=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 gwx_ngfft(18)=Information about 3D FFT for the oscillator strengths, see ~abinit/doc/variables/vargs.htm#ngfft gwx_nfftot=number of points of the FFT grid for GW wavefunctions Vcp <vcoul_t datatype> containing information on the cutoff technique %vc_sqrt(npwx,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 representations of the little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states. prtvol=Flags governing verbosity level.
OUTPUT
Sr%x_mat(minbnd:maxbnd,minbnd:maxbnd,%nsppol*Sigp%nsig_ab)=Matrix elements of Sigma_x.
NOTES
1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) is included.
2) 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.
3) 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. Here we divide the states
where the QP energies are required into complexes. Note however that this approach is not
based on group theory, and it might lead to spurious results in case of accidental degeneracies.
PARENTS
sigma
CHILDREN
cwtime,esymm_symmetrize_mels,findqg0,get_bz_item,gsph_fft_tabs hermitianize,littlegroup_print,paw_cross_rho_tw_g,paw_rho_tw_g paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawmknhat_psipsi 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 wfd_paw_get_aeur,wrtout,xmpi_sum
SOURCE
97 #if defined HAVE_CONFIG_H 98 #include "config.h" 99 #endif 100 101 #include "abi_common.h" 102 103 subroutine calc_sigx_me(sigmak_ibz,ikcalc,minbnd,maxbnd,Cryst,QP_BSt,Sigp,Sr,Gsph_x,Vcp,Kmesh,Qmesh,& 104 & Ltg_k,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,allQP_sym,gwx_ngfft,ngfftf,& 105 & prtvol,pawcross) 106 107 use defs_basis 108 use defs_datatypes 109 use m_profiling_abi 110 use m_gwdefs 111 use m_xmpi 112 use m_defs_ptgroups 113 use m_errors 114 use m_time 115 116 use m_blas, only : xdotc, xgemv 117 use m_numeric_tools, only : hermitianize 118 use m_geometry, only : normv 119 use m_crystal, only : crystal_t 120 use m_fft_mesh, only : rotate_FFT_mesh, cigfft 121 use m_bz_mesh, only : kmesh_t, get_BZ_item, findqg0, littlegroup_t, littlegroup_print 122 use m_gsphere, only : gsphere_t, gsph_fft_tabs 123 use m_vcoul, only : vcoul_t 124 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g 125 use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t 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, pawcprj_copy, paw_overlap 130 use m_wfd, only : wfd_t, wfd_get_ur, wfd_get_cprj, wfd_change_ngfft, wfd_paw_get_aeur, wfd_get_many_ur,& 131 & wfd_sym_ur 132 use m_sigma, only : sigma_t, sigma_distribute_bks 133 use m_oscillators, only : rho_tw_g 134 use m_esymm, only : esymm_t, esymm_symmetrize_mels, esymm_failed 135 use m_ptgroups, only : sum_irreps 136 137 !This section has been created automatically by the script Abilint (TD). 138 !Do not modify the following lines by hand. 139 #undef ABI_FUNC 140 #define ABI_FUNC 'calc_sigx_me' 141 use interfaces_14_hidewrite 142 use interfaces_18_timing 143 use interfaces_65_paw 144 !End of the abilint section 145 146 implicit none 147 148 !Arguments ------------------------------------ 149 !scalars 150 integer,intent(in) :: sigmak_ibz,ikcalc,prtvol,minbnd,maxbnd,pawcross 151 type(crystal_t),intent(in) :: Cryst 152 type(ebands_t),target,intent(in) :: QP_BSt 153 type(kmesh_t),intent(in) :: Kmesh,Qmesh 154 type(vcoul_t),intent(in) :: Vcp 155 type(gsphere_t),intent(in) :: Gsph_x 156 type(littlegroup_t),intent(in) :: Ltg_k 157 type(Pseudopotential_type),intent(in) :: Psps 158 type(sigparams_t),target,intent(in) :: Sigp 159 type(sigma_t),intent(inout) :: Sr 160 type(pawang_type),intent(in) :: Pawang 161 type(wfd_t),target,intent(inout) :: Wfd,Wfdf 162 !arrays 163 integer,intent(in) :: gwx_ngfft(18),ngfftf(18) 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 type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw) 168 type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom*Psps%usepaw) 169 170 !Local variables ------------------------------ 171 !scalars 172 integer,parameter :: tim_fourdp=2,ndat1=1 173 integer,parameter :: use_pawnhat=0,ider0=0 174 integer :: gwcalctyp,izero,iab,ib_sum,ib,ib1,ib2,ierr,ig,ig_rot,ii,iik,itim_q,i2 175 integer :: ik_bz,ik_ibz,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx 176 integer :: jik,jk_bz,jk_ibz,kb,nspinor,nsppol,ifft 177 integer :: nq_summed,ibsp,dimcprj_gw,dim_rtwg 178 integer :: spad,spadx1,spadx2,irow,npw_k,ndegs,wtqm,wtqp,my_nbks 179 integer :: isym_kgw,isym_ki,gwx_mgfft,use_padfft,use_padfftf,gwx_fftalga,gwx_fftalgb 180 integer :: gwx_nfftot,nfftf,mgfftf,nhat12_grdim,npwx 181 real(dp) :: cpu_time,wall_time,gflops 182 real(dp) :: fact_sp,theta_mu_minus_esum,tol_empty 183 complex(dpc) :: ctmp,ph_mkgwt,ph_mkt 184 complex(gwpc) :: gwpc_sigxme 185 logical :: iscompatibleFFT,q_is_gamma 186 character(len=500) :: msg 187 !arrays 188 integer :: g0(3),spinor_padx(2,4) 189 integer,allocatable :: igfftxg0(:),igfftfxg0(:) 190 integer,allocatable :: degtab(:,:,:) 191 integer,allocatable :: gwx_gfft(:,:),gwx_gbound(:,:),gboundf(:,:) 192 integer,allocatable :: ktabr(:,:),irottb(:,:),ktabrf(:,:) 193 integer,allocatable :: proc_distrb(:,:,:) 194 real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),qbz(3),q0(3),tsec(2) 195 real(dp) :: spinrot_kbz(4),spinrot_kgw(4) 196 real(dp),ABI_CONTIGUOUS pointer :: qp_ene(:,:,:),qp_occ(:,:,:) 197 real(dp),allocatable :: nhat12(:,:,:),grnhat12(:,:,:,:) 198 complex(gwpc),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:) 199 complex(gwpc),allocatable :: rhotwg_ki(:,:) 200 complex(gwpc),allocatable :: wfr_bdgw(:,:),ur_ibz(:) 201 complex(gwpc),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:) 202 complex(gwpc),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:) 203 complex(dpc),allocatable :: sigxme_tmp(:,:,:),sym_sigx(:,:,:),sigx(:,:,:,:) 204 complex(gwpc),ABI_CONTIGUOUS pointer :: cg_jb(:),cg_sum(:) 205 logical :: can_symmetrize(Wfd%nsppol) 206 logical,allocatable :: bks_mask(:,:,:) 207 type(esymm_t),pointer :: QP_sym(:) 208 type(sigijtab_t),pointer :: Sigxij_tab(:) 209 type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:) 210 type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:) 211 212 !************************************************************************ 213 214 DBG_ENTER("COLL") 215 216 call timab(430,1,tsec) ! csigme (SigX) 217 call cwtime(cpu_time,wall_time,gflops,"start") 218 219 ! Initialize some values. 220 gwcalctyp=Sigp%gwcalctyp 221 nspinor = Wfd%nspinor; nsppol = Wfd%nsppol; npwx = sigp%npwx 222 dim_rtwg = 1; if (nspinor == 2) dim_rtwg = 2 223 spinor_padx = RESHAPE([0, 0, npwx, npwx, 0, npwx, npwx, 0], [2, 4]) 224 ABI_CHECK(Sigp%npwx==Gsph_x%ng,"") 225 226 qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ 227 228 ! Exctract the symmetries of the bands for this k-point 229 QP_sym => allQP_sym(sigmak_ibz,1:nsppol) 230 231 ib1=minbnd; ib2=maxbnd 232 233 ! Index of the GW point in the BZ array, its image in IBZ and time-reversal 234 jk_bz = Sigp%kptgw2bz(ikcalc) 235 call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt) 236 237 ! TODO: the new version based of get_uug won't support kptgw vectors that are not in 238 ! the IBZ since one should perform the rotation before entering the band loop 239 ! In the old version, the rotation was done in rho_tw_g 240 !ABI_CHECK(jik==1,"jik!=1") 241 !ABI_CHECK(isym_kgw==1,"isym_kgw!=1") 242 !ABI_CHECK((ABS(ph_mkgwt - cone) < tol12),"ph_mkgwt!") 243 !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk) 244 spinrot_kgw(:)=Cryst%spinrot(:,isym_kgw) 245 write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,& 246 & ' Calculating <nk|Sigma_x|nk> at k= ',kgw,ch10,& 247 & ' bands from ',ib1,' to ',ib2,ch10 248 call wrtout(std_out,msg,'COLL') 249 250 if (ANY(gwx_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,gwx_ngfft) 251 gwx_mgfft = MAXVAL(gwx_ngfft(1:3)) 252 gwx_fftalga = gwx_ngfft(7)/100; gwx_fftalgb = MOD(gwx_ngfft(7),100)/10 253 254 if (pawcross==1) mgfftf = MAXVAL(ngfftf(1:3)) 255 256 can_symmetrize = .FALSE. 257 if (Sigp%symsigma>0) then 258 can_symmetrize = .TRUE. 259 if (gwcalctyp >= 20) then 260 do spin=1,Wfd%nsppol 261 can_symmetrize(spin) = .not. esymm_failed(QP_sym(spin)) 262 if (.not.can_symmetrize(spin)) then 263 write(msg,'(a,i0,4a)')& 264 & "Symmetrization cannot be performed for spin: ",spin,ch10,& 265 & "band classification encountered the following problem: ",ch10,& 266 & TRIM(QP_sym(spin)%err_msg) 267 MSG_WARNING(msg) 268 end if 269 end do 270 end if 271 if (nspinor == 2) then 272 MSG_WARNING('Symmetrization with nspinor=2 not implemented') 273 end if 274 end if 275 276 ABI_MALLOC(rhotwg_ki, (npwx*nspinor, minbnd:maxbnd)) 277 rhotwg_ki=czero_gw 278 ABI_MALLOC(rhotwg, (npwx*nspinor)) 279 ABI_MALLOC(rhotwgp, (npwx*nspinor)) 280 ABI_MALLOC(vc_sqrt_qbz, (npwx)) 281 282 ! Normalization of theta_mu_minus_esum 283 ! If nsppol==2, qp_occ $\in [0,1]$ 284 SELECT CASE (nsppol) 285 CASE (1) 286 fact_sp=half; tol_empty=0.01 ! below this value the state is assumed empty 287 if (Sigp%nspinor==2) then 288 fact_sp=one; tol_empty=0.005 ! below this value the state is assumed empty 289 end if 290 CASE (2) 291 fact_sp=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic 292 CASE DEFAULT ! spin unpolarized system is treated using nsppol==2 293 MSG_BUG('Wrong nsppol') 294 END SELECT 295 296 ! Table for \Sigmax_ij matrix elements. 297 Sigxij_tab => Sigp%Sigxij_tab(ikcalc,1:nsppol) 298 299 ! Remove empty states from the list of states that will be distributed. 300 ABI_MALLOC(bks_mask,(Wfd%mband,Kmesh%nbz,nsppol)) 301 bks_mask=.FALSE. 302 do spin=1,nsppol 303 do ik_bz=1,Kmesh%nbz 304 ik_ibz = Kmesh%tab(ik_bz) 305 do ib_sum=1,Sigp%nbnds 306 bks_mask(ib_sum,ik_bz,spin) = qp_occ(ib_sum,ik_ibz,spin) >= tol_empty 307 end do 308 end do 309 end do 310 311 ! Distribute tasks. 312 ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,nsppol)) 313 call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,proc_distrb,bks_mask=bks_mask) 314 ABI_FREE(bks_mask) 315 write(msg,'(2(a,i4),a)')" Will sum ",my_nbks ," (b, k, s) occupied states in Sigma_x." 316 call wrtout(std_out,msg) 317 318 ! The index of G-G0 in the FFT mesh the oscillators 319 ! Sigp%mG0 gives the MAX G0 component to account for umklapp. 320 ! Note the size MAX(npwx, Sigp%npwc). 321 ABI_MALLOC(igfftxg0, (Gsph_x%ng)) 322 323 ! Precalculate the FFT index of $ R^{-1}(r-\tau)$ 324 ! S = \transpose R^{-1} and k_BZ = S k_IBZ 325 ! irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk. 326 gwx_nfftot = PRODUCT(gwx_ngfft(1:3)) 327 ABI_MALLOC(irottb,(gwx_nfftot,Cryst%nsym)) 328 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwx_ngfft,irottb,iscompatibleFFT) 329 if (.not. iscompatibleFFT) then 330 MSG_WARNING("FFT mesh is not compatible with symmetries. Results might be affected by large errors!") 331 end if 332 333 ABI_MALLOC(ktabr,(gwx_nfftot,Kmesh%nbz)) 334 do ik_bz=1,Kmesh%nbz 335 isym=Kmesh%tabo(ik_bz) 336 do ifft=1,gwx_nfftot 337 ktabr(ifft,ik_bz)=irottb(ifft,isym) 338 end do 339 end do 340 ABI_FREE(irottb) 341 342 if (Psps%usepaw==1 .and. pawcross==1) then 343 nfftf = PRODUCT(ngfftf(1:3)) 344 ABI_MALLOC(irottb,(nfftf,Cryst%nsym)) 345 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfftf,irottb,iscompatibleFFT) 346 347 ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz)) 348 do ik_bz=1,Kmesh%nbz 349 isym=Kmesh%tabo(ik_bz) 350 do ifft=1,nfftf 351 ktabrf(ifft,ik_bz)=irottb(ifft,isym) 352 end do 353 end do 354 ABI_FREE(irottb) 355 end if 356 357 ! Additional allocations for PAW. 358 if (Psps%usepaw==1) then 359 ABI_DT_MALLOC(Cprj_ksum,(Cryst%natom,nspinor)) 360 call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm) 361 362 nhat12_grdim=0 363 if (use_pawnhat==1) then 364 ! Compensation charge for \phi_a^*\phi_b 365 call wrtout(std_out,"Using nhat12","COLL") 366 ABI_MALLOC(nhat12 ,(2,gwx_nfftot,nspinor**2)) 367 ABI_MALLOC(grnhat12,(2,gwx_nfftot,nspinor**2,3*nhat12_grdim)) 368 end if 369 end if 370 371 ABI_MALLOC(sigxme_tmp, (minbnd:maxbnd, minbnd:maxbnd, Wfd%nsppol*Sigp%nsig_ab)) 372 ABI_MALLOC(sigx,(2, ib1:ib2, ib1:ib2, nsppol*Sigp%nsig_ab)) 373 sigxme_tmp = czero; sigx = czero 374 375 nq_summed=Kmesh%nbz 376 if (Sigp%symsigma > 0) then 377 call littlegroup_print(Ltg_k,std_out,prtvol,'COLL') 378 nq_summed=SUM(Ltg_k%ibzq(:)) 379 380 ! Find number of degenerates subspaces and number of bands in each subspace. 381 ! The tolerance is a little bit arbitrary (0.001 eV) 382 ! It could be reduced, in particular in case of nearly accidental degeneracies 383 ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,nsppol)) 384 degtab=0 385 do spin=1,nsppol 386 do ib=ib1,ib2 387 do jb=ib1,ib2 388 if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin))<0.001/Ha_ev) degtab(ib,jb,spin)=1 389 end do 390 end do 391 end do 392 end if !symsigma 393 394 write(msg,'(2a,i0,a)')ch10,' calc_sigx_me: calculation status (', nq_summed, ' to be completed):' 395 call wrtout(std_out,msg,'COLL') 396 397 ABI_MALLOC(ur_ibz,(gwx_nfftot*nspinor)) 398 if (pawcross==1) then 399 ABI_MALLOC(ur_ae_sum,(nfftf*nspinor)) 400 ABI_MALLOC(ur_ae_onsite_sum,(nfftf*nspinor)) 401 ABI_MALLOC(ur_ps_onsite_sum,(nfftf*nspinor)) 402 end if 403 404 ! ======================================= 405 ! ==== Begin loop over k_i in the BZ ==== 406 ! ======================================= 407 do spin=1,nsppol 408 if (ALL(proc_distrb(:,:,spin) /= Wfd%my_rank)) CYCLE 409 410 ! Load wavefunctions for GW corrections. 411 ABI_STAT_MALLOC(wfr_bdgw, (gwx_nfftot*nspinor,ib1:ib2), ierr) 412 ABI_CHECK(ierr==0, "Out of memory in wfr_bdgw") 413 call wfd_get_many_ur(Wfd, [(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw) 414 415 if (Wfd%usepaw==1) then 416 ! Load cprj for GW states, note the indexing. 417 dimcprj_gw=nspinor*(ib2-ib1+1) 418 ABI_DT_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1)) 419 call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm) 420 ibsp=ib1 421 do jb=ib1,ib2 422 call wfd_get_cprj(Wfd,jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 423 call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 424 call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1))) 425 ibsp=ibsp+nspinor 426 end do 427 if (pawcross==1) then 428 ABI_MALLOC(ur_ae_bdgw,(nfftf*nspinor,ib1:ib2)) 429 ABI_MALLOC(ur_ae_onsite_bdgw,(nfftf*nspinor,ib1:ib2)) 430 ABI_MALLOC(ur_ps_onsite_bdgw,(nfftf*nspinor,ib1:ib2)) 431 do jb=ib1,ib2 432 call wfd_paw_get_aeur(Wfdf,jb,jk_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,& 433 & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum) 434 ur_ae_bdgw(:,jb)=ur_ae_sum 435 ur_ae_onsite_bdgw(:,jb)=ur_ae_onsite_sum 436 ur_ps_onsite_bdgw(:,jb)=ur_ps_onsite_sum 437 end do 438 end if 439 end if 440 441 do ik_bz=1,Kmesh%nbz 442 ! Parallelization over k-points and spin. 443 ! For the spin there is another check in the inner loop 444 if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE 445 446 ! Find the corresponding irreducible k-point 447 call get_BZ_item(Kmesh,ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt) 448 spinrot_kbz = Cryst%spinrot(:,isym_ki) 449 450 ! Identify q and G0 where q+G0=k_GW-k_i 451 kgw_m_ksum = kgw - ksum 452 call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0) 453 454 ! Symmetrize the matrix elements === 455 ! Sum only q"s in IBZ_k. In this case elements are weighted 456 ! according to wtqp and wtqm. wtqm is for time-reversal. 457 wtqp=1; wtqm=0 458 if (can_symmetrize(spin)) then 459 if (Ltg_k%ibzq(iq_bz)/=1) CYCLE 460 wtqp=0; wtqm=0 461 do isym=1,Ltg_k%nsym_sg 462 wtqp = wtqp + Ltg_k%wtksym(1,isym,iq_bz) 463 wtqm = wtqm + Ltg_k%wtksym(2,isym,iq_bz) 464 end do 465 end if 466 467 write(msg,'(2(a,i4),a,i3)')' calc_sigx_me: ik_bz ',ik_bz,'/',Kmesh%nbz,' done by mpi-rank: ',Wfd%my_rank 468 call wrtout(std_out,msg,'PERS') 469 470 ! Find the corresponding irreducible q-point. 471 call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q) 472 q_is_gamma = (normv(qbz,Cryst%gmet,"G") < GW_TOL_W0) 473 474 ! Tables for the FFT of the oscillators. 475 ! a) FFT index of the G-G0. 476 ! b) gwx_gbound table for the zero-padded FFT performed in rhotwg. 477 ABI_MALLOC(gwx_gbound,(2*gwx_mgfft+8,2)) 478 call gsph_fft_tabs(Gsph_x,g0,gwx_mgfft,gwx_ngfft,use_padfft,gwx_gbound,igfftxg0) 479 480 if (ANY(gwx_fftalga == [2, 4])) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 481 !use_padfft=0 482 if (use_padfft==0) then 483 ABI_FREE(gwx_gbound) 484 ABI_MALLOC(gwx_gbound,(2*gwx_mgfft+8,2*use_padfft)) 485 end if 486 487 if (pawcross==1) then 488 ABI_MALLOC(gboundf,(2*mgfftf+8,2)) 489 ABI_MALLOC(igfftfxg0,(Gsph_x%ng)) 490 call gsph_fft_tabs(Gsph_x,g0,mgfftf,ngfftf,use_padfftf,gboundf,igfftfxg0) 491 if ( ANY(gwx_fftalga == [2, 4]) ) use_padfftf=0 492 if (use_padfftf==0) then 493 ABI_FREE(gboundf) 494 ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf)) 495 end if 496 end if 497 498 if (Psps%usepaw==1 .and. use_pawnhat==0) then 499 ! Evaluate oscillator matrix elements 500 ! $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form 501 q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT 502 ABI_DT_MALLOC(Pwij_qg,(Psps%ntypat)) 503 call pawpwij_init(Pwij_qg, npwx, q0, Gsph_x%gvec, Cryst%rprimd, Psps, Pawtab, Paw_pwff) 504 end if 505 506 ! Get Fourier components of the Coulomb interaction in the BZ 507 ! In 3D systems, neglecting umklapp, vc(Sq,sG)=vc(q,G)=4pi/|q+G| 508 ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S. 509 do ig=1,npwx 510 ig_rot = Gsph_x%rottb(ig,itim_q,isym_q) 511 vc_sqrt_qbz(ig_rot)=Vcp%vc_sqrt_resid(ig,iq_ibz) 512 end do 513 514 ! Sum over (occupied) bands. 515 do ib_sum=1,Sigp%nbnds 516 ! Parallelism over bands 517 ! This processor has this k-point but what about spin? 518 if (proc_distrb(ib_sum,ik_bz,spin)/=Wfd%my_rank) CYCLE 519 520 ! Skip empty states. 521 if (qp_occ(ib_sum,ik_ibz,spin)<tol_empty) CYCLE 522 523 call wfd_get_ur(Wfd,ib_sum,ik_ibz,spin,ur_ibz) 524 525 if (Psps%usepaw==1) then 526 ! Load cprj for point ksum, this spin or spinor and *THIS* band. 527 ! TODO MG I could avoid doing this but I have to exchange spin and bands ??? 528 ! For sure there is a better way to do this! 529 call wfd_get_cprj(Wfd,ib_sum,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 530 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 531 if (pawcross==1) then 532 call wfd_paw_get_aeur(Wfdf,ib_sum,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,& 533 & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum) 534 end if 535 end if 536 537 ! Get all <k-q,ib_sum,s|e^{-i(q+G).r}|s,jb,k> 538 do jb=ib1,ib2 539 540 if (Psps%usepaw==1 .and. use_pawnhat==1) then 541 MSG_ERROR("use_pawnhat is disabled") 542 i2=jb; if (nspinor==2) i2=(2*jb-1) 543 spad=(nspinor-1) 544 545 izero=0 546 call pawmknhat_psipsi(Cprj_ksum,Cprj_kgw(:,i2:i2+spad),ider0,izero,Cryst%natom,& 547 & Cryst%natom,gwx_nfftot,gwx_ngfft,nhat12_grdim,nspinor,Cryst%ntypat,Pawang,Pawfgrtab,& 548 & grnhat12,nhat12,pawtab) 549 550 else 551 call rho_tw_g(nspinor,npwx,gwx_nfftot,ndat1,gwx_ngfft,1,use_padfft,igfftxg0,gwx_gbound, & 552 & ur_ibz ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, & 553 & wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw, & 554 & nspinor,rhotwg_ki(:,jb)) 555 556 if (Psps%usepaw==1.and.use_pawnhat==0) then 557 ! Add on-site contribution, projectors are already in BZ. 558 i2=jb; if (nspinor==2) i2=(2*jb-1) 559 spad=(nspinor-1) 560 call paw_rho_tw_g(npwx,nspinor,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_x%gvec,& 561 & Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),Pwij_qg,rhotwg_ki(:,jb)) 562 end if 563 if (Psps%usepaw==1.and.pawcross==1) then ! Add paw cross term 564 call paw_cross_rho_tw_g(nspinor,npwx,nfftf,ngfftf,1,use_padfftf,igfftfxg0,gboundf,& 565 & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum,iik,ktabrf(:,ik_bz),ph_mkt,spinrot_kbz,& 566 & ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,& 567 & nspinor,rhotwg_ki(:,jb)) 568 end if 569 end if 570 571 ! Multiply by the square root of the Coulomb term 572 ! In 3-D systems, the factor sqrt(4pi) is included 573 do ii=1,nspinor 574 spad = (ii-1) * npwx 575 rhotwg_ki(spad+1:spad+npwx,jb) = rhotwg_ki(spad+1:spad + npwx,jb) * vc_sqrt_qbz(1:npwx) 576 end do 577 578 if (ik_bz == jk_bz) then 579 ! Treat analytically the case q --> 0 580 ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma 581 ! while the Colulomb term is integrated out. 582 ! * In the scalar case we have nonzero contribution only if ib_sum==jb 583 ! * For nspinor==2 evalute <ib_sum,up|jb,up> and <ib_sum,dwn|jb,dwn>, 584 ! impose orthonormalization since npwwfn might be < npwvec. 585 586 if (nspinor==1) then 587 rhotwg_ki(1,jb) = czero_gw 588 !Note the use of i_sz_resid and not i_sz, to account for the possibility to have generalized KS basis set from hybrid 589 if (ib_sum == jb) rhotwg_ki(1,jb) = CMPLX(SQRT(Vcp%i_sz_resid), 0.0_gwp) 590 !rhotwg_ki(1,jb) = czero_gw ! DEBUG 591 else 592 npw_k = Wfd%npwarr(ik_ibz) 593 rhotwg_ki(1, jb) = zero; rhotwg_ki(npwx+1, jb) = zero 594 if (ib_sum == jb) then 595 cg_sum => Wfd%Wave(ib_sum,ik_ibz,spin)%ug 596 cg_jb => Wfd%Wave(jb,jk_ibz,spin)%ug 597 ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1) 598 rhotwg_ki(1, jb) = CMPLX(SQRT(Vcp%i_sz_resid), 0.0_gwp) * real(ctmp) 599 ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1) 600 rhotwg_ki(npwx+1, jb) = CMPLX(SQRT(Vcp%i_sz_resid), 0.0_gwp) * real(ctmp) 601 end if 602 !rhotwg_ki(1, jb) = zero; rhotwg_ki(npwx+1, jb) = zero 603 ! PAW is missing 604 end if 605 end if 606 end do ! jb Got all matrix elements from minbnd up to maxbnd. 607 608 theta_mu_minus_esum = fact_sp * qp_occ(ib_sum,ik_ibz,spin) 609 610 do kb=ib1,ib2 611 ! Compute the ket Sigma_x |phi_{k,kb}>. 612 rhotwgp(:) = rhotwg_ki(:,kb) 613 614 ! Loop over the non-zero row elements of this column. 615 ! If gwcalctyp < 20: only diagonal elements since QP == KS. 616 ! If gwcalctyp >= 20: 617 ! * Only off-diagonal elements connecting states with same character. 618 ! * Only the upper triangle if HF, SEX, or COHSEX. 619 do irow=1,Sigxij_tab(spin)%col(kb)%size1 620 jb = Sigxij_tab(spin)%col(kb)%bidx(irow) 621 rhotwg = rhotwg_ki(:,jb) 622 623 ! Calculate bare exchange <phi_j|Sigma_x|phi_k>. 624 ! Do the scalar product only if ib_sum is occupied. 625 if (theta_mu_minus_esum/fact_sp >= tol_empty) then 626 do iab=1,Sigp%nsig_ab 627 spadx1 = spinor_padx(1, iab); spadx2 = spinor_padx(2, iab) 628 gwpc_sigxme = -XDOTC(npwx, rhotwg(spadx1+1:), 1, rhotwgp(spadx2+1:), 1) * theta_mu_minus_esum 629 630 ! Accumulate and symmetrize Sigma_x 631 ! -wtqm comes from time-reversal (exchange of band indeces) 632 is_idx = spin; if (nspinor == 2) is_idx = iab 633 sigxme_tmp(jb, kb, is_idx) = sigxme_tmp(jb, kb, is_idx) + & 634 & (wtqp + wtqm)*DBLE(gwpc_sigxme) + (wtqp - wtqm)*j_gw*AIMAG(gwpc_sigxme) 635 636 sigx(1, jb, kb, is_idx) = sigx(1, jb, kb, is_idx) + wtqp * gwpc_sigxme 637 sigx(2, jb, kb, is_idx) = sigx(2, jb, kb, is_idx) + wtqm *CONJG(gwpc_sigxme) 638 end do 639 end if 640 end do ! jb used to calculate matrix elements of Sigma_x 641 642 end do ! kb to calculate matrix elements of Sigma_x 643 end do ! ib_sum 644 645 ! Deallocate k-dependent quantities. 646 ABI_FREE(gwx_gbound) 647 if (pawcross==1) then 648 ABI_FREE(gboundf) 649 end if 650 651 if (Psps%usepaw==1.and.use_pawnhat==0) then 652 call pawpwij_free(Pwij_qg) 653 ABI_DT_FREE(Pwij_qg) 654 end if 655 end do ! ik_bz Got all diagonal (off-diagonal) matrix elements. 656 657 ABI_FREE(wfr_bdgw) 658 if (Wfd%usepaw==1) then 659 call pawcprj_free(Cprj_kgw) 660 ABI_DT_FREE(Cprj_kgw) 661 if (pawcross==1) then 662 ABI_FREE(ur_ae_bdgw) 663 ABI_FREE(ur_ae_onsite_bdgw) 664 ABI_FREE(ur_ps_onsite_bdgw) 665 end if 666 end if 667 end do !spin 668 669 ABI_FREE(igfftxg0) 670 if (pawcross==1) then 671 ABI_FREE(igfftfxg0) 672 end if 673 674 ! Gather contributions from all the CPUs. 675 call xmpi_sum(sigxme_tmp, wfd%comm, ierr) 676 call xmpi_sum(sigx, wfd%comm, ierr) 677 678 ! Multiply by constants. For 3D systems sqrt(4pi) is included in vc_sqrt_qbz. 679 sigxme_tmp = (one/(Cryst%ucvol*Kmesh%nbz)) * sigxme_tmp * Sigp%sigma_mixing 680 sigx = (one/(Cryst%ucvol*Kmesh%nbz)) * sigx * Sigp%sigma_mixing 681 682 ! 683 ! If we have summed over the IBZ_q now we have to average over degenerate states. 684 ! NOTE: Presently only diagonal terms are considered 685 ! TODO QP-SCGW required a more involved approach, there is a check in sigma 686 ! TODO it does not work if spinor==2. 687 do spin=1,nsppol 688 if (can_symmetrize(spin)) then 689 ABI_MALLOC(sym_sigx, (ib1:ib2, ib1:ib2, sigp%nsig_ab)) 690 sym_sigx = czero 691 692 ! Average over degenerate diagonal elements. 693 ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results. 694 ! another good reason to use a strict criterion for the tollerance on eigenvalues. 695 do ib=ib1,ib2 696 ndegs=0 697 do jb=ib1,ib2 698 if (degtab(ib,jb,spin)==1) then 699 if (nspinor == 1) then 700 sym_sigx(ib, ib, 1) = sym_sigx(ib, ib, 1) + SUM(sigx(:,jb,jb,spin)) 701 else 702 do ii=1,sigp%nsig_ab 703 sym_sigx(ib, ib, ii) = sym_sigx(ib, ib, ii) + SUM(sigx(:,jb,jb,ii)) 704 end do 705 end if 706 end if 707 ndegs = ndegs + degtab(ib,jb,spin) 708 end do 709 sym_sigx(ib,ib,:) = sym_sigx(ib,ib,:) / ndegs 710 end do 711 712 if (gwcalctyp >= 20) then 713 call esymm_symmetrize_mels(QP_sym(spin),ib1,ib2,sigx(:,:,:,spin),sym_sigx(:,:,1)) 714 end if 715 716 ! Copy symmetrized values. 717 do ib=ib1,ib2 718 do jb=ib1,ib2 719 if (nspinor == 1) then 720 sigxme_tmp(ib,jb,spin) = sym_sigx(ib,jb,1) 721 else 722 do ii=1,sigp%nsig_ab 723 sigxme_tmp(ib,jb,ii) = sym_sigx(ib,jb,ii) 724 end do 725 end if 726 end do 727 end do 728 729 ABI_FREE(sym_sigx) 730 end if 731 end do 732 733 if (gwcalctyp>=20) then 734 ! Reconstruct the full sigma_x matrix from the upper triangle. 735 if (nspinor == 1) then 736 do spin=1,nsppol 737 call hermitianize(sigxme_tmp(:,:,spin), "Upper") 738 end do 739 else 740 MSG_WARNING("Should hermitianize non-collinear sigma!") 741 end if 742 end if 743 744 ! Save diagonal elements or ab components of Sigma_x (hermitian) 745 ! TODO It should be hermitian also if nspinor==2 746 do spin=1,nsppol 747 do jb=ib1,ib2 748 do iab=1,Sigp%nsig_ab 749 is_idx = spin; if (Sigp%nsig_ab > 1) is_idx = iab 750 if (is_idx <= 2) then 751 Sr%sigxme(jb,jk_ibz,is_idx) = DBLE(sigxme_tmp(jb,jb,is_idx)) 752 else 753 Sr%sigxme(jb,jk_ibz,is_idx) = sigxme_tmp(jb,jb,is_idx) 754 end if 755 end do 756 !if (Sigp%nsig_ab>1) then 757 ! write(std_out,'(i3,4f8.3,a,f8.3)')jb,Sr%sigxme(jb,jk_ibz,:)*Ha_eV,' Tot ',SUM(Sr%sigxme(jb,jk_ibz,:))*Ha_eV 758 !end if 759 end do 760 end do 761 762 ! Save full exchange matrix in Sr% 763 Sr%x_mat(minbnd:maxbnd,minbnd:maxbnd,jk_ibz,:) = sigxme_tmp(minbnd:maxbnd,minbnd:maxbnd,:) 764 ABI_FREE(sigxme_tmp) 765 766 ! =========================== 767 ! ==== Deallocate memory ==== 768 ! =========================== 769 if (Psps%usepaw==1) then 770 if (allocated(gwx_gfft)) then 771 ABI_FREE(gwx_gfft) 772 end if 773 call pawcprj_free(Cprj_ksum) 774 ABI_DT_FREE(Cprj_ksum) 775 if (allocated(Pwij_fft)) then 776 call pawpwij_free(Pwij_fft) 777 ABI_DT_FREE(Pwij_fft) 778 end if 779 if (use_pawnhat==1) then 780 ABI_FREE(nhat12) 781 ABI_FREE(grnhat12) 782 end if 783 if (pawcross==1) then 784 ABI_FREE(ur_ae_sum) 785 ABI_FREE(ur_ae_onsite_sum) 786 ABI_FREE(ur_ps_onsite_sum) 787 ABI_FREE(ktabrf) 788 end if 789 end if 790 791 ABI_FREE(ur_ibz) 792 ABI_FREE(rhotwg_ki) 793 ABI_FREE(rhotwg) 794 ABI_FREE(rhotwgp) 795 ABI_FREE(vc_sqrt_qbz) 796 ABI_FREE(ktabr) 797 ABI_FREE(sigx) 798 ABI_FREE(proc_distrb) 799 800 if (allocated(degtab)) then 801 ABI_FREE(degtab) 802 end if 803 804 call timab(430,2,tsec) ! csigme (SigX) 805 call cwtime(cpu_time,wall_time,gflops,"stop") 806 write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time 807 808 DBG_EXIT("COLL") 809 810 end subroutine calc_sigx_me