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-2022 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

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 
25 #include "abi_common.h"
26 
27 MODULE m_prep_calc_ucrpa
28 
29  use defs_basis
30  use m_abicore
31  use m_gwdefs!,        only : czero_gw, cone_gw, j_gw, sigparams_t
32  use m_xmpi
33  use m_defs_ptgroups
34  use m_errors
35 
36  use defs_datatypes,  only : pseudopotential_type, ebands_t
37  use m_time,          only : timab
38  use m_hide_blas,     only : xdotc
39  use m_geometry,      only : normv
40  use m_crystal,       only : crystal_t
41  use m_fft_mesh,      only : rotate_FFT_mesh
42  use m_bz_mesh,       only : kmesh_t, findqg0
43  use m_gsphere,       only : gsphere_t
44  use m_io_tools,      only : flush_unit, open_file
45  use m_vcoul,         only : vcoul_t
46  use m_pawpwij,       only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
47  use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t
48  use m_pawang,        only : pawang_type
49  use m_pawtab,        only : pawtab_type
50  use m_pawfgrtab,     only : pawfgrtab_type
51  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
52  use m_paw_nhat,      only : pawmknhat_psipsi
53  use m_paw_sym,       only : paw_symcprj
54  use m_wfd,           only : wfd_t, wave_t
55  use m_oscillators,   only : rho_tw_g
56  use m_esymm,         only : esymm_t, esymm_failed
57  use m_read_plowannier, only : read_plowannier
58  use m_plowannier, only : plowannier_type,operwan_realspace_type
59 
60  implicit none
61 
62  private
63 
64  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-2022 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.

SOURCE

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