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-2018 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 .

PARENTS

CHILDREN

SOURCE

 23 #if defined HAVE_CONFIG_H
 24 #include "config.h"
 25 #endif
 26 
 27 #include "abi_common.h"
 28 
 29 module m_bethe_salpeter
 30 
 31  use defs_basis
 32  use defs_datatypes
 33  use defs_abitypes
 34  use defs_wvltypes
 35  use m_bs_defs
 36  use m_abicore
 37  use m_xmpi
 38  use m_errors
 39  use m_screen
 40  use m_nctk
 41 #ifdef HAVE_NETCDF
 42  use netcdf
 43 #endif
 44  use m_hdr
 45 
 46  use m_gwdefs,          only : GW_Q0_DEFAULT
 47  use m_time,            only : timab
 48  use m_fstrings,        only : strcat, sjoin, endswith
 49  use m_io_tools,        only : file_exists, iomode_from_fname
 50  use m_geometry,        only : mkrdim, metric, normv
 51  use m_hide_lapack,     only : matrginv
 52  use m_mpinfo,          only : destroy_mpi_enreg, initmpi_seq
 53  use m_fftcore,         only : print_ngfft
 54  use m_fft_mesh,        only : rotate_FFT_mesh, get_gftt, setmesh
 55  use m_fft,             only : fourdp
 56  use m_crystal,         only : crystal_t, crystal_free, crystal_print, idx_spatial_inversion
 57  use m_crystal_io,      only : crystal_ncwrite, crystal_from_hdr
 58  use m_bz_mesh,         only : kmesh_t, kmesh_init, kmesh_free, get_ng0sh, kmesh_print, get_BZ_item, find_qmesh, make_mesh
 59  use m_double_grid,     only : double_grid_t, double_grid_init, double_grid_free
 60  use m_ebands,          only : ebands_init, ebands_print, ebands_copy, ebands_free, &
 61                                ebands_update_occ, get_valence_idx, apply_scissor, ebands_report_gap
 62  use m_kg,              only : getph
 63  use m_gsphere,         only : gsphere_t, gsph_free, gsph_init, print_gsphere, gsph_extend
 64  use m_vcoul,           only : vcoul_t, vcoul_init, vcoul_free
 65  use m_qparticles,      only : rdqps, rdgw  !, show_QP , rdgw
 66  use m_wfd,             only : wfd_init, wfd_free, wfd_print, wfd_t, wfd_test_ortho,&
 67                                wfd_read_wfk, wfd_wave_free, wfd_rotate, wfd_reset_ur_cprj, wfd_mkrho, test_charge
 68  use m_wfk,             only : wfk_read_eigenvalues
 69  use m_energies,        only : energies_type, energies_init
 70  use m_io_screening,    only : hscr_t, hscr_free, hscr_io, hscr_bcast, hscr_from_file, hscr_print
 71  use m_haydock,         only : exc_haydock_driver
 72  use m_exc_diago,       only : exc_diago_driver
 73  use m_exc_analyze,     only : exc_den
 74  use m_eprenorms,       only : eprenorms_t, eprenorms_free, eprenorms_from_epnc, eprenorms_bcast
 75  use m_pawang,          only : pawang_type
 76  use m_pawrad,          only : pawrad_type
 77  use m_pawtab,          only : pawtab_type, pawtab_print, pawtab_get_lsize
 78  use m_paw_an,          only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
 79  use m_paw_ij,          only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
 80  use m_pawfgrtab,       only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init
 81  use m_pawrhoij,        only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy,&
 82 &                              pawrhoij_free, pawrhoij_get_nspden, symrhoij
 83  use m_pawdij,          only : pawdij, symdij
 84  use m_pawfgr,          only : pawfgr_type, pawfgr_init, pawfgr_destroy
 85  use m_paw_hr,          only : pawhur_t, pawhur_free, pawhur_init
 86  use m_pawpwij,         only : pawpwff_t, pawpwff_init, pawpwff_free
 87  use m_paw_sphharm,     only : setsym_ylm
 88  use m_paw_denpot,      only : pawdenpot
 89  use m_paw_init,        only : pawinit,paw_gencond
 90  use m_paw_onsite,      only : pawnabla_init
 91  use m_paw_dmft,        only : paw_dmft_type
 92  use m_paw_mkrho,       only : denfgr
 93  use m_paw_nhat,        only : nhatgrid,pawmknhat
 94  use m_paw_tools,       only : chkpawovlp,pawprt
 95  use m_paw_correlations,only : pawpuxinit
 96  use m_exc_build,       only : exc_build_ham
 97  use m_setvtr,          only : setvtr
 98  use m_mkrho,           only : prtrhomxmn
 99  use m_pspini,          only : pspini
