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_ebands<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<wfdgw_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

SOURCE

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

SOURCE

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

SOURCE

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

ABINIT/m_chi0 [ Modules ]

[ Top ] [ Modules ]

NAME

  m_chi0

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2022 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 .

SOURCE

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