TABLE OF CONTENTS


ABINIT/cohsex_me [ Functions ]

[ Top ] [ Functions ]

NAME

 cohsex_me

FUNCTION

 Calculate diagonal and off-diagonal matrix elements of the SEX or COHSEX self-energy operator.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (FB, GMR, VO, LR, RWG, MG, RShaltaf)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

 sigmak_ibz=Index of the k-point in the IBZ.
 minbnd, maxbnd= min and Max band index for GW correction (for this k-point)
 iomode=Option defining the file format of the SCR file (Fortran, NETCDF)
 Er <Epsilonm1_results> (see the definition of this structured datatype)
    %mqmem=if 0 use out-of-core method in which a single q-slice of espilon is read inside the loop over k
    %nomega_i=Number of imaginary frequencies.
    %nomega_r=Number of real frequencies.
    %nomega=Total number of frequencies.
 Gsph_c<gsphere_t>= info on the G-sphere for Sigma_x
    %nsym=number of symmetry operations
    %rottb(ng,timrev,nsym)=index of (IS) G where I is the identity or the inversion
      operation and G is one of the ng vectors in reciprocal space
    %timrev=2 if time-reversal symmetry is used, 1 otherwise
    %gvec(3,Sigp%npwc)=integer coordinates of each plane wave in reciprocal space
 ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
 Ltg_k datatype containing information on the little group
 Kmesh <kmesh_t>
    %nbz=Number of points in the BZ
    %nibz=Number of points in IBZ
    %kibz(3,nibz)=k-point coordinates, irreducible Brillouin zone
    %kbz(3,nbz)=k-point coordinates, full Brillouin zone
    %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
    %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered
    %ktabp(nbz)= phase factor associated to tnons
 gwc_ngfft(18)=Information about 3D FFT for the oscillator strengths used for the correlation part,
 Vcp <vcoul_t datatype> containing information on the cutoff technique
    %vc_sqrt(npwc,nqibz)= square-root of the coulombian potential for q-points in the IBZ
 Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
 Pawang <type(pawang_type)>=paw angular mesh and related data
 Psps <type(pseudopotential_type)>=variables related to pseudopotentials
    %usepaw=1 for PAW, 0 for NC pseudopotentials.
 Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
    %ibz=q points where $\tilde\epsilon^{-1}$ has been computed
    %bz(3,nqbz)=coordinates of all q-points in BZ
 Sigp <sigparams_t> (see the definition of this structured datatype)
 Cryst<crystal_t>=Info on unit cell and symmetries
    %natom=number of atoms in unit cell
    %ucvol=unit cell volume
    %nsym=number of symmetry operations
    %typat(natom)=type of each atom
  much slower but it requires less memory
 QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
  eig(Sigp%nbnds,Kmesh%nibz,%nsppol)=KS or QP energies for k-points, bands and spin
  occ(Sigp%nbnds,Kmesh%nibz,%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
  Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
 allQP_sym(%nkibz,%nsppol)<esymm_t>=Datatype collecting data on the irreducible representaions of the
  little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
  Sr=sigma_t (see the definition of this structured datatype)

OUTPUT

NOTES

  1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) is included.
  2) The calculation of energy derivative is based on finite elements.
  3) On the symmetrization of Sigma matrix elements ***/
        If  Sk = k+G0 then  M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G}   (k,q)
        If -Sk = k+G0 then  M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)

     Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out
     that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a
     similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase
     has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ.

  4) The unitary transformation relating wavefunctions
     at symmetric k-points should be taken into account during the symmetrization
     of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
     it is possible to make an invariant by just including all the degenerate states and
     averaging the final results over the degenerate subset.

PARENTS

      sigma

CHILDREN

      calc_coh,calc_wfwfg,epsm1_symmetrizer,esymm_symmetrize_mels,findqg0
      get_bz_item,get_epsm1,get_gftt,gsph_fft_tabs,hermitianize
      littlegroup_print,paw_rho_tw_g,paw_symcprj,pawcprj_alloc,pawcprj_copy
      pawcprj_free,pawpwij_free,pawpwij_init,rho_tw_g,rotate_fft_mesh
      sigma_distribute_bks,timab,wfd_change_ngfft,wfd_get_cprj
      wfd_get_many_ur,wfd_get_ur,wrtout,xgemv,xmpi_sum

SOURCE

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