TABLE OF CONTENTS


ABINIT/m_bethe_salpeter [ Modules ]

[ Top ] [ Modules ]

NAME

  m_bethe_salpeter

FUNCTION

  Main routine to calculate dielectric properties by solving the Bethe-Salpeter equation in
  Frequency-Reciprocal space on a transition (electron-hole) basis set.

COPYRIGHT

 Copyright (C) 1992-2009 EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida)
 Copyright (C) 2009-2022 ABINIT group (MG, YG)
  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

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 module m_bethe_salpeter
26 
27  use defs_basis
28  use defs_wvltypes
29  use m_bs_defs
30  use m_abicore
31  use m_xmpi
32  use m_errors
33  use m_screen
34  use m_nctk
35  use m_distribfft
36  use netcdf
37  use m_hdr
38  use m_dtset
39  use m_dtfil
40  use m_crystal
41 
42  use defs_datatypes,    only : pseudopotential_type, ebands_t
43  use defs_abitypes,     only : MPI_type
44  use m_gwdefs,          only : GW_Q0_DEFAULT
45  use m_time,            only : timab
46  use m_fstrings,        only : strcat, sjoin, endswith, itoa
47  use m_io_tools,        only : file_exists, iomode_from_fname
48  use m_geometry,        only : mkrdim, metric, normv
49  use m_hide_lapack,     only : matrginv
50  use m_mpinfo,          only : destroy_mpi_enreg, initmpi_seq
51  use m_fftcore,         only : print_ngfft
52  use m_fft_mesh,        only : rotate_FFT_mesh, get_gftt, setmesh
53  use m_fft,             only : fourdp
54  use m_bz_mesh,         only : kmesh_t, get_ng0sh, find_qmesh, make_mesh
55  use m_double_grid,     only : double_grid_t, double_grid_init, double_grid_free
56  use m_ebands,          only : ebands_init, ebands_print, ebands_copy, ebands_free, &
57                                ebands_update_occ, ebands_get_valence_idx, ebands_apply_scissors, ebands_report_gap
58  use m_kg,              only : getph
59  use m_gsphere,         only : gsphere_t
60  use m_vcoul,           only : vcoul_t
61  use m_qparticles,      only : rdqps, rdgw  !, show_QP , rdgw
62  use m_wfd,             only : wfd_init, wfdgw_t, test_charge
63  use m_wfk,             only : wfk_read_eigenvalues
64  use m_energies,        only : energies_type, energies_init
65  use m_io_screening,    only : hscr_t, hscr_io
66  use m_haydock,         only : exc_haydock_driver
67  use m_exc_diago,       only : exc_diago_driver
68  use m_exc_analyze,     only : exc_den
69  use m_eprenorms,       only : eprenorms_t, eprenorms_free, eprenorms_from_epnc, eprenorms_bcast
70  use m_pawang,          only : pawang_type
71  use m_pawrad,          only : pawrad_type
72  use m_pawtab,          only : pawtab_type, pawtab_print, pawtab_get_lsize
73  use m_paw_an,          only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
74  use m_paw_ij,          only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
75  use m_pawfgrtab,       only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init
76  use m_pawrhoij,        only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free,&
77 &                              pawrhoij_inquire_dim, pawrhoij_symrhoij
78  use m_pawdij,          only : pawdij, symdij
79  use m_pawfgr,          only : pawfgr_type, pawfgr_init, pawfgr_destroy
80  use m_paw_hr,          only : pawhur_t, pawhur_free, pawhur_init
81  use m_pawpwij,         only : pawpwff_t, pawpwff_init, pawpwff_free
82  use m_paw_sphharm,     only : setsym_ylm
83  use m_paw_denpot,      only : pawdenpot
84  use m_paw_init,        only : pawinit,paw_gencond
85  use m_paw_onsite,      only : pawnabla_init
86  use m_paw_dmft,        only : paw_dmft_type
87  use m_paw_mkrho,       only : denfgr
88  use m_paw_nhat,        only : nhatgrid,pawmknhat
89  use m_paw_tools,       only : chkpawovlp, pawprt
90  use m_paw_correlations,only : pawpuxinit
91  use m_exc_build,       only : exc_build_ham
92  use m_setvtr,          only : setvtr
93  use m_mkrho,           only : prtrhomxmn
94  use m_pspini,          only : pspini
95  use m_drivexc,         only : mkdenpos
96 
97  implicit none
98 
99  private

m_bethe_salpeter/bethe_salpeter [ Functions ]

[ Top ] [ m_bethe_salpeter ] [ Functions ]

NAME

  bethe_salpeter

FUNCTION

  Main routine to calculate dielectric properties by solving the Bethe-Salpeter equation in
  Frequency-Reciprocal space on a transition (electron-hole) basis set.

INPUTS

 acell(3)=Length scales of primitive translations (bohr)
 codvsn=Code version
 Dtfil<datafiles_type>=Variables related to files.
 Dtset<dataset_type>=All input variables for this dataset.
 Pawang<pawang_type)>=PAW angular mesh and related data.
 Pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
 Pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
 Psps<pseudopotential_type>=Variables related to pseudopotentials.
   Before entering the first time in the routine, a significant part of Psps has been initialized :
   the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm,
   and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in
   the call to pspini. The next time the code enters bethe_salpeter, Psps might be identical to the
   one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90.
 rprim(3,3)=Dimensionless real space primitive translations.
 xred(3,natom)=Reduced atomic coordinates.

 Input files used during the calculation.
  KSS        : Kohn Sham electronic structure file.
  SCR (SUSC) : Files containing the symmetrized epsilon^-1 or the irreducible RPA polarizability,
               respectively. Used to construct the screening W.
  GW file    : Optional file with the GW QP corrections.

OUTPUT

  Output is written on the main output file and on the following external files:
   * _RPA_NLF_MDF: macroscopic RPA dielectric function without non-local field effects.
   * _GW_NLF_MDF: macroscopic RPA dielectric function without non-local field effects calculated
                 with GW energies or the scissors operator.
   * _EXC_MDF: macroscopic dielectric function with excitonic effects obtained by solving the
              Bethe-Salpeter problem at different level of sophistication.

NOTES

 ON THE USE OF FFT GRIDS:
 =================
 In case of PAW:
 ---------------
    Two FFT grids are used:
    - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
      It is defined by nfft, ngfft, mgfft, ...
      Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
    - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
      It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
 In case of norm-conserving:
 ---------------------------
    - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
      For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.

