TABLE OF CONTENTS


ABINIT/cchi0 [ Functions ]

[ Top ] [ Functions ]

NAME

 cchi0

FUNCTION

 Main calculation of the independent-particle susceptibility chi0 for qpoint!=0

INPUTS

 use_tr=If .TRUE. valence states are stored in Wfs_val and only resonant transitions are calculated
  (time reversal is assumed)
 Dtset <type(dataset_type)>=all input variables in this dataset
 Cryst<crystal_t>= data type gathering info on symmetries and unit cell
    %natom=number of atoms
    %nsym=number of symmetries
    %xred(3,natom)=reduced coordinated of atoms
    %typat(natom)=type of each atom
    %rprimd(3,3)=dimensional primitive translations in real space (bohr)
    %timrev= 2 if time reversal can be used, 1 otherwise
 qpoint(3)=reciprocal space coordinates of the q wavevector
 Ep<type(em1params_t_type)>= Parameters related to the calculation of the inverse dielectric matrix.
    %nbnds=number of bands summed over
    %npwe=number of planewaves for the irreducible polarizability X^0_GGp
    %npwvec=maximum number of G vectors
     used to define the dimension of some arrays e.g igfft
    %nsppol=1 for unpolarized, 2 for spin-polarized
    %nomega=total number of frequencies in X^0 (both real and imaginary)
    %nomegasf=number of real frequencies used to sample the imaginary part of X^0 (spectral method)
    %spmeth=1 if we use the spectral method, 0 for standard Adler-Wiser expression
    %spsmear=gaussian broadening used to approximate the delta distribution
    %zcut=small imaginary shift to avoid poles in X^0
 Psps <type(pseudopotential_type)>=variables related to pseudopotentials
 Kmesh <kmesh_t>= datatype gathering parameters related to the k-point sampling
    %nibz=number of k-points in the IBZ
    %nbz=number of k-points in the BZ
    %bz(3,nbz)=reduced coordinates for k-points in the full Brillouin zone
    %ibz(3,nibz)=reduced coordinates for k-points in the irreducible wedge
    %tab(nbz)=mapping between a kpt in the BZ (array bz) and the irred point in the array ibz
    %tabi(nbz)= -1 if inversion is needed to obtain this particular kpt in the BZ, 1 means identity
    %tabo(nbz)= for each point in the BZ, the index of the symmetry operation S in reciprocal
      space which rotates k_IBZ onto \pm k_BZ (depending on tabi)
    %tabp(nbz)= For each k_BZ, it gives the phase factors associated to non-symmorphic operations, i.e
      e^{-i 2 \pi k_IBZ \cdot R{^-1}t} == e{-i 2\pi k_BZ cdot t} where :
      \transpose R{-1}=S and (S k_IBZ) = \pm k_BZ (depending on ktabi)
    %tabr(nfftot,nbz) For each point r on the real mesh and for each k-point in the BZ, tabr
      gives the index of (R^-1 (r-t)) in the FFT array where R=\transpose S^{-1} and k_BZ=S k_IBZ.
      t is the fractional translation associated to R
 Gsph_epsG0<gsphere_t data type> The G-sphere used to describe chi0/eps. (including umklapp G0 vectors)
    %ng=number of G vectors for chi0
    %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G
    %phmGt(Ep%npwe,nsym)=phase factors e^{-iG \cdot t} needed to symmetrize oscillator matrix elements and epsilon
    %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
    %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
 nbvw=number of bands in the arrays wfrv
 ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths (see input variable)
 nfftot_gw=Total number of points in the GW FFT grid
 Ltg_q<Little group>=Data type gathering information on the little group of the q-points.
 Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
 Pawang<pawang_type> angular mesh discretization and related data:
 QP_BSt<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities)
   %mband=MAX number of bands over k-points and spin (==Ep%nbnds)
   %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band
   %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes
  Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
  Wfd<wfd_t>=Object used to access the wavefunctions

OUTPUT

  chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix at wavevector qpoint and
   each frequeny defined by Ep%omega and Ep%nomega

PARENTS

      screening

CHILDREN

      accumulate_chi0sumrule,approxdelta,assemblychi0_sym,assemblychi0sf
      calc_wfwfg,chi0_bbp_mask,completechi0_deltapart,cwtime,flush_unit
      get_bz_diff,get_bz_item,get_gftt,gsph_fft_tabs,gsph_free,gsph_in_fftbox
      hilbert_transform,littlegroup_print,make_transitions,paw_cross_rho_tw_g
      paw_rho_tw_g,paw_symcprj,pawcprj_alloc,pawcprj_free,pawpwij_free
      pawpwij_init,read_plowannier,rho_tw_g,setup_spectral
      symmetrize_afm_chi0,timab,wfd_change_ngfft,wfd_distribute_kb_kpbp
      wfd_get_cprj,wfd_get_ur,wfd_paw_get_aeur,wrtout,xmpi_sum

SOURCE

