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

SOURCE

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)
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
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
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(:,:,:,:)
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.
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  ! ===========================================
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,&
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
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,&
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
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,&
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.
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
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
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
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.
555    do spin=1,Bsp%nsppol
557    end do
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,&
569 &   Gsph_x_dense%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut)
570
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)
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,&
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,&
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,&
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,&
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!
876
878  do spin=1,Bsp%nsppol
881  end do
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
905    W_info%use_mdf = BSp%mdlf_type
906    !W_info%use_ppm
907    !W_info%vtx_test
908    !W_info%wint_method
909    !
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)
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")
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/)
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

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
1208  wfk_fname = dtfil%fnamewffk
1209  if (.not. file_exists(wfk_fname)) then
1210    wfk_fname = nctk_ncify(wfk_fname)
1212  end if
1213
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)
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)
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
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
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)
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:'
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

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