SOURCE

 166 subroutine bethe_salpeter(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim,xred)
 167 
 168 !Arguments ------------------------------------
 169 !scalars
 170  character(len=8),intent(in) :: codvsn
 171  type(datafiles_type),intent(inout) :: Dtfil
 172  type(dataset_type),intent(inout) :: Dtset
 173  type(pawang_type),intent(inout) :: Pawang
 174  type(pseudopotential_type),intent(inout) :: Psps
 175 !arrays
 176  real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,Dtset%natom)
 177  type(pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw)
 178  type(pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw)
 179 
 180 !Local variables ------------------------------
 181 !scalars
 182  integer,parameter :: tim_fourdp0=0,level=40,ipert0=0,idir0=0,cplex1=1,master=0,option1=1
 183  integer :: band,cplex_rhoij,spin,ik_ibz,mqmem,iwarn
 184  integer :: has_dijU,has_dijso,gnt_option
 185  integer :: ik_bz,mband
 186  integer :: choice
 187  integer :: ider
 188  integer :: usexcnhat,nfft_osc,mgfft_osc
 189  integer :: isym,izero
 190  integer :: optcut,optgr0,optgr1,optgr2,option,optrad,optrhoij,psp_gencond
 191  integer :: ngrvdw,nhatgrdim,nkxc1,nprocs,nspden_rhoij,nzlmopt,ifft
 192  integer :: my_rank,rhoxsp_method,comm
 193  integer :: mgfftf,spin_opt,which_fixed
 194  integer :: nscf,nbsc,nkxc,n3xccc
 195  integer :: nfftf,nfftf_tot,nfftot_osc,my_minb,my_maxb
 196  integer :: optene,moved_atm_inside,moved_rhor,initialized,istep,ierr
 197  real(dp) :: ucvol,drude_plsmf,ecore,ecut_eff,ecutdg_eff,norm
 198  real(dp) :: gsqcutc_eff,gsqcutf_eff,gsqcut_shp
 199  real(dp) :: compch_fft,compch_sph,gsq_osc
 200  real(dp) :: vxcavg
 201  logical :: iscompatibleFFT,is_dfpt=.false.,paw_add_onsite,call_pawinit
 202  character(len=500) :: msg
 203  character(len=fnlen) :: wfk_fname,w_fname
 204  type(Pawfgr_type) :: Pawfgr
 205  type(excfiles) :: BS_files
 206  type(excparam) :: BSp
 207  type(paw_dmft_type) :: Paw_dmft
 208  type(MPI_type) :: MPI_enreg_seq
 209  type(crystal_t) :: Cryst
 210  type(kmesh_t) :: Kmesh,Qmesh
 211  type(gsphere_t) :: Gsph_x,Gsph_c,Gsph_x_dense,Gsph_c_dense
 212  type(Hdr_type) :: Hdr_wfk,Hdr_bse
 213  type(ebands_t) :: KS_BSt,QP_BSt,KS_BSt_dense,QP_BSt_dense
 214  type(Energies_type) :: KS_energies
 215  type(vcoul_t) :: Vcp, Vcp_dense
 216  type(wfdgw_t) :: Wfd, Wfd_dense
 217  type(screen_t) :: W
 218  type(screen_info_t) :: W_info
 219  type(wvl_data) :: wvl
 220  type(kmesh_t) :: Kmesh_dense,Qmesh_dense
 221  type(Hdr_type) :: Hdr_wfk_dense
 222  type(double_grid_t) :: grid
 223  type(eprenorms_t) :: Epren
 224 !arrays
 225  integer :: ngfft_osc(18),ngfftc(18),ngfftf(18),nrcell(3)
 226  integer,allocatable :: ktabr(:,:),l_size_atm(:)
 227  integer,allocatable :: nband(:,:),nq_spl(:),irottb(:,:)
 228  integer,allocatable :: qp_vbik(:,:)
 229  integer,allocatable :: gfft_osc(:,:)
 230  real(dp),parameter :: k0(3)=zero
 231  real(dp) :: tsec(2),gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3),rprimd(3,3),eh_rcoord(3),strsxc(6)
 232  real(dp),allocatable :: ph1df(:,:),prev_rhor(:,:),ph1d(:,:)
 233  real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:),ks_rhor(:,:),qp_aerhor(:,:)
 234  real(dp),allocatable :: qp_rhor(:,:),qp_rhog(:,:) !,qp_vhartr(:),qp_vtrial(:,:),qp_vxc(:,:)
 235  real(dp),allocatable :: qp_rhor_paw(:,:),qp_rhor_n_one(:,:),qp_rhor_nt_one(:,:),qp_nhat(:,:)
 236  real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:)
 237  real(dp),allocatable :: vpsp(:),xccc3d(:)
 238  real(dp),allocatable :: ks_vhartr(:),ks_vtrial(:,:),ks_vxc(:,:)
 239  real(dp),allocatable :: kxc(:,:) !,qp_kxc(:,:)
 240  complex(dpc),allocatable :: m_ks_to_qp(:,:,:,:)
 241  logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
 242  type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:)
 243  type(Pawrhoij_type),allocatable :: prev_Pawrhoij(:) !QP_pawrhoij(:),
 244  type(pawpwff_t),allocatable :: Paw_pwff(:)
 245  type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
 246  type(pawhur_t),allocatable :: Hur(:)
 247  type(Paw_ij_type),allocatable :: KS_paw_ij(:)
 248  type(Paw_an_type),allocatable :: KS_paw_an(:)
 249 
 250 !************************************************************************
 251 
 252  DBG_ENTER('COLL')
 253 
 254  call timab(650,1,tsec) ! bse(Total)
 255  call timab(651,1,tsec) ! bse(Init1)
 256 
 257  comm = xmpi_world; nprocs  = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
 258 
 259  wfk_fname = dtfil%fnamewffk
 260 
 261  if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then
 262    ABI_ERROR(msg)
 263  end if
 264  call xmpi_bcast(wfk_fname, master, comm, ierr)
 265 
 266  write(msg,'(8a)')&
 267   ' Exciton: Calculation of dielectric properties by solving the Bethe-Salpeter equation ',ch10,&
 268   ' in frequency domain and reciprocal space on a transitions basis set. ',ch10,&
 269   ' Based on a program developed by L. Reining, V. Olevano, F. Sottile, ',ch10,&
 270   ' S. Albrecht, and G. Onida. Incorporated in ABINIT by M. Giantomassi. ',ch10
 271  call wrtout([std_out, ab_out], msg)
 272 
 273 #ifdef HAVE_GW_DPC
 274  if (gwpc/=8) then
 275    write(msg,'(6a)')ch10,&
 276     ' Number of bytes for double precision complex /=8 ',ch10,&
 277     ' Cannot continue due to kind mismatch in BLAS library ',ch10,&
 278     ' Some BLAS interfaces are not generated by abilint '
 279    ABI_ERROR(msg)
 280  end if
 281  write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10
 282 #else
 283  write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10
 284 #endif
 285  call wrtout([std_out, ab_out], msg)
 286 
 287  !=== Some variables need to be initialized/nullify at start ===
 288  call energies_init(KS_energies)
 289  usexcnhat=0
 290  call mkrdim(acell,rprim,rprimd)
 291  call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
 292  !
 293  !=== Define FFT grid(s) sizes ===
 294  !* Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the
 295  !(usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90.
 296  !See also NOTES in the comments at the beginning of this file.
 297  !NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn.
 298 
 299 !TODO Recheck getng, should use same trick as that used in screening and sigma.
 300  call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,&
 301   gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0)
 302 
 303  call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials')
 304  nfftf_tot=PRODUCT(ngfftf(1:3))
 305 
 306  ! Fake MPI_type for the sequential part.
 307  call initmpi_seq(MPI_enreg_seq)
 308  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all')
 309  call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all')
 310 
 311  ! ===========================================
 312  ! === Open and read pseudopotential files ===
 313  ! ===========================================
 314  call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm)
 315 
 316  ! === Initialization of basic objects including the BSp structure that defines the parameters of the run ===
 317  call setup_bse(codvsn,acell,rprim,ngfft_osc,Dtset,Dtfil,BS_files,Psps,Pawtab,BSp,&
 318    Cryst,Kmesh,Qmesh,KS_BSt,QP_BSt,Hdr_wfk,Gsph_x,Gsph_c,Vcp,Hdr_bse,w_fname,Epren,comm,wvl%descr)
 319 
 320  if (BSp%use_interp) then
 321    call setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh,Kmesh_dense,&
 322      Qmesh_dense,KS_BSt_dense,QP_BSt_dense,Gsph_x_dense,Gsph_c_dense,&
 323      Vcp_dense,Hdr_wfk_dense,grid,comm)
 324  end if
 325 
 326  !call timab(652,2,tsec) ! setup_bse
 327 
 328  nfftot_osc=PRODUCT(ngfft_osc(1:3))
 329  nfft_osc  =nfftot_osc  !no FFT //
 330  mgfft_osc =MAXVAL(ngfft_osc(1:3))
 331 
 332  call print_ngfft(ngfft_osc,header='FFT mesh used for oscillator strengths')
 333 
 334  !TRYING TO RECREATE AN "ABINIT ENVIRONMENT"
 335  KS_energies%e_corepsp=ecore/Cryst%ucvol
 336 
 337  !
 338  !============================
 339  !==== PAW initialization ====
 340  !============================
 341  if (Dtset%usepaw==1) then
 342    call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,xred)
 343 
 344    ABI_MALLOC(KS_Pawrhoij,(Cryst%natom))
 345    call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
 346                nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc)
 347    call pawrhoij_alloc(KS_Pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab)
 348 
 349    ! Initialize values for several basic arrays ===
 350    gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
 351 
 352    ! Test if we have to call pawinit
 353    call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
 354 
 355    if (psp_gencond==1.or.call_pawinit) then
 356      gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
 357      call pawinit(Dtset%effmass_free,gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,&
 358        Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,&
 359        Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%ixc,Dtset%usepotzero)
 360 
 361      ! Update internal values
 362      call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
 363    else
 364      if (Pawtab(1)%has_kij  ==1) Pawtab(1:Cryst%ntypat)%has_kij  =2
 365      if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2
 366    end if
 367    Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore)
 368 
 369    ! Initialize optional flags in Pawtab to zero
 370    ! (Cannot be done in Pawinit since the routine is called only if some parameters are changed)
 371    Pawtab(:)%has_nabla = 0
 372    Pawtab(:)%usepawu   = 0
 373    Pawtab(:)%useexexch = 0
 374    Pawtab(:)%exchmix   =zero
 375    Pawtab(:)%lamb_shielding = zero
 376 
 377    ! * Evaluate <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for the long wavelength limit.
 378    ! TODO solve problem with memory leak and clean this part as well as the associated flag
 379    call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab)
 380 
 381    call setsym_ylm(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,Cryst%rprimd,Cryst%symrec,Pawang%zarot)
 382 
 383    ! Initialize and compute data for DFT+U
 384    Paw_dmft%use_dmft=Dtset%usedmft
 385    call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,&
 386       is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,Dtset%nspinor,Cryst%ntypat,Dtset%optdcmagpawu,Pawang,Dtset%pawprtvol,&
 387       Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu)
 388    if (Dtset%usepawu>0.or.Dtset%useexexch>0) then
 389      ABI_ERROR('BS equation with DFT+U not completely coded!')
 390    end if
 391    if (my_rank == master) call pawtab_print(Pawtab)
 392 
 393    ! Get Pawrhoij from the header of the WFK file.
 394    call pawrhoij_copy(Hdr_wfk%pawrhoij,KS_Pawrhoij)
 395 
 396    ! Re-symmetrize rhoij ===
 397    ! this call leads to a SIGFAULT, likely some pointer is not initialized correctly
 398    choice=1; optrhoij=1
 399 !  call pawrhoij_symrhoij(KS_Pawrhoij,KS_Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,&
 400 !  &             Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,Pawtab,&
 401 !  &             Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
 402 
 403    ! Evaluate form factor of radial part of phi.phj-tphi.tphj ===
 404    rhoxsp_method=1 ! Arnaud-Alouani
 405    if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc
 406 
 407    ABI_MALLOC(gfft_osc,(3,nfftot_osc))
 408    call get_gftt(ngfft_osc,k0,gmet,gsq_osc,gfft_osc)
 409    ABI_FREE(gfft_osc)
 410 
 411    ! Set up q grids, make qmax 20% larger than largest expected:
 412    ABI_MALLOC(nq_spl,(Psps%ntypat))
 413    ABI_MALLOC(qmax,(Psps%ntypat))
 414    nq_spl = Psps%mqgrid_ff
 415    qmax = SQRT(gsq_osc)*1.2d0 ! qmax=Psps%qgrid_ff(Psps%mqgrid_ff)
 416    ABI_MALLOC(Paw_pwff,(Psps%ntypat))
 417 
 418    call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps)
 419 
 420    ABI_FREE(nq_spl)
 421    ABI_FREE(qmax)
 422    !
 423    ! Variables/arrays related to the fine FFT grid ===
 424    ABI_MALLOC(ks_nhat,(nfftf,Dtset%nspden))
 425    ks_nhat=zero
 426    ABI_MALLOC(Pawfgrtab,(Cryst%natom))
 427    call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
 428    call pawfgrtab_init(Pawfgrtab,cplex1,l_size_atm,Dtset%nspden,Dtset%typat)
 429    ABI_FREE(l_size_atm)
 430    compch_fft=greatest_real
 431    usexcnhat=MAXVAL(Pawtab(:)%usexcnhat)
 432    ! * 0 if Vloc in atomic data is Vbare    (Blochl s formulation)
 433    ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation)
 434    write(msg,'(a,i0)')' bethe_salpeter : using usexcnhat = ',usexcnhat
 435    call wrtout(std_out,msg)
 436    !
 437    ! Identify parts of the rectangular grid where the density has to be calculated ===
 438    optcut=0; optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm
 439    if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm
 440 
 441    call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,&
 442    optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
 443  else
 444    ABI_MALLOC(Paw_pwff,(0))
 445  end if !End of PAW Initialization
 446 
 447  ! Consistency check and additional stuff done only for GW with PAW.
 448  if (Dtset%usepaw==1) then
 449    if (Dtset%ecutwfn < Dtset%ecut) then
 450      write(msg,"(5a)")&
 451       "WARNING - ",ch10,&
 452       "  It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,&
 453       "  an excessive truncation of the planewave basis set can lead to unphysical results."
 454      call wrtout(ab_out,msg,'COLL')
 455    end if
 456 
 457    ABI_CHECK(Dtset%usedmft==0,"DMFT + BSE not allowed")
 458    ABI_CHECK(Dtset%useexexch==0,"LEXX + BSE not allowed")
 459  end if
 460 
 461  ! Allocate these arrays anyway, since they are passed to subroutines.
 462  if (.not.allocated(ks_nhat))  then
 463    ABI_MALLOC(ks_nhat,(nfftf,0))
 464  end if
 465 
 466 !==================================================
 467 !==== Read KS band structure from the KSS file ====
 468 !==================================================
 469 
 470  ! Initialize wave function handler, allocate wavefunctions.
 471  my_minb=1; my_maxb=BSp%nbnds; mband=BSp%nbnds
 472  ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol))
 473  nband=mband
 474 
 475 !At present, no memory distribution, each node has the full set of states.
 476  ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol))
 477  bks_mask=.TRUE.
 478 
 479  ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol))
 480  keep_ur=.FALSE.; if (MODULO(Dtset%gwmem,10)==1) keep_ur = .TRUE.
 481 
 482  call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Dtset%nsppol,bks_mask,&
 483   Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_osc,&
 484   Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm)
 485 
 486  ABI_FREE(bks_mask)
 487  ABI_FREE(nband)
 488  ABI_FREE(keep_ur)
 489 
 490  call wfd%print(header="Wavefunctions used to construct the e-h basis set",mode_paral='PERS')
 491 
 492  call timab(651,2,tsec) ! bse(Init1)
 493  call timab(653,1,tsec) ! bse(rdkss)
 494 
 495  call wfd%read_wfk(wfk_fname,iomode_from_fname(wfk_fname))
 496 
 497  ! This test has been disabled (too expensive!)
 498  if (.False.) call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL")
 499 
 500  call timab(653,2,tsec) ! bse(rdkss)
 501  call timab(655,1,tsec) ! bse(mkrho)
 502 
 503  !TODO : check the consistency of Wfd with Wfd_dense !!!
 504  if (BSp%use_interp) then
 505    ! Initialize wave function handler, allocate wavefunctions.
 506    my_minb=1; my_maxb=BSp%nbnds; mband=BSp%nbnds
 507    ABI_MALLOC(nband,(Kmesh_dense%nibz,Dtset%nsppol))
 508    nband=mband
 509 
 510    ! At present, no memory distribution, each node has the full set of states.
 511    ! albeit we allocate only the states that are used.
 512    ABI_MALLOC(bks_mask,(mband,Kmesh_dense%nibz,Dtset%nsppol))
 513    bks_mask=.False.
 514    do spin=1,Bsp%nsppol
 515      bks_mask(Bsp%lomo_spin(spin):Bsp%humo_spin(spin),:,spin) = .True.
 516    end do
 517    !bks_mask=.TRUE.
 518 
 519    ABI_MALLOC(keep_ur,(mband,Kmesh_dense%nibz,Dtset%nsppol))
 520    keep_ur=.FALSE.; if (MODULO(Dtset%gwmem,10)==1) keep_ur = .TRUE.
 521 
 522    call wfd_init(Wfd_dense,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh_dense%nibz,Dtset%nsppol,&
 523     bks_mask,Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,ngfft_osc,&
 524     Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm)
 525 
 526    ABI_FREE(bks_mask)
 527    ABI_FREE(nband)
 528    ABI_FREE(keep_ur)
 529 
 530    call wfd_dense%print(header="Wavefunctions on the dense K-mesh used for interpolation",mode_paral='PERS')
 531 
 532    call wfd_dense%read_wfk(Dtfil%fnameabi_wfkfine, iomode_from_fname(dtfil%fnameabi_wfkfine))
 533    !call wfd_dense%update_bkstab()
 534 
 535    ! This test has been disabled (too expensive!)
 536    if (.False.) call wfd_dense%test_ortho(Cryst,Pawtab,unit=std_out,mode_paral="COLL")
 537  end if
 538 
 539  !=== Calculate the FFT index of $(R^{-1}(r-\tau))$ ===
 540  !* S=\transpose R^{-1} and k_BZ = S k_IBZ
 541  !* irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
 542  ABI_MALLOC(irottb,(nfftot_osc,Cryst%nsym))
 543  call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfft_osc,irottb,iscompatibleFFT)
 544 
 545  ABI_MALLOC(ktabr,(nfftot_osc,Kmesh%nbz))
 546  do ik_bz=1,Kmesh%nbz
 547    isym=Kmesh%tabo(ik_bz)
 548    do ifft=1,nfftot_osc
 549      ktabr(ifft,ik_bz)=irottb(ifft,isym)
 550    end do
 551  end do
 552  ABI_FREE(irottb)
 553  !
 554  !===========================
 555  !=== COMPUTE THE DENSITY ===
 556  !===========================
 557  !* Evaluate Planewave part (complete charge in case of NC pseudos).
 558  !
 559  ABI_MALLOC(ks_rhor,(nfftf,Wfd%nspden))
 560  call wfd%mkrho(Cryst,Psps,Kmesh,KS_BSt,ngfftf,nfftf,ks_rhor)
 561  !
 562  !=== Additional computation for PAW ===
 563  nhatgrdim=0
 564  if (Dtset%usepaw==1) then
 565    !
 566    ! Calculate the compensation charge nhat.
 567    if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc
 568    ider=2*nhatgrdim; izero=0; qphon(:)=zero
 569    if (nhatgrdim>0)  then
 570      ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3))
 571    end if
 572 
 573    call pawmknhat(compch_fft,cplex1,ider,idir0,ipert0,izero,Cryst%gprimd,&
 574      Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
 575      Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,qphon,Cryst%rprimd,&
 576      Cryst%ucvol,Dtset%usewvl,Cryst%xred)
 577 
 578    ! Evaluate onsite energies, potentials, densities ===
 579    !   * Initialize variables/arrays related to the PAW spheres.
 580    !   * Initialize also lmselect (index of non-zero LM-moments of densities).
 581    ABI_MALLOC(KS_paw_ij,(Cryst%natom))
 582    call paw_ij_nullify(KS_paw_ij)
 583 
 584    has_dijso=Dtset%pawspnorb
 585    has_dijU=merge(0,1,Dtset%usepawu==0)
 586 
 587    call paw_ij_init(KS_paw_ij,cplex1,Dtset%nspinor,Dtset%nsppol,&
 588      Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
 589      has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=0,has_dijxc_hat=0,has_dijxc_val=0,&
 590      has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1)
 591 
 592    ABI_MALLOC(KS_paw_an,(Cryst%natom))
 593    call paw_an_nullify(KS_paw_an)
 594 
 595    nkxc1=0
 596    call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
 597      cplex1,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0)
 598 
 599    ! Calculate onsite vxc with and without core charge ===
 600    nzlmopt=-1; option=0; compch_sph=greatest_real
 601 
 602    call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert0,&
 603      Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
 604      Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,&
 605      Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,&
 606      Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp)
 607  end if !PAW
 608 
 609  if (.not.allocated(ks_nhatgr))  then
 610    ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,0))
 611  end if
 612 
 613  call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,&
 614    Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
 615 
 616  ! === For PAW, add the compensation charge on the FFT mesh, then get rho(G) ===
 617  if (Dtset%usepaw==1) ks_rhor(:,:)=ks_rhor(:,:)+ks_nhat(:,:)
 618  call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_rhor,ucvol=ucvol)
 619 
 620  ABI_MALLOC(ks_rhog,(2,nfftf))
 621 
 622  call fourdp(1,ks_rhog,ks_rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,tim_fourdp0)
 623  call timab(655,2,tsec) ! bse(mkrho)
 624 
 625  !
 626  ! The following steps have been gathered in the setvtr routine:
 627  ! - get Ewald energy and Ewald forces
 628  ! - compute local ionic pseudopotential vpsp
 629  ! - eventually compute 3D core electron density xccc3d
 630  ! - eventually compute vxc and vhartr
 631  ! - set up ks_vtrial
 632  !
 633  ! *******************************************************************
 634  ! **** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
 635  ! *******************************************************************
 636  ngrvdw=0
 637  ABI_MALLOC(grvdw,(3,ngrvdw))
 638  ABI_MALLOC(grchempottn,(3,Cryst%natom))
 639  ABI_MALLOC(grewtn,(3,Cryst%natom))
 640  nkxc=0
 641 !if (Wfd%nspden==1) nkxc=2
 642 !if (Wfd%nspden>=2) nkxc=3 ! check GGA and spinor, quite a messy part!!!
 643  ABI_MALLOC(kxc,(nfftf,nkxc))
 644 
 645  n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf
 646  ABI_MALLOC(xccc3d,(n3xccc))
 647  ABI_MALLOC(ks_vhartr,(nfftf))
 648  ABI_MALLOC(ks_vtrial,(nfftf,Wfd%nspden))
 649  ABI_MALLOC(vpsp,(nfftf))
 650  ABI_MALLOC(ks_vxc,(nfftf,Wfd%nspden))
 651 
 652  optene=4; moved_atm_inside=0; moved_rhor=0; initialized=1; istep=1
 653 !
 654 !=== Compute structure factor phases and large sphere cut-off ===
 655 !WARNING cannot use Dtset%mgfft, this has to be checked better
 656 !mgfft=MAXVAL(ngfftc(:))
 657 !allocate(ph1d(2,3*(2*mgfft+1)*Cryst%natom),ph1df(2,3*(2*mgfftf+1)*Cryst%natom))
 658  write(std_out,*)' CHECK ',Dtset%mgfftdg,mgfftf
 659  if (Dtset%mgfftdg/=mgfftf) then
 660 !  write(std_out,*)"WARNING Dtset%mgfftf /= mgfftf"
 661 !  write(std_out,*)'HACKING Dtset%mgfftf'
 662 !  Dtset%mgfftdg=mgfftf
 663  end if
 664  ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom))
 665  ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom))
 666 
 667  call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred)
 668 
 669  if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then
 670    call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred)
 671  else
 672    ph1df(:,:)=ph1d(:,:)
 673  end if
 674 
 675  ABI_FREE(ph1d)
 676 
 677  call setvtr(Cryst%atindx1,Dtset,KS_energies,Cryst%gmet,Cryst%gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
 678 & istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
 679 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
 680 & optene,Pawrad,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,Cryst%rprimd,strsxc,&
 681 & Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,wvl,xccc3d,Cryst%xred)
 682 
 683  ABI_FREE(ph1df)
 684  ABI_FREE(vpsp)
 685 
 686  ! ============================
 687  ! ==== Compute KS PAW Dij ====
 688  ! ============================
 689  if (Wfd%usepaw==1) then
 690    call timab(561,1,tsec)
 691    !
 692    ! Calculate the unsymmetrized Dij.
 693    call pawdij(cplex1,Dtset%enunit,Cryst%gprimd,ipert0,&
 694      Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
 695      Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,&
 696      Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
 697      k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,ks_vxc,Cryst%xred,&
 698      nucdipmom=Dtset%nucdipmom)
 699 
 700    ! Symmetrize KS Dij
 701    call symdij(Cryst%gprimd,Cryst%indsym,ipert0,&
 702      Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,0,KS_paw_ij,Pawang,&
 703      Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
 704 
 705    ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
 706    call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab)
 707    call timab(561,2,tsec)
 708  end if
 709 
 710  ABI_FREE(kxc)
 711  ABI_FREE(xccc3d)
 712  ABI_FREE(grchempottn)
 713  ABI_FREE(grewtn)
 714  ABI_FREE(grvdw)
 715 
 716  !=== QP_BSt stores energies and occ. used for the calculation ===
 717  !* Initialize QP_BSt with KS values.
 718  !* In case of SC update QP_BSt using the QPS file.
 719  ABI_MALLOC(qp_rhor,(nfftf,Dtset%nspden))
 720  qp_rhor=ks_rhor
 721 
 722  ! AE density used for the model dielectric function.
 723  ABI_MALLOC(qp_aerhor,(nfftf,Dtset%nspden))
 724  qp_aerhor=ks_rhor
 725 
 726  ! PAW: Compute AE rhor. Under testing
 727  if (Wfd%usepaw==1 .and. BSp%mdlf_type/=0) then
 728    ABI_WARNING("Entering qp_aerhor with PAW")
 729 
 730    ABI_MALLOC(qp_rhor_paw   ,(nfftf,Wfd%nspden))
 731    ABI_MALLOC(qp_rhor_n_one ,(nfftf,Wfd%nspden))
 732    ABI_MALLOC(qp_rhor_nt_one,(nfftf,Wfd%nspden))
 733 
 734    ABI_MALLOC(qp_nhat,(nfftf,Wfd%nspden))
 735    qp_nhat = ks_nhat
 736    ! TODO: I pass KS_pawrhoij instead of QP_pawrhoij but in the present version there's no difference.
 737 
 738    call denfgr(Cryst%atindx1,Cryst%gmet,Wfd%comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,qp_nhat,&
 739 &   Wfd%nspinor,Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,KS_pawrhoij,Pawtab,Dtset%prtvol,&
 740 &   qp_rhor,qp_rhor_paw,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
 741 
 742    norm = SUM(qp_rhor_paw(:,1))*Cryst%ucvol/PRODUCT(Pawfgr%ngfft(1:3))
 743    write(msg,'(a,F8.4)') '  QUASIPARTICLE DENSITY CALCULATED - NORM OF DENSITY: ',norm
 744    call wrtout(std_out,msg)
 745    write(std_out,*)"MAX", MAXVAL(qp_rhor_paw(:,1))
 746    write(std_out,*)"MIN", MINVAL(qp_rhor_paw(:,1))
 747 
 748    ABI_FREE(qp_nhat)
 749    ABI_FREE(qp_rhor_n_one)
 750    ABI_FREE(qp_rhor_nt_one)
 751 
 752    ! Use ae density for the model dielectric function.
 753    iwarn=0
 754    call mkdenpos(iwarn,nfftf,Wfd%nspden,option1,qp_rhor_paw,dtset%xc_denpos)
 755    qp_aerhor = qp_rhor_paw
 756    ABI_FREE(qp_rhor_paw)
 757  end if
 758 
 759 !call copy_bandstructure(KS_BSt,QP_BSt)
 760 
 761  if (.FALSE.) then
 762 !  $ m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}> $
 763    ABI_MALLOC(m_ks_to_qp,(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol))
 764    m_ks_to_qp=czero
 765    do spin=1,Wfd%nsppol
 766      do ik_ibz=1,Wfd%nkibz
 767        do band=1,Wfd%nband(ik_ibz,spin)
 768          m_ks_to_qp(band,band,ik_ibz,spin)=cone ! Initialize the QP amplitudes with KS wavefunctions.
 769        end do
 770      end do
 771    end do
 772    !
 773    ! Now read m_ks_to_qp and update the energies in QP_BSt.
 774    ! TODO switch on the renormalization of n in sigma.
 775    ABI_MALLOC(prev_rhor,(nfftf,Wfd%nspden))
 776    ABI_MALLOC(prev_Pawrhoij,(Cryst%natom*Wfd%usepaw))
 777 
 778    call rdqps(QP_BSt,Dtfil%fnameabi_qps,Wfd%usepaw,Wfd%nspden,1,nscf,&
 779     nfftf,ngfftf,Cryst%ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_ks_to_qp,prev_rhor,prev_Pawrhoij)
 780 
 781    ABI_FREE(prev_rhor)
 782    if (Psps%usepaw==1.and.nscf>0) then
 783      call pawrhoij_free(prev_pawrhoij)
 784    end if
 785    ABI_FREE(prev_pawrhoij)
 786    !
 787    !  if (nscf>0.and.wfd_iam_master(Wfd)) then ! Print the unitary transformation on std_out.
 788    !  call show_QP(QP_BSt,m_ks_to_qp,fromb=Sigp%minbdgw,tob=Sigp%maxbdgw,unit=std_out,tolmat=0.001_dp)
 789    !  end if
 790    !
 791    !  === Compute QP wfg as linear combination of KS states ===
 792    !  * Wfd%ug is modified inside calc_wf_qp
 793    !  * For PAW, update also the on-site projections.
 794    !  * WARNING the first dimension of MPI_enreg MUST be Kmesh%nibz
 795    !  TODO here we should use nbsc instead of nbnds
 796 
 797    call wfd%rotate(Cryst,m_ks_to_qp)
 798 
 799    ABI_FREE(m_ks_to_qp)
 800    !
 801    !  === Reinit the storage mode of Wfd as ug have been changed ===
 802    !  * Update also the wavefunctions for GW corrections on each processor
 803    call wfd%reset_ur_cprj()
 804 
 805    call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL")
 806 
 807    ! Compute QP occupation numbers.
 808    call wrtout(std_out,'bethe_salpeter: calculating QP occupation numbers')
 809 
 810    call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=0)
 811    ABI_MALLOC(qp_vbik,(QP_BSt%nkpt,QP_BSt%nsppol))
 812    qp_vbik(:,:) = ebands_get_valence_idx(QP_BSt)
 813    ABI_FREE(qp_vbik)
 814 
 815    call wfd%mkrho(Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,qp_rhor)
 816  end if
 817 
 818  ABI_MALLOC(qp_rhog,(2,nfftf))
 819  call fourdp(1,qp_rhog,qp_rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,0)
 820 
 821  ! States up to lomo-1 are useless now since only the states close to the gap are
 822  ! needed to construct the EXC Hamiltonian. Here we deallocate the wavefunctions
 823  ! to make room for the excitonic Hamiltonian that is allocated in exc_build_ham.
 824  ! and for the screening that is allocated below.
 825  ! Hereafter bands from 1 up to lomo-1 and all bands above humo+1 should not be accessed!
 826  ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol))
 827 
 828  bks_mask=.FALSE.
 829  do spin=1,Bsp%nsppol
 830    if (Bsp%lomo_spin(spin)>1) bks_mask(1:Bsp%lomo_spin(spin)-1,:,spin)=.TRUE.
 831    if (Bsp%humo_spin(spin)+1<=Wfd%mband) bks_mask(Bsp%humo_spin(spin)+1:,:,spin)=.TRUE.
 832  end do
 833  call wfd%wave_free(what="All", bks_mask=bks_mask)
 834  ABI_FREE(bks_mask)
 835  !
 836  ! ================================================================
 837  ! Build the screened interaction W in the irreducible wedge.
 838  ! * W(q,G1,G2) = vc^{1/2} (q,G1) e^{-1}(q,G1,G2) vc^{1/2) (q,G2)
 839  ! * Use Coulomb term for q-->0,
 840  ! * Only the first small Q is used, shall we average if nqlwl>1?
 841  ! ================================================================
 842  ! TODO clean this part and add an option to retrieve a single frequency to save memory.
 843  call timab(654,1,tsec) ! bse(rdmkeps^-1)
 844 
 845  call screen_nullify(W)
 846  if (BSp%use_coulomb_term) then !  Init W.
 847    ! Incore or out-of-core solution?
 848    mqmem=0; if (Dtset%gwmem/10==1) mqmem=Qmesh%nibz
 849 
 850    W_info%invalid_freq = Dtset%gw_invalid_freq
 851    W_info%mat_type = MAT_INV_EPSILON
 852    !W_info%mat_type = MAT_W
 853    !W_info%vtx_family
 854    !W_info%ixc
 855    !W_info%use_ada
 856    W_info%use_mdf = BSp%mdlf_type
 857    !W_info%use_ppm
 858    !W_info%vtx_test
 859    !W_info%wint_method
 860    !
 861    !W_info%ada_kappa
 862    W_info%eps_inf = BSp%eps_inf
 863    !W_info%drude_plsmf
 864 
 865    call screen_init(W,W_Info,Cryst,Qmesh,Gsph_c,Vcp,w_fname,mqmem,Dtset%npweps,&
 866     Dtset%iomode,ngfftf,nfftf_tot,Wfd%nsppol,Wfd%nspden,qp_aerhor,Wfd%prtvol,Wfd%comm)
 867  end if
 868  call timab(654,2,tsec) ! bse(rdmkeps^-1)
 869  !
 870  ! =============================================
 871  ! ==== Build of the excitonic Hamiltonian =====
 872  ! =============================================
 873  !
 874  call timab(656,1,tsec) ! bse(mkexcham)
 875 
 876  call exc_build_ham(BSp,BS_files,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,&
 877 & Wfd,W,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff)
 878  !
 879  ! Free Em1 to make room for the full excitonic Hamiltonian.
 880  call screen_free(W)
 881 
 882  call timab(656,2,tsec) ! bse(mkexcham)
 883  !
 884  ! =========================================
 885  ! ==== Macroscopic dielectric function ====
 886  ! =========================================
 887  call timab(657,1,tsec) ! bse(mkexceps)
 888  !
 889  ! First deallocate the internal %ur buffers to make room for the excitonic Hamiltonian.
 890  call timab(658,1,tsec) ! bse(wfd_wave_free)
 891  ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol))
 892  bks_mask=.TRUE.
 893  call wfd%wave_free(what="Real_space", bks_mask=bks_mask)
 894  ABI_FREE(bks_mask)
 895  call timab(658,2,tsec) ! bse(wfd_wave_free)
 896 
 897  ! Compute the commutator [r,Vu] (PAW only).
 898  ABI_MALLOC(HUr,(Cryst%natom*Wfd%usepaw))
 899 
 900  call timab(659,1,tsec) ! bse(make_pawhur_t)
 901  if (Bsp%inclvkb/=0 .and. Wfd%usepaw==1 .and. Dtset%usepawu/=0) then !TODO here I need KS_Paw_ij
 902    ABI_WARNING("Commutator for DFT+U not tested")
 903    call pawhur_init(hur,Wfd%nsppol,Wfd%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,KS_Paw_ij)
 904  end if
 905  call timab(659,2,tsec) ! bse(make_pawhur_t)
 906 
 907  select case (BSp%algorithm)
 908  case (BSE_ALGO_NONE)
 909    ABI_COMMENT("Skipping solution of the BSE equation")
 910 
 911  case (BSE_ALGO_DDIAGO, BSE_ALGO_CG)
 912    call timab(660,1,tsec) ! bse(exc_diago_driver)
 913    call exc_diago_driver(Wfd,Bsp,BS_files,KS_BSt,QP_BSt,Cryst,Kmesh,Psps, Pawtab,Hur,Hdr_bse,drude_plsmf,Epren)
 914    call timab(660,2,tsec) ! bse(exc_diago_driver)
 915 
 916    if (.FALSE.) then ! Calculate electron-hole excited state density. Not tested at all.
 917      call exc_den(BSp,BS_files,ngfftf,nfftf,Kmesh,ktabr,Wfd)
 918    end if
 919 
 920    if (.FALSE.) then
 921      paw_add_onsite=.FALSE.; spin_opt=1; which_fixed=1; eh_rcoord=(/zero,zero,zero/); nrcell=(/2,2,2/)
 922 !    call exc_plot(Bsp,Bs_files,Wfd,Kmesh,Cryst,Psps,Pawtab,Pawrad,paw_add_onsite,spin_opt,which_fixed,eh_rcoord,nrcell,ngfftf)
 923    end if
 924 
 925    if (BSp%use_interp) then
 926      ABI_ERROR("Interpolation technique not coded for diagonalization and CG")
 927    end if
 928 
 929  case (BSE_ALGO_Haydock)
 930    call timab(661,1,tsec) ! bse(exc_haydock_driver)
 931 
 932    if (BSp%use_interp) then
 933 
 934      call exc_haydock_driver(BSp,BS_files,Cryst,Kmesh,Hdr_bse,KS_BSt,QP_BSt,Wfd,Psps,Pawtab,Hur,Epren, &
 935        kmesh_dense=Kmesh_dense, ks_bst_dense=KS_BSt_dense, qp_bst_dense=QP_BSt_dense,wfd_dense=Wfd_dense, &
 936        vcp_dense=Vcp_dense, grid=grid)
 937    else
 938      call exc_haydock_driver(BSp,BS_files,Cryst,Kmesh,Hdr_bse,KS_BSt,QP_BSt,Wfd,Psps,Pawtab,Hur,Epren)
 939    end if
 940 
 941    call timab(661,2,tsec) ! bse(exc_haydock_driver)
 942 
 943  case default
 944    write(msg,'(a,i0)')" Wrong BSE algorithm: ",BSp%algorithm
 945    ABI_ERROR(msg)
 946  end select
 947 
 948  call timab(657,2,tsec) ! bse(mkexceps)
 949 
 950  !=====================
 951  !==== Free memory ====
 952  !=====================
 953  ABI_FREE(ktabr)
 954  ABI_FREE(ks_vhartr)
 955  ABI_FREE(ks_vtrial)
 956  ABI_FREE(ks_vxc)
 957  ABI_FREE(ks_nhat)
 958  ABI_FREE(ks_nhatgr)
 959  ABI_FREE(ks_rhog)
 960  ABI_FREE(ks_rhor)
 961  ABI_FREE(qp_rhog)
 962  ABI_FREE(qp_rhor)
 963  ABI_FREE(qp_aerhor)
 964  !
 965  ! Free local data structures.
 966  call destroy_mpi_enreg(MPI_enreg_seq)
 967  call cryst%free()
 968  call Gsph_x%free()
 969  call Gsph_c%free()
 970  call Kmesh%free()
 971  call Qmesh%free()
 972  call Hdr_wfk%free()
 973  call Hdr_bse%free()
 974  call ebands_free(KS_BSt)
 975  call ebands_free(QP_BSt)
 976  call Vcp%free()
 977  call pawhur_free(Hur)
 978  ABI_FREE(Hur)
 979  call bs_parameters_free(BSp)
 980  call wfd%free()
 981  call pawfgr_destroy(Pawfgr)
 982  call eprenorms_free(Epren)
 983 
 984  ! Free memory used for interpolation.
 985  if (BSp%use_interp) then
 986    call double_grid_free(grid)
 987    call wfd_dense%free()
 988    call Gsph_x_dense%free()
 989    call Gsph_c_dense%free()
 990    call Kmesh_dense%free()
 991    call Qmesh_dense%free()
 992    call ebands_free(KS_BSt_dense)
 993    call ebands_free(QP_BSt_dense)
 994    call Vcp_dense%free()
 995    call Hdr_wfk_dense%free()
 996  end if
 997 
 998  ! Optional deallocation for PAW.
 999  if (Dtset%usepaw==1) then