1142 subroutine cchi0(use_tr,Dtset,Cryst,qpoint,Ep,Psps,Kmesh,QP_BSt,Gsph_epsG0,&
1143 & Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,nbvw,ngfft_gw,nfftot_gw,ngfftf,nfftf_tot,&
1144 & chi0,ktabr,ktabrf,Ltg_q,chi0_sumrule,Wfd,Wfdf)
1145 
1146 
1147 !This section has been created automatically by the script Abilint (TD).
1148 !Do not modify the following lines by hand.
1149 #undef ABI_FUNC
1150 #define ABI_FUNC 'cchi0'
1151 !End of the abilint section
1152 
1153  implicit none
1154 
1155 !Arguments ------------------------------------
1156 !scalars
1157  integer,intent(in) :: nbvw,nfftot_gw,nfftf_tot
1158  logical,intent(in) :: use_tr
1159  type(ebands_t),target,intent(in) :: QP_BSt
1160  type(kmesh_t),intent(in) :: Kmesh
1161  type(crystal_t),intent(in) :: Cryst
1162  type(Dataset_type),intent(in) :: Dtset
1163  type(em1params_t),intent(in) :: Ep
1164  type(gsphere_t),intent(in) :: Gsph_epsG0
1165  type(littlegroup_t),intent(in) :: Ltg_q
1166  type(Pawang_type),intent(in) :: Pawang
1167  type(Pseudopotential_type),intent(in) :: Psps
1168  type(wfd_t),target,intent(inout) :: Wfd,Wfdf
1169 !arrays
1170  integer,intent(in) :: ktabr(nfftot_gw,Kmesh%nbz),ktabrf(nfftf_tot*Dtset%pawcross,Kmesh%nbz)
1171  integer,intent(in) :: ngfft_gw(18),ngfftf(18)
1172  real(dp),intent(in) :: qpoint(3)
1173  real(dp),intent(out) :: chi0_sumrule(Ep%npwe)
1174  complex(gwpc),intent(out) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega)
1175  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
1176  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
1177  type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
1178  type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
1179 
1180 !Local variables ------------------------------
1181 !scalars
1182  integer,parameter :: tim_fourdp1=1,two_poles=2,one_pole=1,ndat1=1
1183  integer :: bandinf,bandsup,dim_rtwg,band1,band2,ierr
1184  integer :: ig1,ig2,iat,ik_bz,ik_ibz,ikmq_bz,ikmq_ibz
1185  integer :: io,iomegal,iomegar,ispinor1,ispinor2,isym_k,itypatcor,nfft
1186  integer :: isym_kmq,itim_k,itim_kmq,m1,m2,my_wl,my_wr,size_chi0
1187  integer :: nfound,nkpt_summed,nspinor,nsppol,mband
1188  integer :: comm,gw_mgfft,use_padfft,gw_fftalga,lcor,mgfftf,use_padfftf
1189  integer :: my_nbbp,my_nbbpks,spin,nbmax,dummy
1190  real(dp) :: cpu_time,wall_time,gflops
1191  real(dp) :: deltaeGW_b1kmq_b2k,deltaeGW_enhigh_b2k,deltaf_b1kmq_b2k
1192  real(dp) :: e_b1_kmq,en_high,fac,fac2,fac3,f_b1_kmq,factor,max_rest,min_rest,my_max_rest
1193  real(dp) :: my_min_rest,numerator,spin_fact,weight,wl,wr
1194  real(dp) :: gw_gsq,memreq
1195  complex(dpc) :: ph_mkmqt,ph_mkt
1196  complex(gwpc) :: local_czero_gw
1197  logical :: qzero,isirred_k,isirred_kmq,luwindow
1198  character(len=500) :: msg,allup
1199  type(gsphere_t) :: Gsph_FFT
1200 !arrays
1201  integer :: G0(3),umklp_k(3),umklp_kmq(3)
1202  integer :: ucrpa_bands(2)
1203  integer :: wtk_ltg(Kmesh%nbz),got(Wfd%nproc)
1204  integer,allocatable :: tabr_k(:),tabr_kmq(:),tabrf_k(:),tabrf_kmq(:)
1205  integer,allocatable :: igfftepsG0(:),gspfft_igfft(:),igfftepsG0f(:)
1206  integer,allocatable :: gw_gfft(:,:),gw_gbound(:,:),dummy_gbound(:,:),gboundf(:,:)
1207  integer,allocatable :: bbp_ks_distrb(:,:,:,:)
1208  real(dp) :: kbz(3),kmq_bz(3),spinrot_k(4),spinrot_kmq(4),q0(3),tsec(2)
1209  real(dp),ABI_CONTIGUOUS pointer :: qp_energy(:,:,:),qp_occ(:,:,:)
1210  real(dp),allocatable :: omegasf(:)
1211  complex(dpc),allocatable :: green_enhigh_w(:),green_w(:),kkweight(:,:)
1212  complex(gwpc),allocatable :: sf_chi0(:,:,:),rhotwg(:)
1213  complex(gwpc),allocatable :: ur1_kmq_ibz(:),ur2_k_ibz(:),wfwfg(:)
1214  complex(gwpc),allocatable :: usr1_kmq(:),ur2_k(:)
1215  complex(gwpc),allocatable :: ur_ae1(:),ur_ae_onsite1(:),ur_ps_onsite1(:)
1216  complex(gwpc),allocatable :: ur_ae2(:),ur_ae_onsite2(:),ur_ps_onsite2(:)
1217  complex(dpc), allocatable :: coeffW_BZ(:,:,:,:,:,:)
1218  logical,allocatable :: bbp_mask(:,:)
1219  type(pawcprj_type),allocatable :: Cprj1_kmq(:,:),Cprj2_k(:,:)
1220  type(pawpwij_t),allocatable :: Pwij(:),Pwij_fft(:)
1221 !************************************************************************
1222 
1223  DBG_ENTER("COLL")
1224 
1225  call timab(331,1,tsec) ! cchi0
1226  call cwtime(cpu_time,wall_time,gflops,"start")
1227 
1228  nsppol = Wfd%nsppol; nspinor = Wfd%nspinor
1229  ucrpa_bands(1)=dtset%ucrpa_bands(1)
1230  ucrpa_bands(2)=dtset%ucrpa_bands(2)
1231  luwindow=.false.
1232  if(abs(dtset%ucrpa_window(1)+1_dp)>tol8.and.(abs(dtset%ucrpa_window(2)+1_dp)>tol8)) then
1233    luwindow=.true.
1234  endif
1235 ! write(6,*)"ucrpa_bands",ucrpa_bands
1236 ! write(6,*)"ucrpa_window",dtset%ucrpa_window
1237 ! write(6,*)"luwindow",luwindow
1238 
1239 !  For cRPA calculation of U: read forlb.ovlp
1240  if(dtset%ucrpa>=1) then
1241    call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor,Kmesh,lcor,luwindow,&
1242 & nspinor,nsppol,pawang,dtset%prtvol,ucrpa_bands)
1243  endif
1244 ! End of reading forlb.ovlp
1245 
1246  if ( ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,ngfft_gw)
1247  gw_mgfft = MAXVAL(ngfft_gw(1:3))
1248  gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10)
1249 
1250  if (Dtset%pawcross==1) mgfftf = MAXVAL(ngfftf(1:3))
1251 
1252  ! == Copy some values ===
1253  comm = Wfd%comm
1254  mband   = Wfd%mband
1255  nfft    = Wfd%nfft
1256  ABI_CHECK(Wfd%nfftot==nfftot_gw,"Wrong nfftot_gw")
1257 
1258  dim_rtwg=1 !; if (nspinor==2) dim_rtwg=2  ! can reduce size depending on Ep%nI and Ep%nj
1259  size_chi0 = Ep%npwe*Ep%nI*Ep%npwe*Ep%nJ*Ep%nomega
1260 
1261  qp_energy => QP_BSt%eig; qp_occ => QP_BSt%occ
1262 
1263  ! Initialize the completeness correction
1264  if (Ep%gwcomp==1) then
1265    en_high=MAXVAL(qp_energy(Ep%nbnds,:,:)) + Ep%gwencomp
1266    write(msg,'(a,f8.2,a)')' Using completeness correction with the energy ',en_high*Ha_eV,' [eV]'
1267    call wrtout(std_out,msg,'COLL')
1268 
1269    ! Allocation of wfwfg and green_enhigh_w moved inside openmp loop
1270    ! Init the largest G-sphere contained in the FFT box for the wavefunctions.
1271    call gsph_in_fftbox(Gsph_FFT,Cryst,Wfd%ngfft)
1272 
1273    !call print_gsphere(Gsph_FFT,unit=std_out,prtvol=10)
1274 
1275    ABI_MALLOC(gspfft_igfft,(Gsph_FFT%ng))
1276    ABI_MALLOC(dummy_gbound,(2*gw_mgfft+8,2))
1277 
1278    ! Mapping between G-sphere and FFT box.
1279    call gsph_fft_tabs(Gsph_FFT,(/0,0,0/),Wfd%mgfft,Wfd%ngfft,dummy,dummy_gbound,gspfft_igfft)
1280    ABI_FREE(dummy_gbound)
1281 
1282    if (Psps%usepaw==1) then  ! * Prepare the onsite contributions on the GW FFT mesh.
1283      ABI_MALLOC(gw_gfft,(3,nfft))
1284      q0=zero
1285      call get_gftt(ngfft_gw,q0,Cryst%gmet,gw_gsq,gw_gfft) ! Get the set of plane waves in the FFT Box.
1286      ABI_DT_MALLOC(Pwij_fft,(Psps%ntypat))
1287      call pawpwij_init(Pwij_fft,nfft,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
1288    end if
1289  end if
1290 
1291  ! Setup weights (2 for spin unpolarized sistem, 1 for polarized).
1292  ! spin_fact is used to normalize the occupation factors to one. Consider also the AFM case.
1293  SELECT CASE (nsppol)
1294  CASE (1)
1295    weight=two/Kmesh%nbz; spin_fact=half
1296    if (Wfd%nspden==2) then
1297     weight=one/Kmesh%nbz; spin_fact=half
1298    end if
1299    if (nspinor==2) then
1300     weight=one/Kmesh%nbz; spin_fact=one
1301    end if
1302  CASE (2)
1303    weight=one/Kmesh%nbz; spin_fact=one
1304  CASE DEFAULT
1305    MSG_BUG("Wrong nsppol")
1306  END SELECT
1307 
1308  ! Weight for points in the IBZ_q.
1309  wtk_ltg(:)=1
1310  if (Ep%symchi==1) then
1311    do ik_bz=1,Ltg_q%nbz
1312      wtk_ltg(ik_bz)=0
1313      if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only k-points in the IBZ_q.
1314      wtk_ltg(ik_bz)=SUM(Ltg_q%wtksym(:,:,ik_bz))
1315    end do
1316  end if
1317 
1318  write(msg,'(a,i2,2a,i2)')&
1319 &  ' Using spectral method for the imaginary part = ',Ep%spmeth,ch10,&
1320 &  ' Using symmetries to sum only over the IBZ_q  = ',Ep%symchi
1321  call wrtout(std_out,msg,'COLL')
1322 
1323  if (use_tr) then
1324    ! Special care has to be taken in metals and/or spin dependent systems
1325    ! as Wfs_val might contain unoccupied states.
1326    call wrtout(std_out,' Using faster algorithm based on time reversal symmetry.','COLL')
1327  else
1328    call wrtout(std_out,' Using slow algorithm without time reversal symmetry.')
1329  end if
1330 
1331  ! TODO this table can be calculated for each k-point
1332  my_nbbpks=0; allup="All"; got=0
1333  ABI_MALLOC(bbp_ks_distrb,(mband,mband,Kmesh%nbz,nsppol))
1334  ABI_MALLOC(bbp_mask,(mband,mband))
1335 
1336  do spin=1,nsppol
1337    do ik_bz=1,Kmesh%nbz
1338      if (Ep%symchi==1) then
1339        if (Ltg_q%ibzq(ik_bz)/=1) CYCLE  ! Only IBZ_q
1340      end if
1341 
1342      ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz.
1343      call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k)
1344 
1345      ! Get index of k-q in the BZ, stop if not found as the weight=one/nkbz is not correct.
1346      call get_BZ_diff(Kmesh,kbz,qpoint,ikmq_bz,g0,nfound)
1347      ABI_CHECK(nfound==1,"Check kmesh")
1348 
1349      ! Get ikmq_ibz, non-symmorphic phase, ph_mkmqt, and symmetries from ikmq_bz.
1350      call get_BZ_item(Kmesh,ikmq_bz,kmq_bz,ikmq_ibz,isym_kmq,itim_kmq)
1351 
1352      call chi0_bbp_mask(Ep,use_tr,QP_BSt,mband,ikmq_ibz,ik_ibz,spin,spin_fact,bbp_mask)
1353 
1354      call wfd_distribute_kb_kpbp(Wfd,ikmq_ibz,ik_ibz,spin,allup,my_nbbp,bbp_ks_distrb(:,:,ik_bz,spin),got,bbp_mask)
1355      my_nbbpks = my_nbbpks + my_nbbp
1356    end do
1357  end do
1358 
1359  ABI_FREE(bbp_mask)
1360 
1361  write(msg,'(a,i0,a)')" Will sum ",my_nbbpks," (b,b',k,s) states in chi0."
1362  call wrtout(std_out,msg,'PERS')
1363 
1364  if (Psps%usepaw==1) then
1365    ABI_DT_MALLOC(Pwij,(Psps%ntypat))
1366    call pawpwij_init(Pwij,Ep%npwepG0,qpoint,Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
1367    ! Allocate statements moved to inside openmp loop
1368  end if
1369 
1370  SELECT CASE (Ep%spmeth)
1371  CASE (0)
1372    call wrtout(std_out,' Calculating chi0(q,omega,G,G")','COLL')
1373    ! Allocation of green_w moved inside openmp loop
1374 
1375  CASE (1, 2)
1376    call wrtout(std_out,' Calculating Im chi0(q,omega,G,G")','COLL')
1377 
1378    ! Find Max and min resonant transitions for this q, report also treated by this proc.
1379    call make_transitions(Wfd,1,Ep%nbnds,nbvw,nsppol,Ep%symchi,Cryst%timrev,GW_TOL_DOCC,&
1380 &    max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,qp_energy,qp_occ,qpoint,bbp_ks_distrb)
1381    !
1382    ! Calculate frequency dependent weights for Hilbert transform.
1383    ABI_MALLOC(omegasf,(Ep%nomegasf))
1384    ABI_MALLOC(kkweight,(Ep%nomegasf,Ep%nomega))
1385    !my_wl=1; my_wr=Ep%nomegasf
1386    call setup_spectral(Ep%nomega,Ep%omega,Ep%nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,&
1387 &    0,Ep%zcut,zero,my_wl,my_wr,kkweight)
1388 
1389    if (.not.use_tr) then
1390      MSG_BUG('spectral method requires time-reversal')
1391    end if
1392 
1393    memreq = two*gwpc*Ep%npwe**2*(my_wr-my_wl+1)*b2Gb
1394    write(msg,'(a,f10.4,a)')' memory required per spectral point: ',two*gwpc*Ep%npwe**2*b2Mb,' [Mb]'
1395    call wrtout(std_out,msg,'PERS')
1396    write(msg,'(a,f10.4,a)')' memory required by sf_chi0: ',memreq,' [Gb]'
1397    call wrtout(std_out,msg,'PERS')
1398    if (memreq > two) then
1399      MSG_WARNING(' Memory required for sf_chi0 is larger than 2.0 Gb!')
1400    end if
1401    ABI_STAT_MALLOC(sf_chi0,(Ep%npwe,Ep%npwe,my_wl:my_wr), ierr)
1402    ABI_CHECK(ierr==0, 'out-of-memory in sf_chi0')
1403    sf_chi0=czero_gw
1404 
1405  CASE DEFAULT
1406    MSG_BUG("Wrong spmeth")
1407  END SELECT
1408 
1409  nkpt_summed=Kmesh%nbz
1410  if (Ep%symchi==1) then
1411    nkpt_summed=Ltg_q%nibz_ltg
1412    call littlegroup_print(Ltg_q,std_out,Dtset%prtvol,'COLL')
1413  end if
1414 
1415  write(msg,'(a,i6,a)')' Calculation status : ',nkpt_summed,' to be completed '
1416  call wrtout(std_out,msg,'COLL')
1417 
1418  ! ============================================
1419  ! === Begin big fat loop over transitions ===
1420  ! ============================================
1421  chi0=czero_gw; chi0_sumrule=zero
1422 
1423  ! === Loop on spin to calculate trace $\chi_{up,up}+\chi_{down,down}$ ===
1424  ! Only $\chi_{up,up} for AFM.
1425  do spin=1,nsppol
1426    if (ALL(bbp_ks_distrb(:,:,:,spin) /= Wfd%my_rank)) CYCLE
1427 
1428    ! Allocation of arrays that are private to loop
1429    if (Ep%gwcomp==1)  then
1430      ABI_MALLOC(wfwfg,(nfft*nspinor**2))
1431    end if
1432    if (Ep%gwcomp==1)  then
1433      ABI_MALLOC(green_enhigh_w,(Ep%nomega))
1434    end if
1435    if (Ep%spmeth==0)  then
1436      ABI_MALLOC(green_w,(Ep%nomega))
1437    end if
1438    if (Psps%usepaw==1) then
1439      ABI_DT_MALLOC(Cprj2_k  ,(Cryst%natom,nspinor))
1440      call pawcprj_alloc(Cprj2_k,  0,Wfd%nlmn_atm)
1441      ABI_DT_MALLOC(Cprj1_kmq,(Cryst%natom,nspinor))
1442      call pawcprj_alloc(Cprj1_kmq,0,Wfd%nlmn_atm)
1443      if (Dtset%pawcross==1) then
1444        ABI_MALLOC(ur_ae1,(nfftf_tot*nspinor))
1445        ABI_MALLOC(ur_ae_onsite1,(nfftf_tot*nspinor))
1446        ABI_MALLOC(ur_ps_onsite1,(nfftf_tot*nspinor))
1447        ABI_MALLOC(ur_ae2,(nfftf_tot*nspinor))
1448        ABI_MALLOC(ur_ae_onsite2,(nfftf_tot*nspinor))
1449        ABI_MALLOC(ur_ps_onsite2,(nfftf_tot*nspinor))
1450        ABI_MALLOC(igfftepsG0f,(Ep%npwepG0))
1451        ABI_MALLOC(tabrf_k,(nfftf_tot))
1452        ABI_MALLOC(tabrf_kmq,(nfftf_tot))
1453      end if
1454    end if
1455 
1456    ABI_MALLOC(rhotwg,(Ep%npwepG0*dim_rtwg))
1457    ABI_MALLOC(tabr_k,(nfft))
1458    ABI_MALLOC(tabr_kmq,(nfft))
1459    ABI_MALLOC(ur1_kmq_ibz,(nfft*nspinor))
1460    ABI_MALLOC(ur2_k_ibz,(nfft*nspinor))
1461    ABI_MALLOC(usr1_kmq,(nfft*nspinor))
1462    ABI_MALLOC(ur2_k,   (nfft*nspinor))
1463    ABI_MALLOC(igfftepsG0,(Ep%npwepG0))
1464 
1465    ! Loop over k-points in the BZ.
1466    do ik_bz=1,Kmesh%nbz
1467 
1468      if (Ep%symchi==1) then
1469        if (Ltg_q%ibzq(ik_bz)/=1) CYCLE  ! Only IBZ_q
1470      end if
1471 
1472      if (ALL(bbp_ks_distrb(:,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE
1473 
1474      write(msg,'(2(a,i4),a,i2,a,i3)')' ik= ',ik_bz,'/',Kmesh%nbz,' spin= ',spin,' done by mpi rank:',Wfd%my_rank
1475      call wrtout(std_out,msg,'PERS')
1476 
1477      ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz.
1478      call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt,umklp_k,isirred_k)
1479 
1480      call get_BZ_diff(Kmesh,kbz,qpoint,ikmq_bz,G0,nfound)
1481      if (nfound==0) then
1482        MSG_ERROR("Cannot find kbz - qpoint in Kmesh")
1483      end if
1484 
1485      ! Get ikmq_ibz, non-symmorphic phase, ph_mkmqt, and symmetries from ikmq_bz.
1486      call get_BZ_item(Kmesh,ikmq_bz,kmq_bz,ikmq_ibz,isym_kmq,itim_kmq,ph_mkmqt,umklp_kmq,isirred_kmq)
1487 
1488 !BEGIN DEBUG
1489      !if (ANY(umklp_k /=0)) then
1490      !  write(msg,'(a,3i2)')" umklp_k /= 0 ",umklp_k
1491      !  MSG_ERROR(msg)
1492      !end if
1493      !if (ANY( g0 /= -umklp_kmq + umklp_k) ) then
1494      !if (ANY( g0 /= -umklp_kmq ) ) then
1495      !  write(msg,'(a,3(1x,3i2))')" g0 /= -umklp_kmq + umklp_k ",g0, umklp_kmq, umklp_k
1496      !  MSG_ERROR(msg)
1497      !end if
1498      !g0 = -umklp_k + umklp_kmq
1499      !g0 = +umklp_k - umklp_kmq
1500      !if (ANY (ABS(g0) > Ep%mg0) ) then
1501      !  write(msg,'(a,6(1x,i0))')"  ABS(g0) > Ep%mg0 ",g0,Ep%mg0
1502      !  MSG_ERROR(msg)
1503      !end if
1504 !END DEBUG
1505 
1506      ! Copy tables for rotated FFT points
1507      tabr_k(:)  =ktabr(:,ik_bz)
1508      spinrot_k(:)=Cryst%spinrot(:,isym_k)
1509 
1510      tabr_kmq(:)=ktabr(:,ikmq_bz)
1511      spinrot_kmq(:)=Cryst%spinrot(:,isym_kmq)
1512 
1513      if (Dtset%pawcross==1) then
1514        tabrf_k(:)  =ktabrf(:,ik_bz)
1515        tabrf_kmq(:)=ktabrf(:,ikmq_bz)
1516      end if
1517      !
1518      ! Tables for the FFT of the oscillators.
1519      !  a) FFT index of G-G0.
1520      !  b) gw_gbound table for the zero-padded FFT performed in rhotwg.
1521      ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2))
1522      call gsph_fft_tabs(Gsph_epsG0,g0,gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igfftepsG0)
1523      if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
1524      if (use_padfft==0) then
1525        ABI_FREE(gw_gbound)
1526        ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft))
1527      end if
1528 
1529      if (Dtset%pawcross==1) then
1530         ABI_MALLOC(gboundf,(2*mgfftf+8,2))
1531        call gsph_fft_tabs(Gsph_epsG0,g0,mgfftf,ngfftf,use_padfftf,gboundf,igfftepsG0f)
1532        if (ANY(gw_fftalga == [2, 4])) use_padfftf=0
1533        if (use_padfftf==0) then
1534          ABI_FREE(gboundf)
1535          ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf))
1536        end if
1537      end if
1538 
1539      nbmax=Ep%nbnds
1540      do band1=1,nbmax ! Loop over "conduction" states.
1541        if (ALL(bbp_ks_distrb(band1,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE
1542 
1543        call wfd_get_ur(Wfd,band1,ikmq_ibz,spin,ur1_kmq_ibz)
1544 
1545        if (Psps%usepaw==1) then
1546          call wfd_get_cprj(Wfd,band1,ikmq_ibz,spin,Cryst,Cprj1_kmq,sorted=.FALSE.)
1547          call paw_symcprj(ikmq_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_kmq)
1548          if (Dtset%pawcross==1) then
1549            call wfd_paw_get_aeur(Wfdf,band1,ikmq_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae1,ur_ae_onsite1,ur_ps_onsite1)
1550          end if
1551        end if
1552 
1553        e_b1_kmq=qp_energy(band1,ikmq_ibz,spin)
1554        f_b1_kmq=   qp_occ(band1,ikmq_ibz,spin)
1555 
1556        do band2=1,nbmax ! Loop over "valence" states.
1557 !debug         if (.not.luwindow.AND.dtset%ucrpa==1.AND.band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)&
1558 !debug&                                            .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) CYClE
1559        !write(6,*) "ik,band1,band2",ik_bz,band1,band2
1560          if (luwindow.AND.dtset%ucrpa==1.AND.((QP_Bst%eig(band1,ik_ibz   ,spin)-QP_Bst%fermie)<=dtset%ucrpa_window(2)) &
1561 &                                       .AND.((QP_Bst%eig(band1,ik_ibz   ,spin)-QP_Bst%fermie)>=dtset%ucrpa_window(1)) &
1562 &                                       .AND.((QP_Bst%eig(band2,ikmq_ibz,spin)-QP_Bst%fermie)<=dtset%ucrpa_window(2)) &
1563 &                                       .AND.((QP_Bst%eig(band2,ikmq_ibz,spin)-QP_Bst%fermie)>=dtset%ucrpa_window(1))) CYCLE
1564 
1565          if (bbp_ks_distrb(band1,band2,ik_bz,spin) /= Wfd%my_rank) CYCLE
1566 
1567          deltaf_b1kmq_b2k=spin_fact*(f_b1_kmq-qp_occ(band2,ik_ibz,spin))
1568 
1569          if (Ep%gwcomp==0) then ! Skip negligible transitions.
1570            if (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC) CYCLE
1571          else
1572            ! When the completeness correction is used,
1573            ! we need to also consider transitions with vanishing deltaf
1574            !if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC) CYCLE
1575            !
1576            ! Rangel This is to compute chi correctly when using the extrapolar method
1577            if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC .and. (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC .or. band1<band2)) CYCLE
1578          end if
1579 
1580          deltaeGW_b1kmq_b2k=e_b1_kmq-qp_energy(band2,ik_ibz,spin)
1581 
1582          call wfd_get_ur(Wfd,band2,ik_ibz,spin,ur2_k_ibz)
1583 
1584          if (Psps%usepaw==1) then
1585            call wfd_get_cprj(Wfd,band2,ik_ibz,spin,Cryst,Cprj2_k,sorted=.FALSE.)
1586            call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj2_k)
1587            if (Dtset%pawcross==1) then
1588              call wfd_paw_get_aeur(Wfdf,band2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae2,ur_ae_onsite2,ur_ps_onsite2)
1589            end if
1590          end if
1591 
1592          SELECT CASE (Ep%spmeth)
1593          CASE (0)
1594            ! Standard Adler-Wiser expression.
1595            ! Add the small imaginary of the Time-Ordered RF only for non-zero real omega ! FIXME What about metals?
1596            if (.not.use_tr) then
1597              ! Have to sum over all possible resonant and anti-resonant transitions.
1598              do io=1,Ep%nomega
1599                green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,one_pole)
1600              end do
1601 
1602            else
1603              if (Ep%gwcomp==0) then ! cannot be completely skipped in case of completeness correction
1604                if (band1<band2) CYCLE ! Here we GAIN a factor ~2
1605              end if
1606 
1607              do io=1,Ep%nomega
1608                !Rangel: In metals, the intra-band transitions term does not contain the antiresonant part
1609                !green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0)
1610                if (band1==band2) then
1611                  green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,one_pole)
1612                else
1613                  green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,two_poles)
1614                end if
1615 
1616                if (Ep%gwcomp==1) then ! Calculate the completeness correction
1617                  numerator= -spin_fact*qp_occ(band2,ik_ibz,spin)
1618                  deltaeGW_enhigh_b2k = en_high-qp_energy(band2,ik_ibz,spin)
1619 
1620                  if (REAL(Ep%omega(io))<GW_TOL_W0) then ! Completeness correction is NOT valid for real frequencies
1621                    green_enhigh_w(io) = g0g0w(Ep%omega(io),numerator,deltaeGW_enhigh_b2k,Ep%zcut,GW_TOL_W0,two_poles)
1622                  else
1623                    green_enhigh_w(io) = local_czero_gw
1624                  end if
1625                  !
1626                  !Rangel Correction for metals
1627                  !if (deltaf_b1kmq_b2k<0.d0) then
1628                  if (band1>=band2 .and. ABS(deltaf_b1kmq_b2k) > GW_TOL_DOCC ) then
1629                    green_w(io)= green_w(io) - green_enhigh_w(io)
1630                  else ! Disregard green_w, since it is already accounted for through the time-reversal
1631                    green_w(io)=             - green_enhigh_w(io)
1632                  end if
1633                end if !gwcomp==1
1634              end do !io
1635 
1636              if (Ep%gwcomp==1.and.band1==band2) then
1637                ! Add the "delta part" of the extrapolar method. TODO doesnt work for spinor
1638                call calc_wfwfg(tabr_k,itim_k,spinrot_k,nfft,nspinor,ngfft_gw,ur2_k_ibz,ur2_k_ibz,wfwfg)
1639 
1640                if (Psps%usepaw==1) then
1641                  call paw_rho_tw_g(nfft,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,gw_gfft,&
1642 &                  Cprj2_k,Cprj2_k,Pwij_fft,wfwfg)
1643 
1644                  ! Add PAW cross term
1645                  if (Dtset%pawcross==1) then
1646                    call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,&
1647 &                   ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_kmq,tabrf_kmq,ph_mkmqt,spinrot_kmq,&
1648 &                   ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k  ,tabrf_k  ,ph_mkt  ,spinrot_k,dim_rtwg,wfwfg)
1649                  end if
1650                end if
1651 
1652                qzero=.FALSE.
1653                call completechi0_deltapart(ik_bz,qzero,Ep%symchi,Ep%npwe,Gsph_FFT%ng,Ep%nomega,nspinor,&
1654 &                nfft,ngfft_gw,gspfft_igfft,gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0)
1655 
1656              end if
1657            end if ! use_tr
1658 
1659          CASE (1, 2)
1660            ! Spectral method, WARNING time-reversal here is always assumed!
1661            if (deltaeGW_b1kmq_b2k<0) CYCLE
1662            call approxdelta(Ep%nomegasf,omegasf,deltaeGW_b1kmq_b2k,Ep%spsmear,iomegal,iomegar,wl,wr,Ep%spmeth)
1663          END SELECT
1664 
1665          ! Form rho-twiddle(r)=u^*_{b1,kmq_bz}(r) u_{b2,kbz}(r) and its FFT transform.
1666          call rho_tw_g(nspinor,Ep%npwepG0,nfft,ndat1,ngfft_gw,1,use_padfft,igfftepsG0,gw_gbound,&
1667 &          ur1_kmq_ibz,itim_kmq,tabr_kmq,ph_mkmqt,spinrot_kmq,&
1668 &          ur2_k_ibz,  itim_k  ,tabr_k  ,ph_mkt  ,spinrot_k,dim_rtwg,rhotwg)
1669 
1670          if (Psps%usepaw==1) then
1671            ! Add PAW on-site contribution, projectors are already in the BZ.
1672            call paw_rho_tw_g(Ep%npwepG0,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,&
1673 &           Cprj1_kmq,Cprj2_k,Pwij,rhotwg)
1674 
1675            ! Add PAW cross term
1676            if (Dtset%pawcross==1) then
1677              call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,&
1678 &             ur_ae1,ur_ae_onsite1,ur_ps_onsite1,itim_kmq,tabrf_kmq,ph_mkmqt,spinrot_kmq,&
1679 &             ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k  ,tabrf_k  ,ph_mkt  ,spinrot_k,dim_rtwg,rhotwg)
1680            end if
1681          end if
1682 
1683          SELECT CASE (Ep%spmeth)
1684 
1685          CASE (0) ! Adler-Wiser.
1686 !debug           if(dtset%ucrpa==2)  then
1687            if(dtset%ucrpa>=1.and..not.luwindow)  then
1688              fac=one
1689              fac2=one
1690              fac3=one
1691              m1=-1
1692              m2=-1
1693              call flush_unit(std_out)
1694              if(dtset%ucrpa<=2) then
1695                if (       band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)&
1696 &                    .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then
1697                  do iat=1, cryst%nattyp(itypatcor)
1698                    do ispinor1=1,nspinor
1699                      do ispinor2=1,nspinor
1700                        do m1=1,2*lcor+1
1701                          do m2=1,2*lcor+1
1702                            fac=fac - real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
1703 &                                     conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* &
1704 &                                     coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor2,m2)*&
1705 &                                     conjg(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor2,m2)))
1706                          enddo
1707                        enddo
1708                      enddo
1709                    enddo
1710                  enddo
1711                  if(dtset%ucrpa==1) fac=zero
1712                endif
1713              else if (dtset%ucrpa>=3) then
1714                if (band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)) then
1715                  do iat=1, cryst%nattyp(itypatcor)
1716                    do ispinor1=1,nspinor
1717                      do m1=1,2*lcor+1
1718                        fac2=fac2-real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
1719 &                                 conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)))
1720                      enddo
1721                    enddo
1722                  enddo
1723                  if(dtset%ucrpa==4) fac2=zero
1724                endif
1725                if (band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then
1726                  do iat=1, cryst%nattyp(itypatcor)
1727                    do ispinor1=1,nspinor
1728                      do m1=1,2*lcor+1
1729                        fac3=fac3-real(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor1,m1)*&
1730 &                                 conjg(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor1,m1)))
1731                      enddo
1732                    enddo
1733                  enddo
1734                  if(dtset%ucrpa==4) fac3=zero
1735                endif
1736                fac=real(fac2*fac3)
1737              endif
1738 
1739 !             if(dtset%prtvol>=10) write(6,'(6i3,e15.5,a)') ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q=/0"
1740 !             if(dtset%prtvol>=10.and.abs(fac-one)>0.00001) &
1741 !&             write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q/=0"
1742 !             if(dtset%prtvol>=10) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q/=0"
1743              green_w=green_w*fac
1744            endif
1745 
1746            call assemblychi0_sym(ik_bz,nspinor,Ep,Ltg_q,green_w,Ep%npwepG0,rhotwg,Gsph_epsG0,chi0)
1747 
1748          CASE (1, 2)
1749            ! Spectral method (not yet adapted for nspinor=2)
1750            call assemblychi0sf(ik_bz,Ep%symchi,Ltg_q,Ep%npwepG0,Ep%npwe,rhotwg,Gsph_epsG0,&
1751 &            deltaf_b1kmq_b2k,my_wl,iomegal,wl,my_wr,iomegar,wr,Ep%nomegasf,sf_chi0)
1752 
1753          CASE DEFAULT
1754            MSG_BUG("Wrong spmeth")
1755          END SELECT
1756 
1757          ! Accumulating the sum rule on chi0. Eq. (5.284) in G.D. Mahan Many-Particle Physics 3rd edition. [[cite:Mahan2000]]
1758          ! TODO Does not work with spinor
1759          factor=spin_fact*qp_occ(band2,ik_ibz,spin)
1760          call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_b1kmq_b2k,&
1761 &          Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule)
1762 
1763          ! Include also the completeness correction in the sum rule
1764          if (Ep%gwcomp==1) then
1765            factor=-spin_fact*qp_occ(band2,ik_ibz,spin)
1766            call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_enhigh_b2k,&
1767 &            Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule)
1768            if (band1==Ep%nbnds) then
1769              chi0_sumrule(:)=chi0_sumrule(:) + wtk_ltg(ik_bz)*spin_fact*qp_occ(band2,ik_ibz,spin)*deltaeGW_enhigh_b2k
1770            end if
1771          end if
1772 
1773        end do !band2
1774      end do !band1
1775 
1776      ABI_FREE(gw_gbound)
1777      if (Dtset%pawcross==1) then
1778        ABI_FREE(gboundf)
1779      end if
1780    end do !ik_bz
1781 
1782    ! Deallocation of arrays private to the spin loop.
1783    ABI_FREE(igfftepsG0)
1784    ABI_FREE(ur1_kmq_ibz)
1785    ABI_FREE(ur2_k_ibz)
1786    ABI_FREE(usr1_kmq)
1787    ABI_FREE(ur2_k)
1788    ABI_FREE(rhotwg)
1789    ABI_FREE(tabr_k)
1790    ABI_FREE(tabr_kmq)
1791    if (allocated(green_w)) then
1792      ABI_FREE(green_w)
1793    end if
1794    if (allocated(wfwfg)) then
1795      ABI_FREE(wfwfg)
1796    end if
1797    if (allocated(green_enhigh_w)) then
1798      ABI_FREE(green_enhigh_w)
1799    end if
1800    if (Psps%usepaw==1) then
1801      call pawcprj_free(Cprj2_k)
1802      ABI_DT_FREE(Cprj2_k)
1803      call pawcprj_free(Cprj1_kmq)
1804      ABI_DT_FREE(Cprj1_kmq)
1805      if (Dtset%pawcross==1) then
1806        ABI_FREE(ur_ae1)
1807        ABI_FREE(ur_ae_onsite1)
1808        ABI_FREE(ur_ps_onsite1)
1809        ABI_FREE(ur_ae2)
1810        ABI_FREE(ur_ae_onsite2)
1811        ABI_FREE(ur_ps_onsite2)
1812        ABI_FREE(tabrf_k)
1813        ABI_FREE(tabrf_kmq)
1814        ABI_FREE(gboundf)
1815        ABI_FREE(igfftepsG0f)
1816      end if
1817    end if
1818  end do !spin
1819 
1820  ! After big loop over transitions, now MPI
1821  ! Master took care of the contribution in case of metallic|spin polarized systems.
1822  SELECT CASE (Ep%spmeth)
1823  CASE (0)
1824    ! Adler-Wiser
1825    ! Collective sum of the contributions of each node.
1826    ! Looping on frequencies to avoid problems with the size of the MPI packet
1827    do io=1,Ep%nomega
1828      call xmpi_sum(chi0(:,:,io),comm,ierr)
1829    end do
1830 
1831  CASE (1, 2)
1832    ! Spectral method.
1833    call hilbert_transform(Ep%npwe,Ep%nomega,Ep%nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,Ep%spmeth)
1834 
1835    ! Deallocate here before xmpi_sum
1836    if (allocated(sf_chi0)) then
1837      ABI_FREE(sf_chi0)
1838    end if
1839 
1840    ! Collective sum of the contributions.
1841    ! Looping over frequencies to avoid problems with the size of the MPI packet
1842    do io=1,Ep%nomega
1843      call xmpi_sum(chi0(:,:,io),comm,ierr)
1844    end do
1845 
1846  CASE DEFAULT
1847    MSG_BUG("Wrong spmeth")
1848  END SELECT
1849 
1850  ! Divide by the volume
1851 !$OMP PARALLEL WORKSHARE
1852    chi0=chi0*weight/Cryst%ucvol
1853 !$OMP END PARALLEL WORKSHARE
1854 
1855  ! === Collect the sum rule ===
1856  ! * The pi factor comes from Im[1/(x-ieta)] = pi delta(x)
1857  call xmpi_sum(chi0_sumrule,comm,ierr)
1858  chi0_sumrule=chi0_sumrule*pi*weight/Cryst%ucvol
1859  !
1860  ! *************************************************
1861  ! **** Now each node has chi0(q,G,Gp,Ep%omega) ****
1862  ! *************************************************
1863 
1864  ! Impose Hermiticity (valid only for zero or purely imaginary frequencies)
1865  ! MG what about metals, where we have poles around zero?
1866  do io=1,Ep%nomega
1867    if (ABS(REAL(Ep%omega(io))) <0.00001) then
1868      do ig2=1,Ep%npwe
1869        do ig1=1,ig2-1
1870          chi0(ig2,ig1,io) = GWPC_CONJG(chi0(ig1,ig2,io))
1871        end do
1872      end do
1873    end if
1874  end do
1875 
1876  ! === Symmetrize chi0 in case of AFM system ===
1877  ! Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$.
1878  ! Works only in case of magnetic group Shubnikov type IV.
1879  if (Cryst%use_antiferro) call symmetrize_afm_chi0(Cryst,Gsph_epsG0,Ltg_q,Ep%npwe,Ep%nomega,chi0)
1880 
1881  ! =====================
1882  ! ==== Free memory ====
1883  ! =====================
1884  ABI_FREE(bbp_ks_distrb)
1885 
1886  if (allocated(gw_gfft)) then
1887    ABI_FREE(gw_gfft)
1888  end if
1889  if (allocated(kkweight)) then
1890    ABI_FREE(kkweight)
1891  end if
1892  if (allocated(omegasf)) then
1893    ABI_FREE(omegasf)
1894  end if
1895  if (allocated(gspfft_igfft)) then
1896    ABI_FREE(gspfft_igfft)
1897  end if
1898 
1899  call gsph_free(Gsph_FFT)
1900 
1901  ! deallocation for PAW.
1902  if (Psps%usepaw==1) then
1903    call pawpwij_free(Pwij)
1904    ABI_DT_FREE(Pwij)
1905    if (allocated(Pwij_fft)) then
1906      call pawpwij_free(Pwij_fft)
1907      ABI_DT_FREE(Pwij_fft)
1908    end if
1909  end if
1910 
1911  if(dtset%ucrpa>=1) then
1912    ABI_DEALLOCATE(coeffW_BZ)
1913  endif
1914 
1915  call timab(331,2,tsec)
1916  call cwtime(cpu_time,wall_time,gflops,"stop")
1917  write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time
1918 
1919  DBG_EXIT("COLL")
1920 
1921 end subroutine cchi0

