TABLE OF CONTENTS


ABINIT/m_prep_calc_ucrpa [ Modules ]

[ Top ] [ Modules ]

NAME

  m_prep_calc_ucrpa

FUNCTION

 Prepare data for the calculation of U with the CRPA method: oscillators strenghs and k-points.

COPYRIGHT

 Copyright (C) 2006-2018 ABINIT group (BAmadon)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 
29 #include "abi_common.h"
30 
31 MODULE m_prep_calc_ucrpa
32 
33  use defs_basis
34  implicit none
35 
36  private
37 
38  public :: prep_calc_ucrpa

ABINIT/prep_calc_ucrpa [ Functions ]

[ Top ] [ Functions ]

NAME

 prep_calc_ucrpa

FUNCTION

 Prepare data for the calculation of U with the CRPA method: oscillators strenghs and k-points.

COPYRIGHT

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

INPUTS

 sigmak_ibz=Index of the k-point in the IBZ.
 minbnd, maxbnd= min and Max band index for GW correction (for this k-point)
 Gsph_x<gsphere_t>= Info on the G-sphere used for Sigma_x
    %nsym=number of symmetry operations
    %rottb(ng,timrev,nsym)=index of (IS) G where I is the identity or the inversion
      operation and G is one of the ng vectors in reciprocal space
    %timrev=2 if time-reversal symmetry is used, 1 otherwise
    %gvec(3,ng)=integer coordinates of each plane wave in reciprocal space
 ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
 Kmesh <kmesh_t>
    %nbz=Number of points in the BZ
    %nibz=Number of points in IBZ
    %kibz(3,nibz)=k-point coordinates, irreducible Brillouin zone
    %kbz(3,nbz)=k-point coordinates, full Brillouin zone
    %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
    %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered
    %ktabp(nbz)= phase factor associated to tnons
 gwx_ngfft(18)=Information about 3D FFT for the oscillator strengths, see ~abinit/doc/variables/vargs.htm#ngfft
 gwx_nfftot=number of points of the FFT grid for GW wavefunctions
 Vcp <vcoul_t datatype> containing information on the cutoff technique
    %vc_sqrt(npwx,nqibz)= square-root of the coulombian potential for q-points in the IBZ
 Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
 Pawang <type(pawang_type)>=paw angular mesh and related data
 Psps <type(pseudopotential_type)>=variables related to pseudopotentials
    %usepaw=1 for PAW, 0 for NC pseudopotentials.
 Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
    %ibz=q points where $\tilde\epsilon^{-1}$ has been computed
    %bz(3,nqbz)=coordinates of all q-points in BZ
 Sigp <sigparams_t> (see the definition of this structured datatype)
 Cryst<crystal_t>=Info on unit cell and symmetries
    %natom=number of atoms in unit cell
    %ucvol=unit cell volume
    %nsym=number of symmetry operations
    %typat(natom)=type of each atom
  much slower but it requires less memory
 QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
  eig(Sigp%nbnds,Kmesh%nibz,%nsppol)=KS or QP energies for k-points, bands and spin
  occ(Sigp%nbnds,Kmesh%nibz,nsppol)=occupation numbers, for each k point in IBZ, each band and spin
  Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
  allQP_sym(%nkibz,%nsppol)<esymm_t>=Datatype collecting data on the irreducible representaions of the
    little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
 prtvol=Flags governing verbosity level.

OUTPUT

NOTES

  1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) [[cite:Gigy1986]] is included.

  2) On the symmetrization of Sigma matrix elements
     If  Sk = k+G0 then  M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G}   (k,q)
     If -Sk = k+G0 then  M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)

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

  3) the unitary transformation relating wavefunctions
     at symmetric k-points should be taken into account during the symmetrization
     of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
     it is possible to make an invariant by just including all the degenerate states and
     averaging the final results over the degenerate subset. Here we divide the states
     where the QP energies are required into complexes. Note however that this approach is not
     based on group theory, and it might lead to spurious results in case of accidental degeneracies.

PARENTS

      sigma

CHILDREN

      findqg0,flush_unit,get_bz_item,gsph_fft_tabs,paw_cross_rho_tw_g
      paw_rho_tw_g,paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free
      pawmknhat_psipsi,pawpwij_free,pawpwij_init,read_plowannier,rho_tw_g
      rotate_fft_mesh,timab,wfd_change_ngfft,wfd_get_cprj,wfd_get_ur
      wfd_paw_get_aeur,wrtout

SOURCE