1000    call pawrhoij_free(KS_Pawrhoij)
1001    ABI_FREE(KS_Pawrhoij)
1002    call pawfgrtab_free(Pawfgrtab)
1003    ABI_FREE(Pawfgrtab)
1004    call paw_ij_free(KS_paw_ij)
1005    ABI_FREE(KS_paw_ij)
1006    call paw_an_free(KS_paw_an)
1007    ABI_FREE(KS_paw_an)
1008    call pawpwff_free(Paw_pwff)
1009  end if
1010  ABI_FREE(Paw_pwff)
1011 
1012  call timab(650,2,tsec) ! bse(Total)
1013 
1014  DBG_EXIT('COLL')
1015 
1016 end subroutine bethe_salpeter

m_bethe_salpeter/setup_bse [ Functions ]

[ Top ] [ m_bethe_salpeter ] [ Functions ]

NAME

  setup_bse

FUNCTION

  This routine performs the initialization of basic objects and quantities used for Bethe-Salpeter calculations.
  In particular the excparam data type that defines the parameters of the calculation is completely
  initialized starting from the content of Dtset and the parameters read from the external WFK and SCR (SUSC) file.

INPUTS

 codvsn=Code version
 ngfft_gw(18)=Information about 3D FFT for density and potentials, see ~abinit/doc/variables/vargs.htm#ngfft
 acell(3)=Length scales of primitive translations (bohr)
 rprim(3,3)=Dimensionless real space primitive translations.
 Dtset<dataset_type>=All input variables for this dataset.
  Some of them might be redefined here TODO
 Dtfil=filenames and unit numbers used in abinit.
 Psps <pseudopotential_type>=variables related to pseudopotentials
 Pawtab(Psps%ntypat*Dtset%usepaw)<pawtab_type>=PAW tabulated starting data