ABINIT/cchi0q0 [ Functions ]

[ Top ] [ Functions ]

NAME

 cchi0q0

FUNCTION

 Calculate chi0 in the limit q-->0

INPUTS

  use_tr=If .TRUE. Wfs_val are allocate and only resonant transitions are evaluated (assumes time reversal symmetry)
  Dtset <type(dataset_type)>=all input variables in this dataset
  Ep= datatype gathering differening parameters related to the calculation of the inverse dielectric matrix
  Gsph_epsG0<gvectors_data_type>: Info on the G-sphere used to describe chi0/espilon (including umklapp)
    %ng=number of G vectors
    %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G  in the array gvec
    %phmGt(ng,nsym)=phase factor e^{-iG.\tau} needed to symmetrize oscillator matrix elements and chi0
    %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
    %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
  Ep%inclvkb=flag to include (or not) the grad of Vkb
  Ltg_q= little group datatype
  nbvw=number of bands in the arrays wfrv,wfgv
  Kmesh<kmesh_t> The k-point mesh
   %kbz(3,nbz)=k-point coordinates, full Brillouin zone
   %tab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
   irreducible point (kIBZ), where kBZ= (IS) kIBZ and I is either the inversion or the identity
   %tabi(nbzx)= for each point in the BZ defines whether inversion  has to be
   considered in the relation kBZ=(IS) kIBZ (1 => only S; -1 => -S)
   %tabo(nbzx)= the symmetry operation S that takes kIBZ to each kBZ
   %tabp(nbzx)= phase factor associated to tnons e^{-i 2 \pi k\cdot R{^-1}t}
  ktabr(nfftot_gw,Kmesh%nbz) index of R^-(r-t) in the FFT array, where k_BZ = (IS) k_IBZ and S = \transpose R^{-1}
  Ep%nbnds=number of bands
  ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths.
  Ep%nomega=number of frequencies
  Cryst<crystal_t>= data type gathering info on symmetries and unit cell
   %natom=number of atoms
   %nsym=number of symmetry operations
   %symrec(3,3,nsym)=symmetry operations in reciprocal space
   %typat(natom)=type of each atom
   %xred(3,natom)=reduced coordinated of atoms
   %rprimd(3,3)=dimensional primitive translations in real space (bohr)
   %timrev=2 if time-reversal symmetry can be used, 1 otherwise
  Ep%npwe=number of planewaves for sigma exchange (input variable)
  nfftot_gw=Total number of points in the GW FFT grid
  Ep%omega(Ep%nomega)=frequencies
  Psps <type(pseudopotential_type)>=variables related to pseudopotentials
     %mpsang=1+maximum angular momentum for nonlocal pseudopotential
  Pawang<pawang_type> angular mesh discretization and related data:
  Pawrad(ntypat*usepaw)<Pawrad_type>=paw radial mesh and related data
  Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels
  QP_BSt<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities)
    %mband=MAX number of bands over k-points and spin (==Ep%nbnds)
    %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band
    %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes
  KS_BSt<ebands_t>=KS energies and occupations.
    %eig(mband,nkpt,nsppol)=KS energies
  Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
  Wfd<wfd_t>=Object used to access the wavefunctions