136 subroutine prep_calc_ucrpa(sigmak_ibz,ikcalc,itypatcor,minbnd,maxbnd,Cryst,QP_BSt,Sigp,Gsph_x,Vcp,Kmesh,Qmesh,lpawu,&
137 & M1_q_m,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,&
138 & Psps,Wfd,Wfdf,allQP_sym,gwx_ngfft,ngfftf,&
139 & prtvol,pawcross,rhot1_q_m)
140 
141  use defs_basis
142  use defs_datatypes
143  use m_abicore
144  use m_gwdefs!,        only : czero_gw, cone_gw, j_gw, sigparams_t
145  use m_xmpi
146  use m_defs_ptgroups
147  use m_errors
148 
149  use m_time,          only : timab
150  use m_hide_blas,     only : xdotc
151  use m_geometry,      only : normv
152  use m_crystal,       only : crystal_t
153  use m_fft_mesh,      only : rotate_FFT_mesh
154  use m_bz_mesh,       only : kmesh_t, get_BZ_item, findqg0, has_IBZ_item
155  use m_gsphere,       only : gsphere_t, gsph_fft_tabs
156  use m_io_tools,      only : flush_unit, open_file
157  use m_vcoul,         only : vcoul_t
158  use m_pawpwij,       only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
159  use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t
160  use m_pawang,        only : pawang_type
161  use m_pawtab,        only : pawtab_type
162  use m_pawfgrtab,     only : pawfgrtab_type
163  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
164  use m_paw_nhat,      only : pawmknhat_psipsi
165  use m_paw_sym,       only : paw_symcprj
166  use m_wfd,           only : wfd_t, wfd_get_ur, wfd_get_cprj, wfd_change_ngfft, wfd_paw_get_aeur
167  use m_oscillators,   only : rho_tw_g
168  use m_esymm,         only : esymm_t, esymm_failed
169  use m_read_plowannier, only : read_plowannier
170 
171 !This section has been created automatically by the script Abilint (TD).
172 !Do not modify the following lines by hand.
173 #undef ABI_FUNC
174 #define ABI_FUNC 'prep_calc_ucrpa'
175 !End of the abilint section
176 
177  implicit none
178 
179 !Arguments ------------------------------------
180 !scalars
181  integer,intent(in) :: sigmak_ibz,ikcalc,itypatcor,prtvol,lpawu,minbnd,maxbnd,pawcross
182  type(crystal_t),intent(in) :: Cryst
183  type(ebands_t),target,intent(in) :: QP_BSt
184  type(kmesh_t),intent(in) :: Kmesh,Qmesh
185  type(vcoul_t),intent(in) :: Vcp
186  type(gsphere_t),intent(in) :: Gsph_x
187 ! type(littlegroup_t),intent(in) :: Ltg_k
188  type(Pseudopotential_type),intent(in) :: Psps
189  type(sigparams_t),target,intent(in) :: Sigp
190  type(pawang_type),intent(in) :: Pawang
191  type(wfd_t),target,intent(inout) :: Wfd,Wfdf
192 !arrays
193  complex(dpc), intent(out) :: rhot1_q_m(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,2*lpawu+1,2*lpawu+1,sigp%npwx,Qmesh%nibz)
194  complex(dpc), intent(out) :: M1_q_m(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,2*lpawu+1,2*lpawu+1,sigp%npwx,Qmesh%nibz)
195  integer,intent(in) :: gwx_ngfft(18),ngfftf(18)
196  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
197  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
198  type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol)
199  type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
200  type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
201 
202 !Local variables ------------------------------
203 !scalars
204  integer,parameter :: use_pawnhat=0,ider0=0,ndat1=1
205  integer :: bandinf,bandsup
206  integer :: gwcalctyp,izero,ib_sum,ib,ib1,ib2,ig,ig_rot,ii,iik,itim_q,i2
207  integer :: ik_bz,ik_ibz,isym_q,iq_bz,iq_ibz,spin,isym,itypatcor_read,jb,iat
208  integer :: jik,jk_bz,jk_ibz,lcor,m1,m3,nspinor,nsppol,ifft
209  integer :: ibsp,dimcprj_gw
210  integer :: spad
211  integer :: comm
212  integer :: ispinor1,ispinor3,isym_kgw,isym_ki,gwx_mgfft,use_padfft,use_padfftf,gwx_fftalga,gwx_fftalgb
213  integer :: gwx_nfftot,nfftf,mgfftf,ierr
214  integer :: nhat12_grdim
215  real(dp) :: fact_sp,theta_mu_minus_esum,tol_empty,norm,weight
216  complex(dpc) :: ctmp,scprod,ph_mkgwt,ph_mkt
217  logical :: iscompatibleFFT,q_is_gamma
218  character(len=500) :: msg
219 !arrays
220  integer :: g0(3),spinor_padx(2,4)
221  integer,pointer :: igfftxg0(:),igfftfxg0(:)
222  integer,allocatable :: gwx_gfft(:,:),gwx_gbound(:,:),gboundf(:,:)
223  integer,allocatable ::  ktabr(:,:),irottb(:,:),ktabrf(:,:)
224  real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),qbz(3),q0(3),tsec(2)
225  real(dp) :: spinrot_kbz(4),spinrot_kgw(4)
226  real(dp),pointer :: qp_ene(:,:,:),qp_occ(:,:,:)
227  real(dp),allocatable :: nhat12(:,:,:),grnhat12(:,:,:,:)
228  complex(gwpc),allocatable :: vc_sqrt_qbz(:)
229  complex(gwpc),allocatable :: rhotwg_ki(:,:)
230  complex(gwpc),allocatable :: wfr_bdgw(:,:),wfr_sum(:)
231  complex(gwpc),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:)
232  complex(gwpc),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:)
233  complex(gwpc),pointer :: cg_jb(:),cg_sum(:)
234  complex(dpc) :: ovlp(2)
235  complex(dpc),allocatable :: coeffW_BZ(:,:,:,:,:,:)
236  logical :: can_symmetrize(Wfd%nsppol)
237  logical,allocatable :: bks_mask(:,:,:)
238  type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
239  type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
240  type(esymm_t),pointer :: QP_sym(:)
241  logical     :: ecriture=.FALSE.
242  logical     :: l_ucrpa,luwindow
243  integer     :: g0_dump(3),iq_ibz_dump,dumint(2)
244 
245 !************************************************************************
246 
247  l_ucrpa=.true.
248 
249  DBG_ENTER("COLL")
250 
251  !
252  ! === Initial check ===
253  ABI_CHECK(Sigp%npwx==Gsph_x%ng,'')
254 
255  call timab(430,1,tsec) ! csigme (SigX)
256 
257  gwcalctyp=Sigp%gwcalctyp
258  !
259  ! === Initialize MPI variables ===
260  comm = Wfd%comm
261 
262  !
263  ! === Initialize some values ===
264  nspinor = Wfd%nspinor
265  nsppol  = Wfd%nsppol
266  spinor_padx(:,:)=RESHAPE((/0,0,Sigp%npwx,Sigp%npwx,0,Sigp%npwx,Sigp%npwx,0/),(/2,4/))
267 
268  qp_ene => QP_BSt%eig(:,:,:)
269  qp_occ => QP_BSt%occ(:,:,:)
270 
271  ! Exctract the symmetries of the bands for this k-point
272  QP_sym => allQP_sym(sigmak_ibz,1:nsppol)
273 
274  ib1=minbnd
275  ib2=maxbnd
276 
277  ! === Read Wannier function coefficients for Ucrpa
278  ! === for future computation of rhot_m_q directly in this routine.
279 
280  dumint=0
281  luwindow=.true.
282 ! write(6,*) "cc",allocated(coeffW_BZ)
283  call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor_read,Kmesh,lcor,luwindow,&
284 & nspinor,nsppol,pawang,prtvol,dumint)
285 
286  if(lcor/=lpawu) then
287      msg = "lcor and lpawu differ in prep_calc_ucrpa"
288      MSG_ERROR(msg)
289  endif
290 
291  ! === End of read Wannier function coefficients for Ucrpa
292 
293 
294  !
295  ! === Index of the GW point in the BZ array, its image in IBZ and time-reversal ===
296  jk_bz=Sigp%kptgw2bz(ikcalc)
297  !write(6,*) "ikcalc,jk_bz",ikcalc,jk_bz
298  !write(6,*) "ikcalc",Kmesh%bz(:,ikcalc)
299  !write(6,*) "jk_bz",Kmesh%bz(:,jk_bz)
300 ! jk_bz=ikcalc
301  call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
302 ! write(6,*) "jk_ibz",Kmesh%ibz(:,jk_ibz)
303 ! write(6,*) "jk_bz,jk_ibz",jk_bz,jk_ibz,isym_kgw,itim
304  !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
305  spinrot_kgw(:)=Cryst%spinrot(:,isym_kgw)
306  !
307  write(msg,'(2a,3f8.3,a,i4,a,2(i3,a))')ch10,&
308 &  ' Calculating Oscillator element at k= ',kgw, "k-point number",ikcalc,&
309 &  ' bands n = from ',ib1,' to ',ib2,ch10
310  call wrtout(std_out,msg,'COLL')
311 
312 
313  if (ANY(gwx_ngfft(1:3) /= Wfd%ngfft(1:3)) ) then
314    call wfd_change_ngfft(Wfd,Cryst,Psps,gwx_ngfft)
315  end if
316  gwx_mgfft   = MAXVAL(gwx_ngfft(1:3))
317  gwx_fftalga = gwx_ngfft(7)/100
318  gwx_fftalgb = MOD(gwx_ngfft(7),100)/10
319 
320  if (pawcross==1) then
321    mgfftf = MAXVAL(ngfftf(1:3))
322  end if
323 
324  can_symmetrize = .FALSE.
325  if (Sigp%symsigma>0) then
326    can_symmetrize = .TRUE.
327    if (gwcalctyp >= 20) then
328     do spin=1,Wfd%nsppol
329       can_symmetrize(spin) = .not.esymm_failed(QP_sym(spin))
330       if (.not.can_symmetrize(spin)) then
331         write(msg,'(a,i0,4a)')&
332 &         " Symmetrization cannot be performed for spin: ",spin,ch10,&
333 &         " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg)
334         MSG_WARNING(msg)
335       end if
336     end do
337    end if
338    ABI_CHECK(nspinor==1,'Symmetrization with nspinor=2 not implemented')
339  end if
340 
341  ABI_ALLOCATE(rhotwg_ki,(Sigp%npwx*nspinor,minbnd:maxbnd))
342  rhotwg_ki=czero_gw
343  ABI_ALLOCATE(vc_sqrt_qbz,(Sigp%npwx))
344  !
345  ! === Normalization of theta_mu_minus_esum ===
346  ! * If nsppol==2, qp_occ $\in [0,1]$
347  SELECT CASE (nsppol)
348  CASE (1)
349    fact_sp=half; tol_empty=0.01   ! below this value the state is assumed empty
350    if (Sigp%nspinor==2) then
351     fact_sp=one; tol_empty=0.005  ! below this value the state is assumed empty
352    end if
353  CASE (2)
354    fact_sp=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic
355  CASE DEFAULT                    ! spin unpolarized system is treated using nsppol==2
356    MSG_BUG('Wrong nsppol')
357  END SELECT
358 
359  ! Remove empty states from the list of states that will be distributed.
360  ABI_ALLOCATE(bks_mask,(Wfd%mband,Kmesh%nbz,nsppol))
361  bks_mask=.FALSE.
362  do spin=1,nsppol
363    do ik_bz=1,Kmesh%nbz
364      ik_ibz = Kmesh%tab(ik_bz)
365      do ib_sum=1,Sigp%nbnds
366        bks_mask(ib_sum,ik_bz,spin) = (qp_occ(ib_sum,ik_ibz,spin)>=tol_empty)
367      end do
368    end do
369  end do
370 
371 ! ABI_ALLOCATE(proc_distrb,(Wfd%mband,Kmesh%nbz,nsppol))
372 ! call sigma_distribution(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,proc_distrb,bks_mask=bks_mask)
373 ! call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,proc_distrb,bks_mask=bks_mask)
374 ! write(6,*)"lim", ib1,ib2
375 ! do ib_sum= ib1,ib2
376 !   do ik_bz=1, Kmesh%nbz
377 !     write(6,*) ib_sum,ik_bz, proc_distrb(ib_sum,ik_bz,1),Wfd%my_rank
378 !   enddo
379 ! enddo
380 
381  ABI_DEALLOCATE(bks_mask)
382 
383  write(msg,'(a,i8)')" Will sum all (b,k,s) occupied states in Sigma_x for k-point",ikcalc
384  call wrtout(std_out,msg,'PERS')
385  !
386  ! The index of G-G0 in the FFT mesh the oscillators ===
387  ! * Sigp%mG0 gives the MAX G0 component to account for umklapp.
388  ! * Note the size MAX(Sigp%npwx,Sigp%npwc).
389  ABI_ALLOCATE(igfftxg0,(Gsph_x%ng))
390  !
391  ! === Precalculate the FFT index of $ R^{-1}(r-\tau) $ ===
392  ! * S=\transpose R^{-1} and k_BZ = S k_IBZ
393  ! * irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
394  gwx_nfftot = PRODUCT(gwx_ngfft(1:3))
395  ABI_ALLOCATE(irottb,(gwx_nfftot,Cryst%nsym))
396  call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwx_ngfft,irottb,iscompatibleFFT)
397  if (.not.iscompatibleFFT) then
398    msg = "FFT mesh is not compatible with symmetries. Results might be affected by large errors!"
399    MSG_WARNING(msg)
400  end if
401 
402  ABI_ALLOCATE(ktabr,(gwx_nfftot,Kmesh%nbz))
403  do ik_bz=1,Kmesh%nbz
404    isym=Kmesh%tabo(ik_bz)
405    do ifft=1,gwx_nfftot
406      ktabr(ifft,ik_bz)=irottb(ifft,isym)
407    end do
408  end do
409  ABI_DEALLOCATE(irottb)
410 
411  if (Psps%usepaw==1 .and. pawcross==1) then
412    nfftf = PRODUCT(ngfftf(1:3))
413    ABI_ALLOCATE(irottb,(nfftf,Cryst%nsym))
414    call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfftf,irottb,iscompatibleFFT)
415 
416    ABI_ALLOCATE(ktabrf,(nfftf,Kmesh%nbz))
417    do ik_bz=1,Kmesh%nbz
418      isym=Kmesh%tabo(ik_bz)
419      do ifft=1,nfftf
420        ktabrf(ifft,ik_bz)=irottb(ifft,isym)
421      end do
422    end do
423    ABI_DEALLOCATE(irottb)
424  end if
425  !
426  ! === Additional allocations for PAW ===
427  if (Psps%usepaw==1) then
428    ABI_DATATYPE_ALLOCATE(Cprj_ksum,(Cryst%natom,nspinor))
429    call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm)
430 
431    nhat12_grdim=0
432    if (use_pawnhat==1) then ! Compensation charge for \phi_a^*\phi_b
433      call wrtout(std_out,"Using nhat12","COLL")
434      ABI_ALLOCATE(nhat12  ,(2,gwx_nfftot,nspinor**2))
435      ABI_ALLOCATE(grnhat12,(2,gwx_nfftot,nspinor**2,3*nhat12_grdim))
436    end if
437  end if ! usepaw==1
438  !
439 
440  if (Sigp%symsigma>0) then
441    !call littlegroup_print(Ltg_k,std_out,prtvol,'COLL')
442    !
443    ! === Find number of complexes and number of bands in each complex ===
444    ! The tolerance is a little bit arbitrary (0.001 eV)
445    ! It could be reduced, in particular in case of nearly accidental degeneracies
446 !   if (ANY(degtab/=0)) then ! If two states do not belong to the same complex => matrix elements of v_xc differ
447 !     write(msg,'(a,3f8.3,a)')' Degenerate states at k-point = ( ',kgw(:),' ).'
448 !     call wrtout(std_out,msg,'COLL')
449 !     do spin=1,nsppol
450 !       do ib=ib1,ib2
451 !         do jb=ib+1,ib2
452 !           if (degtab(ib,jb,spin)==1) then
453 !             write(msg,'(a,i2,a,i4,a,i4)')' (spin ',spin,')',ib,' <====> ',jb
454 !             call wrtout(std_out,msg,'COLL')
455 !             if (ABS(Sr%vxcme(ib,jk_ibz,spin)-Sr%vxcme(jb,jk_ibz,spin))>ABS(tol6*Sr%vxcme(jb,jk_ibz,spin))) then
456 !               write(msg,'(7a)')&
457 !&                ' It seems that an accidental degeneracy is occurring at this k-point ',ch10,&
458 !&                ' In this case, using symsigma=1 might lead to spurious results as the algorithm ',ch10,&
459 !&                ' will treat these states as degenerate, and it won''t be able to remove the degeneracy. ',ch10,&
460 !&                ' In order to avoid this deficiency, run the calculation using symsigma=0'
461 !               MSG_WARNING(msg)
462 !             end if
463 !           end if
464 !         end do
465 !       end do
466 !     end do
467 !   end if
468  end if !symsigma
469 
470  ABI_ALLOCATE(wfr_sum,(gwx_nfftot*nspinor))
471  if (pawcross==1) then
472    ABI_ALLOCATE(ur_ae_sum,(nfftf*nspinor))
473    ABI_ALLOCATE(ur_ae_onsite_sum,(nfftf*nspinor))
474    ABI_ALLOCATE(ur_ps_onsite_sum,(nfftf*nspinor))
475  end if