TABLE OF CONTENTS


ABINIT/m_sigma_driver [ Modules ]

[ Top ] [ Modules ]

NAME

  m_sigma_driver

FUNCTION

 Calculate the matrix elements of the self-energy operator.

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (MG, GMR, VO, LR, RWG, MT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

 21 #if defined HAVE_CONFIG_H
 22 #include "config.h"
 23 #endif
 24 
 25 #include "abi_common.h"
 26 
 27 module m_sigma_driver
 28 
 29  use defs_basis
 30  use m_gwdefs
 31  use defs_datatypes
 32  use defs_abitypes
 33  use defs_wvltypes
 34  use m_xmpi
 35  use m_xomp
 36  use m_errors
 37  use m_abicore
 38  use m_ab7_mixing
 39  use m_nctk
 40  use m_kxc
 41 #ifdef HAVE_NETCDF
 42  use netcdf
 43 #endif
 44  use m_hdr
 45  use libxc_functionals
 46  use m_wfd
 47 
 48  use m_time,          only : timab
 49  use m_numeric_tools, only : imax_loc
 50  use m_fstrings,      only : strcat, sjoin, itoa, basename, ktoa, ltoa
 51  use m_hide_blas,     only : xdotc
 52  use m_io_tools,      only : open_file, file_exists, iomode_from_fname
 53  use m_mpinfo,        only : destroy_mpi_enreg, initmpi_seq
 54  use m_geometry,      only : normv, mkrdim, metric
 55  use m_fftcore,       only : print_ngfft
 56  use m_fft_mesh,      only : get_gftt, setmesh
 57  use m_fft,           only : fourdp
 58  use m_ioarr,         only : fftdatar_write, read_rhor
 59  use m_crystal,       only : crystal_free, crystal_t, crystal_print, idx_spatial_inversion
 60  use m_crystal_io,    only : crystal_ncwrite, crystal_from_hdr
 61  use m_ebands,        only : ebands_update_occ, ebands_copy, ebands_report_gap, get_valence_idx, get_bandenergy, &
 62 &                            ebands_free, ebands_init, ebands_ncwrite, ebands_interpolate_kpath, get_eneocc_vect, &
 63                              enclose_degbands, get_gaps, gaps_free, gaps_t, gaps_print
 64  use m_energies,      only : energies_type, energies_init
 65  use m_bz_mesh,       only : kmesh_t, kmesh_free, littlegroup_t, littlegroup_init, littlegroup_free, &
 66                              kmesh_init, has_BZ_item, isamek, get_ng0sh, kmesh_print, &
 67                              get_bz_item, has_IBZ_item, find_qmesh
 68  use m_gsphere,       only : gsphere_t, gsph_init, gsph_free, merge_and_sort_kg, gsph_extend, setshells
 69  use m_kg,            only : getph
 70  use m_xcdata,        only : get_xclevel
 71  use m_vcoul,         only : vcoul_t, vcoul_init, vcoul_free
 72  use m_qparticles,    only : wrqps, rdqps, rdgw, show_QP, updt_m_lda_to_qp
 73  use m_screening,     only : mkdump_er, em1results_free, epsilonm1_results, init_er_from_file
 74  use m_ppmodel,       only : ppm_init, ppm_free, setup_ppmodel, getem1_from_PPm, ppmodel_t
 75  use m_sigma,         only : sigma_init, sigma_free, sigma_ncwrite, sigma_t, sigma_get_exene, &
 76 &                            write_sigma_header, write_sigma_results
 77  use m_dyson_solver,  only : solve_dyson
 78  use m_esymm,         only : esymm_t, esymm_free, esymm_failed
 79  use m_melemts,       only : melflags_reset, melements_print, melements_free, melflags_t, melements_t, melements_zero
 80  use m_pawang,        only : pawang_type
 81  use m_pawrad,        only : pawrad_type
 82  use m_pawtab,        only : pawtab_type, pawtab_print, pawtab_get_lsize
 83  use m_paw_an,        only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
 84  use m_paw_ij,        only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_print
 85  use m_pawfgrtab,     only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print
 86  use m_pawrhoij,      only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_get_nspden, symrhoij, &
 87                              pawrhoij_unpack
 88  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_free, paw_overlap
 89  use m_pawdij,        only : pawdij, symdij_all
 90  use m_pawfgr,        only : pawfgr_type, pawfgr_init, pawfgr_destroy
 91  use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
 92  use m_pawpwij,       only : pawpwff_t, pawpwff_init, pawpwff_free
 93  use m_paw_slater,    only : paw_mkdijexc_core, paw_dijhf
 94  use m_paw_dmft,      only : paw_dmft_type
 95  use m_paw_sphharm,   only : setsym_ylm
 96  use m_paw_mkrho,     only : denfgr
 97  use m_paw_nhat,      only : nhatgrid,pawmknhat
 98  use m_paw_tools,     only : chkpawovlp,pawprt
 99  use m_paw_denpot,    only : pawdenpot
100  use m_paw_init,      only : pawinit,paw_gencond
101  use m_classify_bands,only : classify_bands
102  use m_wfk,           only : wfk_read_eigenvalues
103  use m_io_kss,        only : make_gvec_kss
104  use m_vhxc_me,       only : calc_vhxc_me
105  use m_cohsex,        only : cohsex_me
106  use m_sigx,          only : calc_sigx_me
107  use m_sigc,          only : calc_sigc_me
108  use m_setvtr,        only : setvtr
109  use m_mkrho,         only : prtrhomxmn
110  use m_pspini,        only : pspini
111  use m_calc_ucrpa,    only : calc_ucrpa
112  use m_prep_calc_ucrpa,only : prep_calc_ucrpa
113 
114  use m_paw_correlations,only : pawpuxinit
115 
116  implicit none
117 
118  private

ABINIT/paw_qpscgw [ Functions ]

[ Top ] [ Functions ]

NAME

 paw_qpscgw

FUNCTION

  This routine is called during QP self-consistent GW calculations. It calculates the new QP on-site quantities
  using the QP amplitudes read from the QPS file.

INPUTS

  Wfd<wfd_t>=Datatype gathering data on QP amplitudes.
  nscf=Number of QPSCGW iterations done so far (read from the QPS file).
  nfftf=Number of points in the fine FFT grid.
  ngfft(18)=information on the fine FFT grid used for densities and potentials.
  Dtset<dataset_type>=All input variables for this dataset.
  Cryst<crystal_t>=Info on unit cell and symmetries.
  Kmesh<kmesh_t>=Structure describing the k-point sampling.
  Psps<Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the KSS file
  Pawang<pawang_type>=PAW angular mesh and related data.
  Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data
  Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data
  Pawfgrtab(natom)<Pawfgrtab_type>= For PAW, various arrays giving data related to fine grid for a given atom.
  prev_Pawrhoij(Cryst%natom))<Pawrhoij_type>=Previous QP rhoij used for mixing if nscf>0 and rhoqpmix/=one.
  MPI_enreg=Information about MPI parallelization.
  QP_BSt<ebands_t>=QP band structure.
  QP_energies<Energies_type>=Simple datastructure to gather all part of total energy.
  nhatgrdim= 0 if pawgrnhat array is not used ; 1 otherwise

OUTPUT

  QP_pawrhoij(Cryst%natom))<Pawrhoij_type>=on-site densities calculated from the QP amplitudes.
  qp_nhat(nfftf,Dtset%nspden)=Compensation charge density calculated from the QP amplitudes.
  qp_nhatgr(nfftf,Dtset%nspden,3*nhatgrdim)=Derivatives of the QP nhat on fine rectangular grid (and derivatives).
  qp_compch_sph=QP compensation charge integral inside spheres computed over spherical meshes.
  qp_compch_fft=QP compensation charge inside spheres computed over fine fft grid.
  QP_paw_ij(Cryst%natom)<Paw_ij_type>=Non-local D_ij strengths of the QP Hamiltonian.
  QP_paw_an(Cryst%natom)<Paw_an_type>=Various arrays related to the Hamiltonian given
   on ANgular mesh or ANgular moments.

SIDE EFFECTS

PARENTS

      sigma

CHILDREN

      paw_an_init,paw_an_nullify,paw_ij_init,paw_ij_nullify,pawdenpot
      pawmknhat,pawrhoij_alloc,pawrhoij_unpack,symrhoij,wfd_pawrhoij,wrtout

SOURCE

