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