100  use m_drivexc,         only : mkdenpos
101 
102  implicit none
103 
104  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.

PARENTS

      driver

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.

CHILDREN

      bs_parameters_free,chkpawovlp,crystal_free,denfgr,destroy_mpi_enreg
      double_grid_free,ebands_free,ebands_update_occ,energies_init
      eprenorms_free,exc_build_ham,exc_den,exc_diago_driver
      exc_haydock_driver,fourdp,get_gftt,getph,gsph_free,hdr_free
      init_distribfft_seq,initmpi_seq,kmesh_free,metric,mkdenpos,mkrdim
      nhatgrid,paw_an_free,paw_an_init,paw_an_nullify,paw_gencond,paw_ij_free
      paw_ij_init,paw_ij_nullify,pawdenpot,pawdij,pawfgr_destroy,pawfgr_init
      pawfgrtab_free,pawfgrtab_init,pawhur_free,pawhur_init,pawinit,pawmknhat
      pawnabla_init,pawprt,pawpuxinit,pawpwff_free,pawpwff_init
      pawrhoij_alloc,pawrhoij_copy,pawrhoij_free,pawtab_get_lsize
      pawtab_print,print_ngfft,prtrhomxmn,pspini,rdqps,rotate_fft_mesh
      screen_free,screen_init,screen_nullify,setsym_ylm,setup_bse
      setup_bse_interp,setvtr,symdij,test_charge,timab,vcoul_free,wfd_free
      wfd_init,wfd_mkrho,wfd_print,wfd_read_wfk,wfd_reset_ur_cprj,wfd_rotate
      wfd_test_ortho,wfd_wave_free,wrtout,xmpi_bcast

SOURCE

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

PARENTS

      bethe_salpeter

CHILDREN

      apply_scissor,bsp_calctype2str,crystal_from_hdr,crystal_print
      ebands_copy,ebands_init,ebands_print,ebands_report_gap
      ebands_update_occ,eprenorms_bcast,eprenorms_from_epnc,find_qmesh
      get_bz_item,get_ng0sh,gsph_extend,gsph_init,hdr_init,hdr_update
      hdr_vs_dtset,hscr_bcast,hscr_free,hscr_from_file,hscr_print
      init_transitions,kmesh_init,kmesh_print,make_mesh,matrginv,metric
      mkrdim,pawrhoij_alloc,pawrhoij_copy,pawrhoij_free,print_bs_files
      print_bs_parameters,print_gsphere,print_ngfft,rdgw,setmesh,vcoul_init
      wfk_read_eigenvalues,wrtout,xmpi_bcast,xmpi_max,xmpi_split_work

SOURCE

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

PARENTS

      bethe_salpeter

CHILDREN

      apply_scissor,double_grid_init,ebands_copy,ebands_init,ebands_print
      ebands_report_gap,ebands_update_occ,find_qmesh,gsph_extend,gsph_init
      init_transitions,kmesh_init,kmesh_print,make_mesh,print_gsphere
      vcoul_init,wfk_read_eigenvalues,wrtout

SOURCE