3839 subroutine paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,QP_BSt,&
3840 &  Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij,&
3841 &  QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim,qp_nhatgr,qp_compch_sph,qp_compch_fft)
3842 
3843 
3844 !This section has been created automatically by the script Abilint (TD).
3845 !Do not modify the following lines by hand.
3846 #undef ABI_FUNC
3847 #define ABI_FUNC 'paw_qpscgw'
3848 !End of the abilint section
3849 
3850  implicit none
3851 
3852 !Arguments ------------------------------------
3853 !scalars
3854  integer,intent(in) :: nfftf,nscf,nhatgrdim
3855  real(dp),intent(out) :: qp_compch_fft,qp_compch_sph
3856  type(kmesh_t),intent(in) :: Kmesh
3857  type(crystal_t),intent(in) :: Cryst
3858  type(Dataset_type),intent(in) :: Dtset
3859  type(Pseudopotential_type),intent(in) :: Psps
3860  type(Pawang_type),intent(in) :: Pawang
3861  type(ebands_t),intent(in) :: QP_BSt
3862  type(Energies_type),intent(inout) :: QP_energies
3863  type(wfd_t),intent(inout) :: Wfd
3864 !arrays
3865  integer,intent(in) :: ngfftf(18)
3866  real(dp),intent(out) :: qp_nhat(nfftf,Dtset%nspden)
3867  real(dp),intent(out) :: qp_nhatgr(nfftf,Dtset%nspden,3*nhatgrdim)
3868  type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom)
3869  type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw)
3870  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
3871  type(Pawrhoij_type),intent(inout) :: prev_Pawrhoij(Cryst%natom)
3872  type(Pawrhoij_type),intent(out) :: QP_pawrhoij(Cryst%natom)
3873  type(Paw_ij_type),intent(out) :: QP_paw_ij(Cryst%natom)
3874  type(Paw_an_type),intent(inout) :: QP_paw_an(Cryst%natom)
3875 
3876 !Local variables ------------------------------
3877 !scalars
3878  integer :: choice,cplex,has_dijU,has_dijso,iat,ider,idir,ipert
3879  integer :: izero,nkxc1,nspden_rhoij,nzlmopt
3880  integer :: option,optrhoij,usexcnhat
3881  character(len=500) :: msg
3882 !arrays
3883  real(dp) :: k0(3)
3884 
3885 !************************************************************************
3886 
3887  ABI_UNUSED(Kmesh%nibz)
3888  !
3889  !  * 0 if Vloc in atomic data is Vbare    (Blochl s formulation)
3890  !  * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation)
3891  usexcnhat=MAXVAL(Pawtab(:)%usexcnhat)
3892  !
3893  ! Calculate new rhoij_qp from updated Cprj_ibz, note use_rhoij_=1.
3894  nspden_rhoij=pawrhoij_get_nspden(Dtset%nspden,Dtset%nspinor,Dtset%pawspnorb)
3895 
3896  call pawrhoij_alloc(QP_pawrhoij,Dtset%pawcpxocc,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,&
3897 &                 pawtab=Pawtab,use_rhoij_=1,use_rhoijres=1)
3898 
3899  ! FIXME kptop should be passed via Kmesh, in GW time reversal is always assumed.
3900  call wfd_pawrhoij(Wfd,Cryst,QP_Bst,Dtset%kptopt,QP_pawrhoij,Dtset%pawprtvol)
3901  !
3902  ! * Symmetrize QP $\rho_{ij}$.
3903  choice=1; optrhoij=1; ipert=0
3904  call symrhoij(QP_pawrhoij,QP_pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert,&
3905 &              Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,&
3906 &              Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
3907 
3908  ! ======================
3909  ! ==== Make QP nhat ====
3910  ! ======================
3911  cplex=1; ider=2*nhatgrdim; idir=0; ipert=0; izero=0; k0(:)=zero
3912 
3913  call pawmknhat(qp_compch_fft,cplex,ider,idir,ipert,izero,Cryst%gprimd,&
3914 &  Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
3915 &  Pawfgrtab,qp_nhatgr,qp_nhat,QP_Pawrhoij,QP_Pawrhoij,Pawtab,k0,Cryst%rprimd,&
3916 &  Cryst%ucvol,Dtset%usewvl,Cryst%xred)
3917 
3918  ! Allocate quantities related to the PAW spheres for the QP Hamiltonian.
3919  ! TODO call paw_ij_init in scfcv and respfn, fix small issues
3920  has_dijso=Dtset%pawspnorb; has_dijU=Dtset%usepawu
3921 
3922  call paw_ij_nullify(QP_paw_ij)
3923  call paw_ij_init(QP_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,&
3924 &  Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
3925 &  has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,&
3926 &  has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1)
3927 
3928  call paw_an_nullify(QP_paw_an); nkxc1=0 ! No kernel
3929  call paw_an_init(QP_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
3930 &     cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1)
3931 
3932  ! =====================================================
3933  ! ==== Optional mixing of the PAW onsite densities ====
3934  ! =====================================================
3935  if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then
3936    write(msg,'(a,f6.3)')' sigma: mixing on-site QP rho_ij densities using rhoqpmix= ',Dtset%rhoqpmix
3937    call wrtout(std_out,msg,'COLL')
3938    ! qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor)
3939 
3940    call pawrhoij_unpack(QP_Pawrhoij)   ! Unpack new QP %rhoijp
3941    call pawrhoij_unpack(prev_Pawrhoij) ! Unpack previous QP %rhoijp
3942 
3943    do iat=1,Cryst%natom
3944      QP_pawrhoij(iat)%rhoij_ = prev_Pawrhoij(iat)%rhoij_ &
3945 &      + Dtset%rhoqpmix * (QP_pawrhoij(iat)%rhoij_ - prev_pawrhoij(iat)%rhoij_)
3946 
3947      prev_pawrhoij(iat)%use_rhoij_=0
3948      ABI_DEALLOCATE(prev_pawrhoij(iat)%rhoij_)
3949    end do
3950    !
3951    ! * Re-Symmetrize mixed QP $\rho_{ij}$.
3952    choice=1; optrhoij=1; ipert=0
3953    call symrhoij(QP_pawrhoij,QP_pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert,&
3954 &                Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,&
3955 &                Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
3956  end if
3957 
3958  do iat=1,Cryst%natom
3959    QP_pawrhoij(iat)%use_rhoij_=0
3960    ABI_DEALLOCATE(QP_pawrhoij(iat)%rhoij_)
3961  end do
3962  !
3963  ! =================================================================================
3964  ! ==== Evaluate on-site energies, potentials, densities using (mixed) QP rhoij ====
3965  ! =================================================================================
3966  ! * Initialize also "lmselect" (index of non-zero LM-moments of densities).
3967 
3968  nzlmopt=-1; option=0; qp_compch_sph=greatest_real
3969 
3970  call pawdenpot(qp_compch_sph,QP_energies%e_paw,QP_energies%e_pawdc,&
3971 &  ipert,Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
3972 &  Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,QP_paw_an,QP_paw_an,&
3973 &  QP_paw_ij,Pawang,Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,&
3974 &  Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp)
3975 
3976 end subroutine paw_qpscgw

ABINIT/setup_sigma [ Functions ]

[ Top ] [ Functions ]

NAME

 setup_sigma

FUNCTION

  Initialize the data type containing parameters for a sigma calculation.

INPUTS

 acell(3)=length scales of primitive translations (bohr)
 wfk_fname=Name of the WFK file.
 Dtset<type(dataset_type)>=all input variables for this dataset
 Dtfil<type(datafiles_type)>=variables related to files
 rprim(3,3)=dimensionless real space primitive translations
 ngfft(18)=information on the (fine) FFT grid used for the density.
 Psps <Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the WFK file

OUTPUT

 Sigp<sigparams_t>=Parameters governing the self-energy calculation.
 Kmesh <kmesh_t>=Structure describing the k-point sampling.
 Qmesh <kmesh_t>=Structure describing the q-point sampling.
 Cryst<crystal_t>=Info on unit cell and symmetries.
 Gsph_Max<gsphere_t>=Info on the G-sphere
 Gsph_c<gsphere_t>=Info on the G-sphere for W and Sigma_c
 Gsph_x<gsphere_t>=Info on the G-sphere for and Sigma_x
 Hdr_wfk<hdr_type>=The header of the WFK file
 Hdr_out<hdr_type>=The header to be used for the results of sigma calculations.
 Vcp<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique.
 Er<Epsilonm1_results>=Datatype storing data used to construct the screening (partially Initialized in OUTPUT)
 KS_BSt<ebands_t>=The KS energies and occupation factors.
 gwc_ngfft(18), gwx_ngfft(18)= FFT meshes for the oscillator strengths used for the correlated and the
   exchange part of the self-energy, respectively.
 comm=MPI communicator.

PARENTS

      sigma

CHILDREN

SOURCE

2507 subroutine setup_sigma(codvsn,wfk_fname,acell,rprim,ngfftf,Dtset,Dtfil,Psps,Pawtab,&
2508 & gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,KS_BSt,Gsph_Max,Gsph_x,Gsph_c,Vcp,Er,Sigp,comm)
2509 
2510  !use m_gwdefs,        only : GW_Q0_DEFAULT, SIG_GW_AC, sigparams_t, sigma_is_herm, sigma_needs_w
2511 
2512 !This section has been created automatically by the script Abilint (TD).
2513 !Do not modify the following lines by hand.
2514 #undef ABI_FUNC
2515 #define ABI_FUNC 'setup_sigma'
2516 !End of the abilint section
2517 
2518  implicit none
2519 
2520 !Arguments ------------------------------------
2521 !scalars
2522  integer,intent(in) :: comm
2523  character(len=6),intent(in) :: codvsn
2524  character(len=*),intent(in) :: wfk_fname
2525  type(Datafiles_type),intent(in) :: Dtfil
2526  type(Dataset_type),intent(inout) :: Dtset
2527  type(Pseudopotential_type),intent(in) :: Psps
2528  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw)
2529  type(sigparams_t),intent(out) :: Sigp
2530  type(Epsilonm1_results),intent(out) :: Er
2531  type(ebands_t),intent(out) :: KS_BSt
2532  type(kmesh_t),intent(out) :: Kmesh,Qmesh
2533  type(crystal_t),intent(out) :: Cryst
2534  type(gsphere_t),intent(out) :: Gsph_Max,Gsph_x,Gsph_c
2535  type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out
2536  type(vcoul_t),intent(out) :: Vcp
2537 !arrays
2538  integer,intent(in) :: ngfftf(18)
2539  integer,intent(out) :: gwc_ngfft(18),gwx_ngfft(18)
2540  real(dp),intent(in) :: acell(3),rprim(3,3)
2541 
2542 !Local variables-------------------------------
2543 !scalars
2544  integer,parameter :: pertcase0=0,master=0
2545  integer :: bantot,enforce_sym,gwcalctyp,ib,ibtot,icutcoul_eff,ii,ikcalc,ikibz,io,isppol,itypat,jj,method
2546  integer :: mod10,mqmem,mband,ng_kss,nsheps,ikcalc2bz,ierr,gap_err,ng
2547  integer :: gwc_nfftot,gwx_nfftot,nqlwl,test_npwkss,my_rank,nprocs,ik,nk_found,ifo,timrev,usefock_ixc
2548  integer :: iqbz,isym,iq_ibz,itim,ic,pinv,ig1,ng_sigx,spin,gw_qprange,ivcoul_init,nvcoul_init,xclevel_ixc
2549  real(dp),parameter :: OMEGASIMIN=0.01d0,tol_enediff=0.001_dp*eV_Ha
2550  real(dp) :: domegas,domegasi,ucvol,rcut
2551  logical,parameter :: linear_imag_mesh=.TRUE.
2552  logical :: ltest,remove_inv,changed,found
2553  character(len=500) :: msg
2554  character(len=fnlen) :: fname,fcore,string
2555  type(wvl_internal_type) :: wvl
2556  type(gaps_t) :: gaps
2557 !arrays
2558  integer :: ng0sh_opt(3),G0(3),q_umklp(3),kpos(6)
2559  integer,allocatable :: npwarr(:),val_indeces(:,:)
2560  integer,pointer :: gvec_kss(:,:),gsphere_sigx_p(:,:)
2561  integer,pointer :: test_gvec_kss(:,:)
2562  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3),q_bz(3),gamma_point(3,1)
2563  real(dp),pointer :: energies_p(:,:,:)
2564  real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:)
2565  type(Pawrhoij_type),allocatable :: Pawrhoij(:)
2566  type(vcoul_t) :: Vcp_ks
2567 
2568 ! *************************************************************************
2569 
2570  DBG_ENTER('COLL')
2571 
2572  ! Check for calculations that are not implemented
2573  ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1))
2574  ABI_CHECK(ltest,'Dtset%nband(:) must be constant')
2575 
2576  my_rank = xmpi_comm_rank(comm); nprocs  = xmpi_comm_size(comm)
2577 
2578  ! Basic parameters
2579  Sigp%ppmodel    = Dtset%ppmodel
2580  Sigp%gwcalctyp  = Dtset%gwcalctyp
2581  Sigp%nbnds      = Dtset%nband(1)
2582  Sigp%symsigma   = Dtset%symsigma
2583  Sigp%zcut       = Dtset%zcut
2584  Sigp%mbpt_sciss = Dtset%mbpt_sciss
2585 
2586  timrev=  2 ! This information is not reported in the header
2587             ! 1 => do not use time-reversal symmetry
2588             ! 2 => take advantage of time-reversal symmetry
2589  if (any(dtset%kptopt == [3, 4])) timrev = 1
2590 
2591  ! === For HF, SEX or COHSEX use Hybertsen-Louie PPM (only $\omega=0$) ===
2592  ! * Use fake screening for HF.
2593  ! FIXME Why, we should not redefine Sigp%ppmodel
2594  gwcalctyp=Sigp%gwcalctyp
2595  mod10 =MOD(Sigp%gwcalctyp,10)
2596  if (mod10==5.or.mod10==6.or.mod10==7) Sigp%ppmodel=2
2597  if (mod10<5.and.MOD(Sigp%gwcalctyp,1)/=1) then ! * One shot GW (PPM or contour deformation).
2598    if (Dtset%nomegasrd==1) then ! avoid division by zero!
2599      Sigp%nomegasrd  =1
2600      Sigp%maxomega4sd=zero
2601      Sigp%deltae     =zero
2602    else
2603      Sigp%nomegasrd   = Dtset%nomegasrd
2604      Sigp%maxomega4sd = Dtset%omegasrdmax
2605      Sigp%deltae     = (2*Sigp%maxomega4sd)/(Sigp%nomegasrd-1)
2606    endif
2607  else
2608    ! For AC no need to evaluate derivative by finite differences.
2609    Sigp%nomegasrd  =1
2610    Sigp%maxomega4sd=zero
2611    Sigp%deltae     =zero
2612  end if
2613 
2614  ! For analytic continuation define the number of imaginary frequencies for Sigma
2615  ! Tests show than more than 12 freqs in the Pade approximant worsen the results!
2616  Sigp%nomegasi=0
2617 
2618  if (mod10==1) then
2619    Sigp%nomegasi  =Dtset%nomegasi
2620    Sigp%omegasimax=Dtset%omegasimax
2621    Sigp%omegasimin=OMEGASIMIN
2622    write(msg,'(4a,i3,2(2a,f8.3),a)')ch10,&
2623 &    ' Parameters for analytic continuation : ',ch10,&
2624 &    '  number of imaginary frequencies for sigma =  ',Sigp%nomegasi,ch10,&
2625 &    '  min frequency for sigma on imag axis [eV] =  ',Sigp%omegasimin*Ha_eV,ch10,&
2626 &    '  max frequency for sigma on imag axis [eV] =  ',Sigp%omegasimax*Ha_eV,ch10
2627    call wrtout(std_out,msg,'COLL')
2628 
2629    !TODO this should not be done here but in init_sigma_t
2630    ABI_MALLOC(Sigp%omegasi,(Sigp%nomegasi))
2631 
2632    if (linear_imag_mesh) then  ! * Linear mesh along the imaginary axis.
2633      domegasi=Sigp%omegasimax/(Sigp%nomegasi-1)
2634      do io=1,Sigp%nomegasi
2635        Sigp%omegasi(io)=CMPLX(zero,(io-1)*domegasi)
2636      end do
2637    else ! * Logarithmic mesh along the imaginary axis.
2638      MSG_ERROR("AC + log mesh not implemented")
2639      !domegasi=(Sigp%omegasimax/Sigp%omegasimin)**(one/(Sigp%nomegasi-1))
2640      !Sigp%omegasi(1)=czero; ldi=domegasi
2641      !do io=2,Sigp%nomegasi
2642      ! omega(io)=CMPLX(zero,ldi*Sigp%omegasimin)
2643      ! Sigp%omegasi(io)=ldi*domegasi
2644      !end do
2645    end if
2646 
2647    write(msg,'(4a)')ch10,&
2648 &    ' setup_sigma : calculating Sigma(iw)',&
2649 &    ' at imaginary frequencies [eV] (Fermi Level set to 0) ',ch10
2650    call wrtout(std_out,msg,'COLL')
2651    call wrtout(ab_out,msg,'COLL')
2652    do io=1,Sigp%nomegasi
2653      write(msg,'(2(f10.3,2x))')Sigp%omegasi(io)*Ha_eV
2654      call wrtout(std_out,msg,'COLL')
2655      call wrtout(ab_out,msg,'COLL')
2656    end do
2657 
2658    ltest=(Sigp%omegasimax>0.1d-4.and.Sigp%nomegasi>0)
2659    ABI_CHECK(ltest,'Wrong value of omegasimax or nomegasi')
2660    if (Sigp%gwcalctyp/=1) then ! only one shot GW is allowed for AC.
2661      MSG_ERROR("SC-GW with analytic continuation is not coded")
2662    end if
2663  end if
2664 
2665  if (Sigp%symsigma/=0.and.gwcalctyp>=20) then
2666    MSG_WARNING("SC-GW with symmetries is still under development. Use at your own risk!")
2667  end if
2668 
2669  ! Setup parameters for Spectral function.
2670  if (Dtset%gw_customnfreqsp/=0) then
2671    Sigp%nomegasr = Dtset%gw_customnfreqsp
2672    MSG_WARNING('Custom grid for spectral function specified. Assuming experienced user.')
2673    if (Dtset%gw_customnfreqsp/=0) then
2674      Dtset%nfreqsp = Dtset%gw_customnfreqsp
2675      MSG_WARNING('nfreqsp has been set to the same number as gw_customnfreqsp')
2676    end if
2677  else
2678    Sigp%nomegasr  =Dtset%nfreqsp
2679    Sigp%minomega_r=Dtset%freqspmin
2680    Sigp%maxomega_r=Dtset%freqspmax
2681  end if
2682 
2683  if (Sigp%nomegasr>0) then
2684    if (Dtset%gw_customnfreqsp==0) then
2685      ! Check
2686      if (Sigp%minomega_r >= Sigp%maxomega_r) then
2687        MSG_ERROR('freqspmin must be smaller than freqspmax!')
2688      end if
2689      if(Sigp%nomegasr==1) then
2690       domegas=0.d0
2691      else
2692       domegas=(Sigp%maxomega_r-Sigp%minomega_r)/(Sigp%nomegasr-1)
2693      endif
2694      !TODO this should be moved to Sr% and done in init_sigma_t
2695      ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr))
2696      do io=1,Sigp%nomegasr
2697        Sigp%omega_r(io) = CMPLX(Sigp%minomega_r + domegas*(io-1),zero)
2698      end do
2699      write(msg,'(4a,i8,3(2a,f8.3),a)')ch10,&
2700 &      ' Parameters for the calculation of the spectral function : ',ch10,&
2701 &      '  Number of points    = ',Sigp%nomegasr,ch10,&
2702 &      '  Min frequency  [eV] = ',Sigp%minomega_r*Ha_eV,ch10,&
2703 &      '  Max frequency  [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,&
2704 &      '  Frequency step [eV] = ',domegas*Ha_eV,ch10
2705      call wrtout(std_out,msg,'COLL')
2706    else
2707      Sigp%minomega_r = MINVAL(Dtset%gw_freqsp(:))
2708      Sigp%maxomega_r = MAXVAL(Dtset%gw_freqsp(:))
2709      !TODO this should be moved to Sr% and done in init_sigma_t
2710      ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr))
2711      do io=1,Sigp%nomegasr
2712        Sigp%omega_r(io) = CMPLX(Dtset%gw_freqsp(io),zero)
2713      end do
2714      write(msg,'(4a,i8,2(2a,f8.3),3a)')ch10,&
2715 &      ' Parameters for the calculation of the spectral function : ',ch10,&
2716 &      '  Number of points    = ',Sigp%nomegasr,ch10,&
2717 &      '  Min frequency  [eV] = ',Sigp%minomega_r*Ha_eV,ch10,&
2718 &      '  Max frequency  [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,&
2719 &      '  A custom set of frequencies is used! See the input file for values.',ch10
2720      call wrtout(std_out,msg,'COLL')
2721    end if
2722  else
2723    !In indefo all these quantities are set to zero
2724    !Sigp%nomegasr=1
2725    !allocate(Sigp%omega_r(Sigp%nomegasr))
2726    !Sigp%omega_r(1)=0
2727  end if
2728 
2729  ! Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume
2730  call mkrdim(acell,rprim,rprimd)
2731  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2732 
2733  Sigp%npwwfn=Dtset%npwwfn
2734  Sigp%npwx  =Dtset%npwsigx
2735 
2736  ! Read parameters of the WFK, verifify them and retrieve all G-vectors.
2737  call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm)
2738  mband = MAXVAL(Hdr_wfk%nband)
2739 
2740  remove_inv = .FALSE.
2741  call hdr_vs_dtset(Hdr_wfk,Dtset)
2742 
2743  test_npwkss = 0
2744  call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,&
2745 &  gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr)
2746  ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss")
2747 
2748  ABI_MALLOC(gvec_kss,(3,test_npwkss))
2749  gvec_kss = test_gvec_kss
2750  ng_kss = test_npwkss
2751 
2752  ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2))
2753  ierr = 0
2754  do ig1=1,ng
2755    if (ANY(gvec_kss(:,ig1)/=test_gvec_kss(:,ig1))) then
2756      ierr=ierr+1
2757      write(std_out,*)" gvec_kss ",ig1,"/",ng,gvec_kss(:,ig1),test_gvec_kss(:,ig1)
2758    end if
2759  end do
2760  ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss")
2761 
2762  ABI_FREE(test_gvec_kss)
2763 
2764  ! Get important dimensions from the WFK header
2765  Sigp%nsppol =Hdr_wfk%nsppol
2766  Sigp%nspinor=Hdr_wfk%nspinor
2767  Sigp%nsig_ab=Hdr_wfk%nspinor**2  ! TODO Is it useful calculating only diagonal terms?
2768 
2769  if (Sigp%nbnds>mband) then
2770    Sigp%nbnds    =mband
2771    Dtset%nband(:)=mband
2772    Dtset%mband   =MAXVAL(Dtset%nband)
2773    write(msg,'(3a,i4,a)')&
2774 &    'Number of bands found less then required',ch10,&
2775 &    'calculation will proceed with nbnds = ',mband,ch10
2776    MSG_WARNING(msg)
2777  end if
2778 
2779  ! Check input
2780  if (Sigp%ppmodel==3.or.Sigp%ppmodel==4) then
2781    if (gwcalctyp>=10) then
2782      write(msg,'(a,i3,a)')' The ppmodel chosen and gwcalctyp ',Dtset%gwcalctyp,' are not compatible. '
2783      MSG_ERROR(msg)
2784    end if
2785    if (Sigp%nspinor==2) then
2786      write(msg,'(a,i3,a)')' The ppmodel chosen and nspinor ',Sigp%nspinor,' are not compatible. '
2787      MSG_ERROR(msg)
2788    end if
2789  end if
2790 
2791  ! Create crystal_t data type
2792  call crystal_from_hdr(Cryst,Hdr_wfk,timrev,remove_inv)
2793  call crystal_print(Cryst)
2794 
2795  if (Sigp%npwwfn>ng_kss) then ! cannot use more G"s for the wfs than those stored on file
2796    Sigp%npwwfn  =ng_kss
2797    Dtset%npwwfn =ng_kss
2798    write(msg,'(2a,(a,i8,a))')&
2799 &    'Number of G-vectors for WFS found in the KSS file is less than required',ch10,&
2800 &    'calculation will proceed with npwwfn  = ',Sigp%npwwfn,ch10
2801    MSG_WARNING(msg)
2802  end if
2803 
2804  if (Sigp%npwx>ng_kss) then
2805    ! Have to recalcuate the (large) sphere for Sigma_x.
2806    pinv=1; if (remove_inv.and.Cryst%timrev==2) pinv=-1
2807    gamma_point(:,1) = (/zero,zero,zero/); nullify(gsphere_sigx_p)
2808 
2809    call merge_and_sort_kg(1,gamma_point,Dtset%ecutsigx,Cryst%nsym,pinv,Cryst%symrel,&
2810 &   Cryst%gprimd,gsphere_sigx_p,Dtset%prtvol)
2811 
2812    ng_sigx=SIZE(gsphere_sigx_p,DIM=2)
2813    Sigp%npwx     = ng_sigx
2814    Dtset%npwsigx = ng_sigx
2815 
2816    write(msg,'(2a,(a,i8,a))')&
2817 &    'Number of G-vectors for Sigma_x found in the KSS file is less than required',ch10,&
2818 &    'calculation will proceed with npwsigx = ',Sigp%npwx,ch10
2819    MSG_WARNING(msg)
2820 
2821    ltest = (Sigp%npwx >= ng_kss)
2822    ABI_CHECK(ltest,"Sigp%npwx<ng_kss!")
2823 
2824    ! * Fill gvec_kss with larger sphere.
2825    ABI_FREE(gvec_kss)
2826    ABI_MALLOC(gvec_kss,(3,Sigp%npwx))
2827    gvec_kss = gsphere_sigx_p
2828    ABI_FREE(gsphere_sigx_p)
2829  end if
2830 
2831  ! Set up of the k-points and tables in the whole BZ ===
2832  ! TODO Recheck symmorphy and inversion
2833  call kmesh_init(Kmesh,Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.FALSE.)
2834  !call kmesh_init(Kmesh,Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.TRUE.)
2835 
2836  ! Some required information are not filled up inside kmesh_init
2837  ! So doing it here, even though it is not clean
2838  Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:)
2839  Kmesh%nshift        =Dtset%nshiftk
2840  ABI_ALLOCATE(Kmesh%shift,(3,Kmesh%nshift))
2841  Kmesh%shift(:,:)    =Dtset%shiftk(:,1:Dtset%nshiftk)
2842 
2843  call kmesh_print(Kmesh,"K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL")
2844  call kmesh_print(Kmesh,"K-mesh for the wavefunctions",ab_out, 0,           "COLL")
2845 
2846  ! === Initialize the band structure datatype ===
2847  ! * Copy WFK energies and occupations up to Sigp%nbnds==Dtset%nband(:)
2848  bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol))
2849  ABI_MALLOC(doccde,(bantot))
2850  ABI_MALLOC(eigen,(bantot))
2851  ABI_MALLOC(occfact,(bantot))
2852  doccde(:)=zero; eigen(:)=zero; occfact(:)=zero
2853 
2854  jj=0; ibtot=0
2855  do isppol=1,Dtset%nsppol
2856    do ikibz=1,Dtset%nkpt
2857      do ib=1,Hdr_wfk%nband(ikibz+(isppol-1)*Dtset%nkpt)
2858        ibtot=ibtot+1
2859        if (ib<=Sigp%nbnds) then
2860          jj=jj+1
2861          occfact(jj)=Hdr_wfk%occ(ibtot)
2862          eigen  (jj)=energies_p(ib,ikibz,isppol)
2863        end if
2864      end do
2865    end do
2866  end do
2867  ABI_FREE(energies_p)
2868 
2869  ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of
2870  ! symmetry operations in the old GW code (symmorphy and inversion)
2871  ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6))
2872  ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt')
2873 
2874  ABI_MALLOC(npwarr,(Dtset%nkpt))
2875  npwarr(:)=Sigp%npwwfn
2876 
2877  call ebands_init(bantot,KS_BSt,Dtset%nelect,doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,&
2878 &  Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,&
2879 &  dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig,&
2880 &  dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
2881 
2882  ABI_FREE(doccde)
2883  ABI_FREE(eigen)
2884  ABI_FREE(npwarr)
2885 
2886  ! Calculate KS occupation numbers and ks_vbk(nkibz,nsppol) ====
2887  ! ks_vbk gives the (valence|last Fermi band) index for each k and spin.
2888  ! spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals
2889  call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=0)
2890 
2891  gap_err = get_gaps(KS_BSt, gaps)
2892  call gaps_print(gaps, unit=std_out)
2893  call ebands_report_gap(KS_BSt, unit=std_out)
2894 
2895  ABI_MALLOC(val_indeces,(KS_BSt%nkpt,KS_BSt%nsppol))
2896  val_indeces = get_valence_idx(KS_BSt)
2897 
2898  ! Create Sigma header
2899  ! TODO Fix problems with symmorphy and k-points
2900  call hdr_init(KS_BSt,codvsn,Dtset,Hdr_out,Pawtab,pertcase0,Psps,wvl)
2901 
2902  ! Get Pawrhoij from the header of the WFK file
2903  ABI_DT_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw))
2904  if (Dtset%usepaw==1) then
2905    call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab)
2906    call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij)
2907  end if
2908  call hdr_update(hdr_out,bantot,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
2909 
2910  ABI_FREE(occfact)
2911  call pawrhoij_free(Pawrhoij)
2912  ABI_DT_FREE(Pawrhoij)
2913 
2914  ! ===========================================================
2915  ! ==== Setup of k-points and bands for the GW corrections ====
2916  ! ===========================================================
2917  ! * maxbdgw and minbdgw are the Max and min band index for GW corrections over k-points.
2918  !   They are used to dimension the wavefunctions and to calculate the matrix elements.
2919  !
2920  if (dtset%nkptgw==0) then
2921    ! Use qp_range to select the interesting k-points and the corresponing bands.
2922    !
2923    !    0 --> Compute the QP corrections only for the fundamental and the optical gap.
2924    ! +num --> Compute the QP corrections for all the k-points in the irreducible zone and include `num`
2925    !           bands above and below the Fermi level.
2926    ! -num --> Compute the QP corrections for all the k-points in the irreducible zone.
2927    !          Include all occupied states and `num` empty states.
2928 
2929    call wrtout(std_out, "nkptgw == 0 ==> Automatic selection of k-points and bands for the corrections.")
2930    if (gap_err /=0 .and. dtset%gw_qprange==0) then
2931      msg = "Problem while computing the fundamental and optical gap (likely metal). Will replace gw_qprange=0 with gw_qprange=1"
2932      MSG_WARNING(msg)
2933      dtset%gw_qprange = 1
2934    end if
2935    gw_qprange = dtset%gw_qprange
2936 
2937    if (dtset%ucrpa>0) then
2938      dtset%nkptgw=Kmesh%nbz
2939      Sigp%nkptgw =dtset%nkptgw
2940      ABI_MALLOC(Sigp%kptgw,(3,Sigp%nkptgw))
2941      ABI_MALLOC(Sigp%minbnd,(Sigp%nkptgw,Sigp%nsppol))
2942      ABI_MALLOC(Sigp%maxbnd,(Sigp%nkptgw,Sigp%nsppol))
2943      Sigp%kptgw(:,:)=Kmesh%bz(:,:)
2944      Sigp%minbnd=1
2945      Sigp%maxbnd=Sigp%nbnds
2946 
2947    else if (gw_qprange/=0) then
2948      ! Include all the k-points in the IBZ.
2949      ! Note that kptgw == ebands%kptns so we can use a single ik index in the loop over k-points
2950      ! No need to map kptgw onto ebands%kptns.
2951      dtset%nkptgw=Kmesh%nibz
2952      Sigp%nkptgw =dtset%nkptgw
2953      ABI_MALLOC(Sigp%kptgw,(3,Sigp%nkptgw))
2954      ABI_MALLOC(Sigp%minbnd,(Sigp%nkptgw,Sigp%nsppol))
2955      ABI_MALLOC(Sigp%maxbnd,(Sigp%nkptgw,Sigp%nsppol))
2956      Sigp%kptgw(:,:)=Kmesh%ibz(:,:)
2957      Sigp%minbnd=1
2958      Sigp%maxbnd=Sigp%nbnds
2959 
2960      if (gw_qprange>0) then
2961        ! All k-points: Add buffer of bands above and below the Fermi level.
2962        do spin=1,Sigp%nsppol
2963          do ik=1,Sigp%nkptgw
2964            Sigp%minbnd(ik,spin) = MAX(val_indeces(ik,spin) - gw_qprange, 1)
2965            Sigp%maxbnd(ik,spin) = MIN(val_indeces(ik,spin) + gw_qprange + 1, Sigp%nbnds)
2966          end do
2967        end do
2968 
2969      else
2970        ! All k-points: include all occupied states and -gw_qprange empty states.
2971        Sigp%minbnd = 1
2972        do spin=1,Sigp%nsppol
2973          do ik=1,Sigp%nkptgw
2974            Sigp%maxbnd(ik,spin) = MIN(val_indeces(ik,spin) - gw_qprange, Sigp%nbnds)
2975          end do
2976        end do
2977      end if
2978 
2979    else
2980      ! gw_qprange is not specified in the input.
2981      ! Include the optical and the fundamental KS gap.
2982      ! The main problem here is that kptgw and nkptgw do not depend on the spin and therefore
2983      ! we have compute the union of the k-points where the fundamental and the optical gaps are located.
2984      !
2985      ! Find the list of `interesting` kpoints.
2986      ABI_CHECK(gap_err == 0, "gw_qprange 0 cannot be used because I cannot find the gap (gap_err !=0)")
2987      nk_found = 1; kpos(1) = gaps%fo_kpos(1,1)
2988 
2989      do spin=1,Sigp%nsppol
2990        do ifo=1,3
2991          ik = gaps%fo_kpos(ifo, spin)
2992          found = .FALSE.; jj = 0
2993          do while (.not. found .and. jj < nk_found)
2994            jj = jj + 1; found = (kpos(jj) == ik)
2995          end do
2996          if (.not. found) then
2997            nk_found = nk_found + 1; kpos(nk_found) = ik
2998          end if
2999        end do
3000      end do
3001 
3002      ! Now we can define the list of k-points and the bands range.
3003      dtset%nkptgw=nk_found
3004      Sigp%nkptgw =dtset%nkptgw
3005 
3006      ABI_MALLOC(Sigp%kptgw,(3,Sigp%nkptgw))
3007      ABI_MALLOC(Sigp%minbnd,(Sigp%nkptgw,Sigp%nsppol))
3008      ABI_MALLOC(Sigp%maxbnd,(Sigp%nkptgw,Sigp%nsppol))
3009 
3010      do ii=1,Sigp%nkptgw
3011        ik = kpos(ii)
3012        Sigp%kptgw(:,ii)=Kmesh%ibz(:,ik)
3013        do spin=1,Sigp%nsppol
3014          Sigp%minbnd(ii,spin) = val_indeces(ik,spin)
3015          Sigp%maxbnd(ii,spin) = val_indeces(ik,spin) + 1
3016        end do
3017      end do
3018    end if
3019 
3020  else
3021    ! Treat only the k-points and bands specified in the input file.
3022    Sigp%nkptgw=dtset%nkptgw
3023    ABI_MALLOC(Sigp%kptgw,(3,Sigp%nkptgw))
3024    ABI_MALLOC(Sigp%minbnd,(Sigp%nkptgw,Sigp%nsppol))
3025    ABI_MALLOC(Sigp%maxbnd,(Sigp%nkptgw,Sigp%nsppol))
3026 
3027    do spin=1,Sigp%nsppol
3028      Sigp%minbnd(:,spin)=dtset%bdgw(1,:,spin)
3029      Sigp%maxbnd(:,spin)=dtset%bdgw(2,:,spin)
3030    end do
3031 
3032    do ii=1,3
3033      do ikcalc=1,Sigp%nkptgw
3034        Sigp%kptgw(ii,ikcalc)=Dtset%kptgw(ii,ikcalc)
3035      end do
3036    end do
3037 
3038    do spin=1,Sigp%nsppol
3039      do ikcalc=1,Sigp%nkptgw
3040        if (Dtset%bdgw(2,ikcalc,spin)>Sigp%nbnds) then
3041          write(msg,'(a,2i0,2(a,i0),2a,i0)')&
3042 &          "For (k,s) ",ikcalc,spin," bdgw= ",Dtset%bdgw(2,ikcalc,spin), " > nbnds=",Sigp%nbnds,ch10,&
3043 &          "Calculation will continue with bdgw =",Sigp%nbnds
3044          MSG_COMMENT(msg)
3045          Dtset%bdgw(2,ikcalc,spin)=Sigp%nbnds
3046        end if
3047      end do
3048    end do
3049 
3050  end if
3051 
3052  ! Make sure that all the degenerate states are included.
3053  ! * We will have to average the GW corrections over degenerate states if symsigma=1 is used.
3054  ! * KS states belonging to the same irreducible representation should be included in the basis set used for SCGW.
3055  if (Sigp%symsigma/=0 .or. gwcalctyp>=10) then
3056    do isppol=1,Sigp%nsppol
3057      do ikcalc=1,Sigp%nkptgw
3058 
3059        if (has_IBZ_item(Kmesh,Sigp%kptgw(:,ikcalc),ikibz,G0)) then
3060          call enclose_degbands(KS_BSt,ikibz,isppol,Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol),changed,tol_enediff)
3061          if (changed) then
3062            write(msg,'(2(a,i0),2a,2(1x,i0))')&
3063 &            "Not all the degenerate states at ikcalc= ",ikcalc,", spin= ",isppol,ch10,&
3064 &            "were included in the bdgw set. bdgw has been changed to: ",Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol)
3065            MSG_COMMENT(msg)
3066          end if
3067        else
3068          MSG_ERROR(sjoin('k-point', ktoa(Sigp%kptgw(:,ikcalc)), 'not in IBZ'))
3069        end if
3070 
3071      end do
3072    end do
3073  end if
3074 
3075  !if (.not. associated(Dtset%bdgw)) then
3076  !  ABI_MALLOC(Dtset%bdgw, (2,Sigp%nkptgw,Sigp%nsppol))
3077  !end if
3078  !do spin=1,Sigp%nsppol
3079  !  Dtset%bdgw(1,:,spin) = Sigp%minbnd(:,spin)
3080  !  Dtset%bdgw(2,:,spin) = Sigp%maxbnd(:,spin)
3081  !end do
3082 
3083  Sigp%minbdgw=MINVAL(Sigp%minbnd)
3084  Sigp%maxbdgw=MAXVAL(Sigp%maxbnd)
3085 
3086  ABI_MALLOC(Sigp%kptgw2bz,(Sigp%nkptgw))
3087  !
3088  !=== Check if the k-points are in the BZ ===
3089  !FB TODO Honestly the code is not able to treat k-points, which are not in the IBZ.
3090  !This extension should require to change the code in different places.
3091  !Therefore, one should by now prevent the user from calculating sigma for a k-point not in the IBZ.
3092 
3093  do ikcalc=1,Sigp%nkptgw
3094    if (has_BZ_item(Kmesh,Sigp%kptgw(:,ikcalc),ikcalc2bz,G0)) then
3095      !found = has_IBZ_item(Kmesh,Sigp%kptgw(:,ikcalc),ikcalc2bz,G0)
3096      Sigp%kptgw2bz(ikcalc) = ikcalc2bz
3097    else
3098      MSG_ERROR(sjoin('k-point:', ktoa(Sigp%kptgw(:,ikcalc)), 'not in the kbz set'))
3099    end if
3100  end do
3101 
3102  ! Check if there are duplicated k-point in Sigp%
3103  do ii=1,Sigp%nkptgw
3104    do jj=ii+1,Sigp%nkptgw
3105      if (isamek(Sigp%kptgw(:,ii),Sigp%kptgw(:,jj),G0)) then
3106        write(msg,'(5a)')&
3107 &        'kptgw contains duplicated k-points. This is not allowed since ',ch10,&
3108 &        'the QP corrections for this k-point will be calculated more than once. ',ch10,&
3109 &        'Check your input file. '
3110        MSG_ERROR(msg)
3111      end if
3112    end do
3113  end do
3114  !
3115  ! Warn the user if SCGW run and not all the k-points are included.
3116  if (gwcalctyp>=10 .and. Sigp%nkptgw/=Hdr_wfk%nkpt) then
3117    write(msg,'(3a,2(a,i0),2a)')ch10,&
3118 &    " COMMENT: In a self-consistent GW run, the QP corrections should be calculated for all the k-points of the KSS file ",ch10,&
3119 &    " but nkptgw= ",Sigp%nkptgw," and WFK nkpt= ",Hdr_wfk%nkpt,ch10,&
3120 &    " Assuming expert user. Execution will continue. "
3121    call wrtout(ab_out,msg,"COLL")
3122  end if
3123 
3124  ! Setup of the table used in the case of SCGW on wavefunctions to reduce the number
3125  ! of elements <i,kgw,s|\Sigma|j,kgw,s> that have to be calculated. No use of symmetries, except for Hermiticity.
3126  call sigma_tables(Sigp,Kmesh)
3127 
3128  ! === Read external file and initialize basic dimension of Er% ===
3129  ! TODO use mqmem as input variable instead of gwmem
3130 
3131  ! === If required, use a matrix for $\Sigma_c$ which is smaller than that stored on file ===
3132  ! * By default the entire matrix is read and used,
3133  ! * Define consistently npweps and ecuteps for \Sigma_c according the input
3134  if (Dtset%npweps>0.or.Dtset%ecuteps>0) then
3135    ! This should not happen: the Dtset array should not be modified after having been initialized.
3136    if (Dtset%npweps>0) Dtset%ecuteps=zero
3137    nsheps=0
3138    call setshells(Dtset%ecuteps,Dtset%npweps,nsheps,Dtset%nsym,gmet,gprimd,Dtset%symrel,'eps',ucvol)
3139  end if
3140 
3141  mqmem=0; if (Dtset%gwmem/10==1) mqmem=1
3142 
3143  if (Dtset%getscr/=0.or.Dtset%irdscr/=0) then
3144    fname=Dtfil%fnameabi_scr
3145  else if (Dtset%getsuscep/=0.or.Dtset%irdsuscep/=0) then
3146    fname=Dtfil%fnameabi_sus
3147  else
3148    fname=Dtfil%fnameabi_scr
3149    !FIXME this has to be cleaned, in tgw2_3 Dtset%get* and Dtset%ird* are  not defined
3150    !MSG_ERROR("getsuscep or irdsuscep are not defined")
3151  end if
3152  !
3153  ! === Setup of q-mesh in the whole BZ ===
3154  ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed,
3155  !   epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed)
3156  !
3157  if (sigma_needs_w(Sigp)) then
3158    if (.not. file_exists(fname)) then
3159      fname = nctk_ncify(fname)
3160      MSG_COMMENT(sjoin("File not found. Will try netcdf file:", fname))
3161    end if
3162 
3163    call init_Er_from_file(Er,fname,mqmem,Dtset%npweps,comm)
3164 
3165    Sigp%npwc=Er%npwe
3166    if (Sigp%npwc>Sigp%npwx) then
3167      Sigp%npwc=Sigp%npwx
3168      MSG_COMMENT("Found npw_correlation > npw_exchange, Imposing npwc=npwx")
3169      ! There is a good reason for doing so, see csigme.F90 and the size of the arrays
3170      ! rhotwgp and rhotwgp: we need to define a max size and we opt for Sigp%npwx.
3171    end if
3172    Er%npwe=Sigp%npwc
3173    Dtset%npweps=Er%npwe
3174    call kmesh_init(Qmesh,Cryst,Er%nqibz,Er%qibz,Dtset%kptopt)
3175 
3176  else
3177    Er%npwe     =1
3178    Sigp%npwc   =1
3179    Dtset%npweps=1
3180    call find_qmesh(Qmesh,Cryst,Kmesh)
3181    ABI_MALLOC(Er%gvec,(3,1))
3182    Er%gvec(:,1) = (/0,0,0/)
3183  end if
3184 
3185  call kmesh_print(Qmesh,"Q-mesh for screening function",std_out,Dtset%prtvol,"COLL")
3186  call kmesh_print(Qmesh,"Q-mesh for screening function",ab_out ,0           ,"COLL")
3187 
3188  do iqbz=1,Qmesh%nbz
3189    call get_BZ_item(Qmesh,iqbz,q_bz,iq_ibz,isym,itim,umklp=q_umklp)
3190 
3191    if (ANY(q_umklp/=0)) then
3192      sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz))
3193      write(std_out,*) sq,Qmesh%bz(:,iqbz)
3194      write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')&
3195 &      'qpoint ',Qmesh%bz(:,iqbz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,&
3196 &      'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,&
3197 &      'however a non zero umklapp G_o vector is required and this is not yet allowed'
3198      MSG_ERROR(msg)
3199    end if
3200  end do
3201  !
3202  ! === Find optimal value for G-sphere enlargment due to oscillator matrix elements ===
3203  ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true size bz is buggy
3204  ! * -one is used because we loop over all the possibile differences, unlike screening
3205 
3206  call get_ng0sh(Sigp%nkptgw,Sigp%kptgw,Kmesh%nbz,Kmesh%bz,Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt)
3207  call wrtout(std_out, sjoin(' Optimal value for ng0sh ', ltoa(ng0sh_opt)), "COLL")
3208  Sigp%mG0=ng0sh_opt
3209 
3210  ! G-sphere for W and Sigma_c is initialized from the SCR file.
3211  call gsph_init(Gsph_c,Cryst,Er%npwe,gvec=Er%gvec)
3212  call gsph_init(Gsph_x,Cryst,Sigp%npwx,gvec=gvec_kss)
3213 
3214  Sigp%ecuteps = Gsph_c%ecut
3215  Dtset%ecuteps = Sigp%ecuteps
3216 
3217  ! === Make biggest G-sphere of Sigp%npwvec vectors ===
3218  Sigp%npwvec=MAX(Sigp%npwwfn,Sigp%npwx)
3219  call gsph_init(Gsph_Max,Cryst,Sigp%npwvec,gvec=gvec_kss)
3220 
3221 !BEGINDEBUG
3222  ! Make sure that the two G-spheres are equivalent.
3223  ierr=0
3224  if (sigma_needs_w(Sigp)) then
3225    ng = MIN(SIZE(Gsph_c%gvec,DIM=2),SIZE(gvec_kss,DIM=2))
3226    do ig1=1,ng
3227      if (ANY(Gsph_c%gvec(:,ig1)/=gvec_kss(:,ig1))) then
3228        ierr=ierr+1
3229        write(std_out,*)" Gsph_c, gvec_kss ",ig1,"/",ng,Gsph_c%gvec(:,ig1),gvec_kss(:,ig1)
3230      end if
3231    end do
3232    ABI_CHECK(ierr==0,"Mismatch between Gsph_c and gvec_kss")
3233  end if
3234  ierr=0
3235  ng = MIN(SIZE(Gsph_x%gvec,DIM=2),SIZE(gvec_kss,DIM=2))
3236  do ig1=1,ng
3237    if (ANY(Gsph_x%gvec(:,ig1)/=gvec_kss(:,ig1))) then
3238      ierr=ierr+1
3239      write(std_out,*)" Gsph_x, gvec_kss ",ig1,"/",ng,Gsph_x%gvec(:,ig1),gvec_kss(:,ig1)
3240    end if
3241  end do
3242  ABI_CHECK(ierr==0,"Mismatch between Gsph_x and gvec_kss")
3243 !ENDDEBUG
3244 
3245  ABI_FREE(gvec_kss)
3246  !
3247  ! === Get Fourier components of the Coulombian for all q-points in the IBZ ===
3248  ! * If required, use a cutoff in the interaction
3249  ! * Pcv%vc_sqrt contains Vc^{-1/2}
3250  ! * Setup also the analytical calculation of the q->0 component
3251  ! FIXME recheck ngfftf since I got different charge outside the cutoff region
3252 
3253  if (Dtset%gw_nqlwl==0) then
3254    nqlwl=1
3255    ABI_MALLOC(qlwl,(3,nqlwl))
3256    qlwl(:,1)= GW_Q0_DEFAULT
3257  else
3258    nqlwl=Dtset%gw_nqlwl
3259    ABI_MALLOC(qlwl,(3,nqlwl))
3260    qlwl(:,:)=Dtset%gw_qlwl(:,1:nqlwl)
3261  end if
3262 
3263 !The Coulomb interaction used here might have two terms :
3264 !the first term generates the usual sigma self-energy, but possibly, one should subtract
3265 !from it the Coulomb interaction already present in the Kohn-Sham basis,
3266 !if the usefock associated to ixc is one.
3267 !The latter excludes (in the present implementation) mod(Dtset%gwcalctyp,10)==5
3268  nvcoul_init=1
3269  call get_xclevel(Dtset%ixc,xclevel_ixc,usefock_ixc)
3270  if(usefock_ixc==1)then
3271    nvcoul_init=2
3272    if(mod(Dtset%gwcalctyp,10)==5)then
3273      write(msg,'(4a,i3,a,i3,4a,i5)')ch10,&
3274 &     ' The starting wavefunctions were obtained from self-consistent calculations in the planewave basis set',ch10,&
3275 &     ' with ixc = ',Dtset%ixc,' associated with usefock =',usefock_ixc,ch10,&
3276 &     ' In this case, the present implementation does not allow that the self-energy for sigma corresponds to',ch10,&
3277 &     '  mod(gwcalctyp,10)==5, while your gwcalctyp= ',Dtset%gwcalctyp
3278      MSG_ERROR(msg)
3279    endif
3280  endif
3281 
3282  do ivcoul_init=1,nvcoul_init
3283    rcut = Dtset%rcut
3284    icutcoul_eff=Dtset%icutcoul
3285    Sigp%sigma_mixing=one
3286    if( mod(Dtset%gwcalctyp,10)==5 .or. ivcoul_init==2)then
3287      if(abs(Dtset%hyb_mixing)>tol8)then
3288 !      Warning : the absolute value is needed, because of the singular way used to define the default for this input variable
3289        Sigp%sigma_mixing=abs(Dtset%hyb_mixing)
3290      else if(abs(Dtset%hyb_mixing_sr)>tol8)then
3291        Sigp%sigma_mixing=abs(Dtset%hyb_mixing_sr)
3292        icutcoul_eff=5
3293      endif
3294      if(abs(rcut)<tol6 .and. abs(Dtset%hyb_range_fock)>tol8)rcut=one/Dtset%hyb_range_fock
3295    endif
3296 
3297 !#if 1
3298    if(ivcoul_init==1)then
3299 
3300      if (Gsph_x%ng > Gsph_c%ng) then
3301        call vcoul_init(Vcp,Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3302 &        Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,ngfftf,comm)
3303      else
3304        call vcoul_init(Vcp,Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3305 &        Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,ngfftf,comm)
3306      end if
3307 
3308    else
3309 
3310 !    Use a temporary Vcp_ks to compute the Coulomb interaction already present in the Fock part of the Kohn-Sham Hamiltonian
3311      if (Gsph_x%ng > Gsph_c%ng) then
3312        call vcoul_init(Vcp_ks,Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3313 &        Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,ngfftf,comm)
3314      else
3315        call vcoul_init(Vcp_ks,Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3316 &        Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,ngfftf,comm)
3317      end if
3318 
3319 !    Now compute the residual Coulomb interaction
3320      Vcp%vc_sqrt_resid=sqrt(Vcp%vc_sqrt**2-Sigp%sigma_mixing*Vcp_ks%vc_sqrt**2)
3321      Vcp%i_sz_resid=Vcp%i_sz-Sigp%sigma_mixing*Vcp_ks%i_sz
3322 !    The mixing factor has already been accounted for, so set it back to one
3323      Sigp%sigma_mixing=one
3324      call vcoul_free(Vcp_ks)
3325 
3326 !DEBUG
3327      write(std_out,'(a)')' setup_sigma : the residual Coulomb interaction has been computed'
3328 !ENDDEBUG
3329 
3330    endif
3331 !#else
3332 !   call vcoul_init(Vcp,Gsph_Max,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,ivcoul_init,Dtset%vcutgeo,&
3333 !&    Dtset%ecutsigx,Sigp%npwx,nqlwl,qlwl,ngfftf,comm)
3334 !#endif
3335 
3336  enddo
3337 
3338 #if 0
3339  ! Using the random q for the optical limit is one of the reasons
3340  ! why sigma breaks the initial energy degeneracies.
3341  Vcp%i_sz=zero
3342  Vcp%vc_sqrt(1,:)=czero
3343  Vcp%vcqlwl_sqrt(1,:)=czero
3344 #endif
3345 
3346  ABI_FREE(qlwl)
3347 
3348  Sigp%ecuteps = Dtset%ecuteps
3349  Sigp%ecutwfn = Dtset%ecutwfn
3350  Sigp%ecutsigx = Dtset%ecutsigx
3351 
3352  ! === Setup of the FFT mesh for the oscilator strengths ===
3353  ! * gwc_ngfft(7:18)==Dtset%ngfft(7:18) which is initialized before entering screening.
3354  ! * Here we redefine gwc_ngfft(1:6) according to the following options :
3355  !
3356  ! method==0 --> FFT grid read from fft.in (debugging purpose)
3357  ! method==1 --> Normal FFT mesh
3358  ! method==2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
3359  ! method==3 --> Doubled FFT grid, same as the the FFT for the density,
3360  !
3361  ! enforce_sym==1 ==> Enforce a FFT mesh compatible with all the symmetry operation and FFT library
3362  ! enforce_sym==0 ==> Find the smallest FFT grid compatbile with the library, do not care about symmetries
3363  !
3364  gwc_ngfft(1:18)=Dtset%ngfft(1:18)
3365  gwx_ngfft(1:18)=Dtset%ngfft(1:18)
3366 
3367  method=2
3368  if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0
3369  if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1
3370  if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2
3371  if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3
3372  enforce_sym=MOD(Dtset%fftgw,10)
3373 
3374  ! FFT mesh for sigma_x.
3375  call setmesh(gmet,Gsph_Max%gvec,gwx_ngfft,Sigp%npwvec,Sigp%npwx,Sigp%npwwfn,&
3376 &  gwx_nfftot,method,Sigp%mG0,Cryst,enforce_sym)
3377 
3378  ! FFT mesh for sigma_c.
3379  call setmesh(gmet,Gsph_Max%gvec,gwc_ngfft,Sigp%npwvec,Er%npwe,Sigp%npwwfn,&
3380 &  gwc_nfftot,method,Sigp%mG0,Cryst,enforce_sym,unit=dev_null)
3381 
3382  !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Sigp%mG0,enforce_sym,gwx_ngfft,gwx_nfftot)
3383  !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Sigp%mG0,enforce_sym,gwc_ngfft,gwc_nfftot)
3384 
3385  ! ======================================================================
3386  ! ==== Check for presence of files with core orbitals, for PAW only ====
3387  ! ======================================================================
3388  Sigp%use_sigxcore=0
3389  if (Dtset%usepaw==1.and.Dtset%gw_sigxcore==1) then
3390    ii = 0
3391    do itypat=1,Cryst%ntypat
3392      string = Psps%filpsp(itypat)
3393      fcore = "CORE_"//TRIM(basename(string))
3394      ic = INDEX (TRIM(string), "/" , back=.TRUE.)
3395      if (ic>0 .and. ic<LEN_TRIM(string)) then
3396        ! string defines a path, prepend path to fcore
3397        fcore = Psps%filpsp(itypat)(1:ic)//TRIM(fcore)
3398      end if
3399      if (file_exists(fcore)) then
3400        ii = ii+1
3401      else
3402        MSG_WARNING(sjoin("HF decoupling is required but cannot find file:", fcore))
3403      end if
3404    end do
3405 
3406    Sigp%use_sigxcore=1
3407    if (ii/=Cryst%ntypat) then
3408      MSG_ERROR("Files with core orbitals not found")
3409    end if
3410  end if ! PAW+HF decoupling
3411  !
3412  ! ==============================
3413  ! ==== Extrapolar technique ====
3414  ! ==============================
3415  Sigp%gwcomp   = Dtset%gwcomp
3416  Sigp%gwencomp = Dtset%gwencomp
3417 
3418  if (Sigp%gwcomp==1) then
3419    write(msg,'(6a,e11.4,a)')ch10,&
3420 &    'Using the extrapolar approximation to accelerate convergence',ch10,&
3421 &    'with respect to the number of bands included',ch10,&
3422 &    'with gwencomp: ',Sigp%gwencomp*Ha_eV,' [eV]'
3423    call wrtout(std_out,msg,'COLL')
3424  end if
3425  !
3426  ! ===================================
3427  ! ==== Final compatibility tests ====
3428  ! ===================================
3429  if (ANY(KS_BSt%istwfk/=1)) then
3430    MSG_WARNING('istwfk/=1 is still under development')
3431  end if
3432 
3433  ltest=(KS_BSt%mband==Sigp%nbnds.and.ALL(KS_BSt%nband==Sigp%nbnds))
3434  ABI_CHECK(ltest,'BUG in definition of KS_BSt%nband')
3435 
3436  ! FIXME
3437  if (Dtset%symsigma/=0 .and. Sigp%nomegasr/=0) then
3438    if (idx_spatial_inversion(Cryst) == 0) then
3439      write(msg,'(5a)')' setup_sigma : BUG :',ch10,&
3440 &      'It is not possible to use symsigma/=0 to calculate the spectral function ',ch10,&
3441 &      'when the system does not have the spatial inversion. Please use symsigma=0 '
3442      MSG_WARNING(msg)
3443    end if
3444  end if
3445 
3446  if (mod10==SIG_GW_AC) then
3447    if (Sigp%gwcalctyp/=1) MSG_ERROR("Self-consistency with AC not implemented")
3448    if (Sigp%gwcomp==1) MSG_ERROR("AC with extrapolar technique not implemented")
3449  end if
3450 
3451  call gaps_free(gaps)
3452 
3453  ABI_FREE(val_indeces)
3454 
3455  DBG_EXIT('COLL')
3456 
3457 end subroutine setup_sigma