OUTPUT

  chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix for wavevector qq,
   and frequencies defined by Ep%omega
  chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)= Lower wings
  chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)= Upper wings
  chi0_head(3,3,Ep%nomega)=Head of chi0.

NOTES

  *) The terms "head", "wings" and "body" of chi(G,Gp) refer to
     G=Gp=0, either G or Gp=0, and neither=0 respectively

  *) Symmetry conventions:
      1) symmetry in real space is defined as: R_t f(r) = f(R^-1(r-t))
      2) S=\transpose R^-1
      3) kbz=S kibz

  The wavefunctions for the k-point in the BZ are (assuming nondegenerate states):

  u(G,b, Sk) = u ( S^-1G,b,k)* e^{-i(Sk+G)*t)
  u(G,b,-Sk) = u*(-S^-1G,b,k)* e^{ i(Sk-G)*t)

  u(r,b, Sk) = u (R^-1(r-t),b,k) e^{-iSk*t}
  u(r,b,-Sk) = u*(R^-1(r-t),b,k) e^{ iSK*t}

  The gradient of Vnl(K,Kp) for the k-point in the BZ should be:
   gradvnl(SG,SGp,Sk)=S gradvnl(G,Gp,kibz)
 /***********************************************************************/

TODO

  Check npwepG0 before Switching on umklapp

PARENTS

      screening

CHILDREN

      accumulate_chi0_q0,accumulate_chi0sumrule,accumulate_sfchi0_q0
      approxdelta,calc_wfwfg,chi0_bbp_mask,completechi0_deltapart,cwtime
      flush_unit,get_bz_item,get_gftt,gsph_fft_tabs,gsph_free,gsph_in_fftbox
      hilbert_transform,hilbert_transform_headwings,littlegroup_print
      make_transitions,paw_cross_ihr_comm,paw_cross_rho_tw_g,paw_rho_tw_g
      paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawhur_free
      pawhur_init,pawpwij_free,pawpwij_init,print_gsphere,read_plowannier
      rho_tw_g,setup_spectral,symmetrize_afm_chi0,vkbr_free,vkbr_init
      wfd_change_ngfft,wfd_distribute_bbp,wfd_get_cprj,wfd_get_ur
      wfd_paw_get_aeur,wrtout,xmpi_sum