2079 subroutine setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh,&
2080 & Kmesh_dense,Qmesh_dense,KS_BSt_dense,QP_bst_dense,Gsph_x,Gsph_c,Vcp_dense,Hdr_wfk_dense,ngfftf,grid,comm)
2081 
2082 
2083 !This section has been created automatically by the script Abilint (TD).
2084 !Do not modify the following lines by hand.
2085 #undef ABI_FUNC
2086 #define ABI_FUNC 'setup_bse_interp'
2087 !End of the abilint section
2088 
2089  implicit none
2090 
2091 !Arguments ------------------------------------
2092 !scalars
2093  integer,intent(in) :: comm
2094  type(dataset_type),intent(in) :: Dtset
2095  type(datafiles_type),intent(inout) :: Dtfil
2096  type(excparam),intent(inout) :: Bsp
2097  type(hdr_type),intent(out) :: Hdr_wfk_dense
2098  type(crystal_t),intent(in) :: Cryst
2099  type(kmesh_t),intent(in) :: Kmesh
2100  type(kmesh_t),intent(out) :: Kmesh_dense,Qmesh_dense
2101  type(ebands_t),intent(out) :: KS_BSt_dense,QP_Bst_dense
2102  type(double_grid_t),intent(out) :: grid
2103  type(vcoul_t),intent(out) :: Vcp_dense
2104  type(gsphere_t),intent(out) :: Gsph_x,Gsph_c
2105 !arrays
2106  integer,intent(in) :: ngfftf(18)
2107 
2108 !Local variables ------------------------------
2109 !scalars
2110  integer,parameter :: pertcase0=0,master=0
2111  integer :: bantot_dense,ib,ibtot,ik_ibz,isppol,jj
2112  integer :: nbnds_kss_dense
2113  integer :: spin,hexc_size
2114  integer :: my_rank
2115  integer :: it
2116  integer :: nprocs
2117  integer :: is1,is2,is3,is4
2118  real(dp) :: nelect_hdr_dense
2119  logical,parameter :: remove_inv=.FALSE.
2120  character(len=500) :: msg
2121  character(len=fnlen) :: wfk_fname_dense
2122  integer :: nqlwl
2123 !arrays
2124  integer,allocatable :: npwarr(:)
2125  real(dp),allocatable :: shiftk(:,:)
2126  real(dp),allocatable :: doccde(:),eigen(:),occfact(:)
2127  real(dp),pointer :: energies_p_dense(:,:,:)
2128  complex(dpc),allocatable :: gw_energy(:,:,:)
2129  integer,allocatable :: nbands_temp(:)
2130  integer :: kptrlatt_dense(3,3)
2131  real(dp),allocatable :: qlwl(:,:)
2132  real(dp) :: minmax_tene(2)
2133 
2134 !************************************************************************
2135 
2136  DBG_ENTER("COLL")
2137 
2138  kptrlatt_dense = zero
2139 
2140  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
2141 
2142  SELECT CASE(BSp%interp_mode)
2143  CASE (1,2,3,4)
2144    nbnds_kss_dense = -1
2145    wfk_fname_dense = Dtfil%fnameabi_wfkfine
2146    call wrtout(std_out,"BSE Interpolation: will read energies from: "//trim(wfk_fname_dense),"COLL")
2147 
2148    if (nctk_try_fort_or_ncfile(wfk_fname_dense, msg) /= 0) then
2149      MSG_ERROR(msg)
2150    end if
2151 
2152    Dtfil%fnameabi_wfkfine = wfk_fname_dense
2153 
2154    call wfk_read_eigenvalues(wfk_fname_dense,energies_p_dense,Hdr_wfk_dense,comm)
2155    nbnds_kss_dense = MAXVAL(Hdr_wfk_dense%nband)
2156  CASE DEFAULT
2157    MSG_ERROR("Not yet implemented")
2158  END SELECT
2159 
2160  nelect_hdr_dense = Hdr_wfk_dense%nelect
2161 
2162  if (ABS(Dtset%nelect-nelect_hdr_dense)>tol6) then
2163    write(msg,'(2(a,f8.2))')&
2164 &   "File contains ", nelect_hdr_dense," electrons but nelect initialized from input is ",Dtset%nelect
2165    MSG_ERROR(msg)
2166  end if
2167 
2168  ! Setup of the k-point list and symmetry tables in the  BZ
2169  SELECT CASE(BSp%interp_mode)
2170  CASE (1,2,3,4)
2171    if(Dtset%chksymbreak == 0) then
2172      ABI_MALLOC(shiftk,(3,Dtset%nshiftk))
2173      kptrlatt_dense(:,1) = BSp%interp_kmult(1)*Dtset%kptrlatt(:,1)
2174      kptrlatt_dense(:,2) = BSp%interp_kmult(2)*Dtset%kptrlatt(:,2)
2175      kptrlatt_dense(:,3) = BSp%interp_kmult(3)*Dtset%kptrlatt(:,3)
2176      do jj = 1,Dtset%nshiftk
2177        shiftk(:,jj) = Bsp%interp_kmult(:)*Dtset%shiftk(:,jj)
2178      end do
2179      call make_mesh(Kmesh_dense,Cryst,Dtset%kptopt,kptrlatt_dense,Dtset%nshiftk,shiftk,break_symmetry=.TRUE.)
2180      ABI_FREE(shiftk)
2181    else
2182      !Initialize Kmesh with no wrapping inside ]-0.5;0.5]
2183      call kmesh_init(Kmesh_dense,Cryst,Hdr_wfk_dense%nkpt,Hdr_wfk_dense%kptns,Dtset%kptopt)
2184    end if
2185  CASE DEFAULT
2186    MSG_ERROR("Not yet implemented")
2187  END SELECT
2188 
2189  ! Init Qmesh
2190  call find_qmesh(Qmesh_dense,Cryst,Kmesh_dense)
2191 
2192  call gsph_init(Gsph_c,Cryst,0,ecut=Dtset%ecuteps)
2193 
2194  call double_grid_init(Kmesh,Kmesh_dense,Dtset%kptrlatt,BSp%interp_kmult,grid)
2195 
2196  BSp%nkibz_interp = Kmesh_dense%nibz  !We might allow for a smaller number of points....
2197 
2198  call kmesh_print(Kmesh_dense,"Interpolated K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL")
2199  call kmesh_print(Kmesh_dense,"Interpolated K-mesh for the wavefunctions",ab_out, 0,           "COLL")
2200 
2201  if (nbnds_kss_dense < Dtset%nband(1)) then
2202    write(msg,'(2(a,i0),3a,i0)')&
2203 &    'Interpolated WFK file contains only ', nbnds_kss_dense,' levels instead of ',Dtset%nband(1),' required;',ch10,&
2204 &    'The calculation will be done with nbands= ',nbnds_kss_dense
2205    MSG_WARNING(msg)
2206    MSG_ERROR("Not supported yet !")
2207  end if
2208 
2209  ABI_MALLOC(nbands_temp,(Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol))
2210  do isppol=1,Hdr_wfk_dense%nsppol
2211    do ik_ibz=1,Hdr_wfk_dense%nkpt
2212      nbands_temp(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt) = Dtset%nband(1)
2213    end do
2214  end do
2215 
2216  call gsph_extend(Gsph_c,Cryst,Dtset%ecutwfn,Gsph_x)
2217  call print_gsphere(Gsph_x,unit=std_out,prtvol=Dtset%prtvol)
2218 
2219  nqlwl=1
2220  ABI_MALLOC(qlwl,(3,nqlwl))
2221  qlwl(:,nqlwl)= GW_Q0_DEFAULT
2222 
2223  ! Compute Coulomb term on the largest G-sphere.
2224  if (Gsph_x%ng > Gsph_c%ng ) then
2225    call vcoul_init(Vcp_dense,Gsph_x,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%icutcoul,&
2226 &    Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,ngfftf,comm)
2227  else
2228    call vcoul_init(Vcp_dense,Gsph_c,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%icutcoul,&
2229 &    Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,ngfftf,comm)
2230  end if
2231 
2232  ABI_FREE(qlwl)
2233 
2234  bantot_dense=SUM(Hdr_wfk_dense%nband(1:Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol))
2235  ABI_ALLOCATE(doccde,(bantot_dense))
2236  ABI_ALLOCATE(eigen,(bantot_dense))
2237  ABI_ALLOCATE(occfact,(bantot_dense))
2238  doccde=zero; eigen=zero; occfact=zero
2239 
2240  jj=0; ibtot=0
2241  do isppol=1,Hdr_wfk_dense%nsppol
2242    do ik_ibz=1,Hdr_wfk_dense%nkpt
2243      do ib=1,Hdr_wfk_dense%nband(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt)
2244        ibtot=ibtot+1
2245        if (ib<=BSP%nbnds) then
2246          jj=jj+1
2247          occfact(jj)=Hdr_wfk_dense%occ(ibtot)
2248          eigen  (jj)=energies_p_dense(ib,ik_ibz,isppol)
2249        end if
2250      end do
2251    end do
2252  end do
2253 
2254  ABI_FREE(energies_p_dense)
2255 
2256  ABI_MALLOC(npwarr,(kmesh_dense%nibz))
2257  npwarr=BSP%npwwfn
2258 
2259  call ebands_init(bantot_dense,KS_BSt_dense,Dtset%nelect,doccde,eigen,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,nbands_temp,&
2260 &  Kmesh_dense%nibz,npwarr,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,Hdr_wfk_dense%tphysel,Hdr_wfk_dense%tsmear,&
2261 &  Hdr_wfk_dense%occopt,occfact,Kmesh_dense%wt,&
2262 &  hdr_wfk_dense%charge, hdr_wfk_dense%kptopt, hdr_wfk_dense%kptrlatt_orig, hdr_wfk_dense%nshiftk_orig, &
2263 &  hdr_wfk_dense%shiftk_orig, hdr_wfk_dense%kptrlatt, hdr_wfk_dense%nshiftk, hdr_wfk_dense%shiftk)
2264 
2265  ABI_DEALLOCATE(doccde)
2266  ABI_DEALLOCATE(eigen)
2267  ABI_DEALLOCATE(npwarr)
2268 
2269  ABI_FREE(nbands_temp)
2270 
2271  ABI_FREE(occfact)
2272 
2273  !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when
2274  ! the occupation scheme for semiconductors is used.
2275  call ebands_update_occ(KS_BSt_dense,Dtset%spinmagntarget,prtvol=Dtset%prtvol)
2276 
2277  call ebands_print(KS_BSt_dense,"Interpolated band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol)
2278 
2279  call ebands_report_gap(KS_BSt_dense,header="Interpolated KS band structure",unit=std_out,mode_paral="COLL")
2280 
2281  BSp%nkbz_interp = Kmesh_dense%nbz
2282 
2283  call ebands_copy(KS_BSt_dense,QP_bst_dense)
2284 
2285  SELECT CASE (Bsp%calc_type)
2286  CASE (BSE_HTYPE_RPA_KS)
2287    if (ABS(BSp%mbpt_sciss)>tol6) then
2288      write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies."
2289      call wrtout(std_out,msg,"COLL")
2290      call apply_scissor(QP_BSt_dense,BSp%mbpt_sciss)
2291    else
2292      write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]."
2293      call wrtout(std_out,msg,"COLL")
2294    end if
2295    !
2296  CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw
2297    MSG_ERROR("Not yet implemented with interpolation !")
2298  CASE (BSE_HTYPE_RPA_QP)
2299    MSG_ERROR("Not implemented error!")
2300  CASE DEFAULT
2301    write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type
2302    MSG_ERROR(msg)
2303  END SELECT
2304 
2305  call ebands_report_gap(QP_BSt_dense,header=" Interpolated QP band structure",unit=std_out,mode_paral="COLL")
2306 
2307  ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index.
2308  ! FIXME: linewidths not coded.
2309  ABI_ALLOCATE(gw_energy,(BSp%nbnds,Kmesh_dense%nibz,Dtset%nsppol))
2310  gw_energy = QP_BSt_dense%eig
2311 
2312  ABI_ALLOCATE(Bsp%nreh_interp,(Hdr_wfk_dense%nsppol))
2313  Bsp%nreh_interp=zero
2314 
2315  call init_transitions(BSp%Trans_interp,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz_interp,Bsp%nbnds,&
2316 &  Bsp%nkibz_interp,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,gw_energy,QP_BSt_dense%occ,Kmesh_dense%tab,minmax_tene,&
2317 &  Bsp%nreh_interp)
2318 
2319  ABI_DEALLOCATE(gw_energy)
2320 
2321  do spin=1,Dtset%nsppol
2322    write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh_interp(spin)
2323    call wrtout(std_out,msg,"COLL")
2324  end do
2325 
2326  if (ANY(Bsp%nreh_interp/=Bsp%nreh_interp(1))) then
2327    write(msg,'(a,(i0))')" BSE code does not support different number of transitions for the two spin channels",Bsp%nreh
2328    MSG_ERROR(msg)
2329  end if
2330  !
2331  ! Create transition table vcks2t
2332  is1=BSp%lomo_min;is2=BSp%homo_max;is3=BSp%lumo_min;is4=BSp%humo_max
2333  ABI_ALLOCATE(Bsp%vcks2t_interp,(is1:is2,is3:is4,BSp%nkbz_interp,Dtset%nsppol))
2334  Bsp%vcks2t_interp = 0
2335 
2336  do spin=1,Dtset%nsppol
2337    do it=1,BSp%nreh_interp(spin)
2338      BSp%vcks2t_interp(BSp%Trans_interp(it,spin)%v,BSp%Trans_interp(it,spin)%c,&
2339 & BSp%Trans_interp(it,spin)%k,spin) = it
2340    end do
2341  end do
2342 
2343  hexc_size = SUM(Bsp%nreh_interp); if (Bsp%use_coupling>0) hexc_size=2*hexc_size
2344  if (Bsp%nstates_interp<=0) then
2345    Bsp%nstates_interp=hexc_size
2346  else
2347    if (Bsp%nstates_interp>hexc_size) then
2348       Bsp%nstates_interp=hexc_size
2349       write(msg,'(2(a,i0),2a)')&
2350 &      "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates_interp,ch10,&
2351 &      "the number of excitonic states nstates has been modified"
2352      MSG_WARNING(msg)
2353    end if
2354  end if
2355 
2356  DBG_EXIT("COLL")
2357 
2358 end subroutine setup_bse_interp