ABINIT/sigma_bksmask [ Functions ]

[ Top ] [ Functions ]

NAME

 sigma_bksmask

FUNCTION

  Compute tables for the distribution and the storage of the wavefunctions in the SIGMA code.

INPUTS

 Dtset<type(dataset_type)>=all input variables for this dataset
 Sigp<sigparams_t>=Parameters governing the self-energy calculation.
 Kmesh <kmesh_t>=Structure describing the k-point sampling.
 nprocs=Total number of MPI processors
 my_rank=Rank of this this processor.

OUTPUT

 my_spins(:)=Pointer to NULL in input. In output: list of spins treated by this node.
 bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state.
 keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space.
 ierr=Exit status.

PARENTS

      sigma

CHILDREN

SOURCE

3679 subroutine sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr)
3680 
3681 
3682 !This section has been created automatically by the script Abilint (TD).
3683 !Do not modify the following lines by hand.
3684 #undef ABI_FUNC
3685 #define ABI_FUNC 'sigma_bksmask'
3686 !End of the abilint section
3687 
3688  implicit none
3689 
3690 !Arguments ------------------------------------
3691 !scalars
3692  integer,intent(in) :: my_rank,nprocs
3693  integer,intent(out) :: ierr
3694  type(Dataset_type),intent(in) :: Dtset
3695  type(sigparams_t),intent(in) :: Sigp
3696  type(kmesh_t),intent(in) :: Kmesh
3697 !arrays
3698  integer,allocatable,intent(out) :: my_spins(:)
3699  logical,intent(out) :: bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
3700  logical,intent(out) :: keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
3701 
3702 !Local variables-------------------------------
3703 !scalars
3704  integer :: my_nspins,my_maxb,my_minb,isp,spin,band,nsppol,rank_spin
3705  logical :: store_ur
3706 !arrays
3707  integer :: tmp_spins(Sigp%nsppol),nprocs_spin(Sigp%nsppol)
3708 
3709 ! *************************************************************************
3710 
3711  ierr=0; nsppol=Sigp%nsppol
3712 
3713  ! List of spins for each node, number of processors per each spin
3714  ! and the MPI rank in the "spin" communicator.
3715  my_nspins=nsppol
3716  ABI_MALLOC(my_spins, (nsppol))
3717  my_spins= [(isp, isp=1,nsppol)]
3718  nprocs_spin = nprocs; rank_spin = my_rank
3719 
3720  if (nsppol==2 .and. nprocs>1) then
3721    ! Distribute spins (optimal distribution if nprocs is even)
3722    nprocs_spin(1) = nprocs/2
3723    nprocs_spin(2) = nprocs - nprocs/2
3724    my_nspins=1
3725    my_spins(1)=1
3726    if (my_rank+1>nprocs/2) then
3727      ! I will treate spin=2, compute shifted rank.
3728      my_spins(1)=2
3729      rank_spin = my_rank - nprocs/2
3730    end if
3731  end if
3732 
3733  store_ur = (MODULO(Dtset%gwmem,10)==1)
3734  keep_ur=.FALSE.; bks_mask=.FALSE.
3735 
3736  select case (Dtset%gwpara)
3737  case (1)
3738    ! Parallelization over transitions **without** memory distributions (Except for the spin).
3739    my_minb=1; my_maxb=Sigp%nbnds
3740    if (dtset%ucrpa>0) then
3741      bks_mask(my_minb:my_maxb,:,:)=.TRUE.
3742      if (store_ur) keep_ur(my_minb:my_maxb,:,:)=.TRUE.
3743    else
3744      do isp=1,my_nspins
3745        spin = my_spins(isp)
3746        bks_mask(my_minb:my_maxb,:,spin)=.TRUE.
3747        if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
3748      end do
3749    end if
3750 
3751  case (2)
3752    ! Distribute bands and spin (alternating planes of bands)
3753    do isp=1,my_nspins
3754      spin = my_spins(isp)
3755      do band=1,Sigp%nbnds
3756        if (xmpi_distrib_with_replicas(band,Sigp%nbnds,rank_spin,nprocs_spin(spin))) then
3757        !if (MODULO(band-1,nprocs_spin(spin))==rank_spin) then
3758          bks_mask(band,:,spin)=.TRUE.
3759          if (store_ur) keep_ur(band,:,spin)=.TRUE.
3760        end if
3761      end do
3762    end do
3763 
3764 #if 0
3765    ! Each node store the full set of occupied states to speed up Sigma_x.
3766    do isp=1,my_nspins
3767      spin = my_spins(isp)
3768      do ik_ibz=1,Kmesh%nibz
3769        ks_iv=ks_vbik(ik_ibz,spin) ! Valence index for this (k,s)
3770        bks_mask(1:ks_iv,:,spin)=.TRUE.
3771        if (store_ur) keep_ur(1:ks_iv,:,spin)=.TRUE.
3772      end do
3773    end do
3774 #endif
3775 
3776  case default
3777    MSG_WARNING("Wrong value for gwpara")
3778    ierr = 1
3779  end select
3780 
3781  ! Return my_spins with correct size.
3782  tmp_spins(1:my_nspins) = my_spins(1:my_nspins)
3783 
3784  ABI_FREE(my_spins)
3785  ABI_MALLOC(my_spins, (my_nspins))
3786  my_spins = tmp_spins(1:my_nspins)
3787 
3788 end subroutine sigma_bksmask

ABINIT/sigma_tables [ Functions ]

[ Top ] [ Functions ]

NAME

 sigma_tables

FUNCTION

  Build symmetry tables used to speedup self-consistent GW calculations.

INPUTS

 Kmesh <kmesh_t>=Structure describing the k-point sampling.
 [Bnd_sym(Kmesh%nibz,Sigp%nsppol)] <type(Bands_Symmetries)>

 SiDE EFFECTS
 Sigp<sigparams_t>=This routine initializes the tables:
   %Sigcij_tab
   %Sigxij_tab
  that are used to select the matrix elements of the self-energy that have to be calculated.

PARENTS

      setup_sigma,sigma

CHILDREN

SOURCE