OUTPUT

 Cryst<crystal_t>=Info on the crystalline Structure.
 Kmesh<kmesh_t>=Structure defining the k-sampling for the wavefunctions.
 Qmesh<kmesh_t>=Structure defining the q-sampling for the symmetrized inverse dielectric matrix.
 Gsph_x<gsphere_t=Data type gathering info on the G-sphere for wave functions and e^{-1},
 KS_BSt<ebands_t>=The KS band structure (energies, occupancies, k-weights...)
 Vcp<vcoul_t>=Structure gathering information on the Coulomb interaction in reciprocal space,
   including a possible cutoff in real space.
 ngfft_osc(18)=Contain all needed information about the 3D FFT for the oscillator matrix elements.
   See ~abinit/doc/variables/vargs.htm#ngfft
 Bsp<excparam>=Basic parameters defining the Bethe-Salpeter run. Completely initialed in output.
 Hdr_wfk<Hdr_type>=The header of the WFK file.
 Hdr_bse<Hdr_type>=Local header initialized from the parameters used for the Bethe-Salpeter calculation.
 BS_files<excfiles>=Files used in the calculation.
 w_file=File name used to construct W. Set to ABI_NOFILE if no external file is used.

SOURCE

1057 subroutine setup_bse(codvsn,acell,rprim,ngfft_osc,Dtset,Dtfil,BS_files,Psps,Pawtab,BSp,&
1058 & Cryst,Kmesh,Qmesh,KS_BSt,QP_bst,Hdr_wfk,Gsph_x,Gsph_c,Vcp,Hdr_bse,w_fname,Epren,comm,Wvl)
1059 
1060 !Arguments ------------------------------------
1061 !scalars
1062  integer,intent(in) :: comm
1063  character(len=8),intent(in) :: codvsn
1064  character(len=fnlen),intent(out) :: w_fname
1065  type(dataset_type),intent(inout) :: Dtset
1066  type(datafiles_type),intent(in) :: Dtfil
1067  type(pseudopotential_type),intent(in) :: Psps
1068  type(excparam),intent(inout) :: Bsp
1069  type(hdr_type),intent(out) :: Hdr_wfk,Hdr_bse
1070  type(crystal_t),intent(out) :: Cryst
1071  type(kmesh_t),intent(out) :: Kmesh,Qmesh
1072  type(gsphere_t),intent(out) :: Gsph_x,Gsph_c
1073  type(ebands_t),intent(out) :: KS_BSt,QP_Bst
1074  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw)
1075  type(vcoul_t),intent(out) :: Vcp
1076  type(excfiles),intent(out) :: BS_files
1077  type(wvl_internal_type), intent(in) :: Wvl
1078  type(eprenorms_t),intent(out) :: Epren
1079 !arrays
1080  integer,intent(out) :: ngfft_osc(18)
1081  real(dp),intent(in) :: acell(3),rprim(3,3)
1082 
1083 !Local variables ------------------------------
1084 !scalars
1085  integer,parameter :: pertcase0=0,master=0
1086  integer(i8b) :: work_size,tot_nreh,neh_per_proc,il
1087  integer :: bantot,enforce_sym,ib,ibtot,ik_ibz,isppol,jj,method,iat,ount !ii,
1088  integer :: mband,io,nfftot_osc,spin,hexc_size,nqlwl,iq
1089  integer :: timrev,iq_bz,isym,iq_ibz,itim
1090  integer :: my_rank,nprocs,fform,npwe_file,ierr,my_k1, my_k2,my_nbks
1091  integer :: first_dig,second_dig,it
1092  real(dp) :: ucvol,qnorm
1093  real(dp):: eff,mempercpu_mb,wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem
1094  logical,parameter :: remove_inv=.FALSE.
1095  logical :: ltest,occ_from_dtset
1096  character(len=500) :: msg
1097  character(len=fnlen) :: gw_fname,test_file,wfk_fname
1098  character(len=fnlen) :: ep_nc_fname
1099  type(hscr_t) :: Hscr
1100 !arrays
1101  integer :: ng0sh_opt(3),val_idx(Dtset%nsppol)
1102  integer,allocatable :: npwarr(:),val_indeces(:,:),nlmn_atm(:)
1103  real(dp) :: qpt_bz(3),minmax_tene(2)
1104  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3)
1105  real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:)
1106  real(dp),allocatable :: igwene(:,:,:)
1107  real(dp),pointer :: energies_p(:,:,:)
1108  complex(dpc),allocatable :: gw_energy(:,:,:)
1109  type(Pawrhoij_type),allocatable :: Pawrhoij(:)
1110 
1111 !************************************************************************
1112 
1113  DBG_ENTER("COLL")
1114 
1115  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1116 
1117  ! === Check for calculations that are not implemented ===
1118  ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1))
1119  ABI_CHECK(ltest,'Dtset%nband must be constant')
1120  ABI_CHECK(Dtset%nspinor==1,"nspinor==2 not coded")
1121 
1122  ! === Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume ===
1123  call mkrdim(acell,rprim,rprimd)
1124  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1125 
1126  ! Read energies and header from the WFK file.
1127  wfk_fname = dtfil%fnamewffk
1128  if (.not. file_exists(wfk_fname)) then
1129    wfk_fname = nctk_ncify(wfk_fname)
1130    ABI_COMMENT(sjoin("File not found. Will try netcdf file: ", wfk_fname))
1131  end if
1132 
1133  call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm)
1134  mband = MAXVAL(Hdr_wfk%nband)
1135 
1136  call hdr_wfk%vs_dtset(dtset)
1137 
1138  ! === Create crystal_t data type ===
1139  !remove_inv= .FALSE. !(nsym_kss/=Hdr_wfk%nsym)
1140  timrev=  2 ! This information is not reported in the header
1141             ! 1 => do not use time-reversal symmetry
1142             ! 2 => take advantage of time-reversal symmetry
1143 
1144  cryst = Hdr_wfk%get_crystal(gw_timrev=timrev, remove_inv=remove_inv)
1145  call cryst%print()
1146  !
1147  ! Setup of the k-point list and symmetry tables in the  BZ -----------------------------------
1148  if (Dtset%chksymbreak==0) then
1149    ABI_WARNING("Calling make_mesh")
1150    call make_mesh(Kmesh,Cryst,Dtset%kptopt,Dtset%kptrlatt,Dtset%nshiftk,Dtset%shiftk,break_symmetry=.TRUE.)
1151    ! TODO
1152    !Check if kibz from KSS file corresponds to the one returned by make_mesh.
1153  else
1154    call Kmesh%init(Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt)
1155  end if
1156  BSp%nkibz = Kmesh%nibz  !We might allow for a smaller number of points....
1157 
1158  call Kmesh%print("K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL")
1159  call Kmesh%print("K-mesh for the wavefunctions",ab_out, 0,           "COLL")
1160 
1161  nqlwl = 0
1162  w_fname = ABI_NOFILE
1163  if (Dtset%getscr/=0 .or. Dtset%irdscr/=0 .or. dtset%getscr_filepath /= ABI_NOFILE) then
1164    w_fname=Dtfil%fnameabi_scr
1165  else if (Dtset%getsuscep/=0.or.Dtset%irdsuscep/=0) then
1166    w_fname=Dtfil%fnameabi_sus
1167    ABI_ERROR("(get|ird)suscep not implemented")
1168  end if
1169 
1170  if (w_fname /= ABI_NOFILE) then ! Read dimensions from the external file
1171 
1172    if (.not. file_exists(w_fname)) then
1173      w_fname = nctk_ncify(w_fname)
1174      ABI_COMMENT(sjoin("File not found. Will try netcdf file: ", w_fname))
1175    end if
1176 
1177    if (my_rank==master) then
1178      ! Master reads npw and nqlwl from SCR file.
1179      call wrtout(std_out,sjoin('Testing file: ', w_fname))
1180 
1181      call hscr%from_file(w_fname, fform, xmpi_comm_self)
1182      ! Echo the header.
1183      if (Dtset%prtvol>0) call Hscr%print()
1184 
1185      npwe_file = Hscr%npwe ! Have to change %npweps if it was larger than dim on disk.
1186      nqlwl     = Hscr%nqlwl
1187 
1188      if (Dtset%npweps>npwe_file) then
1189        write(msg,'(2(a,i0),2a,i0)')&
1190         "The number of G-vectors stored on file (",npwe_file,") is smaller than Dtset%npweps = ",Dtset%npweps,ch10,&
1191         "Calculation will proceed with the maximum available set, npwe_file = ",npwe_file
1192        ABI_WARNING(msg)
1193        Dtset%npweps = npwe_file
1194      else  if (Dtset%npweps<npwe_file) then
1195        write(msg,'(2(a,i0),2a,i0)')&
1196         "The number of G-vectors stored on file (",npwe_file,") is larger than Dtset%npweps = ",Dtset%npweps,ch10,&
1197         "Calculation will proceed with Dtset%npweps = ",Dtset%npweps
1198        ABI_COMMENT(msg)
1199      end if
1200    end if
1201 
1202    call Hscr%bcast(master, my_rank, comm)
1203    call xmpi_bcast(Dtset%npweps,master,comm,ierr)
1204    call xmpi_bcast(nqlwl,master,comm,ierr)
1205 
1206    if (nqlwl>0) then
1207      ABI_MALLOC(qlwl,(3,nqlwl))
1208      qlwl = Hscr%qlwl
1209    end if
1210    !
1211    ! Init Qmesh from the SCR file.
1212    call Qmesh%init(Cryst,Hscr%nqibz,Hscr%qibz,Dtset%kptopt)
1213 
1214    ! The G-sphere for W and Sigma_c is initialized from the gvectors found in the SCR file.
1215    call Gsph_c%init(Cryst, Dtset%npweps, gvec=Hscr%gvec)
1216 
1217    call Hscr%free()
1218  else
1219    ! Init Qmesh from the K-mesh reported in the WFK file.
1220    call find_qmesh(Qmesh,Cryst,Kmesh)
1221 
1222    ! The G-sphere for W and Sigma_c is initialized from ecuteps.
1223    call Gsph_c%init(Cryst, 0, ecut=Dtset%ecuteps)
1224    Dtset%npweps = Gsph_c%ng
1225  end if
1226 
1227  BSp%npweps = Dtset%npweps
1228  BSp%ecuteps = Dtset%ecuteps
1229 
1230  if (nqlwl==0) then
1231    nqlwl=1
1232    ABI_MALLOC(qlwl,(3,nqlwl))
1233    qlwl(:,nqlwl)= GW_Q0_DEFAULT
1234    write(msg,'(3a,i2,a,3f9.6)')&
1235      "The Header of the screening file does not contain the list of q-point for the optical limit ",ch10,&
1236      "Using nqlwl= ",nqlwl," and qlwl = ",qlwl(:,1)
1237    ABI_COMMENT(msg)
1238  end if
1239  write(std_out,*)"nqlwl and qlwl for Coulomb singularity and e^-1",nqlwl,qlwl
1240 
1241  ! === Setup of q-mesh in the whole BZ ===
1242  ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed,
1243  !   epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed)
1244  !
1245  call Qmesh%print("Q-mesh for the screening function",std_out,Dtset%prtvol,"COLL")
1246  call Qmesh%print("Q-mesh for the screening function",ab_out ,0           ,"COLL")
1247 
1248  do iq_bz=1,Qmesh%nbz
1249    call qmesh%get_BZ_item(iq_bz,qpt_bz,iq_ibz,isym,itim)
1250    sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz))
1251    if (ANY(ABS(Qmesh%bz(:,iq_bz)-sq )>1.0d-4)) then
1252      write(std_out,*) sq,Qmesh%bz(:,iq_bz)
1253      write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')&
1254       'qpoint ',Qmesh%bz(:,iq_bz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,&
1255       'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,&
1256       'however a non zero umklapp G_o vector is required and this is not yet allowed'
1257      ABI_ERROR(msg)
1258    end if
1259  end do
1260 
1261  BSp%algorithm = Dtset%bs_algorithm
1262  BSp%nstates   = Dtset%bs_nstates
1263  Bsp%nsppol    = Dtset%nsppol
1264  Bsp%hayd_term = Dtset%bs_hayd_term
1265 
1266  ! Define the algorithm for solving the BSE.
1267  if (BSp%algorithm == BSE_ALGO_HAYDOCK) then
1268    BSp%niter       = Dtset%bs_haydock_niter
1269    BSp%haydock_tol = Dtset%bs_haydock_tol
1270 
1271  else if (BSp%algorithm == BSE_ALGO_CG) then
1272    ! FIXME For the time being use an hardcoded value.
1273    ! TODO change name in Dtset%
1274    BSp%niter       = Dtset%nstep !100
1275    BSp%cg_tolwfr   = Dtset%tolwfr
1276    BSp%nline       = Dtset%nline
1277    BSp%nbdbuf      = Dtset%nbdbuf
1278    BSp%nstates     = Dtset%bs_nstates
1279    ABI_WARNING("Check CG setup")
1280  else
1281    !BSp%niter       = 0
1282    !BSp%tol_iter    = HUGE(one)
1283  end if
1284  !
1285  ! Shall we include Local field effects?
1286  SELECT CASE (Dtset%bs_exchange_term)
1287  CASE (0,1)
1288    BSp%exchange_term = Dtset%bs_exchange_term
1289  CASE DEFAULT
1290    write(msg,'(a,i0)')" Wrong bs_exchange_term: ",Dtset%bs_exchange_term
1291    ABI_ERROR(msg)
1292  END SELECT
1293  !
1294  ! Treatment of the off-diagonal coupling block.
1295  SELECT CASE (Dtset%bs_coupling)
1296  CASE (0)
1297    BSp%use_coupling = 0
1298    msg = 'RESONANT ONLY CALCULATION'
1299  CASE (1)
1300    BSp%use_coupling = 1
1301    msg = ' RESONANT+COUPLING CALCULATION '
1302  CASE DEFAULT
1303    write(msg,'(a,i0)')" Wrong bs_coupling: ",Dtset%bs_coupling
1304    ABI_ERROR(msg)
1305  END SELECT
1306  call wrtout(std_out,msg)
1307 
1308  BSp%use_diagonal_Wgg = .FALSE.
1309  Bsp%use_coulomb_term = .TRUE.
1310  BSp%eps_inf=zero
1311  Bsp%mdlf_type=0
1312 
1313  first_dig =MOD(Dtset%bs_coulomb_term,10)
1314  second_dig=Dtset%bs_coulomb_term/10
1315 
1316  Bsp%wtype = second_dig
1317  SELECT CASE (second_dig)
1318  CASE (BSE_WTYPE_NONE)
1319    call wrtout(std_out,"Coulomb term won't be calculated")
1320    Bsp%use_coulomb_term = .FALSE.
1321 
1322  CASE (BSE_WTYPE_FROM_SCR)
1323    call wrtout(std_out,"W is read from an external SCR file")
1324    Bsp%use_coulomb_term = .TRUE.
1325 
1326  CASE (BSE_WTYPE_FROM_MDL)
1327    call wrtout(std_out,"W is approximated with the model dielectric function")
1328    Bsp%use_coulomb_term = .TRUE.
1329    BSp%mdlf_type = MDL_BECHSTEDT
1330    BSp%eps_inf = Dtset%mdf_epsinf
1331    ABI_CHECK(Bsp%eps_inf > zero, "mdf_epsinf <= 0")
1332 
1333  CASE DEFAULT
1334    write(msg,'(a,i0)')" Wrong second digit in bs_coulomb_term: ",Dtset%bs_coulomb_term
1335    ABI_ERROR(msg)
1336  END SELECT
1337  !
1338  ! Diagonal approximation or full matrix?
1339  BSp%use_diagonal_Wgg = .TRUE.
1340  if (Bsp%wtype /= BSE_WTYPE_NONE) then
1341    SELECT CASE (first_dig)
1342    CASE (0)
1343      call wrtout(std_out,"Using diagonal approximation W_GG")
1344      BSp%use_diagonal_Wgg = .TRUE.
1345    CASE (1)
1346      call wrtout(std_out,"Using full W_GG' matrix ")
1347      BSp%use_diagonal_Wgg = .FALSE.
1348    CASE DEFAULT
1349      write(msg,'(a,i0)')" Wrong first digit in bs_coulomb_term: ",Dtset%bs_coulomb_term
1350      ABI_ERROR(msg)
1351    END SELECT
1352  end if
1353 
1354  !TODO move the initialization of the parameters for the interpolation in setup_bse_interp
1355 
1356  BSp%use_interp = .FALSE.
1357  BSp%interp_mode = BSE_INTERP_YG
1358  BSp%interp_kmult(1:3) = 0
1359  BSp%prep_interp = .FALSE.
1360  BSp%sum_overlaps = .TRUE. ! Sum over the overlaps
1361 
1362  ! Printing ncham
1363  BSp%prt_ncham = .FALSE.
1364 
1365  ! Deactivate Interpolation Technique by default
1366 ! if (.FALSE.) then
1367 
1368  ! Reading parameters from the input file
1369  BSp%use_interp = (dtset%bs_interp_mode /= 0)
1370  BSp%prep_interp = (dtset%bs_interp_prep == 1)
1371 
1372  SELECT CASE (dtset%bs_interp_mode)
1373  CASE (0)
1374    ! No interpolation, do not print anything !
1375  CASE (1)
1376    call wrtout(std_out,"Using interpolation technique with energies and wavefunctions from dense WFK")
1377  CASE (2)
1378    call wrtout(std_out,"Interpolation technique with energies and wfn on dense WFK + treatment ABC of divergence")
1379  CASE (3)
1380    call wrtout(std_out,"Interpolation technique + divergence ABC along diagonal")
1381  CASE (4)
1382    call wrtout(std_out,"Using interpolation technique mode 1 with full computation of hamiltonian")
1383  CASE DEFAULT
1384    ABI_ERROR(sjoin("Wrong interpolation mode for bs_interp_mode:", itoa(dtset%bs_interp_mode)))
1385  END SELECT
1386 
1387  ! Read from dtset
1388  if(BSp%use_interp) then
1389    BSp%interp_method = dtset%bs_interp_method
1390    BSp%rl_nb = dtset%bs_interp_rl_nb
1391    BSp%interp_m3_width = dtset%bs_interp_m3_width
1392    BSp%interp_kmult(1:3) = dtset%bs_interp_kmult(1:3)
1393    BSp%interp_mode = dtset%bs_interp_mode
1394  end if
1395 
1396  ! Dimensions and parameters of the calculation.
1397  ! TODO one should add npwx as well
1398  !BSp%npweps=Dtset%npweps
1399  !BSp%npwwfn=Dtset%npwwfn
1400 
1401  ABI_MALLOC(Bsp%lomo_spin, (Bsp%nsppol))
1402  ABI_MALLOC(Bsp%homo_spin, (Bsp%nsppol))
1403  ABI_MALLOC(Bsp%lumo_spin, (Bsp%nsppol))
1404  ABI_MALLOC(Bsp%humo_spin, (Bsp%nsppol))
1405  ABI_MALLOC(Bsp%nbndv_spin, (Bsp%nsppol))
1406  ABI_MALLOC(Bsp%nbndc_spin, (Bsp%nsppol))
1407 
1408  ! FIXME use bs_loband(nsppol)
1409  Bsp%lomo_spin = Dtset%bs_loband
1410  write(std_out,*)"bs_loband",Dtset%bs_loband
1411  !if (Bsp%nsppol == 2) Bsp%lomo_spin(2) = Dtset%bs_loband
1412 
1413  ! Check lomo correct only for unpolarized semiconductors
1414  !if (Dtset%nsppol == 1 .and. Bsp%lomo > Dtset%nelect/2) then
1415  !  write(msg,'(a,i0,a,f8.3)') " Bsp%lomo = ",Bsp%lomo," cannot be greater than nelect/2 = ",Dtset%nelect/2
1416  !  ABI_ERROR(msg)
1417  !end if
1418  !
1419  ! ==============================================
1420  ! ==== Setup of the q for the optical limit ====
1421  ! ==============================================
1422  Bsp%inclvkb = Dtset%inclvkb
1423 
1424  if (Dtset%gw_nqlwl == 0) then
1425    ! Predefined list of 6 q-versors (b vectors and Cart axis)
1426    call cryst%get_redcart_qdirs(Bsp%nq, Bsp%q)
1427  else
1428    BSp%nq = Dtset%gw_nqlwl
1429    ABI_MALLOC(BSp%q, (3,BSp%nq))
1430    BSp%q = Dtset%gw_qlwl
1431    do iq=1,BSp%nq ! normalization
1432      qnorm = normv(BSp%q(:,iq), Cryst%gmet,"G")
1433      BSp%q(:,iq) = BSp%q(:,iq) / qnorm
1434    end do
1435  end if
1436 
1437  ! ======================================================
1438  ! === Define the flags defining the calculation type ===
1439  ! ======================================================
1440  Bsp%calc_type = Dtset%bs_calctype
1441 
1442  BSp%mbpt_sciss = zero ! Shall we use the scissors operator to open the gap?
1443  if (ABS(Dtset%mbpt_sciss)>tol6) BSp%mbpt_sciss = Dtset%mbpt_sciss
1444 
1445 !now test input parameters from input and WFK file and assume some defaults
1446 !
1447 ! TODO Add the possibility of using a randomly shifted k-mesh with nsym>1.
1448 ! so that densities and potentials are correctly symmetrized but
1449 ! the list of the k-point in the IBZ is not expanded.
1450 
1451  if (mband < Dtset%nband(1)) then
1452    write(msg,'(2(a,i0),3a,i0)')&
1453     'WFK file contains only ', mband,' levels instead of ',Dtset%nband(1),' required;',ch10,&
1454     'The calculation will be done with nbands= ',mband
1455    ABI_WARNING(msg)
1456    Dtset%nband(:) = mband
1457  end if
1458 
1459  BSp%nbnds = Dtset%nband(1) ! TODO Note the change in the meaning of input variables
1460 
1461  if (BSp%nbnds<=Dtset%nelect/2) then
1462    write(msg,'(2a,a,i0,a,f8.2)')&
1463     'BSp%nbnds cannot be smaller than homo ',ch10,&
1464     'while BSp%nbnds = ',BSp%nbnds,' and Dtset%nelect = ',Dtset%nelect
1465    ABI_ERROR(msg)
1466  end if
1467 
1468 !TODO add new dim for exchange part and consider the possibility of having npwsigx > npwwfn (see setup_sigma).
1469 
1470  ! === Build enlarged G-sphere for the exchange part ===
1471  call Gsph_c%extend(Cryst, Dtset%ecutwfn, Gsph_x)
1472  call Gsph_x%print(unit=std_out,prtvol=Dtset%prtvol)
1473 
1474  ! NPWVEC as the biggest between npweps and npwwfn. MG RECHECK this part.
1475  !BSp%npwwfn = Dtset%npwwfn
1476  Bsp%npwwfn = Gsph_x%ng  ! FIXME temporary hack
1477  BSp%npwvec=MAX(BSp%npwwfn,BSp%npweps)
1478  Bsp%ecutwfn = Dtset%ecutwfn
1479 
1480  ! Compute Coulomb term on the largest G-sphere.
1481  if (Gsph_x%ng > Gsph_c%ng ) then
1482    call Vcp%init(Gsph_x,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%gw_icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,&
1483      nqlwl,qlwl,comm)
1484  else
1485    call Vcp%init(Gsph_c,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%gw_icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,&
1486      nqlwl,qlwl,comm)
1487  end if
1488 
1489  ABI_FREE(qlwl)
1490 
1491  bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol))
1492  ABI_MALLOC(doccde,(bantot))
1493  ABI_MALLOC(eigen,(bantot))
1494  ABI_MALLOC(occfact,(bantot))
1495  doccde=zero; eigen=zero; occfact=zero
1496 
1497  ! Get occupation from input if occopt == 2
1498  occ_from_dtset = (Dtset%occopt == 2)
1499 
1500  jj=0; ibtot=0
1501  do isppol=1,Dtset%nsppol
1502    do ik_ibz=1,Dtset%nkpt
1503      do ib=1,Hdr_wfk%nband(ik_ibz+(isppol-1)*Dtset%nkpt)
1504        ibtot=ibtot+1
1505        if (ib<=BSP%nbnds) then
1506          jj=jj+1
1507          eigen  (jj)=energies_p(ib,ik_ibz,isppol)
1508          if (occ_from_dtset) then
1509            !Not occupations must be the same for different images
1510            occfact(jj)=Dtset%occ_orig(ibtot,1)
1511          else
1512            occfact(jj)=Hdr_wfk%occ(ibtot)
1513          end if
1514        end if
1515      end do
1516    end do
1517  end do
1518 
1519  ABI_FREE(energies_p)
1520  !
1521  ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of
1522  ! symmetry operations in the old GW code (symmorphy and inversion)
1523  ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6))
1524  ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt')
1525 
1526  ABI_MALLOC(npwarr,(Dtset%nkpt))
1527  npwarr=BSP%npwwfn
1528 
1529  call ebands_init(bantot,KS_BSt,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,&
1530    doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,&
1531    Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,&
1532    dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
1533    dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1534 
1535  ABI_FREE(doccde)
1536  ABI_FREE(eigen)
1537  ABI_FREE(npwarr)
1538 
1539  !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when
1540  ! the occupation scheme for semiconductors is used.
1541  call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=Dtset%prtvol)
1542 
1543  call ebands_print(KS_BSt,"Band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol)
1544 
1545  call ebands_report_gap(KS_BSt,header=" KS band structure",unit=std_out,mode_paral="COLL")
1546 
1547  ABI_MALLOC(val_indeces,(KS_BSt%nkpt,KS_BSt%nsppol))
1548  val_indeces = ebands_get_valence_idx(KS_BSt)
1549 
1550  do spin=1,KS_BSt%nsppol
1551    val_idx(spin) = val_indeces(1,spin)
1552    write(msg,'(a,i2,a,i0)')" For spin : ",spin," val_idx ",val_idx(spin)
1553    call wrtout(std_out,msg)
1554    if ( ANY(val_indeces(1,spin) /= val_indeces(:,spin)) ) then
1555      ABI_ERROR("BSE code does not support metals")
1556    end if
1557  end do
1558 
1559  ABI_FREE(val_indeces)
1560  !
1561  ! === Create the BSE header ===
1562  call hdr_init(KS_BSt,codvsn,Dtset,Hdr_bse,Pawtab,pertcase0,Psps,wvl)
1563 
1564  ! === Get Pawrhoij from the header of the WFK file ===
1565  ABI_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw))
1566  if (Dtset%usepaw==1) then
1567    call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab)
1568    call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij)
1569  end if
1570 
1571  ! CP modified
1572  !call hdr_bse%update(bantot,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
1573  call hdr_bse%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
1574  ! End CP modified
1575 
1576  ABI_FREE(occfact)
1577 
1578  if (Dtset%usepaw==1) call pawrhoij_free(Pawrhoij)
1579  ABI_FREE(Pawrhoij)
1580 
1581  ! Find optimal value for G-sphere enlargment due to oscillator matrix elements
1582  ! We will split k-points over processors
1583  call xmpi_split_work(Kmesh%nbz, comm, my_k1, my_k2)
1584 
1585  ! If there is no work to do, just skip the computation
1586  if (my_k2-my_k1+1 <= 0) then
1587    ng0sh_opt(:)=(/zero,zero,zero/)
1588  else
1589    ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true since bz is buggy
1590    ! * -one is used because we loop over all the possibile differences, unlike screening
1591    call get_ng0sh(my_k2-my_k1+1,Kmesh%bz(:,my_k1:my_k2),Kmesh%nbz,Kmesh%bz,&
1592 &    Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt)
1593  end if
1594 
1595  call xmpi_max(ng0sh_opt,BSp%mg0,comm,ierr)
1596 
1597  write(msg,'(a,3(i0,1x))') ' optimal value for ng0sh = ',BSp%mg0
1598  call wrtout(std_out,msg)
1599 
1600  ! === Setup of the FFT mesh for the oscilator strengths ===
1601  ! * ngfft_osc(7:18)==Dtset%ngfft(7:18) which is initialized before entering screening.
1602  ! * Here we redefine ngfft_osc(1:6) according to the following options :
1603  !
1604  ! method==0 --> FFT grid read from fft.in (debugging purpose)
1605  ! method==1 --> Normal FFT mesh
1606  ! method==2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
1607  ! method==3 --> Doubled FFT grid, same as the the FFT for the density,
1608  !
1609  ! enforce_sym==1 ==> Enforce a FFT mesh compatible with all the symmetry operation and FFT library
1610  ! enforce_sym==0 ==> Find the smallest FFT grid compatbile with the library, do not care about symmetries
1611  !
1612  ngfft_osc(1:18)=Dtset%ngfft(1:18); method=2
1613  if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0
1614  if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1
1615  if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2
1616  if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3
1617  enforce_sym=MOD(Dtset%fftgw,10)
1618 
1619  call setmesh(gmet,Gsph_x%gvec,ngfft_osc,BSp%npwvec,BSp%npweps,BSp%npwwfn,nfftot_osc,method,BSp%mg0,Cryst,enforce_sym)
1620  nfftot_osc=PRODUCT(ngfft_osc(1:3))
1621 
1622  call print_ngfft(ngfft_osc,"FFT mesh for oscillator matrix elements",std_out,"COLL",prtvol=Dtset%prtvol)
1623  !
1624  ! BSp%homo gives the
1625  !BSp%homo  = val_idx(1)
1626  ! highest occupied band for each spin
1627  BSp%homo_spin = val_idx
1628 
1629  ! TODO generalize the code to account for this unlikely case.
1630  !if (Dtset%nsppol==2) then
1631  !  ABI_CHECK(BSp%homo == val_idx(2),"Different valence indeces for spin up and down")
1632  !end if
1633 
1634  !BSp%lumo = BSp%homo + 1
1635  !BSp%humo = BSp%nbnds
1636  !BSp%nbndv = BSp%homo  - BSp%lomo + 1
1637  !BSp%nbndc = BSp%nbnds - BSp%homo
1638 
1639  BSp%lumo_spin = BSp%homo_spin + 1
1640  BSp%humo_spin = BSp%nbnds
1641  BSp%nbndv_spin = BSp%homo_spin  - BSp%lomo_spin + 1
1642  BSp%nbndc_spin = BSp%nbnds - BSp%homo_spin
1643  BSp%maxnbndv = MAXVAL(BSp%nbndv_spin(:))
1644  BSp%maxnbndc = MAXVAL(BSp%nbndc_spin(:))
1645 
1646  BSp%nkbz = Kmesh%nbz
1647 
1648  call ebands_copy(KS_BSt,QP_bst)
1649  ABI_MALLOC(igwene,(QP_bst%mband,QP_bst%nkpt,QP_bst%nsppol))
1650  igwene=zero
1651 
1652  call bsp_calctype2str(Bsp,msg)
1653  call wrtout(std_out,"Calculation type: "//TRIM(msg))
1654 
1655  SELECT CASE (Bsp%calc_type)
1656  CASE (BSE_HTYPE_RPA_KS)
1657    if (ABS(BSp%mbpt_sciss)>tol6) then
1658      write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies."
1659      call wrtout(std_out,msg)
1660      call ebands_apply_scissors(QP_BSt,BSp%mbpt_sciss)
1661    else
1662      write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]."
1663      call wrtout(std_out,msg)
1664    end if
1665 
1666  CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw
1667    gw_fname=TRIM(Dtfil%filnam_ds(4))//'_GW'
1668    gw_fname="__in.gw__"
1669    if (.not.file_exists(gw_fname)) then
1670      msg = " File "//TRIM(gw_fname)//" not found. Aborting now"
1671      ABI_ERROR(msg)
1672    end if
1673 
1674    call rdgw(QP_Bst,gw_fname,igwene,extrapolate=.FALSE.) ! here gwenergy is real
1675 
1676    do isppol=1,Dtset%nsppol
1677      write(std_out,*) ' k       GW energies [eV]'
1678      do ik_ibz=1,Kmesh%nibz
1679        write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(QP_bst%eig(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds)
1680      end do
1681      write(std_out,*) ' k       Im GW energies [eV]'
1682      do ik_ibz=1,Kmesh%nibz
1683        write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(igwene(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds)
1684      end do
1685    end do
1686    !
1687    ! If required apply the scissors operator on top of the QP bands structure (!)
1688    if (ABS(BSp%mbpt_sciss)>tol6) then
1689      write(msg,'(a,f8.2,a)')' Applying a scissors operator ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the QP energies!"
1690      ABI_COMMENT(msg)
1691      call ebands_apply_scissors(QP_BSt,BSp%mbpt_sciss)
1692    end if
1693 
1694  CASE (BSE_HTYPE_RPA_QP)
1695    ABI_ERROR("Not implemented error!")
1696 
1697  CASE DEFAULT
1698    write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type
1699    ABI_ERROR(msg)
1700  END SELECT
1701 
1702  call ebands_report_gap(QP_BSt,header=" QP band structure",unit=std_out,mode_paral="COLL")
1703 
1704  ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index.
1705  ! FIXME: linewidths not coded.
1706  ABI_MALLOC(gw_energy,(BSp%nbnds,Kmesh%nibz,Dtset%nsppol))
1707  gw_energy = QP_BSt%eig
1708 
1709  BSp%have_complex_ene = ANY(igwene > tol16)
1710 
1711  ! Compute the number of resonant transitions, nreh, for the two spin channels and initialize BSp%Trans.
1712  ABI_MALLOC(Bsp%nreh,(Bsp%nsppol))
1713 
1714  ! Possible cutoff on the transitions.
1715  BSp%ircut = Dtset%bs_eh_cutoff(1)
1716  BSp%uvcut = Dtset%bs_eh_cutoff(2)
1717 
1718  call init_transitions(BSp%Trans,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz,Bsp%nbnds,Bsp%nkibz,&
1719 &  BSp%nsppol,Dtset%nspinor,gw_energy,QP_BSt%occ,Kmesh%tab,minmax_tene,Bsp%nreh)
1720 
1721  ! Setup of the frequency mesh for the absorption spectrum.
1722  ! If not specified, use the min-max resonant transition energy and make it 10% smaller|larger.
1723 
1724  !if (ABS(Dtset%bs_freq_mesh(1)) < tol6) then
1725  !   Dtset%bs_freq_mesh(1) = MAX(minmax_tene(1) - minmax_tene(1) * 0.1, zero)
1726  !end if
1727 
1728  if (ABS(Dtset%bs_freq_mesh(2)) < tol6) then
1729     Dtset%bs_freq_mesh(2) = minmax_tene(2) + minmax_tene(2) * 0.1
1730  end if
1731 
1732  Bsp%omegai = Dtset%bs_freq_mesh(1)
1733  Bsp%omegae = Dtset%bs_freq_mesh(2)
1734  Bsp%domega = Dtset%bs_freq_mesh(3)
1735  BSp%broad  = Dtset%zcut
1736 
1737  ! The frequency mesh (including the complex imaginary shift)
1738  BSp%nomega = (BSp%omegae - BSp%omegai)/BSp%domega + 1
1739  ABI_MALLOC(BSp%omega,(BSp%nomega))
1740  do io=1,BSp%nomega
1741    BSp%omega(io) = (BSp%omegai + (io-1)*BSp%domega)  + j_dpc*BSp%broad
1742  end do
1743 
1744  ABI_FREE(gw_energy)
1745  ABI_FREE(igwene)
1746 
1747  do spin=1,Bsp%nsppol
1748    write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh(spin)
1749    call wrtout(std_out,msg)
1750  end do
1751 
1752  if (ANY(Bsp%nreh/=Bsp%nreh(1))) then
1753    write(msg,'(a,2(i0,1x))')" BSE code with different number of transitions for the two spin channels: ",Bsp%nreh
1754    ABI_WARNING(msg)
1755  end if
1756  !
1757  ! Create transition table vcks2t
1758  Bsp%lomo_min = MINVAL(BSp%lomo_spin)
1759  Bsp%homo_max = MAXVAL(BSp%homo_spin)
1760  Bsp%lumo_min = MINVAL(BSp%lumo_spin)
1761  Bsp%humo_max = MAXVAL(BSp%humo_spin)
1762 
1763  ABI_MALLOC(Bsp%vcks2t,(BSp%lomo_min:BSp%homo_max,BSp%lumo_min:BSp%humo_max,BSp%nkbz,Dtset%nsppol))
1764  Bsp%vcks2t = 0
1765 
1766  do spin=1,BSp%nsppol
1767    do it=1,BSp%nreh(spin)
1768      BSp%vcks2t(BSp%Trans(it,spin)%v,BSp%Trans(it,spin)%c,BSp%Trans(it,spin)%k,spin) = it
1769    end do
1770  end do
1771 
1772  hexc_size = SUM(Bsp%nreh); if (Bsp%use_coupling>0) hexc_size=2*hexc_size
1773  if (Bsp%nstates<=0) then
1774    Bsp%nstates=hexc_size
1775  else
1776    if (Bsp%nstates>hexc_size) then
1777       Bsp%nstates=hexc_size
1778       write(msg,'(2(a,i0),2a)')&
1779        "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates,ch10,&
1780        "the number of excitonic states nstates has been modified"
1781      ABI_WARNING(msg)
1782    end if
1783  end if
1784 
1785  msg=' Fundamental parameters for the solution of the Bethe-Salpeter equation:'
1786  call print_bs_parameters(BSp,unit=std_out,header=msg,mode_paral="COLL",prtvol=Dtset%prtvol)
1787  call print_bs_parameters(BSp,unit=ab_out, header=msg,mode_paral="COLL")
1788 
1789  if (ANY (Cryst%symrec(:,:,1) /= RESHAPE ( (/1,0,0,0,1,0,0,0,1/),(/3,3/) )) .or. ANY( ABS(Cryst%tnons(:,1)) > tol6) ) then
1790    write(msg,'(3a,9i2,2a,3f6.3,2a)')&
1791      "The first symmetry operation should be the Identity with zero tnons while ",ch10,&
1792      "symrec(:,:,1) = ",Cryst%symrec(:,:,1),ch10,&
1793      "tnons(:,1)    = ",Cryst%tnons(:,1),ch10,&
1794      "This is not allowed, sym_rhotwgq0 should be changed."
1795    ABI_ERROR(msg)
1796  end if
1797  !
1798  ! Prefix for generic output files.
1799  BS_files%out_basename = TRIM(Dtfil%filnam_ds(4))
1800  !
1801  ! Search for files to restart from.
1802  if (Dtset%gethaydock/=0 .or. Dtset%irdhaydock/=0) then
1803    BS_files%in_haydock_basename = TRIM(Dtfil%fnameabi_haydock)
1804  end if
1805 
1806  test_file = Dtfil%fnameabi_bsham_reso
1807  if (file_exists(test_file)) then
1808    BS_files%in_hreso = test_file
1809  else
1810    BS_files%out_hreso = TRIM(Dtfil%filnam_ds(4))//'_BSR'
1811  end if
1812 
1813  test_file = Dtfil%fnameabi_bsham_coup
1814  if (file_exists(test_file) ) then
1815    BS_files%in_hcoup = test_file
1816  else
1817    BS_files%out_hcoup = TRIM(Dtfil%filnam_ds(4))//'_BSC'
1818  end if
1819  !
1820  ! in_eig is the name of the input file with eigenvalues and eigenvectors
1821  ! constructed from getbseig or irdbseig. out_eig is the name of the output file
1822  ! produced by this dataset. in_eig_exists checks for the presence of the input file.
1823  !
1824  if (file_exists(Dtfil%fnameabi_bseig)) then
1825    BS_files%in_eig = Dtfil%fnameabi_bseig
1826  else
1827    BS_files%out_eig = TRIM(BS_files%out_basename)//"_BSEIG"
1828  end if
1829 
1830  call print_bs_files(BS_files,unit=std_out)
1831  !
1832  ! ==========================================================
1833  ! ==== Temperature dependence of the spectrum ==============
1834  ! ==========================================================
1835  BSp%do_ep_renorm = .FALSE.
1836  BSp%do_lifetime = .FALSE. ! Not yet implemented
1837 
1838  ep_nc_fname = 'test_EP.nc'
1839  if(file_exists(ep_nc_fname)) then
1840    BSp%do_ep_renorm = .TRUE.
1841 
1842    if(my_rank == master) then
1843      call eprenorms_from_epnc(Epren,ep_nc_fname)
1844    end if
1845    call eprenorms_bcast(Epren,master,comm)
1846  end if
1847  !
1848  ! ==========================================================
1849  ! ==== Final check on the parameters of the calculation ====
1850  ! ==========================================================
1851  if ( Bsp%use_coupling>0 .and. ALL(Bsp%algorithm /= [BSE_ALGO_DDIAGO, BSE_ALGO_HAYDOCK]) ) then
1852    ABI_ERROR("Resonant+Coupling is only available with the direct diagonalization or the haydock method.")
1853  end if
1854 
1855  ! autoparal section
1856  if (dtset%max_ncpus /=0 .and. dtset%autoparal /=0 ) then
1857    ount = ab_out
1858    ! TODO:
1859    ! nsppol and calculation with coupling!
1860 
1861    ! Temporary table needed to estimate memory
1862    ABI_MALLOC(nlmn_atm,(Cryst%natom))
1863    if (Dtset%usepaw==1) then
1864      do iat=1,Cryst%natom
1865        nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size
1866      end do
1867    end if
1868 
1869    tot_nreh = SUM(BSp%nreh)
1870    work_size = tot_nreh * (tot_nreh + 1) / 2
1871 
1872    write(ount,'(a)')"--- !Autoparal"
1873    write(ount,"(a)")'#Autoparal section for Bethe-Salpeter runs.'
1874 
1875    write(ount,"(a)")   "info:"
1876    write(ount,"(a,i0)")"    autoparal: ",dtset%autoparal
1877    write(ount,"(a,i0)")"    max_ncpus: ",dtset%max_ncpus
1878    write(ount,"(a,i0)")"    nkibz: ",Bsp%nkibz
1879    write(ount,"(a,i0)")"    nkbz: ",Bsp%nkbz
1880    write(ount,"(a,i0)")"    nsppol: ",dtset%nsppol
1881    write(ount,"(a,i0)")"    nspinor: ",dtset%nspinor
1882    write(ount,"(a,i0)")"    lomo_min: ",Bsp%lomo_min
1883    write(ount,"(a,i0)")"    humo_max: ",Bsp%humo_max
1884    write(ount,"(a,i0)")"    tot_nreh: ",tot_nreh
1885    !write(ount,"(a,i0)")"    nbnds: ",Ep%nbnds
1886 
1887    ! Wavefunctions are not distributed. We read all the bands
1888    ! from 1 up to Bsp%nbnds because we have to recompute rhor
1889    ! but then we deallocate all the states that are not used for the construction of the e-h
1890    ! before allocating the EXC hamiltonian. Hence we can safely use  (humo - lomo + 1) instead of Bsp%nbnds.
1891    !my_nbks = (Bsp%humo - Bsp%lomo +1) * Bsp%nkibz * Dtset%nsppol
1892 
1893    ! This one overestimates the memory but it seems to be safer.
1894    my_nbks = Bsp%nbnds * Dtset%nkpt * Dtset%nsppol
1895 
1896    ! Memory needed for Fourier components ug.
1897    ug_mem = two*gwpc*Dtset%nspinor*Bsp%npwwfn*my_nbks*b2Mb
1898 
1899    ! Memory needed for real space ur.
1900    ur_mem = zero
1901    if (MODULO(Dtset%gwmem,10)==1) then
1902      ur_mem = two*gwpc*Dtset%nspinor*nfftot_osc*my_nbks*b2Mb
1903    end if
1904 
1905    ! Memory needed for PAW projections Cprj
1906    cprj_mem = zero
1907    if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb
1908 
1909    wfsmem_mb = ug_mem + ur_mem + cprj_mem
1910 
1911    ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI:  wavefunctions + W
1912    nonscal_mem = (wfsmem_mb + two*gwpc*BSp%npweps**2*b2Mb) * 1.1_dp
1913 
1914    ! List of configurations.
1915    write(ount,"(a)")"configurations:"
1916    do il=1,dtset%max_ncpus
1917      if (il > work_size) cycle
1918      neh_per_proc = work_size / il
1919      neh_per_proc = neh_per_proc + MOD(work_size, il)
1920      eff = (one * work_size) / (il * neh_per_proc)
1921 
1922      ! EXC matrix is distributed.
1923      mempercpu_mb = nonscal_mem + two * dpc * neh_per_proc * b2Mb
1924 
1925      write(ount,"(a,i0)")"    - tot_ncpus: ",il
1926      write(ount,"(a,i0)")"      mpi_ncpus: ",il
1927      !write(ount,"(a,i0)")"      omp_ncpus: ",omp_ncpus
1928      write(ount,"(a,f12.9)")"      efficiency: ",eff
1929      write(ount,"(a,f12.2)")"      mem_per_cpu: ",mempercpu_mb
1930    end do
1931 
1932    write(ount,'(a)')"..."
1933 
1934    ABI_FREE(nlmn_atm)
1935    ABI_ERROR_NODUMP("aborting now")
1936  end if
1937 
1938  DBG_EXIT("COLL")
1939 
1940 end subroutine setup_bse

m_bethe_salpeter/setup_bse_interp [ Functions ]

[ Top ] [ m_bethe_salpeter ] [ Functions ]

NAME

  setup_bse_interp

FUNCTION

INPUTS

 ngfft_gw(18)=Information about 3D FFT for density and potentials, see ~abinit/doc/variables/vargs.htm#ngfft
 acell(3)=Length scales of primitive translations (bohr)
 rprim(3,3)=Dimensionless real space primitive translations.
 Dtset<dataset_type>=All input variables for this dataset.
  Some of them might be redefined here TODO
 Dtfil=filenames and unit numbers used in abinit. fnameabi_wfkfile is changed is Fortran file is not
 found but a netcdf version with similar name is available.

OUTPUT

 Cryst<crystal_structure>=Info on the crystalline Structure.
 Kmesh<BZ_mesh_type>=Structure defining the k-sampling for the wavefunctions.
 Qmesh<BZ_mesh_type>=Structure defining the q-sampling for the symmetrized inverse dielectric matrix.
 Gsph_x<gsphere_t=Data type gathering info on the G-sphere for wave functions and e^{-1},
 KS_BSt<Bandstructure_type>=The KS band structure (energies, occupancies, k-weights...)
 Vcp<vcoul_t>=Structure gathering information on the Coulomb interaction in reciprocal space,
   including a possible cutoff in real space.
 ngfft_osc(18)=Contain all needed information about the 3D FFT for the oscillator matrix elements.
   See ~abinit/doc/variables/vargs.htm#ngfft
 Bsp<excparam>=Basic parameters defining the Bethe-Salpeter run. Completely initialed in output.
 Hdr_wfk<Hdr_type>=The header of the WFK file.
 Hdr_bse<Hdr_type>=Local header initialized from the parameters used for the Bethe-Salpeter calculation.
 w_file=File name used to construct W. Set to ABI_NOFILE if no external file is used.

SOURCE

1975 subroutine setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh, &
1976    Kmesh_dense,Qmesh_dense,KS_BSt_dense,QP_bst_dense,Gsph_x,Gsph_c,Vcp_dense,Hdr_wfk_dense,grid,comm)
1977 
1978 !Arguments ------------------------------------
1979 !scalars
1980  integer,intent(in) :: comm
1981  type(dataset_type),intent(in) :: Dtset
1982  type(datafiles_type),intent(inout) :: Dtfil
1983  type(excparam),intent(inout) :: Bsp
1984  type(hdr_type),intent(out) :: Hdr_wfk_dense
1985  type(crystal_t),intent(in) :: Cryst
1986  type(kmesh_t),intent(in) :: Kmesh
1987  type(kmesh_t),intent(out) :: Kmesh_dense,Qmesh_dense
1988  type(ebands_t),intent(out) :: KS_BSt_dense,QP_Bst_dense
1989  type(double_grid_t),intent(out) :: grid
1990  type(vcoul_t),intent(out) :: Vcp_dense
1991  type(gsphere_t),intent(out) :: Gsph_x,Gsph_c
1992 !arrays
1993 
1994 !Local variables ------------------------------
1995 !scalars
1996  integer,parameter :: pertcase0=0,master=0
1997  integer :: bantot_dense,ib,ibtot,ik_ibz,isppol,jj
1998  integer :: nbnds_kss_dense
1999  integer :: spin,hexc_size
2000  integer :: my_rank
2001  integer :: it
2002  integer :: nprocs
2003  integer :: is1,is2,is3,is4
2004  real(dp) :: nelect_hdr_dense
2005  logical,parameter :: remove_inv=.FALSE.
2006  character(len=500) :: msg
2007  character(len=fnlen) :: wfk_fname_dense
2008  integer :: nqlwl
2009 !arrays
2010  integer,allocatable :: npwarr(:)
2011  real(dp),allocatable :: shiftk(:,:)
2012  real(dp),allocatable :: doccde(:),eigen(:),occfact(:)
2013  real(dp),pointer :: energies_p_dense(:,:,:)
2014  complex(dpc),allocatable :: gw_energy(:,:,:)
2015  integer,allocatable :: nbands_temp(:)
2016  integer :: kptrlatt_dense(3,3)
2017  real(dp),allocatable :: qlwl(:,:)
2018  real(dp) :: minmax_tene(2)
2019 
2020 !************************************************************************
2021 
2022  DBG_ENTER("COLL")
2023 
2024  kptrlatt_dense = zero
2025 
2026  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
2027 
2028  SELECT CASE(BSp%interp_mode)
2029  CASE (1,2,3,4)
2030    nbnds_kss_dense = -1
2031    wfk_fname_dense = Dtfil%fnameabi_wfkfine
2032    call wrtout(std_out,"BSE Interpolation: will read energies from: "//trim(wfk_fname_dense),"COLL")
2033 
2034    if (nctk_try_fort_or_ncfile(wfk_fname_dense, msg) /= 0) then
2035      ABI_ERROR(msg)
2036    end if
2037 
2038    Dtfil%fnameabi_wfkfine = wfk_fname_dense
2039 
2040    call wfk_read_eigenvalues(wfk_fname_dense,energies_p_dense,Hdr_wfk_dense,comm)
2041    nbnds_kss_dense = MAXVAL(Hdr_wfk_dense%nband)
2042  CASE DEFAULT
2043    ABI_ERROR("Not yet implemented")
2044  END SELECT
2045 
2046  nelect_hdr_dense = Hdr_wfk_dense%nelect
2047 
2048  if (ABS(Dtset%nelect-nelect_hdr_dense)>tol6) then
2049    write(msg,'(2(a,f8.2))')&
2050 &   "File contains ", nelect_hdr_dense," electrons but nelect initialized from input is ",Dtset%nelect
2051    ABI_ERROR(msg)
2052  end if
2053 
2054  ! Setup of the k-point list and symmetry tables in the  BZ
2055  SELECT CASE(BSp%interp_mode)
2056  CASE (1,2,3,4)
2057    if(Dtset%chksymbreak == 0) then
2058      ABI_MALLOC(shiftk,(3,Dtset%nshiftk))
2059      kptrlatt_dense(:,1) = BSp%interp_kmult(1)*Dtset%kptrlatt(:,1)
2060      kptrlatt_dense(:,2) = BSp%interp_kmult(2)*Dtset%kptrlatt(:,2)
2061      kptrlatt_dense(:,3) = BSp%interp_kmult(3)*Dtset%kptrlatt(:,3)
2062      do jj = 1,Dtset%nshiftk
2063        shiftk(:,jj) = Bsp%interp_kmult(:)*Dtset%shiftk(:,jj)
2064      end do
2065      call make_mesh(Kmesh_dense,Cryst,Dtset%kptopt,kptrlatt_dense,Dtset%nshiftk,shiftk,break_symmetry=.TRUE.)
2066      ABI_FREE(shiftk)
2067    else
2068      !Initialize Kmesh with no wrapping inside ]-0.5;0.5]
2069      call Kmesh_dense%init(Cryst,Hdr_wfk_dense%nkpt,Hdr_wfk_dense%kptns,Dtset%kptopt)
2070    end if
2071  CASE DEFAULT
2072    ABI_ERROR("Not yet implemented")
2073  END SELECT
2074 
2075  ! Init Qmesh
2076  call find_qmesh(Qmesh_dense,Cryst,Kmesh_dense)
2077 
2078  call Gsph_c%init(Cryst, 0, ecut=Dtset%ecuteps)
2079 
2080  call double_grid_init(Kmesh,Kmesh_dense,Dtset%kptrlatt,BSp%interp_kmult,grid)
2081 
2082  BSp%nkibz_interp = Kmesh_dense%nibz  !We might allow for a smaller number of points....
2083 
2084  call Kmesh_dense%print("Interpolated K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL")
2085  call Kmesh_dense%print("Interpolated K-mesh for the wavefunctions",ab_out, 0,           "COLL")
2086 
2087  if (nbnds_kss_dense < Dtset%nband(1)) then
2088    write(msg,'(2(a,i0),3a,i0)')&
2089     'Interpolated WFK file contains only ', nbnds_kss_dense,' levels instead of ',Dtset%nband(1),' required;',ch10,&
2090     'The calculation will be done with nbands= ',nbnds_kss_dense
2091    ABI_WARNING(msg)
2092    ABI_ERROR("Not supported yet !")
2093  end if
2094 
2095  ABI_MALLOC(nbands_temp,(Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol))
2096  do isppol=1,Hdr_wfk_dense%nsppol
2097    do ik_ibz=1,Hdr_wfk_dense%nkpt
2098      nbands_temp(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt) = Dtset%nband(1)
2099    end do
2100  end do
2101 
2102  call Gsph_c%extend(Cryst, Dtset%ecutwfn, Gsph_x)
2103  call Gsph_x%print(unit=std_out,prtvol=Dtset%prtvol)
2104 
2105  nqlwl=1
2106  ABI_MALLOC(qlwl,(3,nqlwl))
2107  qlwl(:,nqlwl)= GW_Q0_DEFAULT
2108 
2109  ! Compute Coulomb term on the largest G-sphere.
2110  if (Gsph_x%ng > Gsph_c%ng ) then
2111    call Vcp_dense%init(Gsph_x,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%gw_icutcoul,&
2112                        Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm)
2113  else
2114    call Vcp_dense%init(Gsph_c,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%gw_icutcoul,&
2115                        Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm)
2116  end if
2117 
2118  ABI_FREE(qlwl)
2119 
2120  bantot_dense=SUM(Hdr_wfk_dense%nband(1:Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol))
2121  ABI_MALLOC(doccde,(bantot_dense))
2122  ABI_MALLOC(eigen,(bantot_dense))
2123  ABI_MALLOC(occfact,(bantot_dense))
2124  doccde=zero; eigen=zero; occfact=zero
2125 
2126  jj=0; ibtot=0
2127  do isppol=1,Hdr_wfk_dense%nsppol
2128    do ik_ibz=1,Hdr_wfk_dense%nkpt
2129      do ib=1,Hdr_wfk_dense%nband(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt)
2130        ibtot=ibtot+1
2131        if (ib<=BSP%nbnds) then
2132          jj=jj+1
2133          occfact(jj)=Hdr_wfk_dense%occ(ibtot)
2134          eigen  (jj)=energies_p_dense(ib,ik_ibz,isppol)
2135        end if
2136      end do
2137    end do
2138  end do
2139 
2140  ABI_FREE(energies_p_dense)
2141 
2142  ABI_MALLOC(npwarr,(kmesh_dense%nibz))
2143  npwarr=BSP%npwwfn
2144 
2145  ! CP modified
2146  !call ebands_init(bantot_dense,KS_BSt_dense,Dtset%nelect,doccde,eigen,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,nbands_temp,&
2147 !&  Kmesh_dense%nibz,npwarr,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,Hdr_wfk_dense%tphysel,Hdr_wfk_dense%tsmear,&
2148 !&  Hdr_wfk_dense%occopt,occfact,Kmesh_dense%wt,&
2149 !&  hdr_wfk_dense%cellcharge, hdr_wfk_dense%kptopt, hdr_wfk_dense%kptrlatt_orig, hdr_wfk_dense%nshiftk_orig, &
2150 !&  hdr_wfk_dense%shiftk_orig, hdr_wfk_dense%kptrlatt, hdr_wfk_dense%nshiftk, hdr_wfk_dense%shiftk)
2151  call ebands_init(bantot_dense,KS_BSt_dense,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,& ! CP added
2152 &  doccde,eigen,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,nbands_temp,&
2153 &  Kmesh_dense%nibz,npwarr,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,Hdr_wfk_dense%tphysel,Hdr_wfk_dense%tsmear,&
2154 &  Hdr_wfk_dense%occopt,occfact,Kmesh_dense%wt,&
2155 &  hdr_wfk_dense%cellcharge, hdr_wfk_dense%kptopt, hdr_wfk_dense%kptrlatt_orig, hdr_wfk_dense%nshiftk_orig, &
2156 &  hdr_wfk_dense%shiftk_orig, hdr_wfk_dense%kptrlatt, hdr_wfk_dense%nshiftk, hdr_wfk_dense%shiftk)
2157  ! End CP modified
2158 
2159  ABI_FREE(doccde)
2160  ABI_FREE(eigen)
2161  ABI_FREE(npwarr)
2162 
2163  ABI_FREE(nbands_temp)
2164 
2165  ABI_FREE(occfact)
2166 
2167  !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when
2168  ! the occupation scheme for semiconductors is used.
2169  call ebands_update_occ(KS_BSt_dense,Dtset%spinmagntarget,prtvol=Dtset%prtvol)
2170 
2171  call ebands_print(KS_BSt_dense,"Interpolated band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol)
2172 
2173  call ebands_report_gap(KS_BSt_dense,header="Interpolated KS band structure",unit=std_out,mode_paral="COLL")
2174 
2175  BSp%nkbz_interp = Kmesh_dense%nbz
2176 
2177  call ebands_copy(KS_BSt_dense,QP_bst_dense)
2178 
2179  SELECT CASE (Bsp%calc_type)
2180  CASE (BSE_HTYPE_RPA_KS)
2181    if (ABS(BSp%mbpt_sciss)>tol6) then
2182      write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies."
2183      call wrtout(std_out,msg)
2184      call ebands_apply_scissors(QP_BSt_dense,BSp%mbpt_sciss)
2185    else
2186      write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]."
2187      call wrtout(std_out,msg)
2188    end if
2189    !
2190  CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw
2191    ABI_ERROR("Not yet implemented with interpolation !")
2192  CASE (BSE_HTYPE_RPA_QP)
2193    ABI_ERROR("Not implemented error!")
2194  CASE DEFAULT
2195    write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type
2196    ABI_ERROR(msg)
2197  END SELECT
2198 
2199  call ebands_report_gap(QP_BSt_dense,header=" Interpolated QP band structure",unit=std_out,mode_paral="COLL")
2200 
2201  ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index.
2202  ! FIXME: linewidths not coded.
2203  ABI_MALLOC(gw_energy, (BSp%nbnds,Kmesh_dense%nibz,Dtset%nsppol))
2204  gw_energy = QP_BSt_dense%eig
2205 
2206  ABI_MALLOC(Bsp%nreh_interp,(Hdr_wfk_dense%nsppol))
2207  Bsp%nreh_interp=zero
2208 
2209  call init_transitions(BSp%Trans_interp,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz_interp,Bsp%nbnds,&
2210 &  Bsp%nkibz_interp,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,gw_energy,QP_BSt_dense%occ,Kmesh_dense%tab,minmax_tene,&
2211 &  Bsp%nreh_interp)
2212 
2213  ABI_FREE(gw_energy)
2214 
2215  do spin=1,Dtset%nsppol
2216    write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh_interp(spin)
2217    call wrtout(std_out,msg)
2218  end do
2219 
2220  if (ANY(Bsp%nreh_interp/=Bsp%nreh_interp(1))) then
2221    write(msg,'(a,(i0))')" BSE code does not support different number of transitions for the two spin channels",Bsp%nreh
2222    ABI_ERROR(msg)
2223  end if
2224  !
2225  ! Create transition table vcks2t
2226  is1=BSp%lomo_min;is2=BSp%homo_max;is3=BSp%lumo_min;is4=BSp%humo_max
2227  ABI_MALLOC(Bsp%vcks2t_interp, (is1:is2,is3:is4,BSp%nkbz_interp,Dtset%nsppol))
2228  Bsp%vcks2t_interp = 0
2229 
2230  do spin=1,Dtset%nsppol
2231    do it=1,BSp%nreh_interp(spin)
2232      BSp%vcks2t_interp(BSp%Trans_interp(it,spin)%v,BSp%Trans_interp(it,spin)%c, BSp%Trans_interp(it,spin)%k,spin) = it
2233    end do
2234  end do
2235 
2236  hexc_size = SUM(Bsp%nreh_interp); if (Bsp%use_coupling>0) hexc_size=2*hexc_size
2237  if (Bsp%nstates_interp<=0) then
2238    Bsp%nstates_interp=hexc_size
2239  else
2240    if (Bsp%nstates_interp>hexc_size) then
2241       Bsp%nstates_interp=hexc_size
2242       write(msg,'(2(a,i0),2a)')&
2243        "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates_interp,ch10,&
2244        "the number of excitonic states nstates has been modified"
2245      ABI_WARNING(msg)
2246    end if
2247  end if
2248 
2249  DBG_EXIT("COLL")
2250 
2251 end subroutine setup_bse_interp