TABLE OF CONTENTS
ABINIT/calc_coh [ Functions ]
NAME
calc_coh
FUNCTION
Calculates the partial contribution to the COH part of the COHSEX self-energy for a given q-point.
INPUTS
iqibz=index of the irreducible q-point in the array qibz, point which is related by a symmetry operation to the point q summed over (see csigme). This index is also used to treat the integrable coulombian singularity at q=0 ngfft(18)=contain all needed information about 3D FFT for GW wavefuntions, see ~abinit/doc/variables/vargs.htm#ngfft nsig_ab=Number of components in the self-energy operator (1 for collinear magnetism) npwc=number of plane waves in $\tilde epsilon^{-1}$ nspinor=Number of spinorial components. nfftot=number of points in real space i_sz=contribution arising from the integrable coulombian singularity at q==0 (see csigme for the method used), note that in case of 3-D systems the factor 4pi in the coulombian potential is included in the definition of i_sz gvec(3,npwc)=G vectors in reduced coordinates vc_sqrt(npwc)= square root of the coulombian matrix elements for this q-point epsm1q_o(npwc,npwc)= contains $\tilde epsilon^{-1}(q,w=0) - \delta_{G Gp}$ for the particular q-point considered in the sum wfg2_jk(nsig_ab*nfftot)= Fourier Transform of $\u_{jb k}^*(r) u_{kb k}$ jb,kb=left and righ band indices definining the left and right states where the partial contribution to the matrix element of $\Sigma_{COH}$ is evaluated
OUTPUT
sigcohme=partial contribution to the matrix element of $<jb k \sigma|\Sigma_{COH} | kb k \sigma>$ coming from this single q-point
SOURCE
871 subroutine calc_coh(nspinor,nsig_ab,nfftot,ngfft,npwc,gvec,wfg2_jk,epsm1q_o,vc_sqrt,i_sz,iqibz,same_band,sigcohme) 872 873 !Arguments ------------------------------------ 874 !scalars 875 integer,intent(in) :: iqibz,nfftot,npwc,nsig_ab,nspinor 876 real(dp),intent(in) :: i_sz 877 logical,intent(in) :: same_band 878 !arrays 879 integer,intent(in) :: gvec(3,npwc),ngfft(18) 880 complex(gwpc),intent(in) :: epsm1q_o(npwc,npwc),vc_sqrt(npwc) 881 complex(gwpc),intent(in) :: wfg2_jk(nfftot*nsig_ab) 882 complex(gwpc),intent(out) :: sigcohme(nsig_ab) 883 884 !Local variables------------------------------- 885 !scalars 886 integer,save :: enough=0 887 integer :: ig,ig4,ig4x,ig4y,ig4z,igp,igmin,ispinor,spad,outofbox 888 !arrays 889 integer :: g2mg1(3) 890 891 ! ************************************************************************* 892 893 DBG_ENTER("COLL") 894 895 ! === Partial contribution to the matrix element of Sigma_c === 896 ! * For nspinor==2, the closure relation reads: 897 ! $\sum_s \psi_a^*(1)\psi_b(2) = \delta_{ab} \delta(1-2)$ 898 ! where a,b are the spinor components. As a consequence, Sigma_{COH} is always 899 ! diagonal in spin-space and only diagonal matrix elements have to be calculated. 900 ! MG TODO wfg2_jk should be calculated on an augmented FFT box to avoid spurious wrapping of G1-G2. 901 ! MG: One has to make sure G1-G2 is still in the FFT mesh for each G1 and G2 in chi0 (not always true) 902 ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic! 903 904 ! * Treat the case q --> 0 adequately. 905 ! TODO Better treatment of wings, check cutoff in the coulombian interaction. 906 igmin=1; if (iqibz==1) igmin=2 907 908 sigcohme(:)=czero_gw 909 910 do ispinor=1,nspinor 911 spad=(ispinor-1)*nfftot 912 outofbox=0 913 914 do igp=igmin,npwc 915 do ig=igmin,npwc 916 917 g2mg1 = gvec(:,igp)-gvec(:,ig) 918 if (ANY(g2mg1(:)>ngfft(1:3)/2) .or. ANY(g2mg1(:)<-(ngfft(1:3)-1)/2)) then 919 outofbox = outofbox+1; CYCLE 920 end if 921 922 ig4x=MODULO(g2mg1(1),ngfft(1)) 923 ig4y=MODULO(g2mg1(2),ngfft(2)) 924 ig4z=MODULO(g2mg1(3),ngfft(3)) 925 ig4= 1+ig4x+ig4y*ngfft(1)+ig4z*ngfft(1)*ngfft(2) 926 927 sigcohme(ispinor) = sigcohme(ispinor) + & 928 & half*wfg2_jk(spad+ig4)*epsm1q_o(ig,igp)*vc_sqrt(ig)*vc_sqrt(igp) 929 end do !ig 930 end do !igp 931 932 if (iqibz==1.and.same_band) then 933 sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*epsm1q_o(1,1)*i_sz 934 end if 935 end do !ispinor 936 937 if (outofbox/=0) then 938 enough=enough+1 939 if (enough<=50) then 940 ABI_WARNING(sjoin(' Number of G1-G2 pairs outside the G-sphere for Wfns:', itoa(outofbox))) 941 if (enough==50) then 942 call wrtout(std_out,' ========== Stop writing Warnings ==========') 943 end if 944 end if 945 end if 946 947 DBG_EXIT("COLL") 948 949 end subroutine calc_coh
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-2024 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) [[cite:Gigy1986]] 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.
SOURCE
150 subroutine cohsex_me(sigmak_ibz,ikcalc,nomega_sigc,minbnd,maxbnd,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_c,Vcp,& 151 & Kmesh,Qmesh,Ltg_k,Pawtab,Pawang,Paw_pwff,Psps,Wfd,allQP_sym,gwc_ngfft,iomode,prtvol,sigcme_tmp) 152 153 !Arguments ------------------------------------ 154 !scalars 155 integer,intent(in) :: sigmak_ibz,ikcalc,prtvol,iomode,nomega_sigc,minbnd,maxbnd 156 type(crystal_t),intent(in) :: Cryst 157 type(ebands_t),target,intent(in) :: QP_BSt 158 type(kmesh_t),intent(in) :: Kmesh,Qmesh 159 type(vcoul_t),intent(in) :: Vcp 160 type(Epsilonm1_results),intent(inout) :: Er 161 type(gsphere_t),intent(in) :: Gsph_c 162 type(littlegroup_t),intent(in) :: Ltg_k 163 type(Pseudopotential_type),intent(in) :: Psps 164 type(pawang_type),intent(in) :: pawang 165 type(sigparams_t),target,intent(in) :: Sigp 166 type(sigma_t),intent(in) :: Sr 167 type(wfdgw_t),target,intent(inout) :: Wfd 168 !arrays 169 integer,intent(in) :: gwc_ngfft(18) 170 complex(dpc),intent(out) :: sigcme_tmp(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Wfd%nsppol*Sigp%nsig_ab) 171 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat) 172 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 173 type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol) 174 175 !Local variables ------------------------------ 176 !scalars 177 integer,parameter :: tim_fourdp=2,ndat1=1 178 integer :: iab,ib,ib1,ib2,ierr,ig,ii,iik,itim_q,i1,i2,npwc 179 integer :: ik_bz,ik_ibz,io,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx 180 integer :: band,band1,band2,idle,rank 181 integer :: jik,jk_bz,jk_ibz,kb,nspinor,nsppol 182 integer :: nomega_tot,nq_summed,ispinor,ibsp,dimcprj_gw 183 integer :: spad,spadc,spadc1,spadc2,irow,my_nbks 184 integer :: ndegs,wtqm,wtqp,mod10 185 integer :: isym_kgw,isym_ki,gwc_mgfft,use_padfft,gwc_fftalga,gwc_nfftot,ifft,npw_k 186 real(dp) :: fact_sp,theta_mu_minus_e0i,tol_empty,gw_gsq 187 complex(dpc) :: ctmp,ph_mkgwt,ph_mkt 188 logical :: iscompatibleFFT,q_is_gamma 189 character(len=500) :: msg 190 type(wave_t),pointer :: wave_sum, wave_jb 191 !arrays 192 integer :: g0(3),spinor_padc(2,4),nbv_ks(Kmesh%nibz,Wfd%nsppol) 193 integer,allocatable :: proc_distrb(:,:,:),coh_distrb(:,:,:,:),degtab(:,:,:) 194 integer,allocatable :: igfftcg0(:),gw_gfft(:,:),gw_gbound(:,:),irottb(:,:),ktabr(:,:) 195 integer :: got(Wfd%nproc) 196 real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),q0(3),tsec(2),qbz(3),spinrot_kbz(4),spinrot_kgw(4) 197 real(dp),ABI_CONTIGUOUS pointer :: qp_ene(:,:,:),qp_occ(:,:,:) 198 complex(gwpc) :: sigcohme(Sigp%nsig_ab) 199 complex(gwpc),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:),sigsex(:) 200 complex(gwpc),allocatable :: epsm1_qbz(:,:,:) 201 complex(gwpc),allocatable :: sigc_ket(:,:) 202 complex(gwpc),allocatable :: rhotwg_ki(:,:) 203 complex(gwpc),allocatable :: sigctmp(:,:) 204 complex(gwpc),allocatable :: wfr_bdgw(:,:),ur_sum(:),wf1swf2_g(:) 205 complex(gwpc),ABI_CONTIGUOUS pointer :: cg_jb(:),cg_sum(:) 206 complex(dpc),allocatable :: sym_cme(:,:,:,:),sigc(:,:,:,:,:) 207 logical :: rank_mask(Wfd%nproc),can_symmetrize(Wfd%nsppol) 208 logical,allocatable :: bks_mask(:,:,:) 209 type(sigijtab_t),pointer :: Sigcij_tab(:) 210 type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:) 211 type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:) 212 type(esymm_t),pointer :: QP_sym(:) 213 214 !************************************************************************ 215 216 DBG_ENTER("COLL") 217 218 call timab(423,1,tsec) ! cohsex_me 219 220 ! Initial check 221 ABI_CHECK(Sr%nomega_r == Sigp%nomegasr,"") 222 ABI_CHECK(Sr%nomega4sd == Sigp%nomegasrd,"") 223 !ABI_CHECK(Sigp%npwc==Gsph_c%ng,"") 224 225 ! Initialize some values 226 nspinor = Wfd%nspinor; nsppol = Wfd%nsppol 227 npwc = sigp%npwc 228 spinor_padc = RESHAPE([0, 0, npwc, npwc, 0, npwc, npwc,0], [2, 4]) 229 230 qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ 231 232 ! Extract the symmetries of the bands for this k-point 233 QP_sym => allQP_sym(sigmak_ibz,1:nsppol) 234 235 ! Index of the GW point in the BZ array, its image in IBZ and time-reversal === 236 jk_bz=Sigp%kptgw2bz(ikcalc) 237 call kmesh%get_BZ_item(jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt) 238 !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk) 239 spinrot_kgw=Cryst%spinrot(:,isym_kgw) 240 ib1 = minbnd; ib2 = maxbnd 241 242 write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,& 243 ' Calculating <nk|Sigma_c(omega)|nk> at k = ',kgw(:),ch10,& 244 ' bands n = from ',ib1,' to ',ib2,ch10 245 call wrtout(std_out,msg) 246 247 if (ANY(gwc_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst,Psps,gwc_ngfft) 248 gwc_mgfft = MAXVAL(gwc_ngfft(1:3)) 249 gwc_fftalga = gwc_ngfft(7)/100 !; gwc_fftalgc=MOD(gwc_ngfft(7),10) 250 251 can_symmetrize = .FALSE. 252 if (Sigp%symsigma>0) then 253 can_symmetrize = .TRUE. 254 if (Sigp%gwcalctyp >= 20) then 255 do spin=1,Wfd%nsppol 256 can_symmetrize(spin) = .not. esymm_failed(QP_sym(spin)) 257 if (.not.can_symmetrize(spin)) then 258 write(msg,'(a,i0,4a)')" Symmetrization cannot be performed for spin: ",spin,ch10,& 259 & " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg) 260 ABI_WARNING(msg) 261 end if 262 end do 263 end if 264 if (nspinor == 2) ABI_WARNING('Symmetrization with nspinor=2 not implemented') 265 end if 266 267 mod10=MOD(Sigp%gwcalctyp, 10) 268 269 call timab(491,1,tsec) ! csigme(tot) Overall clock. TODO check this 270 call timab(495,1,tsec) ! csigme (SigC) 271 272 ! Normalization of theta_mu_minus_e0i 273 ! If nsppol==2, qp_occ $\in [0,1]$ 274 SELECT CASE (nsppol) 275 CASE (1) 276 fact_sp=half; tol_empty=0.01 ! below this value the state is assumed empty 277 if (nspinor==2) then 278 fact_sp=one; tol_empty=0.005 ! below this value the state is assumed empty 279 end if 280 CASE (2) 281 fact_sp=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic 282 CASE DEFAULT ! spin unpolarized system is treated using nsppol==2 283 ABI_BUG('Wrong nsppol') 284 END SELECT 285 286 call timab(442,1,tsec) ! csigme(init0) 287 288 ! Precalculate the FFT index of $(R^{-1}(r-\tau))$ === 289 ! S=\transpose R^{-1} and k_BZ = S k_IBZ 290 ! irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk. 291 gwc_nfftot = PRODUCT(gwc_ngfft(1:3)) 292 ABI_MALLOC(irottb,(gwc_nfftot,Cryst%nsym)) 293 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwc_ngfft,irottb,iscompatibleFFT) 294 if (.not.iscompatibleFFT) then 295 ABI_WARNING("FFT mesh is not compatible with symmetries. Results might be affected by large errors!") 296 end if 297 298 ABI_MALLOC(ktabr,(gwc_nfftot, Kmesh%nbz)) 299 do ik_bz=1,Kmesh%nbz 300 isym=Kmesh%tabo(ik_bz) 301 do ifft=1,gwc_nfftot 302 ktabr(ifft,ik_bz)=irottb(ifft,isym) 303 end do 304 end do 305 ABI_FREE(irottb) 306 307 ! The number of occupied states for each point in the IBZ and spin. 308 do spin=1,nsppol 309 do ik_ibz=1,Kmesh%nibz 310 nbv_ks(ik_ibz,spin) = COUNT(qp_occ(:,ik_ibz,spin)>=tol_empty) 311 end do 312 end do 313 314 ! (b,k,s) mask for MPI distribution of the sum over occupied states in the BZ. 315 ABI_MALLOC(bks_mask,(Wfd%mband,Kmesh%nbz,nsppol)) 316 bks_mask=.FALSE. 317 do spin=1,nsppol 318 do ik_bz=1,Kmesh%nbz 319 ik_ibz = Kmesh%tab(ik_bz) 320 bks_mask(1:nbv_ks(ik_ibz,spin),ik_bz,spin) = .TRUE. 321 end do 322 end do 323 324 ! Distribute the individual terms of the sum over the BZ taking into account symmetries and MPI memory distribution. 325 ! got is used to optimize the distribution if more than one node can calculate the same (b,k,s) element. 326 got=0 327 ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,nsppol)) 328 call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,& 329 & proc_distrb,got,bks_mask,global=.TRUE.) 330 331 ABI_FREE(bks_mask) 332 333 write(msg,'(a,i0,a)')" Will sum ",my_nbks," (b,k,s) occupied states in (COHSEX|SEX)." 334 call wrtout(std_out,msg,'PERS') 335 336 Sigcij_tab => Sigp%Sigcij_tab(ikcalc,1:nsppol) 337 338 if (mod10==SIG_COHSEX) then 339 ! Distribute the COHSEX terms, taking into account the symmetries of the Sigma_ij matrix. 340 ABI_MALLOC(coh_distrb,(ib1:ib2,ib1:ib2,Kmesh%nbz,nsppol)) 341 342 coh_distrb = xmpi_undefined_rank 343 do spin=1,nsppol 344 do ik_bz=1,Kmesh%nbz 345 if (ANY(proc_distrb(:,ik_bz,spin) /= xmpi_undefined_rank) ) then ! This BZ point will be calculated. 346 rank_mask = .FALSE. ! To select only those nodes that will treat (k,s). 347 do band=1,Wfd%mband 348 rank = proc_distrb(band,ik_bz,spin) 349 if (rank /= xmpi_undefined_rank) rank_mask(rank+1)=.TRUE. 350 end do 351 do band2=ib1,ib2 352 do irow=1,Sigcij_tab(spin)%col(band2)%size1 ! Looping over the upper triangle of sigma_ij with non-zero elements. 353 band1 = Sigcij_tab(spin)%col(band2)%bidx(irow) 354 idle = imin_loc(got,mask=rank_mask) 355 got(idle) = got(idle)+1 356 coh_distrb(band1,band2,ik_bz,spin) = idle-1 357 end do 358 end do 359 end if 360 end do 361 end do 362 363 write(msg,'(a,i0,a)')" will treat ",COUNT(coh_distrb==Wfd%my_rank)," COH terms." 364 call wrtout(std_out,msg,'PERS') 365 end if 366 367 ABI_MALLOC(rhotwg_ki, (npwc * nspinor, minbnd:maxbnd)) 368 rhotwg_ki=czero_gw 369 ABI_MALLOC(rhotwg, (npwc * nspinor)) 370 ABI_MALLOC(rhotwgp ,(npwc * nspinor)) 371 ABI_MALLOC(vc_sqrt_qbz, (npwc)) 372 373 ! Additional allocations for PAW 374 if (Psps%usepaw==1) then 375 ABI_MALLOC(Cprj_ksum,(Cryst%natom,nspinor)) 376 call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm) 377 378 ! For COHSEX we need the onsite terms of the PW on the FFT mesh. 379 ! gw_gfft is the set of plane waves in the FFT Box for the oscillators. 380 if (mod10==SIG_COHSEX) then 381 ABI_MALLOC(gw_gfft,(3,gwc_nfftot)) 382 q0=zero 383 call get_gftt(gwc_ngfft,q0,Cryst%gmet,gw_gsq,gw_gfft) 384 ABI_MALLOC(Pwij_fft,(Psps%ntypat)) 385 call pawpwij_init(Pwij_fft,gwc_nfftot,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 386 end if 387 end if ! usepaw==1 388 389 ! === Calculate total number of frequencies and allocate related arrays === 390 ! sigcme2 is used to accumulate the diagonal matrix elements over k-points and 391 ! GW bands, used only in case of ppmodel 3 and 4 (TODO save memory) 392 nomega_tot=Sr%nomega_r+Sr%nomega4sd 393 394 ABI_MALLOC(sigctmp, (nomega_sigc,Sigp%nsig_ab)) 395 sigctmp = czero_gw 396 ABI_MALLOC(sigc_ket, (npwc*nspinor, nomega_sigc)) 397 398 if (mod10==SIG_COHSEX) then 399 ABI_MALLOC(wf1swf2_g,(gwc_nfftot*nspinor)) 400 end if 401 402 ! Arrays storing the contribution given by the Hermitian/anti-Hermitian part of \Sigma_c 403 !allocate(aherm_sigc_ket(npwc*nspinor,nomega_sigc)) 404 !allocate( herm_sigc_ket(npwc*nspinor,nomega_sigc)) 405 ABI_MALLOC(sigsex,(npwc)) 406 sigcme_tmp=czero 407 408 ABI_MALLOC(sigc,(2,nomega_sigc,ib1:ib2,ib1:ib2,nsppol*Sigp%nsig_ab)) 409 sigc=czero 410 411 ! Here we divide the states where the QP energies are required into complexes. Note however that this approach is not 412 ! based on group theory, and it might lead to spurious results in case of accidental degeneracies. 413 nq_summed=Kmesh%nbz 414 if (Sigp%symsigma > 0) then 415 call Ltg_k%print(std_out, prtvol, mode_paral='COLL') 416 nq_summed=SUM(Ltg_k%ibzq(:)) 417 ! 418 ! Find number of degenerate states and number of bands in each subspace 419 ! The tolerance is a little bit arbitrary (0.001 eV) 420 ! It could be reduced, in particular in case of nearly accidental degeneracies 421 ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,nsppol)) 422 degtab=0 423 do spin=1,nsppol 424 do ib=ib1,ib2 425 do jb=ib1,ib2 426 if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin))<0.001/Ha_ev) then 427 degtab(ib,jb,spin)=1 428 end if 429 end do 430 end do 431 end do 432 end if !symsigma 433 434 write(msg,'(2a,i6,a)')ch10,' calculation status ( ',nq_summed,' to be completed):' 435 call wrtout(std_out,msg) 436 437 ! TODO if single q (ex molecule) dont allocate epsm1q, avoid waste of memory 438 ABI_MALLOC_OR_DIE(epsm1_qbz, (npwc, npwc, 1), ierr) 439 ABI_MALLOC(igfftcg0,(Gsph_c%ng)) 440 441 ! Out-of-core solution for epsilon. 442 if (Er%mqmem==0) then 443 ABI_COMMENT('Reading q-slices from file. Slower but less memory.') 444 end if 445 446 call timab(442,2,tsec) 447 448 ! ========================================== 449 ! ==== Fat loop over k_i in the full BZ ==== 450 ! ========================================== 451 ABI_MALLOC(ur_sum,(gwc_nfftot*nspinor)) 452 453 do spin=1,nsppol 454 if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE 455 456 ABI_MALLOC(wfr_bdgw,(gwc_nfftot*nspinor,ib1:ib2)) 457 call wfd%get_many_ur([(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw) 458 459 if (Wfd%usepaw==1) then 460 ! Load cprj for GW states, note the indexing. 461 dimcprj_gw=nspinor*(ib2-ib1+1) 462 ABI_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1)) 463 call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm) 464 ibsp=ib1 465 do jb=ib1,ib2 466 call wfd%get_cprj(jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 467 call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 468 call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1))) 469 ibsp=ibsp+nspinor 470 end do 471 end if 472 473 do ik_bz=1,Kmesh%nbz 474 ! Parallelization over k-points and spin 475 ! For the spin there is another check in the inner loop. 476 if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE 477 478 call timab(443,1,tsec) ! csigme (initq) 479 480 ! Find the corresponding irreducible k-point 481 call kmesh%get_BZ_item(ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt) 482 spinrot_kbz(:)=Cryst%spinrot(:,isym_ki) 483 484 ! Identify q and G0 where q+G0=k_GW-k_i 485 kgw_m_ksum=kgw-ksum 486 call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0) 487 488 ! Symmetrize the matrix elements. 489 ! Sum only q"s in IBZ_k. In this case elements are weighted 490 ! according to wtqp and wtqm. wtqm is for time-reversal. 491 wtqp=1; wtqm=0 492 if (can_symmetrize(spin)) then 493 if (Ltg_k%ibzq(iq_bz)/=1) CYCLE 494 wtqp=0; wtqm=0 495 do isym=1,Ltg_k%nsym_sg 496 wtqp=wtqp+Ltg_k%wtksym(1,isym,iq_bz) 497 wtqm=wtqm+Ltg_k%wtksym(2,isym,iq_bz) 498 end do 499 end if 500 501 !%write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank 502 !%call wrtout(std_out,msg,'PERS') 503 504 ! Find the corresponding irred q-point. 505 call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q) 506 q_is_gamma = (normv(qbz, Cryst%gmet, "G") < GW_TOLQ0) 507 508 ! Tables for the FFT of the oscillators. 509 ! a) FFT index of the G-G0. 510 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 511 ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2)) 512 call Gsph_c%fft_tabs(g0,gwc_mgfft,gwc_ngfft,use_padfft,gw_gbound,igfftcg0) 513 if ( ANY(gwc_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 514 #ifdef FC_IBM 515 ! XLF does not deserve this optimization (problem with [v67mbpt][t03]) 516 use_padfft = 0 517 #endif 518 if (use_padfft==0) then 519 ABI_FREE(gw_gbound) 520 ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2*use_padfft)) 521 end if 522 523 if (Psps%usepaw==1) then 524 ! Get PAW oscillator matrix elements $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form. 525 ABI_MALLOC(Pwij_qg,(Psps%ntypat)) 526 q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT 527 call pawpwij_init(Pwij_qg, npwc, q0, Gsph_c%gvec, Cryst%rprimd, Psps, Pawtab, Paw_pwff) 528 end if 529 530 if (Er%mqmem==0) then 531 ! Read q-slice of epsilon^{-1}|chi0 in Er%epsm1(:,:,:,1) (much slower but less memory). 532 call get_epsm1(Er,Vcp,0,0,iomode,xmpi_comm_self,iqibzA=iq_ibz) 533 end if 534 535 ! Only omega==0 for SEX or COHSEX 536 call Epsm1_symmetrizer(iq_bz, 1, npwc, Er, Gsph_c, Qmesh, .True., epsm1_qbz) 537 538 ! Get Fourier components of the Coulomb interaction in the BZ. 539 ! In 3D systems, neglecting umklapp, vc(Sq,sG)=vc(q,G)=4pi/|q+G| 540 ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S. 541 do ig=1,npwc 542 vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q)) = Vcp%vc_sqrt(ig,iq_ibz) 543 end do 544 545 call timab(443,2,tsec) ! csigme (initq) 546 547 ! Sum over bands. 548 do ib=1,Sigp%nbnds 549 ! Parallelism over spin 550 ! This processor has this k-point but what about spin? 551 if (proc_distrb(ib,ik_bz,spin)/=Wfd%my_rank) CYCLE 552 553 ! Skip empty state ib for HF, SEX, and COHSEX. 554 if (qp_occ(ib,ik_ibz,spin)<tol_empty) CYCLE 555 556 theta_mu_minus_e0i=fact_sp*qp_occ(ib,ik_ibz,spin) 557 558 call wfd%get_ur(ib,ik_ibz,spin,ur_sum) 559 560 if (Psps%usepaw==1) then 561 ! Load cprj for point ksum, this spin or spinor and *THIS* band. 562 ! TODO MG I could avoid doing this but I have to exchange spin and bands ??? 563 ! For sure there is a better way to do this! 564 call wfd%get_cprj(ib,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 565 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 566 end if 567 568 do jb=ib1,ib2 569 ! Get all <k-q,ib,s|e^{-i(q+G).r}|s,jb,k>, at once. 570 call rho_tw_g(nspinor,npwc,gwc_nfftot,ndat1,gwc_ngfft,1,use_padfft,igfftcg0,gw_gbound,& 571 & ur_sum ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, & 572 & wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,& 573 & nspinor,rhotwg_ki(:,jb)) 574 575 if (Psps%usepaw==1) then 576 ! Add on-site contribution, projectors are already in BZ !TODO Recheck this! 577 i2=jb; if (nspinor==2) i2=(2*jb-1) 578 spad=(nspinor-1) 579 call paw_rho_tw_g(cryst,Pwij_qg, npwc,nspinor,nspinor,Gsph_c%gvec,& 580 Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),rhotwg_ki(:,jb)) 581 end if 582 583 ! Multiply by the square root of the Coulomb term. 584 ! In 3-D systems, the factor sqrt(4pi) is included) 585 do ii=1,nspinor 586 spad = (ii-1) * npwc 587 rhotwg_ki(spad+1:spad+npwc,jb) = rhotwg_ki(spad+1:spad+npwc,jb)*vc_sqrt_qbz(1:npwc) 588 end do 589 590 ! === Treat analytically the case q --> 0 === 591 ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma 592 ! while the Colulomb term is integrated out 593 ! * In the scalar case we have nonzero contribution only if ib==jb 594 ! * For nspinor==2 evalute <ib,up|jb,up> and <ib,dwn|jb,dwn>, 595 ! impose orthonormalization since npwwfn might be < npwvec. 596 if (ik_bz == jk_bz) then 597 if (nspinor == 1) then 598 rhotwg_ki(1, jb) = czero_gw 599 if (ib==jb) rhotwg_ki(1, jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) 600 else 601 npw_k = Wfd%npwarr(ik_ibz) 602 rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero 603 if (ib == jb) then 604 ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, spin, wave_sum, msg) == 0, msg) 605 cg_sum => wave_sum%ug 606 ABI_CHECK(wfd%get_wave_ptr(jb, jk_ibz, spin, wave_jb, msg) == 0, msg) 607 cg_jb => wave_jb%ug 608 ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1) 609 rhotwg_ki(1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp) 610 ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1) 611 rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp) 612 ! PAW is missing 613 614 !rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half) 615 !rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half) 616 end if 617 end if 618 end if 619 end do !jb Got all matrix elements from minbnd up to maxbnd. 620 621 do kb=ib1,ib2 622 ! Get the ket \Sigma|\phi_{k,kb}> according to the method. 623 rhotwgp(:) = rhotwg_ki(:,kb) 624 sigc_ket = czero_gw 625 626 ! SEX part. TODO add check on theta_mu_minus_e0i 627 do ispinor=1,nspinor 628 spadc = (ispinor-1) * npwc 629 call XGEMV('N',npwc,npwc,cone_gw,epsm1_qbz(:,:,1),npwc,rhotwgp(1+spadc:),1,czero_gw,sigsex,1) 630 631 sigsex(:)= -theta_mu_minus_e0i*sigsex(:) 632 633 do io=1,nomega_tot ! nomega==1 as SEX is energy independent. 634 sigc_ket(spadc+1:spadc+npwc,io) = sigsex(:) 635 end do 636 end do 637 638 ! Loop over the non-zero row elements of this column. 639 ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS. 640 ! 2) If gwcalctyp>=20: 641 ! * Only off-diagonal elements connecting states with same character. 642 ! * Only the upper triangle if HF, SEX, or COHSEX. 643 do irow=1,Sigcij_tab(spin)%col(kb)%size1 644 jb = Sigcij_tab(spin)%col(kb)%bidx(irow) 645 rhotwg = rhotwg_ki(:,jb) 646 647 ! Calculate <\phi_j|\Sigma_c|\phi_k> 648 ! Different freqs according to method (AC or Perturbative), see nomega_sigc. 649 do iab=1,Sigp%nsig_ab 650 spadc1=spinor_padc(1,iab); spadc2=spinor_padc(2,iab) 651 do io=1,nomega_sigc 652 sigctmp(io,iab) = XDOTC(npwc,rhotwg(spadc1+1:),1,sigc_ket(spadc2+1:,io),1) 653 end do 654 end do 655 656 ! TODO: save wf1swf2_g to avoid having to recalculate it at each q-point. 657 if (mod10==SIG_COHSEX) then 658 ! Evaluate Static COH. TODO add spinor. 659 if (coh_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank) then 660 ! COH term is done only once for each k-point. 661 ! It does not depend on the index ib summed over. 662 coh_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank 663 664 #if 1 665 call calc_wfwfg(ktabr(:,jk_ibz),jik, spinrot_kgw, & ! TODO why jk_ibz? 666 & gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g) 667 #else 668 ABI_CHECK(jik==1,"jik") 669 call calc_wfwfg(ktabr(:,jk_bz),jik, spinrot_kgw, & 670 gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g) 671 #endif 672 673 if (Psps%usepaw==1) then 674 i1=jb; i2=kb 675 if (nspinor==2) then 676 i1=(2*jb-1); i2=(2*kb-1) 677 end if 678 spad=(nspinor-1) 679 call paw_rho_tw_g(cryst,Pwij_fft,gwc_nfftot,Sigp%nsig_ab,nspinor,& 680 gw_gfft,Cprj_kgw(:,i1:i1+spad),Cprj_kgw(:,i2:i2+spad),wf1swf2_g) 681 end if 682 683 call calc_coh(nspinor,Sigp%nsig_ab,gwc_nfftot,gwc_ngfft,npwc,Gsph_c%gvec,wf1swf2_g,epsm1_qbz(:,:,1),& 684 vc_sqrt_qbz,Vcp%i_sz,iq_ibz,(jb==kb),sigcohme) 685 686 do io=1,nomega_sigc ! Should be 1 687 sigctmp(io,:) = sigctmp(io,:)+sigcohme(:) 688 end do 689 690 end if 691 end if ! COHSEX 692 693 ! Accumulate and, in case, symmetrize matrix elements of Sigma_c. 694 do iab=1,Sigp%nsig_ab 695 is_idx = spin; if (nspinor==2) is_idx=iab 696 697 sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + & 698 & (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab)) 699 700 sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp* sigctmp(:,iab) 701 sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab)) 702 ! TODO this should be the contribution coming from the anti-hermitian part. 703 end do 704 end do !jb used to calculate matrix elements of $\Sigma$ 705 706 end do !kb to calculate matrix elements of $\Sigma$ 707 end do !ib 708 709 ! Deallocate k-dependent quantities. 710 ABI_FREE(gw_gbound) 711 if (Psps%usepaw==1) then 712 call pawpwij_free(Pwij_qg) 713 ABI_FREE(Pwij_qg) 714 end if 715 end do !ik_bz 716 717 ABI_FREE(wfr_bdgw) 718 if (Wfd%usepaw==1) then 719 call pawcprj_free(Cprj_kgw) 720 ABI_FREE(Cprj_kgw) 721 end if 722 end do !spin 723 724 ABI_FREE(igfftcg0) 725 726 ! Gather contributions from all the CPUs. 727 call xmpi_sum(sigcme_tmp, wfd%comm, ierr) 728 call xmpi_sum(sigc, wfd%comm, ierr) 729 730 ! Multiply by constants 731 ! For 3D systems sqrt(4pi) is included in vc_sqrt_qbz. 732 sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz) 733 sigc = sigc /(Cryst%ucvol*Kmesh%nbz) 734 735 ! If we have summed over the IBZ_q now we have to average over degenerate states. 736 ! Presently only diagonal terms are considered 737 ! TODO it does not work if nspinor==2. 738 do spin=1,nsppol 739 if (can_symmetrize(spin)) then 740 ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, sigp%nsig_ab)) 741 sym_cme=czero 742 743 ! Average over degenerate diagonal elements. 744 ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results. 745 ! another good reason to use a strict criterion for the tolerance on eigenvalues. 746 do ib=ib1,ib2 747 ndegs=0 748 do jb=ib1,ib2 749 if (degtab(ib,jb,spin)==1) then 750 if (nspinor == 1) then 751 sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:, :, jb, jb, spin), dim=1) 752 else 753 do ii=1,sigp%nsig_ab 754 sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:, :, jb, jb, ii), dim=1) 755 end do 756 end if 757 758 end if 759 ndegs = ndegs + degtab(ib,jb,spin) 760 end do 761 sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs 762 end do 763 764 if (Sigp%gwcalctyp >= 20) then 765 call esymm_symmetrize_mels(QP_sym(spin),ib1,ib2,sigc(:,1,:,:,spin),sym_cme(1,:,:,1)) 766 end if 767 768 ! Copy symmetrized values. 769 do ib=ib1,ib2 770 do jb=ib1,ib2 771 if (nspinor == 1) then 772 sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1) 773 else 774 sigcme_tmp(:,ib,jb,:) = sym_cme(:,ib,jb,:) 775 end if 776 end do 777 end do 778 ABI_FREE(sym_cme) 779 end if 780 end do 781 782 ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX) 783 if (Sigp%gwcalctyp>=20 .and. sigma_is_herm(Sigp) ) then 784 ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!") 785 do spin=1,nsppol 786 do io=1,nomega_sigc 787 call hermitianize(sigcme_tmp(io,:,:,spin),"Upper") 788 end do 789 end do 790 end if 791 792 ! =========================== 793 ! ==== Deallocate memory ==== 794 ! =========================== 795 if (Psps%usepaw==1) then 796 if (allocated(gw_gfft)) then 797 ABI_FREE(gw_gfft) 798 end if 799 call pawcprj_free(Cprj_ksum) 800 ABI_FREE(Cprj_ksum) 801 if (allocated(Pwij_fft)) then 802 call pawpwij_free(Pwij_fft) 803 ABI_FREE(Pwij_fft) 804 end if 805 end if 806 807 ABI_FREE(ktabr) 808 ABI_FREE(ur_sum) 809 ABI_FREE(rhotwg_ki) 810 ABI_FREE(rhotwg) 811 ABI_FREE(rhotwgp) 812 ABI_FREE(vc_sqrt_qbz) 813 ABI_FREE(sigc_ket) 814 ABI_FREE(epsm1_qbz) 815 ABI_FREE(sigctmp) 816 ABI_FREE(sigc) 817 ABI_FREE(sigsex) 818 ABI_FREE(proc_distrb) 819 if (mod10==SIG_COHSEX) then 820 ABI_FREE(wf1swf2_g) 821 ABI_FREE(coh_distrb) 822 end if 823 824 if (allocated(degtab)) then 825 ABI_FREE(degtab) 826 end if 827 828 call timab(495,2,tsec) ! csigme(SigC) 829 call timab(491,2,tsec) 830 call timab(423,2,tsec) ! cohsex_me 831 832 DBG_EXIT("COLL") 833 834 end subroutine cohsex_me
ABINIT/m_cohsex [ Modules ]
NAME
m_cohsex
FUNCTION
Calculate diagonal and off-diagonal matrix elements of the SEX or COHSEX self-energy operator.
COPYRIGHT
Copyright (C) 1999-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_cohsex 23 24 use defs_basis 25 use m_defs_ptgroups 26 use m_gwdefs 27 use m_xmpi 28 use m_errors 29 use m_abicore 30 31 use defs_datatypes, only : pseudopotential_type, ebands_t 32 use m_time, only : timab 33 use m_fstrings, only : sjoin, itoa 34 use m_hide_blas, only : xdotc, xgemv 35 use m_numeric_tools, only : hermitianize, imin_loc 36 use m_geometry, only : normv 37 use m_crystal, only : crystal_t 38 use m_bz_mesh, only : kmesh_t, findqg0, littlegroup_t 39 use m_gsphere, only : gsphere_t 40 use m_fft_mesh, only : get_gftt, rotate_fft_mesh, cigfft 41 use m_vcoul, only : vcoul_t 42 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g 43 use m_wfd, only : wfdgw_t, wave_t 44 use m_oscillators, only : rho_tw_g, calc_wfwfg 45 use m_screening, only : epsm1_symmetrizer, get_epsm1, epsilonm1_results 46 use m_esymm, only : esymm_t, esymm_symmetrize_mels, esymm_failed 47 use m_sigma, only : sigma_t, sigma_distribute_bks 48 use m_pawang, only : pawang_type 49 use m_pawtab, only : pawtab_type 50 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap 51 use m_paw_sym, only : paw_symcprj 52 53 implicit none 54 55 private