3487 subroutine sigma_tables(Sigp,Kmesh,Bnd_sym)
3488 
3489 
3490 !This section has been created automatically by the script Abilint (TD).
3491 !Do not modify the following lines by hand.
3492 #undef ABI_FUNC
3493 #define ABI_FUNC 'sigma_tables'
3494 !End of the abilint section
3495 
3496  implicit none
3497 
3498 !Arguments ------------------------------------
3499 !scalars
3500  type(sigparams_t),intent(inout) :: Sigp
3501  type(kmesh_t),intent(in) :: Kmesh
3502 !arrays
3503  type(esymm_t),optional,intent(in) :: Bnd_sym(Kmesh%nibz,Sigp%nsppol)
3504 
3505 !Local variables-------------------------------
3506 !scalars
3507  integer :: gwcalctyp,spin,ikcalc,ik_ibz,bmin,bmax,bcol,brow
3508  integer :: ii,idx_x,idx_c,irr_idx1,irr_idx2
3509 !arrays
3510  integer,allocatable :: sigc_bidx(:),sigx_bidx(:)
3511  logical :: use_sym_at(Kmesh%nibz,Sigp%nsppol)
3512 
3513 ! *************************************************************************
3514 
3515  gwcalctyp=Sigp%gwcalctyp
3516 
3517  ! Recreate the Sig_ij tables taking advantage of the classification of the bands.
3518  if (allocated(Sigp%Sigxij_tab)) then
3519    call sigijtab_free(Sigp%Sigxij_tab)
3520    ABI_DT_FREE(Sigp%Sigxij_tab)
3521  end if
3522  if (allocated(Sigp%Sigcij_tab)) then
3523    call sigijtab_free(Sigp%Sigcij_tab)
3524    ABI_DT_FREE(Sigp%Sigcij_tab)
3525  end if
3526 
3527  ABI_DT_MALLOC(Sigp%Sigcij_tab,(Sigp%nkptgw,Sigp%nsppol))
3528  ABI_DT_MALLOC(Sigp%Sigxij_tab,(Sigp%nkptgw,Sigp%nsppol))
3529 
3530  use_sym_at=.FALSE.
3531  if (PRESENT(Bnd_sym)) then
3532    do spin=1,Sigp%nsppol
3533      do ikcalc=1,Sigp%nkptgw
3534       ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
3535       use_sym_at(ik_ibz,spin) = ( .not.esymm_failed(Bnd_sym(ik_ibz,spin)) )
3536      end do
3537    end do
3538  end if
3539 
3540  do spin=1,Sigp%nsppol
3541    do ikcalc=1,Sigp%nkptgw
3542      ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
3543 
3544      if (use_sym_at(ik_ibz,spin)) then
3545        if (gwcalctyp<20) then
3546          MSG_ERROR("You should not be here!")
3547        end if
3548 
3549        bmin=Sigp%minbnd(ikcalc,spin); bmax=Sigp%maxbnd(ikcalc,spin)
3550        ABI_DT_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col,(bmin:bmax))
3551        ABI_DT_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col,(bmin:bmax))
3552 
3553        do bcol=bmin,bmax
3554          ABI_MALLOC(sigc_bidx,(bmax-bmin+1))
3555          ABI_MALLOC(sigx_bidx,(bmax-bmin+1))
3556 
3557          if (Bnd_sym(ik_ibz,spin)%err_status/=0) then   ! Band classification failed.
3558            sigc_bidx = (/(ii,ii=bmin,bmax)/)
3559            idx_c = bmax-bmin+1
3560            sigx_bidx = (/(ii,ii=bmin,bcol)/) ! Hermitian
3561            idx_x = bcol-bmin+1
3562          else
3563            irr_idx2 = Bnd_sym(ik_ibz,spin)%b2irrep(bcol)
3564            idx_c = 0
3565            do brow=bmin,bmax
3566              irr_idx1 = Bnd_sym(ik_ibz,spin)%b2irrep(brow)
3567              if (sigma_is_herm(Sigp).and.bcol<brow) CYCLE  ! Only the upper triangle for HF, SEX, or COHSEX.
3568              if (irr_idx1 == irr_idx2) then ! same character, add this row to the list.
3569                idx_c = idx_c +1
3570                sigc_bidx(idx_c) = brow
3571              end if
3572            end do
3573            idx_x = 0
3574            do brow=bmin,bcol
3575              irr_idx1 = Bnd_sym(ik_ibz,spin)%b2irrep(brow)
3576              if (bcol<brow) CYCLE  ! Sig_x is always Hermitian.
3577              if (irr_idx1 == irr_idx2) then ! same character, add this row to the list.
3578                idx_x = idx_x +1
3579                sigx_bidx(idx_x) = brow
3580              end if
3581            end do
3582          end if
3583          !
3584          ! Table for Sigma_x matrix elements taking into account symmetries of the bands.
3585          ABI_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx,(idx_x))
3586 
3587          Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%size1= idx_x
3588          Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(:) = sigx_bidx(1:idx_x)
3589          !write(std_out,*)" Sigxij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol
3590          !write(std_out,*)" size: ",idx_x,(Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(ii),ii=1,idx_x)
3591          !
3592          ! Table for Sigma_c matrix elements taking into account symmetries of the bands.
3593          ABI_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx,(idx_c))
3594 
3595          Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%size1= idx_c
3596          Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(:) = sigc_bidx(1:idx_c)
3597          !write(std_out,*)" Sigcij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol
3598          !write(std_out,*)" size: ",idx_c,(Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(ii), ii=1,idx_c)
3599 
3600          ABI_FREE(sigx_bidx)
3601          ABI_FREE(sigc_bidx)
3602        end do ! bcol
3603 
3604      else  ! Symmetries cannot be used for this (k,s).
3605 
3606        bmin=Sigp%minbnd(ikcalc,spin); bmax=Sigp%maxbnd(ikcalc,spin)
3607        ABI_DT_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col,(bmin:bmax))
3608        ABI_DT_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col,(bmin:bmax))
3609 
3610        if (gwcalctyp<20) then  ! QP wavefunctions == KS, therefore only diagonal elements are calculated.
3611          do bcol=bmin,bmax
3612            ABI_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx,(1:1))
3613            Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%size1= 1
3614            Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(1) = bcol
3615            ABI_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx,(1:1))
3616            Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%size1= 1
3617            Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(1) = bcol
3618          end do
3619        else
3620          ! Use QP wavefunctions, Sigma_ij matrix is sparse but we have to classify the states in sigma.
3621          ! The only thing we can do here is filling the entire matrix taking advantage of Hermiticity (if any).
3622          do bcol=bmin,bmax
3623            ABI_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx,(bcol-bmin+1))
3624            Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%size1= bcol-bmin+1
3625            Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(:) = (/(ii,ii=bmin,bcol)/) ! Sigma_x is Hermitian.
3626            !write(std_out,*)"Sigxij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol,Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(:)
3627 
3628            ABI_MALLOC(sigc_bidx,(bmax-bmin+1))
3629            idx_c = 0
3630            do brow=bmin,bmax
3631              if (sigma_is_herm(Sigp).and.bcol<brow) CYCLE  ! Only the upper triangle of Sigc_ij is needed (SEX, COHSEX).
3632              idx_c = idx_c +1
3633              sigc_bidx(idx_c) = brow
3634            end do
3635            ABI_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx,(idx_c))
3636            Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%size1= idx_c
3637            Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(:) = sigc_bidx(1:idx_c)
3638            ABI_FREE(sigc_bidx)
3639            !write(std_out,*)"Sigcij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol,Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(:)
3640          end do
3641        end if
3642      end if
3643 
3644    end do !ikcalc
3645  end do !spin
3646 
3647 end subroutine sigma_tables

m_sigma_driver/sigma [ Functions ]

[ Top ] [ m_sigma_driver ] [ Functions ]

NAME

 sigma

FUNCTION

 Calculate the matrix elements of the self-energy operator.

INPUTS

 acell(3)=length scales of primitive translations (bohr)
 codvsn=code version
 Dtfil<type(datafiles_type)>=variables related to files
 Dtset<type(dataset_type)>=all input variables for this dataset
 Pawang<type(pawang_type)>=paw angular mesh and related data
 Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data
 Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data
 Psps<type(pseudopotential_type)>=variables related to pseudopotentials
   Before entering the first time in sigma, 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 screening, 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

OUTPUT

  Converged=.TRUE. if degw are within the user-specified tolerance.
  Output is written on the main abinit output file. Some results are stored in external files

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

      calc_sigc_me,calc_sigx_me,calc_ucrpa,calc_vhxc_me,chkpawovlp
      classify_bands,cohsex_me,crystal_free,denfgr,destroy_mpi_enreg
      ebands_copy,ebands_free,ebands_interpolate_kpath,ebands_report_gap
      ebands_update_occ,em1results_free,energies_init,esymm_free
      fftdatar_write,fourdp,get_gftt,getph,gsph_free,hdr_free
      init_distribfft_seq,initmpi_seq,kmesh_free,kxc_ada,kxc_driver
      littlegroup_free,littlegroup_init,melements_free,melements_print
      melements_zero,melflags_reset,metric,mkdump_er,mkrdim,nhatgrid
      paw_an_free,paw_an_init,paw_an_nullify,paw_check_symcprj,paw_dijhf
      paw_gencond,paw_ij_free,paw_ij_init,paw_ij_nullify,paw_ij_print
      paw_mkdijexc_core,paw_pwaves_lmn_free,paw_pwaves_lmn_init,paw_qpscgw
      pawcprj_alloc,pawcprj_free,pawdenpot,pawdij,pawfgr_destroy,pawfgr_init
      pawfgrtab_free,pawfgrtab_init,pawfgrtab_print,pawinit,pawmknhat,pawprt
      pawpuxinit,pawpwff_free,pawpwff_init,pawrhoij_alloc,pawrhoij_copy
      pawrhoij_free,pawtab_get_lsize,pawtab_print,ppm_free,ppm_init
      prep_calc_ucrpa,print_ngfft,prtrhomxmn,pspini,rdgw,rdqps,read_rhor
      setsym_ylm,setup_ppmodel,setup_sigma,setvtr,show_qp,sigma_bksmask
      sigma_free,sigma_init,sigma_tables,sigparams_free,solve_dyson,symdij
      symdij_all,test_charge,timab,updt_m_lda_to_qp,vcoul_free
      wfd_change_ngfft,wfd_copy,wfd_distribute_bands,wfd_free,wfd_get_cprj
      wfd_init,wfd_mkrho,wfd_print,wfd_read_wfk,wfd_reset_ur_cprj,wfd_rotate
      wfd_test_ortho,write_sigma_header,write_sigma_results,wrqps,wrtout
      xmpi_barrier,xmpi_bcast,xmpi_sum

