TABLE OF CONTENTS


ABINIT/setup_bse [ Functions ]

[ Top ] [ 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.

COPYRIGHT

 Copyright (C) 1992-2009 EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida)
 Copyright (C) 2009-2018 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

 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

  62 #if defined HAVE_CONFIG_H
  63 #include "config.h"
  64 #endif
  65 
  66 #include "abi_common.h"
  67 
  68 subroutine setup_bse(codvsn,acell,rprim,ngfftf,ngfft_osc,Dtset,Dtfil,BS_files,Psps,Pawtab,BSp,&
  69 & Cryst,Kmesh,Qmesh,KS_BSt,QP_bst,Hdr_wfk,Gsph_x,Gsph_c,Vcp,Hdr_bse,w_fname,Epren,comm,Wvl)
  70 
  71  use defs_basis
  72  use defs_datatypes
  73  use defs_abitypes
  74  use defs_wvltypes
  75  use m_bs_defs
  76  use m_profiling_abi
  77  use m_errors
  78  use m_xmpi
  79  use m_nctk
  80  use m_hdr
  81 
  82  use m_gwdefs,        only : GW_Q0_DEFAULT
  83  use m_fstrings,      only : toupper, sjoin
  84  use m_io_tools,      only : file_exists, open_file
  85  use m_geometry,      only : normv, mkrdim, metric
  86  use m_abilasi,       only : matrginv
  87  use m_crystal,       only : crystal_print, idx_spatial_inversion, crystal_t
  88  use m_crystal_io,    only : crystal_from_hdr
  89  use m_bz_mesh,       only : kmesh_t, kmesh_init, get_ng0sh, kmesh_print, get_BZ_item, find_qmesh, make_mesh
  90  use m_ebands,        only : ebands_init, ebands_print, ebands_copy, ebands_free, &
  91 &                            ebands_update_occ, get_valence_idx, apply_scissor, ebands_report_gap
  92  use m_eprenorms,     only : eprenorms_t, eprenorms_from_epnc, eprenorms_bcast
  93  use m_vcoul,         only : vcoul_t, vcoul_init
  94  use m_fftcore,       only : print_ngfft
  95  use m_fft_mesh,      only : setmesh
  96  use m_gsphere,       only : gsphere_t, gsph_init, print_gsphere, merge_and_sort_kg, gsph_extend
  97  use m_io_screening,  only : hscr_t, hscr_free, hscr_io, hscr_bcast, hscr_from_file, hscr_print
  98  use m_pawtab,        only : pawtab_type
  99  use m_pawrhoij,      only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free
 100  use m_qparticles,    only : rdgw
 101  use m_screen,        only : MDL_BECHSTEDT
 102  use m_wfk,           only : wfk_read_eigenvalues
 103 
 104 !This section has been created automatically by the script Abilint (TD).
 105 !Do not modify the following lines by hand.
 106 #undef ABI_FUNC
 107 #define ABI_FUNC 'setup_bse'
 108  use interfaces_14_hidewrite
 109  use interfaces_56_io_mpi
 110 !End of the abilint section
 111 
 112  implicit none
 113 
 114 !Arguments ------------------------------------
 115 !scalars
 116  integer,intent(in) :: comm
 117  character(len=6),intent(in) :: codvsn
 118  character(len=fnlen),intent(out) :: w_fname
 119  type(dataset_type),intent(inout) :: Dtset
 120  type(datafiles_type),intent(in) :: Dtfil
 121  type(pseudopotential_type),intent(in) :: Psps
 122  type(excparam),intent(inout) :: Bsp
 123  type(hdr_type),intent(out) :: Hdr_wfk,Hdr_bse
 124  type(crystal_t),intent(out) :: Cryst
 125  type(kmesh_t),intent(out) :: Kmesh,Qmesh
 126  type(gsphere_t),intent(out) :: Gsph_x,Gsph_c
 127  type(ebands_t),intent(out) :: KS_BSt,QP_Bst
 128  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw)
 129  type(vcoul_t),intent(out) :: Vcp
 130  type(excfiles),intent(out) :: BS_files
 131  type(wvl_internal_type), intent(in) :: Wvl
 132  type(eprenorms_t),intent(out) :: Epren
 133 !arrays
 134  integer,intent(in) :: ngfftf(18)
 135  integer,intent(out) :: ngfft_osc(18)
 136  real(dp),intent(in) :: acell(3),rprim(3,3)
 137 
 138 !Local variables ------------------------------
 139 !scalars
 140  integer,parameter :: pertcase0=0,master=0
 141  integer(i8b) :: work_size,tot_nreh,neh_per_proc,il
 142  integer :: bantot,enforce_sym,ib,ibtot,ik_ibz,isppol,jj,method,iat,ount !ii,
 143  integer :: mband,io,nfftot_osc,spin,hexc_size,nqlwl,iq
 144  integer :: timrev,iq_bz,isym,iq_ibz,itim
 145  integer :: my_rank,nprocs,fform,npwe_file,ierr,my_k1, my_k2,my_nbks
 146  integer :: first_dig,second_dig,it
 147  real(dp) :: ucvol,qnorm
 148  real(dp):: eff,mempercpu_mb,wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem
 149  logical,parameter :: remove_inv=.FALSE.
 150  logical :: ltest,occ_from_dtset
 151  character(len=500) :: msg
 152  character(len=fnlen) :: gw_fname,test_file,wfk_fname
 153  character(len=fnlen) :: ep_nc_fname
 154  type(hscr_t) :: Hscr
 155 !arrays
 156  integer :: ng0sh_opt(3),val_idx(Dtset%nsppol)
 157  integer,allocatable :: npwarr(:),val_indeces(:,:),nlmn_atm(:)
 158  real(dp) :: qpt_bz(3),minmax_tene(2)
 159  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3)
 160  real(dp) :: qred2cart(3,3),qcart2red(3,3)
 161  real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:)
 162  real(dp),allocatable :: igwene(:,:,:)
 163  real(dp),pointer :: energies_p(:,:,:)
 164  complex(dpc),allocatable :: gw_energy(:,:,:)
 165  type(Pawrhoij_type),allocatable :: Pawrhoij(:)
 166 !Interp@BSE
 167  !integer :: mode
 168  !integer :: kmult(3)
 169  !integer :: unt
 170  !integer :: rl_nb
 171  !logical :: interp_params_exists, prepare, sum_overlaps
 172  !namelist /interp_params/ mode,kmult,prepare,rl_nb,sum_overlaps
 173  !character(len=fnlen) :: tmp_fname
 174 
 175 !************************************************************************
 176 
 177  DBG_ENTER("COLL")
 178 
 179  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
 180 
 181  ! === Check for calculations that are not implemented ===
 182  ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1))
 183  ABI_CHECK(ltest,'Dtset%nband must be constant')
 184  ABI_CHECK(Dtset%nspinor==1,"nspinor==2 not coded")
 185 
 186  ! === Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume ===
 187  call mkrdim(acell,rprim,rprimd)
 188  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 189 
 190  ! Read energies and header from the WFK file.
 191  wfk_fname = dtfil%fnamewffk
 192  if (.not. file_exists(wfk_fname)) then
 193    wfk_fname = nctk_ncify(wfk_fname)
 194    MSG_COMMENT(sjoin("File not found. Will try netcdf file: ", wfk_fname))
 195  end if
 196 
 197  call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm)
 198  mband = MAXVAL(Hdr_wfk%nband)
 199 
 200  call hdr_vs_dtset(Hdr_wfk,Dtset)
 201 
 202  ! === Create crystal_t data type ===
 203  !remove_inv= .FALSE. !(nsym_kss/=Hdr_wfk%nsym)
 204  timrev=  2 ! This information is not reported in the header
 205             ! 1 => do not use time-reversal symmetry
 206             ! 2 => take advantage of time-reversal symmetry
 207 
 208  call crystal_from_hdr(Cryst,Hdr_wfk,timrev,remove_inv)
 209  call crystal_print(Cryst)
 210  !
 211  ! Setup of the k-point list and symmetry tables in the  BZ -----------------------------------
 212  if (Dtset%chksymbreak==0) then
 213    MSG_WARNING("Calling make_mesh")
 214    call make_mesh(Kmesh,Cryst,Dtset%kptopt,Dtset%kptrlatt,Dtset%nshiftk,Dtset%shiftk,break_symmetry=.TRUE.)
 215    ! TODO
 216    !Check if kibz from KSS file corresponds to the one returned by make_mesh.
 217  else
 218    call kmesh_init(Kmesh,Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt)
 219  end if
 220  BSp%nkibz = Kmesh%nibz  !We might allow for a smaller number of points....
 221 
 222  call kmesh_print(Kmesh,"K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL")
 223  call kmesh_print(Kmesh,"K-mesh for the wavefunctions",ab_out, 0,           "COLL")
 224 
 225  nqlwl = 0
 226  w_fname = ABI_NOFILE
 227  if (Dtset%getscr/=0.or.Dtset%irdscr/=0) then
 228    w_fname=Dtfil%fnameabi_scr
 229  else if (Dtset%getsuscep/=0.or.Dtset%irdsuscep/=0) then
 230    w_fname=Dtfil%fnameabi_sus
 231    MSG_ERROR("(get|ird)suscep not implemented")
 232  end if
 233 
 234  if (w_fname /= ABI_NOFILE) then ! Read dimensions from the external file
 235 
 236    if (.not. file_exists(w_fname)) then
 237      w_fname = nctk_ncify(w_fname)
 238      MSG_COMMENT(sjoin("File not found. Will try netcdf file: ", w_fname))
 239    end if
 240 
 241    if (my_rank==master) then
 242      ! Master reads npw and nqlwl from SCR file.
 243      call wrtout(std_out,sjoin('Testing file: ', w_fname),"COLL")
 244 
 245      call hscr_from_file(hscr,w_fname,fform,xmpi_comm_self)
 246      ! Echo the header.
 247      if (Dtset%prtvol>0) call hscr_print(Hscr)
 248 
 249      npwe_file = Hscr%npwe ! Have to change %npweps if it was larger than dim on disk.
 250      nqlwl     = Hscr%nqlwl
 251 
 252      if (Dtset%npweps>npwe_file) then
 253        write(msg,'(2(a,i0),2a,i0)')&
 254 &       "The number of G-vectors stored on file (",npwe_file,") is smaller than Dtset%npweps = ",Dtset%npweps,ch10,&
 255 &       "Calculation will proceed with the maximum available set, npwe_file = ",npwe_file
 256        MSG_WARNING(msg)
 257        Dtset%npweps = npwe_file
 258      else  if (Dtset%npweps<npwe_file) then
 259        write(msg,'(2(a,i0),2a,i0)')&
 260 &       "The number of G-vectors stored on file (",npwe_file,") is larger than Dtset%npweps = ",Dtset%npweps,ch10,&
 261 &       "Calculation will proceed with Dtset%npweps = ",Dtset%npweps
 262        MSG_COMMENT(msg)
 263      end if
 264    end if
 265 
 266    call hscr_bcast(Hscr,master,my_rank,comm)
 267    call xmpi_bcast(Dtset%npweps,master,comm,ierr)
 268    call xmpi_bcast(nqlwl,master,comm,ierr)
 269 
 270    if (nqlwl>0) then
 271      ABI_MALLOC(qlwl,(3,nqlwl))
 272      qlwl = Hscr%qlwl
 273    end if
 274    !
 275    ! Init Qmesh from the SCR file.
 276    call kmesh_init(Qmesh,Cryst,Hscr%nqibz,Hscr%qibz,Dtset%kptopt)
 277 
 278    ! The G-sphere for W and Sigma_c is initialized from the gvectors found in the SCR file.
 279    call gsph_init(Gsph_c,Cryst,Dtset%npweps,gvec=Hscr%gvec)
 280 
 281    call hscr_free(Hscr)
 282  else
 283    ! Init Qmesh from the K-mesh reported in the WFK file.
 284    call find_qmesh(Qmesh,Cryst,Kmesh)
 285 
 286    ! The G-sphere for W and Sigma_c is initialized from ecuteps.
 287    call gsph_init(Gsph_c,Cryst,0,ecut=Dtset%ecuteps)
 288    Dtset%npweps = Gsph_c%ng
 289  end if
 290 
 291  BSp%npweps = Dtset%npweps
 292  BSp%ecuteps = Dtset%ecuteps
 293 
 294  if (nqlwl==0) then
 295    nqlwl=1
 296    ABI_MALLOC(qlwl,(3,nqlwl))
 297    qlwl(:,nqlwl)= GW_Q0_DEFAULT
 298    write(msg,'(3a,i2,a,3f9.6)')&
 299 &    "The Header of the screening file does not contain the list of q-point for the optical limit ",ch10,&
 300 &    "Using nqlwl= ",nqlwl," and qlwl = ",qlwl(:,1)
 301    MSG_COMMENT(msg)
 302  end if
 303  write(std_out,*)"nqlwl and qlwl for Coulomb singularity and e^-1",nqlwl,qlwl
 304 
 305  ! === Setup of q-mesh in the whole BZ ===
 306  ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed,
 307  !   epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed)
 308  !
 309  call kmesh_print(Qmesh,"Q-mesh for the screening function",std_out,Dtset%prtvol,"COLL")
 310  call kmesh_print(Qmesh,"Q-mesh for the screening function",ab_out ,0           ,"COLL")
 311 
 312  do iq_bz=1,Qmesh%nbz
 313    call get_BZ_item(Qmesh,iq_bz,qpt_bz,iq_ibz,isym,itim)
 314    sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz))
 315    if (ANY(ABS(Qmesh%bz(:,iq_bz)-sq )>1.0d-4)) then
 316      write(std_out,*) sq,Qmesh%bz(:,iq_bz)
 317      write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')&
 318 &      'qpoint ',Qmesh%bz(:,iq_bz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,&
 319 &      'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,&
 320 &      'however a non zero umklapp G_o vector is required and this is not yet allowed'
 321      MSG_ERROR(msg)
 322    end if
 323  end do
 324 
 325  BSp%algorithm = Dtset%bs_algorithm
 326  BSp%nstates   = Dtset%bs_nstates
 327  Bsp%nsppol    = Dtset%nsppol
 328  Bsp%hayd_term = Dtset%bs_hayd_term
 329  !
 330  ! Define the algorithm for solving the BSE.
 331  if (BSp%algorithm == BSE_ALGO_HAYDOCK) then
 332    BSp%niter       = Dtset%bs_haydock_niter
 333    BSp%haydock_tol = Dtset%bs_haydock_tol
 334 
 335  else if (BSp%algorithm == BSE_ALGO_CG) then
 336    ! FIXME For the time being use an hardcoded value.
 337    ! TODO change name in Dtset%
 338    BSp%niter       = Dtset%nstep !100
 339    BSp%cg_tolwfr   = Dtset%tolwfr
 340    BSp%nline       = Dtset%nline
 341    BSp%nbdbuf      = Dtset%nbdbuf
 342    BSp%nstates     = Dtset%bs_nstates
 343    MSG_WARNING("Check CG setup")
 344  else
 345    !BSp%niter       = 0
 346    !BSp%tol_iter    = HUGE(one)
 347  end if
 348  !
 349  ! Shall we include Local field effects?
 350  SELECT CASE (Dtset%bs_exchange_term)
 351  CASE (0,1)
 352    BSp%exchange_term = Dtset%bs_exchange_term
 353  CASE DEFAULT
 354    write(msg,'(a,i0)')" Wrong bs_exchange_term: ",Dtset%bs_exchange_term
 355    MSG_ERROR(msg)
 356  END SELECT
 357  !
 358  ! Treatment of the off-diagonal coupling block.
 359  SELECT CASE (Dtset%bs_coupling)
 360  CASE (0)
 361    BSp%use_coupling = 0
 362    msg = 'RESONANT ONLY CALCULATION'
 363  CASE (1)
 364    BSp%use_coupling = 1
 365    msg = ' RESONANT+COUPLING CALCULATION '
 366  CASE DEFAULT
 367    write(msg,'(a,i0)')" Wrong bs_coupling: ",Dtset%bs_coupling
 368    MSG_ERROR(msg)
 369  END SELECT
 370  call wrtout(std_out,msg,"COLL")
 371 
 372  BSp%use_diagonal_Wgg = .FALSE.
 373  Bsp%use_coulomb_term = .TRUE.
 374  BSp%eps_inf=zero
 375  Bsp%mdlf_type=0
 376 
 377  first_dig =MOD(Dtset%bs_coulomb_term,10)
 378  second_dig=Dtset%bs_coulomb_term/10
 379 
 380  Bsp%wtype = second_dig
 381  SELECT CASE (second_dig)
 382  CASE (BSE_WTYPE_NONE)
 383    call wrtout(std_out,"Coulomb term won't be calculated","COLL")
 384    Bsp%use_coulomb_term = .FALSE.
 385 
 386  CASE (BSE_WTYPE_FROM_SCR)
 387    call wrtout(std_out,"W is read from an external SCR file","COLL")
 388    Bsp%use_coulomb_term = .TRUE.
 389 
 390  CASE (BSE_WTYPE_FROM_MDL)
 391    call wrtout(std_out,"W is approximated with the model dielectric function","COLL")
 392    Bsp%use_coulomb_term = .TRUE.
 393    BSp%mdlf_type = MDL_BECHSTEDT
 394    BSp%eps_inf = Dtset%mdf_epsinf
 395    ABI_CHECK(Bsp%eps_inf > zero, "mdf_epsinf <= 0")
 396 
 397  CASE DEFAULT
 398    write(msg,'(a,i0)')" Wrong second digit in bs_coulomb_term: ",Dtset%bs_coulomb_term
 399    MSG_ERROR(msg)
 400  END SELECT
 401  !
 402  ! Diagonal approximation or full matrix?
 403  BSp%use_diagonal_Wgg = .TRUE.
 404  if (Bsp%wtype /= BSE_WTYPE_NONE) then
 405    SELECT CASE (first_dig)
 406    CASE (0)
 407      call wrtout(std_out,"Using diagonal approximation W_GG","COLL")
 408      BSp%use_diagonal_Wgg = .TRUE.
 409    CASE (1)
 410      call wrtout(std_out,"Using full W_GG' matrix ","COLL")
 411      BSp%use_diagonal_Wgg = .FALSE.
 412    CASE DEFAULT
 413      write(msg,'(a,i0)')" Wrong first digit in bs_coulomb_term: ",Dtset%bs_coulomb_term
 414      MSG_ERROR(msg)
 415    END SELECT
 416  end if
 417 
 418  !TODO move the initialization of the parameters for the interpolation in setup_bse_interp
 419 
 420  BSp%use_interp = .FALSE.
 421  BSp%interp_mode = BSE_INTERP_YG
 422  BSp%interp_kmult(1:3) = 0
 423  BSp%prep_interp = .FALSE.
 424  BSp%sum_overlaps = .TRUE. ! Sum over the overlaps
 425 
 426  ! Printing ncham
 427  BSp%prt_ncham = .FALSE.
 428 
 429  ! Deactivate Interpolation Technique by default
 430 ! if (.FALSE.) then
 431 
 432  ! Reading parameters from the input file
 433  BSp%use_interp = (dtset%bs_interp_mode /= 0)
 434  BSp%prep_interp = (dtset%bs_interp_prep == 1)
 435 
 436  SELECT CASE (dtset%bs_interp_mode)
 437  CASE (0)
 438    ! No interpolation, do not print anything !
 439  CASE (1)
 440    call wrtout(std_out,"Using interpolation technique with energies and wavefunctions from dense WFK","COLL")
 441  CASE (2)
 442    call wrtout(std_out,"Interpolation technique with energies and wfn on dense WFK + treatment ABC of divergence","COLL")
 443  CASE (3)
 444    call wrtout(std_out,"Interpolation technique + divergence ABC along diagonal","COLL")
 445  CASE (4)
 446    call wrtout(std_out,"Using interpolation technique mode 1 with full computation of hamiltonian","COLL")
 447  CASE DEFAULT
 448    write(msg,'(a,i0)')" Wrong interpolation mode for bs_interp_mode: ",dtset%bs_interp_mode
 449    MSG_ERROR(msg)
 450  END SELECT
 451 
 452  ! Read from dtset
 453  if(BSp%use_interp) then
 454    BSp%interp_method = dtset%bs_interp_method
 455    BSp%rl_nb = dtset%bs_interp_rl_nb
 456    BSp%interp_m3_width = dtset%bs_interp_m3_width
 457    BSp%interp_kmult(1:3) = dtset%bs_interp_kmult(1:3)
 458    BSp%interp_mode = dtset%bs_interp_mode
 459  end if
 460 
 461  ! Dimensions and parameters of the calculation.
 462  ! TODO one should add npwx as well
 463  !BSp%npweps=Dtset%npweps
 464  !BSp%npwwfn=Dtset%npwwfn
 465 
 466  ABI_MALLOC(Bsp%lomo_spin, (Bsp%nsppol))
 467  ABI_MALLOC(Bsp%homo_spin, (Bsp%nsppol))
 468  ABI_MALLOC(Bsp%lumo_spin, (Bsp%nsppol))
 469  ABI_MALLOC(Bsp%humo_spin, (Bsp%nsppol))
 470  ABI_MALLOC(Bsp%nbndv_spin, (Bsp%nsppol))
 471  ABI_MALLOC(Bsp%nbndc_spin, (Bsp%nsppol))
 472 
 473  ! FIXME use bs_loband(nsppol)
 474  Bsp%lomo_spin = Dtset%bs_loband
 475  write(std_out,*)"bs_loband",Dtset%bs_loband
 476  !if (Bsp%nsppol == 2) Bsp%lomo_spin(2) = Dtset%bs_loband
 477 
 478  ! Check lomo correct only for unpolarized semiconductors
 479  !if (Dtset%nsppol == 1 .and. Bsp%lomo > Dtset%nelect/2) then
 480  !  write(msg,'(a,i0,a,f8.3)') " Bsp%lomo = ",Bsp%lomo," cannot be greater than nelect/2 = ",Dtset%nelect/2
 481  !  MSG_ERROR(msg)
 482  !end if
 483  !
 484  ! ==============================================
 485  ! ==== Setup of the q for the optical limit ====
 486  ! ==============================================
 487  Bsp%inclvkb = Dtset%inclvkb
 488 
 489  qred2cart = two_pi*Cryst%gprimd
 490  qcart2red = qred2cart
 491  call matrginv(qcart2red,3,3)
 492 
 493  if (Dtset%gw_nqlwl==0) then
 494    BSp%nq = 6
 495    ABI_MALLOC(BSp%q,(3,BSp%nq))
 496    BSp%q(:,1) = (/one,zero,zero/)  ! (100)
 497    BSp%q(:,2) = (/zero,one,zero/)  ! (010)
 498    BSp%q(:,3) = (/zero,zero,one/)  ! (001)
 499    BSp%q(:,4) = MATMUL(qcart2red,(/one,zero,zero/)) ! (x)
 500    BSp%q(:,5) = MATMUL(qcart2red,(/zero,one,zero/)) ! (y)
 501    BSp%q(:,6) = MATMUL(qcart2red,(/zero,zero,one/)) ! (z)
 502  else
 503    BSp%nq = Dtset%gw_nqlwl
 504    ABI_MALLOC(BSp%q,(3,BSp%nq))
 505    BSp%q = Dtset%gw_qlwl
 506  end if
 507 
 508  do iq=1,BSp%nq ! normalization
 509    qnorm = normv(BSp%q(:,iq),Cryst%gmet,"G")
 510    BSp%q(:,iq) = BSp%q(:,iq)/qnorm
 511  end do
 512  !
 513  ! ======================================================
 514  ! === Define the flags defining the calculation type ===
 515  ! ======================================================
 516  Bsp%calc_type = Dtset%bs_calctype
 517 
 518  BSp%mbpt_sciss = zero ! Shall we use the scissors operator to open the gap?
 519  if (ABS(Dtset%mbpt_sciss)>tol6) BSp%mbpt_sciss = Dtset%mbpt_sciss
 520 
 521 !now test input parameters from input and WFK file and assume some defaults
 522 !
 523 ! TODO Add the possibility of using a randomly shifted k-mesh with nsym>1.
 524 ! so that densities and potentials are correctly symmetrized but
 525 ! the list of the k-point in the IBZ is not expanded.
 526 
 527  if (mband < Dtset%nband(1)) then
 528    write(msg,'(2(a,i0),3a,i0)')&
 529 &    'WFK file contains only ', mband,' levels instead of ',Dtset%nband(1),' required;',ch10,&
 530 &    'The calculation will be done with nbands= ',mband
 531    MSG_WARNING(msg)
 532    Dtset%nband(:) = mband
 533  end if
 534 
 535  BSp%nbnds = Dtset%nband(1) ! TODO Note the change in the meaning of input variables
 536 
 537  if (BSp%nbnds<=Dtset%nelect/2) then
 538    write(msg,'(2a,a,i0,a,f8.2)')&
 539 &    'BSp%nbnds cannot be smaller than homo ',ch10,&
 540 &    'while BSp%nbnds = ',BSp%nbnds,' and Dtset%nelect = ',Dtset%nelect
 541    MSG_ERROR(msg)
 542  end if
 543 
 544 !TODO add new dim for exchange part and consider the possibility of having npwsigx > npwwfn (see setup_sigma).
 545 
 546  ! === Build enlarged G-sphere for the exchange part ===
 547  call gsph_extend(Gsph_c,Cryst,Dtset%ecutwfn,Gsph_x)
 548  call print_gsphere(Gsph_x,unit=std_out,prtvol=Dtset%prtvol)
 549 
 550  ! NPWVEC as the biggest between npweps and npwwfn. MG RECHECK this part.
 551  !BSp%npwwfn = Dtset%npwwfn
 552  Bsp%npwwfn = Gsph_x%ng  ! FIXME temporary hack
 553  BSp%npwvec=MAX(BSp%npwwfn,BSp%npweps)
 554  Bsp%ecutwfn = Dtset%ecutwfn
 555 
 556  ! Compute Coulomb term on the largest G-sphere.
 557  if (Gsph_x%ng > Gsph_c%ng ) then
 558    call vcoul_init(Vcp,Gsph_x,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,&
 559 &    nqlwl,qlwl,ngfftf,comm)
 560  else
 561    call vcoul_init(Vcp,Gsph_c,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,&
 562 &    nqlwl,qlwl,ngfftf,comm)
 563  end if
 564 
 565  ABI_FREE(qlwl)
 566 
 567  bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol))
 568  ABI_MALLOC(doccde,(bantot))
 569  ABI_MALLOC(eigen,(bantot))
 570  ABI_MALLOC(occfact,(bantot))
 571  doccde=zero; eigen=zero; occfact=zero
 572 
 573  ! Get occupation from input if occopt == 2
 574  occ_from_dtset = (Dtset%occopt == 2)
 575 
 576  jj=0; ibtot=0
 577  do isppol=1,Dtset%nsppol
 578    do ik_ibz=1,Dtset%nkpt
 579      do ib=1,Hdr_wfk%nband(ik_ibz+(isppol-1)*Dtset%nkpt)
 580        ibtot=ibtot+1
 581        if (ib<=BSP%nbnds) then
 582          jj=jj+1
 583          eigen  (jj)=energies_p(ib,ik_ibz,isppol)
 584          if (occ_from_dtset) then
 585            occfact(jj)=Dtset%occ_orig(ibtot)
 586          else
 587            occfact(jj)=Hdr_wfk%occ(ibtot)
 588          end if
 589        end if
 590      end do
 591    end do
 592  end do
 593 
 594  ABI_FREE(energies_p)
 595  !
 596  ! * Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of
 597  !   symmetry operations in the old GW code (symmorphy and inversion)
 598  ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6))
 599  ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt')
 600 
 601  ABI_MALLOC(npwarr,(Dtset%nkpt))
 602  npwarr=BSP%npwwfn
 603 
 604  call ebands_init(bantot,KS_BSt,Dtset%nelect,doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,&
 605 &  Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,&
 606 &  dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
 607 &  dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
 608 
 609  ABI_FREE(doccde)
 610  ABI_FREE(eigen)
 611  ABI_FREE(npwarr)
 612 
 613  !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when
 614  ! the occupation scheme for semiconductors is used.
 615  call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=Dtset%prtvol)
 616 
 617  call ebands_print(KS_BSt,"Band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol)
 618 
 619  call ebands_report_gap(KS_BSt,header=" KS band structure",unit=std_out,mode_paral="COLL")
 620 
 621  ABI_MALLOC(val_indeces,(KS_BSt%nkpt,KS_BSt%nsppol))
 622  val_indeces = get_valence_idx(KS_BSt)
 623 
 624  do spin=1,KS_BSt%nsppol
 625    val_idx(spin) = val_indeces(1,spin)
 626    write(msg,'(a,i2,a,i0)')" For spin : ",spin," val_idx ",val_idx(spin)
 627    call wrtout(std_out,msg,"COLL")
 628    if ( ANY(val_indeces(1,spin) /= val_indeces(:,spin)) ) then
 629      MSG_ERROR("BSE code does not support metals")
 630    end if
 631  end do
 632 
 633  ABI_FREE(val_indeces)
 634  !
 635  ! === Create the BSE header ===
 636  call hdr_init(KS_BSt,codvsn,Dtset,Hdr_bse,Pawtab,pertcase0,Psps,wvl)
 637 
 638  ! === Get Pawrhoij from the header of the WFK file ===
 639  ABI_DT_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw))
 640  if (Dtset%usepaw==1) then
 641    call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab)
 642    call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij)
 643  end if
 644 
 645  call hdr_update(hdr_bse,bantot,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
 646 
 647  ABI_FREE(occfact)
 648 
 649  if (Dtset%usepaw==1) call pawrhoij_free(Pawrhoij)
 650  ABI_DT_FREE(Pawrhoij)
 651 
 652  ! === Find optimal value for G-sphere enlargment due to oscillator matrix elements ===
 653 
 654  ! We will split k-points over processors
 655  call xmpi_split_work(Kmesh%nbz,comm,my_k1,my_k2,msg,ierr)
 656  if (ierr/=0) then
 657    MSG_WARNING(msg)
 658  end if
 659 
 660  ! If there is no work to do, just skip the computation
 661  if (my_k2-my_k1+1 <= 0) then
 662    ng0sh_opt(:)=(/zero,zero,zero/)
 663  else
 664    ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true since bz is buggy
 665    ! * -one is used because we loop over all the possibile differences, unlike screening
 666    call get_ng0sh(my_k2-my_k1+1,Kmesh%bz(:,my_k1:my_k2),Kmesh%nbz,Kmesh%bz,&
 667 &    Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt)
 668  end if
 669 
 670  call xmpi_max(ng0sh_opt,BSp%mg0,comm,ierr)
 671 
 672  write(msg,'(a,3(i0,1x))') ' optimal value for ng0sh = ',BSp%mg0
 673  call wrtout(std_out,msg,"COLL")
 674 
 675  ! === Setup of the FFT mesh for the oscilator strengths ===
 676  ! * ngfft_osc(7:18)==Dtset%ngfft(7:18) which is initialized before entering screening.
 677  ! * Here we redefine ngfft_osc(1:6) according to the following options :
 678  !
 679  ! method==0 --> FFT grid read from fft.in (debugging purpose)
 680  ! method==1 --> Normal FFT mesh
 681  ! method==2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
 682  ! method==3 --> Doubled FFT grid, same as the the FFT for the density,
 683  !
 684  ! enforce_sym==1 ==> Enforce a FFT mesh compatible with all the symmetry operation and FFT library
 685  ! enforce_sym==0 ==> Find the smallest FFT grid compatbile with the library, do not care about symmetries
 686  !
 687  ngfft_osc(1:18)=Dtset%ngfft(1:18); method=2
 688  if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0
 689  if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1
 690  if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2
 691  if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3
 692  enforce_sym=MOD(Dtset%fftgw,10)
 693 
 694  call setmesh(gmet,Gsph_x%gvec,ngfft_osc,BSp%npwvec,BSp%npweps,BSp%npwwfn,nfftot_osc,method,BSp%mg0,Cryst,enforce_sym)
 695  nfftot_osc=PRODUCT(ngfft_osc(1:3))
 696 
 697  call print_ngfft(ngfft_osc,"FFT mesh for oscillator matrix elements",std_out,"COLL",prtvol=Dtset%prtvol)
 698  !
 699  ! BSp%homo gives the
 700  !BSp%homo  = val_idx(1)
 701  ! highest occupied band for each spin
 702  BSp%homo_spin = val_idx
 703 
 704  ! TODO generalize the code to account for this unlikely case.
 705  !if (Dtset%nsppol==2) then
 706  !  ABI_CHECK(BSp%homo == val_idx(2),"Different valence indeces for spin up and down")
 707  !end if
 708 
 709  !BSp%lumo = BSp%homo + 1
 710  !BSp%humo = BSp%nbnds
 711  !BSp%nbndv = BSp%homo  - BSp%lomo + 1
 712  !BSp%nbndc = BSp%nbnds - BSp%homo
 713 
 714  BSp%lumo_spin = BSp%homo_spin + 1
 715  BSp%humo_spin = BSp%nbnds
 716  BSp%nbndv_spin = BSp%homo_spin  - BSp%lomo_spin + 1
 717  BSp%nbndc_spin = BSp%nbnds - BSp%homo_spin
 718  BSp%maxnbndv = MAXVAL(BSp%nbndv_spin(:))
 719  BSp%maxnbndc = MAXVAL(BSp%nbndc_spin(:))
 720 
 721  BSp%nkbz = Kmesh%nbz
 722 
 723  call ebands_copy(KS_BSt,QP_bst)
 724  ABI_MALLOC(igwene,(QP_bst%mband,QP_bst%nkpt,QP_bst%nsppol))
 725  igwene=zero
 726 
 727  call bsp_calctype2str(Bsp,msg)
 728  call wrtout(std_out,"Calculation type: "//TRIM(msg))
 729 
 730  SELECT CASE (Bsp%calc_type)
 731  CASE (BSE_HTYPE_RPA_KS)
 732    if (ABS(BSp%mbpt_sciss)>tol6) then
 733      write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies."
 734      call wrtout(std_out,msg,"COLL")
 735      call apply_scissor(QP_BSt,BSp%mbpt_sciss)
 736    else
 737      write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]."
 738      call wrtout(std_out,msg,"COLL")
 739    end if
 740 
 741  CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw
 742    gw_fname=TRIM(Dtfil%filnam_ds(4))//'_GW'
 743    gw_fname="__in.gw__"
 744    if (.not.file_exists(gw_fname)) then
 745      msg = " File "//TRIM(gw_fname)//" not found. Aborting now"
 746      MSG_ERROR(msg)
 747    end if
 748 
 749    call rdgw(QP_Bst,gw_fname,igwene,extrapolate=.FALSE.) ! here gwenergy is real
 750 
 751    do isppol=1,Dtset%nsppol
 752      write(std_out,*) ' k       GW energies [eV]'
 753      do ik_ibz=1,Kmesh%nibz
 754        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)
 755      end do
 756      write(std_out,*) ' k       Im GW energies [eV]'
 757      do ik_ibz=1,Kmesh%nibz
 758        write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(igwene(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds)
 759      end do
 760    end do
 761    !
 762    ! If required apply the scissors operator on top of the QP bands structure (!)
 763    if (ABS(BSp%mbpt_sciss)>tol6) then
 764      write(msg,'(a,f8.2,a)')' Applying a scissors operator ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the QP energies!"
 765      MSG_COMMENT(msg)
 766      call apply_scissor(QP_BSt,BSp%mbpt_sciss)
 767    end if
 768 
 769  CASE (BSE_HTYPE_RPA_QP)
 770    MSG_ERROR("Not implemented error!")
 771 
 772  CASE DEFAULT
 773    write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type
 774    MSG_ERROR(msg)
 775  END SELECT
 776 
 777  call ebands_report_gap(QP_BSt,header=" QP band structure",unit=std_out,mode_paral="COLL")
 778 
 779  ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index.
 780  ! FIXME: linewidths not coded.
 781  ABI_MALLOC(gw_energy,(BSp%nbnds,Kmesh%nibz,Dtset%nsppol))
 782  gw_energy = QP_BSt%eig
 783 
 784  BSp%have_complex_ene = ANY(igwene > tol16)
 785 
 786  ! Compute the number of resonant transitions, nreh, for the two spin channels and initialize BSp%Trans.
 787  ABI_MALLOC(Bsp%nreh,(Bsp%nsppol))
 788 
 789  ! Possible cutoff on the transitions.
 790  BSp%ircut = Dtset%bs_eh_cutoff(1)
 791  BSp%uvcut = Dtset%bs_eh_cutoff(2)
 792 
 793  call init_transitions(BSp%Trans,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz,Bsp%nbnds,Bsp%nkibz,&
 794 &  BSp%nsppol,Dtset%nspinor,gw_energy,QP_BSt%occ,Kmesh%tab,minmax_tene,Bsp%nreh)
 795 
 796  ! Setup of the frequency mesh for the absorption spectrum.
 797  ! If not specified, use the min-max resonant transition energy and make it 10% smaller|larger.
 798 
 799  !if (ABS(Dtset%bs_freq_mesh(1)) < tol6) then
 800  !   Dtset%bs_freq_mesh(1) = MAX(minmax_tene(1) - minmax_tene(1) * 0.1, zero)
 801  !end if
 802 
 803  if (ABS(Dtset%bs_freq_mesh(2)) < tol6) then
 804     Dtset%bs_freq_mesh(2) = minmax_tene(2) + minmax_tene(2) * 0.1
 805  end if
 806 
 807  Bsp%omegai = Dtset%bs_freq_mesh(1)
 808  Bsp%omegae = Dtset%bs_freq_mesh(2)
 809  Bsp%domega = Dtset%bs_freq_mesh(3)
 810  BSp%broad  = Dtset%zcut
 811 
 812  ! The frequency mesh (including the complex imaginary shift)
 813  BSp%nomega = (BSp%omegae - BSp%omegai)/BSp%domega + 1
 814  ABI_MALLOC(BSp%omega,(BSp%nomega))
 815  do io=1,BSp%nomega
 816    BSp%omega(io) = (BSp%omegai + (io-1)*BSp%domega)  + j_dpc*BSp%broad
 817  end do
 818 
 819  ABI_FREE(gw_energy)
 820  ABI_FREE(igwene)
 821 
 822  do spin=1,Bsp%nsppol
 823    write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh(spin)
 824    call wrtout(std_out,msg,"COLL")
 825  end do
 826 
 827  if (ANY(Bsp%nreh/=Bsp%nreh(1))) then
 828    write(msg,'(a,2(i0,1x))')" BSE code with different number of transitions for the two spin channels: ",Bsp%nreh
 829    MSG_WARNING(msg)
 830  end if
 831  !
 832  ! Create transition table vcks2t
 833  Bsp%lomo_min = MINVAL(BSp%lomo_spin)
 834  Bsp%homo_max = MAXVAL(BSp%homo_spin)
 835  Bsp%lumo_min = MINVAL(BSp%lumo_spin)
 836  Bsp%humo_max = MAXVAL(BSp%humo_spin)
 837 
 838  ABI_MALLOC(Bsp%vcks2t,(BSp%lomo_min:BSp%homo_max,BSp%lumo_min:BSp%humo_max,BSp%nkbz,Dtset%nsppol))
 839  Bsp%vcks2t = 0
 840 
 841  do spin=1,BSp%nsppol
 842    do it=1,BSp%nreh(spin)
 843      BSp%vcks2t(BSp%Trans(it,spin)%v,BSp%Trans(it,spin)%c,BSp%Trans(it,spin)%k,spin) = it
 844    end do
 845  end do
 846 
 847  hexc_size = SUM(Bsp%nreh); if (Bsp%use_coupling>0) hexc_size=2*hexc_size
 848  if (Bsp%nstates<=0) then
 849    Bsp%nstates=hexc_size
 850  else
 851    if (Bsp%nstates>hexc_size) then
 852       Bsp%nstates=hexc_size
 853       write(msg,'(2(a,i0),2a)')&
 854 &      "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates,ch10,&
 855 &      "the number of excitonic states nstates has been modified"
 856      MSG_WARNING(msg)
 857    end if
 858  end if
 859 
 860  msg=' Fundamental parameters for the solution of the Bethe-Salpeter equation:'
 861  call print_bs_parameters(BSp,unit=std_out,header=msg,mode_paral="COLL",prtvol=Dtset%prtvol)
 862  call print_bs_parameters(BSp,unit=ab_out, header=msg,mode_paral="COLL")
 863 
 864  if (ANY (Cryst%symrec(:,:,1) /= RESHAPE ( (/1,0,0,0,1,0,0,0,1/),(/3,3/) )) .or. &
 865 &    ANY( ABS(Cryst%tnons(:,1)) > tol6) ) then
 866    write(msg,'(3a,9i2,2a,3f6.3,2a)')&
 867 &    "The first symmetry operation should be the Identity with zero tnons while ",ch10,&
 868 &    "symrec(:,:,1) = ",Cryst%symrec(:,:,1),ch10,&
 869 &    "tnons(:,1)    = ",Cryst%tnons(:,1),ch10,&
 870 &    "This is not allowed, sym_rhotwgq0 should be changed."
 871    MSG_ERROR(msg)
 872  end if
 873  !
 874  ! Prefix for generic output files.
 875  BS_files%out_basename = TRIM(Dtfil%filnam_ds(4))
 876  !
 877  ! Search for files to restart from.
 878  if (Dtset%gethaydock/=0 .or. Dtset%irdhaydock/=0) then
 879    BS_files%in_haydock_basename = TRIM(Dtfil%fnameabi_haydock)
 880  end if
 881 
 882  test_file = Dtfil%fnameabi_bsham_reso
 883  if (file_exists(test_file)) then
 884    BS_files%in_hreso = test_file
 885  else
 886    BS_files%out_hreso = TRIM(Dtfil%filnam_ds(4))//'_BSR'
 887  end if
 888 
 889  test_file = Dtfil%fnameabi_bsham_coup
 890  if (file_exists(test_file) ) then
 891    BS_files%in_hcoup = test_file
 892  else
 893    BS_files%out_hcoup = TRIM(Dtfil%filnam_ds(4))//'_BSC'
 894  end if
 895  !
 896  ! in_eig is the name of the input file with eigenvalues and eigenvectors
 897  ! constructed from getbseig or irdbseig. out_eig is the name of the output file
 898  ! produced by this dataset. in_eig_exists checks for the presence of the input file.
 899  !
 900  if (file_exists(Dtfil%fnameabi_bseig)) then
 901    BS_files%in_eig = Dtfil%fnameabi_bseig
 902  else
 903    BS_files%out_eig = TRIM(BS_files%out_basename)//"_BSEIG"
 904  end if
 905 
 906  call print_bs_files(BS_files,unit=std_out)
 907 
 908 
 909  !
 910  ! ==========================================================
 911  ! ==== Temperature dependence of the spectrum ==============
 912  ! ==========================================================
 913  BSp%do_ep_renorm = .FALSE.
 914  BSp%do_lifetime = .FALSE. ! Not yet implemented
 915 
 916  ep_nc_fname = 'test_EP.nc'
 917  if(file_exists(ep_nc_fname)) then
 918    BSp%do_ep_renorm = .TRUE.
 919 
 920    if(my_rank == master) then
 921      call eprenorms_from_epnc(Epren,ep_nc_fname)
 922    end if
 923    call eprenorms_bcast(Epren,master,comm)
 924  end if
 925 
 926 
 927  !
 928  ! ==========================================================
 929  ! ==== Final check on the parameters of the calculation ====
 930  ! ==========================================================
 931  if ( Bsp%use_coupling>0 .and. ALL(Bsp%algorithm /= [BSE_ALGO_DDIAGO, BSE_ALGO_HAYDOCK]) ) then
 932    msg = "Resonant+Coupling is only available with the direct diagonalization or the haydock method."
 933    MSG_ERROR(msg)
 934  end if
 935 
 936  ! autoparal section
 937  if (dtset%max_ncpus /=0 .and. dtset%autoparal /=0 ) then
 938    ount = ab_out
 939    ! TODO:
 940    ! nsppol and calculation with coupling!
 941 
 942    ! Temporary table needed to estimate memory
 943    ABI_MALLOC(nlmn_atm,(Cryst%natom))
 944    if (Dtset%usepaw==1) then
 945      do iat=1,Cryst%natom
 946        nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size
 947      end do
 948    end if
 949 
 950    tot_nreh = SUM(BSp%nreh)
 951    work_size = tot_nreh * (tot_nreh + 1) / 2
 952 
 953    write(ount,'(a)')"--- !Autoparal"
 954    write(ount,"(a)")'#Autoparal section for Bethe-Salpeter runs.'
 955 
 956    write(ount,"(a)")   "info:"
 957    write(ount,"(a,i0)")"    autoparal: ",dtset%autoparal
 958    write(ount,"(a,i0)")"    max_ncpus: ",dtset%max_ncpus
 959    write(ount,"(a,i0)")"    nkibz: ",Bsp%nkibz
 960    write(ount,"(a,i0)")"    nkbz: ",Bsp%nkbz
 961    write(ount,"(a,i0)")"    nsppol: ",dtset%nsppol
 962    write(ount,"(a,i0)")"    nspinor: ",dtset%nspinor
 963    write(ount,"(a,i0)")"    lomo_min: ",Bsp%lomo_min
 964    write(ount,"(a,i0)")"    humo_max: ",Bsp%humo_max
 965    write(ount,"(a,i0)")"    tot_nreh: ",tot_nreh
 966    !write(ount,"(a,i0)")"    nbnds: ",Ep%nbnds
 967 
 968    ! Wavefunctions are not distributed. We read all the bands
 969    ! from 1 up to Bsp%nbnds because we have to recompute rhor
 970    ! but then we deallocate all the states that are not used for the construction of the e-h
 971    ! before allocating the EXC hamiltonian. Hence we can safely use  (humo - lomo + 1) instead of Bsp%nbnds.
 972    !my_nbks = (Bsp%humo - Bsp%lomo +1) * Bsp%nkibz * Dtset%nsppol
 973 
 974    ! This one overestimates the memory but it seems to be safer.
 975    my_nbks = Bsp%nbnds * Dtset%nkpt * Dtset%nsppol
 976 
 977    ! Memory needed for Fourier components ug.
 978    ug_mem = two*gwpc*Dtset%nspinor*Bsp%npwwfn*my_nbks*b2Mb
 979 
 980    ! Memory needed for real space ur.
 981    ur_mem = zero
 982    if (MODULO(Dtset%gwmem,10)==1) then
 983      ur_mem = two*gwpc*Dtset%nspinor*nfftot_osc*my_nbks*b2Mb
 984    end if
 985 
 986    ! Memory needed for PAW projections Cprj
 987    cprj_mem = zero
 988    if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb
 989 
 990    wfsmem_mb = ug_mem + ur_mem + cprj_mem
 991 
 992    ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI:  wavefunctions + W
 993    nonscal_mem = (wfsmem_mb + two*gwpc*BSp%npweps**2*b2Mb) * 1.1_dp
 994 
 995    ! List of configurations.
 996    write(ount,"(a)")"configurations:"
 997    do il=1,dtset%max_ncpus
 998      if (il > work_size) cycle
 999      neh_per_proc = work_size / il
1000      neh_per_proc = neh_per_proc + MOD(work_size, il)
1001      eff = (one * work_size) / (il * neh_per_proc)
1002 
1003      ! EXC matrix is distributed.
1004      mempercpu_mb = nonscal_mem + two * dpc * neh_per_proc * b2Mb
1005 
1006      write(ount,"(a,i0)")"    - tot_ncpus: ",il
1007      write(ount,"(a,i0)")"      mpi_ncpus: ",il
1008      !write(ount,"(a,i0)")"      omp_ncpus: ",omp_ncpus
1009      write(ount,"(a,f12.9)")"      efficiency: ",eff
1010      write(ount,"(a,f12.2)")"      mem_per_cpu: ",mempercpu_mb
1011    end do
1012 
1013    write(ount,'(a)')"..."
1014 
1015    ABI_FREE(nlmn_atm)
1016    MSG_ERROR_NODUMP("aborting now")
1017  end if
1018 
1019  DBG_EXIT("COLL")
1020 
1021 end subroutine setup_bse