TABLE OF CONTENTS
ABINIT/calc_sigc_me [ Functions ]
NAME
calc_sigc_me
FUNCTION
Calculate diagonal and off-diagonal matrix elements 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) Dtset <type(dataset_type)>=all input variables in this dataset %iomode %paral_kgb %nspinor=Number of spinorial components. %gwcomp=If 1 use an extrapolar approximation to accelerate convergence. %gwencomp=Extrapolar energy. 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 much slower but it requires less memory %nomega_i=Number of imaginary frequencies. %nomega_r=Number of real frequencies. %nomega=Total number of frequencies. Gsph_c<gsphere_t>= info on G-sphere for Sigma_c Gsph_Max<gsphere_t>= info on biggest G-sphere %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 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 PPm<ppmodel_t>= Datatype gathering information on the Plasmonpole technique (see also ppm_get_qbz). %model=type of ppmodel %npwc=number of G-vectors for the correlation part. %dm2_otq =size of second dimension of otq array %dm2_bots=size of second dimension of botsq arrays %dm_eig =size of second dimension of eig arrays QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,Kmesh%nibz,Wfd%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(Wfd%nkibz,Wfd%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) use_aerhor=1 is aepaw_rhor is used, 0 otherwise. aepaw_rhor(rho_nfftot,Wfd%nspden*use_aerhor)=AE PAW density used to generate PPmodel paramenters if mqmem==0
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_comp,calc_sig_ppm,calc_sig_ppm_comp,calc_sigc_cd,calc_wfwfg coeffs_gausslegint,cwtime,epsm1_symmetrizer,epsm1_symmetrizer_inplace esymm_symmetrize_mels,findqg0,get_bz_item,get_epsm1,get_gftt gsph_fft_tabs,littlegroup_print,paw_cross_rho_tw_g,paw_rho_tw_g paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawpwij_free pawpwij_init,ppm_get_qbz,rho_tw_g,rotate_fft_mesh,setup_ppmodel sigma_distribute_bks,timab,wfd_change_ngfft,wfd_get_cprj wfd_get_many_ur,wfd_get_ur,wfd_paw_get_aeur,wrtout,xmpi_max,xmpi_sum
SOURCE
115 #if defined HAVE_CONFIG_H 116 #include "config.h" 117 #endif 118 119 #include "abi_common.h" 120 121 122 subroutine calc_sigc_me(sigmak_ibz,ikcalc,nomega_sigc,minbnd,maxbnd,& 123 & Dtset,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,Ltg_k,& 124 & PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,allQP_sym,& 125 & gwc_ngfft,rho_ngfft,rho_nfftot,rhor,use_aerhor,aepaw_rhor,sigcme_tmp) 126 127 use defs_basis 128 use defs_datatypes 129 use defs_abitypes 130 use m_gwdefs 131 use m_profiling_abi 132 use m_xmpi 133 use m_defs_ptgroups 134 use m_errors 135 use m_time 136 137 use m_blas, only : xdotc, xgemv 138 use m_numeric_tools, only : hermitianize, imin_loc, coeffs_gausslegint 139 use m_fstrings, only : sjoin, itoa 140 use m_geometry, only : normv 141 use m_crystal, only : crystal_t 142 use m_bz_mesh, only : kmesh_t, get_BZ_item, findqg0, littlegroup_t, littlegroup_print 143 use m_gsphere, only : gsphere_t, gsph_fft_tabs 144 use m_fft_mesh, only : get_gftt, rotate_fft_mesh, cigfft 145 use m_vcoul, only : vcoul_t 146 use m_wfd, only : wfd_get_ur, wfd_t, wfd_get_cprj, wfd_barrier, wfd_change_ngfft, wfd_paw_get_aeur, & 147 & wfd_get_many_ur,wfd_sym_ur 148 use m_oscillators, only : rho_tw_g, calc_wfwfg 149 use m_screening, only : epsilonm1_results, epsm1_symmetrizer, epsm1_symmetrizer_inplace, get_epsm1 150 use m_ppmodel, only : setup_ppmodel, ppm_get_qbz, ppmodel_t, calc_sig_ppm 151 use m_sigma, only : sigma_t, sigma_distribute_bks 152 use m_esymm, only : esymm_t, esymm_symmetrize_mels, esymm_failed 153 use m_pawang, only : pawang_type 154 use m_pawtab, only : pawtab_type 155 use m_pawfgrtab, only : pawfgrtab_type 156 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap 157 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g 158 use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t 159 160 !This section has been created automatically by the script Abilint (TD). 161 !Do not modify the following lines by hand. 162 #undef ABI_FUNC 163 #define ABI_FUNC 'calc_sigc_me' 164 use interfaces_14_hidewrite 165 use interfaces_18_timing 166 use interfaces_65_paw 167 use interfaces_70_gw, except_this_one => calc_sigc_me 168 !End of the abilint section 169 170 implicit none 171 172 !Arguments ------------------------------------ 173 !scalars 174 integer,intent(in) :: sigmak_ibz,ikcalc,rho_nfftot,nomega_sigc,minbnd,maxbnd 175 integer,intent(in) :: use_aerhor 176 type(crystal_t),intent(in) :: Cryst 177 type(ebands_t),target,intent(in) :: QP_BSt 178 type(kmesh_t),intent(in) :: Kmesh,Qmesh 179 type(vcoul_t),intent(in) :: Vcp 180 type(dataset_type),intent(in) :: Dtset 181 type(epsilonm1_results),intent(inout) :: Er 182 type(gsphere_t),intent(in) :: Gsph_Max,Gsph_c 183 type(littlegroup_t),intent(in) :: Ltg_k 184 type(ppmodel_t),intent(inout) :: PPm 185 type(Pseudopotential_type),intent(in) :: Psps 186 type(pawang_type),intent(in) :: pawang 187 type(sigparams_t),target,intent(in) :: Sigp 188 type(sigma_t),intent(in) :: Sr 189 type(wfd_t),target,intent(inout) :: Wfd,Wfdf 190 !arrays 191 integer,intent(in) :: gwc_ngfft(18),rho_ngfft(18) 192 real(dp),intent(in) :: rhor(rho_nfftot,Wfd%nspden) 193 real(dp),intent(in) :: aepaw_rhor(rho_nfftot,Wfd%nspden*use_aerhor) 194 complex(dpc),intent(out) :: sigcme_tmp(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Wfd%nsppol*Sigp%nsig_ab) 195 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat) 196 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 197 type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol) 198 type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw) 199 type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom*Psps%usepaw) 200 201 !Local variables ------------------------------ 202 !scalars 203 integer,parameter :: tim_fourdp2=2,ndat1=1 204 integer :: npw_k,iab,ib,ib1,ib2,ierr,ig,iggp,igp,ii,iik,itim_q,i1,i2,npls 205 integer :: ik_bz,ik_ibz,io,iiw,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx 206 integer :: band,band1,band2,idle,rank,jik,jk_bz,jk_ibz,kb,nspinor 207 integer :: nomega_tot,nq_summed,ispinor,ibsp,dimcprj_gw,npwc 208 integer :: spad,spadc,spadc1,spadc2,irow,my_nbks,ndegs,wtqm,wtqp,mod10 209 integer :: isym_kgw,isym_ki,gwc_mgfft,use_padfft,gwc_fftalga,gwc_nfftot,nfftf,mgfftf,use_padfftf 210 integer :: iwc,ifft 211 real(dp) :: cpu_time,wall_time,gflops 212 real(dp) :: e0i,fact_sp,theta_mu_minus_e0i,tol_empty,z2,en_high,gw_gsq,w_localmax,w_max 213 complex(dpc) :: ctmp,omegame0i2_ac,omegame0i_ac,ph_mkgwt,ph_mkt 214 logical :: iscompatibleFFT,q_is_gamma 215 character(len=500) :: msg,sigma_type 216 complex(gwpc),allocatable :: botsq(:,:),otq(:,:),eig(:,:) 217 !arrays 218 integer :: g0(3),spinor_padc(2,4),got(Wfd%nproc) 219 integer,allocatable :: proc_distrb(:,:,:),extrapolar_distrb(:,:,:,:),degtab(:,:,:) 220 integer,allocatable :: igfftcg0(:),gw_gfft(:,:),gw_gbound(:,:),irottb(:,:),ktabr(:,:) 221 integer,allocatable :: igfftfcg0(:),gboundf(:,:),ktabrf(:,:),npoles_missing(:) 222 real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),omegap(Er%nomega_i),omegap2(Er%nomega_i),q0(3),tsec(2),qbz(3) 223 real(dp) :: gl_knots(Er%nomega_i),gl_wts(Er%nomega_i) 224 real(dp) :: spinrot_kbz(4),spinrot_kgw(4) 225 real(dp),ABI_CONTIGUOUS pointer :: qp_ene(:,:,:),qp_occ(:,:,:) 226 real(dp),allocatable :: omegame0i(:), w_maxval(:) 227 complex(gwpc) :: sigcohme(Sigp%nsig_ab) 228 complex(gwpc),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:) 229 complex(gwpc),allocatable :: botsq_conjg_transp(:,:),ac_epsm1cqwz2(:,:,:) 230 complex(gwpc),allocatable :: epsm1_qbz(:,:,:),epsm1_trcc_qbz(:,:,:), epsm1_tmp(:,:) 231 complex(gwpc),allocatable :: ac_integr(:,:,:),sigc_ket(:,:),ket1(:,:),ket2(:,:) 232 complex(gwpc),allocatable :: herm_sigc_ket(:,:),aherm_sigc_ket(:,:) 233 complex(gwpc),allocatable :: rhotwg_ki(:,:) 234 complex(gwpc),allocatable :: sigcme2(:,:),sigcme_3(:),sigcme_new(:),sigctmp(:,:) 235 complex(gwpc),allocatable :: wfr_bdgw(:,:),ur_ibz(:),wf1swf2_g(:),usr_bz(:) 236 complex(gwpc),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:) 237 complex(gwpc),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:) 238 complex(gwpc),allocatable :: otq_transp(:,:) 239 complex(gwpc),ABI_CONTIGUOUS pointer :: cg_jb(:),cg_sum(:) 240 complex(dpc),allocatable :: sym_cme(:,:,:,:),sigc(:,:,:,:,:) 241 logical :: rank_mask(Wfd%nproc),can_symmetrize(Wfd%nsppol) 242 !logical :: me_calc_poles(Sr%nomega_r+Sr%nomega4sd) 243 type(sigijtab_t),pointer :: Sigcij_tab(:) 244 type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:) 245 type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:) 246 type(esymm_t),pointer :: QP_sym(:) 247 248 !************************************************************************ 249 250 DBG_ENTER("COLL") 251 252 ! Initial check 253 ABI_CHECK(Sr%nomega_r==Sigp%nomegasr,"") 254 ABI_CHECK(Sr%nomega4sd==Sigp%nomegasrd,"") 255 ABI_CHECK(Sigp%npwc==Gsph_c%ng,"") 256 ABI_CHECK(Sigp%npwvec==Gsph_Max%ng,"") 257 258 call timab(424,1,tsec) ! calc_sigc_me 259 call timab(431,1,tsec) ! calc_sigc_me 260 call timab(432,1,tsec) ! Init 261 call cwtime(cpu_time,wall_time,gflops,"start") 262 263 ! Initialize MPI variables 264 qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ 265 266 ! Extract the symmetries of the bands for this k-point 267 QP_sym => allQP_sym(sigmak_ibz,1:Wfd%nsppol) 268 269 ! Index of the GW point in the BZ array, its image in IBZ and time-reversal === 270 jk_bz=Sigp%kptgw2bz(ikcalc) 271 call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt) 272 !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk) 273 274 ! TODO: the new version based of get_uug won't suppporte kptgw vectors that are not in 275 ! the IBZ since one should perform the rotation before entering the band loop 276 ! In the old version, the rotation was done in rho_tw_g 277 !ABI_CHECK(jik==1,"jik!=1") 278 !ABI_CHECK(isym_kgw==1,"isym_kgw!=1") 279 !ABI_CHECK((ABS(ph_mkgwt - cone) < tol12),"ph_mkgwt!") 280 281 spinrot_kgw=Cryst%spinrot(:,isym_kgw) 282 ib1=minbnd; ib2=maxbnd 283 284 write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,& 285 & ' Calculating <nk|Sigma_c(omega)|nk> at k = ',kgw(:),ch10,& 286 & ' bands n = from ',ib1,' to ',ib2,ch10 287 call wrtout(std_out,msg,'COLL') 288 289 ABI_MALLOC(w_maxval,(minbnd:maxbnd)) 290 w_maxval = zero 291 292 if ( ANY(gwc_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,gwc_ngfft) 293 gwc_mgfft = MAXVAL(gwc_ngfft(1:3)) 294 gwc_fftalga = gwc_ngfft(7)/100 !; gwc_fftalgc=MOD(gwc_ngfft(7),10) 295 296 if (Dtset%pawcross==1) mgfftf = MAXVAL(rho_ngfft(1:3)) 297 298 can_symmetrize = .FALSE. 299 if (Sigp%symsigma>0) then 300 can_symmetrize = .TRUE. 301 if (Sigp%gwcalctyp >= 20) then 302 do spin=1,Wfd%nsppol 303 can_symmetrize(spin) = .not. esymm_failed(QP_sym(spin)) 304 if (.not.can_symmetrize(spin)) then 305 write(msg,'(a,i0,4a)')& 306 & " Symmetrization cannot be performed for spin: ",spin,ch10,& 307 & " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg) 308 MSG_WARNING(msg) 309 end if 310 end do 311 end if 312 if (wfd%nspinor == 2) MSG_WARNING("Symmetrization with nspinor=2 not implemented") 313 end if 314 315 ABI_UNUSED(Pawang%l_max) 316 317 ! Print type of calculation. 318 mod10=MOD(Sigp%gwcalctyp,10); sigma_type = sigma_type_from_key(mod10) 319 call wrtout(std_out,sigma_type,'COLL') 320 321 ! Set up logical flags for Sigma calculation 322 if (mod10==SIG_GW_AC.and.Sigp%gwcalctyp/=1) then 323 MSG_ERROR("not implemented") 324 end if 325 326 ! Initialize some values 327 nspinor = Wfd%nspinor 328 npwc = sigp%npwc 329 spinor_padc(:,:)=RESHAPE([0, 0, npwc, npwc, 0, npwc, npwc, 0], [2, 4]) 330 ABI_MALLOC(npoles_missing,(minbnd:maxbnd)) 331 npoles_missing=0 332 333 ! Normalization of theta_mu_minus_e0i 334 ! If nsppol==2, qp_occ $\in [0,1]$ 335 SELECT CASE (Wfd%nsppol) 336 CASE (1) 337 fact_sp=half; tol_empty=0.01 ! below this value the state is assumed empty 338 if (Wfd%nspinor==2) then 339 fact_sp=one; tol_empty=0.005 ! below this value the state is assumed empty 340 end if 341 CASE (2) 342 fact_sp=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic 343 CASE DEFAULT ! spin unpolarized system is treated using nsppol==2 344 MSG_BUG('Wrong nsppol') 345 END SELECT 346 347 ! Allocate arrays used to accumulate the matrix elements of \Sigma_c over 348 ! k-points and bands. Note that for AC requires only the imaginary frequencies 349 !nomega_sigc=Sr%nomega_r+Sr%nomega4sd 350 !if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i 351 ! 352 ! === Define the G-G0 shifts for the FFT of the oscillators === 353 ! * Sigp%mG0 gives the MAX G0 component to account for umklapp. 354 ! * Note the size MAX(Sigp%npwx,npwc). 355 ! 356 ! === Precalculate the FFT index of $(R^{-1}(r-\tau))$ === 357 ! * S=\transpose R^{-1} and k_BZ = S k_IBZ 358 ! * irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk. 359 gwc_nfftot = PRODUCT(gwc_ngfft(1:3)) 360 ABI_MALLOC(irottb,(gwc_nfftot,Cryst%nsym)) 361 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwc_ngfft,irottb,iscompatibleFFT) 362 if (.not.iscompatibleFFT) then 363 msg = "FFT mesh is not compatible with symmetries. Results might be affected by large errors!" 364 MSG_WARNING(msg) 365 end if 366 367 ABI_MALLOC(ktabr,(gwc_nfftot,Kmesh%nbz)) 368 do ik_bz=1,Kmesh%nbz 369 isym=Kmesh%tabo(ik_bz) 370 do ifft=1,gwc_nfftot 371 ktabr(ifft,ik_bz)=irottb(ifft,isym) 372 end do 373 end do 374 ABI_FREE(irottb) 375 376 if (Psps%usepaw==1 .and. Dtset%pawcross==1) then 377 nfftf = PRODUCT(rho_ngfft(1:3)) 378 ABI_MALLOC(irottb,(nfftf,Cryst%nsym)) 379 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,rho_ngfft,irottb,iscompatibleFFT) 380 381 ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz)) 382 do ik_bz=1,Kmesh%nbz 383 isym=Kmesh%tabo(ik_bz) 384 do ifft=1,nfftf 385 ktabrf(ifft,ik_bz)=irottb(ifft,isym) 386 end do 387 end do 388 ABI_FREE(irottb) 389 end if 390 391 Sigcij_tab => Sigp%Sigcij_tab(ikcalc,1:Wfd%nsppol) 392 393 got=0 394 ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,Wfd%nsppol)) 395 call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,Wfd%nsppol,can_symmetrize,kgw,Sigp%mg0,& 396 & my_nbks,proc_distrb,got,global=.TRUE.) 397 398 write(msg,'(a,i0,a)')" Will sum ",my_nbks," (b,k,s) states in Sigma_c." 399 call wrtout(std_out,msg,'PERS') 400 401 if (Sigp%gwcomp==1) then 402 en_high=MAXVAL(qp_ene(Sigp%nbnds,:,:)) + Sigp%gwencomp 403 write(msg,'(6a,e11.4,a)')ch10,& 404 & ' Using the extrapolar approximation to accelerate convergence',ch10,& 405 & ' with respect to the number of bands included',ch10,& 406 & ' with extrapolar energy: ',en_high*Ha_eV,' [eV]' 407 call wrtout(std_out,msg,'COLL') 408 ABI_MALLOC(wf1swf2_g,(gwc_nfftot*nspinor)) 409 endif 410 411 if (Sigp%gwcomp == 1) then 412 ! Setup of MPI table for extrapolar contributions. 413 ABI_MALLOC(extrapolar_distrb,(ib1:ib2,ib1:ib2,Kmesh%nbz,Wfd%nsppol)) 414 extrapolar_distrb = xmpi_undefined_rank 415 416 do spin=1,Wfd%nsppol 417 do ik_bz=1,Kmesh%nbz 418 if (ANY(proc_distrb(:,ik_bz,spin) /= xmpi_undefined_rank) ) then ! This BZ point will be calculated. 419 rank_mask = .FALSE. ! The set of node that will treat (k,s). 420 do band=1,Wfd%mband 421 rank = proc_distrb(band,ik_bz,spin) 422 if (rank /= xmpi_undefined_rank) rank_mask(rank+1)=.TRUE. 423 end do 424 do band2=ib1,ib2 425 do irow=1,Sigcij_tab(spin)%col(band2)%size1 ! Looping over the non-zero elements of sigma_ij. 426 band1 = Sigcij_tab(spin)%col(band2)%bidx(irow) 427 idle = imin_loc(got,mask=rank_mask) 428 got(idle) = got(idle)+1 429 extrapolar_distrb(band1,band2,ik_bz,spin) = idle-1 430 end do 431 end do 432 end if 433 end do 434 end do 435 436 write(msg,'(a,i0,a)')" Will treat ",COUNT(extrapolar_distrb==Wfd%my_rank)," extrapolar terms." 437 call wrtout(std_out,msg,'PERS') 438 end if 439 440 ABI_MALLOC(rhotwg_ki, (npwc*nspinor, minbnd:maxbnd)) 441 rhotwg_ki=czero_gw 442 ABI_MALLOC(rhotwg, (npwc*nspinor)) 443 ABI_MALLOC(rhotwgp, (npwc*nspinor)) 444 ABI_MALLOC(vc_sqrt_qbz, (npwc)) 445 446 if (Er%mqmem==0) then ! Use out-of-core solution for epsilon. 447 MSG_COMMENT('Reading q-slices from file. Slower but less memory.') 448 end if ! 449 450 ! Additional allocations for PAW. 451 if (Psps%usepaw==1) then 452 ABI_DT_MALLOC(Cprj_ksum,(Cryst%natom,nspinor)) 453 call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm) 454 ! 455 ! For the extrapolar method we need the onsite terms of the PW in the FT mesh. 456 ! * gw_gfft is the set of plane waves in the FFT Box for the oscillators. 457 if (Sigp%gwcomp==1) then 458 ABI_MALLOC(gw_gfft,(3,gwc_nfftot)) 459 q0=zero 460 call get_gftt(gwc_ngfft,q0,Cryst%gmet,gw_gsq,gw_gfft) 461 ABI_DT_MALLOC(Pwij_fft,(Psps%ntypat)) 462 call pawpwij_init(Pwij_fft,gwc_nfftot,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 463 end if 464 end if ! usepaw==1 465 466 if (mod10==SIG_GW_AC) then ! Calculate Gauss-Legendre quadrature knots and weights for analytic continuation 467 call coeffs_gausslegint(zero,one,gl_knots,gl_wts,Er%nomega_i) 468 469 do io=1,Er%nomega_i ! First frequencies are always real 470 if (ABS(AIMAG(one*Er%omega(Er%nomega_r+io))-(one/gl_knots(io)-one)) > 0.0001) then 471 write(msg,'(3a)')& 472 & ' Frequencies in the SCR file are not compatible with the analytic continuation. ',ch10,& 473 & ' Verify the frequencies in the SCR file. ' 474 MSG_WARNING(msg) 475 if (Wfd%my_rank==Wfd%master) write(std_out,*)AIMAG(Er%omega(Er%nomega_r+io)),(one/gl_knots(io)-one) 476 MSG_ERROR("Cannot continue!") 477 end if 478 end do 479 480 ! To calculate \int_0^\infty domegap f(omegap), we calculate \int_0^1 dz f(1/z-1)/z^2. 481 omegap(:)=one/gl_knots(:)-one 482 omegap2(:)=omegap(:)*omegap(:) 483 ABI_MALLOC(ac_epsm1cqwz2, (npwc, npwc, Er%nomega_i)) 484 ABI_MALLOC(ac_integr, (npwc, npwc, Sr%nomega_i)) 485 end if 486 487 ! Calculate total number of frequencies and allocate related arrays. 488 ! sigcme2 is used to accumulate the diagonal matrix elements over k-points and 489 ! GW bands, used only in case of ppmodel 3 and 4 (TODO save memory) 490 nomega_tot=Sr%nomega_r+Sr%nomega4sd 491 ABI_MALLOC(sigcme2,(nomega_tot,ib1:ib2)) 492 ABI_MALLOC(sigcme_3,(nomega_tot)) 493 sigcme2=czero_gw; sigcme_3=czero_gw 494 495 ABI_MALLOC(sigctmp,(nomega_sigc,Sigp%nsig_ab)) 496 sigctmp=czero_gw 497 ABI_MALLOC(sigc_ket, (npwc*nspinor, nomega_sigc)) 498 499 ! Arrays storing the contribution given by the Hermitian/anti-Hermitian part of \Sigma_c 500 ABI_MALLOC(aherm_sigc_ket, (npwc*nspinor, nomega_sigc)) 501 ABI_MALLOC( herm_sigc_ket, (npwc*nspinor, nomega_sigc)) 502 503 sigcme_tmp=czero 504 505 ABI_MALLOC(sigc,(2,nomega_sigc,ib1:ib2,ib1:ib2,Wfd%nsppol*Sigp%nsig_ab)) 506 sigc=czero 507 508 !FIXME This quantities are only used for model GW if I am not wrong 509 ABI_MALLOC(ket1, (npwc*nspinor, nomega_tot)) 510 ABI_MALLOC(ket2, (npwc*nspinor, nomega_tot)) 511 ABI_MALLOC(omegame0i,(nomega_tot)) 512 513 ! Here we divide the states where the QP energies are required into complexes. Note however that this approach is not 514 ! based on group theory, and it might lead to spurious results in case of accidental degeneracies. 515 nq_summed=Kmesh%nbz 516 if (Sigp%symsigma>0) then 517 call littlegroup_print(Ltg_k,std_out,Dtset%prtvol,'COLL') 518 nq_summed=SUM(Ltg_k%ibzq(:)) 519 ! 520 ! Find number of degenerate subspaces and number of bands in each subspace 521 ! The tolerance is a little bit arbitrary (0.001 eV) 522 ! It could be reduced, in particular in case of nearly accidental degeneracies 523 ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,Wfd%nsppol)) 524 degtab=0 525 do spin=1,Wfd%nsppol 526 do ib=ib1,ib2 527 do jb=ib1,ib2 528 if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin))<0.001/Ha_ev) then 529 degtab(ib,jb,spin)=1 530 end if 531 end do 532 end do 533 end do 534 end if !symsigma 535 536 write(msg,'(2a,i6,a)')ch10,' calculation status ( ',nq_summed,' to be completed):' 537 call wrtout(std_out,msg,'COLL') 538 539 ! Here we have a problem in case of CD since epsm1q might be huge 540 ! TODO if single q (ex molecule) dont allocate epsm1q, avoid waste of memory 541 if ( ANY(mod10== [SIG_GW_AC, SIG_GW_CD, SIG_QPGW_CD])) then 542 if (.not.(mod10==SIG_GW_CD.and.Er%mqmem==0)) then 543 ABI_STAT_MALLOC(epsm1_qbz, (npwc, npwc, Er%nomega), ierr) 544 ABI_CHECK(ierr==0, "out-of-memory in epsm1_qbz") 545 end if 546 end if 547 548 ! TODO epsm1_trcc_qbz is needed for SIG_GW_CD with symmetries since 549 ! the Hermitian and the anti-Hermitian part have to be symmetrized in a different way. 550 !if (mod10==SIG_QPGW_CD) then 551 if (mod10==SIG_QPGW_CD) then 552 ABI_STAT_MALLOC(epsm1_trcc_qbz, (npwc, npwc, Er%nomega), ierr) 553 ABI_CHECK(ierr==0, "out-of-memory in epsm1_trcc_qbz") 554 ABI_MALLOC(epsm1_tmp, (npwc, npwc)) 555 end if 556 557 ABI_MALLOC(igfftcg0,(Gsph_Max%ng)) 558 ABI_MALLOC(ur_ibz,(gwc_nfftot*nspinor)) 559 ABI_MALLOC(usr_bz,(gwc_nfftot*nspinor)) 560 561 if (Dtset%pawcross==1) then 562 ABI_MALLOC(igfftfcg0,(Gsph_c%ng)) 563 ABI_MALLOC(ur_ae_sum,(nfftf*nspinor)) 564 ABI_MALLOC(ur_ae_onsite_sum,(nfftf*nspinor)) 565 ABI_MALLOC(ur_ps_onsite_sum,(nfftf*nspinor)) 566 end if 567 call timab(432,2,tsec) ! Init 568 569 ! ========================================== 570 ! ==== Fat loop over k_i in the full BZ ==== 571 ! ========================================== 572 do spin=1,Wfd%nsppol 573 if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE 574 call timab(433,1,tsec) ! Init spin 575 576 ! Load wavefunctions for GW corrections 577 ! TODO: Rotate the functions here instead of calling rho_tw_g 578 ABI_MALLOC(wfr_bdgw,(gwc_nfftot*nspinor,ib1:ib2)) 579 call wfd_get_many_ur(Wfd, [(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw) 580 581 if (Wfd%usepaw==1) then 582 ! Load cprj for GW states, note the indexing. 583 dimcprj_gw=nspinor*(ib2-ib1+1) 584 ABI_DT_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1)) 585 call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm) 586 ibsp=ib1 587 do jb=ib1,ib2 588 call wfd_get_cprj(Wfd,jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 589 call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 590 call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1))) 591 ibsp=ibsp+nspinor 592 end do 593 if (Dtset%pawcross==1) then 594 ABI_MALLOC(ur_ae_bdgw,(nfftf*nspinor,ib1:ib2)) 595 ABI_MALLOC(ur_ae_onsite_bdgw,(nfftf*nspinor,ib1:ib2)) 596 ABI_MALLOC(ur_ps_onsite_bdgw,(nfftf*nspinor,ib1:ib2)) 597 do jb=ib1,ib2 598 call wfd_paw_get_aeur(Wfdf,jb,jk_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,& 599 & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum) 600 ur_ae_bdgw(:,jb)=ur_ae_sum 601 ur_ae_onsite_bdgw(:,jb)=ur_ae_onsite_sum 602 ur_ps_onsite_bdgw(:,jb)=ur_ps_onsite_sum 603 end do 604 end if 605 end if ! usepaw 606 607 call timab(433,2,tsec) ! Init spin 608 609 do ik_bz=1,Kmesh%nbz 610 ! Parallelization over k-points and spin 611 ! For the spin there is another check in the inner loop 612 if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE 613 614 call timab(434,1,tsec) ! initq 615 616 ! Find the corresponding irreducible k-point 617 call get_BZ_item(Kmesh,ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt) 618 spinrot_kbz(:)=Cryst%spinrot(:,isym_ki) 619 620 ! Identify q and G0 where q + G0 = k_GW - k_i 621 kgw_m_ksum=kgw-ksum 622 call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0) 623 624 ! If symsigma, symmetrize the matrix elements. 625 ! Sum only q"s in IBZ_k. In this case elements are weighted 626 ! according to wtqp and wtqm. wtqm is for time-reversal. 627 wtqp=1; wtqm=0 628 if (can_symmetrize(spin)) then 629 if (Ltg_k%ibzq(iq_bz)/=1) CYCLE 630 wtqp=0; wtqm=0 631 do isym=1,Ltg_k%nsym_sg 632 wtqp=wtqp+Ltg_k%wtksym(1,isym,iq_bz) 633 wtqm=wtqm+Ltg_k%wtksym(2,isym,iq_bz) 634 end do 635 end if 636 637 write(msg,'(3(a,i0),a,i0)')' Sigma_c: ik_bz ',ik_bz,'/',Kmesh%nbz,", spin: ",spin,' done by mpi-rank: ',Wfd%my_rank 638 call wrtout(std_out,msg,'PERS') 639 640 ! Find the corresponding irred q-point. 641 call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q) 642 q_is_gamma = (normv(qbz,Cryst%gmet,"G") < GW_TOL_W0) 643 644 !q_is_gamma = (normv(qbz,Cryst%gmet,"G") < 0.7) 645 !if (iq_ibz/=2.and.iq_ibz/=1) CYCLE 646 !if (ANY(qbz<=-(half-tol16)) .or. ANY(qbz>(half+tol16))) CYCLE 647 !if (q_is_gamma) then; write(std_out,*)"skipping q=Gamma"; CYCLE; end if 648 ! 649 ! Tables for the FFT of the oscillators. 650 ! a) FFT index of the G-G0. 651 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 652 ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2)) 653 call gsph_fft_tabs(Gsph_c,g0,gwc_mgfft,gwc_ngfft,use_padfft,gw_gbound,igfftcg0) 654 655 if (ANY(gwc_fftalga == [2, 4])) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 656 #ifdef FC_IBM 657 ! XLF does not deserve this optimization (problem with [v67mbpt][t03]) 658 use_padfft = 0 659 #endif 660 if (use_padfft==0) then 661 ABI_FREE(gw_gbound) 662 ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2*use_padfft)) 663 end if 664 665 if (Dtset%pawcross==1) then 666 ABI_MALLOC(gboundf,(2*mgfftf+8,2)) 667 call gsph_fft_tabs(Gsph_c,g0,mgfftf,rho_ngfft,use_padfftf,gboundf,igfftfcg0) 668 if ( ANY(gwc_fftalga == (/2,4/)) ) use_padfftf=0 669 if (use_padfftf==0) then 670 ABI_FREE(gboundf) 671 ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf)) 672 end if 673 end if 674 675 ! Evaluate oscillator matrix elements 676 ! $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form. 677 if (Psps%usepaw==1) then 678 ABI_DT_MALLOC(Pwij_qg,(Psps%ntypat)) 679 q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT 680 call pawpwij_init(Pwij_qg,npwc,q0,Gsph_c%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 681 end if 682 683 if (Er%mqmem==0) then 684 ! Read q-slice of epsilon^{-1}|chi0 in Er%epsm1(:,:,:,1) (much slower but less memory). 685 call get_epsm1(Er,Vcp,0,0,Dtset%iomode,xmpi_comm_self,iqibzA=iq_ibz) 686 if (sigma_needs_ppm(Sigp)) then 687 if (Wfd%usepaw==1.and.PPm%userho==1) then 688 ! Use PAW AE rhor. 689 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,rho_nfftot,Gsph_c%gvec,& 690 & rho_ngfft,aepaw_rhor(:,1),iqiA=iq_ibz) 691 else 692 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,rho_nfftot,Gsph_c%gvec,& 693 & rho_ngfft,rhor(:,1),iqiA=iq_ibz) 694 end if 695 end if 696 end if 697 698 ! Symmetrize PPM parameters and epsm1 (q_IBZ --> q_BZ): 699 ! NOTE: 700 ! - We are not considering umklapp with G0/=0. In this case, 701 ! indeed the equation is different since we have to use G-G0. 702 ! A check, however, is performed in sigma. 703 ! - If gwcomp==1 and mod10 in [1,2,9], one needs both to set up botsq and epsm1_q 704 if (sigma_needs_ppm(Sigp)) then 705 ABI_MALLOC(botsq,(PPm%npwc,PPm%dm2_botsq)) 706 ABI_MALLOC(otq,(PPm%npwc,PPm%dm2_otq)) 707 ABI_MALLOC(eig,(PPm%dm_eig,PPm%dm_eig)) 708 call ppm_get_qbz(PPm,Gsph_c,Qmesh,iq_bz,botsq,otq,eig) 709 end if 710 711 if ( ANY(mod10 == [SIG_GW_AC,SIG_GW_CD,SIG_QPGW_CD])) then 712 ! Numerical integration or model GW with contour deformation or Analytic Continuation 713 ! TODO In case of AC we should symmetrize only the imaginary frequencies 714 if (mod10==SIG_GW_CD.and.Er%mqmem==0) then 715 ! Do in-place symmetrisation 716 call Epsm1_symmetrizer_inplace(iq_bz,Er%nomega, npwc,Er,Gsph_c,Qmesh,.TRUE.) 717 else 718 call Epsm1_symmetrizer(iq_bz,Er%nomega, npwc,Er,Gsph_c,Qmesh,.TRUE.,epsm1_qbz) 719 end if 720 721 if (mod10==SIG_GW_AC) then 722 ! Prepare the first term: \Sum w_i 1/z_i^2 f(1/z_i-1).. 723 ! The first frequencies are always real, skip them. 724 ! Memory is not optimized. 725 do iiw=1,Er%nomega_i 726 z2=gl_knots(iiw)*gl_knots(iiw) 727 ac_epsm1cqwz2(:,:,iiw)= gl_wts(iiw)*epsm1_qbz(:,:,Er%nomega_r+iiw)/z2 728 end do 729 end if 730 731 if (mod10==SIG_QPGW_CD) then 732 ! For model GW we need transpose(conjg(epsm1_qbz)) === 733 do io=1,Er%nomega 734 epsm1_tmp(:,:) = GWPC_CONJG(epsm1_qbz(:,:,io)) 735 epsm1_trcc_qbz(:,:,io) = TRANSPOSE(epsm1_tmp) 736 end do 737 end if 738 end if !gwcalctyp 739 740 ! Get Fourier components of the Coulombian interaction in the BZ === 741 ! In 3D systems, neglecting umklapp: vc(Sq,sG) = vc(q,G) = 4pi/|q+G|**2 742 ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S. 743 do ig=1,npwc 744 vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q))=Vcp%vc_sqrt(ig,iq_ibz) 745 end do 746 747 call timab(434,2,tsec) ! initq 748 749 ! Sum over band 750 call timab(445,1,tsec) ! loop 751 do ib=1,Sigp%nbnds 752 call timab(436,1,tsec) ! (1) 753 754 ! Parallelism over spin 755 ! This processor has this k-point but what about spin? 756 if (proc_distrb(ib,ik_bz,spin)/=Wfd%my_rank) CYCLE 757 758 call wfd_get_ur(Wfd,ib,ik_ibz,spin,ur_ibz) 759 760 if (Psps%usepaw==1) then 761 ! Load cprj for point ksum, this spin or spinor and *THIS* band. 762 ! TODO MG I could avoid doing this but I have to exchange spin and bands ??? 763 ! For sure there is a better way to do this! 764 call wfd_get_cprj(Wfd,ib,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 765 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 766 if (Dtset%pawcross==1) then 767 call wfd_paw_get_aeur(Wfdf,ib,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,& 768 & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum) 769 end if 770 end if 771 772 if (mod10==SIG_GW_AC) then 773 ! Calculate integral over omegap with Gauss-Legendre quadrature. 774 ! * -1/pi \int domegap epsm1c*(omega-e0i) / ( (omega-e0i)^2 + omegap^2) 775 ! * Note that energies are calculated wrt the Fermi level. 776 ac_integr(:,:,:)=czero_gw 777 do io=1,Sr%nomega_i 778 omegame0i_ac = Sr%omega_i(io)-qp_ene(ib,ik_ibz,spin) 779 omegame0i2_ac = omegame0i_ac*omegame0i_ac 780 do iiw=1,Er%nomega_i 781 do iggp=0,npwc*npwc-1 782 ig=iggp/npwc+1 783 igp= iggp-(ig-1)*npwc+1 ! \int domegap epsm1c/((omega-e0i)^2 + omegap^2) 784 ac_integr(ig,igp,io)= ac_integr(ig,igp,io) + ac_epsm1cqwz2(ig,igp,iiw)/(omegame0i2_ac + omegap2(iiw)) 785 end do 786 end do 787 ac_integr(:,:,io)=ac_integr(:,:,io)*omegame0i_ac 788 end do 789 ac_integr(:,:,:)=-ac_integr(:,:,:)*piinv 790 end if 791 792 call timab(436,2,tsec) ! (1) 793 call timab(437,1,tsec) ! rho_tw_g 794 795 ! Get all <k-q,ib,s|e^{-i(q+G).r}|s,jb,k>, at once 796 do jb=ib1,ib2 797 798 call rho_tw_g(nspinor,npwc,gwc_nfftot,ndat1,gwc_ngfft,1,use_padfft,igfftcg0,gw_gbound,& 799 & ur_ibz ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, & 800 & wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,& 801 & nspinor,rhotwg_ki(:,jb)) 802 803 if (Psps%usepaw==1) then 804 ! Add on-site contribution, projectors are already in BZ !TODO Recheck this! 805 i2=jb; if (nspinor==2) i2=(2*jb-1) 806 spad=(nspinor-1) 807 call paw_rho_tw_g(npwc,nspinor,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_c%gvec,& 808 & Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),Pwij_qg,rhotwg_ki(:,jb)) 809 if (Dtset%pawcross==1) then ! Add paw cross term 810 call paw_cross_rho_tw_g(nspinor,npwc,nfftf,rho_ngfft,1,use_padfftf,igfftfcg0,gboundf,& 811 & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum,iik,ktabrf(:,ik_bz),ph_mkt,spinrot_kbz,& 812 & ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,& 813 & nspinor,rhotwg_ki(:,jb)) 814 end if 815 end if 816 817 ! Multiply by the square root of the Coulomb term 818 ! In 3-D systems, the factor sqrt(4pi) is included) 819 do ii=1,nspinor 820 spad=(ii-1)*npwc 821 rhotwg_ki(spad+1:spad+npwc,jb) = rhotwg_ki(spad+1:spad+npwc,jb)*vc_sqrt_qbz(1:npwc) 822 end do 823 824 ! === Treat analytically the case q --> 0 === 825 ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma 826 ! while the Colulomb term is integrated out. 827 ! * In the scalar case we have nonzero contribution only if ib==jb 828 ! * For nspinor==2 evalute <ib,up|jb,up> and <ib,dwn|jb,dwn>, 829 ! impose orthonormalization since npwwfn might be < npwvec. 830 if (ik_bz==jk_bz) then 831 if (nspinor==1) then 832 rhotwg_ki(1,jb)=czero_gw 833 if (ib==jb) rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) 834 835 else 836 npw_k = Wfd%npwarr(ik_ibz) 837 rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero 838 if (ib==jb) then 839 cg_sum => Wfd%Wave(ib,ik_ibz,spin)%ug 840 cg_jb => Wfd%Wave(jb,jk_ibz,spin)%ug 841 ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1) 842 rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp) 843 ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1) 844 rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp) 845 ! PAW is missing 846 end if 847 end if 848 end if 849 end do !jb Got all matrix elements from minbnd up to maxbnd. 850 851 theta_mu_minus_e0i=fact_sp*qp_occ(ib,ik_ibz,spin) 852 853 ! Starting point to evaluate the derivative of Sigma and the Spectral function 854 e0i=qp_ene(ib,ik_ibz,spin) 855 856 ! Frequencies for the spectral function, e0i=qp_ene(ib,ik_ibz,spin) 857 ! FIXME the interval is not centered on eoi ! WHY? 858 if (Sr%nomega_r>0) omegame0i(1:Sr%nomega_r)=DBLE(Sr%omega_r(1:Sr%nomega_r))-e0i 859 860 call timab(437,2,tsec) ! rho_tw_g 861 862 do kb=ib1,ib2 863 call timab(438,1,tsec) ! (2) 864 865 ! Get frequencies $\omega$-\epsilon_in$ to evaluate $d\Sigma/dE$, note the spin 866 ! subtract e_KS since we have stored e_KS+ Delta \omega in Sr%omega4sd, not required for AC 867 do io=Sr%nomega_r+1,nomega_tot 868 omegame0i(io)=DBLE(Sr%omega4sd(kb,jk_ibz,io-Sr%nomega_r,spin))-e0i 869 end do 870 871 ! Get the ket \Sigma|\phi_{k,kb}> according to the method. 872 rhotwgp(:)=rhotwg_ki(:,kb) 873 sigc_ket = czero_gw 874 ket1 = czero_gw 875 ket2 = czero_gw 876 !aherm_sigc_ket = czero_gw 877 ! herm_sigc_ket = czero_gw 878 879 SELECT CASE (mod10) 880 CASE (SIG_GW_PPM) 881 ! GW WITH Plasmon-Pole Model. 882 ! Note that ppmodel 3 or 4 work only in case of standard perturbative approach! 883 ! Moreover, for ppmodel 3 and 4, spinorial case is not allowed 884 call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq,otq,& 885 & omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,sigc_ket,sigcme_3) 886 887 if (PPm%model==3.or.PPm%model==4) then 888 sigcme2(:,kb)=sigcme2(:,kb) + (wtqp+wtqm)*DBLE(sigcme_3(:)) + (wtqp-wtqm)*j_gw*AIMAG(sigcme_3(:)) 889 end if 890 891 CASE (SIG_GW_AC) 892 ! GW with Analytic continuation. 893 ! Evaluate \sum_Gp integr_GGp(omegasi) rhotw_Gp TODO this part can be optimized 894 do io=1,Sr%nomega_i 895 do ispinor=1,nspinor 896 spadc=(ispinor-1)*npwc 897 do ig=1,npwc 898 ctmp=czero 899 do igp=1,npwc 900 ctmp=ctmp+ac_integr(ig,igp,io)*rhotwgp(igp+spadc) 901 end do 902 sigc_ket(ig+spadc,io)=ctmp 903 end do 904 end do !ispinor 905 end do !io 906 907 CASE (SIG_GW_CD) 908 ! GW with contour deformation. 909 ! Check if pole contributions need to be summed 910 ! (this avoids unnecessary splint calls and saves time) 911 !me_calc_poles = .TRUE. 912 do io=1,nomega_tot 913 if (omegame0i(io)>=zero.AND.(ABS(one-theta_mu_minus_e0i)>zero)) then 914 !me_calc_poles(io) = .TRUE. 915 if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io)) 916 else if (omegame0i(io)<zero.AND.(ABS(theta_mu_minus_e0i)>zero)) then 917 !me_calc_poles(io) = .TRUE. 918 if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io)) 919 end if 920 end do 921 922 ! Check memory saving 923 if (Er%mqmem == 0) then 924 call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 925 & Er%omega,Er%epsm1(:,:,:,1),omegame0i,theta_mu_minus_e0i,sigc_ket,Dtset%ppmfrq,npoles_missing(kb),& 926 & method=Dtset%cd_frqim_method) 927 else 928 call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 929 & Er%omega,epsm1_qbz,omegame0i,theta_mu_minus_e0i,sigc_ket,Dtset%ppmfrq,npoles_missing(kb),& 930 & method=Dtset%cd_frqim_method) 931 end if 932 933 #if 0 934 if (wtqm/=0) then 935 call calc_sigc_cd(npwc,npwc,nspinor,,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 936 & Er%omega,epsm1_trcc_qbz,omegame0i,theta_mu_minus_e0i,aherm_sigc_ket,Dtset%ppmfrq,npoles_missing(kb),& 937 & method=Dtset%cd_frqim_method) 938 939 herm_sigc_ket = half*(sigc_ket + aherm_sigc_ket) 940 aherm_sigc_ket = half*(sigc_ket - aherm_sigc_ket) 941 else 942 herm_sigc_ket = sigc_ket 943 aherm_sigc_ket = czero_gw 944 end if 945 #endif 946 947 CASE (SIG_QPGW_PPM) 948 ! MODEL GW calculation WITH PPm TODO Spinor not tested. 949 ! Calculate \Sigma(E_k) |k> to obtain <j|\Sigma(E_k)|k> 950 ABI_MALLOC(sigcme_new,(nomega_tot)) 951 952 call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq,otq,& 953 & omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,ket1,sigcme_new) 954 955 if (Sigp%gwcalctyp==28) then 956 if (PPm%model/=1.and.PPm%model/=2) then 957 ! This is needed to have npwc=PPm%dm2_botsq=PPm%dm2_otq 958 write(msg,'(3a)')& 959 & 'For the time being, gwcalctyp=28 cannot be used with ppmodel=3,4.',ch10,& 960 & 'Use another Plasmon Pole Model when gwcalctyp=28.' 961 MSG_ERROR(msg) 962 end if 963 ABI_MALLOC(botsq_conjg_transp,(PPm%dm2_botsq,npwc)) 964 botsq_conjg_transp=TRANSPOSE(botsq) ! Keep these two lines separated, otherwise gfortran messes up 965 botsq_conjg_transp=CONJG(botsq_conjg_transp) 966 ABI_MALLOC(otq_transp,(PPm%dm2_otq,PPm%npwc)) 967 otq_transp=TRANSPOSE(otq) 968 969 call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq_conjg_transp,otq_transp,& 970 & omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,ket2,sigcme_3) 971 972 ABI_FREE(botsq_conjg_transp) 973 ABI_FREE(otq_transp) 974 sigc_ket= half*(ket1+ket2) 975 else 976 sigc_ket= ket1 977 end if 978 979 ABI_FREE(sigcme_new) 980 981 CASE (SIG_QPGW_CD) 982 ! MODEL GW with numerical integration. 983 ! Check if pole contributions need to be summed 984 ! (this avoids unnecessary splint calls and saves time) 985 !me_calc_poles = .TRUE. 986 do io=1,nomega_tot 987 if (omegame0i(io)>=zero.AND.(ABS(one-theta_mu_minus_e0i)>zero)) then 988 !me_calc_poles(io) = .TRUE. 989 if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io)) 990 else if (omegame0i(io)<zero.AND.(ABS(theta_mu_minus_e0i)>zero)) then 991 !me_calc_poles(io) = .TRUE. 992 if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io)) 993 end if 994 end do 995 996 ! Calculate \Sigma(E_k)|k> to obtain <j|\Sigma(E_k)|k> 997 call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 998 & Er%omega,epsm1_qbz,omegame0i,theta_mu_minus_e0i,ket1,Dtset%ppmfrq,npoles_missing(kb),& 999 & method=Dtset%cd_frqim_method) 1000 1001 if (Sigp%gwcalctyp==29) then 1002 ! Calculate \Sigma^*(E_k)|k> to obtain <k|\Sigma(E_k)|j>^* 1003 call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 1004 & Er%omega,epsm1_trcc_qbz,omegame0i,theta_mu_minus_e0i,ket2,Dtset%ppmfrq,npoles_missing(kb),& 1005 & method=Dtset%cd_frqim_method) 1006 sigc_ket = half*(ket1+ket2) 1007 else 1008 sigc_ket = ket1 1009 end if 1010 1011 CASE DEFAULT 1012 MSG_ERROR(sjoin("Unsupported value for mod10:", itoa(mod10))) 1013 END SELECT 1014 1015 if (Sigp%gwcomp==1) then 1016 ! TODO spinor not implemented 1017 call calc_sig_ppm_comp(npwc,nomega_tot,rhotwgp,botsq,otq,DBLE(Sr%egw(kb,jk_ibz,spin)-en_high),& 1018 & Sigp%zcut,theta_mu_minus_e0i,sigc_ket,PPm%model,npwc,PPm%dm2_botsq,PPm%dm2_otq) 1019 end if 1020 1021 call timab(438,2,tsec) ! 1022 call timab(439,1,tsec) ! sigma_me 1023 1024 ! Loop over the non-zero row elements of this column. 1025 ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS. 1026 ! 2) If gwcalctyp>=20: only off-diagonal elements connecting states with same character. 1027 do irow=1,Sigcij_tab(spin)%col(kb)%size1 1028 jb = Sigcij_tab(spin)%col(kb)%bidx(irow) 1029 rhotwg=rhotwg_ki(:,jb) 1030 1031 ! Calculate <\phi_j|\Sigma_c|\phi_k> 1032 ! Different freqs according to method (AC or Perturbative), see nomega_sigc. 1033 do iab=1,Sigp%nsig_ab 1034 spadc1 = spinor_padc(1, iab); spadc2 = spinor_padc(2, iab) 1035 do io=1,nomega_sigc 1036 sigctmp(io,iab) = XDOTC(npwc,rhotwg(spadc1+1:),1,sigc_ket(spadc2+1:,io),1) 1037 end do 1038 end do 1039 1040 if (Sigp%gwcomp==1) then 1041 ! Evaluate Extrapolar term TODO this does not work with spinor 1042 if (extrapolar_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank ) then 1043 ! Do it once as it does not depend on the ib index summed over. 1044 extrapolar_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank 1045 #if 1 1046 call calc_wfwfg(ktabr(:,jk_ibz),jik, spinrot_kgw, & ! TODO: why jk_ibz? 1047 gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g) 1048 #else 1049 call calc_wfwfg(ktabr(:,jk_bz),jik,spinrot_kgw,& 1050 gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g) 1051 #endif 1052 1053 if (Psps%usepaw==1) then 1054 i1=jb; i2=kb 1055 if (nspinor==2) then 1056 i1=(2*jb-1); i2=(2*kb-1) 1057 end if 1058 spad=(nspinor-1) 1059 call paw_rho_tw_g(gwc_nfftot,Sigp%nsig_ab,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,& 1060 & gw_gfft,Cprj_kgw(:,i1:i1+spad),Cprj_kgw(:,i2:i2+spad),Pwij_fft,wf1swf2_g) 1061 1062 if (Dtset%pawcross==1) then ! Add paw cross term 1063 call paw_cross_rho_tw_g(nspinor,npwc,nfftf,rho_ngfft,1,use_padfftf,igfftfcg0,gboundf,& 1064 & ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,& 1065 & ur_ae_bdgw(:,kb),ur_ae_onsite_bdgw(:,kb),ur_ps_onsite_bdgw(:,kb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,& 1066 & nspinor,wf1swf2_g) 1067 end if 1068 end if 1069 1070 ! The static contribution from completeness relation is calculated once === 1071 call calc_coh_comp(iq_ibz,Vcp%i_sz,(jb==kb),nspinor,Sigp%nsig_ab,DBLE(Sr%egw(kb,jk_ibz,spin)-en_high),& 1072 & npwc,Gsph_c%gvec,gwc_ngfft,gwc_nfftot,wf1swf2_g,vc_sqrt_qbz,botsq,otq,sigcohme) 1073 1074 do io=1,nomega_sigc 1075 sigctmp(io,:) = sigctmp(io,:)+sigcohme(:) 1076 end do 1077 end if ! gwcomp==1 1078 end if ! gwcom==1 1079 1080 ! Accumulate and, in case, symmetrize matrix elements of Sigma_c 1081 do iab=1,Sigp%nsig_ab 1082 is_idx=spin; if (nspinor==2) is_idx=iab 1083 1084 sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + & 1085 & (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab)) 1086 1087 sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp* sigctmp(:,iab) 1088 sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab)) 1089 ! TODO this should be the contribution coming from the anti-hermitian part. 1090 end do 1091 end do ! irow used to calculate matrix elements of $\Sigma$ 1092 1093 ! shaltaf (030406): this has to be done in a clean way later. 1094 ! TODO does not work with spinor. 1095 if (mod10==SIG_GW_PPM.and.(PPm%model==3.or.PPm%model==4)) then 1096 sigcme_tmp(:,kb,kb,spin)= sigcme2(:,kb) 1097 sigc(1,:,kb,kb,spin)= sigcme2(:,kb) 1098 sigc(2,:,kb,kb,spin)= czero 1099 end if 1100 1101 call timab(439,2,tsec) ! csigme(SigC) 1102 end do !kb to calculate matrix elements of $\Sigma$ 1103 end do !ib 1104 1105 call timab(445,2,tsec) ! csigme(SigC) 1106 1107 ! Deallocate k-dependent quantities. 1108 ABI_FREE(gw_gbound) 1109 if (Dtset%pawcross==1) then 1110 ABI_FREE(gboundf) 1111 end if 1112 1113 if (sigma_needs_ppm(Sigp)) then 1114 ABI_FREE(botsq) 1115 ABI_FREE(otq) 1116 ABI_FREE(eig) 1117 end if 1118 if (Psps%usepaw==1) then 1119 call pawpwij_free(Pwij_qg) 1120 ABI_DT_FREE(Pwij_qg) 1121 end if 1122 1123 end do ! ik_bz 1124 1125 ABI_FREE(wfr_bdgw) 1126 if (Wfd%usepaw==1) then 1127 call pawcprj_free(Cprj_kgw) 1128 ABI_DT_FREE(Cprj_kgw) 1129 if (Dtset%pawcross==1) then 1130 ABI_FREE(ur_ae_bdgw) 1131 ABI_FREE(ur_ae_onsite_bdgw) 1132 ABI_FREE(ur_ps_onsite_bdgw) 1133 end if 1134 end if 1135 end do ! spin 1136 1137 ABI_FREE(sigcme2) 1138 ABI_FREE(sigcme_3) 1139 ABI_FREE(igfftcg0) 1140 if (Dtset%pawcross==1) then 1141 ABI_FREE(igfftfcg0) 1142 end if 1143 1144 ! Gather contributions from all the CPUs 1145 call timab(440,1,tsec) ! wfd_barrier 1146 call timab(440,2,tsec) ! wfd_barrier 1147 call timab(441,1,tsec) ! xmpi_sum 1148 1149 call xmpi_sum(sigcme_tmp, wfd%comm, ierr) 1150 call xmpi_sum(sigc, wfd%comm, ierr) 1151 call timab(441,2,tsec) ! xmpi_sum 1152 1153 ! Multiply by constants. In 3D systems sqrt(4pi) is included in vc_sqrt_qbz === 1154 sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz) 1155 sigc = sigc /(Cryst%ucvol*Kmesh%nbz) 1156 1157 ! If we have summed over the IBZ_q now we have to average over complexes === 1158 ! Presently only diagonal terms are considered 1159 ! TODO QP-SCGW required a more involved approach, there is a check in sigma 1160 ! TODO it does not work if nspinor==2. 1161 call timab(442,1,tsec) ! final ops 1162 1163 do spin=1,Wfd%nsppol 1164 if (can_symmetrize(spin)) then 1165 if (mod10==SIG_GW_AC) then ! FIXME here there is a problem in case of AC with symmetries 1166 ABI_MALLOC(sym_cme, (Sr%nomega_i, ib1:ib2, ib1:ib2, sigp%nsig_ab)) 1167 else 1168 ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, sigp%nsig_ab)) 1169 end if 1170 sym_cme=czero 1171 1172 ! Average over degenerate diagonal elements 1173 ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results. 1174 ! another good reason to use a strict criterion for the tolerance on eigenvalues. 1175 do ib=ib1,ib2 1176 ndegs=0 1177 do jb=ib1,ib2 1178 if (degtab(ib,jb,spin)==1) then 1179 if (nspinor == 1) then 1180 sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:,:,jb,jb,spin), DIM=1) 1181 else 1182 do ii=1,sigp%nsig_ab 1183 sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:,:,jb,jb,ii), dim=1) 1184 end do 1185 end if 1186 end if 1187 ndegs = ndegs + degtab(ib,jb,spin) 1188 end do 1189 sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs 1190 end do 1191 1192 if (Sigp%gwcalctyp >= 20) then 1193 do iwc=1,nomega_sigc 1194 call esymm_symmetrize_mels(QP_sym(spin),ib1,ib2,sigc(:,iwc,:,:,spin),sym_cme(iwc,:,:,1)) 1195 end do 1196 end if 1197 1198 ! Copy symmetrized values 1199 do ib=ib1,ib2 1200 do jb=ib1,ib2 1201 !if (mod10==SIG_GW_AC.and.average_real) CYCLE ! this is to check another scheme in case of AC 1202 if (nspinor == 1) then 1203 sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1) 1204 else 1205 do ii=1,sigp%nsig_ab 1206 sigcme_tmp(:,ib,jb,ii) = sym_cme(:,ib,jb,ii) 1207 end do 1208 end if 1209 end do 1210 end do 1211 ABI_FREE(sym_cme) 1212 end if 1213 end do 1214 1215 ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX) 1216 !if (Sigp%gwcalctyp>=20 .and. sigma_is_herm(Sigp) ) then 1217 ! ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!") 1218 ! do spin=1,Wfd%nsppol 1219 ! do io=1,nomega_sigc 1220 ! call hermitianize(sigcme_tmp(io,:,:,spin),"Upper") 1221 ! end do 1222 ! end do 1223 !end if 1224 1225 ! GW with contour deformation: check on the number of poles not included. 1226 if (ANY(mod10 == [SIG_GW_CD,SIG_QPGW_CD])) then 1227 call xmpi_sum(npoles_missing, wfd%comm, ierr) 1228 npls = SUM(npoles_missing) 1229 if (npls>0) then 1230 MSG_WARNING(sjoin("Total number of missing poles for contour deformation method:", itoa(npls))) 1231 do band=minbnd,maxbnd 1232 npls = npoles_missing(band) 1233 if (npls>0) then 1234 write(msg,'(a,2(i0,a))')" For band ",band," there are ",npls," missing poles" 1235 call wrtout(std_out,msg,"COLL") 1236 end if 1237 end do 1238 end if 1239 ! Print data on the maximum value needed for the screening along the real axis 1240 w_localmax = MAXVAL(w_maxval) 1241 call xmpi_max(w_localmax,w_max, wfd%comm, ierr) 1242 write(msg,'(a,f12.5,a)') ' Max omega value used in W(omega): ',w_max*Ha_eV,' [eV]' 1243 call wrtout(std_out,msg,"COLL") 1244 end if 1245 call timab(442,2,tsec) ! final ops 1246 1247 ! =========================== 1248 ! ==== Deallocate memory ==== 1249 ! =========================== 1250 if (Psps%usepaw==1) then 1251 if (allocated(gw_gfft)) then 1252 ABI_FREE(gw_gfft) 1253 end if 1254 call pawcprj_free(Cprj_ksum) 1255 ABI_DT_FREE(Cprj_ksum) 1256 if (allocated(Pwij_fft)) then 1257 call pawpwij_free(Pwij_fft) 1258 ABI_DT_FREE(Pwij_fft) 1259 end if 1260 if (Dtset%pawcross==1) then 1261 ABI_FREE(ur_ae_sum) 1262 ABI_FREE(ur_ae_onsite_sum) 1263 ABI_FREE(ur_ps_onsite_sum) 1264 ABI_FREE(ktabrf) 1265 end if 1266 end if 1267 1268 ABI_FREE(npoles_missing) 1269 ABI_FREE(ur_ibz) 1270 ABI_FREE(usr_bz) 1271 ABI_FREE(ktabr) 1272 ABI_FREE(sigc_ket) 1273 ABI_FREE(ket1) 1274 ABI_FREE(ket2) 1275 ABI_FREE(rhotwg_ki) 1276 ABI_FREE(rhotwg) 1277 ABI_FREE(rhotwgp) 1278 ABI_FREE(vc_sqrt_qbz) 1279 ABI_FREE(omegame0i) 1280 ABI_FREE(sigctmp) 1281 ABI_FREE(sigc) 1282 ABI_FREE(w_maxval) 1283 1284 if (allocated(epsm1_qbz)) then 1285 ABI_FREE(epsm1_qbz) 1286 end if 1287 if (allocated(epsm1_trcc_qbz)) then 1288 ABI_FREE(epsm1_trcc_qbz) 1289 end if 1290 if (allocated(epsm1_tmp)) then 1291 ABI_FREE(epsm1_tmp) 1292 end if 1293 if (allocated(degtab)) then 1294 ABI_FREE(degtab) 1295 end if 1296 if (allocated(ac_epsm1cqwz2)) then 1297 ABI_FREE(ac_epsm1cqwz2) 1298 end if 1299 if (allocated(ac_integr)) then 1300 ABI_FREE(ac_integr) 1301 end if 1302 if (allocated(aherm_sigc_ket)) then 1303 ABI_FREE(aherm_sigc_ket) 1304 end if 1305 if (allocated(herm_sigc_ket)) then 1306 ABI_FREE(herm_sigc_ket) 1307 end if 1308 if (Sigp%gwcomp==1) then 1309 ABI_FREE(wf1swf2_g) 1310 endif 1311 if (Sigp%gwcomp==1) then 1312 ABI_FREE(extrapolar_distrb) 1313 end if 1314 ABI_FREE(proc_distrb) 1315 1316 call timab(431,2,tsec) 1317 call timab(424,2,tsec) ! calc_sigc_me 1318 1319 call cwtime(cpu_time,wall_time,gflops,"stop") 1320 write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time 1321 1322 DBG_EXIT("COLL") 1323 1324 end subroutine calc_sigc_me