SOURCE

 185 subroutine cchi0q0(use_tr,Dtset,Cryst,Ep,Psps,Kmesh,QP_BSt,KS_BSt,Gsph_epsG0,&
 186 &  Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,Pawfgrtab,Paw_onsite,ktabr,ktabrf,nbvw,ngfft_gw,&
 187 &  nfftot_gw,ngfftf,nfftf_tot,chi0,chi0_head,chi0_lwing,chi0_uwing,Ltg_q,chi0_sumrule,Wfd,Wfdf)
 188 
 189 
 190 !This section has been created automatically by the script Abilint (TD).
 191 !Do not modify the following lines by hand.
 192 #undef ABI_FUNC
 193 #define ABI_FUNC 'cchi0q0'
 194 !End of the abilint section
 195 
 196  implicit none
 197 
 198 !Arguments ------------------------------------
 199 !scalars
 200  integer,intent(in) :: nbvw,nfftot_gw,nfftf_tot
 201  logical,intent(in) :: use_tr
 202  type(ebands_t),target,intent(in) :: QP_BSt,KS_BSt
 203  type(crystal_t),intent(in) :: Cryst
 204  type(Dataset_type),intent(in) :: Dtset
 205  type(littlegroup_t),intent(in) :: Ltg_q
 206  type(em1params_t),intent(in) :: Ep
 207  type(kmesh_t),intent(in) :: Kmesh
 208  type(gsphere_t),intent(in) :: Gsph_epsG0
 209  type(Pseudopotential_type),intent(in) :: Psps
 210  type(Pawang_type),intent(in) :: Pawang
 211  type(wfd_t),target,intent(inout) :: Wfd,Wfdf
 212 !arrays
 213  integer,intent(in) :: ktabr(nfftot_gw,Kmesh%nbz),ktabrf(nfftf_tot*Dtset%pawcross,Kmesh%nbz)
 214  integer,intent(in) :: ngfft_gw(18),ngfftf(18)
 215  real(dp),intent(out) :: chi0_sumrule(Ep%npwe)
 216  complex(gwpc),intent(out) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega)
 217  complex(dpc),intent(out) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)
 218  complex(dpc),intent(out) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)
 219  complex(dpc),intent(out) :: chi0_head(3,3,Ep%nomega)
 220  type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw)
 221  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
 222  type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw)
 223  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
 224  type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
 225  type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
 226 
 227 !Local variables ------------------------------
 228 !scalars
 229  integer,parameter :: tim_fourdp=1,enough=10,two_poles=2,one_pole=1,ndat1=1
 230  integer :: bandinf,bandsup,lcor,nspinor,npw_k,istwf_k,mband,nfft
 231  integer :: band1,band2,iat,ig,itim_k,ik_bz,ik_ibz,io,iqlwl,ispinor1,ispinor2,isym_k
 232  integer :: itypatcor,m1,m2,nkpt_summed,dim_rtwg,use_padfft,gw_fftalga,use_padfftf,mgfftf
 233  integer :: my_nbbp,my_nbbpks,spin,nsppol !ig1,ig2,
 234  integer :: comm,ierr,my_wl,my_wr,iomegal,iomegar,gw_mgfft,dummy
 235  real(dp) :: cpu_time,wall_time,gflops
 236  real(dp) :: fac,fac1,fac2,fac3,fac4,spin_fact,deltaf_b1b2,weight,factor
 237  real(dp) :: max_rest,min_rest,my_max_rest,my_min_rest
 238  real(dp) :: en_high,deltaeGW_enhigh_b2
 239  real(dp) :: wl,wr,numerator,deltaeGW_b1b2
 240  real(dp) :: gw_gsq,memreq
 241  complex(dpc) :: deltaeKS_b1b2
 242  logical :: qzero,luwindow
 243  character(len=500) :: msg_tmp,msg,allup
 244  type(gsphere_t) :: Gsph_FFT
 245 !arrays
 246  integer,ABI_CONTIGUOUS pointer :: kg_k(:,:)
 247  integer :: ucrpa_bands(2)
 248  integer :: wtk_ltg(Kmesh%nbz)
 249  integer :: got(Wfd%nproc)
 250  integer,allocatable :: tabr_k(:),tabrf_k(:)
 251  integer,allocatable :: igffteps0(:),gspfft_igfft(:),igfftepsG0f(:)
 252  integer,allocatable :: gw_gfft(:,:),gw_gbound(:,:),dummy_gbound(:,:),gboundf(:,:)
 253  integer,allocatable :: bbp_ks_distrb(:,:,:,:)
 254  real(dp) :: kbz(3),spinrot_kbz(4),q0(3)
 255  real(dp),ABI_CONTIGUOUS pointer :: ks_energy(:,:,:),qp_energy(:,:,:),qp_occ(:,:,:)
 256  real(dp),allocatable :: omegasf(:)
 257  complex(gwpc) :: rhotwx(3,Wfd%nspinor**2)
 258  complex(gwpc),allocatable :: rhotwg(:)
 259  complex(dpc),allocatable :: green_w(:),green_enhigh_w(:)
 260  complex(dpc),allocatable :: sf_lwing(:,:,:),sf_uwing(:,:,:),sf_head(:,:,:)
 261  complex(dpc) :: wng(3),chq(3)
 262  complex(dpc) :: ph_mkt
 263  complex(gwpc),allocatable :: sf_chi0(:,:,:)
 264  complex(dpc),allocatable :: kkweight(:,:)
 265  complex(gwpc),allocatable :: ur1_kibz(:),ur2_kibz(:)
 266  complex(gwpc),allocatable :: usr1_k(:),ur2_k(:)
 267  complex(gwpc),allocatable :: wfwfg(:)
 268  complex(gwpc),allocatable :: ur_ae1(:),ur_ae_onsite1(:),ur_ps_onsite1(:)
 269  complex(gwpc),allocatable :: ur_ae2(:),ur_ae_onsite2(:),ur_ps_onsite2(:)
 270  complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:),ug2(:)
 271  complex(dpc), allocatable :: coeffW_BZ(:,:,:,:,:,:)
 272  logical :: gradk_not_done(Kmesh%nibz)
 273  logical,allocatable :: bbp_mask(:,:)
 274  type(pawcprj_type),allocatable :: Cprj1_bz(:,:),Cprj2_bz(:,:)
 275  type(pawcprj_type),allocatable :: Cprj1_ibz(:,:),Cprj2_ibz(:,:)
 276  type(pawpwij_t),allocatable :: Pwij(:),Pwij_fft(:)
 277  type(pawhur_t),allocatable :: Hur(:)
 278  type(vkbr_t),allocatable :: vkbr(:)
 279 !************************************************************************
 280 
 281  DBG_ENTER("COLL")
 282  call cwtime(cpu_time,wall_time,gflops,"start")
 283 
 284  ! Change FFT mesh if needed
 285  if (ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3))) call wfd_change_ngfft(Wfd,Cryst,Psps,ngfft_gw)
 286 
 287  gw_mgfft = MAXVAL(ngfft_gw(1:3)); gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10)
 288  if (Dtset%pawcross==1) mgfftf = MAXVAL(ngfftf(1:3))
 289 
 290  ! Copy important variables.
 291  comm = Wfd%comm; nsppol = Wfd%nsppol; nspinor = Wfd%nspinor; mband = Wfd%mband
 292  nfft = Wfd%nfft
 293  ABI_CHECK(Wfd%nfftot==nfftot_gw,"Wrong nfftot_gw")
 294  dim_rtwg=1 !; if (nspinor==2) dim_rtwg=2 ! Can reduce size depending on Ep%nI and Ep%nj
 295 
 296  ucrpa_bands(1)=dtset%ucrpa_bands(1)
 297  ucrpa_bands(2)=dtset%ucrpa_bands(2)
 298  luwindow=.false.
 299  if(abs(dtset%ucrpa_window(1)+1_dp)>tol8.and.(abs(dtset%ucrpa_window(2)+1_dp)>tol8)) then
 300    luwindow=.true.
 301  endif
 302 
 303  ! For cRPA calculation of U: read forlb.ovlp
 304  if(dtset%ucrpa>=1) then
 305    call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor,Kmesh,lcor,luwindow,&
 306 & nspinor,nsppol,pawang,dtset%prtvol,ucrpa_bands)
 307  endif
 308 
 309  ks_energy => KS_BSt%eig
 310  qp_energy => QP_BSt%eig; qp_occ => QP_BSt%occ
 311 
 312  chi0_lwing = czero; chi0_uwing= czero; chi0_head = czero
 313 
 314  if (Psps%usepaw==0) then
 315    if (Ep%inclvkb/=0) then
 316      ! Include the term <n,k|[Vnl,iqr]|n"k>' for q->0.
 317      ABI_CHECK(nspinor==1,"nspinor+inclvkb not coded")
 318    else
 319      MSG_WARNING('Neglecting <n,k|[Vnl,iqr]|m,k>')
 320    end if
 321  else
 322    ! For PAW+LDA+U, precalculate <\phi_i|[Hu,r]|phi_j\>
 323    ABI_DT_MALLOC(HUr,(Cryst%natom))
 324    if (Dtset%usepawu/=0) then
 325      call pawhur_init(hur,nsppol,Dtset%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij)
 326    end if
 327  end if
 328 
 329  ! Initialize the completeness correction.
 330  ABI_MALLOC(green_enhigh_w,(Ep%nomega))
 331  green_enhigh_w=czero
 332 
 333  if (Ep%gwcomp==1) then
 334    en_high=MAXVAL(qp_energy(Ep%nbnds,:,:))+Ep%gwencomp
 335    write(msg,'(a,f8.2,a)')' Using completeness correction with energy ',en_high*Ha_eV,' [eV] '
 336    call wrtout(std_out,msg,'COLL')
 337    ABI_MALLOC(wfwfg,(nfft*nspinor**2))
 338 
 339    ! Init the largest G-sphere contained in the FFT box for the wavefunctions.
 340    call gsph_in_fftbox(Gsph_FFT,Cryst,Wfd%ngfft)
 341    call print_gsphere(Gsph_FFT,unit=std_out,prtvol=10)
 342 
 343    ABI_MALLOC(gspfft_igfft,(Gsph_FFT%ng))
 344    ABI_MALLOC(dummy_gbound,(2*gw_mgfft+8,2))
 345 
 346    ! Mapping between G-sphere and FFT box.
 347    call gsph_fft_tabs(Gsph_FFT, [0, 0, 0],Wfd%mgfft,Wfd%ngfft,dummy,dummy_gbound,gspfft_igfft)
 348    ABI_FREE(dummy_gbound)
 349 
 350    if (Psps%usepaw==1) then
 351      ! Prepare the onsite contributions on the GW FFT mesh.
 352      ABI_MALLOC(gw_gfft,(3,nfft))
 353      q0=zero
 354      call get_gftt(ngfft_gw,q0,Cryst%gmet,gw_gsq,gw_gfft) ! The set of plane waves in the FFT Box.
 355      ABI_DT_MALLOC(Pwij_fft,(Psps%ntypat))
 356      call pawpwij_init(Pwij_fft,nfft,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
 357    end if
 358  end if
 359 
 360  ! Setup weight (2 for spin unpolarized systems, 1 for polarized).
 361  ! spin_fact is used to normalize the occupation factors to one.
 362  ! Consider also the AFM case.
 363  SELECT CASE (nsppol)
 364  CASE (1)
 365    weight=two/Kmesh%nbz; spin_fact=half
 366    if (Wfd%nspden==2) then
 367      weight=one/Kmesh%nbz; spin_fact=half
 368    end if
 369    if (nspinor==2) then
 370      weight=one/Kmesh%nbz; spin_fact=one
 371    end if
 372 
 373  CASE (2)
 374    weight=one/Kmesh%nbz; spin_fact=one
 375 
 376  CASE DEFAULT
 377    MSG_BUG("Wrong nsppol")
 378  END SELECT
 379 
 380  ! Weight for points in the IBZ_q
 381  wtk_ltg(:)=1
 382  if (Ep%symchi==1) then
 383    do ik_bz=1,Ltg_q%nbz
 384      wtk_ltg(ik_bz)=0
 385      if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only k in IBZ_q
 386      wtk_ltg(ik_bz)=SUM(Ltg_q%wtksym(:,:,ik_bz))
 387    end do
 388  end if
 389 
 390  write(msg,'(a,i3,a)')' Q-points for long wave-length limit. # ',Ep%nqlwl,ch10
 391  do iqlwl=1,Ep%nqlwl
 392    write(msg_tmp,'(1x,i5,a,2x,3f12.6,a)') iqlwl,')',Ep%qlwl(:,iqlwl),ch10
 393    msg=TRIM(msg)//msg_tmp
 394  end do
 395  call wrtout(std_out,msg,'COLL')
 396 
 397  write(msg,'(a,i2,2a,i2)')&
 398 &  ' Using spectral method for the imaginary part = ',Ep%spmeth,ch10,&
 399 &  ' Using symmetries to sum only over the IBZ_q  = ',Ep%symchi
 400  call wrtout(std_out,msg,'COLL')
 401 
 402  if (use_tr) then
 403    ! Special care has to be taken in metals and/or spin dependent systems
 404    ! as Wfs_val might contain unoccupied states.
 405    call wrtout(std_out,' Using faster algorithm based on time reversal symmetry.')
 406  else
 407    call wrtout(std_out,' Using slow algorithm without time reversal symmetry.')
 408  end if
 409 
 410  ! Evaluate oscillator matrix elements btw partial waves. Note q=Gamma
 411  if (Psps%usepaw==1) then
 412    ABI_DT_MALLOC(Pwij,(Psps%ntypat))
 413    call pawpwij_init(Pwij,Ep%npwepG0,(/zero,zero,zero/),Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
 414 
 415    ABI_DT_MALLOC(Cprj1_bz,(Cryst%natom,nspinor))
 416    call pawcprj_alloc(Cprj1_bz,0,Wfd%nlmn_atm)
 417    ABI_DT_MALLOC(Cprj2_bz,(Cryst%natom,nspinor))
 418    call pawcprj_alloc(Cprj2_bz,0,Wfd%nlmn_atm)
 419 
 420    ABI_DT_MALLOC(Cprj1_ibz,(Cryst%natom,nspinor))
 421    call pawcprj_alloc(Cprj1_ibz,0,Wfd%nlmn_atm)
 422    ABI_DT_MALLOC(Cprj2_ibz,(Cryst%natom,nspinor))
 423    call pawcprj_alloc(Cprj2_ibz,0,Wfd%nlmn_atm)
 424    if (Dtset%pawcross==1) then
 425      ABI_MALLOC(ur_ae1,(nfftf_tot*nspinor))
 426      ABI_MALLOC(ur_ae_onsite1,(nfftf_tot*nspinor))
 427      ABI_MALLOC(ur_ps_onsite1,(nfftf_tot*nspinor))
 428      ABI_MALLOC(ur_ae2,(nfftf_tot*nspinor))
 429      ABI_MALLOC(ur_ae_onsite2,(nfftf_tot*nspinor))
 430      ABI_MALLOC(ur_ps_onsite2,(nfftf_tot*nspinor))
 431      ABI_MALLOC(igfftepsG0f,(Ep%npwepG0))
 432      ABI_MALLOC(tabrf_k,(nfftf_tot))
 433    end if
 434  end if
 435 
 436  ABI_MALLOC(rhotwg,(Ep%npwe*dim_rtwg))
 437  ABI_MALLOC(tabr_k,(nfft))
 438  ABI_MALLOC(ur1_kibz,(nfft*nspinor))
 439  ABI_MALLOC(ur2_kibz,(nfft*nspinor))
 440 
 441  ABI_MALLOC(usr1_k,(nfft*nspinor))
 442  ABI_MALLOC(ur2_k,(nfft*nspinor))
 443  !
 444  ! Tables for the FFT of the oscillators.
 445  !  a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere).
 446  !  b) gw_gbound table for the zero-padded FFT performed in rhotwg.
 447  ABI_MALLOC(igffteps0,(Gsph_epsG0%ng))
 448  ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2))
 449  call gsph_fft_tabs(Gsph_epsG0, [0, 0, 0], gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igffteps0)
 450  if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
 451 #ifdef FC_IBM
 452  ! XLF does not deserve this optimization (problem with [v67mbpt][t03])
 453  use_padfft = 0
 454 #endif
 455  if (use_padfft==0) then
 456    ABI_FREE(gw_gbound)
 457    ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft))
 458  end if
 459  if (Dtset%pawcross==1) then
 460     ABI_MALLOC(gboundf,(2*mgfftf+8,2))
 461    call gsph_fft_tabs(Gsph_epsG0,(/0,0,0/),mgfftf,ngfftf,use_padfftf,gboundf,igfftepsG0f)
 462    if ( ANY(gw_fftalga == (/2,4/)) ) use_padfftf=0
 463    if (use_padfftf==0) then
 464      ABI_FREE(gboundf)
 465      ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf))
 466    end if
 467  end if
 468 
 469  ! TODO this table can be calculated for each k-point
 470  my_nbbpks=0; allup="All"; got=0
 471  ABI_MALLOC(bbp_ks_distrb,(mband,mband,Kmesh%nbz,nsppol))
 472  ABI_MALLOC(bbp_mask,(mband,mband))
 473 
 474  do spin=1,nsppol
 475    do ik_bz=1,Kmesh%nbz
 476 
 477      if (Ep%symchi==1) then
 478        if (Ltg_q%ibzq(ik_bz)/=1) CYCLE  ! Only IBZ_q
 479      end if
 480      ik_ibz=Kmesh%tab(ik_bz)
 481 
 482      call chi0_bbp_mask(Ep,use_tr,QP_BSt,mband,ik_ibz,ik_ibz,spin,spin_fact,bbp_mask)
 483 
 484      call wfd_distribute_bbp(Wfd,ik_ibz,spin,allup,my_nbbp,bbp_ks_distrb(:,:,ik_bz,spin),got,bbp_mask)
 485      my_nbbpks = my_nbbpks + my_nbbp
 486    end do
 487  end do
 488 
 489  ABI_FREE(bbp_mask)
 490 
 491  write(msg,'(a,i0,a)')" Will sum ",my_nbbpks," (b,b',k,s) states in chi0q0."
 492  call wrtout(std_out,msg,'PERS')
 493 
 494  SELECT CASE (Ep%spmeth)
 495  CASE (0)
 496    call wrtout(std_out,' Calculating chi0(q=(0,0,0),omega,G,G")',"COLL")
 497    ABI_MALLOC(green_w,(Ep%nomega))
 498 
 499  CASE (1, 2)
 500    call wrtout(std_out,' Calculating Im chi0(q=(0,0,0),omega,G,G")','COLL')
 501    !
 502    ! === Find max and min resonant transitions for this q, report values for this processor ===
 503    call make_transitions(Wfd,1,Ep%nbnds,nbvw,nsppol,Ep%symchi,Cryst%timrev,GW_TOL_DOCC,&
 504 &    max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,qp_energy,qp_occ,(/zero,zero,zero/),bbp_ks_distrb)
 505 
 506    !FIXME there is a problem in make_transitions due to MPI_enreg
 507    !ltest=(MPI_enreg%gwpara==0.or.MPI_enreg%gwpara==2)
 508    !ABI_CHECK(ltest,"spectral method with gwpara==1 not coded")
 509    !
 510    ! === Calculate frequency dependent weights for Kramers Kronig transform ===
 511    ABI_MALLOC(omegasf,(Ep%nomegasf))
 512    ABI_MALLOC(kkweight,(Ep%nomegasf,Ep%nomega))
 513    !my_wl=1; my_wr=Ep%nomegasf
 514    call setup_spectral(Ep%nomega,Ep%omega,Ep%nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,&
 515 &     0,Ep%zcut,zero,my_wl,my_wr,kkweight)
 516 
 517    if (.not.use_tr) then
 518      MSG_BUG('Hilbert transform requires time-reversal')
 519    end if
 520 
 521    ! allocate heads and wings of the spectral function.
 522    ABI_MALLOC(sf_head,(3,3,my_wl:my_wr))
 523    ABI_MALLOC(sf_lwing,(Ep%npwe,my_wl:my_wr,3))
 524    ABI_MALLOC(sf_uwing,(Ep%npwe,my_wl:my_wr,3))
 525    sf_head=czero; sf_lwing=czero; sf_uwing=czero
 526 
 527    memreq = two*gwpc*Ep%npwe**2*(my_wr-my_wl+1)*b2Gb
 528    write(msg,'(a,f10.4,a)')' memory required per spectral point: ',two*gwpc*Ep%npwe**2*b2Mb,' [Mb]'
 529    call wrtout(std_out,msg,'PERS')
 530    write(msg,'(a,f10.4,a)')' memory required by sf_chi0q0:       ',memreq,' [Gb]'
 531    call wrtout(std_out,msg,'PERS')
 532    if (memreq > two) then
 533      MSG_WARNING(' Memory required for sf_chi0q0 is larger than 2.0 Gb!')
 534    end if
 535    ABI_STAT_MALLOC(sf_chi0,(Ep%npwe,Ep%npwe,my_wl:my_wr), ierr)
 536    ABI_CHECK(ierr==0, 'out-of-memory in sf_chi0q0')
 537    sf_chi0=czero_gw
 538 
 539  CASE DEFAULT
 540    MSG_BUG("Wrong spmeth")
 541  END SELECT
 542 
 543  nkpt_summed=Kmesh%nbz
 544  if (Ep%symchi/=0) then
 545    nkpt_summed=Ltg_q%nibz_ltg
 546    call littlegroup_print(Ltg_q,std_out,Dtset%prtvol,'COLL')
 547  end if
 548 
 549  ABI_DT_MALLOC(vkbr,(Kmesh%nibz))
 550  gradk_not_done=.TRUE.
 551 
 552  write(msg,'(a,i6,a)')' Calculation status ( ',nkpt_summed,' to be completed):'
 553  call wrtout(std_out,msg,'COLL')
 554  !
 555  ! ============================================
 556  ! === Begin big fat loop over transitions ====
 557  ! ============================================
 558  chi0=czero_gw; chi0_sumrule =zero
 559 
 560  ! Loop on spin to calculate $\chi_{\up,\up} + \chi_{\down,\down}$
 561  do spin=1,nsppol
 562    if (ALL(bbp_ks_distrb(:,:,:,spin) /= Wfd%my_rank)) CYCLE
 563 
 564    ! Loop over k-points in the BZ.
 565    do ik_bz=1,Kmesh%nbz
 566      if (Ep%symchi==1) then
 567        if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q
 568      end if
 569 
 570      if (ALL(bbp_ks_distrb(:,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE
 571 
 572      write(msg,'(2(a,i4),a,i2,a,i3)')' ik= ',ik_bz,'/',Kmesh%nbz,' spin=',spin,' done by mpi rank:',Wfd%my_rank
 573      call wrtout(std_out,msg,'PERS')
 574 
 575      ! Get ik_ibz, non-symmorphic phase and symmetries from ik_bz.
 576      call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt)
 577      tabr_k=ktabr(:,ik_bz) ! Table for rotated FFT points
 578      spinrot_kbz(:)=Cryst%spinrot(:,isym_k)
 579      if (Dtset%pawcross==1) tabrf_k(:) = ktabrf(:,ik_bz)
 580 
 581      istwf_k =  Wfd%istwfk(ik_ibz)
 582      npw_k   =  Wfd%npwarr(ik_ibz)
 583      kg_k    => Wfd%Kdata(ik_ibz)%kg_k
 584 
 585      if (Psps%usepaw==0.and.Ep%inclvkb/=0.and.gradk_not_done(ik_ibz)) then
 586        ! Include term <n,k|[Vnl,iqr]|n"k>' for q->0.
 587        call vkbr_init(vkbr(ik_ibz),Cryst,Psps,Ep%inclvkb,istwf_k,npw_k,Kmesh%ibz(:,ik_ibz),kg_k)
 588        gradk_not_done(ik_ibz)=.FALSE.
 589      end if
 590 
 591      ! Loop over "conduction" states.
 592      do band1=1,Ep%nbnds
 593        if (ALL(bbp_ks_distrb(band1,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE
 594 
 595        ug1 => Wfd%Wave(band1,ik_ibz,spin)%ug
 596        call wfd_get_ur(Wfd,band1,ik_ibz,spin,ur1_kibz)
 597 
 598        if (Psps%usepaw==1) then
 599          call wfd_get_cprj(Wfd,band1,ik_ibz,spin,Cryst,Cprj1_ibz,sorted=.FALSE.)
 600          call pawcprj_copy(Cprj1_ibz,Cprj1_bz)
 601          call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_bz)
 602          if (Dtset%pawcross==1) then
 603            call wfd_paw_get_aeur(Wfdf,band1,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae1,ur_ae_onsite1,ur_ps_onsite1)
 604          end if
 605        end if
 606 
 607        do band2=1,Ep%nbnds ! Loop over "valence" states.
 608 !       write(6,*) "ik,band1,band2",ik_bz,band1,band2
 609 
 610 !         -----------------  cRPA for U
 611 !debug         if (.not.luwindow.AND.dtset%ucrpa==1.AND.band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)&
 612 !debug&                                            .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) CYCLE
 613          if (luwindow.AND.dtset%ucrpa==1.AND.((KS_Bst%eig(band1,ik_ibz,spin)-KS_Bst%fermie)<=dtset%ucrpa_window(2)) &
 614 &                                       .AND.((KS_Bst%eig(band1,ik_ibz,spin)-KS_Bst%fermie)>=dtset%ucrpa_window(1)) &
 615 &                                       .AND.((KS_Bst%eig(band2,ik_ibz,spin)-KS_Bst%fermie)<=dtset%ucrpa_window(2)) &
 616 &                                       .AND.((KS_Bst%eig(band2,ik_ibz,spin)-KS_Bst%fermie)>=dtset%ucrpa_window(1))) CYCLE
 617 !         -----------------  cRPA for U
 618 
 619          if (bbp_ks_distrb(band1,band2,ik_bz,spin) /= Wfd%my_rank) CYCLE
 620 
 621          deltaf_b1b2  =spin_fact*(qp_occ(band1,ik_ibz,spin)-qp_occ(band2,ik_ibz,spin))
 622          deltaeKS_b1b2= ks_energy(band1,ik_ibz,spin) - ks_energy(band2,ik_ibz,spin)
 623          deltaeGW_b1b2= qp_energy(band1,ik_ibz,spin) - qp_energy(band2,ik_ibz,spin)
 624 
 625          if (Ep%gwcomp==0) then ! Skip negligible transitions.
 626            if (ABS(deltaf_b1b2) < GW_TOL_DOCC) CYCLE
 627          else
 628            ! when the completeness trick is used, we need to also consider transitions with vanishing deltaf
 629            !Rangel Correction for metals
 630            !if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC) CYCLE
 631            if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC .and. ( ABS(deltaf_b1b2)< GW_TOL_DOCC .or. band1<band2)) CYCLE
 632          end if
 633 
 634          ug2 => Wfd%Wave(band2,ik_ibz,spin)%ug
 635          call wfd_get_ur(Wfd,band2,ik_ibz,spin,ur2_kibz)
 636 
 637          if (Psps%usepaw==1) then
 638            call wfd_get_cprj(Wfd,band2,ik_ibz,spin,Cryst,Cprj2_ibz,sorted=.FALSE.)
 639            call pawcprj_copy(Cprj2_ibz,Cprj2_bz)
 640            call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj2_bz)
 641            if (Dtset%pawcross==1) then
 642              call wfd_paw_get_aeur(Wfdf,band2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae2,ur_ae_onsite2,ur_ps_onsite2)
 643            end if
 644          end if
 645 
 646          SELECT CASE (Ep%spmeth)
 647          CASE (0)
 648            ! Adler-Wiser expression.
 649            ! Add small imaginary of the Time-Ordered resp function but only for non-zero real omega  FIXME What about metals?
 650 
 651            if (.not.use_tr) then
 652              ! Adler-Wiser without time-reversal.
 653              do io=1,Ep%nomega
 654                green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0,one_pole)
 655              end do
 656            else
 657              if (Ep%gwcomp==0) then ! cannot be completely skipped in case of completeness correction
 658                if (band1<band2) CYCLE ! Here we GAIN a factor ~2
 659              end if
 660 
 661              do io=1,Ep%nomega
 662                !Rangel: In metals, the intra-band transitions term does not contain the antiresonant part
 663                !if(abs(deltaeGW_b1b2)>GW_TOL_W0) green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0)
 664                if (band1==band2) green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0,one_pole)
 665                if (band1/=band2) green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0,two_poles)
 666 
 667                if (Ep%gwcomp==1) then ! Calculate the completeness correction
 668                  numerator= -spin_fact*qp_occ(band2,ik_ibz,spin)
 669                  deltaeGW_enhigh_b2=en_high-qp_energy(band2,ik_ibz,spin)
 670                  ! Completeness correction is NOT valid for real frequencies
 671                  if (REAL(Ep%omega(io))<GW_TOL_W0) then
 672                    green_enhigh_w(io) = g0g0w(Ep%omega(io),numerator,deltaeGW_enhigh_b2,Ep%zcut,GW_TOL_W0,two_poles)
 673                  else
 674                    green_enhigh_w(io) = czero_gw
 675                  endif
 676                  !
 677                  !Rangel Correction for metals
 678                  !if (deltaf_b1b2<0.d0) then
 679                  if (band1>=band2 .and. abs(deltaf_b1b2) > GW_TOL_DOCC) then
 680                    green_w(io)= green_w(io) - green_enhigh_w(io)
 681                  else ! Disregard green_w, since it is already accounted for through the time-reversal
 682                    green_w(io)=             - green_enhigh_w(io)
 683                  end if
 684                end if !gwcomp==1
 685              end do !io
 686 
 687              if (Ep%gwcomp==1.and.band1==band2) then
 688                ! Add the "delta part", symmetrization is done inside the routine.
 689                call calc_wfwfg(tabr_k,itim_k,spinrot_kbz,nfft,nspinor,ngfft_gw,ur2_kibz,ur2_kibz,wfwfg)
 690 
 691                if (Psps%usepaw==1) then
 692                  call paw_rho_tw_g(nfft,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,gw_gfft,&
 693 &                  Cprj2_bz,Cprj2_bz,Pwij_fft,wfwfg)
 694 
 695                 ! Add PAW cross term
 696                 if (Dtset%pawcross==1) then
 697                   call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,&
 698 &                  ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,&
 699 &                  ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,&
 700 &                  dim_rtwg,wfwfg)
 701                 end if
 702                end if
 703 
 704                qzero=.TRUE.
 705                call completechi0_deltapart(ik_bz,qzero,Ep%symchi,Ep%npwe,Gsph_FFT%ng,Ep%nomega,nspinor,&
 706 &                nfft,ngfft_gw,gspfft_igfft,Gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0)
 707              end if
 708            end if ! use_tr
 709 
 710          CASE (1, 2)
 711            ! Spectral method, here time-reversal is always assumed.
 712            if (deltaeGW_b1b2<0) CYCLE
 713            call approxdelta(Ep%nomegasf,omegasf,deltaeGW_b1b2,Ep%spsmear,iomegal,iomegar,wl,wr,Ep%spmeth)
 714          END SELECT
 715 
 716          ! FFT of u^*_{b1,k}(r) u_{b2,k}(r) and (q,G=0) limit using small q and k.p perturbation theory
 717          call rho_tw_g(nspinor,Ep%npwe,nfft,ndat1,ngfft_gw,1,use_padfft,igffteps0,gw_gbound,&
 718 &          ur1_kibz,itim_k,tabr_k,ph_mkt,spinrot_kbz,&
 719 &          ur2_kibz,itim_k,tabr_k,ph_mkt,spinrot_kbz,&
 720 &          dim_rtwg,rhotwg)
 721 
 722          if (Psps%usepaw==0) then
 723            ! Matrix elements of i[H,r] for NC pseudopotentials.
 724            rhotwx = nc_ihr_comm(vkbr(ik_ibz),cryst,psps,npw_k,nspinor,istwf_k,Ep%inclvkb,Kmesh%ibz(:,ik_ibz),ug1,ug2,kg_k)
 725 
 726          else
 727            ! 1) Add PAW onsite contribution, projectors are already in the BZ.
 728            call paw_rho_tw_g(Ep%npwe,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,&
 729 &            Cprj1_bz,Cprj2_bz,Pwij,rhotwg)
 730 
 731            ! 2) Matrix elements of i[H,r] for PAW.
 732            rhotwx = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug1,ug2,kg_k,Cprj1_ibz,Cprj2_ibz,HUr)
 733 
 734            ! Add PAW cross term
 735            if (Dtset%pawcross==1) then
 736              call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,&
 737 &             ur_ae1,ur_ae_onsite1,ur_ps_onsite1,itim_k,tabrf_k,ph_mkt,spinrot_kbz,&
 738 &             ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,&
 739 &             dim_rtwg,rhotwg)
 740               ! Add cross-term contribution to the commutator
 741              if (Dtset%userib/=111) then
 742                call paw_cross_ihr_comm(rhotwx,nspinor,nfftf_tot,Cryst,Pawfgrtab,Paw_onsite,&
 743 &                   ur_ae1,ur_ae2,ur_ae_onsite1,ur_ae_onsite2,Cprj1_ibz,Cprj2_ibz)
 744              end if
 745            end if
 746          end if
 747 
 748          ! Treat a possible degeneracy between v and c.
 749          if (ABS(deltaeKS_b1b2)>GW_TOL_W0) then
 750            rhotwx=-rhotwx/deltaeKS_b1b2
 751          else
 752            rhotwx=czero_gw
 753          end if
 754 
 755          SELECT CASE (Ep%spmeth)
 756          CASE (0)
 757 
 758 !          ---------------- Ucrpa (begin)
 759            if(dtset%ucrpa>=1.and..not.luwindow)  then
 760              fac=one
 761              fac1=zero
 762              fac2=zero
 763              fac3=zero
 764              fac4=one
 765              m1=-1
 766              m2=-1
 767              if(dtset%ucrpa<=2) then
 768                call flush_unit(std_out)
 769                if (       band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)&
 770 &                    .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then
 771 !               if(dtset%prtvol>=10)write(6,*)"calculation is in progress",band1,band2,ucrpa_bands(1),ucrpa_bands(2)
 772                  do iat=1, cryst%nattyp(itypatcor)
 773                    do ispinor1=1,nspinor
 774                      do m1=1,2*lcor+1
 775                        fac1=fac1+ real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
 776 &                                 conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)))
 777                        fac2=fac2+ real(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)*&
 778 &                                 conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)))
 779                        do ispinor2=1,nspinor
 780                          do m2=1,2*lcor+1
 781                            fac=fac - real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
 782 &                                    conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* &
 783 &                                    coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)*&
 784 &                                    conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)))
 785                            fac3=fac3 + real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
 786 &                                      conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* &
 787 &                                      coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)*&
 788 &                                      conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)))
 789 !                         if(dtset%prtvol>=10)write(6,*) fac,fac3
 790                          enddo
 791                        enddo
 792 !                       if(dtset%prtvol>=10)write(6,*) fac,fac3,fac1,fac2,fac1*fac2
 793                      enddo
 794                    enddo
 795                  enddo
 796                  fac4=fac
 797 !                 fac=zero
 798                  if(dtset%ucrpa==1) fac=zero
 799 !                 write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0"
 800                endif
 801              else if (dtset%ucrpa==3) then
 802                if        (band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)) then
 803                  do iat=1, cryst%nattyp(itypatcor)
 804                    do ispinor1=1,nspinor
 805                      do m1=1,2*lcor+1
 806                        fac2=fac2-real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
 807 &                                conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)))
 808                      enddo
 809                    enddo
 810                  enddo
 811                  if(dtset%ucrpa==4) fac2=zero
 812                endif
 813                if        (band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then
 814                  do iat=1, cryst%nattyp(itypatcor)
 815                    do ispinor1=1,nspinor
 816                      do m1=1,2*lcor+1
 817                        fac3=fac3-real(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)*&
 818 &                                conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)))
 819                      enddo
 820                    enddo
 821                  enddo
 822                  if(dtset%ucrpa==4) fac3=zero
 823                endif
 824                fac=real(fac2*fac3)
 825              endif
 826 !             if(dtset%prtvol>=10) write(6,'(6i4,e15.5,a)') ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0"
 827 !             if(abs(fac-one)>0.00001) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0"
 828 !             if(dtset%prtvol>=10) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac4," q==0"
 829              green_w=green_w*fac
 830            endif
 831 !          ---------------- Ucrpa (end)
 832 
 833            ! Adler-Wiser expression, to be consistent here we use the KS eigenvalues (?)
 834            call accumulate_chi0_q0(ik_bz,isym_k,itim_k,Ep%gwcomp,nspinor,Ep%npwepG0,Ep,&
 835 &           Cryst,Ltg_q,Gsph_epsG0,chi0,rhotwx,rhotwg,green_w,green_enhigh_w,deltaf_b1b2,chi0_head,chi0_lwing,chi0_uwing)
 836 
 837          CASE (1, 2)
 838            ! Spectral method, to be consistent here we use the KS eigenvalues.
 839            call accumulate_sfchi0_q0(ik_bz,isym_k,itim_k,nspinor,Ep%symchi,Ep%npwepG0,Ep%npwe,Cryst,Ltg_q,&
 840 &            Gsph_epsG0,deltaf_b1b2,my_wl,iomegal,wl,my_wr,iomegar,wr,rhotwx,rhotwg,Ep%nomegasf,sf_chi0,sf_head,sf_lwing,sf_uwing)
 841 
 842          CASE DEFAULT
 843            MSG_BUG("Wrong spmeth")
 844          END SELECT
 845 
 846          ! Accumulating the sum rule on chi0. Eq. (5.284) in G.D. Mahan Many-Particle Physics 3rd edition. [[cite:Mahan2000]]
 847          factor=spin_fact*qp_occ(band2,ik_ibz,spin)
 848 
 849          call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_b1b2,&
 850 &          Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule)
 851 
 852          if (Ep%gwcomp==1) then
 853            ! Include also the completeness correction in the sum rule.
 854            factor=-spin_fact*qp_occ(band2,ik_ibz,spin)
 855            call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_enhigh_b2,&
 856 &            Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule)
 857            if (band1==Ep%nbnds) then
 858              chi0_sumrule(:)=chi0_sumrule(:) + wtk_ltg(ik_bz)*spin_fact*qp_occ(band2,ik_ibz,spin)*deltaeGW_enhigh_b2
 859            end if
 860          end if
 861 
 862        end do !band2
 863      end do !band1
 864 
 865      if (Psps%usepaw==0.and.Ep%inclvkb/=0.and.Ep%symchi==1) then
 866        call vkbr_free(vkbr(ik_ibz)) ! Not need anymore as we loop only over IBZ.
 867      end if
 868    end do !ik_bz
 869  end do !spin
 870 
 871  ABI_FREE(igffteps0)
 872 
 873  call vkbr_free(vkbr)
 874  ABI_DT_FREE(vkbr)
 875 
 876  ! === After big fat loop over transitions, now MPI ===
 877  ! * Master took care of the contribution in case of (metallic|spin) polarized systems.
 878  SELECT CASE (Ep%spmeth)
 879  CASE (0)
 880    ! Adler-Wiser expression
 881    ! Sum contributions from each proc.
 882    ! Looping on frequencies to avoid problems with the size of the MPI packet.
 883    do io=1,Ep%nomega
 884      call xmpi_sum(chi0(:,:,io),comm,ierr)
 885    end do
 886 
 887  CASE (1, 2)
 888    ! Spectral method.
 889    call hilbert_transform(Ep%npwe,Ep%nomega,Ep%nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,Ep%spmeth)
 890 
 891    if (allocated(sf_chi0)) then
 892      ABI_FREE(sf_chi0)
 893    end if
 894 
 895    ! Sum contributions from each proc ===
 896    ! Looping on frequencies to avoid problems with the size of the MPI packet
 897    do io=1,Ep%nomega
 898      call xmpi_sum(chi0(:,:,io),comm,ierr)
 899    end do
 900 
 901    call hilbert_transform_headwings(Ep%npwe,Ep%nomega,Ep%nomegasf,&
 902 &   my_wl,my_wr,kkweight,sf_lwing,sf_uwing,sf_head,chi0_lwing,&
 903 &   chi0_uwing,chi0_head,Ep%spmeth)
 904 
 905  CASE DEFAULT
 906    MSG_BUG("Wrong spmeth")
 907  END SELECT
 908 
 909  ! Divide by the volume
 910 !$OMP PARALLEL WORKSHARE
 911    chi0=chi0*weight/Cryst%ucvol
 912 !$OMP END PARALLEL WORKSHARE
 913 
 914  ! Collect sum rule. pi comes from Im[1/(x-ieta)] = pi delta(x)
 915  call xmpi_sum(chi0_sumrule,comm,ierr)
 916  chi0_sumrule=chi0_sumrule * pi * weight / Cryst%ucvol
 917 
 918  ! Collect heads and wings.
 919  call xmpi_sum(chi0_head,comm,ierr)
 920  call xmpi_sum(chi0_lwing,comm,ierr)
 921  call xmpi_sum(chi0_uwing,comm,ierr)
 922 
 923  chi0_head  = chi0_head * weight/Cryst%ucvol
 924  do io=1,Ep%nomega ! Tensor in the basis of the reciprocal lattice vectors.
 925    chi0_head(:,:,io) = MATMUL(chi0_head(:,:,io),Cryst%gmet) * (two_pi**2)
 926  end do
 927  chi0_lwing = chi0_lwing * weight/Cryst%ucvol
 928  chi0_uwing = chi0_uwing * weight/Cryst%ucvol
 929 
 930  ! ===============================================
 931  ! ==== Symmetrize chi0 in case of AFM system ====
 932  ! ===============================================
 933  ! Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$.
 934  ! Works only in the case of magnetic group Shubnikov type IV.
 935  if (Cryst%use_antiferro) then
 936    call symmetrize_afm_chi0(Cryst,Gsph_epsG0,Ltg_q,Ep%npwe,Ep%nomega,chi0,chi0_head,chi0_lwing,chi0_uwing)
 937  end if
 938 
 939  ! ===================================================
 940  ! ==== Construct heads and wings from the tensor ====
 941  ! ===================================================
 942  do io=1,Ep%nomega
 943    do ig=2,Ep%npwe
 944      wng = chi0_uwing(ig,io,:)
 945      chi0(1,ig,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G")
 946      wng = chi0_lwing(ig,io,:)
 947      chi0(ig,1,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G")
 948    end do
 949    chq = MATMUL(chi0_head(:,:,io), Ep%qlwl(:,1))
 950    chi0(1,1,io) = vdotw(Ep%qlwl(:,1),chq,Cryst%gmet,"G")  ! Use user-defined small q
 951  end do
 952 
 953  ! Impose Hermiticity (valid only for zero or purely imaginary frequencies)
 954  ! MG what about metals, where we have poles around zero?
 955  !  do io=1,Ep%nomega
 956  !    if (ABS(REAL(Ep%omega(io)))<0.00001) then
 957  !      do ig2=1,Ep%npwe
 958  !        do ig1=1,ig2-1
 959  !         chi0(ig2,ig1,io)=GWPC_CONJG(chi0(ig1,ig2,io))
 960  !        end do
 961  !      end do
 962  !    end if
 963  !  end do
 964  !
 965  ! =====================
 966  ! ==== Free memory ====
 967  ! =====================
 968  ABI_FREE(bbp_ks_distrb)
 969  ABI_FREE(rhotwg)
 970  ABI_FREE(tabr_k)
 971  ABI_FREE(ur1_kibz)
 972  ABI_FREE(ur2_kibz)
 973  ABI_FREE(usr1_k)
 974  ABI_FREE(ur2_k)
 975  ABI_FREE(gw_gbound)
 976 
 977  if (Dtset%pawcross==1) then
 978    ABI_FREE(gboundf)
 979  end if
 980 
 981  if (allocated(green_enhigh_w))  then
 982    ABI_FREE(green_enhigh_w)
 983  end if
 984  if (allocated(gw_gfft))  then
 985    ABI_FREE(gw_gfft)
 986  end if
 987  if (allocated(wfwfg))  then
 988    ABI_FREE(wfwfg)
 989  end if
 990  if (allocated(kkweight))  then
 991    ABI_FREE(kkweight)
 992  end if
 993  if (allocated(omegasf))  then
 994    ABI_FREE(omegasf)
 995  end if
 996  if (allocated(green_w))  then
 997    ABI_FREE(green_w)
 998  end if
 999 
1000  if (allocated(sf_head))  then
1001    ABI_FREE(sf_head)
1002  end if
1003  if (allocated(sf_lwing))  then
1004    ABI_FREE(sf_lwing)
1005  end if
1006  if (allocated(sf_uwing))  then
1007    ABI_FREE(sf_uwing)
1008  end if
1009  if (allocated(gspfft_igfft))  then
1010    ABI_FREE(gspfft_igfft)
1011  end if
1012 
1013  call gsph_free(Gsph_FFT)
1014 
1015  if (Psps%usepaw==1) then ! deallocation for PAW.
1016    call pawcprj_free(Cprj1_bz )
1017    ABI_DT_FREE(Cprj1_bz)
1018    call pawcprj_free(Cprj2_bz )
1019    ABI_DT_FREE(Cprj2_bz)
1020    call pawcprj_free(Cprj1_ibz )
1021    ABI_DT_FREE(Cprj1_ibz)
1022    call pawcprj_free(Cprj2_ibz )
1023    ABI_DT_FREE(Cprj2_ibz)
1024    call pawpwij_free(Pwij)
1025    ABI_DT_FREE(Pwij)
1026    if (allocated(Pwij_fft)) then
1027      call pawpwij_free(Pwij_fft)
1028      ABI_DT_FREE(Pwij_fft)
1029    end if
1030    call pawhur_free(Hur)
1031    ABI_DT_FREE(Hur)
1032    if (Dtset%pawcross==1) then
1033      ABI_FREE(ur_ae1)
1034      ABI_FREE(ur_ae_onsite1)
1035      ABI_FREE(ur_ps_onsite1)
1036      ABI_FREE(ur_ae2)
1037      ABI_FREE(ur_ae_onsite2)
1038      ABI_FREE(ur_ps_onsite2)
1039      ABI_FREE(tabrf_k)
1040      ABI_FREE(gboundf)
1041      ABI_FREE(igfftepsG0f)
1042    end if
1043  end if
1044 
1045  if(dtset%ucrpa>=1) then
1046    ABI_DEALLOCATE(coeffW_BZ)
1047  endif
1048 
1049  call cwtime(cpu_time,wall_time,gflops,"stop")
1050  write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time
1051 
1052  DBG_EXIT("COLL")
1053 
1054 end subroutine cchi0q0

ABINIT/chi0q0_intraband [ Functions ]

[ Top ] [ Functions ]

NAME

 chi0q0_intraband

FUNCTION

 Calculate chi0 in the limit q-->0

INPUTS

  use_tr=If .TRUE. Wfs_val are allocate and only resonant transitions are evaluated (assumes time reversal symmetry)
  Ep= datatype gathering differening parameters related to the calculation of the inverse dielectric matrix
  Gsph_epsG0<gvectors_data_type>: Info on the G-sphere used to describe chi0/espilon (including umklapp)
    %ng=number of G vectors
    %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G  in the array gvec
    %phmGt(ng,nsym)=phase factor e^{-iG.\tau} needed to symmetrize oscillator matrix elements and chi0
    %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
    %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
  Ep%nbnds=number of bands
  ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths.
  Ep%nomega=number of frequencies
  Cryst<crystal_t>= data type gathering info on symmetries and unit cell
   %natom=number of atoms
   %nsym=number of symmetry operations
   %symrec(3,3,nsym)=symmetry operations in reciprocal space
   %typat(natom)=type of each atom
   %xred(3,natom)=reduced coordinated of atoms
   %rprimd(3,3)=dimensional primitive translations in real space (bohr)
   %timrev=2 if time-reversal symmetry can be used, 1 otherwise
  Ep%npwe=number of planewaves for sigma exchange (input variable)
  Ep%nsppol=1 for unpolarized, 2 for spin-polarized
  Ep%omega(Ep%nomega)=frequencies
  Psps <type(pseudopotential_type)>=variables related to pseudopotentials
     %mpsang=1+maximum angular momentum for nonlocal pseudopotential
  Pawang<pawang_type> angular mesh discretization and related data:
  Pawrad(ntypat*usepaw)<Pawrad_type>=paw radial mesh and related data
  Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels
  BSt<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities)
    %mband=MAX number of bands over k-points and spin (==Ep%nbnds)
    %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band
    %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes
  Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.

OUTPUT

  chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix for wavevector qq,
   and frequencies defined by Ep%omega

NOTES

  *) The terms "head", "wings" and "body" of chi(G,Gp) refer to
     G=Gp=0, either G or Gp=0, and neither=0 respectively