SOURCE

 201 subroutine sigma(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim,converged)
 202 
 203 
 204 !This section has been created automatically by the script Abilint (TD).
 205 !Do not modify the following lines by hand.
 206 #undef ABI_FUNC
 207 #define ABI_FUNC 'sigma'
 208 !End of the abilint section
 209 
 210  implicit none
 211 
 212 !Arguments ------------------------------------
 213 !scalars
 214  logical,intent(out) :: converged
 215  character(len=6),intent(in) :: codvsn
 216  type(Datafiles_type),intent(in) :: Dtfil
 217  type(Dataset_type),intent(inout) :: Dtset
 218  type(Pawang_type),intent(inout) :: Pawang
 219  type(Pseudopotential_type),intent(inout) :: Psps
 220 !arrays
 221  real(dp),intent(in) :: acell(3),rprim(3,3)
 222  type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw)
 223  type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw)
 224 
 225 !Local variables-------------------------------
 226 !scalars
 227  integer,parameter :: tim_fourdp5=5,master=0,cplex1=1
 228  integer :: approx_type,b1gw,b2gw,choice,cplex,cplex_dij,band
 229  integer :: dim_kxcg,gwcalctyp,gnt_option,has_dijU,has_dijso,iab,bmin,bmax,irr_idx1,irr_idx2
 230  integer :: iat,ib,ib1,ib2,ic,id_required,ider,idir,ii,ik,ierr,ount
 231  integer :: ik_bz,ikcalc,ik_ibz,ikxc,ipert,npw_k,omp_ncpus
 232  integer :: isp,is_idx,istep,itypat,itypatcor,izero,jj,first_band,last_band
 233  integer :: ks_iv,lcor,lmn2_size_max,mband,my_nband
 234  integer :: mgfftf,mod10,moved_atm_inside,moved_rhor,n3xccc !,mgfft
 235  integer :: nbsc,ndij,ndim,nfftf,nfftf_tot,nkcalc,gwc_nfft,gwc_nfftot,gwx_nfft,gwx_nfftot
 236  integer :: ngrvdw,nhatgrdim,nkxc,nkxc1,nprocs,nscf,nspden_rhoij,nzlmopt,optene
 237  integer :: optcut,optgr0,optgr1,optgr2,option,option_test,option_dij,optrad,optrhoij,psp_gencond
 238  integer :: my_rank,rhoxsp_method,comm,use_aerhor,use_umklp,usexcnhat
 239  integer :: ioe0j,spin,io,jb,nomega_sigc
 240  integer :: temp_unt,ncid
 241  integer :: work_size,nstates_per_proc,my_nbks
 242  !integer :: jb_qp,ib_ks,ks_irr
 243  real(dp) :: compch_fft,compch_sph,r_s,rhoav,alpha,opt_ecut
 244  real(dp) :: drude_plsmf,my_plsmf,ecore,ecut_eff,ecutdg_eff,ehartree
 245  real(dp) :: ex_energy,gsqcutc_eff,gsqcutf_eff,gsqcut_shp,norm,oldefermi
 246  real(dp) :: ucvol,vxcavg,vxcavg_qp
 247  real(dp) :: gwc_gsq,gwx_gsq,gw_gsq
 248  real(dp):: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem
 249  complex(dpc) :: max_degw,cdummy
 250  logical :: use_paw_aeur,dbg_mode,pole_screening,call_pawinit
 251  character(len=500) :: msg
 252  character(len=fnlen) :: wfk_fname,pawden_fname
 253  type(kmesh_t) :: Kmesh,Qmesh
 254  type(ebands_t) :: KS_BSt,QP_BSt
 255  type(vcoul_t) :: Vcp
 256  type(crystal_t) :: Cryst
 257  type(Energies_type) :: KS_energies,QP_energies
 258  type(Epsilonm1_results) :: Er
 259  type(gsphere_t) :: Gsph_Max,Gsph_x,Gsph_c
 260  type(hdr_type) :: Hdr_wfk,Hdr_sigma,hdr_rhor
 261  type(melflags_t) :: KS_mflags,QP_mflags
 262  type(melements_t) :: KS_me,QP_me
 263  type(MPI_type) :: MPI_enreg_seq
 264  type(paw_dmft_type) :: Paw_dmft
 265  type(pawfgr_type) :: Pawfgr
 266  type(ppmodel_t) :: PPm
 267  type(sigparams_t) :: Sigp
 268  type(sigma_t) :: Sr
 269  type(wfd_t),target :: Wfd,Wfdf
 270  type(wvl_data) :: Wvl
 271 !arrays
 272  integer :: gwc_ngfft(18),ngfftc(18),ngfftf(18),gwx_ngfft(18)
 273  integer,allocatable :: nq_spl(:),nlmn_atm(:),my_spins(:)
 274  integer,allocatable :: tmp_gfft(:,:),ks_vbik(:,:),nband(:,:),l_size_atm(:),qp_vbik(:,:)
 275  integer,allocatable :: tmp_kstab(:,:,:),ks_irreptab(:,:,:),qp_irreptab(:,:,:),my_band_list(:)
 276  real(dp),parameter ::  k0(3)=zero
 277  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),strsxc(6),tsec(2)
 278  real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:)
 279  real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:)
 280  real(dp),allocatable :: ks_rhor(:,:),ks_vhartr(:),ks_vtrial(:,:),ks_vxc(:,:)
 281  real(dp),allocatable :: ks_taug(:,:),ks_taur(:,:)
 282  real(dp),allocatable :: kxc(:,:),qp_kxc(:,:),ph1d(:,:),ph1df(:,:)
 283  real(dp),allocatable :: prev_rhor(:,:),prev_taur(:,:),qp_nhat(:,:)
 284  real(dp),allocatable :: qp_nhatgr(:,:,:),qp_rhog(:,:),qp_rhor_paw(:,:)
 285  real(dp),allocatable :: qp_rhor_n_one(:,:),qp_rhor_nt_one(:,:)
 286  real(dp),allocatable :: qp_rhor(:,:),qp_vhartr(:),qp_vtrial(:,:),qp_vxc(:,:)
 287  real(dp),allocatable :: qp_taur(:,:),qp_taug(:,:),igwene(:,:,:)
 288  real(dp),allocatable :: vpsp(:),xccc3d(:),dijexc_core(:,:,:),dij_hf(:,:,:)
 289  !real(dp),allocatable :: osoc_bks(:, :, :)
 290  real(dp),allocatable :: ks_aepaw_rhor(:,:) !,ks_n_one_rhor(:,:),ks_nt_one_rhor(:,:)
 291  complex(dpc) :: ovlp(2)
 292  complex(dpc),allocatable :: ctmp(:,:),hbare(:,:,:,:)
 293  complex(dpc),target,allocatable :: sigcme(:,:,:,:,:)
 294  complex(dpc),allocatable :: hlda(:,:,:,:),htmp(:,:,:,:),uks2qp(:,:)
 295  complex(gwpc),allocatable :: kxcg(:,:),fxc_ADA(:,:,:)
 296  complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:)
 297 !complex(dpc),pointer :: sigcme_p(:,:,:,:)
 298  complex(dpc),allocatable :: sigcme_k(:,:,:,:)
 299  complex(dpc), allocatable :: rhot1_q_m(:,:,:,:,:,:,:)
 300  complex(dpc), allocatable :: M1_q_m(:,:,:,:,:,:,:)
 301  logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:),bmask(:)
 302  type(esymm_t),target,allocatable :: KS_sym(:,:)
 303  type(esymm_t),pointer :: QP_sym(:,:)
 304  type(pawcprj_type),allocatable :: Cp1(:,:),Cp2(:,:)
 305  type(littlegroup_t),allocatable :: Ltg_k(:)
 306  type(Paw_an_type),allocatable :: KS_paw_an(:),QP_paw_an(:)
 307  type(Paw_ij_type),allocatable :: KS_paw_ij(:),QP_paw_ij(:)
 308  type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
 309  type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:),QP_pawrhoij(:),prev_Pawrhoij(:),tmp_pawrhoij(:)
 310  type(pawpwff_t),allocatable :: Paw_pwff(:)
 311  type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
 312 
 313 !************************************************************************
 314 
 315  DBG_ENTER('COLL')
 316 
 317  call timab(401,1,tsec) ! sigma(Total)
 318  call timab(402,1,tsec) ! sigma(Init1)
 319 
 320  write(msg,'(7a)')&
 321 & ' SIGMA: Calculation of the GW corrections ',ch10,ch10,&
 322 & ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,&
 323 & ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.'
 324  call wrtout(std_out,msg,'COLL')
 325  call wrtout(ab_out,msg,'COLL')
 326 
 327  if(dtset%ucrpa>0) then
 328    write(msg,'(6a)')ch10,&
 329 &   ' cRPA Calculation: Calculation of the screened Coulomb interaction (ucrpa/=0) ',ch10
 330    call wrtout(ab_out, msg,'COLL')
 331    call wrtout(std_out,msg,'COLL')
 332  end if
 333 
 334 #if defined HAVE_GW_DPC
 335  if (gwpc/=8) then
 336    write(msg,'(6a)')ch10,&
 337 &   'Number of bytes for double precision complex /=8 ',ch10,&
 338 &   'Cannot continue due to kind mismatch in BLAS library ',ch10,&
 339 &   'Some BLAS interfaces are not generated by abilint '
 340    MSG_ERROR(msg)
 341  end if
 342  write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10
 343 #else
 344  write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10
 345 #endif
 346  call wrtout(std_out,msg,'COLL')
 347  call wrtout(ab_out,msg,'COLL')
 348 
 349  gwcalctyp=Dtset%gwcalctyp
 350  mod10 =MOD(Dtset%gwcalctyp,10)
 351 
 352  ! Perform some additional checks for hybrid functional calculations
 353  if(mod10==5) then
 354    if (Dtset%ixc_sigma<0 .and. .not.libxc_functionals_check()) then
 355      msg='Hybrid functional calculations require the compilation with LIBXC library'
 356      MSG_ERROR(msg)
 357    end if
 358 !  XG 20171116 : I do not agree with this condition, as one might like to do a one-shot hybrid functional calculation
 359 !  on top of a LDA/GGA calculation ... give the power (and risks) to the user !
 360 !  if(gwcalctyp<10) then
 361 !    msg='gwcalctyp should enforce update of the energies and/or wavefunctions when performing hybrid functional calculation'
 362 !    MSG_ERROR(msg)
 363 !  end if
 364    if(Dtset%usepaw==1) then
 365      msg='PAW version of hybrid functional calculations is not implemented'
 366      MSG_ERROR(msg)
 367    end if
 368  end if
 369 
 370  !=== Initialize MPI variables, and parallelization level ===
 371  !* gwpara 0--> sequential run, 1--> parallelism over k-points, 2--> parallelism over bands.
 372  !* In case of gwpara==1 memory is not parallelized.
 373  !* If gwpara==2, bands are divided among processors but each proc has all the states where GW corrections are required.
 374  comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
 375  converged = .FALSE.
 376 
 377  if (my_rank == master) then
 378    wfk_fname = dtfil%fnamewffk
 379    if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then
 380      MSG_ERROR(msg)
 381    end if
 382  end if
 383  call xmpi_bcast(wfk_fname, master, comm, ierr)
 384 
 385  ! === Some variables need to be initialized/nullify at start ===
 386  call energies_init(KS_energies)
 387  usexcnhat=0
 388  call mkrdim(acell,rprim,rprimd)
 389  call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
 390  !
 391  ! === Define FFT grid(s) sizes ===
 392  ! * Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the
 393  ! (usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90.
 394  ! See also NOTES in the comments at the beginning of this file.
 395  ! NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn.
 396 
 397  call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,&
 398 & gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0)
 399 
 400  ! Fake MPI_type for the sequential part.
 401  call initmpi_seq(MPI_enreg_seq)
 402  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all')
 403  call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all')
 404 
 405  call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials')
 406  nfftf_tot=PRODUCT(ngfftf(1:3))
 407 
 408  ! Open and read pseudopotential files ===
 409  call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm)
 410 
 411  call timab(402,2,tsec) ! Init1
 412  !
 413  ! ===============================================
 414  ! ==== Initialize Sigp, Er and basic objects ====
 415  ! ===============================================
 416  ! * Sigp is completetly initialized here.
 417  ! * Er is only initialized with dimensions, (SCR|SUSC) file is read in mkdump_Er
 418  call timab(403,1,tsec) ! setup_sigma
 419 
 420  call setup_sigma(codvsn,wfk_fname,acell,rprim,ngfftf,Dtset,Dtfil,Psps,Pawtab,&
 421 & gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_sigma,Cryst,Kmesh,Qmesh,KS_BSt,Gsph_Max,Gsph_x,Gsph_c,Vcp,Er,Sigp,comm)
 422 
 423  call timab(403,2,tsec) ! setup_sigma
 424  call timab(402,1,tsec) ! Init1
 425 
 426 !XG090617 Please, do not remove this write, unless you have checked
 427 !that the code executes correctly on max+g95 (especially, Tv5#70).
 428 !It is one more a silly write, perhaps needed because the compiler does not treat correctly non-nullified pointers.
 429  if (sigma_needs_w(Sigp) .and. my_rank==master) then
 430    write(std_out,*)' screening after setup_sigma : Er%Hscr%headform=',Er%Hscr%headform
 431  end if
 432 !END XG090617
 433 
 434  pole_screening = .FALSE.
 435  if (Er%fform==2002) then
 436    pole_screening = .TRUE.
 437    MSG_WARNING(' EXPERIMENTAL - Using a pole-fit screening!')
 438  end if
 439 
 440  call print_ngfft(gwc_ngfft,header='FFT mesh for oscillator strengths used for Sigma_c')
 441  call print_ngfft(gwx_ngfft,header='FFT mesh for oscillator strengths used for Sigma_x')
 442 
 443  b1gw=Sigp%minbdgw
 444  b2gw=Sigp%maxbdgw
 445 
 446  gwc_nfftot=PRODUCT(gwc_ngfft(1:3))
 447  gwc_nfft  =gwc_nfftot  !no FFT //
 448 
 449  gwx_nfftot=PRODUCT(gwx_ngfft(1:3))
 450  gwx_nfft  =gwx_nfftot  !no FFT //
 451 !
 452 !TRYING TO RECREATE AN "ABINIT ENVIRONMENT"
 453  KS_energies%e_corepsp=ecore/Cryst%ucvol
 454 
 455 !=== Calculate KS occupation numbers and ks_vbk(nkibz,nsppol) ====
 456 !* ks_vbk gives the (valence|last Fermi band) index for each k and spin.
 457 !* spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals
 458  ABI_MALLOC(ks_vbik,(KS_BSt%nkpt,KS_BSt%nsppol))
 459  ABI_MALLOC(qp_vbik,(KS_BSt%nkpt,KS_BSt%nsppol))
 460 
 461  !call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=0)
 462  ks_vbik(:,:) = get_valence_idx(KS_BSt)
 463 
 464  ! ============================
 465  ! ==== PAW initialization ====
 466  ! ============================
 467  if (Dtset%usepaw==1) then
 468    call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,Cryst%xred)
 469 
 470    cplex_dij=Dtset%nspinor; cplex=1; ndij=1
 471 
 472    ABI_DT_MALLOC(KS_Pawrhoij,(Cryst%natom))
 473    nspden_rhoij=pawrhoij_get_nspden(Dtset%nspden,Dtset%nspinor,Dtset%pawspnorb)
 474    call pawrhoij_alloc(KS_Pawrhoij,Dtset%pawcpxocc,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab)
 475 
 476    ! Initialize values for several basic arrays
 477    gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
 478 
 479    ! Test if we have to call pawinit
 480    call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
 481 
 482    if (psp_gencond==1.or. call_pawinit) then
 483      call timab(553,1,tsec)
 484      gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
 485      call pawinit(gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,&
 486 &     Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,&
 487 &     Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%xclevel,Dtset%usepotzero)
 488      call timab(553,2,tsec)
 489 
 490      ! Update internal values
 491      call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
 492 
 493    else
 494      if (Pawtab(1)%has_kij  ==1) Pawtab(1:Cryst%ntypat)%has_kij  =2
 495      if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2
 496    end if
 497    Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore)
 498 
 499    ! Initialize optional flags in Pawtab to zero
 500    ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed)
 501    Pawtab(:)%has_nabla = 0
 502    Pawtab(:)%usepawu   = 0
 503    Pawtab(:)%useexexch = 0
 504    Pawtab(:)%exchmix   =zero
 505 
 506    call setsym_ylm(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,Cryst%rprimd,Cryst%symrec,Pawang%zarot)
 507 
 508    ! Initialize and compute data for LDA+U
 509    Paw_dmft%use_dmft=Dtset%usedmft
 510    if (Dtset%usepawu>0.or.Dtset%useexexch>0) then
 511      call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,&
 512 &     Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,Cryst%ntypat,Pawang,Dtset%pawprtvol,&
 513 &     Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa)
 514    end if
 515 
 516    if (my_rank == master) call pawtab_print(Pawtab)
 517 
 518    ! Get Pawrhoij from the header of the WFK file.
 519    call pawrhoij_copy(Hdr_wfk%pawrhoij,KS_Pawrhoij)
 520 
 521    ! Re-symmetrize symrhoij.
 522    ! FIXME this call leads to a SIGFAULT, likely some pointer is not initialized correctly
 523    choice=1; optrhoij=1; ipert=0; idir=0
 524 !  call symrhoij(KS_Pawrhoij,KS_Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert,&
 525 !  &             Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,&
 526 !  &             Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
 527 
 528    !  Evaluate form factor of radial part of phi.phj-tphi.tphj.
 529    rhoxsp_method=1 ! Arnaud-Alouani
 530 !  rhoxsp_method=2 ! Shiskin-Kresse
 531    if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc
 532 
 533 !  The q-grid must contain the FFT mesh used for sigma_c and the G-sphere for the exchange part.
 534 !  We use the FFT mesh for sigma_c since COHSEX and the extrapolar method require oscillator
 535 !  strengths on the FFT mesh.
 536    ABI_MALLOC(tmp_gfft,(3,gwc_nfftot))
 537    call get_gftt(gwc_ngfft,k0,gmet,gwc_gsq,tmp_gfft)
 538    ABI_FREE(tmp_gfft)
 539 
 540    gwx_gsq = Dtset%ecutsigx/(two*pi**2)
 541 !  allocate(tmp_gfft(3,gwx_nfftot)); q0=zero
 542 !  call get_gftt(gwx_ngfft,q0,gmet,gwx_gsq,tmp_gfft)
 543 !  deallocate(tmp_gfft)
 544    gw_gsq = MAX(gwx_gsq,gwc_gsq)
 545 
 546 !  * Set up q-grid, make qmax 20% larger than largest expected.
 547    ABI_MALLOC(nq_spl,(Psps%ntypat))
 548    ABI_MALLOC(qmax,(Psps%ntypat))
 549    qmax = SQRT(gw_gsq)*1.2d0 ! qmax=Psps%qgrid_ff(Psps%mqgrid_ff)
 550    nq_spl = Psps%mqgrid_ff
 551 !  write(std_out,*)"using nq_spl",nq_spl,"qmax=",qmax
 552    ABI_DT_MALLOC(Paw_pwff,(Psps%ntypat))
 553 
 554    call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps)
 555 
 556    ABI_FREE(nq_spl)
 557    ABI_FREE(qmax)
 558    !
 559    ! Variables/arrays related to the fine FFT grid ===
 560    ABI_MALLOC(ks_nhat,(nfftf,Dtset%nspden))
 561    ks_nhat=zero
 562    ABI_DT_MALLOC(Pawfgrtab,(Cryst%natom))
 563    call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
 564    cplex=1
 565    call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat)
 566    ABI_FREE(l_size_atm)
 567    compch_fft=greatest_real
 568    usexcnhat=MAXVAL(Pawtab(:)%usexcnhat)
 569    !  * 0 if Vloc in atomic data is Vbare    (Blochl s formulation)
 570    !  * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation)
 571    write(msg,'(a,i2)')' sigma : using usexcnhat = ',usexcnhat
 572    call wrtout(std_out,msg,'COLL')
 573    !
 574    ! Identify parts of the rectangular grid where the density has to be calculated ===
 575    optcut=0; optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm
 576    if (Dtset%pawcross==1) optrad=1
 577    if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm
 578 
 579    call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,&
 580 &   optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
 581 
 582    call pawfgrtab_print(Pawfgrtab,Cryst%natom,unit=std_out,prtvol=Dtset%pawprtvol)
 583  else
 584    ABI_DT_MALLOC(Paw_pwff,(0))
 585    ABI_DT_MALLOC(Pawfgrtab,(0))
 586  end if !End of PAW Initialization
 587 
 588  ! Consistency check and additional stuff done only for GW with PAW.
 589  ABI_DT_MALLOC(Paw_onsite,(0))
 590 
 591  if (Dtset%usepaw==1) then
 592    if (Dtset%ecutwfn < Dtset%ecut) then
 593      write(msg,"(5a)")&
 594 &     "WARNING - ",ch10,&
 595 &     "  It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,&
 596 &     "  an excessive truncation of the planewave basis set can lead to unphysical results."
 597      MSG_WARNING(msg)
 598    end if
 599 
 600    ABI_CHECK(Dtset%useexexch==0,"LEXX not yet implemented in GW")
 601    ABI_CHECK(Paw_dmft%use_dmft==0,"DMFT + GW not available")
 602 
 603    ! Optionally read core orbitals from file and calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for the HF decoupling.
 604    if (Sigp%use_sigxcore==1) then
 605      lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
 606      ABI_MALLOC(dijexc_core,(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat))
 607 
 608      call paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,Dtset%prtvol,Psps%filpsp)
 609    end if ! HF decoupling
 610 
 611    if (Dtset%pawcross==1) then
 612      if (allocated(Paw_onsite) ) then
 613        ABI_DT_FREE(Paw_onsite)
 614      end if
 615      ABI_DT_MALLOC(Paw_onsite,(Cryst%natom))
 616      call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,&
 617 &     Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab)
 618    end if
 619  end if
 620 
 621  ! Allocate these arrays anyway, since they are passed to subroutines.
 622  if (.not.allocated(ks_nhat)) then
 623    ABI_MALLOC(ks_nhat,(nfftf,0))
 624  end if
 625  if (.not.allocated(dijexc_core)) then
 626    ABI_MALLOC(dijexc_core,(1,1,0))
 627  end if
 628 
 629  ! ==================================================
 630  ! ==== Read KS band structure from the KSS file ====
 631  ! ==================================================
 632  !
 633  ! Initialize Wfd, allocate wavefunctions and precalculate tables to do the FFT using the coarse gwc_ngfft.
 634  mband=Sigp%nbnds
 635  ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Sigp%nsppol))
 636  ABI_MALLOC(keep_ur ,(mband,Kmesh%nibz,Sigp%nsppol))
 637  keep_ur=.FALSE.; bks_mask=.FALSE.
 638 
 639  ABI_MALLOC(nband,(Kmesh%nibz,Sigp%nsppol))
 640  nband=mband
 641 
 642  ! autoparal section
 643  if (dtset%max_ncpus/=0) then
 644    ount =ab_out
 645     ! Temporary table needed to estimate memory
 646    ABI_MALLOC(nlmn_atm,(Cryst%natom))
 647    if (Dtset%usepaw==1) then
 648      do iat=1,Cryst%natom
 649        nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size
 650      end do
 651    end if
 652 
 653    write(ount,'(a)')"--- !Autoparal"
 654    write(ount,"(a)")'#Autoparal section for Sigma runs.'
 655    write(ount,"(a)")   "info:"
 656    write(ount,"(a,i0)")"    autoparal: ",dtset%autoparal
 657    write(ount,"(a,i0)")"    max_ncpus: ",dtset%max_ncpus
 658    write(ount,"(a,i0)")"    gwpara: ",dtset%gwpara
 659    write(ount,"(a,i0)")"    nkpt: ",dtset%nkpt
 660    write(ount,"(a,i0)")"    nsppol: ",dtset%nsppol
 661    write(ount,"(a,i0)")"    nspinor: ",dtset%nspinor
 662    write(ount,"(a,i0)")"    nbnds: ",Sigp%nbnds
 663 
 664    work_size = mband * Kmesh%nibz * Sigp%nsppol
 665 
 666    ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI.
 667    nonscal_mem = (two*gwpc*Er%npwe**2*Er%nomega*(Er%mqmem+1)*b2Mb) * 1.1_dp
 668 
 669    ! List of configurations.
 670    ! Assuming an OpenMP implementation with perfect speedup!
 671    write(ount,"(a)")"configurations:"
 672    do ii=1,dtset%max_ncpus
 673      nstates_per_proc = 0
 674      eff = HUGE(one)
 675      max_wfsmem_mb = zero
 676 
 677      do my_rank=0,ii-1
 678        call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,ii,my_spins,bks_mask,keep_ur,ierr)
 679        ABI_FREE(my_spins)
 680        if (ierr /= 0) exit
 681        my_nbks = COUNT(bks_mask)
 682        nstates_per_proc = MAX(nstates_per_proc, my_nbks)
 683        eff = MIN(eff, (one * work_size) / (ii * nstates_per_proc))
 684 
 685        ! Memory needed for Fourier components ug.
 686        ug_mem = two*gwpc*Dtset%nspinor*Sigp%npwwfn*my_nbks*b2Mb
 687        ! Memory needed for real space ur (use gwc_nfft, instead of gwx_nfft)
 688        ur_mem = two*gwpc*Dtset%nspinor*gwc_nfft*COUNT(keep_ur)*b2Mb
 689        ! Memory needed for PAW projections Cprj
 690        cprj_mem = zero
 691        if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb
 692        max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem)
 693      end do
 694      if (ierr /= 0) cycle
 695 
 696      ! Add the non-scalable part and increase by 10% to account for other datastructures.
 697      mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp
 698      do omp_ncpus=1,xomp_get_max_threads()
 699        write(ount,"(a,i0)")"    - tot_ncpus: ",ii * omp_ncpus
 700        write(ount,"(a,i0)")"      mpi_ncpus: ",ii
 701        write(ount,"(a,i0)")"      omp_ncpus: ",omp_ncpus
 702        write(ount,"(a,f12.9)")"      efficiency: ",eff
 703        write(ount,"(a,f12.2)")"      mem_per_cpu: ",mempercpu_mb
 704      end do
 705    end do
 706    write(ount,'(a)')"..."
 707 
 708    ABI_FREE(nlmn_atm)
 709    MSG_ERROR_NODUMP("aborting now")
 710 
 711  else
 712    call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr)
 713    ABI_CHECK(ierr==0, "Error in sigma_bksmask")
 714  end if
 715 
 716  ! Then each node owns the wavefunctions where GW corrections are required.
 717  do isp=1,SIZE(my_spins)
 718    spin = my_spins(isp)
 719    do ikcalc=1,Sigp%nkptgw
 720      ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW
 721      ii=Sigp%minbnd(ikcalc,spin); jj=Sigp%maxbnd(ikcalc,spin)
 722      bks_mask(ii:jj,ik_ibz,spin) = .TRUE.
 723      if (MODULO(Dtset%gwmem,10)==1) keep_ur(ii:jj,ik_ibz,spin)=.TRUE.
 724    end do
 725  end do
 726 
 727  ABI_FREE(my_spins)
 728  opt_ecut=Dtset%ecutwfn
 729 !opt_ecut=zero
 730 
 731  call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,Sigp%npwwfn,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,&
 732 & Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,&
 733 & Gsph_Max%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut)
 734 
 735  if (Dtset%pawcross==1) then
 736    call wfd_init(Wfdf,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,Sigp%npwwfn,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,&
 737 &   Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,&
 738 &   Gsph_Max%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut)
 739  end if
 740 
 741  ABI_FREE(bks_mask)
 742  ABI_FREE(nband)
 743  ABI_FREE(keep_ur)
 744 
 745  call timab(402,2,tsec) ! sigma(Init1)
 746  call timab(404,1,tsec) ! rdkss
 747 
 748  call wfd_read_wfk(Wfd,wfk_fname,iomode_from_fname(wfk_fname))
 749 
 750  if (Dtset%pawcross==1) then
 751    call wfd_copy(Wfd,Wfdf)
 752    call wfd_change_ngfft(Wfdf,Cryst,Psps,ngfftf)
 753  end if
 754 
 755  ! This test has been disabled (too expensive!)
 756  if (.False.) call wfd_test_ortho(Wfd,Cryst,Pawtab,unit=ab_out,mode_paral="COLL")
 757 
 758  call timab(404,2,tsec) ! rdkss
 759  call timab(405,1,tsec) ! Init2
 760 
 761  ! ==============================================================
 762  ! ==== Find little group of the k-points for GW corrections ====
 763  ! ==============================================================
 764  ! * The little group is calculated only if sys_sigma.
 765  ! * If use_umklp==1 then also symmetries requiring an umklapp to preserve k_gw are included.
 766  !
 767  ABI_DT_MALLOC(Ltg_k,(Sigp%nkptgw))
 768  use_umklp=1
 769  do ikcalc=1,Sigp%nkptgw
 770    if (Sigp%symsigma/=0) then
 771      call littlegroup_init(Sigp%kptgw(:,ikcalc),Qmesh,Cryst,use_umklp,Ltg_k(ikcalc),0)
 772    end if
 773  end do
 774 !
 775 !=== Compute structure factor phases and large sphere cut-off ===
 776 !WARNING cannot use Dtset%mgfft, this has to be checked better
 777 !mgfft=MAXVAL(ngfftc(:))
 778 !allocate(ph1d(2,3*(2*mgfft+1)*Cryst%natom),ph1df(2,3*(2*mgfftf+1)*Cryst%natom))
 779 !write(std_out,*)' CHECK ',Dtset%mgfftdg,mgfftf
 780 !if (Dtset%mgfftdg/=mgfftf) then
 781 !write(std_out,*)"WARNING Dtset%mgfftf /= mgfftf"
 782 !write(std_out,*)'HACKING Dtset%mgfftf'
 783 !Dtset%mgfftdg=mgfftf
 784 !end if
 785  ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom))
 786  ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom))
 787 
 788  call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred)
 789 
 790  if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then
 791    call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred)
 792  else
 793    ph1df(:,:)=ph1d(:,:)
 794  end if
 795 
 796  !===================================================================================
 797  !==== Classify the GW wavefunctions according to the irreducible representation ====
 798  !===================================================================================
 799  !* Warning still under development.
 800  !* Only for SCGW.
 801  !bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
 802 
 803  ABI_DT_MALLOC(KS_sym,(Wfd%nkibz,Wfd%nsppol))
 804 
 805  if (Sigp%symsigma==1.and.gwcalctyp>=20) then
 806    ! call check_zarot(Gsph_c%ng,Cryst,gwc_ngfft,Gsph_c%gvec,Psps,Pawang,Gsph_c%rottb,Gsph_c%rottbm1)
 807    use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric
 808    do spin=1,Wfd%nsppol
 809      do ikcalc=1,Sigp%nkptgw
 810        ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
 811        first_band = Sigp%minbnd(ikcalc,spin)
 812        last_band  = Sigp%maxbnd(ikcalc,spin)
 813 !      call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,KS_BSt,Pawtab,Pawrad,Pawang,Psps,&
 814        call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,KS_BSt,Pawtab,Pawrad,Pawang,Psps,&
 815 &       Dtset%tolsym,KS_sym(ik_ibz,spin))
 816      end do
 817    end do
 818    ! Recreate the Sig_ij tables taking advantage of the classification of the bands.
 819    call sigma_tables(Sigp,Kmesh,KS_sym)
 820  end if
 821 
 822  call timab(405,2,tsec) ! Init2
 823  call timab(406,1,tsec) ! make_vhxc
 824 
 825  !===========================
 826  !=== COMPUTE THE DENSITY ===
 827  !===========================
 828  ! * Evaluate the planewave part (complete charge in case of NC pseudos).
 829  ABI_MALLOC(ks_rhor,(nfftf,Dtset%nspden))
 830  ABI_MALLOC(ks_taur,(nfftf,Dtset%nspden*Dtset%usekden))
 831 
 832  call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,KS_BSt,ngfftf,nfftf,ks_rhor)
 833 
 834  if (Dtset%usekden==1) then
 835    call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,KS_BSt,ngfftf,nfftf,ks_taur,optcalc=1)
 836  end if
 837 
 838  !========================================
 839  !==== Additional computation for PAW ====
 840  !========================================
 841  nhatgrdim=0
 842  if (Dtset%usepaw==1) then
 843    ! Calculate the compensation charge nhat.
 844    if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc
 845    cplex=1; ider=2*nhatgrdim; izero=0
 846    if (nhatgrdim>0)  then
 847      ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim))
 848    end if
 849    if (nhatgrdim==0)  then
 850      ABI_MALLOC(ks_nhatgr,(0,0,0))
 851    end if
 852 
 853    call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,Cryst%gprimd,&
 854 &   Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
 855 &   Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,k0,Cryst%rprimd,&
 856 &   Cryst%ucvol,dtset%usewvl,Cryst%xred)
 857 
 858    !if (nhatgrdim==0)  then
 859    !  ABI_FREE(ks_nhatgr)
 860    !end if
 861 
 862    ! === Evaluate onsite energies, potentials, densities ===
 863    !  * Initialize variables/arrays related to the PAW spheres.
 864    !  * Initialize also lmselect (index of non-zero LM-moments of densities).
 865    ABI_DT_MALLOC(KS_paw_ij,(Cryst%natom))
 866 !  cplex=1;cplex_dij=Dtset%nspinor
 867    has_dijso=Dtset%pawspnorb; has_dijU=Dtset%usepawu
 868    call paw_ij_nullify(KS_paw_ij)
 869    call paw_ij_init(KS_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,&
 870 &   Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
 871 &   has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,&
 872 &   has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1)
 873 
 874    nkxc1=0
 875    ABI_DT_MALLOC(KS_paw_an,(Cryst%natom))
 876    call paw_an_nullify(KS_paw_an)
 877    call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
 878 &   cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1)
 879 !
 880 !  Calculate onsite vxc with and without core charge.
 881    nzlmopt=-1; option=0; compch_sph=greatest_real
 882    call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert,&
 883 &   Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
 884 &   Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,&
 885 &   Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,&
 886 &   Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp)
 887 
 888 ! TO BE REMOVED ASAP !!!
 889    _IBM6("Silly write for IBM6")
 890  else
 891    ABI_MALLOC(ks_nhatgr,(0,0,0))
 892    ABI_DT_MALLOC(KS_paw_ij,(0))
 893    ABI_DT_MALLOC(KS_paw_an,(0))
 894  end if !PAW
 895 
 896  !if (.not.allocated(ks_nhatgr))  then
 897  !  ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,0))
 898  !end if
 899 
 900  call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,&
 901 & Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
 902 !
 903 !For PAW, add the compensation charge on the FFT mesh, then get rho(G).
 904  if (Dtset%usepaw==1) ks_rhor=ks_rhor+ks_nhat
 905 
 906  call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_rhor,ucvol=ucvol)
 907  if(Dtset%usekden==1)then
 908    call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=ucvol)
 909  end if
 910 
 911  ABI_MALLOC(ks_rhog,(2,nfftf))
 912  ABI_MALLOC(ks_taug,(2,nfftf*Dtset%usekden))
 913  call fourdp(1,ks_rhog,ks_rhor(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp5)
 914  if(Dtset%usekden==1)then
 915    call fourdp(1,ks_taug,ks_taur(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp5)
 916  end if
 917 
 918  !The following steps have been gathered in the setvtr routine:
 919  !- get Ewald energy and Ewald forces
 920  !- compute local ionic pseudopotential vpsp
 921  !- eventually compute 3D core electron density xccc3d
 922  !- eventually compute vxc and vhartr
 923  !- set up ks_vtrial
 924  !
 925  !*******************************************************************
 926  !**** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
 927  !*******************************************************************
 928 
 929  ngrvdw=0
 930  ABI_MALLOC(grvdw,(3,ngrvdw))
 931  ABI_MALLOC(grchempottn,(3,Cryst%natom))
 932  ABI_MALLOC(grewtn,(3,Cryst%natom))
 933  nkxc=0
 934  if (Dtset%nspden==1) nkxc=2
 935  if (Dtset%nspden>=2) nkxc=3 ! check GGA and spinor, quite a messy part!!!
 936 !In case of MGGA, fxc and kxc are not available and we dont need them
 937 !for the screening part (for now ...)
 938  if (Dtset%ixc<0.and.libxc_functionals_ismgga()) nkxc=0
 939  if (nkxc/=0)  then
 940    ABI_MALLOC(kxc,(nfftf,nkxc))
 941  end if
 942 
 943  n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf
 944  ABI_MALLOC(xccc3d,(n3xccc))
 945  ABI_MALLOC(ks_vhartr,(nfftf))
 946  ABI_MALLOC(ks_vtrial,(nfftf,Dtset%nspden))
 947  ABI_MALLOC(vpsp,(nfftf))
 948  ABI_MALLOC(ks_vxc,(nfftf,Dtset%nspden))
 949 
 950  optene=4; moved_atm_inside=0; moved_rhor=0; istep=1
 951 
 952  call setvtr(Cryst%atindx1,Dtset,KS_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
 953 & istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
 954 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
 955 & optene,pawrad,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,Cryst%rprimd,strsxc,&
 956 & Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,Wvl,xccc3d,Cryst%xred,taug=ks_taug,taur=ks_taur)
 957 
 958 !============================
 959 !==== Compute KS PAW Dij ====
 960 !============================
 961  if (Dtset%usepaw==1) then
 962    !TO BE REMOVED
 963    _IBM6("Another silly write for IBM6")
 964    call timab(561,1,tsec)
 965 
 966    !  Calculate the unsymmetrized Dij.
 967    ipert=0; idir=0
 968    call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert,&
 969 &   Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
 970 &   Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,&
 971 &   Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
 972 &   k0,Dtset%spnorbscl,Cryst%ucvol,dtset%charge,ks_vtrial,ks_vxc,Cryst%xred,&
 973 &   nucdipmom=Dtset%nucdipmom)
 974 
 975     !  Symmetrize KS Dij
 976 #if 0
 977    call symdij(Cryst%gprimd,Cryst%indsym,ipert,&
 978 &   Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,0,KS_paw_ij,Pawang,&
 979 &   Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
 980 #else
 981    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert,&
 982 &   Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,KS_paw_ij,Pawang,&
 983 &   Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
 984 #endif
 985 !
 986 !  Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
 987    call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab)
 988    call timab(561,2,tsec)
 989  end if
 990 
 991  call timab(406,2,tsec) ! make_vhxc
 992 
 993  !=== Calculate Vxc(b1,b2,k,s)=<b1,k,s|v_{xc}|b2,k,s>  for all the states included in GW ===
 994  !  * ks_vxcvalme is calculated without NLCC, ks_vxcme contains NLCC (if any)
 995  !  * This part is parallelized within MPI_COMM_WORD since each node has all GW wavefunctions.
 996  !  * ks_vUme is zero unless we are using LDA+U as starting point, see calc_vHxc_braket
 997  !  * Note that vH matrix elements are calculated using the true uncutted interaction.
 998  call timab(407,1,tsec) ! vHxc_me
 999 
1000  call melflags_reset(KS_mflags)
1001  KS_mflags%has_vhartree=1
1002  KS_mflags%has_vxc     =1
1003  KS_mflags%has_vxcval  =1
1004  if (Dtset%usepawu>0     )  KS_mflags%has_vu     =1
1005  if (Dtset%useexexch>0   )  KS_mflags%has_lexexch=1
1006  if (Sigp%use_sigxcore==1)  KS_mflags%has_sxcore =1
1007  if (gwcalctyp<10           )  KS_mflags%only_diago =1 ! off-diagonal elements only for SC on wavefunctions.
1008 
1009  if (.FALSE.) then ! quick and dirty hack to test HF contribution.
1010    MSG_WARNING("testing on-site HF")
1011    lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
1012    ABI_MALLOC(dij_hf,(cplex_dij*lmn2_size_max,ndij,Cryst%natom))
1013    call paw_dijhf(ndij,cplex_dij,lmn2_size_max,Cryst%natom,Cryst%ntypat,Pawtab,Pawrad,Pawang,&
1014 &   KS_Pawrhoij,dij_hf,Dtset%prtvol)
1015 
1016    do iat=1,Cryst%natom
1017      itypat = Cryst%typat(iat)
1018      ii = Pawtab(itypat)%lmn2_size
1019      KS_Paw_ij(iat)%dijxc(:,:) = dij_hf(1:cplex_dij*ii,:,iat)
1020      KS_Paw_ij(iat)%dijxc_hat(:,:) = zero
1021    end do
1022    ABI_FREE(dij_hf)
1023 
1024 !  option_dij=3
1025 !  call symdij(Cryst%gprimd,Cryst%indsym,ipert,&
1026 !  &   Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
1027 !  &   KS_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,&
1028 !  &   Cryst%symafm,Cryst%symrec)
1029  end if
1030 
1031  ABI_MALLOC(tmp_kstab,(2,Wfd%nkibz,Wfd%nsppol))
1032  tmp_kstab=0
1033  do spin=1,Sigp%nsppol
1034    do ikcalc=1,Sigp%nkptgw ! No spin dependent!
1035      ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1036      tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin)
1037      tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin)
1038    end do
1039  end do
1040 
1041  call calc_vhxc_me(Wfd,KS_mflags,KS_me,Cryst,Dtset,nfftf,ngfftf,&
1042 & ks_vtrial,ks_vhartr,ks_vxc,Psps,Pawtab,KS_paw_an,Pawang,Pawfgrtab,KS_paw_ij,dijexc_core,&
1043 & ks_rhor,usexcnhat,ks_nhat,ks_nhatgr,nhatgrdim,tmp_kstab,taug=ks_taug,taur=ks_taur)
1044  ABI_FREE(tmp_kstab)
1045 
1046 !#ifdef DEV_HAVE_SCGW_SYM
1047 !Set KS matrix elements connecting different irreps to zero. Do not touch unknown bands!.
1048  if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then
1049    bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1050    ABI_MALLOC(ks_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
1051    ks_irreptab=0
1052    do spin=1,Sigp%nsppol
1053      do ikcalc=1,Sigp%nkptgw
1054        ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1055        first_band = Sigp%minbnd(ikcalc,spin)
1056        last_band  = Sigp%maxbnd(ikcalc,spin)
1057        if (.not.esymm_failed(KS_sym(ik_ibz,spin))) then
1058          ks_irreptab(first_band:last_band,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(first_band:last_band)
1059 !        ks_irreptab(bmin:bmax,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(bmin:bmax)
1060        end if
1061      end do
1062    end do
1063    call melements_zero(KS_me,ks_irreptab)
1064    ABI_FREE(ks_irreptab)
1065  end if
1066 !#endif
1067 
1068  call melements_print(KS_me,header="Matrix elements in the KS basis set",prtvol=Dtset%prtvol)
1069  !
1070  ! If possible, calculate the EXX energy from the between the frozen core
1071  ! and the valence electrons using KS wavefunctions
1072  ! MG: BE careful here, since ex_energy is meaningful only is all occupied states are calculated.
1073  if( KS_mflags%has_sxcore ==1 ) then
1074    ! TODO
1075    !ex_energy = mels_get_exene_core(KS_me,kmesh,KS_BSt)
1076    ex_energy=zero
1077    do spin=1,Sigp%nsppol
1078      do ik=1,Kmesh%nibz
1079        do ib=b1gw,b2gw
1080          if (Sigp%nsig_ab==1) then
1081            ex_energy = ex_energy + half*KS_BSt%occ(ib,ik,spin)*Kmesh%wt(ik)*KS_me%sxcore(ib,ib,ik,spin)
1082          else
1083            ex_energy = ex_energy + half*KS_BSt%occ(ib,ik,spin)*Kmesh%wt(ik)*SUM(KS_me%sxcore(ib,ib,ik,:))
1084          end if
1085        end do
1086      end do
1087    end do
1088    write(msg,'(a,2(es16.6,a))')' CORE Exchange energy with KS wavefunctions: ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV'
1089    call wrtout(std_out,msg,'COLL')
1090  end if
1091 
1092  call timab(407,2,tsec) ! vHxc_me
1093 
1094  call timab(408,1,tsec) ! hqp_init
1095 
1096  ! Do not break this coding! When gwcalctyp>10, the order of the bands can be interexchanged after
1097  ! the diagonalization. Therefore, we have to correctly assign the matrix elements to the corresponding
1098  ! bands and we cannot skip the following even though it looks unuseful.
1099  if (gwcalctyp>=10) then
1100    call wrtout(std_out,ch10//' *************** KS Energies *******************','COLL')
1101  end if
1102 
1103  !=== QP_BSt stores energies and occ. used for the calculation ===
1104  !  * Initialize QP_BSt with KS values.
1105  !  * In case of SC update QP_BSt using the QPS file.
1106  call ebands_copy(KS_BSt,QP_BSt)
1107 
1108  ABI_MALLOC(qp_rhor,(nfftf,Dtset%nspden))
1109  ABI_MALLOC(qp_taur,(nfftf,Dtset%nspden*Dtset%usekden))
1110  QP_sym => KS_sym
1111 
1112  if (gwcalctyp<10) then
1113    ! one-shot GW, just do a copy of the KS density.
1114    qp_rhor=ks_rhor
1115    if(Dtset%usekden==1)qp_taur=ks_taur
1116    QP_sym => KS_sym
1117  else
1118    ! Self-consistent GW.
1119    ! Read the unitary matrix and the QP energies of the previous step from the QPS file.
1120    call energies_init(QP_energies)
1121    QP_energies%e_corepsp=ecore/Cryst%ucvol
1122 
1123    !  m_lda_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}>
1124    ! Initialize the QP amplitudes with KS wavefunctions.
1125    ABI_MALLOC(Sr%m_lda_to_qp,(Sigp%nbnds,Sigp%nbnds,Kmesh%nibz,Sigp%nsppol))
1126    Sr%m_lda_to_qp=czero
1127    do ib=1,Sigp%nbnds
1128      Sr%m_lda_to_qp(ib,ib,:,:)=cone
1129    end do
1130 
1131    ! Now read m_lda_to_qp and update the energies in QP_BSt.
1132    ! TODO switch on the renormalization of n in sigma.
1133    ABI_MALLOC(prev_rhor,(nfftf,Dtset%nspden))
1134    ABI_MALLOC(prev_taur,(nfftf,Dtset%nspden*Dtset%usekden))
1135    ABI_DT_MALLOC(prev_Pawrhoij,(Cryst%natom*Psps%usepaw))
1136 
1137    call rdqps(QP_BSt,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,&
1138 &   nfftf,ngfftf,Cryst%ucvol,Dtset%paral_kgb,Cryst,Pawtab,MPI_enreg_seq,nbsc,Sr%m_lda_to_qp,prev_rhor,prev_Pawrhoij)
1139 
1140 !  Find the irreps associated to the QP amplitudes starting from the analogous table for the KS states.
1141 !  bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1142 !  allocate(qp_irreptab(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
1143 !  qp_irreptab=0
1144 !  !qp_irreptab=ks_irreptab
1145 
1146 !  do jb_qp=bmin,bmax
1147 !  do ib_ks=bmin,bmax
1148 !  if (ABS(Sr%m_lda_to_qp(ib_ks,jb_qp,ik_ibz,spin)) > tol12) then ! jb_qp has same the same character as ib_ks.
1149 !  ks_irr = ks_irreptab(ib_ks,ib_ks,ik_ibz,spin)
1150 !  qp_irreptab(jb_qp,jb_qp,ik_ibz,spin) = ks_irr
1151 !  do ii=bmin,bmax
1152 !  if (ks_irr == ks_irreptab(ii,ib_ks,ik_ibz,spin)) then
1153 !  qp_irreptab(jb_qp,ii,ik_ibz,spin) = ks_irr
1154 !  end if
1155 !  end do
1156 !  end if
1157 !  end do
1158 !  end do
1159 
1160    if (nscf==0) prev_rhor=ks_rhor
1161    if (nscf==0 .and. Dtset%usekden==1) prev_taur=ks_taur
1162 
1163    if (nscf>0.and.gwcalctyp>=20.and.wfd_iam_master(Wfd)) then
1164      ! Print the unitary transformation on std_out.
1165      call show_QP(QP_BSt,Sr%m_lda_to_qp,fromb=Sigp%minbdgw,tob=Sigp%maxbdgw,unit=std_out,tolmat=0.001_dp)
1166    end if
1167 
1168    ! Compute QP wfg as linear combination of KS states ===
1169    !  * Wfd%ug is modified inside calc_wf_qp
1170    !  * For PAW, update also the on-site projections.
1171    !  * WARNING the first dimension of MPI_enreg MUST be Kmesh%nibz
1172    !  TODO here we should use nbsc instead of nbnds
1173 
1174    call wfd_rotate(Wfd,Cryst,Sr%m_lda_to_qp)
1175 
1176    ! Reinit the storage mode of Wfd as ug have been changed ===
1177    ! Update also the wavefunctions for GW corrections on each processor
1178    call wfd_reset_ur_cprj(Wfd)
1179 
1180    ! This test has been disabled (too expensive!)
1181    if (.False.) call wfd_test_ortho(Wfd, Cryst, Pawtab, unit=std_out)
1182 
1183    ! Compute QP occupation numbers.
1184    call wrtout(std_out,'sigma: calculating QP occupation numbers:','COLL')
1185    call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=0)
1186    qp_vbik(:,:) = get_valence_idx(QP_BSt)
1187 
1188 !  #ifdef DEV_HAVE_SCGW_SYM
1189    ! Calculate the irreducible representations of the new QP amplitdues.
1190    if (Sigp%symsigma==1.and.gwcalctyp>=20) then
1191      ABI_DT_MALLOC(QP_sym,(Wfd%nkibz,Wfd%nsppol))
1192      use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric
1193      do spin=1,Wfd%nsppol
1194        do ikcalc=1,Sigp%nkptgw
1195          ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1196 !        Quick fix for SCGW+symm TODO fix properly!
1197          first_band = Sigp%minbnd(ikcalc,spin)
1198          last_band  = Sigp%maxbnd(ikcalc,spin)
1199 !        first_band = MINVAL(Sigp%minbnd(:,spin))
1200 !        last_band  = MAXVAL(Sigp%maxbnd(:,spin))
1201 !        call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,QP_BSt,Pawtab,Pawrad,Pawang,Psps,&
1202          call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,QP_BSt,Pawtab,Pawrad,Pawang,Psps,&
1203 &         Dtset%tolsym,QP_sym(ik_ibz,spin))
1204        end do
1205      end do
1206 
1207      ! Recreate the Sig_ij tables taking advantage of the classification of the bands.
1208      call sigma_tables(Sigp,Kmesh,QP_sym)
1209    end if
1210 !  #endif
1211 
1212    ! Compute QP density using the updated wfg.
1213    call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,qp_rhor)
1214    if (Dtset%usekden==1) then
1215      call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,qp_taur,optcalc=1)
1216    end if
1217 
1218    ! ========================================
1219    ! ==== QP self-consistent GW with PAW ====
1220    ! ========================================
1221    if (Dtset%usepaw==1) then
1222      ABI_MALLOC(qp_nhat,(nfftf,Dtset%nspden))
1223      nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat
1224      ABI_MALLOC(qp_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim))
1225 
1226      ABI_DT_MALLOC(QP_pawrhoij,(Cryst%natom))
1227      ABI_DT_MALLOC(QP_paw_ij,(Cryst%natom))
1228      ABI_DT_MALLOC(QP_paw_an,(Cryst%natom))
1229 
1230      ! Calculate new QP quantities: nhat, nhatgr, rho_ij, paw_ij, and paw_an.
1231      call paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,QP_BSt,&
1232 &     Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij,&
1233 &     QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim,qp_nhatgr,compch_sph,compch_fft)
1234    else
1235      ABI_MALLOC(qp_nhat,(0,0))
1236      ABI_MALLOC(qp_nhatgr,(0,0,0))
1237      ABI_DT_MALLOC(QP_pawrhoij,(0))
1238      ABI_DT_MALLOC(QP_paw_ij,(0))
1239      ABI_DT_MALLOC(QP_paw_an,(0))
1240    end if
1241 
1242 !  here I should renormalize the density
1243    call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,qp_rhor,Cryst%ucvol,&
1244 &   Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
1245 
1246    if (Dtset%usepaw==1) qp_rhor(:,:)=qp_rhor(:,:)+qp_nhat(:,:) ! Add the "hat" term.
1247 
1248    call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_rhor,ucvol=ucvol)
1249    if(Dtset%usekden==1) then
1250      call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_taur,optrhor=1,ucvol=ucvol)
1251    end if
1252 
1253    ! Simple mixing of the PW density to damp oscillations in the Hartree potential.
1254    if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then
1255      write(msg,'(2a,f6.3)')ch10,' sigma: mixing QP densities using rhoqpmix= ',Dtset%rhoqpmix
1256      call wrtout(std_out,msg,'COLL')
1257      qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor)
1258      if(Dtset%usekden==1)qp_taur = prev_taur + Dtset%rhoqpmix*(qp_taur-prev_taur) ! mix taur.
1259    end if
1260 
1261    ABI_FREE(prev_rhor)
1262    ABI_FREE(prev_taur)
1263    if (Psps%usepaw==1.and.nscf>0) then
1264      call pawrhoij_free(prev_pawrhoij)
1265    end if
1266    ABI_DT_FREE(prev_pawrhoij)
1267 
1268    ABI_MALLOC(qp_rhog,(2,nfftf))
1269    ABI_MALLOC(qp_taug,(2,nfftf*Dtset%usekden))
1270    call fourdp(1,qp_rhog,qp_rhor(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp5)
1271    if(Dtset%usekden==1)call fourdp(1,qp_taug,qp_taur(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp5)
1272 
1273    ! ===========================================
1274    ! ==== Optional output of the QP density ====
1275    ! ===========================================
1276    if (Dtset%prtden/=0.and.wfd_iam_master(Wfd)) then
1277      call fftdatar_write("qp_rhor",dtfil%fnameabo_qp_den,dtset%iomode,hdr_sigma,&
1278      cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor,mpi_enreg_seq,ebands=qp_bst)
1279    end if
1280 
1281    ! ===========================================
1282    ! === Optional output of the full QP density
1283    ! ===========================================
1284    if (Wfd%usepaw==1.and.Dtset%prtden==2) then
1285      ABI_MALLOC(qp_rhor_paw   ,(nfftf,Wfd%nspden))
1286      ABI_MALLOC(qp_rhor_n_one ,(nfftf,Wfd%nspden))
1287      ABI_MALLOC(qp_rhor_nt_one,(nfftf,Wfd%nspden))
1288 
1289      call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,qp_nhat,&
1290 &     Wfd%nspinor,Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,QP_pawrhoij,Pawtab,Dtset%prtvol,&
1291 &     qp_rhor,qp_rhor_paw,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
1292 
1293      ABI_FREE(qp_rhor_n_one)
1294      ABI_FREE(qp_rhor_nt_one)
1295      if (Dtset%prtvol>9) then ! Print a normalisation check
1296        norm = SUM(qp_rhor_paw(:,1))*Cryst%ucvol/PRODUCT(Pawfgr%ngfft(1:3))
1297        write(msg,'(a,F8.4)') '  QUASIPARTICLE DENSITY CALCULATED - NORM OF DENSITY: ',norm
1298        call wrtout(std_out,msg,'PERS')
1299      end if
1300 
1301      ! Write the density to file
1302      if (my_rank==master) then
1303        call fftdatar_write("qp_pawrhor",dtfil%fnameabo_qp_pawden,dtset%iomode,hdr_sigma,&
1304        cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor_paw,mpi_enreg_seq,ebands=qp_bst)
1305      end if
1306      ABI_FREE(qp_rhor_paw)
1307    end if
1308 
1309    nkxc=0
1310    if (Dtset%nspden==1) nkxc=2
1311    if (Dtset%nspden>=2) nkxc=3 !check GGA and spinor that is messy !!!
1312    !In case of MGGA, fxc and kxc are not available and we dont need them
1313    !for the screening part (for now ...)
1314    if (Dtset%ixc<0.and.libxc_functionals_ismgga()) nkxc=0
1315    if (nkxc/=0)  then
1316      ABI_MALLOC(qp_kxc,(nfftf,nkxc))
1317    end if
1318 
1319    ! **** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
1320    n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf
1321    ABI_MALLOC(qp_vhartr,(nfftf))
1322    ABI_MALLOC(qp_vtrial,(nfftf,Dtset%nspden))
1323    ABI_MALLOC(qp_vxc,(nfftf,Dtset%nspden))
1324 
1325    optene=4; moved_atm_inside=0; moved_rhor=0; istep=1
1326 
1327    call setvtr(Cryst%atindx1,Dtset,QP_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
1328 &   istep,qp_kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
1329 &   Cryst%nattyp,nfftf,ngfftf,ngrvdw,qp_nhat,qp_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
1330 &   optene,pawrad,Pawtab,ph1df,Psps,qp_rhog,qp_rhor,Cryst%rmet,Cryst%rprimd,strsxc,&
1331 &   Cryst%ucvol,usexcnhat,qp_vhartr,vpsp,qp_vtrial,qp_vxc,vxcavg_qp,Wvl,&
1332 &   xccc3d,Cryst%xred,taug=qp_taug,taur=qp_taur)
1333 
1334    if (allocated(qp_kxc)) then
1335      ABI_FREE(qp_kxc)
1336    end if
1337 
1338    if (Dtset%usepaw==1) then
1339      call timab(561,1,tsec)
1340 
1341      ! Compute QP Dij
1342      ipert=0; idir=0
1343      call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert,&
1344 &     Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
1345 &     Dtset%nspden,Cryst%ntypat,QP_paw_an,QP_paw_ij,Pawang,Pawfgrtab,&
1346 &     Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
1347 &     k0,Dtset%spnorbscl,Cryst%ucvol,dtset%charge,qp_vtrial,qp_vxc,Cryst%xred,&
1348 &     nucdipmom=Dtset%nucdipmom)
1349 
1350      ! Symmetrize total Dij
1351      option_dij=0
1352 
1353 #if 0
1354      call symdij(Cryst%gprimd,Cryst%indsym,ipert,&
1355 &     Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
1356 &     QP_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,&
1357 &     Cryst%symafm,Cryst%symrec)
1358 #else
1359      call symdij_all(Cryst%gprimd,Cryst%indsym,ipert,&
1360 &     Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
1361 &     QP_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,&
1362 &     Cryst%symafm,Cryst%symrec)
1363 #endif
1364 
1365      ! Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij.
1366      call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab)
1367      call timab(561,2,tsec)
1368    end if
1369 
1370    ehartree=half*SUM(qp_rhor(:,1)*qp_vhartr(:))/DBLE(nfftf)*Cryst%ucvol
1371 
1372    write(msg,'(a,80a)')ch10,('-',ii=1,80)
1373    call wrtout(ab_out,msg,'COLL')
1374    write(msg,'(5a,f9.4,3a,es21.14,2a,es21.14)')ch10,&
1375 &   ' QP results after the unitary transformation in the KS subspace: ',ch10,ch10,&
1376 &   '  Number of electrons    = ',qp_rhog(1,1)*Cryst%ucvol,ch10,ch10,&
1377 &   '  QP Band energy    [Ha] = ',get_bandenergy(QP_BSt),ch10,&
1378 &   '  QP Hartree energy [Ha] = ',ehartree
1379    call wrtout(ab_out,msg,'COLL')
1380    write(msg,'(a,80a)')ch10,('-',ii=1,80)
1381    call wrtout(ab_out,msg,'COLL')
1382 !
1383 !  TODO Since plasmonpole model 2-3-4 depend on the Fourier components of the density
1384 !  in case of self-consistency we might calculate here the ppm coefficients using qp_rhor
1385  end if ! gwcalctyp>=10
1386 !
1387 !=== KS hamiltonian hlda(b1,b1,k,s)= <b1,k,s|H_s|b1,k,s> ===
1388  ABI_MALLOC(hlda,(b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab))
1389  hlda = czero
1390 
1391  if (Dtset%nspinor==1) then
1392    do spin=1,Sigp%nsppol
1393      do ik=1,Kmesh%nibz
1394        do ib=b1gw,b2gw
1395          hlda(ib,ib,ik,spin) = KS_BSt%eig(ib,ik,spin)
1396        end do
1397      end do
1398    end do
1399  else
1400    ! Spinorial case
1401    !  * Note that here vxc contains the contribution of the core.
1402    !  * Scale ovlp if orthonormalization is not satisfied as npwwfn might be < npwvec.
1403    !  TODO add spin-orbit case and gwpara 2
1404    ABI_MALLOC(my_band_list, (wfd%mband))
1405    ABI_MALLOC(bmask, (wfd%mband))
1406    bmask = .False.; bmask(b1gw:b2gw) = .True.
1407 
1408    if (Wfd%usepaw==1) then
1409      ABI_DT_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor))
1410      call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm)
1411    end if
1412 
1413    do spin=1,Sigp%nsppol
1414      do ik_ibz=1,Kmesh%nibz
1415        ! Distribute bands in [b1gw, b2gw] range
1416        call wfd_distribute_bands(wfd, ik_ibz, spin, my_nband, my_band_list, bmask=bmask)
1417        if (my_nband == 0) cycle
1418        npw_k = Wfd%npwarr(ik_ibz)
1419        do ii=1,my_nband  ! ib=b1gw,b2gw in sequential
1420          ib = my_band_list(ii)
1421          ug1  => Wfd%Wave(ib,ik_ibz,spin)%ug
1422          cdummy = xdotc(npw_k*Wfd%nspinor,ug1,1,ug1,1)
1423          ovlp(1) = REAL(cdummy); ovlp(2) = AIMAG(cdummy)
1424 
1425          if (Psps%usepaw==1) then
1426            call wfd_get_cprj(Wfd,ib,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.)
1427            ovlp = ovlp + paw_overlap(Cp1,Cp1,Cryst%typat,Pawtab)
1428          end if
1429 !        write(std_out,*)ovlp(1),ovlp(2)
1430          norm=DBLE(ovlp(1)+ovlp(2))
1431          ovlp(1)=DBLE(ovlp(1)/norm); ovlp(2)=DBLE(ovlp(2)/norm) ! ovlp(2)=cone-ovlp(1)
1432          hlda(ib,ib,ik_ibz,1) = KS_BSt%eig(ib,ik_ibz,1)*ovlp(1)-KS_me%vxc(ib,ib,ik_ibz,3)
1433          hlda(ib,ib,ik_ibz,2) = KS_BSt%eig(ib,ik_ibz,1)*ovlp(2)-KS_me%vxc(ib,ib,ik_ibz,4)
1434          hlda(ib,ib,ik_ibz,3) = KS_me%vxc(ib,ib,ik_ibz,3)
1435          hlda(ib,ib,ik_ibz,4) = KS_me%vxc(ib,ib,ik_ibz,4)
1436        end do
1437      end do
1438    end do
1439 
1440    call xmpi_sum(hlda, wfd%comm, ierr)
1441 
1442    ABI_FREE(my_band_list)
1443    ABI_FREE(bmask)
1444    if (Wfd%usepaw==1) then
1445      call pawcprj_free(Cp1)
1446      ABI_DT_FREE(Cp1)
1447    end if
1448  end if
1449 
1450 !=== Initialize Sigma results ===
1451 !TODO it is better if we use ragged arrays indexed by the k-point
1452  call sigma_init(Sigp,Kmesh%nibz,Dtset%usepawu,Sr)
1453 
1454  ! Setup of the bare Hamiltonian := T + v_{loc} + v_{nl} + v_H.
1455  ! * The representation depends wheter we are updating the wfs or not.
1456  ! * ks_vUme is zero unless we are using LDA+U as starting point, see calc_vHxc_braket
1457  ! * Note that vH matrix elements are calculated using the true uncutted interaction.
1458 
1459  if (gwcalctyp<10) then
1460    ! For one-shot GW use the KS representation.
1461    Sr%hhartree=hlda-KS_me%vxcval
1462    ! Additional goodies for PAW
1463    !  * LDA +U Hamiltonian
1464    !  * LEXX.
1465    !  * Core contribution estimated using Fock exchange.
1466    if (Dtset%usepaw==1) then
1467      if (Sigp%use_sigxcore==1) Sr%hhartree=hlda - (KS_me%vxc - KS_me%sxcore)
1468      if (Dtset%usepawu>0) Sr%hhartree=Sr%hhartree-KS_me%vu
1469      if (Dtset%useexexch>0) then
1470        MSG_ERROR("useexexch > 0 not implemented")
1471        Sr%hhartree = Sr%hhartree - KS_me%vlexx
1472      end if
1473    end if
1474  else
1475    ! Self-consistent on energies and|or wavefunctions.
1476    !   * For NC get the bare Hamiltonian  $H_{bare}= T+v_{loc}+ v_{nl}$ in the KS representation
1477    !   * For PAW, calculate the matrix elements of h0, store also the new Dij in QP_Paw_ij.
1478    !   * h0 is defined as T+vH[tn+nhat+tnZc] + vxc[tnc] + dij_eff and
1479    !     dij_eff = dij^0 + dij^hartree + dij^xc-dij^xc_val + dijhat - dijhat_val.
1480    !     In the above expression tn, tnhat are QP quantities.
1481    if (Dtset%usepaw==0) then
1482      ABI_MALLOC(hbare,(b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab))
1483      hbare=hlda-KS_me%vhartree-KS_me%vxcval
1484 
1485      ! Change basis from KS to QP, hbare is overwritten: A_{QP} = U^\dagger A_{KS} U
1486      ABI_MALLOC(htmp,(b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab))
1487      ABI_MALLOC(ctmp,(b1gw:b2gw,b1gw:b2gw))
1488      ABI_MALLOC(uks2qp,(b1gw:b2gw,b1gw:b2gw))
1489 
1490      htmp=hbare; hbare=czero
1491 
1492      do spin=1,Sigp%nsppol
1493        do ik=1,Kmesh%nibz
1494          uks2qp(:,:) = Sr%m_lda_to_qp(b1gw:b2gw,b1gw:b2gw,ik,spin)
1495          do iab=1,Sigp%nsig_ab
1496            is_idx=spin; if (Sigp%nsig_ab>1) is_idx=iab
1497            ctmp(:,:)=MATMUL(htmp(:,:,ik,is_idx),uks2qp(:,:))
1498            hbare(:,:,ik,is_idx)=MATMUL(TRANSPOSE(CONJG(uks2qp)),ctmp)
1499          end do
1500        end do
1501      end do
1502      ABI_FREE(htmp)
1503      ABI_FREE(ctmp)
1504      ABI_FREE(uks2qp)
1505    end if ! usepaw==0
1506 
1507    ! Calculate the QP matrix elements
1508    ! This part is parallelized within MPI_COMM_WORD since each node has all GW wavefunctions.
1509    ! For PAW, we have to construct the new bare Hamiltonian.
1510    call wrtout(std_out,ch10//' *************** QP Energies *******************','COLL')
1511 
1512    call melflags_reset(QP_mflags)
1513 !  if (gwcalctyp<20) QP_mflags%only_diago=1 ! For e-only, no need of off-diagonal elements.
1514    QP_mflags%has_vhartree=1
1515    if (Dtset%usepaw==1)    QP_mflags%has_hbare  =1
1516 !  QP_mflags%has_vxc     =1
1517 !  QP_mflags%has_vxcval  =1
1518 !  if (Sigp%gwcalctyp >100) QP_mflags%has_vxcval_hybrid=1
1519    if (mod10==5 .and. &
1520 &   (Dtset%ixc_sigma==-402 .or. Dtset%ixc_sigma==-406 .or. Dtset%ixc_sigma==-427 .or. Dtset%ixc_sigma==-428 .or.&
1521 &   Dtset%ixc_sigma==-456 .or. Dtset%ixc_sigma==41 .or. Dtset%ixc_sigma==42)) then
1522      QP_mflags%has_vxcval_hybrid=1
1523    end if
1524 !  if (Sigp%use_sigxcore==1) QP_mflags%has_sxcore =1
1525 !  if (Dtset%usepawu>0)    QP_mflags%has_vu     =1
1526 !  if (Dtset%useexexch>0)  QP_mflags%has_lexexch=1
1527 
1528    ABI_MALLOC(tmp_kstab,(2,Wfd%nkibz,Wfd%nsppol))
1529    tmp_kstab=0
1530    do spin=1,Sigp%nsppol
1531      do ikcalc=1,Sigp%nkptgw ! No spin dependent!
1532        ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1533        tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin)
1534        tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin)
1535      end do
1536    end do
1537 
1538    call calc_vhxc_me(Wfd,QP_mflags,QP_me,Cryst,Dtset,nfftf,ngfftf,&
1539 &   qp_vtrial,qp_vhartr,qp_vxc,Psps,Pawtab,QP_paw_an,Pawang,Pawfgrtab,QP_paw_ij,dijexc_core,&
1540 &   qp_rhor,usexcnhat,qp_nhat,qp_nhatgr,nhatgrdim,tmp_kstab,taug=qp_taug,taur=qp_taur)
1541    ABI_FREE(tmp_kstab)
1542 
1543 !  #ifdef DEV_HAVE_SCGW_SYM
1544    if (gwcalctyp>=20 .and. Sigp%symsigma>0) then
1545      bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1546      ABI_MALLOC(qp_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
1547      qp_irreptab=0
1548      do spin=1,Sigp%nsppol
1549        do ikcalc=1,Sigp%nkptgw
1550          ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1551          first_band = Sigp%minbnd(ikcalc,spin)
1552          last_band  = Sigp%maxbnd(ikcalc,spin)
1553          if (.not.esymm_failed(QP_sym(ik_ibz,spin))) then
1554            qp_irreptab(first_band:last_band,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(first_band:last_band)
1555 !          qp_irreptab(bmin:bmax,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(bmin:bmax)
1556          end if
1557        end do
1558      end do
1559      call melements_zero(QP_me,qp_irreptab)
1560      ABI_FREE(qp_irreptab)
1561    end if
1562 !  #endif
1563 
1564    call melements_print(QP_me,header="Matrix elements in the QP basis set",prtvol=Dtset%prtvol)
1565 
1566    ! * Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij.
1567    if (Dtset%usepaw==1) then
1568      call wrtout(std_out," *** After calc_vHxc_braket *** ",'COLL')
1569 !    TODO terminate the implementation of this routine.
1570      call paw_ij_print(QP_Paw_ij,unit=std_out,pawprtvol=Dtset%pawprtvol,pawspnorb=Dtset%pawspnorb,mode_paral="COLL")
1571      call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab)
1572    end if
1573 
1574    if (Dtset%usepaw==0) then
1575 !    GA : We have an odd bug here. I have to unroll this loop, otherwise it
1576 !    might cause segfault when running on several nodes.
1577 !
1578 !    Sr%hhartree = hbare + QP_me%vhartree
1579      if(QP_mflags%has_vxcval_hybrid==0) then
1580        do spin=1,Sigp%nsppol*Sr%nsig_ab
1581          do ikcalc=1,Sr%nkibz
1582            do ib1=b1gw,b2gw
1583              do ib2=b1gw,b2gw
1584                Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + QP_me%vhartree(ib2,ib1,ikcalc,spin)
1585              end do
1586            end do
1587          end do
1588        end do
1589      else
1590        do spin=1,Sigp%nsppol*Sr%nsig_ab
1591          do ikcalc=1,Sr%nkibz
1592            do ib1=b1gw,b2gw
1593              do ib2=b1gw,b2gw
1594                Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + &
1595 &               QP_me%vhartree(ib2,ib1,ikcalc,spin) + &
1596 &               QP_me%vxcval_hybrid(ib2,ib1,ikcalc,spin)
1597              end do
1598            end do
1599          end do
1600        end do
1601      end if
1602    else
1603      Sr%hhartree=QP_me%hbare
1604    end if
1605 
1606 !  #ifdef DEV_HAVE_SCGW_SYM
1607    if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then
1608 !    bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1609      do spin=1,Sigp%nsppol
1610        do ik_ibz=1,Kmesh%nibz
1611          if (.not.esymm_failed(QP_sym(ik_ibz,spin))) then
1612            bmin=Sigp%minbnd(ik_ibz,spin); bmax=Sigp%minbnd(ik_ibz,spin)
1613            do ib2=bmin,bmax
1614              irr_idx2 = QP_sym(ik_ibz,spin)%b2irrep(ib2)
1615              do ib1=bmin,bmax
1616                irr_idx1 = QP_sym(ik_ibz,spin)%b2irrep(ib1)
1617                if (irr_idx1/=irr_idx2 .and. ALL((/irr_idx1,irr_idx2/)/=0) ) Sr%hhartree(ib1,ib2,ik_ibz,spin) = czero
1618              end do
1619            end do
1620          end if
1621        end do
1622      end do
1623    end if
1624 !  #endif
1625 
1626    ABI_FREE(qp_rhog)
1627    ABI_FREE(qp_vhartr)
1628    ABI_FREE(qp_vxc)
1629    ABI_FREE(qp_taug)
1630    ABI_FREE(qp_nhat)
1631    ABI_FREE(qp_nhatgr)
1632    call melements_free(QP_me)
1633  end if ! gwcalctyp<10
1634 
1635  ! Free some memory
1636  if (allocated(hbare)) then
1637    ABI_FREE(hbare)
1638  end if
1639  if (allocated(hlda )) then
1640    ABI_FREE(hlda)
1641  end if
1642 
1643  ! Prepare the storage of QP amplitudes and energies ===
1644  ! Initialize with KS wavefunctions and energies.
1645  Sr%eigvec_qp=czero;  Sr%en_qp_diago=zero
1646  do ib=1,Sigp%nbnds
1647    Sr%en_qp_diago(ib,:,:)=KS_BSt%eig(ib,:,:)
1648    Sr%eigvec_qp(ib,ib,:,:)=cone
1649  end do
1650 
1651  ! Store <n,k,s|V_xc[n_val]|n,k,s> and <n,k,s|V_U|n,k,s> ===
1652  ! Note that we store the matrix elements of V_xc in the KS basis set, not in the QP basis set
1653  ! Matrix elements of V_U are zero unless we are using LDA+U as starting point
1654  do ib=b1gw,b2gw
1655    Sr%vxcme(ib,:,:)=KS_me%vxcval(ib,ib,:,:)
1656    if (Dtset%usepawu>0) Sr%vUme (ib,:,:)=KS_me%vu(ib,ib,:,:)
1657  end do
1658 
1659  ! Initial guess for the GW energies
1660  ! Save the energies of the previous iteration.
1661  do spin=1,Sigp%nsppol
1662    do ik=1,Kmesh%nibz
1663      do ib=1,Sigp%nbnds
1664        Sr%e0 (ib,ik,spin)=QP_BSt%eig(ib,ik,spin)
1665        Sr%egw(ib,ik,spin)=QP_BSt%eig(ib,ik,spin)
1666      end do
1667      Sr%e0gap(ik,spin)=zero
1668      ks_iv=ks_vbik(ik,spin)
1669      if (Sigp%nbnds>=ks_iv+1) Sr%e0gap(ik,spin)=Sr%e0(ks_iv+1,ik,spin)-Sr%e0(ks_iv,ik,spin)
1670    end do
1671  end do
1672 !
1673 !=== If required apply a scissor operator or update the energies ===
1674 !TODO check if other Sr entries have to be updated
1675 !moreover this part should be done only in case of semiconductors
1676 !FIXME To me it makes more sense if we apply the scissor to KS_BS but I have to RECHECK csigme
1677  if (ABS(Sigp%mbpt_sciss)>tol6) then
1678    write(msg,'(6a,f10.5,2a)')ch10,&
1679 &   ' sigma : performing a first self-consistency',ch10,&
1680 &   '  update of the energies in G by a scissor operator',ch10, &
1681 &   '  applying a scissor operator of ',Sigp%mbpt_sciss*Ha_eV,' [eV] ',ch10
1682    call wrtout(std_out,msg,'COLL')
1683    call wrtout(ab_out,msg,'COLL')
1684    do spin=1,Sigp%nsppol
1685      do ik=1,Kmesh%nibz
1686        ks_iv=ks_vbik(ik,spin)
1687        if (Sigp%nbnds>=ks_iv+1) then
1688          Sr%egw    (ks_iv+1:Sigp%nbnds,ik,spin) = Sr%egw    (ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss
1689          QP_BSt%eig(ks_iv+1:Sigp%nbnds,ik,spin) = QP_BSt%eig(ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss
1690        end if
1691      end do
1692    end do
1693 !  %call apply_scissor(QP_BSt,Sigp%mbpt_sciss)
1694  else if (.FALSE.) then
1695    write(msg,'(4a)')ch10,&
1696 &   ' sigma : performing a first self-consistency',ch10,&
1697 &   '  update of the energies in G by a previous GW calculation'
1698    call wrtout(std_out,msg,'COLL')
1699    call wrtout(ab_out,msg,'COLL')
1700 !  TODO Recheck this part, is not clear to me!
1701    ABI_MALLOC(igwene,(QP_Bst%mband,QP_Bst%nkpt,QP_Bst%nsppol))
1702    call rdgw(QP_BSt,'__in.gw__',igwene,extrapolate=.TRUE.)
1703    ABI_FREE(igwene)
1704    Sr%egw=QP_BSt%eig
1705    !
1706    ! * Recalculate the new fermi level.
1707    call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=0)
1708  end if
1709 !
1710 !In case of AC refer all the energies wrt to the fermi level
1711 !Take care because results from ppmodel cannot be used for AC
1712 !FIXME check ks_energy or qp_energy (in case of SCGW?)
1713 
1714  if (mod10==SIG_GW_AC) then
1715    ! All these quantities will be passed to csigme
1716    ! if I skipped the self-consistent part then here I have to use fermi
1717    QP_BSt%eig = QP_BSt%eig -QP_BSt%fermie
1718    Sr%egw = Sr%egw-QP_BSt%fermie
1719    Sr%e0  = Sr%e0 -QP_BSt%fermie
1720    oldefermi=QP_BSt%fermie
1721    ! TODO Recheck fermi
1722    ! Clean EVERYTHING in particulare the treatment of E fermi
1723    QP_BSt%fermie=zero
1724  end if
1725  !
1726  ! === Setup frequencies around the KS\QP eigenvalues to compute Sigma derivatives (notice the spin) ===
1727  ! TODO it is better using an odd Sr%nomega4sd so that the KS\QP eigenvalue is in the middle
1728  ioe0j=Sr%nomega4sd/2+1
1729  do spin=1,Sigp%nsppol
1730    do ikcalc=1,Sigp%nkptgw
1731      ib1=Sigp%minbnd(ikcalc,spin)
1732      ib2=Sigp%maxbnd(ikcalc,spin)
1733      ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1734      do jb=ib1,ib2
1735        do io=1,Sr%nomega4sd
1736          Sr%omega4sd(jb,ik_ibz,io,spin)=Sr%egw(jb,ik_ibz,spin)+Sigp%deltae*(io-ioe0j)
1737        end do
1738      end do
1739    end do
1740  end do
1741 
1742  call timab(408,2,tsec) ! hqp_init
1743  call timab(409,1,tsec) ! getW
1744 
1745  ! === Get epsilon^{-1} either from the _SCR or the _SUSC file and store it in Er%epsm1 ===
1746  !   * If Er%mqmem==0, allocate and read a single q-slice inside csigme.
1747  !     TODO Er%nomega should be initialized so that only the frequencies really needed are stored in memory
1748 
1749  ! TODO: The same piece of code is present in screening.
1750  if (sigma_needs_w(Sigp)) then
1751    select case (dtset%gwgamma)
1752    case (0)
1753      id_required=4; ikxc=0; approx_type=0; option_test=0; dim_kxcg=0
1754      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1755 
1756    case (1, 2)
1757      ! ALDA TDDFT kernel vertex
1758      !ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available")
1759      ABI_CHECK(Er%ID==0,"Er%ID should be 0")
1760 
1761      if (Dtset%usepaw==1) then
1762        ! If we have PAW, we need the full density on the fine grid
1763        ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1764        if (Dtset%getpawden==0 .and. Dtset%irdpawden==0) then
1765          MSG_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1766        end if
1767        call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin),"COLL")
1768        if (.not. file_exists(dtfil%filpawdensin)) then
1769          MSG_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1770        end if
1771 
1772        ABI_DT_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1773        call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1774        ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm)
1775 
1776        call hdr_free(hdr_rhor)
1777        call pawrhoij_free(tmp_pawrhoij)
1778        ABI_DT_FREE(tmp_pawrhoij)
1779      end if ! Dtset%usepaw==1
1780 
1781      id_required=4; ikxc=7; approx_type=1; dim_kxcg=1
1782      if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1783      if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1784      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1785 
1786      dbg_mode=.FALSE.
1787      if (Dtset%usepaw==1) then
1788        ! Use PAW all-electron density
1789        call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_aepaw_rhor,&
1790        Er%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode)
1791      else
1792        ! Norm-conserving
1793        call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_rhor,&
1794        Er%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode)
1795      end if
1796 
1797    case (3, 4)
1798      ! ADA non-local kernel vertex
1799      !ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available")
1800      ABI_CHECK(Er%ID==0,"Er%ID should be 0")
1801      ABI_CHECK(Sigp%nsppol==1,"ADA vertex for GWGamma not available yet for spin-polarised cases")
1802 
1803      if (Dtset%usepaw==1) then ! If we have PAW, we need the full density on the fine grid
1804        ABI_MALLOC(ks_aepaw_rhor,(nfftf,Sigp%nsppol))
1805        if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then
1806          MSG_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1807        end if
1808        call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin),"COLL")
1809        if (.not. file_exists(dtfil%filpawdensin)) then
1810          MSG_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1811        end if
1812 
1813        ABI_DT_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1814 
1815        call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1816        ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm)
1817 
1818        call hdr_free(hdr_rhor)
1819        call pawrhoij_free(tmp_pawrhoij)
1820        ABI_DT_FREE(tmp_pawrhoij)
1821      end if ! Dtset%usepaw==1
1822 
1823      id_required=4; ikxc=7; approx_type=2;
1824      if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1825      if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1826      ABI_MALLOC(fxc_ADA,(Er%npwe,Er%npwe,Er%nqibz))
1827 !    Use userrd to set kappa
1828      if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp
1829 !    Set correct value of kappa (should be scaled with alpha*r_s where)
1830 !    r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3)
1831      rhoav = (drude_plsmf*drude_plsmf)/four_pi
1832      r_s = (three/(four_pi*rhoav))**third
1833      alpha = (four*ninth*piinv)**third
1834      Dtset%userrd = Dtset%userrd/(alpha*r_s)
1835 
1836      dbg_mode=.TRUE.
1837      if (Dtset%usepaw==1) then
1838        ! Use PAW all-electron density
1839        call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,&
1840        ks_aepaw_rhor,Er%npwe,Er%nqibz,Er%qibz,&
1841        fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode)
1842      else
1843        ! Norm conserving
1844        call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,&
1845        ks_rhor,Er%npwe,Er%nqibz,Er%qibz,&
1846        fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode)
1847      end if
1848 
1849      dim_kxcg = 0
1850      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1851 
1852 !  @WC: bootstrap --
1853    case (-3, -4, -5, -6, -7, -8)
1854 !    ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available")
1855      ABI_CHECK(Er%ID==0,"Er%ID should be 0")
1856 
1857      if (Dtset%usepaw==1) then
1858        ! If we have PAW, we need the full density on the fine grid
1859        ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1860        if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then
1861          MSG_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1862        end if
1863        call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin),"COLL")
1864        if (.not. file_exists(dtfil%filpawdensin)) then
1865          MSG_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1866        end if
1867 
1868        ABI_DT_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1869 
1870        call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1871        ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm)
1872 
1873        call hdr_free(hdr_rhor)
1874        call pawrhoij_free(tmp_pawrhoij)
1875        ABI_DT_FREE(tmp_pawrhoij)
1876      end if ! Dtset%usepaw==1
1877 
1878      id_required=4; ikxc=7; dim_kxcg=0
1879 
1880      if (dtset%gwgamma>-5) then
1881        approx_type=4  ! full fxc(G,G')
1882      else if (dtset%gwgamma>-7) then
1883        approx_type=5  ! fxc(0,0) one-shot
1884      else
1885        approx_type=6  ! rpa-type bootstrap
1886      end if
1887 
1888      option_test=MOD(Dtset%gwgamma,2)
1889      ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma
1890      ! 0 -> TESTPARTICLE, vertex in chi0 only
1891      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1892      ! --@WC
1893 
1894    case default
1895      MSG_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma)))
1896    end select
1897 
1898 !  Set plasma frequency
1899    my_plsmf = drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf = Dtset%ppmfrq
1900    Dtset%ppmfrq = my_plsmf
1901 
1902 !   if(dtset%ucrpa==0) then
1903    if (Dtset%gwgamma<3) then
1904      call mkdump_Er(Er,Vcp,Er%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,&
1905 &     approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,&
1906 &     nfftf_tot,ngfftf,comm)
1907    else
1908      call mkdump_Er(Er,Vcp,Er%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,&
1909 &     approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,&
1910 &     nfftf_tot,ngfftf,comm,fxc_ADA)
1911    end if
1912 !   end if
1913    if (allocated(kxcg))  then
1914      ABI_FREE(kxcg)
1915    end if
1916    if (allocated(fxc_ADA))  then
1917      ABI_FREE(fxc_ADA)
1918    end if
1919    if (allocated(ks_aepaw_rhor))  then
1920      ABI_FREE(ks_aepaw_rhor)
1921    end if
1922  end if
1923  !
1924  ! ================================================
1925  ! ==== Calculate plasmonpole model parameters ====
1926  ! ================================================
1927  ! TODO Maybe its better if we use mqmem as input variable
1928  use_aerhor=0
1929  ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden*use_aerhor))
1930 
1931  if (sigma_needs_ppm(Sigp)) then
1932    my_plsmf=drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf=Dtset%ppmfrq
1933    call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,Sigp%ppmodel,my_plsmf,Dtset%gw_invalid_freq)
1934 
1935    ! PPm%force_plsmf= force_ppmfrq  ! this line to change the plasme frequency in HL expression.
1936 
1937    if (Wfd%usepaw==1 .and. Ppm%userho==1) then
1938      ! * For PAW and ppmodel 2-3-4 we need the AE rho(G) without compensation charge.
1939      ! * It would be possible to calculate rho(G) using Paw_pwff, though. It should be faster but
1940      !    results will depend on the expression used for the matrix elements. This approach is safer.
1941      use_aerhor=1
1942      ABI_FREE(ks_aepaw_rhor)
1943      ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1944 
1945      ! Check if input density file is available, otherwise compute
1946      pawden_fname = strcat(Dtfil%filnam_ds(3), '_PAWDEN')
1947      call wrtout(std_out,sjoin('Checking for existence of:',pawden_fname))
1948      if (file_exists(pawden_fname)) then
1949        ! Read density from file
1950        ABI_DT_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1951 
1952        call read_rhor(pawden_fname, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1953        ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm)
1954 
1955        call hdr_free(hdr_rhor)
1956        call pawrhoij_free(tmp_pawrhoij)
1957        ABI_DT_FREE(tmp_pawrhoij)
1958      else
1959        ! Have to calculate PAW AW rhor from scratch
1960        ABI_MALLOC(qp_rhor_n_one,(pawfgr%nfft,Dtset%nspden))
1961        ABI_MALLOC(qp_rhor_nt_one,(pawfgr%nfft,Dtset%nspden))
1962 
1963 !      FIXME
1964        MSG_WARNING(" denfgr in sigma seems to produce wrong results")
1965        write(std_out,*)" input tilde ks_rhor integrates: ",SUM(ks_rhor(:,1))*Cryst%ucvol/nfftf
1966 
1967        call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,ks_nhat,Wfd%nspinor,&
1968        Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,KS_Pawrhoij,Pawtab,Dtset%prtvol, &
1969        ks_rhor,ks_aepaw_rhor,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
1970 
1971        ABI_FREE(qp_rhor_n_one)
1972        ABI_FREE(qp_rhor_nt_one)
1973      end if
1974 
1975      write(msg,'(a,f8.4)')' sigma: PAW AE density used for PPmodel integrates to: ',SUM(ks_aepaw_rhor(:,1))*Cryst%ucvol/nfftf
1976      call wrtout(std_out,msg,'PERS')
1977 
1978      if (Er%mqmem/=0) then
1979        ! Calculate ppmodel parameters for all q-points.
1980        call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_aepaw_rhor(:,1))
1981      end if
1982 
1983    else
1984      ! NC or PAW with PPmodel 1.
1985      if (Er%mqmem/=0) then
1986        ! Calculate ppmodel parameters for all q-points
1987        call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_rhor(:,1))
1988      end if
1989    end if ! PAW or NC PPm and/or needs density
1990  end if ! sigma_needs_ppm
1991 
1992  call timab(409,2,tsec) ! getW
1993 
1994  if (wfd_iam_master(Wfd)) then
1995    ! Write info on the run on ab_out, then open files to store final results.
1996    call ebands_report_gap(KS_BSt,header='KS Band Gaps',unit=ab_out)
1997    if(dtset%ucrpa==0) then
1998      call write_sigma_header(Sigp,Er,Cryst,Kmesh,Qmesh)
1999    end if
2000 
2001    if (open_file(Dtfil%fnameabo_gw,msg,unit=unt_gw,status='unknown',form='formatted') /=0) then
2002      MSG_ERROR(msg)
2003    end if
2004    write(unt_gw,*)Sigp%nkptgw,Sigp%nsppol
2005 
2006    if (open_file(Dtfil%fnameabo_gwdiag,msg,unit=unt_gwdiag,status='unknown',form='formatted') /= 0) then
2007      MSG_ERROR(msg)
2008    end if
2009    write(unt_gwdiag,*)Sigp%nkptgw,Sigp%nsppol
2010 
2011    if (open_file(Dtfil%fnameabo_sig,msg,unit=unt_sig,status='unknown',form='formatted') /= 0) then
2012      MSG_ERROR(msg)
2013    end if
2014    if (open_file(Dtfil%fnameabo_sgr,msg,unit=unt_sgr,status='unknown',form='formatted') /= 0) then
2015      MSG_ERROR(msg)
2016    end if
2017 
2018    if (mod10==SIG_GW_AC) then
2019      ! Sigma along the imaginary axis.
2020      if (open_file(Dtfil%fnameabo_sgm,msg,unit=unt_sgm,status='unknown',form='formatted') /= 0) then
2021        MSG_ERROR(msg)
2022      end if
2023    end if
2024  end if
2025 
2026 !=======================================================================
2027 !==== Calculate self-energy and output the results for each k-point ====
2028 !=======================================================================
2029 !* Here it would be possible to calculate the QP correction for the same k-point using a PPmodel
2030 !in the first iteration just to improve the initial guess and CD or AC in the second step. Useful
2031 !if the KS level is far from the expected QP results. Previously it was possible to do such
2032 !calculation by simply listing the same k-point twice in kptgw. Now this trick is not allowed anymore.
2033 !Everything, indeed, should be done in a clean and transparent way inside csigme.
2034 
2035  call wfd_print(Wfd,mode_paral='PERS')
2036 
2037  call wrtout(std_out,sigma_type_from_key(mod10),'COLL')
2038 
2039  if (gwcalctyp<10) then
2040    msg = " Perturbative Calculation"
2041    if (gwcalctyp==1) msg = " Newton Raphson method "
2042  else if (gwcalctyp<20) then
2043    msg = " Self-Consistent on Energies only"
2044  else
2045    msg = " Self-Consistent on Energies and Wavefunctions"
2046  end if
2047  if(Dtset%ucrpa==0) then
2048    call wrtout(std_out,msg,'COLL')
2049  end if
2050 !
2051 !=================================================
2052 !==== Calculate the matrix elements of Sigma =====
2053 !=================================================
2054 
2055  nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i
2056 
2057  ! min and max band indeces for GW corrections.
2058  ib1=Sigp%minbdgw; ib2=Sigp%maxbdgw
2059 
2060  !MG TODO: I don't like the fact that ib1 and ib2 are redefined here because this
2061  ! prevents me from refactoring the code. In particular I want to store the self-energy
2062  ! results inside the sigma_results datatypes hence one needs to know all the dimensions
2063  ! at the beginning of the execution (e.g. in setup_sigma) so that one can easily allocate the arrays in the type.
2064  if(Dtset%ucrpa>=1) then
2065    !Read the band
2066    if (open_file("forlb.ovlp",msg,newunit=temp_unt,form="formatted", status="unknown") /= 0) then
2067      MSG_ERROR(msg)
2068    end if
2069    rewind(temp_unt)
2070    read(temp_unt,*)
2071    read(temp_unt,*) msg, ib1, ib2
2072    close(temp_unt)
2073  end if
2074 
2075  ABI_MALLOC(sigcme,(nomega_sigc,ib1:ib2,ib1:ib2,Sigp%nkptgw,Sigp%nsppol*Sigp%nsig_ab))
2076  sigcme=czero
2077 
2078 ! if (.False. .and. psps%usepaw == 0 .and. wfd%nspinor == 1 .and. any(dtset%so_psp /= 0)) then
2079 !   call wrtout(std_out, "Computing SOC contribution with first-order perturbation theory")
2080 !   ABI_MALLOC(bks_mask, (wfd%mband, wfd%nkibz, wfd%nsppol))
2081 !   bks_mask = .False.
2082 !   do spin=1,wfd%nsppol
2083 !     do ikcalc=1,Sigp%nkptgw
2084 !       ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW
2085 !       ii=Sigp%minbnd(ikcalc, spin); jj=Sigp%maxbnd(ikcalc, spin)
2086 !       bks_mask(ii:jj, ik_ibz, spin) = .True.
2087 !     end do
2088 !   end do
2089 !
2090 !   call wfd_get_socpert(wfd, cryst, psps, pawtab, bks_mask, osoc_bks)
2091 !   ABI_FREE(bks_mask)
2092 !   ABI_FREE(osoc_bks)
2093 ! end if
2094 
2095 !==========================================================
2096 !==== Exchange part using the dense gwx_ngfft FFT mesh ====
2097 !==========================================================
2098 !TODO : load distribution is not optimal if band parallelism is used.
2099 !but this problem was also affecting the old implementation.
2100  call timab(421,1,tsec) ! calc_sigx_me
2101 
2102  ! This routine shows how the wavefunctions are distributed.
2103  !call wfd_show_bkstab(Wfd,unit=std_out)
2104 
2105  if(Dtset%ucrpa>=1) then
2106 !  !Calculation of the Rho_n for calc_ucrpa
2107    call wrtout(std_out,sjoin("begin of Ucrpa calc for a nb of kpoints of: ",itoa(Sigp%nkptgw)),"COLL")
2108 !   Wannier basis: rhot1_q_m will need an atom index in the near future.
2109    ic=0
2110    do  itypat=1,cryst%ntypat
2111      if(pawtab(itypat)%lpawu.ne.-1) then
2112        ndim=2*dtset%lpawu(itypat)+1
2113        ic=ic+1
2114        itypatcor=itypat
2115        lcor=dtset%lpawu(itypat)
2116      end if
2117    end do
2118    if(ic>1) then
2119      MSG_ERROR("number of correlated species is larger than one")
2120    end if
2121    ABI_ALLOCATE(rhot1_q_m,(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz))
2122    ABI_ALLOCATE(M1_q_m,(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz))
2123 
2124    M1_q_m=czero
2125    rhot1_q_m=czero
2126 
2127 !   do ikcalc=1,Sigp%nkptgw
2128 
2129 !   if(cryst%nsym==1) nkcalc=Kmesh%nbz
2130 !   if(cryst%nsym>1)  nkcalc=Kmesh%nibz
2131    nkcalc=Kmesh%nbz
2132    do ikcalc=1,nkcalc ! for the oscillator strengh, spins are identical without SOC
2133 !    if(Sigp%nkptgw/=Kmesh%nbz) then
2134 !      write(msg,'(6a)')ch10,&
2135 !&      ' nkptgw and nbz differs: this is not allowed to compute U in cRPA '
2136 !      MSG_ERROR(msg)
2137 !    endif
2138      write(std_out,*) "ikcalc",ikcalc,size(Sigp%kptgw2bz),nkcalc
2139 
2140      if(mod(ikcalc-1,nprocs)==Wfd%my_rank) then
2141        if(cryst%nsym==1) ik_ibz=Kmesh%tab(ikcalc) ! Index of the irred k-point for GW
2142        if(cryst%nsym>1) ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2143 
2144        call prep_calc_ucrpa(ik_ibz,ikcalc,itypatcor,ib1,ib2,Cryst,QP_bst,Sigp,Gsph_x,Vcp,&
2145 &       Kmesh,Qmesh,lcor,M1_q_m,Pawtab,Pawang,Paw_pwff,&
2146 &       Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,gwx_ngfft,&
2147 &       ngfftf,Dtset%prtvol,Dtset%pawcross,rhot1_q_m)
2148      end if
2149    end do
2150 
2151    call xmpi_sum(rhot1_q_m,Wfd%comm,ierr)
2152    call xmpi_sum(M1_q_m,Wfd%comm,ierr)
2153 !   if(Cryst%nsym==1) then
2154    M1_q_m=M1_q_m/Kmesh%nbz/Wfd%nsppol
2155    rhot1_q_m=rhot1_q_m/Kmesh%nbz/Wfd%nsppol
2156 !   endif
2157 !  Calculation of U in cRPA: need to treat with a different cutoff the
2158 !  bare coulomb and the screened coulomb interaction.
2159    call calc_ucrpa(itypatcor,cryst,Kmesh,lcor,M1_q_m,Qmesh,Er%npwe,sigp%npwx,&
2160 &   Cryst%nsym,rhot1_q_m,Sigp%nomegasr,Sigp%minomega_r,Sigp%maxomega_r,ib1,ib2,'Gsum',Cryst%ucvol,Wfd,Er%fname)
2161 
2162    ABI_DEALLOCATE(rhot1_q_m)
2163    ABI_DEALLOCATE(M1_q_m)
2164 
2165  else
2166    do ikcalc=1,Sigp%nkptgw
2167      ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2168      ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2169      ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2170 
2171      call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,QP_bst,Sigp,Sr,Gsph_x,Vcp,Kmesh,Qmesh,Ltg_k(ikcalc),&
2172 &     Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,&
2173 &     gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross)
2174    end do
2175 
2176    ! for the time being, do not remove this barrier!
2177    call xmpi_barrier(Wfd%comm)
2178 
2179    call timab(421,2,tsec) ! calc_sigx_me
2180 
2181    ! ==========================================================
2182    ! ==== Correlation part using the coarse gwc_ngfft mesh ====
2183    ! ==========================================================
2184    if (mod10/=SIG_HF) then
2185      do ikcalc=1,Sigp%nkptgw
2186        ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2187        ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2188        ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2189 
2190 !      sigcme_p => sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:) ! Causes annoying Fortran runtime warning on abiref.
2191        ABI_ALLOCATE(sigcme_k,(nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab))
2192        sigcme_k=czero
2193        if (any(mod10 == [SIG_SEX, SIG_COHSEX])) then
2194          ! Calculate static COHSEX or SEX using the coarse gwc_ngfft mesh.
2195          call cohsex_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_c,Vcp,Kmesh,Qmesh,&
2196 &         Ltg_k(ikcalc),Pawtab,Pawang,Paw_pwff,Psps,Wfd,QP_sym,&
2197 !&         gwc_ngfft,Dtset%iomode,Dtset%prtvol,sigcme_p)
2198 &         gwc_ngfft,Dtset%iomode,Dtset%prtvol,sigcme_k)
2199        else
2200          ! Compute correlated part using the coarse gwc_ngfft mesh.
2201          call calc_sigc_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Dtset,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,&
2202 &         Ltg_k(ikcalc),PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,&
2203 !&         gwc_ngfft,ngfftf,nfftf,ks_rhor,use_aerhor,ks_aepaw_rhor,sigcme_p)
2204 &         gwc_ngfft,ngfftf,nfftf,ks_rhor,use_aerhor,ks_aepaw_rhor,sigcme_k)
2205        end if
2206        sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:)=sigcme_k
2207        ABI_DEALLOCATE(sigcme_k)
2208 
2209      end do
2210    end if
2211 
2212    call xmpi_barrier(Wfd%comm)
2213 
2214    !  =====================================================
2215    !  ==== Solve Dyson equation storing results in Sr% ====
2216    !  =====================================================
2217    !  * Use perturbative approach or AC to find QP corrections.
2218    !  * If qp-GW, diagonalize also H0+Sigma in the KS basis set to get the
2219    !  new QP amplitudes and energies (Sr%eigvec_qp and Sr%en_qp_diago.
2220    !  TODO AC with spinor not implemented yet.
2221    !  TODO Diagonalization of Sigma+hhartre with AC is wrong.
2222    !
2223    call timab(425,1,tsec) ! solve_dyson
2224    do ikcalc=1,Sigp%nkptgw
2225      ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2226      ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indeces for GW corrections (for this k-point)
2227      ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2228 
2229 !    sigcme_p => sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:)   ! Causes annoying Fortran runtime warning on abiref.
2230      ABI_ALLOCATE(sigcme_k,(nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab))
2231      sigcme_k=sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:)
2232 
2233 !    call solve_dyson(ikcalc,ib1,ib2,nomega_sigc,Sigp,Kmesh,sigcme_p,QP_BSt%eig,Sr,Dtset%prtvol,Dtfil,Wfd%comm)
2234      call solve_dyson(ikcalc,ib1,ib2,nomega_sigc,Sigp,Kmesh,sigcme_k,QP_BSt%eig,Sr,Dtset%prtvol,Dtfil,Wfd%comm)
2235      ABI_DEALLOCATE(sigcme_k)
2236      !
2237      ! Calculate direct gap for each spin and print out final results.
2238      ! We use the valence index of the KS system because we still do not know
2239      ! all the QP corrections. Ideally one should use the QP valence index
2240      do spin=1,Sigp%nsppol
2241        if ( Sigp%maxbnd(ikcalc,spin) >= ks_vbik(ik_ibz,spin)+1 .and. &
2242 &       Sigp%minbnd(ikcalc,spin) <= ks_vbik(ik_ibz,spin)  ) then
2243          ks_iv=ks_vbik(ik_ibz,spin)
2244          Sr%egwgap (ik_ibz,spin)= Sr%egw(ks_iv+1,ik_ibz,spin) -  Sr%egw(ks_iv,ik_ibz,spin)
2245          Sr%degwgap(ik_ibz,spin)= Sr%degw(ks_iv+1,ik_ibz,spin) - Sr%degw(ks_iv,ik_ibz,spin)
2246        else
2247          ! The "gap" cannot be computed
2248          Sr%e0gap(ik_ibz,spin)=zero; Sr%egwgap (ik_ibz,spin)=zero; Sr%degwgap(ik_ibz,spin)=zero
2249        end if
2250      end do
2251 
2252      if (wfd_iam_master(Wfd)) call write_sigma_results(ikcalc,ik_ibz,Sigp,Sr,KS_BSt)
2253    end do !ikcalc
2254 
2255    call timab(425,2,tsec) ! solve_dyson
2256    call timab(426,1,tsec) ! finalize
2257 
2258    ! Update the energies in QP_BSt
2259    ! If QPSCGW, use diagonalized eigenvalues otherwise perturbative results.
2260    if (gwcalctyp>=10) then
2261      do ib=1,Sigp%nbnds
2262        QP_BSt%eig(ib,:,:)=Sr%en_qp_diago(ib,:,:)
2263      end do
2264    else
2265      QP_BSt%eig=Sr%egw
2266    end if
2267 
2268    ! ================================================================================
2269    ! ==== This part is done only if all k-points in the IBZ have been calculated ====
2270    ! ================================================================================
2271    if (Sigp%nkptgw==Kmesh%nibz) then
2272 
2273      ! Recalculate new occupations and Fermi level.
2274      call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=Dtset%prtvol)
2275      qp_vbik(:,:) = get_valence_idx(QP_BSt)
2276 
2277      write(msg,'(2a,3x,2(es16.6,a))')ch10,' New Fermi energy : ',QP_BSt%fermie,' Ha ,',QP_BSt%fermie*Ha_eV,' eV'
2278      call wrtout(std_out,msg,'COLL')
2279      call wrtout(ab_out,msg,'COLL')
2280 
2281      ! === If all k-points and all occupied bands are calculated, output EXX ===
2282      ! FIXME here be careful about the check on  ks_vbik in case of metals
2283      ! if (my_rank==master.and.Sigp%nkptgw==Kmesh%nibz.and.ALL(Sigp%minbnd(:)==1).and.ALL(Sigp%maxbnd(:)>=MAXVAL(nbv(:)))) then
2284      ! if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=MAXVAL(MAXVAL(ks_vbik(:,:),DIM=1))) ) then
2285      if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=ks_vbik) ) then  ! FIXME here the two arrays use a different indexing.
2286 
2287        ex_energy = sigma_get_exene(Sr,Kmesh,QP_BSt)
2288        write(msg,'(a,2(es16.6,a))')' New Exchange energy : ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV'
2289        call wrtout(std_out,msg,'COLL')
2290        call wrtout(ab_out,msg,'COLL')
2291      end if
2292 
2293      ! Report the QP gaps (Fundamental and Optical)
2294      call ebands_report_gap(QP_BSt,header='QP Band Gaps',unit=ab_out)
2295 
2296      ! Band structure interpolation from QP energies computed on the k-mesh.
2297      if (nint(dtset%einterp(1)) /= 0 .and. all(sigp%minbdgw == sigp%minbnd) .and. all(sigp%maxbdgw == sigp%maxbnd)) then
2298        call ebands_interpolate_kpath(QP_BSt, dtset, cryst, [sigp%minbdgw, sigp%maxbdgw], dtfil%filnam_ds(4), comm)
2299      end if
2300    end if ! Sigp%nkptgw==Kmesh%nibz
2301    !
2302    ! Write SCF data in case of self-consistent calculation ===
2303    ! Save Sr%en_qp_diago, Sr%eigvec_qp and m_lda_to_qp in the _QPS file.
2304    ! Note that in the first iteration qp_rhor contains KS rhor, then the mixed rhor.
2305    if (gwcalctyp>=10) then
2306      ! Calculate the new m_lda_to_qp
2307      call updt_m_lda_to_qp(Sigp,Kmesh,nscf,Sr,Sr%m_lda_to_qp)
2308 
2309      if (wfd_iam_master(Wfd)) then
2310        call wrqps(Dtfil%fnameabo_qps,Sigp,Cryst,Kmesh,Psps,Pawtab,QP_Pawrhoij,&
2311 &       Dtset%nspden,nscf,nfftf,ngfftf,Sr,QP_BSt,Sr%m_lda_to_qp,qp_rhor)
2312      end if
2313 
2314      ! Report the MAX variation for each kptgw and spin
2315      call wrtout(ab_out,ch10//' Convergence of QP corrections ','COLL')
2316      converged=.TRUE.
2317      do spin=1,Sigp%nsppol
2318        write(msg,'(a,i2,a)')' >>>>> For spin ',spin,' <<<<< '
2319        call wrtout(ab_out,msg,'COLL')
2320        do ikcalc=1,Sigp%nkptgw
2321          ib1 = Sigp%minbnd(ikcalc,spin); ib2 = Sigp%maxbnd(ikcalc,spin)
2322          ik_bz = Sigp%kptgw2bz(ikcalc); ik_ibz = Kmesh%tab(ik_bz)
2323          ii = imax_loc(ABS(Sr%degw(ib1:ib2,ik_ibz,spin)))
2324          max_degw = Sr%degw(ii,ik_ibz,spin)
2325          write(msg,('(a,i3,a,2f8.3,a,i3)'))'.  kptgw no:',ikcalc,'; Maximum DeltaE = (',max_degw*Ha_eV,') for band index:',ii
2326          call wrtout(ab_out,msg,'COLL')
2327          converged = (converged .and. ABS(max_degw) < Dtset%gw_toldfeig)
2328        end do
2329      end do
2330    end if
2331 
2332    ! ============================================
2333    ! ==== Save the GW results in NETCDF file ====
2334    ! ============================================
2335 #ifdef HAVE_NETCDF
2336    if (wfd_iam_master(Wfd)) then
2337      NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), '_SIGRES.nc'), xmpi_comm_self))
2338      NCF_CHECK(nctk_defnwrite_ivars(ncid, ["sigres_version"], [1]))
2339      NCF_CHECK(crystal_ncwrite(Cryst, ncid))
2340      NCF_CHECK(ebands_ncwrite(KS_Bst, ncid))
2341      NCF_CHECK(sigma_ncwrite(Sigp, Er, Sr, ncid))
2342      ! Add qp_rhor. Note that qp_rhor == ks_rhor if wavefunctions are not updated.
2343      !ncerr = nctk_write_datar("qp_rhor",path,ngfft,cplex,nfft,nspden,&
2344      !                          comm_fft,fftn3_distrib,ffti3_local,datar,action)
2345      NCF_CHECK(nf90_close(ncid))
2346    end if
2347 #endif
2348 
2349  end if ! ucrpa
2350 
2351  !----------------------------- END OF THE CALCULATION ------------------------
2352  !
2353  !=====================
2354  !==== Close Files ====
2355  !=====================
2356  if (wfd_iam_master(Wfd)) then
2357    close(unt_gw )
2358    close(unt_gwdiag)
2359    close(unt_sig)
2360    close(unt_sgr)
2361    if (mod10==SIG_GW_AC) close(unt_sgm)
2362  end if
2363  !
2364  !===============================================
2365  !==== Free arrays and local data structures ====
2366  !===============================================
2367  ABI_FREE(ks_vbik)
2368  ABI_FREE(qp_vbik)
2369  ABI_FREE(ph1d)
2370  ABI_FREE(ph1df)
2371  ABI_FREE(qp_rhor)
2372  ABI_FREE(ks_rhor)
2373  ABI_FREE(ks_rhog)
2374  ABI_FREE(qp_taur)
2375  ABI_FREE(ks_taur)
2376  ABI_FREE(ks_taug)
2377  ABI_FREE(ks_vhartr)
2378  ABI_FREE(ks_vtrial)
2379  ABI_FREE(vpsp)
2380  ABI_FREE(ks_vxc)
2381  ABI_FREE(xccc3d)
2382  ABI_FREE(grchempottn)
2383  ABI_FREE(grewtn)
2384  ABI_FREE(grvdw)
2385  ABI_FREE(sigcme)
2386 
2387  if (allocated(kxc)) then
2388    ABI_FREE(kxc)
2389  end if
2390  if (allocated(qp_vtrial)) then
2391    ABI_FREE(qp_vtrial)
2392  end if
2393 
2394  ABI_FREE(ks_nhat)
2395  ABI_FREE(ks_nhatgr)
2396  ABI_FREE(dijexc_core)
2397  ABI_FREE(ks_aepaw_rhor)
2398  call pawfgr_destroy(Pawfgr)
2399 
2400  ! Deallocation for PAW.
2401  if (Dtset%usepaw==1) then
2402    call pawrhoij_free(KS_Pawrhoij)
2403    ABI_DT_FREE(KS_Pawrhoij)
2404    call pawfgrtab_free(Pawfgrtab)
2405    call paw_ij_free(KS_paw_ij)
2406    call paw_an_free(KS_paw_an)
2407    call pawpwff_free(Paw_pwff)
2408    if (gwcalctyp>=10) then
2409      call pawrhoij_free(QP_pawrhoij)
2410      call paw_ij_free(QP_paw_ij)
2411      call paw_an_free(QP_paw_an)
2412    end if
2413    if (Dtset%pawcross==1) then
2414      call paw_pwaves_lmn_free(Paw_onsite)
2415      Wfdf%bks_comm = xmpi_comm_null
2416      call wfd_free(Wfdf)
2417    end if
2418  end if
2419  ABI_DT_FREE(Paw_onsite)
2420  ABI_DT_FREE(Pawfgrtab)
2421  ABI_DT_FREE(Paw_pwff)
2422  ABI_DT_FREE(KS_paw_ij)
2423  ABI_DT_FREE(KS_paw_an)
2424  if (gwcalctyp>=10) then
2425    ABI_DT_FREE(QP_pawrhoij)
2426    ABI_DT_FREE(QP_paw_an)
2427    ABI_DT_FREE(QP_paw_ij)
2428  end if
2429 
2430  call wfd_free(Wfd)
2431  call destroy_mpi_enreg(MPI_enreg_seq)
2432  call littlegroup_free(Ltg_k)
2433  ABI_DT_FREE(Ltg_k)
2434  call kmesh_free(Kmesh)
2435  call kmesh_free(Qmesh)
2436  call gsph_free(Gsph_Max)
2437  call gsph_free(Gsph_x)
2438  call gsph_free(Gsph_c)
2439  call vcoul_free(Vcp)
2440  call crystal_free(Cryst)
2441  call sigma_free(Sr)
2442  call em1results_free(Er)
2443  call ppm_free(PPm)
2444  call hdr_free(Hdr_sigma)
2445  call hdr_free(Hdr_wfk)
2446  call ebands_free(KS_BSt)
2447  call ebands_free(QP_BSt)
2448  call melements_free(KS_me)
2449  call esymm_free(KS_sym)
2450  ABI_DT_FREE(KS_sym)
2451 
2452  if (Sigp%symsigma==1.and.gwcalctyp>=20) then
2453    call esymm_free(QP_sym)
2454    ABI_DT_FREE(QP_sym)
2455  end if
2456 
2457  call sigparams_free(Sigp)
2458 
2459  call timab(426,2,tsec) ! finalize
2460  call timab(401,2,tsec)
2461 
2462  DBG_EXIT('COLL')
2463 
2464 end subroutine sigma