TABLE OF CONTENTS
ABINIT/m_prep_calc_ucrpa [ 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 ]
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