TODO

  Check npwepG0 before Switching on umklapp

PARENTS

      screening

CHILDREN

      assemblychi0_sym,get_bz_item,getnel,gsph_fft_tabs,kmesh_free,kmesh_init
      littlegroup_free,littlegroup_init,littlegroup_print,pack_eneocc
      paw_rho_tw_g,paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free
      pawhur_free,pawhur_init,pawpwij_free,pawpwij_init,print_arr,rho_tw_g
      rotate_fft_mesh,symmetrize_afm_chi0,unpack_eneocc,vkbr_free,vkbr_init
      wfd_change_ngfft,wfd_distribute_bands,wfd_get_cprj,wfd_get_ur,wrtout
      xmpi_sum

SOURCE

1991 subroutine chi0q0_intraband(Wfd,Cryst,Ep,Psps,BSt,Gsph_epsG0,Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,use_tr,usepawu,&
1992 &  ngfft_gw,chi0,chi0_head,chi0_lwing,chi0_uwing)
1993 
1994 
1995 !This section has been created automatically by the script Abilint (TD).
1996 !Do not modify the following lines by hand.
1997 #undef ABI_FUNC
1998 #define ABI_FUNC 'chi0q0_intraband'
1999 !End of the abilint section
2000 
2001  implicit none
2002 
2003 !Arguments ------------------------------------
2004 !scalars
2005  integer,intent(in) :: usepawu
2006  logical,intent(in) :: use_tr
2007  type(ebands_t),intent(in) :: BSt
2008  type(crystal_t),intent(in) :: Cryst
2009  type(em1params_t),intent(in) :: Ep
2010  type(gsphere_t),intent(in) :: Gsph_epsG0
2011  type(Pseudopotential_type),intent(in) :: Psps
2012  type(Pawang_type),intent(in) :: Pawang
2013  type(wfd_t),target,intent(inout) :: Wfd
2014 !arrays
2015  integer,intent(in) :: ngfft_gw(18)
2016  complex(gwpc),intent(out) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega)
2017  complex(dpc),intent(out) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)
2018  complex(dpc),intent(out) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)
2019  complex(dpc),intent(out) :: chi0_head(3,3,Ep%nomega)
2020  type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw)
2021  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
2022  type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw)
2023  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
2024 
2025 !Local variables ------------------------------
2026 !scalars
2027  integer,parameter :: tim_fourdp1=1,two_poles=2,one_pole=1,ndat1=1
2028  integer,parameter :: unitdos0=0,option1=1,NOMEGA_PRINTED=15
2029  integer :: nqlwl,nband_k,iomega,istwf_k,npw_k,my_nband,lbidx
2030  integer :: band,itim_k,ik_bz,ik_ibz,io,isym_k,spin,iqlwl!,gw_eet !ig,ig1,ig2,my_nbbp,my_nbbpks
2031  integer :: nkpt_summed,dim_rtwg,use_padfft,gw_fftalga,ifft
2032  integer :: kptopt,isym,nsppol,nspinor
2033  integer :: comm,ierr,gw_mgfft,use_umklp,inclvkb
2034  real(dp) :: spin_fact,deltaf_b1b2,weight
2035  real(dp) :: deltaeGW_b1b2,zcut
2036  real(dp),parameter :: dummy_dosdeltae=HUGE(zero)
2037  real(dp) :: o_entropy,o_nelect,maxocc
2038  complex(dpc) :: ph_mkt
2039  logical :: iscompatibleFFT !,ltest
2040  character(len=500) :: msg,msg_tmp !,allup
2041  type(kmesh_t) :: Kmesh
2042  type(littlegroup_t) :: Ltg_q
2043  type(vkbr_t) :: vkbr
2044 !arrays
2045  integer :: my_band_list(Wfd%mband)
2046  integer,ABI_CONTIGUOUS pointer :: kg_k(:,:)
2047  integer,allocatable :: ktabr(:,:),irottb(:,:)
2048  !integer :: got(Wfd%nproc)
2049  integer,allocatable :: tabr_k(:),igffteps0(:),gw_gbound(:,:)
2050  real(dp),parameter :: q0(3)=(/zero,zero,zero/)
2051  real(dp) :: kpt(3),dedk(3),kbz(3),spinrot_kbz(4)
2052  !real(dp),ABI_CONTIGUOUS pointer :: ks_energy(:,:,:),qp_energy(:,:,:),qp_occ(:,:,:)
2053  real(dp) :: shift_ene(BSt%mband,BSt%nkpt,BSt%nsppol)
2054  real(dp) :: delta_occ(BSt%mband,BSt%nkpt,BSt%nsppol)
2055  !real(dp) :: eigen_vec(BSt%bantot)
2056  real(dp) :: o_doccde(BSt%bantot)
2057  real(dp) :: eigen_pdelta_vec(BSt%bantot),eigen_mdelta_vec(BSt%bantot)
2058  real(dp) :: o_occ_pdelta(BSt%bantot),o_occ_mdelta(BSt%bantot)
2059  real(dp) :: delta_ene(BSt%mband,BSt%nkpt,BSt%nsppol)
2060  real(dp) :: test_docc(BSt%mband,BSt%nkpt,BSt%nsppol)
2061  real(dp),allocatable :: qlwl(:,:)
2062  complex(gwpc) :: comm_kbbs(3,Wfd%nspinor**2)
2063  complex(dpc),allocatable :: ihr_comm(:,:,:,:,:)
2064  complex(gwpc),allocatable :: rhotwg(:)
2065  complex(dpc) :: green_w(Ep%nomega)
2066  complex(gwpc),allocatable :: ur1(:)
2067  complex(gwpc),ABI_CONTIGUOUS pointer :: ug(:)
2068  logical :: bmask(Wfd%mband)
2069  type(pawcprj_type),allocatable :: Cprj1_bz(:,:),Cprj1_ibz(:,:),Cp_bks(:,:)
2070  type(pawpwij_t),allocatable :: Pwij(:)
2071  type(pawhur_t),allocatable :: Hur(:)
2072 
2073 !************************************************************************
2074 
2075  DBG_ENTER("COLL")
2076 
2077  nsppol  = Wfd%nsppol
2078  nspinor = Wfd%nspinor
2079 
2080  gw_mgfft = MAXVAL(ngfft_gw(1:3))
2081  gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10)
2082 
2083  ! Calculate <k,b1|i[H,r]|k',b2>.
2084  inclvkb=2; if (Wfd%usepaw==1) inclvkb=0
2085  ABI_MALLOC(ihr_comm,(3,nspinor**2,Wfd%mband,Wfd%nkibz,nsppol))
2086  ihr_comm = czero
2087 
2088  if (Wfd%usepaw==1) then
2089    ABI_DT_MALLOC(Cp_bks,(Cryst%natom,nspinor))
2090    call pawcprj_alloc(Cp_bks,0,Wfd%nlmn_atm)
2091    ABI_DT_MALLOC(HUr,(Cryst%natom))
2092    if (usepawu/=0) then ! For PAW+LDA+U, precalculate <\phi_i|[Hu,r]|phi_j\>.
2093      call pawhur_init(hur,nsppol,Wfd%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij)
2094    end if
2095  end if
2096 
2097  do spin=1,nsppol
2098    do ik_ibz=1,Wfd%nkibz
2099      npw_k  =  Wfd%npwarr(ik_ibz)
2100      nband_k=  Wfd%nband(ik_ibz,spin)
2101      kpt    =  Wfd%kibz(:,ik_ibz)
2102      kg_k   => Wfd%Kdata(ik_ibz)%kg_k
2103      istwf_k = Wfd%istwfk(ik_ibz)
2104      ABI_CHECK(istwf_k==1,"istwf_k/=1 not coded")
2105 
2106      ! Distribute bands.
2107      bmask=.FALSE.; bmask(1:nband_k)=.TRUE. ! TODO only bands around EF should be included.
2108      call wfd_distribute_bands(Wfd,ik_ibz,spin,my_nband,my_band_list,bmask=bmask)
2109      if (my_nband==0) CYCLE
2110 
2111      if (Wfd%usepaw==0.and.inclvkb/=0) then ! Include term <n,k|[Vnl,iqr]|n"k>' for q->0.
2112        call vkbr_init(vkbr,Cryst,Psps,inclvkb,istwf_k,npw_k,kpt,kg_k)
2113      end if
2114 
2115      do lbidx=1,my_nband
2116        band=my_band_list(lbidx)
2117        ug => Wfd%Wave(band,ik_ibz,spin)%ug
2118 
2119        if (Wfd%usepaw==0) then
2120          ! Matrix elements of i[H,r] for NC pseudopotentials.
2121          comm_kbbs = nc_ihr_comm(vkbr,cryst,psps,npw_k,nspinor,istwf_k,inclvkb,Kmesh%ibz(:,ik_ibz),ug,ug,kg_k)
2122        else
2123          ! Matrix elements of i[H,r] for PAW.
2124          call wfd_get_cprj(Wfd,band,ik_ibz,spin,Cryst,Cp_bks,sorted=.FALSE.)
2125          comm_kbbs = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug,ug,kg_k,Cp_bks,Cp_bks,HUr)
2126        end if
2127 
2128        ihr_comm(:,:,band,ik_ibz,spin) = comm_kbbs
2129      end do
2130 
2131      call vkbr_free(vkbr) ! Not need anymore as we loop only over IBZ.
2132    end do
2133  end do
2134  !
2135  ! Gather the commutator on each node.
2136  call xmpi_sum(ihr_comm,Wfd%comm,ierr)
2137 
2138  if (Wfd%usepaw==1) then
2139    call pawcprj_free(Cp_bks)
2140    ABI_DT_FREE(Cp_bks)
2141    call pawhur_free(Hur)
2142    ABI_DT_FREE(Hur)
2143  end if
2144 
2145  nqlwl=1
2146  ABI_MALLOC(qlwl,(3,nqlwl))
2147  !qlwl = GW_Q0_DEFAULT(3)
2148  qlwl(:,1) = (/0.00001_dp, 0.00002_dp, 0.00003_dp/)
2149  !
2150  write(msg,'(a,i3,a)')' Q-points for long wave-length limit in chi0q_intraband. # ',nqlwl,ch10
2151  do iqlwl=1,nqlwl
2152    write(msg_tmp,'(1x,i5,a,2x,3f12.6,a)') iqlwl,')',qlwl(:,iqlwl),ch10
2153    msg=TRIM(msg)//msg_tmp
2154  end do
2155  call wrtout(std_out,msg,'COLL')
2156  !
2157  ! delta_ene =  e_{b,k-q} - e_{b,k} = -q. <b,k| i[H,r] |b,k> + O(q^2).
2158  delta_ene = zero
2159  do spin=1,nsppol
2160    do ik_ibz=1,Wfd%nkibz
2161      do band=1,Wfd%nband(ik_ibz,spin)
2162        dedk = REAL(ihr_comm(:,1,band,ik_ibz,spin))
2163        delta_ene(band,ik_ibz,spin) = -vdotw(qlwl(:,1),dedk,Cryst%gmet,"G")
2164      end do
2165    end do
2166  end do
2167 
2168  maxocc=two/(nsppol*nspinor)
2169 
2170  ! Calculate the occupations at f(e+delta/2).
2171  shift_ene = BSt%eig + half*delta_ene
2172 
2173  call pack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,BSt%bantot,shift_ene,eigen_pdelta_vec)
2174 
2175  call getnel(o_doccde,dummy_dosdeltae,eigen_pdelta_vec,o_entropy,BSt%fermie,maxocc,BSt%mband,BSt%nband,&
2176 &  o_nelect,BSt%nkpt,BSt%nsppol,o_occ_pdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk)
2177  write(std_out,*)"nelect1: ",o_nelect
2178  !
2179  ! Calculate the occupations at f(e-delta/2).
2180  shift_ene = BSt%eig - half*delta_ene
2181 
2182  call pack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,BSt%bantot,shift_ene,eigen_mdelta_vec)
2183 
2184  call getnel(o_doccde,dummy_dosdeltae,eigen_mdelta_vec,o_entropy,BSt%fermie,maxocc,BSt%mband,BSt%nband,&
2185 &  o_nelect,BSt%nkpt,BSt%nsppol,o_occ_mdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk)
2186  write(std_out,*)"nelect2: ",o_nelect
2187  !
2188  ! f(e-delta/2) - f(e+delta/2).
2189  o_occ_pdelta = o_occ_mdelta - o_occ_pdelta
2190 
2191  call unpack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,o_occ_pdelta,delta_occ)
2192  !
2193  ! Expand f(e-delta/2) - f(e+delta/2) up to the first order in the small q.
2194  do spin=1,nsppol
2195    do ik_ibz=1,Wfd%nkibz
2196      do band=1,Wfd%nband(ik_ibz,spin)
2197        dedk = REAL(ihr_comm(:,1,band,ik_ibz,spin))
2198        test_docc(band,ik_ibz,spin) = +vdotw(qlwl(:,1),dedk,Cryst%gmet,"G") * BSt%doccde(band,ik_ibz,spin)
2199        write(std_out,'(a,3(i0,1x),1x,3es16.8)')" spin,ik_ibz,band, delta_occ: ",&
2200 &      spin,ik_ibz,band,delta_occ(band,ik_ibz,spin),&
2201 &      test_docc(band,ik_ibz,spin),delta_occ(band,ik_ibz,spin)-test_docc(band,ik_ibz,spin)
2202      end do
2203    end do
2204  end do
2205 
2206 ! MSG_ERROR("DONE")
2207 ! do spin=1,nsppol
2208 !   do ik_ibz=1,Wfd%nkibz
2209 !     nband_k = Wfd%nband(ik_ibz,spin)
2210 !     do band=1,nband_k
2211 !       write(std_out,'(a,3i3,2es14.6)')" spin, band, ik_ibz, delta_ene, delta_occ ",&
2212 !&        spin,band,ik_ibz,delta_ene(band,ik_ibz,spin),delta_occ(band,ik_ibz,spin)
2213 !     end do
2214 !   end do
2215 ! end do
2216 
2217  ABI_FREE(ihr_comm)
2218  ABI_FREE(qlwl)
2219 
2220  if ( ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,ngfft_gw)
2221 
2222  ! TODO take into account the case of random k-meshes.
2223  kptopt=3
2224  call kmesh_init(Kmesh,Cryst,Wfd%nkibz,Wfd%kibz,kptopt)
2225  !
2226  !=== Get the FFT index of $ (R^{-1}(r-\tau)) $ ===
2227  !* S= $\transpose R^{-1}$ and k_BZ = S k_IBZ
2228  !* irottb is the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Sk.
2229  ABI_MALLOC(irottb,(Wfd%nfftot,Cryst%nsym))
2230 
2231  call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,irottb,iscompatibleFFT)
2232  ABI_CHECK(iscompatibleFFT,"FFT mesh not compatible with symmetries")
2233 
2234  ABI_MALLOC(ktabr,(Wfd%nfftot,Kmesh%nbz))
2235  do ik_bz=1,Kmesh%nbz
2236    isym=Kmesh%tabo(ik_bz)
2237    do ifft=1,Wfd%nfftot
2238      ktabr(ifft,ik_bz)=irottb(ifft,isym)
2239    end do
2240  end do
2241  ABI_FREE(irottb)
2242  !
2243  ! === Setup weight (2 for spin unpolarized systems, 1 for polarized) ===
2244  ! spin_fact is used to normalize the occupation factors to one.
2245  ! Consider also the AFM case.
2246  SELECT CASE (nsppol)
2247  CASE (1)
2248    weight=two/Kmesh%nbz; spin_fact=half
2249    if (Wfd%nspden==2) then
2250      weight=one/Kmesh%nbz; spin_fact=half
2251    end if
2252    if (nspinor==2) then
2253      weight=one/Kmesh%nbz; spin_fact=one
2254    end if
2255  CASE (2)
2256    weight=one/Kmesh%nbz; spin_fact=one
2257  CASE DEFAULT
2258    MSG_BUG("Wrong nsppol")
2259  END SELECT
2260 
2261  use_umklp=0
2262  call littlegroup_init(q0,Kmesh,Cryst,use_umklp,Ltg_q,Ep%npwepG0,gvec=Gsph_epsG0%gvec)
2263 
2264  write(msg,'(a,i2)')' Using symmetries to sum only over the IBZ_q  = ',Ep%symchi
2265  call wrtout(std_out,msg,'COLL')
2266  !
2267  ! === Evaluate oscillator matrix elements btw partial waves. Note that q=Gamma is used.
2268  if (Psps%usepaw==1) then
2269    ABI_DT_MALLOC(Pwij,(Psps%ntypat))
2270    call pawpwij_init(Pwij,Ep%npwepG0, [zero,zero,zero], Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
2271 
2272    ABI_DT_MALLOC(Cprj1_bz ,(Cryst%natom,nspinor))
2273    call pawcprj_alloc(Cprj1_bz, 0,Wfd%nlmn_atm)
2274    ABI_DT_MALLOC(Cprj1_ibz,(Cryst%natom,nspinor))
2275    call pawcprj_alloc(Cprj1_ibz,0,Wfd%nlmn_atm)
2276  end if
2277 
2278  ABI_MALLOC(rhotwg,(Ep%npwe*nspinor**2))
2279  ABI_MALLOC(tabr_k,(Wfd%nfftot))
2280  ABI_MALLOC(ur1,(Wfd%nfft*nspinor))
2281  !
2282  ! Tables for the FFT of the oscillators.
2283  !  a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere).
2284  !  b) gw_gbound table for the zero-padded FFT performed in rhotwg.
2285  ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2))
2286  ABI_MALLOC(igffteps0,(Gsph_epsG0%ng))
2287 
2288  call gsph_fft_tabs(Gsph_epsG0, [0, 0, 0], gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igffteps0)
2289  if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
2290  if (use_padfft==0) then
2291    ABI_FREE(gw_gbound)
2292    ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft))
2293  end if
2294 
2295  nkpt_summed=Kmesh%nbz
2296  if (Ep%symchi/=0) then
2297    nkpt_summed=Ltg_q%nibz_ltg
2298    call littlegroup_print(Ltg_q,std_out,Wfd%prtvol,'COLL')
2299  end if
2300  !
2301  ! ============================================
2302  ! === Begin big fat loop over transitions ====
2303  ! ============================================
2304  chi0 = czero_gw
2305  chi0_head = czero_gw; chi0_lwing = czero_gw; chi0_uwing = czero_gw
2306  dim_rtwg=1; if (nspinor==2) dim_rtwg=2 !can reduce size depending on Ep%nI and Ep%nj
2307 
2308  zcut = Ep%zcut
2309  zcut = 0.1/Ha_eV
2310  write(std_out,*)" using zcut ",zcut*Ha_eV," [eV]"
2311 
2312  ! Loop on spin to calculate $ \chi_{\up,\up} + \chi_{\down,\down}
2313  do spin=1,nsppol
2314    ! Loop over k-points in the BZ.
2315    do ik_bz=1,Kmesh%nbz
2316      if (Ep%symchi==1) then
2317        if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q
2318      end if
2319 
2320      ! Get ik_ibz, non-symmorphic phase and symmetries from ik_bz.
2321      call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt)
2322      tabr_k=ktabr(:,ik_bz) ! Table for rotated FFT points
2323      spinrot_kbz(:)=Cryst%spinrot(:,isym_k)
2324      nband_k=Wfd%nband(ik_ibz,spin)
2325 
2326      ! Distribute bands.
2327      bmask=.FALSE.; bmask(1:nband_k)=.TRUE. ! TODO only bands around EF should be included.
2328      call wfd_distribute_bands(Wfd,ik_ibz,spin,my_nband,my_band_list,bmask=bmask)
2329      if (my_nband==0) CYCLE
2330 
2331      write(msg,'(2(a,i4),a,i2,a,i3)')' ik = ',ik_bz,' / ',Kmesh%nbz,' spin = ',spin,' done by processor ',Wfd%my_rank
2332      call wrtout(std_out,msg,'PERS')
2333 
2334      do lbidx=1,my_nband
2335        ! Loop over bands treated by this node.
2336        band=my_band_list(lbidx)
2337        call wfd_get_ur(Wfd,band,ik_ibz,spin,ur1)
2338 
2339        if (Psps%usepaw==1) then
2340          call wfd_get_cprj(Wfd,band,ik_ibz,spin,Cryst,Cprj1_ibz,sorted=.FALSE.)
2341          call pawcprj_copy(Cprj1_ibz,Cprj1_bz)
2342          call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_bz)
2343        end if
2344 
2345        deltaf_b1b2  = spin_fact*delta_occ(band,ik_ibz,spin)
2346        deltaeGW_b1b2= delta_ene(band,ik_ibz,spin)
2347 
2348        ! Add small imaginary of the Time-Ordered resp function but only for non-zero real omega  FIXME What about metals?
2349        if (.not.use_tr) then
2350          do io=1,Ep%nomega
2351            !green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,-one,one_pole)
2352            green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,GW_TOL_W0,one_pole)
2353          end do
2354        else
2355          do io=1,Ep%nomega ! This expression implements time-reversal even when the input k-mesh breaks it.
2356            !green_w(io) = half * g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,-one,two_poles)
2357            green_w(io) = half * g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,GW_TOL_W0,two_poles)
2358          end do !io
2359        end if ! use_tr
2360 
2361        ! FFT of u^*_{b1,k}(r) u_{b2,k}(r).
2362        call rho_tw_g(nspinor,Ep%npwe,Wfd%nfft,ndat1,ngfft_gw,1,use_padfft,igffteps0,gw_gbound,&
2363 &        ur1,itim_k,tabr_k,ph_mkt,spinrot_kbz,&
2364 &        ur1,itim_k,tabr_k,ph_mkt,spinrot_kbz,&
2365 &        dim_rtwg,rhotwg)
2366 
2367        if (Psps%usepaw==1) then
2368          ! Add PAW onsite contribution, projectors are already in the BZ.
2369          call paw_rho_tw_g(Ep%npwe,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,&
2370 &          Cprj1_bz,Cprj1_bz,Pwij,rhotwg)
2371        end if
2372 
2373        ! ==== Adler-Wiser expression, to be consistent here we use the KS eigenvalues (?) ====
2374        call assemblychi0_sym(ik_bz,nspinor,Ep,Ltg_q,green_w,Ep%npwepG0,rhotwg,Gsph_epsG0,chi0)
2375      end do !band
2376    end do !ik_bz
2377  end do !spin
2378 
2379  ! Collect body, heads and wings within comm
2380  comm=Wfd%comm
2381  do io=1,Ep%nomega
2382    call xmpi_sum(chi0(:,:,io),comm,ierr)
2383  end do
2384  call xmpi_sum(chi0_head,comm,ierr)
2385  call xmpi_sum(chi0_lwing,comm,ierr)
2386  call xmpi_sum(chi0_uwing,comm,ierr)
2387 
2388  ! Divide by the volume
2389  chi0       = chi0       * weight/Cryst%ucvol
2390  chi0_head  = chi0_head  * weight/Cryst%ucvol
2391  do io=1,Ep%nomega ! Tensor in the basis of the reciprocal lattice vectors.
2392    chi0_head(:,:,io) = MATMUL(chi0_head(:,:,io),Cryst%gmet) * (two_pi**2)
2393  end do
2394  chi0_lwing = chi0_lwing * weight/Cryst%ucvol
2395  chi0_uwing = chi0_uwing * weight/Cryst%ucvol
2396  !
2397  ! ===============================================
2398  ! ==== Symmetrize chi0 in case of AFM system ====
2399  ! ===============================================
2400  ! * Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$.
2401  ! * Works only in the case of magnetic group Shubnikov type IV.
2402  if (Cryst%use_antiferro) then
2403    call symmetrize_afm_chi0(Cryst,Gsph_epsG0,Ltg_q,Ep%npwe,Ep%nomega,chi0,chi0_head,chi0_lwing,chi0_uwing)
2404  end if
2405  !
2406  ! ===================================================
2407  ! ==== Construct heads and wings from the tensor ====
2408  ! ===================================================
2409  !do io=1,Ep%nomega
2410  !  do ig=2,Ep%npwe
2411  !    wng = chi0_uwing(ig,io,:)
2412  !    chi0(1,ig,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G")
2413  !    wng = chi0_lwing(ig,io,:)
2414  !    chi0(ig,1,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G")
2415  !  end do
2416  !  chq = MATMUL(chi0_head(:,:,io), Ep%qlwl(:,1))
2417  !  chi0(1,1,io) = vdotw(Ep%qlwl(:,1),chq,Cryst%gmet,"G")  ! Use user-defined small q
2418  !end do
2419  !call wfd_barrier(Wfd)
2420 
2421  ! Impose Hermiticity (valid only for zero or purely imaginary frequencies)
2422  ! MG what about metals, where we have poles around zero?
2423  !if (dtset%gw_eet/=-1) then
2424  !  do io=1,Ep%nomega
2425  !    if (ABS(REAL(Ep%omega(io)))<0.00001) then
2426  !      do ig2=1,Ep%npwe
2427  !        do ig1=1,ig2-1
2428  !         chi0(ig2,ig1,io)=CONJG(chi0(ig1,ig2,io))
2429  !        end do
2430  !      end do
2431  !    end if
2432  !  end do
2433  !end if
2434 
2435  do iomega=1,MIN(Ep%nomega,NOMEGA_PRINTED)
2436    write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0_intra(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]'
2437    call wrtout(std_out,msg,'COLL')
2438    call print_arr(chi0(:,:,iomega),unit=std_out)
2439  end do
2440 
2441  ! =====================
2442  ! ==== Free memory ====
2443  ! =====================
2444  ABI_FREE(rhotwg)
2445  ABI_FREE(tabr_k)
2446  ABI_FREE(ur1)
2447  ABI_FREE(gw_gbound)
2448  ABI_FREE(ktabr)
2449  ABI_FREE(igffteps0)
2450 
2451  ! deallocation for PAW.
2452  if (Psps%usepaw==1) then
2453    call pawcprj_free(Cprj1_bz)
2454    ABI_DT_FREE(Cprj1_bz)
2455    call pawcprj_free(Cprj1_ibz)
2456    ABI_DT_FREE(Cprj1_ibz)
2457    call pawpwij_free(Pwij)
2458    ABI_DT_FREE(Pwij)
2459  end if
2460 
2461  call littlegroup_free(Ltg_q)
2462  call kmesh_free(Kmesh)
2463 
2464  DBG_EXIT("COLL")
2465 
2466 end subroutine chi0q0_intraband

ABINIT/m_chi0 [ Modules ]

[ Top ] [ Modules ]

NAME

  m_chi0

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_chi0
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_abicore
33  use m_xmpi
34  use m_errors
35  use m_hide_blas
36  use m_time
37  use m_wfd
38 
39  use m_gwdefs,          only : GW_TOL_DOCC, GW_TOL_W0, czero_gw, em1params_t, g0g0w
40  use m_numeric_tools,   only : imin_loc, print_arr
41  use m_geometry,        only : normv, vdotw
42  use m_crystal,         only : crystal_t
43  use m_fft_mesh,        only : rotate_FFT_mesh, get_gftt
44  use m_occ,             only : getnel
45  use m_ebands,          only : pack_eneocc, unpack_eneocc
46  use m_bz_mesh,         only : kmesh_t, kmesh_init, kmesh_free, get_BZ_item, get_BZ_diff, &
47 &                              littlegroup_t, littlegroup_print, littlegroup_free, littlegroup_init
48  use m_gsphere,         only : gsphere_t, gsph_fft_tabs, gsph_in_fftbox, gsph_free, print_gsphere
49  use m_io_tools,        only : flush_unit
50  use m_oscillators,     only : rho_tw_g, calc_wfwfg
51  use m_vkbr,            only : vkbr_t, vkbr_free, vkbr_init, nc_ihr_comm
52  use m_chi0tk,          only : hilbert_transform, setup_spectral, assemblychi0_sym, assemblychi0sf, symmetrize_afm_chi0, &
53                                approxdelta, completechi0_deltapart, accumulate_chi0sumrule, make_transitions, &
54                                chi0_bbp_mask, accumulate_chi0_q0, accumulate_sfchi0_q0, hilbert_transform_headwings
55  use m_pawang,          only : pawang_type
56  use m_pawrad,          only : pawrad_type
57  use m_pawtab,          only : pawtab_type
58  use m_paw_ij,          only : paw_ij_type
59  use m_pawfgrtab,       only : pawfgrtab_type
60  use m_pawcprj,         only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
61  use m_pawpwij,         only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
62  use m_paw_sym,         only : paw_symcprj
63  use m_paw_pwaves_lmn,  only : paw_pwaves_lmn_t
64  use m_paw_hr,          only : pawhur_t, pawhur_free, pawhur_init, paw_ihr, paw_cross_ihr_comm
65  use m_read_plowannier, only : read_plowannier
66 
67  implicit none
68 
69  private