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

SOURCE

 16 #if defined HAVE_CONFIG_H
 17 #include "config.h"
 18 #endif
 19 
 20 #include "abi_common.h"
 21 
 22 module m_sigma_driver
 23 
 24  use defs_basis
 25  use m_gwdefs
 26  use defs_wvltypes
 27  use m_dtset
 28  use m_xmpi
 29  use m_xomp
 30  use m_errors
 31  use m_abicore
 32  use m_ab7_mixing
 33  use m_nctk
 34  use m_kxc
 35  use m_distribfft
 36  use netcdf
 37  use m_hdr
 38  use libxc_functionals
 39  use m_dtfil
 40  use m_crystal
 41  use m_cgtools
 42 
 43  use defs_datatypes,  only : pseudopotential_type, ebands_t
 44  use defs_abitypes,   only : MPI_type
 45  use m_time,          only : timab
 46  use m_numeric_tools, only : imax_loc
 47  use m_fstrings,      only : strcat, sjoin, itoa, basename, ktoa, ltoa
 48  use m_hide_blas,     only : xdotc
 49  use m_io_tools,      only : open_file, file_exists, iomode_from_fname
 50  use m_mpinfo,        only : destroy_mpi_enreg, initmpi_seq
 51  use m_geometry,      only : normv, mkrdim, metric
 52  use m_fftcore,       only : print_ngfft
 53  use m_fft_mesh,      only : get_gftt, setmesh
 54  use m_fft,           only : fourdp
 55  use m_ioarr,         only : fftdatar_write, read_rhor
 56  use m_ebands,        only : ebands_update_occ, ebands_copy, ebands_report_gap, ebands_get_valence_idx, ebands_get_bandenergy,&
 57                              ebands_free, ebands_init, ebands_ncwrite, ebands_interpolate_kpath, get_eneocc_vect, &
 58                              ebands_enclose_degbands, ebands_get_gaps, gaps_t
 59  use m_energies,      only : energies_type, energies_init
 60  use m_bz_mesh,       only : kmesh_t, littlegroup_t, littlegroup_free, isamek, get_ng0sh, find_qmesh
 61  use m_gsphere,       only : gsphere_t, merge_and_sort_kg, setshells
 62  use m_kg,            only : getph, getcut
 63  use m_xcdata,        only : get_xclevel
 64  use m_wfd,           only : wfd_init, wfdgw_t, wfdgw_copy, test_charge, wave_t
 65  use m_vcoul,         only : vcoul_t
 66  use m_qparticles,    only : wrqps, rdqps, rdgw, show_QP, updt_m_ks_to_qp
 67  use m_screening,     only : mkdump_er, em1results_free, epsilonm1_results, init_er_from_file
 68  use m_ppmodel,       only : ppm_init, ppm_free, setup_ppmodel, getem1_from_PPm, ppmodel_t
 69  use m_sigma,         only : sigma_init, sigma_free, sigma_ncwrite, sigma_t, sigma_get_exene, &
 70                              mels_get_haene, mels_get_kiene, sigma_get_excene, write_sigma_header, write_sigma_results
 71  use m_dyson_solver,  only : solve_dyson
 72  use m_esymm,         only : esymm_t, esymm_free, esymm_failed
 73  use m_melemts,       only : melflags_t, melements_t
 74  use m_pawang,        only : pawang_type
 75  use m_pawrad,        only : pawrad_type
 76  use m_pawtab,        only : pawtab_type, pawtab_print, pawtab_get_lsize
 77  use m_paw_an,        only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
 78  use m_paw_ij,        only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_print
 79  use m_pawfgrtab,     only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print
 80  use m_pawrhoij,      only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, &
 81                              pawrhoij_inquire_dim, pawrhoij_symrhoij, pawrhoij_unpack
 82  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_free, paw_overlap
 83  use m_pawdij,        only : pawdij, symdij_all
 84  use m_pawfgr,        only : pawfgr_type, pawfgr_init, pawfgr_destroy
 85  use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
 86  use m_pawpwij,       only : pawpwff_t, pawpwff_init, pawpwff_free
 87  use m_paw_slater,    only : paw_mkdijexc_core, paw_dijhf
 88  use m_paw_dmft,      only : paw_dmft_type
 89  use m_paw_sphharm,   only : setsym_ylm
 90  use m_paw_mkrho,     only : denfgr
 91  use m_paw_nhat,      only : nhatgrid, pawmknhat
 92  use m_paw_tools,     only : chkpawovlp, pawprt
 93  use m_paw_denpot,    only : pawdenpot
 94  use m_paw_init,      only : pawinit, paw_gencond
 95  use m_classify_bands,only : classify_bands
 96  use m_wfk,           only : wfk_read_eigenvalues
 97  use m_io_kss,        only : make_gvec_kss
 98  use m_vhxc_me,       only : calc_vhxc_me
 99  use m_cohsex,        only : cohsex_me
100  use m_sigx,          only : calc_sigx_me
101  use m_sigc,          only : calc_sigc_me
102  use m_setvtr,        only : setvtr
103  use m_mkrho,         only : prtrhomxmn
104  use m_pspini,        only : pspini
105  use m_calc_ucrpa,    only : calc_ucrpa
106  use m_prep_calc_ucrpa,only : prep_calc_ucrpa
107  use m_paw_correlations,only : pawpuxinit
108  use m_spacepar,      only : hartre
109  use m_gwrdm,         only : calc_rdmx,calc_rdmc,natoccs,update_hdr_bst,print_tot_occ,get_chkprdm,&
110                              print_chkprdm,change_matrix,print_total_energy,print_band_energies,quadrature_sigma_cw
111  use m_plowannier,    only : operwan_realspace_type,plowannier_type,init_plowannier,get_plowannier,&
112                              fullbz_plowannier,init_operwan_realspace,reduce_operwan_realspace,&
113                              destroy_operwan_realspace,destroy_plowannier,zero_operwan_realspace
114 
115  implicit none
116 
117  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<wfdgw_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_ebands<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.

SOURCE

4139 subroutine paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,qp_ebands, &
4140                       Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij, &
4141                       QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim, &
4142                       qp_nhatgr,qp_compch_sph,qp_compch_fft)
4143 
4144 !Arguments ------------------------------------
4145 !scalars
4146  integer,intent(in) :: nfftf,nscf,nhatgrdim
4147  real(dp),intent(out) :: qp_compch_fft,qp_compch_sph
4148  type(kmesh_t),intent(in) :: Kmesh
4149  type(crystal_t),intent(in) :: Cryst
4150  type(Dataset_type),intent(in) :: Dtset
4151  type(Pseudopotential_type),intent(in) :: Psps
4152  type(Pawang_type),intent(in) :: Pawang
4153  type(ebands_t),intent(in) :: qp_ebands
4154  type(Energies_type),intent(inout) :: QP_energies
4155  type(wfdgw_t),intent(inout) :: Wfd
4156 !arrays
4157  integer,intent(in) :: ngfftf(18)
4158  real(dp),intent(out) :: qp_nhat(nfftf,Dtset%nspden)
4159  real(dp),intent(out) :: qp_nhatgr(nfftf,Dtset%nspden,3*nhatgrdim)
4160  type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom)
4161  type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw)
4162  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
4163  type(Pawrhoij_type),intent(inout) :: prev_Pawrhoij(Cryst%natom)
4164  type(Pawrhoij_type),intent(out) :: QP_pawrhoij(Cryst%natom)
4165  type(Paw_ij_type),intent(out) :: QP_paw_ij(Cryst%natom)
4166  type(Paw_an_type),intent(inout) :: QP_paw_an(Cryst%natom)
4167 
4168 !Local variables ------------------------------
4169 !scalars
4170  integer,parameter :: ipert0 = 0, idir0 = 0, optrhoij1 = 1
4171  integer :: choice,cplex,cplex_rhoij,has_dijU,has_dijso,iat,ider
4172  integer :: izero,nkxc1,nspden_rhoij,nzlmopt, option,usexcnhat
4173  character(len=500) :: msg
4174 !arrays
4175  real(dp) :: k0(3)
4176 
4177 !************************************************************************
4178 
4179  ABI_UNUSED(Kmesh%nibz)
4180  !
4181  !  * 0 if Vloc in atomic data is Vbare    (Blochl s formulation)
4182  !  * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation)
4183  usexcnhat=MAXVAL(Pawtab(:)%usexcnhat)
4184  !
4185  ! Calculate new rhoij_qp from updated Cprj_ibz, note use_rhoij_=1.
4186  call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
4187                            nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc)
4188  call pawrhoij_alloc(QP_pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,&
4189                      pawtab=Pawtab,use_rhoij_=1,use_rhoijres=1)
4190 
4191  ! FIXME kptop should be passed via Kmesh, in GW time reversal is always assumed.
4192  call wfd%pawrhoij(Cryst,qp_ebands,Dtset%kptopt,QP_pawrhoij,Dtset%pawprtvol)
4193  !
4194  ! Symmetrize QP $\rho_{ij}$.
4195  choice=1
4196  call pawrhoij_symrhoij(QP_pawrhoij,QP_pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,&
4197                Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij1,Pawang,Dtset%pawprtvol,&
4198                Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
4199 
4200  ! ======================
4201  ! ==== Make QP nhat ====
4202  ! ======================
4203  cplex=1; ider=2*nhatgrdim; izero=0; k0(:)=zero
4204 
4205  call pawmknhat(qp_compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,&
4206    Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
4207    Pawfgrtab,qp_nhatgr,qp_nhat,QP_Pawrhoij,QP_Pawrhoij,Pawtab,k0,Cryst%rprimd,&
4208    Cryst%ucvol,Dtset%usewvl,Cryst%xred)
4209 
4210  ! Allocate quantities related to the PAW spheres for the QP Hamiltonian.
4211  ! TODO call paw_ij_init in scfcv and respfn, fix small issues
4212  has_dijso=Dtset%pawspnorb; has_dijU=merge(0,1,Dtset%usepawu==0)
4213 
4214  call paw_ij_nullify(QP_paw_ij)
4215  call paw_ij_init(QP_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,&
4216    Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
4217    has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,&
4218    has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1)
4219 
4220  call paw_an_nullify(QP_paw_an); nkxc1=0 ! No kernel
4221  call paw_an_init(QP_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
4222       cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1)
4223 
4224  ! =====================================================
4225  ! ==== Optional mixing of the PAW onsite densities ====
4226  ! =====================================================
4227  if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then
4228    write(msg,'(a,f6.3)')' sigma: mixing on-site QP rho_ij densities using rhoqpmix= ',Dtset%rhoqpmix
4229    call wrtout(std_out, msg)
4230    ! qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor)
4231 
4232    call pawrhoij_unpack(QP_Pawrhoij)   ! Unpack new QP %rhoijp
4233    call pawrhoij_unpack(prev_Pawrhoij) ! Unpack previous QP %rhoijp
4234 
4235    do iat=1,Cryst%natom
4236      QP_pawrhoij(iat)%rhoij_ = prev_Pawrhoij(iat)%rhoij_ &
4237       + Dtset%rhoqpmix * (QP_pawrhoij(iat)%rhoij_ - prev_pawrhoij(iat)%rhoij_)
4238 
4239      prev_pawrhoij(iat)%use_rhoij_=0
4240      ABI_FREE(prev_pawrhoij(iat)%rhoij_)
4241    end do
4242 
4243    ! Re-Symmetrize mixed QP $\rho_{ij}$.
4244    choice=1
4245    call pawrhoij_symrhoij(QP_pawrhoij,QP_pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,&
4246                 Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij1,Pawang,Dtset%pawprtvol,&
4247                 Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
4248  end if
4249 
4250  do iat=1,Cryst%natom
4251    QP_pawrhoij(iat)%use_rhoij_=0
4252    ABI_FREE(QP_pawrhoij(iat)%rhoij_)
4253  end do
4254  !
4255  ! =================================================================================
4256  ! ==== Evaluate on-site energies, potentials, densities using (mixed) QP rhoij ====
4257  ! =================================================================================
4258  ! * Initialize also "lmselect" (index of non-zero LM-moments of densities).
4259 
4260  nzlmopt=-1; option=0; qp_compch_sph=greatest_real
4261 
4262  call pawdenpot(qp_compch_sph,QP_energies%e_paw,QP_energies%e_pawdc,&
4263    ipert0,Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
4264    Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,QP_paw_an,QP_paw_an,&
4265    QP_paw_ij,Pawang,Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,&
4266    Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp)
4267 
4268 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
 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_ebands<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.

SOURCE

2968 subroutine setup_sigma(codvsn,wfk_fname,acell,rprim,Dtset,Dtfil,Psps,Pawtab,&
2969 & gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,ks_ebands,Gsph_Max,Gsph_x,Gsph_c,Vcp,Er,Sigp,comm)
2970 
2971 !Arguments ------------------------------------
2972 !scalars
2973  integer,intent(in) :: comm
2974  character(len=8),intent(in) :: codvsn
2975  character(len=*),intent(in) :: wfk_fname
2976  type(Datafiles_type),intent(in) :: Dtfil
2977  type(Dataset_type),intent(inout) :: Dtset
2978  type(Pseudopotential_type),intent(in) :: Psps
2979  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw)
2980  type(sigparams_t),intent(out) :: Sigp
2981  type(Epsilonm1_results),intent(out) :: Er
2982  type(ebands_t),intent(out) :: ks_ebands
2983  type(kmesh_t),intent(out) :: Kmesh,Qmesh
2984  type(crystal_t),intent(out) :: Cryst
2985  type(gsphere_t),intent(out) :: Gsph_Max,Gsph_x,Gsph_c
2986  type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out
2987  type(vcoul_t),intent(out) :: Vcp
2988 !arrays
2989  integer,intent(out) :: gwc_ngfft(18),gwx_ngfft(18)
2990  real(dp),intent(in) :: acell(3),rprim(3,3)
2991 
2992 !Local variables-------------------------------
2993 !scalars
2994  integer,parameter :: pertcase0=0,master=0
2995  integer :: bantot,enforce_sym,gwcalctyp,ib,ibtot,icutcoul_eff,ii,ikcalc,ikibz,io,isppol,itypat,jj,method
2996  integer :: mod10,mqmem,mband,ng_kss,nsheps,ikcalc2bz,ierr,gap_err,ng, nsppol
2997  integer :: gwc_nfftot,gwx_nfftot,nqlwl,test_npwkss,my_rank,nprocs,ik,nk_found,ifo,timrev,usefock_ixc
2998  integer :: iqbz,isym,iq_ibz,itim,ic,pinv,ig1,ng_sigx,spin,gw_qprange,ivcoul_init,nvcoul_init,xclevel_ixc
2999  real(dp),parameter :: OMEGASIMIN=0.01d0,tol_enediff=0.001_dp*eV_Ha
3000  real(dp) :: domegas,domegasi,ucvol,ucvol_local,rcut
3001  logical,parameter :: linear_imag_mesh=.TRUE.
3002  logical :: ltest,remove_inv,changed,found
3003  character(len=500) :: msg
3004  character(len=fnlen) :: fname,fcore,string
3005  type(wvl_internal_type) :: wvl
3006  type(gaps_t) :: gaps
3007 !arrays
3008  integer :: ng0sh_opt(3),G0(3),q_umklp(3),kpos(6)
3009  integer,allocatable :: npwarr(:),val_indeces(:,:)
3010  integer,pointer :: gvec_kss(:,:),gsphere_sigx_p(:,:)
3011  integer,pointer :: test_gvec_kss(:,:)
3012  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3),q_bz(3),gamma_point(3,1)
3013  real(dp),pointer :: energies_p(:,:,:)
3014  real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:)
3015  type(Pawrhoij_type),allocatable :: Pawrhoij(:)
3016  type(vcoul_t) :: Vcp_ks
3017 
3018 ! *************************************************************************
3019 
3020  DBG_ENTER('COLL')
3021 
3022  ! Check for calculations that are not implemented
3023  nsppol = dtset%nsppol
3024  ltest = ALL(Dtset%nband(1:Dtset%nkpt*nsppol) == Dtset%nband(1))
3025  ABI_CHECK(ltest,'Dtset%nband(:) must be constant')
3026 
3027  my_rank = xmpi_comm_rank(comm); nprocs  = xmpi_comm_size(comm)
3028 
3029  ! Basic parameters
3030  Sigp%ppmodel    = Dtset%ppmodel
3031  Sigp%gwcalctyp  = Dtset%gwcalctyp
3032  Sigp%nbnds      = Dtset%nband(1)
3033  Sigp%symsigma   = Dtset%symsigma
3034  Sigp%zcut       = Dtset%zcut
3035  Sigp%mbpt_sciss = Dtset%mbpt_sciss
3036 
3037  timrev=  2 ! This information is not reported in the header
3038             ! 1 => do not use time-reversal symmetry
3039             ! 2 => take advantage of time-reversal symmetry
3040  if (any(dtset%kptopt == [3, 4])) timrev = 1
3041 
3042  ! === For HF, SEX or COHSEX use Hybertsen-Louie PPM (only $\omega=0$) ===
3043  ! * Use fake screening for HF.
3044  ! FIXME Why, we should not redefine Sigp%ppmodel
3045  gwcalctyp=Sigp%gwcalctyp
3046  mod10 =MOD(Sigp%gwcalctyp,10)
3047  if (mod10==5.or.mod10==6.or.mod10==7) Sigp%ppmodel=2
3048  if (mod10<5.and.MOD(Sigp%gwcalctyp,1)/=1) then ! * One shot GW (PPM or contour deformation).
3049    if (Dtset%nomegasrd==1) then ! avoid division by zero!
3050      Sigp%nomegasrd  =1
3051      Sigp%maxomega4sd=zero
3052      Sigp%deltae     =zero
3053    else
3054      Sigp%nomegasrd   = Dtset%nomegasrd
3055      Sigp%maxomega4sd = Dtset%omegasrdmax
3056      Sigp%deltae     = (2*Sigp%maxomega4sd)/(Sigp%nomegasrd-1)
3057    endif
3058  else
3059    ! For AC no need to evaluate derivative by finite differences.
3060    Sigp%nomegasrd  =1
3061    Sigp%maxomega4sd=zero
3062    Sigp%deltae     =zero
3063  end if
3064 
3065  ! For analytic continuation define the number of imaginary frequencies for Sigma
3066  ! Tests show than more than 12 freqs in the Pade approximant worsen the results!
3067  Sigp%nomegasi=0
3068 
3069  if (mod10==1) then
3070    Sigp%nomegasi  =Dtset%nomegasi
3071    Sigp%omegasimax=Dtset%omegasimax
3072    Sigp%omegasimin=OMEGASIMIN
3073    write(msg,'(4a,i3,2(2a,f8.3),a)')ch10,&
3074     ' Parameters for analytic continuation : ',ch10,&
3075     '  number of imaginary frequencies for sigma =  ',Sigp%nomegasi,ch10,&
3076     '  min frequency for sigma on imag axis [eV] =  ',Sigp%omegasimin*Ha_eV,ch10,&
3077     '  max frequency for sigma on imag axis [eV] =  ',Sigp%omegasimax*Ha_eV,ch10
3078    call wrtout(std_out, msg)
3079 
3080    !TODO this should not be done here but in init_sigma_t
3081    ABI_MALLOC(Sigp%omegasi,(Sigp%nomegasi))
3082 
3083    if (linear_imag_mesh) then
3084      ! Linear mesh along the imaginary axis.
3085      domegasi=Sigp%omegasimax/(Sigp%nomegasi-1)
3086      do io=1,Sigp%nomegasi
3087        Sigp%omegasi(io)=CMPLX(zero,(io-1)*domegasi)
3088      end do
3089    else
3090      ! Logarithmic mesh along the imaginary axis.
3091      ABI_ERROR("AC + log mesh not implemented")
3092      !domegasi=(Sigp%omegasimax/Sigp%omegasimin)**(one/(Sigp%nomegasi-1))
3093      !Sigp%omegasi(1)=czero; ldi=domegasi
3094      !do io=2,Sigp%nomegasi
3095      ! omega(io)=CMPLX(zero,ldi*Sigp%omegasimin)
3096      ! Sigp%omegasi(io)=ldi*domegasi
3097      !end do
3098    end if
3099 
3100    ! MRM: do not print for 1-RDM correction
3101    if(Sigp%gwcalctyp/=21) then
3102     write(msg,'(4a)')ch10,&
3103      ' setup_sigma : calculating Sigma(iw)',&
3104      ' at imaginary frequencies [eV] (Fermi Level set to 0) ',ch10
3105     call wrtout([std_out, ab_out], msg)
3106     do io=1,Sigp%nomegasi
3107       write(msg,'(2(f10.3,2x))')Sigp%omegasi(io)*Ha_eV
3108       call wrtout([std_out, ab_out], msg)
3109     end do
3110    endif
3111 
3112    ltest=(Sigp%omegasimax>0.1d-4.and.Sigp%nomegasi>0)
3113    ABI_CHECK(ltest,'Wrong value of omegasimax or nomegasi')
3114    if (Sigp%gwcalctyp/=1) then ! only one shot GW is allowed for AC.
3115      !ABI_ERROR("SC-GW with analytic continuation is not coded") ! MRM: let's allow it
3116    end if
3117  end if
3118 
3119  if (Sigp%symsigma/=0.and.gwcalctyp>=20) then
3120    ABI_WARNING("SC-GW with symmetries is still under development. Use at your own risk!")
3121    ABI_ERROR("SC-GW requires symsigma == 0 in input. New default in Abinit9 is symsigma 1!")
3122  end if
3123 
3124  ! Setup parameters for Spectral function.
3125  if (Dtset%gw_customnfreqsp/=0) then
3126    Sigp%nomegasr = Dtset%gw_customnfreqsp
3127    ABI_WARNING('Custom grid for spectral function specified. Assuming experienced user.')
3128    if (Dtset%gw_customnfreqsp/=0) then
3129      Dtset%nfreqsp = Dtset%gw_customnfreqsp
3130      ABI_WARNING('nfreqsp has been set to the same number as gw_customnfreqsp')
3131    end if
3132  else
3133    Sigp%nomegasr  =Dtset%nfreqsp
3134    Sigp%minomega_r=Dtset%freqspmin
3135    Sigp%maxomega_r=Dtset%freqspmax
3136  end if
3137 
3138  if (Sigp%nomegasr>0) then
3139    if (Dtset%gw_customnfreqsp==0) then
3140      ! Check
3141      if (Sigp%minomega_r >= Sigp%maxomega_r) then
3142        ABI_ERROR('freqspmin must be smaller than freqspmax!')
3143      end if
3144      if(Sigp%nomegasr==1) then
3145       domegas=0.d0
3146      else
3147       domegas=(Sigp%maxomega_r-Sigp%minomega_r)/(Sigp%nomegasr-1)
3148      endif
3149      !TODO this should be moved to Sr% and done in init_sigma_t
3150      ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr))
3151      do io=1,Sigp%nomegasr
3152        Sigp%omega_r(io) = CMPLX(Sigp%minomega_r + domegas*(io-1),zero)
3153      end do
3154      write(msg,'(4a,i8,3(2a,f8.3),a)')ch10,&
3155        ' Parameters for the calculation of the spectral function : ',ch10,&
3156        '  Number of points    = ',Sigp%nomegasr,ch10,&
3157        '  Min frequency  [eV] = ',Sigp%minomega_r*Ha_eV,ch10,&
3158        '  Max frequency  [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,&
3159        '  Frequency step [eV] = ',domegas*Ha_eV,ch10
3160      call wrtout(std_out, msg)
3161    else
3162      Sigp%minomega_r = MINVAL(Dtset%gw_freqsp(:))
3163      Sigp%maxomega_r = MAXVAL(Dtset%gw_freqsp(:))
3164      !TODO this should be moved to Sr% and done in init_sigma_t
3165      ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr))
3166      do io=1,Sigp%nomegasr
3167        Sigp%omega_r(io) = CMPLX(Dtset%gw_freqsp(io),zero)
3168      end do
3169      write(msg,'(4a,i8,2(2a,f8.3),3a)')ch10,&
3170        ' Parameters for the calculation of the spectral function : ',ch10,&
3171        '  Number of points    = ',Sigp%nomegasr,ch10,&
3172        '  Min frequency  [eV] = ',Sigp%minomega_r*Ha_eV,ch10,&
3173        '  Max frequency  [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,&
3174        '  A custom set of frequencies is used! See the input file for values.',ch10
3175      call wrtout(std_out, msg)
3176    end if
3177  else
3178    !In indefo all these quantities are set to zero
3179    !Sigp%nomegasr=1
3180    !allocate(Sigp%omega_r(Sigp%nomegasr))
3181    !Sigp%omega_r(1)=0
3182  end if
3183 
3184  ! Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume
3185  call mkrdim(acell,rprim,rprimd)
3186  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
3187  ucvol_local = ucvol
3188 
3189  Sigp%npwwfn = Dtset%npwwfn
3190  Sigp%npwx   = Dtset%npwsigx
3191 
3192  ! Read parameters of the WFK, verifify them and retrieve all G-vectors.
3193  call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm)
3194  mband = MAXVAL(Hdr_wfk%nband)
3195 
3196  remove_inv = .FALSE.
3197  call hdr_wfk%vs_dtset(dtset)
3198 
3199  test_npwkss = 0
3200  call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,&
3201                     gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr)
3202  ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss")
3203 
3204  ABI_MALLOC(gvec_kss,(3, test_npwkss))
3205  gvec_kss = test_gvec_kss
3206  ng_kss = test_npwkss
3207 
3208  ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2))
3209  ierr = 0
3210  do ig1=1,ng
3211    if (ANY(gvec_kss(:,ig1)/=test_gvec_kss(:,ig1))) then
3212      ierr=ierr+1
3213      write(std_out,*)" gvec_kss ",ig1,"/",ng,gvec_kss(:,ig1),test_gvec_kss(:,ig1)
3214    end if
3215  end do
3216  ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss")
3217 
3218  ABI_FREE(test_gvec_kss)
3219 
3220  ! Get important dimensions from the WFK header
3221  Sigp%nsppol = Hdr_wfk%nsppol
3222  Sigp%nspinor = Hdr_wfk%nspinor
3223  Sigp%nsig_ab = Hdr_wfk%nspinor**2  ! TODO Is it useful calculating only diagonal terms?
3224 
3225  if (Sigp%nbnds > mband) then
3226    Sigp%nbnds     = mband
3227    Dtset%nband(:) = mband
3228    Dtset%mband    = MAXVAL(Dtset%nband)
3229    write(msg,'(3a,i4,a)')&
3230     'Number of bands found less then required',ch10,&
3231     'calculation will proceed with nbnds = ',mband,ch10
3232    ABI_WARNING(msg)
3233  end if
3234 
3235  ! Check input
3236  if (Sigp%ppmodel==3.or.Sigp%ppmodel==4) then
3237    if (gwcalctyp>=10) then
3238      write(msg,'(a,i3,a)')' The ppmodel chosen and gwcalctyp ',Dtset%gwcalctyp,' are not compatible. '
3239      ABI_ERROR(msg)
3240    end if
3241    if (Sigp%nspinor==2) then
3242      write(msg,'(a,i3,a)')' The ppmodel chosen and nspinor ',Sigp%nspinor,' are not compatible. '
3243      ABI_ERROR(msg)
3244    end if
3245  end if
3246 
3247  ! Create crystal_t data type
3248  cryst = Hdr_wfk%get_crystal(gw_timrev=timrev, remove_inv=remove_inv)
3249  call cryst%print()
3250 
3251  if (Sigp%npwwfn > ng_kss) then ! cannot use more G"s for the wfs than those stored on file
3252    Sigp%npwwfn  =ng_kss
3253    Dtset%npwwfn =ng_kss
3254    write(msg,'(2a,(a,i8,a))')&
3255     'Number of G-vectors for WFS found in the KSS file is less than required',ch10,&
3256     'calculation will proceed with npwwfn  = ',Sigp%npwwfn,ch10
3257    ABI_WARNING(msg)
3258  end if
3259 
3260  if (Sigp%npwx>ng_kss) then
3261    ! Have to recalcuate the (large) sphere for Sigma_x.
3262    pinv=1; if (remove_inv.and.Cryst%timrev==2) pinv=-1
3263    gamma_point(:,1) = (/zero,zero,zero/); nullify(gsphere_sigx_p)
3264 
3265    call merge_and_sort_kg(1,gamma_point,Dtset%ecutsigx,Cryst%nsym,pinv,Cryst%symrel,&
3266      Cryst%gprimd,gsphere_sigx_p,Dtset%prtvol)
3267 
3268    ng_sigx=SIZE(gsphere_sigx_p,DIM=2)
3269    Sigp%npwx     = ng_sigx
3270    Dtset%npwsigx = ng_sigx
3271 
3272    write(msg,'(2a,(a,i8,a))')&
3273      'Number of G-vectors for Sigma_x found in the KSS file is less than required',ch10,&
3274      'calculation will proceed with npwsigx = ',Sigp%npwx,ch10
3275    ABI_WARNING(msg)
3276 
3277    ltest = (Sigp%npwx >= ng_kss)
3278    ABI_CHECK(ltest,"Sigp%npwx<ng_kss!")
3279 
3280    ! * Fill gvec_kss with larger sphere.
3281    ABI_FREE(gvec_kss)
3282    ABI_MALLOC(gvec_kss,(3,Sigp%npwx))
3283    gvec_kss = gsphere_sigx_p
3284    ABI_FREE(gsphere_sigx_p)
3285  end if
3286 
3287  ! Set up of the k-points and tables in the whole BZ ===
3288  ! TODO Recheck symmorphy and inversion
3289  call Kmesh%init(Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.FALSE.)
3290  !call Kmesh%init(Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.TRUE.)
3291 
3292  ! Some required information are not filled up inside kmesh_init
3293  ! So doing it here, even though it is not clean
3294  Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:)
3295  Kmesh%nshift        =Dtset%nshiftk
3296  ABI_MALLOC(Kmesh%shift,(3,Kmesh%nshift))
3297  Kmesh%shift(:,:)    =Dtset%shiftk(:,1:Dtset%nshiftk)
3298 
3299  call Kmesh%print("K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL")
3300  call Kmesh%print("K-mesh for the wavefunctions",ab_out, 0,           "COLL")
3301 
3302  ! === Initialize the band structure datatype ===
3303  ! * Copy WFK energies and occupations up to Sigp%nbnds==Dtset%nband(:)
3304  bantot = SUM(Dtset%nband(1:Dtset%nkpt*nsppol))
3305  ABI_MALLOC(doccde,(bantot))
3306  ABI_MALLOC(eigen,(bantot))
3307  ABI_MALLOC(occfact,(bantot))
3308  doccde(:)=zero; eigen(:)=zero; occfact(:)=zero
3309 
3310  jj=0; ibtot=0
3311  do isppol=1,nsppol
3312    do ikibz=1,Dtset%nkpt
3313      do ib=1,Hdr_wfk%nband(ikibz+(isppol-1)*Dtset%nkpt)
3314        ibtot=ibtot+1
3315        if (ib<=Sigp%nbnds) then
3316          jj=jj+1
3317          occfact(jj)=Hdr_wfk%occ(ibtot)
3318          eigen  (jj)=energies_p(ib,ikibz,isppol)
3319        end if
3320      end do
3321    end do
3322  end do
3323  ABI_FREE(energies_p)
3324 
3325  ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of
3326  ! symmetry operations in the old GW code (symmorphy and inversion)
3327  ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6))
3328  ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt')
3329 
3330  ABI_MALLOC(npwarr,(Dtset%nkpt))
3331  npwarr(:)=Sigp%npwwfn
3332 
3333  call ebands_init(bantot,ks_ebands,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,&
3334    doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,&
3335    Kmesh%nibz,npwarr,nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,&
3336    dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig,&
3337    dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
3338 
3339  ABI_FREE(doccde)
3340  ABI_FREE(eigen)
3341  ABI_FREE(npwarr)
3342 
3343  ! Calculate KS occupation numbers and ks_vbk(nkibz, nsppol) ====
3344  ! ks_vbk gives the (valence|last Fermi band) index for each k and spin.
3345  ! spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals
3346  call ebands_update_occ(ks_ebands, Dtset%spinmagntarget, prtvol=0)
3347 
3348  gaps = ebands_get_gaps(ks_ebands, gap_err)
3349  call gaps%print(unit=std_out)
3350  call ebands_report_gap(ks_ebands, unit=std_out)
3351 
3352  ABI_MALLOC(val_indeces,(ks_ebands%nkpt, nsppol))
3353  val_indeces = ebands_get_valence_idx(ks_ebands)
3354 
3355  ! Create Sigma header
3356  ! TODO Fix problems with symmorphy and k-points
3357  call hdr_init(ks_ebands,codvsn,Dtset,Hdr_out,Pawtab,pertcase0,Psps,wvl)
3358 
3359  ! Get Pawrhoij from the header of the WFK file
3360  ABI_MALLOC(Pawrhoij, (Cryst%natom*Dtset%usepaw))
3361  if (Dtset%usepaw==1) then
3362    call pawrhoij_alloc(Pawrhoij, 1, Dtset%nspden, Dtset%nspinor, nsppol, Cryst%typat, pawtab=Pawtab)
3363    call pawrhoij_copy(Hdr_wfk%Pawrhoij, Pawrhoij)
3364  end if
3365 
3366  call hdr_out%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
3367  ABI_FREE(occfact)
3368  call pawrhoij_free(Pawrhoij)
3369  ABI_FREE(Pawrhoij)
3370 
3371  ! ===========================================================
3372  ! ==== Setup of k-points and bands for the GW corrections ====
3373  ! ===========================================================
3374  ! * maxbdgw and minbdgw are the Max and min band index for GW corrections over k-points.
3375  !   They are used to dimension the wavefunctions and to calculate the matrix elements.
3376  !
3377  if (dtset%nkptgw == 0) then
3378    !
3379    ! Use qp_range to select the interesting k-points and the corresponing bands.
3380    !
3381    !    0 --> Compute the QP corrections only for the fundamental and the direct gap.
3382    ! +num --> Compute the QP corrections for all the k-points in the irreducible zone and include `num`
3383    !           bands above and below the Fermi level.
3384    ! -num --> Compute the QP corrections for all the k-points in the irreducible zone.
3385    !          Include all occupied states and `num` empty states.
3386 
3387    call wrtout(std_out, "nkptgw == 0 ==> Automatic selection of k-points and bands for the corrections.")
3388 
3389    if (gap_err /= 0 .and. dtset%gw_qprange == 0) then
3390      msg = "Problem while computing the fundamental and direct gap (likely metal). Will replace gw_qprange=0 with gw_qprange=1"
3391      ABI_WARNING(msg)
3392      dtset%gw_qprange = 1
3393    end if
3394 
3395    gw_qprange = dtset%gw_qprange
3396 
3397    if (dtset%ucrpa > 0) then
3398      dtset%nkptgw = Kmesh%nbz
3399      Sigp%nkptgw = dtset%nkptgw
3400      ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw))
3401      ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol))
3402      ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol))
3403      Sigp%kptgw(:,:) = Kmesh%bz(:,:)
3404      Sigp%minbnd = 1
3405      Sigp%maxbnd = Sigp%nbnds
3406 
3407    else if (gw_qprange /= 0) then
3408      ! Include all the k-points in the IBZ.
3409      ! Note that kptgw == ebands%kptns so we can use a single ik index in the loop over k-points
3410      ! No need to map kptgw onto ebands%kptns.
3411      dtset%nkptgw = Kmesh%nibz
3412      Sigp%nkptgw  = dtset%nkptgw
3413      ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw))
3414      ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol))
3415      ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol))
3416      Sigp%kptgw(:,:) = Kmesh%ibz(:,:)
3417      Sigp%minbnd = 1
3418      Sigp%maxbnd = Sigp%nbnds
3419 
3420      if (gw_qprange > 0) then
3421        ! All k-points: Add buffer of bands above and below the Fermi level.
3422        do spin=1,nsppol
3423          do ik=1,Sigp%nkptgw
3424            Sigp%minbnd(ik, spin) = MAX(val_indeces(ik,spin) - gw_qprange, 1)
3425            Sigp%maxbnd(ik, spin) = MIN(val_indeces(ik,spin) + gw_qprange + 1, Sigp%nbnds)
3426          end do
3427        end do
3428 
3429      else
3430        ! All k-points: include all occupied states and -gw_qprange empty states.
3431        Sigp%minbnd = 1
3432        do spin=1,nsppol
3433          do ik=1,Sigp%nkptgw
3434            Sigp%maxbnd(ik, spin) = MIN(val_indeces(ik,spin) - gw_qprange, Sigp%nbnds)
3435          end do
3436        end do
3437      end if
3438 
3439    else
3440      ! gw_qprange is not specified in the input.
3441      ! Include the direct and the fundamental KS gap.
3442      ! The main problem here is that kptgw and nkptgw do not depend on the spin and therefore
3443      ! we have compute the union of the k-points where the fundamental and the direct gaps are located.
3444      !
3445      ! Find the list of `interesting` kpoints.
3446      ABI_CHECK(gap_err == 0, "gw_qprange 0 cannot be used because I cannot find the gap (gap_err !=0)")
3447      nk_found = 1; kpos(1) = gaps%fo_kpos(1,1)
3448 
3449      do spin=1,nsppol
3450        do ifo=1,3
3451          ik = gaps%fo_kpos(ifo, spin)
3452          found = .FALSE.; jj = 0
3453          do while (.not. found .and. jj < nk_found)
3454            jj = jj + 1; found = (kpos(jj) == ik)
3455          end do
3456          if (.not. found) then
3457            nk_found = nk_found + 1; kpos(nk_found) = ik
3458          end if
3459        end do
3460      end do
3461 
3462      ! Now we can define the list of k-points and the bands range.
3463      dtset%nkptgw = nk_found
3464      Sigp%nkptgw = dtset%nkptgw
3465 
3466      ABI_MALLOC(Sigp%kptgw,(3, Sigp%nkptgw))
3467      ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol))
3468      ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol))
3469 
3470      do ii=1,Sigp%nkptgw
3471        ik = kpos(ii)
3472        Sigp%kptgw(:,ii) = Kmesh%ibz(:,ik)
3473        do spin=1,nsppol
3474          Sigp%minbnd(ii,spin) = val_indeces(ik, spin)
3475          Sigp%maxbnd(ii,spin) = val_indeces(ik, spin) + 1
3476        end do
3477      end do
3478    end if
3479 
3480  else
3481    ! Treat only the k-points and bands specified in the input file.
3482    Sigp%nkptgw = dtset%nkptgw
3483    ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw))
3484    ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol))
3485    ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol))
3486 
3487    do spin=1,nsppol
3488      Sigp%minbnd(:,spin)= dtset%bdgw(1,:,spin)
3489      Sigp%maxbnd(:,spin)= dtset%bdgw(2,:,spin)
3490    end do
3491 
3492    do ii=1,3
3493      do ikcalc=1,Sigp%nkptgw
3494        Sigp%kptgw(ii,ikcalc) = Dtset%kptgw(ii,ikcalc)
3495      end do
3496    end do
3497 
3498    do spin=1,nsppol
3499      do ikcalc=1,Sigp%nkptgw
3500        if (Dtset%bdgw(2,ikcalc,spin) > Sigp%nbnds) then
3501          write(msg,'(a,2i0,2(a,i0),2a,i0)')&
3502           "For (k,s) ",ikcalc,spin," bdgw= ",Dtset%bdgw(2,ikcalc,spin), " > nbnds=",Sigp%nbnds,ch10,&
3503           "Calculation will continue with bdgw =",Sigp%nbnds
3504          ABI_COMMENT(msg)
3505          Dtset%bdgw(2, ikcalc, spin) = Sigp%nbnds
3506        end if
3507      end do
3508    end do
3509 
3510  end if
3511 
3512  ! Make sure that all the degenerate states are included.
3513  ! * We will have to average the GW corrections over degenerate states if symsigma=1 is used.
3514  ! * KS states belonging to the same irreducible representation should be included in the basis set used for SCGW.
3515  if (Sigp%symsigma /=0 .or. gwcalctyp >= 10) then
3516    do isppol=1,nsppol
3517      do ikcalc=1,Sigp%nkptgw
3518 
3519        if (kmesh%has_IBZ_item(Sigp%kptgw(:,ikcalc), ikibz, G0)) then
3520          call ebands_enclose_degbands(ks_ebands,ikibz,isppol, &
3521                Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol),changed,tol_enediff)
3522 
3523          if (changed) then
3524            write(msg,'(2(a,i0),2a,2(1x,i0))')&
3525             "Not all the degenerate states at ikcalc= ",ikcalc,", spin= ",isppol, ch10, &
3526             "were included in the bdgw set. bdgw has been changed to: ",Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol)
3527            ABI_COMMENT(msg)
3528          end if
3529        else
3530          write(msg,'(3a)')&
3531           ' not in the list of k points treated in the preparatory SCF run.',ch10, &
3532           ' Change kptgw, or shiftk of previous run.'
3533          ABI_ERROR(sjoin('k-point', ktoa(Sigp%kptgw(:,ikcalc)),trim(msg)))
3534        end if
3535 
3536      end do
3537    end do
3538  end if
3539 
3540  Sigp%minbdgw = MINVAL(Sigp%minbnd)
3541  Sigp%maxbdgw = MAXVAL(Sigp%maxbnd)
3542 
3543  ! Check if there are duplicated k-point in Sigp%
3544  do ii=1,Sigp%nkptgw
3545    do jj=ii+1,Sigp%nkptgw
3546      if (isamek(Sigp%kptgw(:,ii), Sigp%kptgw(:,jj), G0)) then
3547        write(msg,'(5a)')&
3548         'kptgw contains duplicated k-points. This is not allowed since ',ch10,&
3549         'the QP corrections for this k-point will be calculated more than once. ',ch10,&
3550         'Check your input file. '
3551        ABI_ERROR(msg)
3552      end if
3553    end do
3554  end do
3555 
3556  !=== Check if the k-points are in the BZ ===
3557  ! FB: Honestly the code is not able to treat k-points, which are not in the input list of points in the IBZ.
3558  ! This extension should require to change the code in different places.
3559  ! Therefore, one should by now prevent the user from calculating sigma for a k-point not in the IBZ.
3560  ABI_MALLOC(Sigp%kptgw2bz, (Sigp%nkptgw))
3561 
3562  do ikcalc=1,Sigp%nkptgw
3563    if (kmesh%has_BZ_item(Sigp%kptgw(:,ikcalc), ikcalc2bz, G0)) then
3564      Sigp%kptgw2bz(ikcalc) = ikcalc2bz
3565    else
3566      write(msg,'(3a)')&
3567       ' not in the list of k points treated in the preparatory SCF run.',ch10,' Change kptgw, or shiftk of previous run.'
3568      ABI_ERROR(sjoin('k-point:', ktoa(Sigp%kptgw(:,ikcalc)),trim(msg)))
3569    end if
3570  end do
3571 
3572  ! Warn the user if SCGW run and not all the k-points are included.
3573  if (gwcalctyp >= 10 .and. Sigp%nkptgw /= Hdr_wfk%nkpt) then
3574    write(msg,'(3a,2(a,i0),2a)')ch10,&
3575     " COMMENT: In a self-consistent GW run, the QP corrections should be calculated for all the k-points of the KSS file ",ch10,&
3576     " but nkptgw= ",Sigp%nkptgw," and WFK nkpt= ",Hdr_wfk%nkpt,ch10,&
3577     " Assuming expert user. Execution will continue. "
3578    call wrtout(ab_out, msg)
3579  end if
3580 
3581  ! Setup of the table used in the case of SCGW on wavefunctions to reduce the number
3582  ! of elements <i,kgw,s|\Sigma|j,kgw,s> that have to be calculated. No use of symmetries, except for Hermiticity.
3583  call sigma_tables(Sigp, Kmesh)
3584 
3585  ! === Read external file and initialize basic dimension of Er% ===
3586  ! TODO use mqmem as input variable instead of gwmem
3587 
3588  ! === If required, use a matrix for $\Sigma_c$ which is smaller than that stored on file ===
3589  ! * By default the entire matrix is read and used,
3590  ! * Define consistently npweps and ecuteps for \Sigma_c according the input
3591  if (Dtset%npweps>0.or.Dtset%ecuteps>0) then
3592    if (Dtset%npweps>0) Dtset%ecuteps=zero
3593    nsheps=0
3594    call setshells(Dtset%ecuteps,Dtset%npweps,nsheps,Dtset%nsym,gmet,gprimd,Dtset%symrel,'eps',ucvol)
3595  end if
3596 
3597  mqmem=0; if (Dtset%gwmem/10==1) mqmem=1
3598 
3599  if (dtset%getscr /=0 .or. dtset%irdscr/=0 .or. dtset%getscr_filepath /= ABI_NOFILE) then
3600    fname=Dtfil%fnameabi_scr
3601  else if (Dtset%getsuscep/=0.or.Dtset%irdsuscep/=0) then
3602    fname=Dtfil%fnameabi_sus
3603  else
3604    fname=Dtfil%fnameabi_scr
3605    !FIXME this has to be cleaned, in tgw2_3 Dtset%get* and Dtset%ird* are  not defined
3606    !ABI_ERROR("getsuscep or irdsuscep are not defined")
3607  end if
3608  !
3609  ! === Setup of q-mesh in the whole BZ ===
3610  ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed,
3611  !   epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed)
3612  !
3613  if (sigma_needs_w(Sigp)) then
3614    if (.not. file_exists(fname)) then
3615      fname = nctk_ncify(fname)
3616      ABI_COMMENT(sjoin("File not found. Will try netcdf file:", fname))
3617    end if
3618 
3619    call init_Er_from_file(Er,fname,mqmem,Dtset%npweps,comm)
3620 
3621    Sigp%npwc=Er%npwe
3622    if (Sigp%npwc>Sigp%npwx) then
3623      Sigp%npwc=Sigp%npwx
3624      ABI_COMMENT("Found npw_correlation > npw_exchange, Imposing npwc=npwx")
3625      ! There is a good reason for doing so, see csigme.F90 and the size of the arrays
3626      ! rhotwgp and rhotwgp: we need to define a max size and we opt for Sigp%npwx.
3627    end if
3628    Er%npwe=Sigp%npwc
3629    Dtset%npweps=Er%npwe
3630    call Qmesh%init(Cryst,Er%nqibz,Er%qibz,Dtset%kptopt)
3631 
3632  else
3633    Er%npwe     =1
3634    Sigp%npwc   =1
3635    Dtset%npweps=1
3636    call find_qmesh(Qmesh,Cryst,Kmesh)
3637    ABI_MALLOC(Er%gvec,(3,1))
3638    Er%gvec(:,1) = [0, 0, 0]
3639  end if
3640 
3641  call Qmesh%print("Q-mesh for screening function",std_out,Dtset%prtvol,"COLL")
3642  call Qmesh%print("Q-mesh for screening function",ab_out ,0           ,"COLL")
3643 
3644 
3645  do iqbz=1,Qmesh%nbz
3646    call qmesh%get_BZ_item(iqbz, q_bz, iq_ibz, isym, itim, umklp=q_umklp)
3647    !print *, "iqbz, q_bz, iq_ibz, isym, itim, umklp"
3648    !print *, iqbz, q_bz, iq_ibz, isym, itim, q_umklp
3649 
3650    if (ANY(q_umklp/=0)) then
3651      sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz))
3652      write(std_out,*) sq,Qmesh%bz(:,iqbz)
3653      write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')&
3654        'qpoint ',Qmesh%bz(:,iqbz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,&
3655        'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,&
3656        'however a non-zero umklapp G_o vector is required and this is not yet allowed'
3657      ABI_ERROR(msg)
3658    end if
3659  end do
3660  !stop
3661  !
3662  ! === Find optimal value for G-sphere enlargment due to oscillator matrix elements ===
3663  ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true size bz is buggy
3664  ! * -one is used because we loop over all the possibile differences, unlike screening
3665 
3666  call get_ng0sh(Sigp%nkptgw,Sigp%kptgw,Kmesh%nbz,Kmesh%bz,Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt)
3667  call wrtout(std_out, sjoin(' Optimal value for ng0sh ', ltoa(ng0sh_opt)))
3668  Sigp%mG0 = ng0sh_opt
3669 
3670 ! G-sphere for W and Sigma_c is initialized from the SCR file.
3671  call Gsph_c%init(Cryst, Er%npwe, gvec=Er%gvec)
3672  call Gsph_x%init(Cryst, Sigp%npwx, gvec=gvec_kss)
3673  Sigp%ecuteps = Gsph_c%ecut
3674  Dtset%ecuteps = Sigp%ecuteps
3675 
3676 ! === Make biggest G-sphere of Sigp%npwvec vectors ===
3677  Sigp%npwvec=MAX(Sigp%npwwfn,Sigp%npwx)
3678  call Gsph_Max%init(Cryst, Sigp%npwvec, gvec=gvec_kss)
3679 !BEGINDEBUG
3680  ! Make sure that the two G-spheres are equivalent.
3681  ierr=0
3682  if (sigma_needs_w(Sigp)) then
3683    ng = MIN(SIZE(Gsph_c%gvec,DIM=2),SIZE(gvec_kss,DIM=2))
3684    do ig1=1,ng
3685      if (ANY(Gsph_c%gvec(:,ig1)/=gvec_kss(:,ig1))) then
3686        ierr=ierr+1
3687        write(std_out,*)" Gsph_c, gvec_kss ",ig1,"/",ng,Gsph_c%gvec(:,ig1),gvec_kss(:,ig1)
3688      end if
3689    end do
3690    ABI_CHECK(ierr==0,"Mismatch between Gsph_c and gvec_kss")
3691  end if
3692  ierr=0
3693  ng = MIN(SIZE(Gsph_x%gvec,DIM=2),SIZE(gvec_kss,DIM=2))
3694  do ig1=1,ng
3695    if (ANY(Gsph_x%gvec(:,ig1)/=gvec_kss(:,ig1))) then
3696      ierr=ierr+1
3697      write(std_out,*)" Gsph_x, gvec_kss ",ig1,"/",ng,Gsph_x%gvec(:,ig1),gvec_kss(:,ig1)
3698    end if
3699  end do
3700  ABI_CHECK(ierr==0,"Mismatch between Gsph_x and gvec_kss")
3701 !ENDDEBUG
3702 
3703  ABI_FREE(gvec_kss)
3704  !
3705  ! === Get Fourier components of the Coulomb term for all q-points in the IBZ ===
3706  ! * If required, use a cutoff in the interaction
3707  ! * Pcv%vc_sqrt contains Vc^{-1/2}
3708  ! * Setup also the analytical calculation of the q->0 component
3709  ! FIXME recheck ngfftf since I got different charge outside the cutoff region
3710 
3711  if (Dtset%gw_nqlwl==0) then
3712    nqlwl=1
3713    ABI_MALLOC(qlwl,(3,nqlwl))
3714    qlwl(:,1)= GW_Q0_DEFAULT
3715  else
3716    nqlwl=Dtset%gw_nqlwl
3717    ABI_MALLOC(qlwl,(3,nqlwl))
3718    qlwl(:,:)=Dtset%gw_qlwl(:,1:nqlwl)
3719  end if
3720 
3721  ! The Coulomb interaction used here might have two terms:
3722  ! the first term generates the usual sigma self-energy, but possibly, one should subtract
3723  ! from it the Coulomb interaction already present in the Kohn-Sham basis,
3724  ! if the usefock associated to ixc is one.
3725  ! The latter excludes (in the present implementation) mod(Dtset%gwcalctyp,10)==5
3726  nvcoul_init=1
3727  call get_xclevel(Dtset%ixc,xclevel_ixc,usefock_ixc)
3728 
3729  if(usefock_ixc==1)then
3730    nvcoul_init=2
3731    if(mod(Dtset%gwcalctyp,10)==5)then
3732      write(msg,'(4a,i3,a,i3,4a,i5)')ch10,&
3733      ' The starting wavefunctions were obtained from self-consistent calculations in the planewave basis set',ch10,&
3734      ' with ixc = ',Dtset%ixc,' associated with usefock =',usefock_ixc,ch10,&
3735      ' In this case, the present implementation does not allow that the self-energy for sigma corresponds to',ch10,&
3736      '  mod(gwcalctyp,10)==5, while your gwcalctyp= ',Dtset%gwcalctyp
3737      ABI_ERROR(msg)
3738    endif
3739  endif
3740  do ivcoul_init=1,nvcoul_init
3741    rcut = Dtset%rcut
3742    icutcoul_eff=Dtset%gw_icutcoul
3743    Sigp%sigma_mixing=one
3744    if( mod(Dtset%gwcalctyp,10)==5 .or. ivcoul_init==2)then
3745      if(abs(Dtset%hyb_mixing)>tol8)then
3746        ! Warning: the absolute value is needed, because of the singular way used to define the default for this input variable
3747        Sigp%sigma_mixing=abs(Dtset%hyb_mixing)
3748      else if(abs(Dtset%hyb_mixing_sr)>tol8)then
3749        Sigp%sigma_mixing=abs(Dtset%hyb_mixing_sr)
3750        icutcoul_eff=5
3751      endif
3752      if(abs(rcut)<tol6 .and. abs(Dtset%hyb_range_fock)>tol8)rcut=one/Dtset%hyb_range_fock
3753    endif
3754 
3755    if (ivcoul_init == 1) then
3756 
3757      if (Gsph_x%ng > Gsph_c%ng) then
3758        call Vcp%init(Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3759         Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm)
3760      else
3761        call Vcp%init(Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3762         Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm)
3763      end if
3764 
3765    else
3766 
3767      ! Use a temporary Vcp_ks to compute the Coulomb interaction already present
3768      ! in the Fock part of the Kohn-Sham Hamiltonian
3769      if (Gsph_x%ng > Gsph_c%ng) then
3770        call Vcp_ks%init(Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3771                         Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm)
3772      else
3773        call Vcp_ks%init(Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3774                         Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm)
3775      end if
3776 
3777      ! Now compute the residual Coulomb interaction.
3778      Vcp%vc_sqrt_resid=sqrt(Vcp%vc_sqrt**2-Sigp%sigma_mixing*Vcp_ks%vc_sqrt**2)
3779      Vcp%i_sz_resid=Vcp%i_sz-Sigp%sigma_mixing*Vcp_ks%i_sz
3780 
3781      ! The mixing factor has already been accounted for, so set it back to one
3782      Sigp%sigma_mixing=one
3783      call Vcp_ks%free()
3784 
3785      !write(std_out,'(a)')' setup_sigma : the residual Coulomb interaction has been computed'
3786    endif
3787 
3788   !#else
3789   !   call Vcp%init(Gsph_Max,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,ivcoul_init,Dtset%vcutgeo,&
3790   !&    Dtset%ecutsigx,Sigp%npwx,nqlwl,qlwl,ngfftf,comm)
3791   !#endif
3792 
3793  enddo
3794 
3795 #if 0
3796  ! Using the random q for the optical limit is one of the reasons
3797  ! why sigma breaks the initial energy degeneracies.
3798  Vcp%i_sz=zero
3799  Vcp%vc_sqrt(1,:)=czero
3800  Vcp%vcqlwl_sqrt(1,:)=czero
3801 #endif
3802 
3803  ABI_FREE(qlwl)
3804 
3805  Sigp%ecuteps = Dtset%ecuteps
3806  Sigp%ecutwfn = Dtset%ecutwfn
3807  Sigp%ecutsigx = Dtset%ecutsigx
3808 
3809  ! === Setup of the FFT mesh for the oscilator strengths ===
3810  ! * Init gwc_ngfft(7:18) and gwx_ngfft(7:18) with Dtset%ngfft(7:18)
3811  ! * Here we redefine gwc_ngfft(1:6) according to the following options:
3812  !
3813  ! method == 0 --> FFT grid read from fft.in (debugging purpose)
3814  ! method == 1 --> Normal FFT mesh
3815  ! method == 2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
3816  ! method == 3 --> Doubled FFT grid, same as the the FFT for the density,
3817  !
3818  ! enforce_sym == 1 --> Enforce a FFT mesh compatible with all the symmetry operation and FFT library
3819  ! enforce_sym == 0 --> Find the smallest FFT grid compatbile with the library, do not care about symmetries
3820  !
3821  gwc_ngfft(1:18) = Dtset%ngfft(1:18)
3822  gwx_ngfft(1:18) = Dtset%ngfft(1:18)
3823 
3824  method = 2
3825  if (Dtset%fftgw == 00 .or. Dtset%fftgw == 01) method = 0
3826  if (Dtset%fftgw == 10 .or. Dtset%fftgw == 11) method = 1
3827  if (Dtset%fftgw == 20 .or. Dtset%fftgw == 21) method = 2
3828  if (Dtset%fftgw == 30 .or. Dtset%fftgw == 31) method = 3
3829  enforce_sym = mod(dtset%fftgw, 10)
3830 
3831  ! FFT mesh for sigma_x.
3832  call setmesh(gmet, Gsph_Max%gvec, gwx_ngfft, Sigp%npwvec, Sigp%npwx, Sigp%npwwfn, &
3833               gwx_nfftot, method, Sigp%mG0, Cryst, enforce_sym)
3834 
3835  ! FFT mesh for sigma_c.
3836  call setmesh(gmet, Gsph_Max%gvec, gwc_ngfft, Sigp%npwvec, Er%npwe, Sigp%npwwfn,&
3837               gwc_nfftot, method, Sigp%mG0, Cryst, enforce_sym, unit=dev_null)
3838 
3839  ! ======================================================================
3840  ! ==== Check for presence of files with core orbitals, for PAW only ====
3841  ! ======================================================================
3842  Sigp%use_sigxcore=0
3843  if (Dtset%usepaw==1.and.Dtset%gw_sigxcore==1) then
3844    ii = 0
3845    do itypat=1,Cryst%ntypat
3846      string = Psps%filpsp(itypat)
3847      fcore = "CORE_"//TRIM(basename(string))
3848      ic = INDEX (TRIM(string), "/" , back=.TRUE.)
3849      if (ic>0 .and. ic<LEN_TRIM(string)) then
3850        ! string defines a path, prepend path to fcore
3851        fcore = Psps%filpsp(itypat)(1:ic)//TRIM(fcore)
3852      end if
3853      if (file_exists(fcore)) then
3854        ii = ii+1
3855      else
3856        ABI_WARNING(sjoin("HF decoupling is required but cannot find file:", fcore))
3857      end if
3858    end do
3859 
3860    Sigp%use_sigxcore=1
3861    if (ii/=Cryst%ntypat) then
3862      ABI_ERROR("Files with core orbitals not found")
3863    end if
3864  end if ! PAW+HF decoupling
3865  !
3866  ! ==============================
3867  ! ==== Extrapolar technique ====
3868  ! ==============================
3869  Sigp%gwcomp   = Dtset%gwcomp
3870  Sigp%gwencomp = Dtset%gwencomp
3871 
3872  if (Sigp%gwcomp==1) then
3873    write(msg,'(6a,e11.4,a)')ch10,&
3874     'Using the extrapolar approximation to accelerate convergence',ch10,&
3875     'with respect to the number of bands included',ch10,&
3876     'with gwencomp: ',Sigp%gwencomp*Ha_eV,' [eV]'
3877    call wrtout(std_out, msg)
3878  end if
3879  !
3880  ! ===================================
3881  ! ==== Final compatibility tests ====
3882  ! ===================================
3883  ltest=(ks_ebands%mband == Sigp%nbnds .and. ALL(ks_ebands%nband == Sigp%nbnds))
3884  ABI_CHECK(ltest,'BUG in definition of ks_ebands%nband')
3885 
3886  ! FIXME
3887  if (Dtset%symsigma/=0 .and. Sigp%nomegasr/=0) then
3888    if (cryst%idx_spatial_inversion() == 0) then
3889      write(msg,'(5a)')' setup_sigma : BUG :',ch10,&
3890       'It is not possible to use symsigma /= 0 to calculate the spectral function ',ch10,&
3891       'when the system does not have the spatial inversion. Please use symsigma=0 '
3892      ABI_WARNING(msg)
3893    end if
3894  end if
3895 
3896  if (mod10==SIG_GW_AC) then
3897    !if (Sigp%gwcalctyp/=1) ABI_ERROR("Self-consistency with AC not implemented") ! MRM: let's allow it
3898    if (Sigp%gwcomp==1) ABI_ERROR("AC with extrapolar technique not implemented")
3899  end if
3900 
3901  call gaps%free()
3902 
3903  ABI_FREE(val_indeces)
3904 
3905  DBG_EXIT('COLL')
3906 
3907 end subroutine setup_sigma

ABINIT/setup_vcp [ Functions ]

[ Top ] [ Functions ]

NAME

 setup_vcp

FUNCTION

  Initialize the Vcp and compute the corresponding residue.

INPUTS

 Dtset<type(dataset_type)>=all input variables for this dataset
 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
 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.
 comm=Information about the xmpi_world

OUTPUT

 Vcp_full<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique.
 Vcp_ks<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique.
 coef_hyb=real variable containing the amount of GLOBAL hybridization

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.

SOURCE

3997 subroutine sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr)
3998 
3999 !Arguments ------------------------------------
4000 !scalars
4001  integer,intent(in) :: my_rank,nprocs
4002  integer,intent(out) :: ierr
4003  type(Dataset_type),intent(in) :: Dtset
4004  type(sigparams_t),intent(in) :: Sigp
4005  type(kmesh_t),intent(in) :: Kmesh
4006 !arrays
4007  integer,allocatable,intent(out) :: my_spins(:)
4008  logical,intent(out) :: bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
4009  logical,intent(out) :: keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
4010 
4011 !Local variables-------------------------------
4012 !scalars
4013  integer :: my_nspins,my_maxb,my_minb,isp,spin,band,nsppol,rank_spin
4014  logical :: store_ur
4015 !arrays
4016  integer :: tmp_spins(Sigp%nsppol),nprocs_spin(Sigp%nsppol)
4017 
4018 ! *************************************************************************
4019 
4020  ierr=0; nsppol=Sigp%nsppol
4021 
4022  ! List of spins for each node, number of processors per each spin
4023  ! and the MPI rank in the "spin" communicator.
4024  my_nspins=nsppol
4025  ABI_MALLOC(my_spins, (nsppol))
4026  my_spins= [(isp, isp=1,nsppol)]
4027  nprocs_spin = nprocs; rank_spin = my_rank
4028 
4029  if (nsppol==2 .and. nprocs>1) then
4030    ! Distribute spins (optimal distribution if nprocs is even)
4031    nprocs_spin(1) = nprocs/2
4032    nprocs_spin(2) = nprocs - nprocs/2
4033    my_nspins=1
4034    my_spins(1)=1
4035    if (my_rank+1>nprocs/2) then
4036      ! I will treate spin=2, compute shifted rank.
4037      my_spins(1)=2
4038      rank_spin = my_rank - nprocs/2
4039    end if
4040  end if
4041 
4042  store_ur = (MODULO(Dtset%gwmem,10)==1)
4043  keep_ur=.FALSE.; bks_mask=.FALSE.
4044 
4045  select case (Dtset%gwpara)
4046  case (1)
4047    ! Parallelization over transitions **without** memory distributions (Except for the spin).
4048    my_minb=1; my_maxb=Sigp%nbnds
4049    if (dtset%ucrpa>0) then
4050      bks_mask(my_minb:my_maxb,:,:)=.TRUE.
4051      if (store_ur) keep_ur(my_minb:my_maxb,:,:)=.TRUE.
4052    else
4053      do isp=1,my_nspins
4054        spin = my_spins(isp)
4055        bks_mask(my_minb:my_maxb,:,spin)=.TRUE.
4056        if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
4057      end do
4058    end if
4059 
4060  case (2)
4061    ! Distribute bands and spin (alternating planes of bands)
4062    do isp=1,my_nspins
4063      spin = my_spins(isp)
4064      do band=1,Sigp%nbnds
4065        if (xmpi_distrib_with_replicas(band,Sigp%nbnds,rank_spin,nprocs_spin(spin))) then
4066        !if (MODULO(band-1,nprocs_spin(spin))==rank_spin) then
4067          bks_mask(band,:,spin)=.TRUE.
4068          if (store_ur) keep_ur(band,:,spin)=.TRUE.
4069        end if
4070      end do
4071    end do
4072 
4073 #if 0
4074    ! Each node store the full set of occupied states to speed up Sigma_x.
4075    do isp=1,my_nspins
4076      spin = my_spins(isp)
4077      do ik_ibz=1,Kmesh%nibz
4078        ks_iv=ks_vbik(ik_ibz,spin) ! Valence index for this (k,s)
4079        bks_mask(1:ks_iv,:,spin)=.TRUE.
4080        if (store_ur) keep_ur(1:ks_iv,:,spin)=.TRUE.
4081      end do
4082    end do
4083 #endif
4084 
4085  case default
4086    ABI_WARNING("Wrong value for gwpara")
4087    ierr = 1
4088  end select
4089 
4090  ! Return my_spins with correct size.
4091  tmp_spins(1:my_nspins) = my_spins(1:my_nspins)
4092 
4093  ABI_FREE(my_spins)
4094  ABI_MALLOC(my_spins, (my_nspins))
4095  my_spins = tmp_spins(1:my_nspins)
4096 
4097 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.
 [esymm(Kmesh%nibz,Sigp%nsppol)]= 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.

SOURCE

3931 subroutine sigma_tables(Sigp, Kmesh, esymm)
3932 
3933   use m_sigtk, only : sigtk_sigma_tables
3934 
3935 !Arguments ------------------------------------
3936 !scalars
3937  type(sigparams_t),target,intent(inout) :: Sigp
3938  type(kmesh_t),intent(in) :: Kmesh
3939 !arrays
3940  type(esymm_t),optional,intent(in) :: esymm(Kmesh%nibz, Sigp%nsppol)
3941 
3942 !Local variables-------------------------------
3943 !scalars
3944  integer :: ikcalc,ik_ibz
3945  logical :: sigc_is_herm, only_diago
3946 !arrays
3947  integer,allocatable :: kcalc2ibz(:)
3948 
3949 ! *************************************************************************
3950 
3951  only_diago = sigp%gwcalctyp < 20
3952  sigc_is_herm = sigma_is_herm(Sigp)
3953 
3954  ABI_MALLOC(kcalc2ibz, (sigp%nkptgw))
3955  do ikcalc=1,sigp%nkptgw
3956    ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
3957    kcalc2ibz(ikcalc) = ik_ibz
3958  end do
3959 
3960  if (present(esymm)) then
3961    call sigtk_sigma_tables(sigp%nkptgw, kmesh%nibz, sigp%nsppol, sigp%minbnd, sigp%maxbnd, kcalc2ibz, &
3962                            only_diago, sigc_is_herm, sigp%sigxij_tab, sigp%sigcij_tab, esymm=esymm)
3963  else
3964    call sigtk_sigma_tables(sigp%nkptgw, kmesh%nibz, sigp%nsppol, sigp%minbnd, sigp%maxbnd, kcalc2ibz, &
3965                            only_diago, sigc_is_herm, sigp%sigxij_tab, sigp%sigcij_tab)
3966  end if
3967 
3968  ABI_FREE(kcalc2ibz)
3969 
3970 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

  Output is written on the main abinit output file. Some results are stored in external files

NOTES

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

SOURCE

 171 subroutine sigma(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim)
 172 
 173 !Arguments ------------------------------------
 174 !scalars
 175  character(len=8),intent(in) :: codvsn
 176  type(Datafiles_type),intent(in) :: Dtfil
 177  type(Dataset_type),intent(inout) :: Dtset
 178  type(Pawang_type),intent(inout) :: Pawang
 179  type(Pseudopotential_type),intent(inout) :: Psps
 180 !arrays
 181  real(dp),intent(in) :: acell(3),rprim(3,3)
 182  type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw)
 183  type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw)
 184 
 185 !Local variables-------------------------------
 186 !scalars
 187  integer,parameter :: tim_fourdp5 = 5, master = 0, cplex1 = 1, ipert0 = 0, idir0 = 0,  optrhoij1 = 1
 188  integer :: approx_type,b1gw,b2gw,cplex,cplex_dij,cplex_rhoij !,band
 189  integer :: dim_kxcg,gwcalctyp,gnt_option,has_dijU,has_dijso,iab,bmin,bmax,irr_idx1,irr_idx2
 190  integer :: iat,ib,ib1,ib2,ic,id_required,ider,ii,ik,ierr,ount
 191  integer :: ik_bz,ikcalc,ik_ibz,ikxc,npw_k,omp_ncpus,pwx,ibz
 192  integer :: isp,is_idx,istep,itypat,itypatcor,izero,jj,first_band,last_band
 193  integer :: ks_iv,lcor,lmn2_size_max,mband,my_nband
 194  integer :: mgfftf,mod10,moved_atm_inside,moved_rhor,n3xccc !,mgfft
 195  integer :: nbsc,ndij,ndim,nfftf,nfftf_tot,nkcalc,gwc_nfft,gwc_nfftot,gwx_nfft,gwx_nfftot
 196  integer :: ngrvdw,nhatgrdim,nkxc,nkxc1,nprocs,nscf,nspden_rhoij,nzlmopt,optene
 197  integer :: optcut,optgr0,optgr1,optgr2,option,option_test,option_dij,optrad,psp_gencond
 198  integer :: my_rank,rhoxsp_method,comm,use_aerhor,use_umklp,usexcnhat
 199  integer :: ioe0j,spin,io,jb,nomega_sigc
 200  integer :: temp_unt,ncid
 201  integer :: work_size,nstates_per_proc,my_nbks
 202  !integer :: jb_qp,ib_ks,ks_irr
 203  integer :: gw1rdm,x1rdm
 204  real(dp) :: compch_fft,compch_sph,r_s,rhoav,alpha
 205  real(dp) :: drude_plsmf,my_plsmf,ecore,ecut_eff,ecutdg_eff,ehartree
 206  real(dp) :: etot_sd,etot_mbb,evextnl_energy,ex_energy,gsqcutc_eff,gsqcutf_eff,gsqcut_shp,norm,oldefermi
 207  real(dp) :: eh_energy,ekin_energy,evext_energy,den_int,coef_hyb,exc_mbb_energy,tol_empty
 208  real(dp) :: ucvol,ucvol_local,vxcavg,vxcavg_qp
 209  real(dp) :: gwc_gsq,gwx_gsq,gw_gsq, gsqcut,boxcut,ecutf
 210  real(dp) :: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem
 211  complex(dpc) :: max_degw,cdummy
 212  logical :: wfknocheck,rdm_update,readchkprdm,prtchkprdm
 213  logical :: use_paw_aeur,dbg_mode,pole_screening,call_pawinit,is_dfpt=.false.
 214  character(len=500) :: msg
 215  character(len=fnlen) :: wfk_fname,pawden_fname,gw1rdm_fname
 216  type(kmesh_t) :: Kmesh,Qmesh
 217  type(ebands_t) :: ks_ebands, qp_ebands
 218  type(vcoul_t) :: Vcp, Vcp_ks, Vcp_full
 219  type(crystal_t) :: Cryst
 220  type(Energies_type) :: KS_energies,QP_energies
 221  type(Epsilonm1_results) :: Er
 222  type(gsphere_t) :: Gsph_Max,Gsph_x,Gsph_c
 223  type(hdr_type) :: Hdr_wfk,Hdr_sigma,Hdr_rhor
 224  type(melflags_t) :: KS_mflags,QP_mflags
 225  type(melements_t) :: KS_me, QP_me, GW1RDM_me
 226  type(MPI_type) :: MPI_enreg_seq
 227  type(paw_dmft_type) :: Paw_dmft
 228  type(pawfgr_type) :: Pawfgr
 229  type(ppmodel_t) :: PPm
 230  type(sigparams_t) :: Sigp
 231  type(sigma_t) :: Sr
 232  type(wfdgw_t),target :: Wfd,Wfdf,Wfd_nato_master
 233  type(wfdgw_t),pointer :: Wfd_nato_all
 234  type(wave_t),pointer :: wave
 235  type(wvl_data) :: Wvl
 236 !arrays
 237  integer :: gwc_ngfft(18),ngfftc(18),ngfftf(18),gwx_ngfft(18)
 238  integer,allocatable :: sigmak_todo(:)
 239  integer,allocatable :: nq_spl(:),nlmn_atm(:),my_spins(:)
 240  integer,allocatable :: tmp_gfft(:,:),ks_vbik(:,:),nband(:,:),l_size_atm(:),qp_vbik(:,:)
 241  integer,allocatable :: tmp_kstab(:,:,:),ks_irreptab(:,:,:),qp_irreptab(:,:,:),my_band_list(:)
 242  real(dp),parameter ::  k0(3)=zero
 243  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),strsxc(6),tsec(2)
 244  real(dp),allocatable :: weights(:),nat_occs(:,:),gw_rhor(:,:),gw_rhog(:,:),gw_vhartr(:)
 245  real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:)
 246  real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:)
 247  real(dp),allocatable :: ks_rhor(:,:),ks_vhartr(:),ks_vtrial(:,:),ks_vxc(:,:), ks_taur(:,:)
 248  real(dp),allocatable :: kxc(:,:),qp_kxc(:,:),ph1d(:,:),ph1df(:,:)
 249  real(dp),allocatable :: prev_rhor(:,:),prev_taur(:,:),qp_nhat(:,:)
 250  real(dp),allocatable :: qp_nhatgr(:,:,:),qp_rhog(:,:),qp_rhor_paw(:,:)
 251  real(dp),allocatable :: qp_rhor_n_one(:,:),qp_rhor_nt_one(:,:)
 252  real(dp),allocatable :: qp_rhor(:,:),qp_vhartr(:),qp_vtrial(:,:),qp_vxc(:,:)
 253  real(dp),allocatable :: qp_taur(:,:),igwene(:,:,:)
 254  real(dp),allocatable :: vpsp(:),xccc3d(:),dijexc_core(:,:,:),dij_hf(:,:,:)
 255  real(dp),allocatable :: nl_bks(:,:,:)
 256  !real(dp),allocatable :: osoc_bks(:, :, :)
 257  real(dp),allocatable :: ks_aepaw_rhor(:,:) !,ks_n_one_rhor(:,:),ks_nt_one_rhor(:,:)
 258  complex(dpc) :: ovlp(2)
 259  complex(dpc),allocatable :: ctmp(:,:),hbare(:,:,:,:)
 260  complex(dpc),target,allocatable :: sigcme(:,:,:,:,:)
 261  complex(dpc),allocatable :: hdft(:,:,:,:),htmp(:,:,:,:),uks2qp(:,:)
 262  complex(dpc),allocatable :: xrdm_k_full(:,:,:),rdm_k(:,:),pot_k(:,:),nateigv(:,:,:,:)
 263  complex(dpc),allocatable :: old_ks_purex(:,:),new_hartr(:,:)
 264  complex(gwpc),allocatable :: kxcg(:,:),fxc_ADA(:,:,:)
 265  complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:)
 266  complex(dpc),allocatable :: sigcme_k(:,:,:,:)
 267  complex(dpc), allocatable :: rhot1_q_m(:,:,:,:,:,:,:)
 268  complex(dpc), allocatable :: M1_q_m(:,:,:,:,:,:,:)
 269  logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:),bmask(:)
 270  logical,allocatable :: bdm_mask(:,:,:),bdm2_mask(:,:,:)
 271  type(esymm_t),target,allocatable :: KS_sym(:,:)
 272  type(esymm_t),pointer :: QP_sym(:,:)
 273  type(pawcprj_type),allocatable :: Cp1(:,:) !,Cp2(:,:)
 274  type(littlegroup_t),allocatable :: Ltg_k(:)
 275  type(Paw_an_type),allocatable :: KS_paw_an(:),QP_paw_an(:)
 276  type(Paw_ij_type),allocatable :: KS_paw_ij(:),QP_paw_ij(:)
 277  type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
 278  type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:),QP_pawrhoij(:),prev_Pawrhoij(:),tmp_pawrhoij(:)
 279  type(pawpwff_t),allocatable :: Paw_pwff(:)
 280  type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
 281  type(plowannier_type) :: wanbz,wanibz,wanibz_in
 282  type(operwan_realspace_type), allocatable :: rhot1(:,:)
 283 
 284 !************************************************************************
 285 
 286  DBG_ENTER('COLL')
 287 
 288  call timab(401,1,tsec) ! sigma(Total)
 289  call timab(402,1,tsec) ! sigma(Init1)
 290 
 291  write(msg,'(7a)')&
 292  ' SIGMA: Calculation of the GW corrections ',ch10,ch10,&
 293  ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,&
 294  ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.'
 295  call wrtout([std_out, ab_out], msg)
 296 
 297  if(dtset%ucrpa>0) then
 298    write(msg,'(6a)')ch10,&
 299     ' cRPA Calculation: Calculation of the screened Coulomb interaction (ucrpa/=0) ',ch10
 300    call wrtout([std_out, ab_out], msg)
 301  end if
 302 
 303 #if defined HAVE_GW_DPC
 304  if (gwpc /= 8) then
 305    write(msg,'(6a)')ch10,&
 306     'Number of bytes for double precision complex /=8 ',ch10,&
 307     'Cannot continue due to kind mismatch in BLAS library ',ch10,&
 308     'Some BLAS interfaces are not generated by abilint '
 309    ABI_ERROR(msg)
 310  end if
 311  write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10
 312 #else
 313  write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10
 314 #endif
 315  call wrtout([std_out, ab_out], msg)
 316 
 317  tol_empty=0.01                            ! Initialize the tolerance used to decide if a band is empty (passed to m_sigx.F90)
 318  gwcalctyp=Dtset%gwcalctyp
 319  gw1rdm=Dtset%gw1rdm                       ! Input variable to decide if updates to the 1-RDM must be performed
 320  x1rdm=Dtset%x1rdm                         ! Input variable to use pure exchange correction on the 1-RDM ( Sigma_x - Vxc )
 321  wfknocheck=.true.                         ! Used for printing WFK file subroutine
 322  rdm_update=(gwcalctyp==21 .and. gw1rdm>0) ! Input variable to decide whether to update GW density matrix
 323  readchkprdm=(Dtset%irdchkprdm==1)         ! Input variable to decide if checkpoint files must be read
 324  prtchkprdm=(Dtset%prtchkprdm==1)          ! Input variable to decide if checkpoint files must be written
 325 
 326  mod10 =MOD(Dtset%gwcalctyp,10)
 327 
 328  ! Perform some additional checks for hybrid functional calculations
 329  if (mod10 == 5) then
 330    !if (Dtset%ixc_sigma<0 .and. .not.libxc_functionals_check()) then
 331    !  ABI_ERROR('Hybrid functional calculations require the compilation with LIBXC library')
 332    !end if
 333    !XG 20171116 : I do not agree with this condition, as one might like to do a one-shot hybrid functional calculation
 334    !on top of a LDA/GGA calculation ... give the power (and risks) to the user !
 335    !if(gwcalctyp<10) then
 336    !  ABI_ERROR('gwcalctyp requires the update of energies and/or wavefunctions when performing hybrid XC calculations')
 337    !end if
 338    if (Dtset%usepaw == 1) then
 339      ABI_ERROR('PAW version of hybrid functional calculations is not implemented')
 340    end if
 341  end if
 342 
 343  !=== Initialize MPI variables, and parallelization level ===
 344  ! gwpara: 1--> parallelism over k-points, 2--> parallelism over bands.
 345  ! In case of gwpara==1 memory is not parallelized.
 346  ! If gwpara==2, bands are divided among processors but each proc has all the states where GW corrections are required.
 347  comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
 348 
 349  if (my_rank == master) then
 350    wfk_fname = dtfil%fnamewffk
 351    if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then
 352      ABI_ERROR(msg)
 353    end if
 354  end if
 355  call xmpi_bcast(wfk_fname, master, comm, ierr)
 356 
 357  ! === Some variables need to be initialized/nullify at start ===
 358  usexcnhat = 0
 359  call energies_init(KS_energies)
 360  call mkrdim(acell, rprim, rprimd)
 361  call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
 362  ucvol_local = ucvol
 363  !
 364  ! === Define FFT grid(s) sizes ===
 365  ! Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the
 366  ! (usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90.
 367  ! See also NOTES in the comments at the beginning of this file.
 368  ! NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn.
 369 
 370  call pawfgr_init(Pawfgr, Dtset, mgfftf, nfftf, ecut_eff, ecutdg_eff, ngfftc, ngfftf, &
 371                   gsqcutc_eff=gsqcutc_eff, gsqcutf_eff=gsqcutf_eff, gmet=gmet, k0=k0)
 372 
 373  ! Fake MPI_type for the sequential part.
 374  call initmpi_seq(MPI_enreg_seq)
 375  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all')
 376  call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all')
 377 
 378  call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials')
 379  nfftf_tot=PRODUCT(ngfftf(1:3))
 380 
 381  ! ===========================================
 382  ! === Open and read pseudopotential files ===
 383  ! ===========================================
 384  call pspini(dtset, dtfil, ecore, psp_gencond, gsqcutc_eff, gsqcutf_eff, pawrad, pawtab, psps, cryst%rprimd, comm_mpi=comm)
 385 
 386  call timab(402,2,tsec) ! Init1
 387  !
 388  ! ===============================================
 389  ! ==== Initialize Sigp, Er and basic objects ====
 390  ! ===============================================
 391  ! * Sigp is completetly initialized here.
 392  ! * Er is only initialized with dimensions, (SCR|SUSC) file is read in mkdump_Er
 393  call timab(403,1,tsec) ! setup_sigma
 394 
 395  call setup_sigma(codvsn,wfk_fname,acell,rprim,Dtset,Dtfil,Psps,Pawtab,&
 396   gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_sigma,Cryst,Kmesh,Qmesh,ks_ebands,Gsph_Max,Gsph_x,Gsph_c,Vcp,Er,Sigp,comm)
 397 
 398  call timab(403,2,tsec) ! setup_sigma
 399  call timab(402,1,tsec) ! Init1
 400 
 401  if (nprocs > Sigp%nbnds) then
 402    write(msg,"(2(a,i0))")"The number of MPI procs: ", nprocs, " is greater than nband: ", sigp%nbnds
 403    ABI_ERROR(msg)
 404  end if
 405 
 406  pole_screening = .FALSE.
 407  if (Er%fform==2002) then
 408    pole_screening = .TRUE.
 409    ABI_WARNING(' EXPERIMENTAL - Using a pole-fit screening!')
 410  end if
 411 
 412  call print_ngfft(gwc_ngfft,header='FFT mesh for oscillator strengths used for Sigma_c')
 413  call print_ngfft(gwx_ngfft,header='FFT mesh for oscillator strengths used for Sigma_x')
 414 
 415  b1gw=Sigp%minbdgw
 416  b2gw=Sigp%maxbdgw
 417 
 418  gwc_nfftot=PRODUCT(gwc_ngfft(1:3))
 419  gwc_nfft  =gwc_nfftot  !no FFT //
 420 
 421  gwx_nfftot=PRODUCT(gwx_ngfft(1:3))
 422  gwx_nfft  =gwx_nfftot  !no FFT //
 423 
 424  ! TRYING TO RECREATE AN "ABINIT ENVIRONMENT"
 425  KS_energies%e_corepsp=ecore/Cryst%ucvol
 426 
 427  ! === Calculate KS occupation numbers and ks_vbk(nkibz,nsppol) ====
 428  ! * ks_vbk gives the (valence|last Fermi band) index for each k and spin.
 429  ! * spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals
 430  ABI_MALLOC(ks_vbik, (ks_ebands%nkpt, ks_ebands%nsppol))
 431  ABI_MALLOC(qp_vbik, (ks_ebands%nkpt, ks_ebands%nsppol))
 432 
 433  !call ebands_update_occ(ks_ebands,Dtset%spinmagntarget,prtvol=0)
 434  ks_vbik(:,:) = ebands_get_valence_idx(ks_ebands)
 435 
 436  ! ============================
 437  ! ==== PAW initialization ====
 438  ! ============================
 439  if (dtset%usepaw == 1) then
 440    call chkpawovlp(cryst%natom, cryst%ntypat, dtset%pawovlp, pawtab, cryst%rmet, cryst%typat, cryst%xred)
 441 
 442    cplex_dij = dtset%nspinor; cplex = 1; ndij = 1
 443 
 444    ABI_MALLOC(ks_pawrhoij, (cryst%natom))
 445    call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij, nspden_rhoij=nspden_rhoij, &
 446                              nspden=dtset%nspden, spnorb=dtset%pawspnorb, cpxocc=Dtset%pawcpxocc)
 447    call pawrhoij_alloc(ks_pawrhoij, cplex_rhoij, nspden_rhoij, dtset%nspinor, dtset%nsppol, cryst%typat, pawtab=pawtab)
 448 
 449    ! Test if we have to call pawinit
 450    gnt_option = 1; if (dtset%pawxcdev == 2 .or. (dtset%pawxcdev == 1 .and. dtset%positron /= 0)) gnt_option = 2
 451    call paw_gencond(dtset, gnt_option, "test", call_pawinit)
 452 
 453    if (psp_gencond == 1 .or. call_pawinit) then
 454      call timab(553,1,tsec)
 455      gsqcut_shp = two * abs(dtset%diecut) * dtset%dilatmx**2 / pi**2
 456      call pawinit(dtset%effmass_free, gnt_option, gsqcut_shp, zero, dtset%pawlcutd, dtset%pawlmix, &
 457                   psps%mpsang, dtset%pawnphi, cryst%nsym, dtset%pawntheta, pawang, pawrad, &
 458                   dtset%pawspnorb, pawtab, dtset%pawxcdev, dtset%ixc, dtset%usepotzero)
 459      call timab(553,2,tsec)
 460 
 461      ! Update internal values
 462      call paw_gencond(dtset, gnt_option, "save", call_pawinit)
 463    else
 464      if (pawtab(1)%has_kij  ==1) pawtab(1:cryst%ntypat)%has_kij   = 2
 465      if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla = 2
 466    end if
 467 
 468    psps%n1xccc = MAXVAL(pawtab(1:cryst%ntypat)%usetcore)
 469 
 470    ! Initialize optional flags in Pawtab to zero
 471    ! Cannot be done in Pawinit since the routine is called only if some pars. are changed
 472    pawtab(:)%has_nabla = 0
 473    Pawtab(:)%lamb_shielding = zero
 474 
 475    call setsym_ylm(gprimd, pawang%l_max-1, cryst%nsym, dtset%pawprtvol, cryst%rprimd, cryst%symrec, pawang%zarot)
 476 
 477    ! Initialize and compute data for DFT+U
 478    Paw_dmft%use_dmft=Dtset%usedmft
 479    call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,&
 480      is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,dtset%nspinor,Cryst%ntypat,dtset%optdcmagpawu,Pawang,Dtset%pawprtvol,&
 481      Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa)
 482 
 483    if (my_rank == master) call pawtab_print(Pawtab)
 484 
 485    ! Get Pawrhoij from the header of the WFK file.
 486    call pawrhoij_copy(Hdr_wfk%pawrhoij, KS_Pawrhoij)
 487 
 488    ! Evaluate form factor of radial part of phi.phj-tphi.tphj.
 489    rhoxsp_method = 1  ! Arnaud-Alouani (default in sigma)
 490    !rhoxsp_method = 2 ! Shiskin-Kresse
 491    if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc
 492 
 493    ! The q-grid must contain the FFT mesh used for sigma_c and the G-sphere for the exchange part.
 494    ! We use the FFT mesh for sigma_c since COHSEX and the extrapolar method require oscillator
 495    ! strengths on the FFT mesh.
 496    ABI_MALLOC(tmp_gfft,(3, gwc_nfftot))
 497    call get_gftt(gwc_ngfft, k0, gmet, gwc_gsq, tmp_gfft)
 498    ABI_FREE(tmp_gfft)
 499 
 500    gwx_gsq = Dtset%ecutsigx/(two*pi**2)
 501    gw_gsq = MAX(gwx_gsq, gwc_gsq)
 502 
 503    ! Set up q-grid, make qmax 20% larger than largest expected.
 504    ABI_MALLOC(nq_spl, (Psps%ntypat))
 505    ABI_MALLOC(qmax, (Psps%ntypat))
 506    qmax = SQRT(gw_gsq)*1.2d0
 507    nq_spl = Psps%mqgrid_ff
 508    ! write(std_out,*)"using nq_spl",nq_spl,"qmax=",qmax
 509    ABI_MALLOC(Paw_pwff, (Psps%ntypat))
 510    call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps)
 511 
 512    ABI_FREE(nq_spl)
 513    ABI_FREE(qmax)
 514 
 515    ! Variables/arrays related to the fine FFT grid
 516    ABI_MALLOC(ks_nhat, (nfftf, Dtset%nspden))
 517    ks_nhat = zero
 518    ABI_MALLOC(Pawfgrtab, (Cryst%natom))
 519    call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
 520 
 521    cplex = 1
 522    call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat)
 523    ABI_FREE(l_size_atm)
 524    compch_fft=greatest_real
 525    usexcnhat = MAXVAL(Pawtab(:)%usexcnhat)
 526    ! * 0 if Vloc in atomic data is Vbare    (Blochl's formulation)
 527    ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse's formulation)
 528    call wrtout(std_out, sjoin(' using usexcnhat: ', itoa(usexcnhat)))
 529    !
 530    ! Identify parts of the rectangular grid where the density has to be calculated ===
 531    optcut = 0; optgr0 = Dtset%pawstgylm; optgr1 = 0; optgr2 = 0; optrad = 1 - Dtset%pawstgylm
 532    if (Dtset%pawcross==1) optrad=1
 533    if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm
 534 
 535    call nhatgrid(Cryst%atindx1, gmet, Cryst%natom, Cryst%natom, Cryst%nattyp, ngfftf, Cryst%ntypat,&
 536     optcut, optgr0, optgr1, optgr2, optrad, Pawfgrtab, Pawtab, Cryst%rprimd, Cryst%typat, Cryst%ucvol, Cryst%xred)
 537 
 538    call pawfgrtab_print(Pawfgrtab,Cryst%natom,unit=std_out,prtvol=Dtset%pawprtvol)
 539 
 540  else
 541    ABI_MALLOC(Paw_pwff, (0))
 542    ABI_MALLOC(Pawfgrtab, (0))
 543  end if ! End of PAW Initialization
 544 
 545  ! Consistency check and additional stuff done only for GW with PAW.
 546  ABI_MALLOC(Paw_onsite, (0))
 547 
 548  if (Dtset%usepaw == 1) then
 549    if (Dtset%ecutwfn < Dtset%ecut) then
 550      write(msg,"(3a)")&
 551      "  It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,&
 552      "  an excessive truncation of the planewave basis set can lead to unphysical results."
 553      ABI_WARNING(msg)
 554    end if
 555 
 556    ABI_CHECK(Dtset%useexexch == 0, "LEXX not yet implemented in GW")
 557    ABI_CHECK(Paw_dmft%use_dmft == 0, "DMFT + GW not available")
 558 
 559    ! Optionally read core orbitals from file and calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for the HF decoupling.
 560    if (Sigp%use_sigxcore == 1) then
 561      lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
 562      ABI_MALLOC(dijexc_core,(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat))
 563 
 564      call paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,Dtset%prtvol,Psps%filpsp)
 565    end if ! HF decoupling
 566 
 567    if (Dtset%pawcross==1) then
 568      ABI_SFREE(Paw_onsite)
 569      ABI_MALLOC(Paw_onsite,(Cryst%natom))
 570      call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat, &
 571                              Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab)
 572    end if
 573  end if
 574 
 575  ! Allocate these arrays anyway, since they are passed to subroutines.
 576  if (.not.allocated(ks_nhat)) then
 577    ABI_MALLOC(ks_nhat, (nfftf, 0))
 578  end if
 579  if (.not.allocated(dijexc_core)) then
 580    ABI_MALLOC(dijexc_core, (1, 1, 0))
 581  end if
 582 
 583  ! ==================================================
 584  ! ==== Read KS band structure from the KSS file ====
 585  ! ==================================================
 586  !
 587  ! Initialize Wfd, allocate wavefunctions and precalculate tables to do the FFT using the coarse gwc_ngfft.
 588  mband=Sigp%nbnds
 589  ABI_MALLOC(bks_mask, (mband, Kmesh%nibz, Sigp%nsppol))
 590  ABI_MALLOC(keep_ur , (mband, Kmesh%nibz, Sigp%nsppol))
 591  keep_ur=.FALSE.; bks_mask=.FALSE.
 592 
 593  if (rdm_update) then
 594    ABI_MALLOC(bdm_mask, (mband, Kmesh%nibz, Sigp%nsppol))
 595    bdm_mask=.FALSE.
 596    if (my_rank==master) then
 597      bdm_mask=.TRUE.
 598    end if
 599    ABI_MALLOC(bdm2_mask, (mband,Kmesh%nibz,Sigp%nsppol))
 600    bdm2_mask=.FALSE.
 601  end if
 602  ABI_MALLOC(nband,(Kmesh%nibz, Sigp%nsppol))
 603  nband=mband
 604 
 605  ! autoparal section
 606  if (dtset%max_ncpus/=0) then
 607    ount =ab_out
 608    ! Temporary table needed to estimate memory
 609    ABI_MALLOC(nlmn_atm,(Cryst%natom))
 610    if (Dtset%usepaw==1) then
 611      do iat=1,Cryst%natom
 612        nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size
 613      end do
 614    end if
 615 
 616    write(ount,'(a)')"--- !Autoparal"
 617    write(ount,"(a)")'#Autoparal section for Sigma runs.'
 618    write(ount,"(a)")   "info:"
 619    write(ount,"(a,i0)")"    autoparal: ",dtset%autoparal
 620    write(ount,"(a,i0)")"    max_ncpus: ",dtset%max_ncpus
 621    write(ount,"(a,i0)")"    gwpara: ",dtset%gwpara
 622    write(ount,"(a,i0)")"    nkpt: ",dtset%nkpt
 623    write(ount,"(a,i0)")"    nsppol: ",dtset%nsppol
 624    write(ount,"(a,i0)")"    nspinor: ",dtset%nspinor
 625    write(ount,"(a,i0)")"    nbnds: ",Sigp%nbnds
 626 
 627    work_size = mband * Kmesh%nibz * Sigp%nsppol
 628 
 629    ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI.
 630    nonscal_mem = (two*gwpc*Er%npwe**2*Er%nomega*(Er%mqmem+1)*b2Mb) * 1.1_dp
 631 
 632    ! List of configurations.
 633    ! Assuming an OpenMP implementation with perfect speedup!
 634    write(ount,"(a)")"configurations:"
 635    do ii=1,dtset%max_ncpus
 636      nstates_per_proc = 0
 637      eff = HUGE(one)
 638      max_wfsmem_mb = zero
 639 
 640      do my_rank=0,ii-1
 641        call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,ii,my_spins,bks_mask,keep_ur,ierr)
 642        ABI_FREE(my_spins)
 643        if (ierr /= 0) exit
 644        my_nbks = COUNT(bks_mask)
 645        nstates_per_proc = MAX(nstates_per_proc, my_nbks)
 646        eff = MIN(eff, (one * work_size) / (ii * nstates_per_proc))
 647 
 648        ! Memory needed for Fourier components ug.
 649        ug_mem = two*gwpc*Dtset%nspinor*Sigp%npwwfn*my_nbks*b2Mb
 650        ! Memory needed for real space ur (use gwc_nfft, instead of gwx_nfft)
 651        ur_mem = two*gwpc*Dtset%nspinor*gwc_nfft*COUNT(keep_ur)*b2Mb
 652        ! Memory needed for PAW projections Cprj
 653        cprj_mem = zero
 654        if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb
 655        max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem)
 656      end do
 657      if (ierr /= 0) cycle
 658 
 659      ! Add the non-scalable part and increase by 10% to account for other datastructures.
 660      mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp
 661      do omp_ncpus=1,xomp_get_max_threads()
 662        write(ount,"(a,i0)")"    - tot_ncpus: ",ii * omp_ncpus
 663        write(ount,"(a,i0)")"      mpi_ncpus: ",ii
 664        write(ount,"(a,i0)")"      omp_ncpus: ",omp_ncpus
 665        write(ount,"(a,f12.9)")"      efficiency: ",eff
 666        write(ount,"(a,f12.2)")"      mem_per_cpu: ",mempercpu_mb
 667      end do
 668    end do
 669    write(ount,'(a)')"..."
 670 
 671    ABI_FREE(nlmn_atm)
 672    ABI_ERROR_NODUMP("aborting now")
 673 
 674  else
 675    call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr)
 676    ABI_CHECK(ierr==0, "Error in sigma_bksmask")
 677  end if
 678 
 679  ! Each core stores the wavefunctions where GW corrections are required.
 680  do isp=1,SIZE(my_spins)
 681    spin = my_spins(isp)
 682    do ikcalc=1,Sigp%nkptgw
 683      ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW
 684      ii=Sigp%minbnd(ikcalc,spin); jj=Sigp%maxbnd(ikcalc,spin)
 685      bks_mask(ii:jj,ik_ibz,spin) = .TRUE.
 686      if (MODULO(Dtset%gwmem,10)==1) keep_ur(ii:jj,ik_ibz,spin)=.TRUE.
 687    end do
 688  end do
 689 
 690  ABI_FREE(my_spins)
 691 
 692  call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,&
 693    Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,&
 694    Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,use_fnl_dir0der0=(Dtset%gw1rdm==2))
 695 
 696  ! MRM: also initialize the Wfd_nato_master for GW 1-RDM if required.
 697  ! Warning, this should be replaced by copy but copy fails due to bands being allocated in different manners. Do it in the future! FIXME
 698  if (rdm_update) then
 699    bdm2_mask=bks_mask ! As bks_mask is going to be removed, save it in bdm2_mask to use it in Evext_nl
 700    call wfd_init(Wfd_nato_master,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bdm_mask,&
 701      Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,&
 702      Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,xmpi_comm_self)!comm)  ! MPI_COMM_SELF
 703    call Wfd_nato_master%read_wfk(wfk_fname,iomode_from_fname(wfk_fname))
 704  end if
 705 
 706  if (Dtset%pawcross==1) then
 707    call wfd_init(Wfdf,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,&
 708      Dtset%nspden,Dtset%nspinor,dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,&
 709      Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm)
 710  end if
 711 
 712  ABI_FREE(bks_mask)
 713  ABI_FREE(nband)
 714  ABI_FREE(keep_ur)
 715 
 716  call timab(402,2,tsec) ! sigma(Init1)
 717  call timab(404,1,tsec) ! rdkss
 718 
 719  call wfd%read_wfk(wfk_fname, iomode_from_fname(wfk_fname))
 720 
 721  if (Dtset%pawcross==1) then
 722    call wfdgw_copy(Wfd, Wfdf)
 723    call wfdf%change_ngfft(Cryst,Psps,ngfftf)
 724  end if
 725 
 726  ! This test has been disabled (too expensive!)
 727  if (.False.) call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL")
 728 
 729  call timab(404,2,tsec) ! rdkss
 730  call timab(405,1,tsec) ! Init2
 731 
 732  ! ==============================================================
 733  ! ==== Find little group of the k-points for GW corrections ====
 734  ! ==============================================================
 735  ! * The little group is used only if symsigma == 1
 736  ! * If use_umklp == 1 then symmetries requiring an umklapp to preserve k_gw are included as well.
 737  !
 738  ABI_MALLOC(Ltg_k, (Sigp%nkptgw))
 739  use_umklp = 1
 740  do ikcalc=1,Sigp%nkptgw
 741    if (Sigp%symsigma /= 0) then
 742      call Ltg_k(ikcalc)%init(Sigp%kptgw(:,ikcalc), Qmesh%nbz, Qmesh%bz, Cryst, use_umklp, npwe=0)
 743    end if
 744  end do
 745 
 746  ! Compute structure factor phases and large sphere cut-off
 747  ABI_MALLOC(ph1d,(2, 3 * (2 * Dtset%mgfft + 1) * Cryst%natom))
 748  ABI_MALLOC(ph1df,(2, 3 * (2 * mgfftf + 1) * Cryst%natom))
 749 
 750  call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred)
 751 
 752  if (Psps%usepaw == 1.and. Pawfgr%usefinegrid == 1) then
 753    call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred)
 754  else
 755    ph1df(:,:)=ph1d(:,:)
 756  end if
 757 
 758  !===================================================================================
 759  !==== Classify the GW wavefunctions according to the irreducible representation ====
 760  !===================================================================================
 761  !* Warning still under development.
 762  !* Only for SCGW.
 763  !bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
 764 
 765  ABI_MALLOC(KS_sym,(Wfd%nkibz,Wfd%nsppol))
 766 
 767  if (Sigp%symsigma==1.and.gwcalctyp>=20) then
 768    ! call check_zarot(Gsph_c%ng,Cryst,gwc_ngfft,Gsph_c%gvec,Psps,Pawang,Gsph_c%rottb,Gsph_c%rottbm1)
 769    use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric
 770    do spin=1,Wfd%nsppol
 771      do ikcalc=1,Sigp%nkptgw
 772        ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
 773        first_band = Sigp%minbnd(ikcalc,spin)
 774        last_band  = Sigp%maxbnd(ikcalc,spin)
 775 !      call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,ks_ebands,Pawtab,Pawrad,Pawang,Psps,&
 776        call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,ks_ebands,Pawtab,Pawrad,Pawang,Psps,&
 777 &       Dtset%tolsym,KS_sym(ik_ibz,spin))
 778      end do
 779    end do
 780    ! Recreate the Sig_ij tables taking advantage of the classification of the bands.
 781    call sigma_tables(Sigp,Kmesh, esymm=KS_sym)
 782  end if
 783 
 784  call timab(405,2,tsec) ! Init2
 785  call timab(406,1,tsec) ! make_vhxc
 786 
 787  !===========================
 788  !=== COMPUTE THE DENSITY ===
 789  !===========================
 790  ! Evaluate the planewave part (complete charge in case of NC pseudos).
 791  ABI_MALLOC(ks_rhor, (nfftf, dtset%nspden))
 792  ABI_MALLOC(ks_taur, (nfftf, dtset%nspden * dtset%usekden))
 793 
 794  call wfd%mkrho(Cryst, Psps, Kmesh, ks_ebands, ngfftf, nfftf, ks_rhor)
 795  if ((rdm_update .and. Dtset%prtden /= 0) .and. Wfd%my_rank == master) then
 796    ! Print initial (KS) density file as read (usefull to compare DEN files, cubes, etc.)
 797    gw1rdm_fname = trim(dtfil%fnameabo_ks_den)  ! and used on Sigma grids
 798    call fftdatar_write("density",gw1rdm_fname,dtset%iomode,hdr_sigma,&
 799                        Cryst,ngfftf,cplex1,nfftf,dtset%nspden,ks_rhor,mpi_enreg_seq,ebands=ks_ebands)
 800  end if
 801 
 802  if (Dtset%usekden == 1) call wfd%mkrho(Cryst,Psps,Kmesh,ks_ebands,ngfftf,nfftf,ks_taur,optcalc=1)
 803 
 804  !========================================
 805  !==== Additional computation for PAW ====
 806  !========================================
 807  nhatgrdim = 0
 808  if (Dtset%usepaw==1) then
 809    ! Calculate the compensation charge nhat.
 810    if (Dtset%xclevel==2) nhatgrdim = usexcnhat * Dtset%pawnhatxc
 811    cplex = 1; ider = 2 * nhatgrdim; izero = 0
 812    if (nhatgrdim > 0) then
 813      ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim))
 814    end if
 815    if (nhatgrdim == 0) then
 816      ABI_MALLOC(ks_nhatgr,(0,0,0))
 817    end if
 818 
 819    call pawmknhat(compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,&
 820                   Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
 821                   Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,k0,Cryst%rprimd,&
 822                   Cryst%ucvol,dtset%usewvl,Cryst%xred)
 823 
 824    ! === Evaluate onsite energies, potentials, densities ===
 825    !  * Initialize variables/arrays related to the PAW spheres.
 826    !  * Initialize also lmselect (index of non-zero LM-moments of densities).
 827    ABI_MALLOC(KS_paw_ij, (Cryst%natom))
 828    has_dijso = Dtset%pawspnorb; has_dijU = merge(0, 1, Dtset%usepawu == 0)
 829 
 830    call paw_ij_nullify(KS_paw_ij)
 831    call paw_ij_init(KS_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,&
 832      Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
 833      has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,&
 834      has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1)
 835 
 836    nkxc1 = 0
 837    ABI_MALLOC(KS_paw_an, (Cryst%natom))
 838    call paw_an_nullify(KS_paw_an)
 839    call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
 840      cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1)
 841 
 842    !  Calculate onsite vxc with and without core charge.
 843    nzlmopt=-1; option=0; compch_sph=greatest_real
 844    call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert0,&
 845      Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
 846      Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,&
 847      Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,&
 848      Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp)
 849 
 850  else
 851    ABI_MALLOC(ks_nhatgr, (0, 0, 0))
 852    ABI_MALLOC(KS_paw_ij, (0))
 853    ABI_MALLOC(KS_paw_an, (0))
 854  end if ! PAW
 855 
 856  call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,&
 857                   Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
 858 
 859  ! For PAW, add the compensation charge on the FFT mesh, then get rho(G).
 860  if (Dtset%usepaw==1) ks_rhor = ks_rhor + ks_nhat
 861 
 862  call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_rhor,ucvol=ucvol)
 863  if (Dtset%usekden==1) call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=ucvol)
 864 
 865  ABI_MALLOC(ks_rhog, (2, nfftf))
 866  call fourdp(1, ks_rhog, ks_rhor(:,1),-1, MPI_enreg_seq, nfftf, 1, ngfftf, tim_fourdp5)
 867 
 868  ! The following steps have been gathered in the setvtr routine:
 869  !  - get Ewald energy and Ewald forces
 870  !  - compute local ionic pseudopotential vpsp
 871  !  - eventually compute 3D core electron density xccc3d
 872  !  - eventually compute vxc and vhartr
 873  !  - set up ks_vtrial
 874  !
 875  !*******************************************************************
 876  !**** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
 877  !*******************************************************************
 878 
 879  ngrvdw = 0
 880  ABI_MALLOC(grvdw, (3, ngrvdw))
 881  ABI_MALLOC(grchempottn, (3, Cryst%natom))
 882  ABI_MALLOC(grewtn, (3, Cryst%natom))
 883  nkxc = 0
 884  if (Dtset%nspden == 1) nkxc = 2
 885  if (Dtset%nspden >= 2) nkxc = 3 ! check GGA and spinor, quite a messy part!!!
 886  ! In case of MGGA, fxc and kxc are not available and we dont need them in sigma (for now ...)
 887  if (Dtset%ixc < 0 .and. libxc_functionals_ismgga()) nkxc = 0
 888  if (nkxc /= 0) then
 889    ABI_MALLOC(kxc, (nfftf, nkxc))
 890  end if
 891 
 892  n3xccc = 0; if (Psps%n1xccc /= 0) n3xccc = nfftf
 893  ABI_MALLOC(xccc3d, (n3xccc))
 894  ABI_MALLOC(ks_vhartr, (nfftf))
 895  ABI_MALLOC(ks_vtrial, (nfftf, Dtset%nspden))
 896  ABI_MALLOC(vpsp, (nfftf))
 897  ABI_MALLOC(ks_vxc, (nfftf, Dtset%nspden))
 898 
 899  optene = 4; moved_atm_inside = 0; moved_rhor = 0; istep = 1
 900 
 901  call setvtr(Cryst%atindx1,Dtset,KS_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
 902              istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
 903              Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
 904              optene,pawrad,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,Cryst%rprimd,strsxc,&
 905              Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,Wvl,xccc3d,Cryst%xred,taur=ks_taur)
 906 
 907   !============================
 908   !==== Compute KS PAW Dij ====
 909   !============================
 910  if (Dtset%usepaw == 1) then
 911    call timab(561,1,tsec)
 912 
 913    ! Calculate the unsymmetrized Dij.
 914    call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,&
 915                Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
 916                Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,&
 917                Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
 918                k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,ks_vxc,Cryst%xred,&
 919                nucdipmom=Dtset%nucdipmom)
 920 
 921    ! Symmetrize KS Dij
 922    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,&
 923                    Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,KS_paw_ij,Pawang,&
 924                    Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
 925 
 926    ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
 927    call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab)
 928    call timab(561,2,tsec)
 929  end if
 930 
 931  call timab(406,2,tsec) ! make_vhxc
 932 
 933  !=== Calculate Vxc(b1,b2,k,s)=<b1,k,s|v_{xc}|b2,k,s> for all the states included in GW ===
 934  !  * This part is parallelized within wfd%comm since each node has all GW wavefunctions.
 935  !  * Note that vH matrix elements are calculated using the true uncutted interaction.
 936  call timab(407,1,tsec) ! vHxc_me
 937 
 938  call KS_mflags%reset()
 939  if (rdm_update) then
 940    KS_mflags%has_hbare=1
 941    KS_mflags%has_kinetic=1
 942  end if
 943  KS_mflags%has_vhartree=1
 944  KS_mflags%has_vxc     =1
 945  KS_mflags%has_vxcval  =1
 946  if (Dtset%usepawu /= 0  )  KS_mflags%has_vu     = 1
 947  if (Dtset%useexexch /= 0)  KS_mflags%has_lexexch= 1
 948  if (Sigp%use_sigxcore == 1) KS_mflags%has_sxcore = 1
 949  if (gwcalctyp<10        )  KS_mflags%only_diago = 1 ! off-diagonal elements only for SC on wavefunctions.
 950 
 951  if (.FALSE.) then ! quick and dirty hack to test HF contribution.
 952    ABI_WARNING("testing on-site HF")
 953    lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
 954    ABI_MALLOC(dij_hf,(cplex_dij*lmn2_size_max,ndij,Cryst%natom))
 955    call paw_dijhf(ndij,cplex_dij,1,lmn2_size_max,Cryst%natom,Cryst%ntypat,Pawtab,Pawrad,Pawang,&
 956                   KS_Pawrhoij,dij_hf,Dtset%prtvol)
 957 
 958    do iat=1,Cryst%natom
 959      itypat = Cryst%typat(iat)
 960      ii = Pawtab(itypat)%lmn2_size
 961      KS_Paw_ij(iat)%dijxc(:,:) = dij_hf(1:cplex_dij*ii,:,iat)
 962      KS_Paw_ij(iat)%dijxc_hat(:,:) = zero
 963    end do
 964    ABI_FREE(dij_hf)
 965 
 966    !option_dij=3
 967    !call symdij(Cryst%gprimd,Cryst%indsym,ipert0,&
 968    !&   Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
 969    !&   KS_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,&
 970    !&   Cryst%symafm,Cryst%symrec)
 971  end if
 972 
 973  ABI_MALLOC(tmp_kstab, (2, Wfd%nkibz, Wfd%nsppol))
 974  tmp_kstab=0
 975  do spin=1,Sigp%nsppol
 976    do ikcalc=1,Sigp%nkptgw ! No spin dependent!
 977      ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
 978      tmp_kstab(1, ik_ibz, spin) = Sigp%minbnd(ikcalc, spin)
 979      tmp_kstab(2, ik_ibz, spin) = Sigp%maxbnd(ikcalc, spin)
 980    end do
 981  end do
 982 
 983  call calc_vhxc_me(Wfd, KS_mflags, KS_me, Cryst, Dtset, nfftf, ngfftf, &
 984                    ks_vtrial, ks_vhartr, ks_vxc, Psps, Pawtab, KS_paw_an, Pawang, Pawfgrtab, KS_paw_ij, dijexc_core, &
 985                    ks_rhor, usexcnhat, ks_nhat, ks_nhatgr, nhatgrdim, tmp_kstab, taur=ks_taur)
 986  ABI_FREE(tmp_kstab)
 987 
 988 !#ifdef DEV_HAVE_SCGW_SYM
 989 !Set KS matrix elements connecting different irreps to zero. Do not touch unknown bands!.
 990  if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then
 991    bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
 992    ABI_MALLOC(ks_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
 993    ks_irreptab=0
 994    do spin=1,Sigp%nsppol
 995      do ikcalc=1,Sigp%nkptgw
 996        ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
 997        first_band = Sigp%minbnd(ikcalc,spin)
 998        last_band  = Sigp%maxbnd(ikcalc,spin)
 999        if (.not.esymm_failed(KS_sym(ik_ibz,spin))) then
1000          ks_irreptab(first_band:last_band,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(first_band:last_band)
1001          !ks_irreptab(bmin:bmax,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(bmin:bmax)
1002        end if
1003      end do
1004    end do
1005    call KS_me%zero(ks_irreptab)
1006    ABI_FREE(ks_irreptab)
1007  end if
1008 !#endif
1009 
1010  call KS_me%print(header="Matrix elements in the KS basis set", prtvol=Dtset%prtvol)
1011  !
1012  ! If possible, calculate the EXX energy between the frozen core
1013  ! and the valence electrons using KS wavefunctions
1014  ! MG: Be careful here, since ex_energy is meaningful only if all occupied states are calculated.
1015  if( KS_mflags%has_sxcore ==1 ) then
1016    ! TODO
1017    !ex_energy = mels_get_exene_core(KS_me,kmesh,ks_ebands)
1018    ex_energy=zero
1019    do spin=1,Sigp%nsppol
1020      do ik=1,Kmesh%nibz
1021        do ib=b1gw,b2gw
1022          if (Sigp%nsig_ab==1) then
1023            ex_energy = ex_energy + half*ks_ebands%occ(ib,ik,spin)*Kmesh%wt(ik)*KS_me%sxcore(ib,ib,ik,spin)
1024          else
1025            ex_energy = ex_energy + half*ks_ebands%occ(ib,ik,spin)*Kmesh%wt(ik)*SUM(KS_me%sxcore(ib,ib,ik,:))
1026          end if
1027        end do
1028      end do
1029    end do
1030    write(msg,'(a,2(es16.6,a))')' CORE Exchange energy with KS wavefunctions: ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV'
1031    call wrtout(std_out, msg)
1032  end if
1033 
1034  call timab(407,2,tsec) ! vHxc_me
1035  call timab(408,1,tsec) ! hqp_init
1036 
1037  ! Do not break this coding!
1038  ! When gwcalctyp>10, the order of the bands can be interexchanged after
1039  ! the diagonalization. Therefore, we have to correctly assign the matrix elements to the corresponding
1040  ! bands and we cannot skip the following even though it looks unuseful.
1041  if (gwcalctyp >= 10) call wrtout(std_out, ch10//' *************** KS Energies *******************')
1042 
1043  !=== qp_ebands stores energies and occ. used for the calculation ===
1044  !  * Initialize qp_ebands with KS values.
1045  !  * In case of SC update qp_ebands using the QPS file.
1046  call ebands_copy(ks_ebands, qp_ebands)
1047 
1048  ABI_MALLOC(qp_rhor, (nfftf, Dtset%nspden))
1049  ABI_MALLOC(qp_taur, (nfftf, Dtset%nspden * Dtset%usekden))
1050  QP_sym => KS_sym
1051 
1052  if (gwcalctyp<10) then
1053    ! one-shot GW, just do a copy of the KS density.
1054    qp_rhor=ks_rhor
1055    if(Dtset%usekden==1)qp_taur=ks_taur
1056    QP_sym => KS_sym
1057  else
1058    ! Self-consistent GW.
1059    ! Read the unitary matrix and the QP energies of the previous step from the QPS file.
1060    call energies_init(QP_energies)
1061    QP_energies%e_corepsp=ecore/Cryst%ucvol
1062 
1063    ! m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}>
1064    ! Initialize the QP amplitudes with KS wavefunctions.
1065    ABI_MALLOC(Sr%m_ks_to_qp, (Sigp%nbnds, Sigp%nbnds, Kmesh%nibz, Sigp%nsppol))
1066    Sr%m_ks_to_qp=czero
1067    do ib=1,Sigp%nbnds
1068      Sr%m_ks_to_qp(ib,ib,:,:)=cone
1069    end do
1070 
1071    ! Now read m_ks_to_qp and update the energies in qp_ebands.
1072    ! TODO switch on the renormalization of n in sigma.
1073    ABI_MALLOC(prev_rhor, (nfftf, Dtset%nspden))
1074    ABI_MALLOC(prev_taur, (nfftf, Dtset%nspden*Dtset%usekden))
1075    ABI_MALLOC(prev_Pawrhoij, (Cryst%natom*Psps%usepaw))
1076 
1077    call rdqps(qp_ebands,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,&
1078               nfftf,ngfftf,Cryst%ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,Sr%m_ks_to_qp,prev_rhor,prev_Pawrhoij)
1079 
1080    !Find the irreps associated to the QP amplitudes starting from the analogous table for the KS states.
1081    !bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1082    !allocate(qp_irreptab(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
1083    !qp_irreptab=0
1084    !!qp_irreptab=ks_irreptab
1085 
1086    !do jb_qp=bmin,bmax
1087    !do ib_ks=bmin,bmax
1088    !if (ABS(Sr%m_ks_to_qp(ib_ks,jb_qp,ik_ibz,spin)) > tol12) then ! jb_qp has same the same character as ib_ks.
1089    !ks_irr = ks_irreptab(ib_ks,ib_ks,ik_ibz,spin)
1090    !qp_irreptab(jb_qp,jb_qp,ik_ibz,spin) = ks_irr
1091    !do ii=bmin,bmax
1092    !if (ks_irr == ks_irreptab(ii,ib_ks,ik_ibz,spin)) then
1093    !qp_irreptab(jb_qp,ii,ik_ibz,spin) = ks_irr
1094    !end if
1095    !end do
1096    !end if
1097    !end do
1098    !end do
1099 
1100    if (nscf==0) prev_rhor=ks_rhor
1101    if (nscf==0 .and. Dtset%usekden==1) prev_taur=ks_taur
1102 
1103    if (nscf>0.and.gwcalctyp>=20.and. wfd%my_rank == master) then
1104      ! Print the unitary transformation on std_out.
1105      call show_QP(qp_ebands,Sr%m_ks_to_qp,fromb=Sigp%minbdgw,tob=Sigp%maxbdgw,unit=std_out,tolmat=0.001_dp)
1106    end if
1107 
1108    ! Compute QP wfg as linear combination of KS states ===
1109    !  * Wfd%ug is modified inside calc_wf_qp
1110    !  * For PAW, update also the on-site projections.
1111    !  * WARNING the first dimension of MPI_enreg MUST be Kmesh%nibz
1112    !  TODO here we should use nbsc instead of nbnds
1113 
1114    call wfd%rotate(Cryst,Sr%m_ks_to_qp)
1115 
1116    ! Reinit the storage mode of Wfd as ug have been changed ===
1117    ! Update also the wavefunctions for GW corrections on each processor
1118    call wfd%reset_ur_cprj()
1119 
1120    ! This test has been disabled (too expensive!)
1121    if (.False.) call wfd%test_ortho(Cryst, Pawtab, unit=std_out)
1122 
1123    ! Compute QP occupation numbers.
1124    call wrtout(std_out,'sigma: calculating QP occupation numbers:')
1125    call ebands_update_occ(qp_ebands, Dtset%spinmagntarget, prtvol=0)
1126    qp_vbik(:,:) = ebands_get_valence_idx(qp_ebands)
1127 
1128 !  #ifdef DEV_HAVE_SCGW_SYM
1129    ! Calculate the irreducible representations of the new QP amplitdues.
1130    if (Sigp%symsigma==1.and.gwcalctyp>=20) then
1131      ABI_MALLOC(QP_sym,(Wfd%nkibz,Wfd%nsppol))
1132      use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric
1133      do spin=1,Wfd%nsppol
1134        do ikcalc=1,Sigp%nkptgw
1135          ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1136 !        Quick fix for SCGW+symm TODO fix properly!
1137          first_band = Sigp%minbnd(ikcalc,spin)
1138          last_band  = Sigp%maxbnd(ikcalc,spin)
1139 !        first_band = MINVAL(Sigp%minbnd(:,spin))
1140 !        last_band  = MAXVAL(Sigp%maxbnd(:,spin))
1141 !        call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,qp_ebands,Pawtab,Pawrad,Pawang,Psps,&
1142          call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,qp_ebands,Pawtab,Pawrad,Pawang,Psps,&
1143 &         Dtset%tolsym,QP_sym(ik_ibz,spin))
1144        end do
1145      end do
1146 
1147      ! Recreate the Sig_ij tables taking advantage of the classification of the bands.
1148      call sigma_tables(Sigp, Kmesh, esymm=QP_sym)
1149    end if
1150 !  #endif
1151 
1152    ! Compute QP density using the updated wfg.
1153    call wfd%mkrho(Cryst,Psps,Kmesh,qp_ebands,ngfftf,nfftf,qp_rhor)
1154    if (Dtset%usekden==1) call wfd%mkrho(Cryst, Psps, Kmesh, qp_ebands, ngfftf, nfftf, qp_taur, optcalc=1)
1155 
1156    ! ========================================
1157    ! ==== QP self-consistent GW with PAW ====
1158    ! ========================================
1159    if (Dtset%usepaw==1) then
1160      ABI_MALLOC(qp_nhat,(nfftf,Dtset%nspden))
1161      nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat
1162      ABI_MALLOC(qp_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim))
1163 
1164      ABI_MALLOC(QP_pawrhoij,(Cryst%natom))
1165      ABI_MALLOC(QP_paw_ij,(Cryst%natom))
1166      ABI_MALLOC(QP_paw_an,(Cryst%natom))
1167 
1168      ! Calculate new QP quantities: nhat, nhatgr, rho_ij, paw_ij, and paw_an.
1169      call paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,qp_ebands,&
1170        Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij,&
1171        QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim,qp_nhatgr,compch_sph,compch_fft)
1172    else
1173      ABI_MALLOC(qp_nhat, (0, 0))
1174      ABI_MALLOC(qp_nhatgr, (0, 0, 0))
1175      ABI_MALLOC(QP_pawrhoij, (0))
1176      ABI_MALLOC(QP_paw_ij, (0))
1177      ABI_MALLOC(QP_paw_an, (0))
1178    end if
1179 
1180    ! here I should renormalize the density
1181    call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,qp_rhor,Cryst%ucvol,&
1182                     Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
1183 
1184    if (Dtset%usepaw==1) qp_rhor(:,:)=qp_rhor(:,:)+qp_nhat(:,:) ! Add the "hat" term.
1185 
1186    call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_rhor,ucvol=ucvol)
1187    if(Dtset%usekden==1) then
1188      call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_taur,optrhor=1,ucvol=ucvol)
1189    end if
1190 
1191    ! Simple mixing of the PW density to damp oscillations in the Hartree potential.
1192    if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then
1193      write(msg,'(2a,f6.3)')ch10,' sigma: mixing QP densities using rhoqpmix= ',Dtset%rhoqpmix
1194      call wrtout(std_out, msg)
1195      qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor)
1196      if(Dtset%usekden==1) qp_taur = prev_taur + Dtset%rhoqpmix*(qp_taur-prev_taur) ! mix taur.
1197    end if
1198 
1199    ABI_FREE(prev_rhor)
1200    ABI_FREE(prev_taur)
1201    if (Psps%usepaw==1.and.nscf>0) then
1202      call pawrhoij_free(prev_pawrhoij)
1203    end if
1204    ABI_FREE(prev_pawrhoij)
1205 
1206    ABI_MALLOC(qp_rhog,(2,nfftf))
1207    call fourdp(1,qp_rhog,qp_rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,tim_fourdp5)
1208 
1209    ! ===========================================
1210    ! ==== Optional output of the QP density ====
1211    ! ===========================================
1212    if (Dtset%prtden/=0 .and. wfd%my_rank == master) then
1213      call fftdatar_write("qp_rhor",dtfil%fnameabo_qp_den,dtset%iomode,hdr_sigma,&
1214                          cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor,mpi_enreg_seq,ebands=qp_ebands)
1215    end if
1216 
1217    ! ===========================================
1218    ! === Optional output of the full QP density
1219    ! ===========================================
1220    if (Wfd%usepaw==1.and.Dtset%prtden==2) then
1221      ABI_MALLOC(qp_rhor_paw   ,(nfftf,Wfd%nspden))
1222      ABI_MALLOC(qp_rhor_n_one ,(nfftf,Wfd%nspden))
1223      ABI_MALLOC(qp_rhor_nt_one,(nfftf,Wfd%nspden))
1224 
1225      call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,qp_nhat,&
1226                  Wfd%nspinor,Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,QP_pawrhoij,Pawtab,Dtset%prtvol,&
1227                  qp_rhor,qp_rhor_paw,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
1228 
1229      ABI_FREE(qp_rhor_n_one)
1230      ABI_FREE(qp_rhor_nt_one)
1231 
1232      if (Dtset%prtvol > 9) then
1233        ! Print a normalisation check
1234        norm = SUM(qp_rhor_paw(:,1))*Cryst%ucvol/PRODUCT(Pawfgr%ngfft(1:3))
1235        write(msg,'(a,F8.4)') '  QUASIPARTICLE DENSITY CALCULATED - NORM OF DENSITY: ',norm
1236        call wrtout(std_out, msg)
1237      end if
1238 
1239      ! Write the density to file
1240      if (my_rank==master) then
1241        call fftdatar_write("qp_pawrhor",dtfil%fnameabo_qp_pawden,dtset%iomode,hdr_sigma,&
1242                            cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor_paw,mpi_enreg_seq,ebands=qp_ebands)
1243      end if
1244      ABI_FREE(qp_rhor_paw)
1245    end if
1246 
1247    nkxc=0
1248    if (Dtset%nspden==1) nkxc=2
1249    if (Dtset%nspden>=2) nkxc=3 !check GGA and spinor that is messy !!!
1250    !In case of MGGA, fxc and kxc are not available and we dont need them
1251    !for the screening part (for now ...)
1252    if (Dtset%ixc<0.and.libxc_functionals_ismgga()) nkxc=0
1253    if (nkxc/=0)  then
1254      ABI_MALLOC(qp_kxc,(nfftf,nkxc))
1255    end if
1256 
1257    ! **** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
1258    n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf
1259    ABI_MALLOC(qp_vhartr,(nfftf))
1260    ABI_MALLOC(qp_vtrial,(nfftf,Dtset%nspden))
1261    ABI_MALLOC(qp_vxc,(nfftf,Dtset%nspden))
1262 
1263    optene=4; moved_atm_inside=0; moved_rhor=0; istep=1
1264 
1265    call setvtr(Cryst%atindx1,Dtset,QP_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
1266                istep,qp_kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
1267                Cryst%nattyp,nfftf,ngfftf,ngrvdw,qp_nhat,qp_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
1268                optene,pawrad,Pawtab,ph1df,Psps,qp_rhog,qp_rhor,Cryst%rmet,Cryst%rprimd,strsxc,&
1269                Cryst%ucvol,usexcnhat,qp_vhartr,vpsp,qp_vtrial,qp_vxc,vxcavg_qp,Wvl,&
1270                xccc3d,Cryst%xred,taur=qp_taur)
1271 
1272    ABI_SFREE(qp_kxc)
1273 
1274    if (Dtset%usepaw==1) then
1275      call timab(561,1,tsec)
1276 
1277      ! Compute QP Dij
1278      call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,&
1279                  Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
1280                  Dtset%nspden,Cryst%ntypat,QP_paw_an,QP_paw_ij,Pawang,Pawfgrtab,&
1281                  Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
1282                  k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),qp_vtrial,qp_vxc,Cryst%xred,&
1283                  nucdipmom=Dtset%nucdipmom)
1284 
1285      ! Symmetrize total Dij
1286      option_dij=0
1287 
1288      call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,&
1289                      Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
1290                      QP_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,&
1291                      Cryst%symafm,Cryst%symrec)
1292 
1293      ! Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij.
1294      call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab)
1295      call timab(561,2,tsec)
1296    end if
1297 
1298    ehartree=half*SUM(qp_rhor(:,1)*qp_vhartr(:))/DBLE(nfftf)*Cryst%ucvol
1299 
1300    write(msg,'(a,80a)')ch10,('-',ii=1,80)
1301    call wrtout(ab_out,msg)
1302    write(msg,'(5a,f9.4,3a,es21.14,2a,es21.14)')ch10,&
1303     ' QP results after the unitary transformation in the KS subspace: ',ch10,ch10,&
1304     '  Number of electrons    = ',qp_rhog(1,1)*Cryst%ucvol,ch10,ch10,&
1305     '  QP Band energy    [Ha] = ',ebands_get_bandenergy(qp_ebands),ch10,&
1306     '  QP Hartree energy [Ha] = ',ehartree
1307    call wrtout(ab_out,msg)
1308    write(msg,'(a,80a)')ch10,('-',ii=1,80)
1309    call wrtout(ab_out,msg)
1310 
1311    ! TODO Since plasmonpole model 2-3-4 depend on the Fourier components of the density
1312    ! in case of self-consistency we might calculate here the ppm coefficients using qp_rhor
1313  end if ! gwcalctyp>=10
1314 
1315  ! KS hamiltonian: hdft(b1,b1,k,s)= <b1,k,s|H_s|b1,k,s>
1316  ABI_MALLOC(hdft, (b1gw:b2gw, b1gw:b2gw, Kmesh%nibz, Sigp%nsppol*Sigp%nsig_ab))
1317  hdft = czero
1318 
1319  if (Dtset%nspinor == 1) then
1320    do spin=1,Sigp%nsppol
1321      do ik=1,Kmesh%nibz
1322        do ib=b1gw,b2gw
1323          hdft(ib, ib, ik, spin) = ks_ebands%eig(ib, ik, spin)
1324        end do
1325      end do
1326    end do
1327  else
1328    ! Spinorial case
1329    !  * Note that here vxc contains the contribution of the core.
1330    !  * Scale ovlp if orthonormalization is not satisfied as npwwfn might be < npwvec.
1331    !  TODO add spin-orbit case and gwpara 2
1332    ABI_MALLOC(my_band_list, (wfd%mband))
1333    ABI_MALLOC(bmask, (wfd%mband))
1334    bmask = .False.; bmask(b1gw:b2gw) = .True.
1335 
1336    if (Wfd%usepaw == 1) then
1337      ABI_MALLOC(Cp1,(Wfd%natom, Wfd%nspinor))
1338      call pawcprj_alloc(Cp1, 0, Wfd%nlmn_atm)
1339    end if
1340 
1341    do spin=1,Sigp%nsppol
1342      do ik_ibz=1,Kmesh%nibz
1343        ! Distribute bands in [b1gw, b2gw] range
1344        call wfd%distribute_bands(ik_ibz, spin, my_nband, my_band_list, bmask=bmask)
1345        if (my_nband == 0) cycle
1346        npw_k = Wfd%npwarr(ik_ibz)
1347        do ii=1,my_nband  ! ib=b1gw,b2gw in sequential
1348          ib = my_band_list(ii)
1349          ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, spin, wave, msg) == 0, msg)
1350          ug1  => wave%ug
1351          cdummy = xdotc(npw_k*Wfd%nspinor,ug1,1,ug1,1)
1352          ovlp(1) = REAL(cdummy); ovlp(2) = AIMAG(cdummy)
1353 
1354          if (Psps%usepaw==1) then
1355            call wfd%get_cprj(ib,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.)
1356            ovlp = ovlp + paw_overlap(Cp1,Cp1,Cryst%typat,Pawtab)
1357          end if
1358          ! write(std_out,*)ovlp(1),ovlp(2)
1359          norm=DBLE(ovlp(1)+ovlp(2))
1360          ovlp(1)=DBLE(ovlp(1)/norm); ovlp(2)=DBLE(ovlp(2)/norm) ! ovlp(2)=cone-ovlp(1)
1361          hdft(ib,ib,ik_ibz,1) = ks_ebands%eig(ib,ik_ibz,1)*ovlp(1) - KS_me%vxc(ib,ib,ik_ibz,3)
1362          hdft(ib,ib,ik_ibz,2) = ks_ebands%eig(ib,ik_ibz,1)*ovlp(2) - KS_me%vxc(ib,ib,ik_ibz,4)
1363          hdft(ib,ib,ik_ibz,3) = KS_me%vxc(ib,ib,ik_ibz,3)
1364          hdft(ib,ib,ik_ibz,4) = KS_me%vxc(ib,ib,ik_ibz,4)
1365        end do
1366      end do
1367    end do
1368 
1369    call xmpi_sum(hdft, wfd%comm, ierr)
1370 
1371    ABI_FREE(my_band_list)
1372    ABI_FREE(bmask)
1373    if (Wfd%usepaw==1) then
1374      call pawcprj_free(Cp1)
1375      ABI_FREE(Cp1)
1376    end if
1377  end if
1378 
1379  ! Initialize Sigma results ===
1380  ! TODO it is better if we use ragged arrays indexed by the k-point
1381  call sigma_init(Sigp,Kmesh%nibz,Dtset%usepawu,Sr)
1382 
1383  ! Setup bare Hamiltonian := T + v_{loc} + v_{nl} + v_H.
1384  !
1385  ! * The representation depends wheter we are updating the wfs or not.
1386  ! * ks_vUme is zero unless we are using DFT+U as starting point, see calc_vHxc_braket
1387  ! * Note that vH matrix elements are calculated using the true uncutted interaction
1388  !   This should be changed if the cutoff is also used in the GS run.
1389 
1390  if (gwcalctyp < 10) then
1391    ! For one-shot GW use the KS representation.
1392    Sr%hhartree = hdft - KS_me%vxcval
1393 
1394    ! Additional stuff for PAW
1395    !  * DFT +U Hamiltonian
1396    !  * LEXX.
1397    !  * Core contribution estimated using Fock exchange.
1398    if (Dtset%usepaw==1) then
1399      if (Sigp%use_sigxcore == 1) Sr%hhartree = hdft - (KS_me%vxc - KS_me%sxcore)
1400      if (Dtset%usepawu /= 0) Sr%hhartree = Sr%hhartree - KS_me%vu
1401      if (Dtset%useexexch /= 0) then
1402        ABI_ERROR("useexexch > 0 not implemented")
1403        Sr%hhartree = Sr%hhartree - KS_me%vlexx
1404      end if
1405    end if
1406 
1407  else
1408    ! Self-consistent on energies and|or wavefunctions.
1409    !   * For NC get the bare Hamiltonian  $H_{bare}= T+v_{loc}+ v_{nl}$ in the KS representation
1410    !   * For PAW, calculate the matrix elements of h0, store also the new Dij in QP_Paw_ij.
1411    !   * h0 is defined as T+vH[tn+nhat+tnZc] + vxc[tnc] + dij_eff and
1412    !     dij_eff = dij^0 + dij^hartree + dij^xc-dij^xc_val + dijhat - dijhat_val.
1413    !     In the above expression tn, tnhat are QP quantities.
1414    if (Dtset%usepaw==0) then
1415      ABI_MALLOC(hbare, (b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab))
1416      hbare = hdft - KS_me%vhartree - KS_me%vxcval
1417 
1418      ! Change basis from KS to QP, hbare is overwritten: A_{QP} = U^\dagger A_{KS} U
1419      ABI_MALLOC(htmp, (b1gw:b2gw,b1gw:b2gw, Kmesh%nibz, Sigp%nsppol*Sigp%nsig_ab))
1420      ABI_MALLOC(ctmp, (b1gw:b2gw, b1gw:b2gw))
1421      ABI_MALLOC(uks2qp, (b1gw:b2gw, b1gw:b2gw))
1422 
1423      htmp = hbare; hbare = czero
1424 
1425      do spin=1,Sigp%nsppol
1426        do ik=1,Kmesh%nibz
1427          uks2qp(:,:) = Sr%m_ks_to_qp(b1gw:b2gw,b1gw:b2gw,ik,spin)
1428          do iab=1,Sigp%nsig_ab
1429            is_idx=spin; if (Sigp%nsig_ab>1) is_idx=iab
1430            ctmp(:,:)=MATMUL(htmp(:,:,ik,is_idx),uks2qp(:,:))
1431            hbare(:,:,ik,is_idx)=MATMUL(TRANSPOSE(CONJG(uks2qp)),ctmp)
1432          end do
1433        end do
1434      end do
1435      ABI_FREE(htmp)
1436      ABI_FREE(ctmp)
1437      ABI_FREE(uks2qp)
1438    end if ! usepaw==0
1439 
1440    ! Calculate the QP matrix elements
1441    ! This part is parallelized within MPI_COMM_WORD since each node has all GW wavefunctions.
1442    ! For PAW, we have to construct the new bare Hamiltonian.
1443    call wrtout(std_out,ch10//' *************** QP Energies *******************')
1444 
1445    call QP_mflags%reset()
1446    ! if (gwcalctyp<20) QP_mflags%only_diago=1 ! For e-only, no need of off-diagonal elements.
1447    QP_mflags%has_vhartree=1
1448    if (Dtset%usepaw==1) then
1449      QP_mflags%has_kinetic  =1
1450      QP_mflags%has_hbare    =1
1451    end if
1452    !QP_mflags%has_vxc     =1
1453    !QP_mflags%has_vxcval  =1
1454    !if (Sigp%gwcalctyp >100) QP_mflags%has_vxcval_hybrid=1
1455    if (mod10==5 .and. &
1456        (Dtset%ixc_sigma==-402 .or. Dtset%ixc_sigma==-406 .or. Dtset%ixc_sigma==-427 .or. Dtset%ixc_sigma==-428 .or. &
1457        Dtset%ixc_sigma==-456 .or. Dtset%ixc_sigma==41 .or. Dtset%ixc_sigma==42)) then
1458      QP_mflags%has_vxcval_hybrid = 1
1459    end if
1460 
1461    !if (Sigp%use_sigxcore==1) QP_mflags%has_sxcore =1
1462    !if (Dtset%usepawu/=0)    QP_mflags%has_vu     =1
1463    !if (Dtset%useexexch/=0)  QP_mflags%has_lexexch=1
1464 
1465    ABI_MALLOC(tmp_kstab, (2, Wfd%nkibz, Wfd%nsppol))
1466    tmp_kstab=0
1467    do spin=1,Sigp%nsppol
1468      do ikcalc=1,Sigp%nkptgw ! No spin dependent!
1469        ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1470        tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin)
1471        tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin)
1472      end do
1473    end do
1474 
1475    call calc_vhxc_me(Wfd,QP_mflags,QP_me,Cryst,Dtset,nfftf,ngfftf,&
1476                      qp_vtrial,qp_vhartr,qp_vxc,Psps,Pawtab,QP_paw_an,Pawang,Pawfgrtab,QP_paw_ij,dijexc_core,&
1477                      qp_rhor,usexcnhat,qp_nhat,qp_nhatgr,nhatgrdim,tmp_kstab,taur=qp_taur)
1478    ABI_FREE(tmp_kstab)
1479 
1480 !  #ifdef DEV_HAVE_SCGW_SYM
1481    if (gwcalctyp>=20 .and. Sigp%symsigma>0) then
1482      bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1483      ABI_MALLOC(qp_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
1484      qp_irreptab=0
1485      do spin=1,Sigp%nsppol
1486        do ikcalc=1,Sigp%nkptgw
1487          ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1488          first_band = Sigp%minbnd(ikcalc,spin)
1489          last_band  = Sigp%maxbnd(ikcalc,spin)
1490          if (.not.esymm_failed(QP_sym(ik_ibz,spin))) then
1491            qp_irreptab(first_band:last_band,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(first_band:last_band)
1492 !          qp_irreptab(bmin:bmax,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(bmin:bmax)
1493          end if
1494        end do
1495      end do
1496      call QP_me%zero(qp_irreptab)
1497      ABI_FREE(qp_irreptab)
1498    end if
1499 !  #endif
1500 
1501    call QP_me%print(header="Matrix elements in the QP basis set", prtvol=Dtset%prtvol)
1502 
1503    ! Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij.
1504    if (Dtset%usepaw==1) then
1505      call wrtout(std_out," *** After calc_vHxc_braket *** ")
1506      ! TODO terminate the implementation of this routine.
1507      call paw_ij_print(QP_Paw_ij,unit=std_out,pawprtvol=Dtset%pawprtvol,pawspnorb=Dtset%pawspnorb,mode_paral="COLL")
1508      call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab)
1509    end if
1510 
1511    if (Dtset%usepaw==0) then
1512      ! GA: We have an odd bug here. I have to unroll this loop, otherwise it
1513      ! might cause segfault when running on several nodes.
1514      !
1515      ! Sr%hhartree = hbare + QP_me%vhartree
1516      if(QP_mflags%has_vxcval_hybrid==0) then
1517        do spin=1,Sigp%nsppol*Sr%nsig_ab
1518          do ikcalc=1,Sr%nkibz
1519            do ib1=b1gw,b2gw
1520              do ib2=b1gw,b2gw
1521                Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + QP_me%vhartree(ib2,ib1,ikcalc,spin)
1522              end do
1523            end do
1524          end do
1525        end do
1526      else
1527        do spin=1,Sigp%nsppol*Sr%nsig_ab
1528          do ikcalc=1,Sr%nkibz
1529            do ib1=b1gw,b2gw
1530              do ib2=b1gw,b2gw
1531                Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + &
1532 &               QP_me%vhartree(ib2,ib1,ikcalc,spin) + &
1533 &               QP_me%vxcval_hybrid(ib2,ib1,ikcalc,spin)
1534              end do
1535            end do
1536          end do
1537        end do
1538      end if
1539    else
1540      Sr%hhartree=QP_me%hbare
1541    end if
1542 
1543 !  #ifdef DEV_HAVE_SCGW_SYM
1544    if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then
1545 !    bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1546      do spin=1,Sigp%nsppol
1547        do ik_ibz=1,Kmesh%nibz
1548          if (.not.esymm_failed(QP_sym(ik_ibz,spin))) then
1549            bmin=Sigp%minbnd(ik_ibz,spin); bmax=Sigp%minbnd(ik_ibz,spin)
1550            do ib2=bmin,bmax
1551              irr_idx2 = QP_sym(ik_ibz,spin)%b2irrep(ib2)
1552              do ib1=bmin,bmax
1553                irr_idx1 = QP_sym(ik_ibz,spin)%b2irrep(ib1)
1554                if (irr_idx1/=irr_idx2 .and. ALL((/irr_idx1,irr_idx2/)/=0) ) Sr%hhartree(ib1,ib2,ik_ibz,spin) = czero
1555              end do
1556            end do
1557          end if
1558        end do
1559      end do
1560    end if
1561 !  #endif
1562 
1563    ABI_FREE(qp_rhog)
1564    ABI_FREE(qp_vhartr)
1565    ABI_FREE(qp_vxc)
1566    ABI_FREE(qp_nhat)
1567    ABI_FREE(qp_nhatgr)
1568    call QP_me%free()
1569  end if ! gwcalctyp<10
1570 
1571  ! Free some memory
1572  ABI_SFREE(hbare)
1573  ABI_SFREE(hdft)
1574 
1575  ! Prepare the storage of QP amplitudes and energies
1576  ! Initialize with KS wavefunctions and energies.
1577  Sr%eigvec_qp=czero; Sr%en_qp_diago=zero
1578  do ib=1,Sigp%nbnds
1579    Sr%en_qp_diago(ib,:,:) = ks_ebands%eig(ib,:,:)
1580    Sr%eigvec_qp(ib,ib,:,:) = cone
1581  end do
1582 
1583  ! Store <n,k,s|V_xc[n_val]|n,k,s> and <n,k,s|V_U|n,k,s> ===
1584  ! Note that we store the matrix elements of V_xc in the KS basis set, not in the QP basis set
1585  ! Matrix elements of V_U are zero unless we are using DFT+U as starting point
1586  do ib=b1gw,b2gw
1587    Sr%vxcme(ib,:,:)=KS_me%vxcval(ib,ib,:,:)
1588    if (Dtset%usepawu/=0) Sr%vUme (ib,:,:)=KS_me%vu(ib,ib,:,:)
1589  end do
1590 
1591  ! Initial guess for the GW energies
1592  ! Save the energies of the previous iteration.
1593  do spin=1,Sigp%nsppol
1594    do ik=1,Kmesh%nibz
1595      do ib=1,Sigp%nbnds
1596        Sr%e0 (ib,ik,spin) = qp_ebands%eig(ib,ik,spin)
1597        Sr%egw(ib,ik,spin) = qp_ebands%eig(ib,ik,spin)
1598      end do
1599      Sr%e0gap(ik,spin) = zero
1600      ks_iv = ks_vbik(ik, spin)
1601      if (Sigp%nbnds>=ks_iv+1) Sr%e0gap(ik,spin)=Sr%e0(ks_iv+1,ik,spin)-Sr%e0(ks_iv,ik,spin)
1602    end do
1603  end do
1604  !
1605  !=== If required apply a scissor operator or update the energies ===
1606  !TODO check if other Sr entries have to be updated
1607  !moreover this part should be done only in case of semiconductors
1608  !FIXME To me it makes more sense if we apply the scissor to ks_ebands but I have to RECHECK csigme
1609  if (ABS(Sigp%mbpt_sciss)>tol6) then
1610    write(msg,'(6a,f10.5,2a)')ch10,&
1611    ' sigma : performing a first self-consistency',ch10,&
1612    '  update of the energies in G by a scissor operator',ch10, &
1613    '  applying a scissor operator of ',Sigp%mbpt_sciss*Ha_eV,' [eV] ',ch10
1614    call wrtout([std_out, ab_out], msg)
1615    do spin=1,Sigp%nsppol
1616      do ik=1,Kmesh%nibz
1617        ks_iv=ks_vbik(ik,spin)
1618        if (Sigp%nbnds>=ks_iv+1) then
1619          Sr%egw    (ks_iv+1:Sigp%nbnds,ik,spin) = Sr%egw    (ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss
1620          qp_ebands%eig(ks_iv+1:Sigp%nbnds,ik,spin) = qp_ebands%eig(ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss
1621        end if
1622      end do
1623    end do
1624    !call apply_scissor(qp_ebands,Sigp%mbpt_sciss)
1625  else if (.FALSE.) then
1626    write(msg,'(4a)')ch10,&
1627     ' sigma : performing a first self-consistency',ch10,&
1628     '  update of the energies in G by a previous GW calculation'
1629    call wrtout([std_out, ab_out], msg)
1630    ! TODO Recheck this part, is not clear to me!
1631    ABI_MALLOC(igwene,(qp_ebands%mband, qp_ebands%nkpt, qp_ebands%nsppol))
1632    call rdgw(qp_ebands, '__in.gw__', igwene, extrapolate=.TRUE.)
1633    ABI_FREE(igwene)
1634    Sr%egw=qp_ebands%eig
1635    !
1636    ! * Recalculate the new fermi level.
1637    call ebands_update_occ(qp_ebands, Dtset%spinmagntarget, prtvol=0)
1638  end if
1639 
1640  ! In case of AC refer all the energies wrt to the fermi level
1641  ! Be careful because results from ppmodel cannot be used for AC
1642  ! FIXME check ks_energy or qp_energy (in case of SCGW?)
1643 
1644  if (mod10 == SIG_GW_AC) then
1645    ! All these quantities will be passed to csigme
1646    ! if I skipped the self-consistent part then here I have to use fermi
1647    qp_ebands%eig = qp_ebands%eig - qp_ebands%fermie
1648    Sr%egw = Sr%egw - qp_ebands%fermie
1649    Sr%e0  = Sr%e0  - qp_ebands%fermie
1650    oldefermi = qp_ebands%fermie
1651    ! TODO Recheck fermi
1652    ! Clean EVERYTHING in particulare the treatment of E fermi
1653    qp_ebands%fermie=zero
1654  end if
1655  !
1656  ! === Setup frequencies around the KS\QP eigenvalues to compute Sigma derivatives (notice the spin) ===
1657  ! TODO it is better using an odd Sr%nomega4sd so that the KS\QP eigenvalue is in the middle
1658  ioe0j=Sr%nomega4sd/2+1
1659  do spin=1,Sigp%nsppol
1660    do ikcalc=1,Sigp%nkptgw
1661      ib1=Sigp%minbnd(ikcalc,spin)
1662      ib2=Sigp%maxbnd(ikcalc,spin)
1663      ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1664      do jb=ib1,ib2
1665        do io=1,Sr%nomega4sd
1666          Sr%omega4sd(jb,ik_ibz,io,spin)=Sr%egw(jb,ik_ibz,spin)+Sigp%deltae*(io-ioe0j)
1667        end do
1668      end do
1669    end do
1670  end do
1671 
1672  call timab(408,2,tsec) ! hqp_init
1673  call timab(409,1,tsec) ! getW
1674 
1675  ! Get epsilon^{-1} either from the _SCR or the _SUSC file and store it in Er%epsm1
1676  ! If Er%mqmem==0, allocate and read a single q-slice inside csigme.
1677  ! TODO Er%nomega should be initialized so that only the frequencies really needed are stored in memory
1678  ! TODO The same piece of code is present in screening.
1679  if (sigma_needs_w(Sigp)) then
1680    select case (dtset%gwgamma)
1681    case (0)
1682      id_required=4; ikxc=0; approx_type=0; option_test=0; dim_kxcg=0
1683      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1684 
1685    case (1, 2)
1686      ! ALDA TDDFT kernel vertex
1687      !ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available")
1688      ABI_CHECK(Er%ID==0,"Er%ID should be 0")
1689 
1690      if (Dtset%usepaw==1) then
1691        ! If we have PAW, we need the full density on the fine grid
1692        ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1693        if (Dtset%getpawden==0 .and. Dtset%irdpawden==0) then
1694          ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1695        end if
1696        call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin))
1697        if (.not. file_exists(dtfil%filpawdensin)) then
1698          ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1699        end if
1700 
1701        ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1702        call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1703                       ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm)
1704 
1705        call Hdr_rhor%free()
1706        call pawrhoij_free(tmp_pawrhoij)
1707        ABI_FREE(tmp_pawrhoij)
1708      end if ! Dtset%usepaw==1
1709 
1710      id_required=4; ikxc=7; approx_type=1; dim_kxcg=1
1711      if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1712      if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1713      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1714 
1715      dbg_mode=.FALSE.
1716      if (Dtset%usepaw==1) then
1717        ! Use PAW all-electron density
1718        call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_aepaw_rhor,&
1719        Er%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode)
1720      else
1721        ! Norm-conserving
1722        call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_rhor,&
1723        Er%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode)
1724      end if
1725 
1726    case (3, 4)
1727      ! ADA non-local kernel vertex
1728      !ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available")
1729      ABI_CHECK(Er%ID==0,"Er%ID should be 0")
1730      ABI_CHECK(Sigp%nsppol==1,"ADA vertex for GWGamma not available yet for spin-polarised cases")
1731 
1732      if (Dtset%usepaw==1) then ! If we have PAW, we need the full density on the fine grid
1733        ABI_MALLOC(ks_aepaw_rhor,(nfftf,Sigp%nsppol))
1734        if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then
1735          ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1736        end if
1737        call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin))
1738        if (.not. file_exists(dtfil%filpawdensin)) then
1739          ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1740        end if
1741 
1742        ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1743 
1744        call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1745                       ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm)
1746 
1747        call Hdr_rhor%free()
1748        call pawrhoij_free(tmp_pawrhoij)
1749        ABI_FREE(tmp_pawrhoij)
1750      end if ! Dtset%usepaw==1
1751 
1752      id_required=4; ikxc=7; approx_type=2;
1753      if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1754      if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1755      ABI_MALLOC(fxc_ADA,(Er%npwe,Er%npwe,Er%nqibz))
1756      ! Use userrd to set kappa
1757      if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp
1758      ! Set correct value of kappa (should be scaled with alpha*r_s where)
1759      ! r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3)
1760      rhoav = (drude_plsmf*drude_plsmf)/four_pi
1761      r_s = (three/(four_pi*rhoav))**third
1762      alpha = (four*ninth*piinv)**third
1763      Dtset%userrd = Dtset%userrd/(alpha*r_s)
1764 
1765      dbg_mode=.TRUE.
1766      if (Dtset%usepaw==1) then
1767        ! Use PAW all-electron density
1768        call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,&
1769        ks_aepaw_rhor,Er%npwe,Er%nqibz,Er%qibz,&
1770        fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode)
1771      else
1772        ! Norm conserving
1773        call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,&
1774        ks_rhor,Er%npwe,Er%nqibz,Er%qibz,&
1775        fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode)
1776      end if
1777 
1778      dim_kxcg = 0
1779      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1780 
1781    case (-3, -4, -5, -6, -7, -8)
1782      ! bootstrap kernels
1783      ! ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available")
1784      ABI_CHECK(Er%ID==0,"Er%ID should be 0")
1785 
1786      if (Dtset%usepaw==1) then
1787        ! If we have PAW, we need the full density on the fine grid
1788        ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1789        if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then
1790          ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1791        end if
1792        call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin))
1793        if (.not. file_exists(dtfil%filpawdensin)) then
1794          ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1795        end if
1796 
1797        ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1798 
1799        call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1800                       ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm)
1801 
1802        call Hdr_rhor%free()
1803        call pawrhoij_free(tmp_pawrhoij)
1804        ABI_FREE(tmp_pawrhoij)
1805      end if ! Dtset%usepaw==1
1806 
1807      id_required=4; ikxc=7; dim_kxcg=0
1808 
1809      if (dtset%gwgamma>-5) then
1810        approx_type=4  ! full fxc(G,G')
1811      else if (dtset%gwgamma>-7) then
1812        approx_type=5  ! fxc(0,0) one-shot
1813      else
1814        approx_type=6  ! rpa-type bootstrap
1815      end if
1816 
1817      option_test=MOD(Dtset%gwgamma,2)
1818      ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma
1819      ! 0 -> TESTPARTICLE, vertex in chi0 only
1820      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1821 
1822    case (-11)
1823      !LR+ALDA kernel
1824      !ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available")
1825      ABI_CHECK(Er%ID==0,"Er%ID should be 0")
1826 
1827      if (Dtset%usepaw==1) then
1828        ! If we have PAW, we need the full density on the fine grid
1829        ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1830        if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then
1831          ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1832        end if
1833        call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin))
1834        if (.not. file_exists(dtfil%filpawdensin)) then
1835          ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1836        end if
1837 
1838        ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1839 
1840        call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1841                       ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm)
1842 
1843        call hdr_rhor%free()
1844        call pawrhoij_free(tmp_pawrhoij)
1845        ABI_FREE(tmp_pawrhoij)
1846      end if ! Dtset%usepaw==1
1847 
1848      id_required=4; ikxc=7; dim_kxcg=1
1849      approx_type=7
1850      option_test=1
1851      ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma
1852      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1853 
1854    case default
1855      ABI_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma)))
1856    end select
1857 
1858    ! Set plasma frequency
1859    my_plsmf = drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf = Dtset%ppmfrq
1860    Dtset%ppmfrq = my_plsmf
1861 
1862    ! if(dtset%ucrpa==0) then
1863    if (Dtset%gwgamma<3) then
1864      call mkdump_Er(Er,Vcp,Er%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,&
1865                     approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,&
1866                     nfftf_tot,ngfftf,comm)
1867    else
1868      call mkdump_Er(Er,Vcp,Er%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,&
1869                     approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,&
1870                     nfftf_tot,ngfftf,comm,fxc_ADA)
1871    end if
1872 !   end if
1873    ABI_SFREE(kxcg)
1874    ABI_SFREE(fxc_ADA)
1875    ABI_SFREE(ks_aepaw_rhor)
1876  end if
1877  !
1878  ! ================================================
1879  ! ==== Calculate plasmonpole model parameters ====
1880  ! ================================================
1881  ! TODO Maybe its better if we use mqmem as input variable
1882  use_aerhor=0
1883  ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden*use_aerhor))
1884 
1885  if (sigma_needs_ppm(Sigp)) then
1886    my_plsmf=drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf=Dtset%ppmfrq
1887    call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,Sigp%ppmodel,my_plsmf,Dtset%gw_invalid_freq)
1888 
1889    ! PPm%force_plsmf= force_ppmfrq  ! this line to change the plasme frequency in HL expression.
1890 
1891    if (Wfd%usepaw==1 .and. Ppm%userho==1) then
1892      ! * For PAW and ppmodel 2-3-4 we need the AE rho(G) without compensation charge.
1893      ! * It would be possible to calculate rho(G) using Paw_pwff, though. It should be faster but
1894      !    results will depend on the expression used for the matrix elements. This approach is safer.
1895      use_aerhor=1
1896      ABI_FREE(ks_aepaw_rhor)
1897      ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1898 
1899      ! Check if input density file is available, otherwise compute
1900      pawden_fname = strcat(Dtfil%filnam_ds(3), '_PAWDEN')
1901      call wrtout(std_out,sjoin('Checking for existence of:',pawden_fname))
1902      if (file_exists(pawden_fname)) then
1903        ! Read density from file
1904        ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1905 
1906        call read_rhor(pawden_fname, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1907                       ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm)
1908 
1909        call Hdr_rhor%free()
1910        call pawrhoij_free(tmp_pawrhoij)
1911        ABI_FREE(tmp_pawrhoij)
1912      else
1913        ! Have to calculate PAW AW rhor from scratch
1914        ABI_MALLOC(qp_rhor_n_one,(pawfgr%nfft,Dtset%nspden))
1915        ABI_MALLOC(qp_rhor_nt_one,(pawfgr%nfft,Dtset%nspden))
1916 
1917        ! FIXME
1918        ABI_WARNING(" denfgr in sigma seems to produce wrong results")
1919        write(std_out,*)" input tilde ks_rhor integrates: ",SUM(ks_rhor(:,1))*Cryst%ucvol/nfftf
1920 
1921        call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,ks_nhat,Wfd%nspinor,&
1922        Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,KS_Pawrhoij,Pawtab,Dtset%prtvol, &
1923        ks_rhor,ks_aepaw_rhor,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
1924 
1925        ABI_FREE(qp_rhor_n_one)
1926        ABI_FREE(qp_rhor_nt_one)
1927      end if
1928 
1929      write(msg,'(a,f8.4)')' sigma: PAW AE density used for PPmodel integrates to: ',SUM(ks_aepaw_rhor(:,1))*Cryst%ucvol/nfftf
1930      call wrtout(std_out, msg)
1931 
1932      if (Er%mqmem/=0) then
1933        ! Calculate ppmodel parameters for all q-points.
1934        call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_aepaw_rhor(:,1))
1935      end if
1936 
1937    else
1938      ! NC or PAW with PPmodel 1.
1939      if (Er%mqmem/=0) then
1940        ! Calculate ppmodel parameters for all q-points
1941        call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_rhor(:,1))
1942      end if
1943    end if ! PAW or NC PPm and/or needs density
1944  end if ! sigma_needs_ppm
1945 
1946  call timab(409,2,tsec) ! getW
1947 
1948  if (wfd%my_rank == master) then
1949    ! Write info on the run on ab_out, then open files to store final results.
1950    call ebands_report_gap(ks_ebands,header='KS Band Gaps',unit=ab_out)
1951    if(dtset%ucrpa == 0) then
1952      call write_sigma_header(Sigp,Er,Cryst,Kmesh,Qmesh)
1953    end if
1954 
1955    ! unt_gw:  File with GW corrections.
1956    ! unt_sig: Self-energy as a function of frequency.
1957    ! unt_sgr: Derivative wrt omega of the Self-energy.
1958    ! unt_sigc: Sigma_c(eik) MRM
1959    ! unt_sgm: Sigma on the Matsubara axis (imag axis)
1960 
1961    if (open_file(Dtfil%fnameabo_gw, msg, unit=unt_gw, status='unknown', form='formatted') /= 0) then
1962      ABI_ERROR(msg)
1963    end if
1964    write(unt_gw,"(a)")"# QP energies E in eV"
1965    write(unt_gw,"(a)")"# Format:"
1966    write(unt_gw,"(a)")"#     kpoint"
1967    write(unt_gw,"(a)")"#     number of bands computed"
1968    write(unt_gw,"(2a)")"#     band index, Re(E), E-E0, Im(E)", ch10
1969    write(unt_gw,"(2(i0,1x),a)")Sigp%nkptgw,Sigp%nsppol, "# nkptgw, nsppol"
1970 
1971    if (open_file(Dtfil%fnameabo_gwdiag, msg, unit=unt_gwdiag, status='unknown', form='formatted') /= 0) then
1972      ABI_ERROR(msg)
1973    end if
1974    write(unt_gwdiag,*)Sigp%nkptgw,Sigp%nsppol
1975 
1976    if (open_file(Dtfil%fnameabo_sig, msg, unit=unt_sig, status='unknown', form='formatted') /= 0) then
1977      ABI_ERROR(msg)
1978    end if
1979    write(unt_sig,"(a)")"# Sigma_xc and spectral function A along the real frequency axis in eV units"
1980    write(unt_sig,"(a)")"# Format:"
1981    write(unt_sig,"(a)")"#   kpoint"
1982    write(unt_sig,"(a)")"#   min_band(k) max_band(k)"
1983    write(unt_sig,"(a)")"#   For each frequency w:"
1984    write(unt_sig,"(2a)")"#       w, {Re(Sigma_b(w)), Im(Sigma_b(w), A_b(w) for b in [min_band, max_band]}",ch10
1985 
1986    if (open_file(Dtfil%fnameabo_sgr, msg, unit=unt_sgr, status='unknown', form='formatted') /= 0) then
1987      ABI_ERROR(msg)
1988    end if
1989    write(unt_sgr,"(a)")"# Derivatives of Sigma_c(omega) wrt omega in eV units"
1990 
1991    !  Sigma_c(eik) MRM
1992    if (open_file(trim(Dtfil%fnameabo_sgr)//'_SIGC', msg, unit=unt_sigc, status='unknown', form='formatted') /= 0) then
1993      ABI_ERROR(msg)
1994    end if
1995 
1996    if (mod10 == SIG_GW_AC) then
1997      ! Sigma along the imaginary axis.
1998      if (open_file(Dtfil%fnameabo_sgm, msg, unit=unt_sgm, status='unknown', form='formatted') /= 0) then
1999        ABI_ERROR(msg)
2000      end if
2001      write(unt_sgm,"(a)")"# Sigma_xc along the imaginary frequency axis in eV units"
2002    end if
2003  end if
2004 
2005   !=======================================================================
2006   !==== Calculate self-energy and output the results for each k-point ====
2007   !=======================================================================
2008   ! Here it would be possible to calculate the QP correction for the same k-point using a PPmodel
2009   ! in the first iteration just to improve the initial guess and CD or AC in the second step. Useful
2010   ! if the KS level is far from the expected QP results. Previously it was possible to do such
2011   ! calculation by simply listing the same k-point twice in kptgw. Now this trick is not allowed anymore.
2012   ! Everything, indeed, should be done in a clean and transparent way inside csigme.
2013 
2014  call wfd%print(mode_paral='PERS')
2015 
2016  call wrtout(std_out,sigma_type_from_key(mod10))
2017 
2018  if (gwcalctyp<10) then
2019    msg = " Perturbative Calculation"
2020    if (gwcalctyp==1) msg = " Newton Raphson method "
2021  else if (gwcalctyp<20) then
2022    msg = " Self-Consistent on Energies only"
2023  else if (gwcalctyp==21) then
2024    msg = " GW 1-RDM Corrections"
2025  else
2026    msg = " Self-Consistent on Energies and Wavefunctions"
2027  end if
2028  if (Dtset%ucrpa==0) call wrtout(std_out,msg)
2029 
2030  !=================================================
2031  !==== Calculate the matrix elements of Sigma =====
2032  !=================================================
2033 
2034  nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i
2035 
2036  ! min and max band indeces for GW corrections.
2037  ib1=Sigp%minbdgw; ib2=Sigp%maxbdgw
2038 
2039  !MG TODO: I don't like the fact that ib1 and ib2 are redefined here because this
2040  ! prevents me from refactoring the code. In particular I want to store the self-energy
2041  ! results inside the sigma_results datatypes hence one needs to know all the dimensions
2042  ! at the beginning of the execution (e.g. in setup_sigma) so that one can easily allocate the arrays in the type.
2043  if(Dtset%ucrpa>=1) then
2044    !Read the band
2045    if (dtset%plowan_compute<10)then
2046      if (open_file("forlb.ovlp",msg,newunit=temp_unt,form="formatted", status="unknown") /= 0) then
2047        ABI_ERROR(msg)
2048      end if
2049      rewind(temp_unt)
2050      read(temp_unt,*)
2051    read(temp_unt,*)
2052      read(temp_unt,*) msg, ib1, ib2
2053      close(temp_unt)
2054    else
2055      if (open_file("data.plowann",msg,newunit=temp_unt,form="formatted", status="unknown") /= 0) then
2056        ABI_ERROR(msg)
2057      end if
2058      rewind(temp_unt)
2059      read(temp_unt,*)
2060      read(temp_unt,*)
2061      read(temp_unt,'(a7,2i4)') msg, ib1, ib2
2062      close(temp_unt)
2063  endif
2064 endif
2065 
2066  ! Do not store it for rdm_update because it might be too large!
2067  if(.not.rdm_update) then
2068    ABI_MALLOC(sigcme,(nomega_sigc,ib1:ib2,ib1:ib2,Sigp%nkptgw,Sigp%nsppol*Sigp%nsig_ab))
2069    sigcme=czero
2070  endif
2071 
2072  ! if (.False. .and. psps%usepaw == 0 .and. wfd%nspinor == 1 .and. any(dtset%so_psp /= 0)) then
2073  !   call wrtout(std_out, "Computing SOC contribution with first-order perturbation theory")
2074  !   ABI_MALLOC(bks_mask, (wfd%mband, wfd%nkibz, wfd%nsppol))
2075  !   bks_mask = .False.
2076  !   do spin=1,wfd%nsppol
2077  !     do ikcalc=1,Sigp%nkptgw
2078  !       ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW
2079  !       ii=Sigp%minbnd(ikcalc, spin); jj=Sigp%maxbnd(ikcalc, spin)
2080  !       bks_mask(ii:jj, ik_ibz, spin) = .True.
2081  !     end do
2082  !   end do
2083  !
2084  !   call wfd_get_socpert(wfd, cryst, psps, pawtab, bks_mask, osoc_bks)
2085  !   ABI_FREE(bks_mask)
2086  !   ABI_FREE(osoc_bks)
2087  ! end if
2088 
2089  !==========================================================
2090  !==== Exchange part using the dense gwx_ngfft FFT mesh ====
2091  !==========================================================
2092  !TODO : load distribution is not optimal if band parallelism is used.
2093  !but this problem was also affecting the old implementation.
2094  call timab(421,1,tsec) ! calc_sigx_me
2095 
2096  ! This routine shows how the wavefunctions are distributed.
2097  !call wfd_show_bkstab(Wfd,unit=std_out)
2098 
2099  if(Dtset%ucrpa>=1) then
2100    ! Calculation of the Rho_n for calc_ucrpa
2101    call wrtout(std_out, sjoin("begin of Ucrpa calc for a nb of kpoints of: ",itoa(Sigp%nkptgw)))
2102    ! Wannier basis: rhot1_q_m will need an atom index in the near future.
2103    ic=0
2104    do  itypat=1,cryst%ntypat
2105      if(pawtab(itypat)%lpawu.ne.-1) then
2106        ndim=2*dtset%lpawu(itypat)+1
2107        ic=ic+1
2108        itypatcor=itypat
2109        lcor=dtset%lpawu(itypat)
2110      end if
2111    end do
2112    if(ic>1) then
2113      ABI_ERROR("number of correlated species is larger than one")
2114    end if
2115    ABI_MALLOC(rhot1_q_m,(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz))
2116    ABI_MALLOC(M1_q_m,(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz))
2117 
2118    M1_q_m=czero
2119    rhot1_q_m=czero
2120 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2121 !Initialization of wan objects and getting psichies
2122 !Allocation of rhot1
2123 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2124    if (dtset%plowan_compute>=10)then
2125      write(msg,'(a)')" cRPA calculations using wannier weights from data.plowann"
2126      call wrtout([std_out, ab_out], msg)
2127      call init_plowannier(dtset%plowan_bandf,dtset%plowan_bandi,dtset%plowan_compute,dtset%plowan_iatom,&
2128        &dtset%plowan_it,dtset%plowan_lcalc,dtset%plowan_natom,dtset%plowan_nbl,dtset%plowan_nt,&
2129        &dtset%plowan_projcalc,dtset%acell_orig,dtset%kptns,dtset%nimage,dtset%nkpt,dtset%nspinor,&
2130        &dtset%nsppol,dtset%wtk,dtset%dmft_t2g,wanibz_in)
2131      call get_plowannier(wanibz_in,wanibz,dtset)
2132      call fullbz_plowannier(dtset,kmesh,cryst,pawang,wanibz,wanbz)
2133      ABI_MALLOC(rhot1,(sigp%npwx,Qmesh%nibz))
2134      do pwx=1,sigp%npwx
2135        do ibz=1,Qmesh%nibz
2136          call init_operwan_realspace(wanbz,rhot1(pwx,ibz))
2137          call zero_operwan_realspace(wanbz,rhot1(pwx,ibz))
2138        enddo
2139      enddo
2140    endif
2141 
2142 !   do ikcalc=1,Sigp%nkptgw
2143 
2144 !   if(cryst%nsym==1) nkcalc=Kmesh%nbz
2145 !   if(cryst%nsym>1)  nkcalc=Kmesh%nibz
2146    nkcalc=Kmesh%nbz
2147    !if(1==1)then!DEBUG
2148    !open(67,file="test.rhot1",status="REPLACE")
2149    do ikcalc=1,nkcalc ! for the oscillator strengh, spins are identical without SOC
2150 !    if(Sigp%nkptgw/=Kmesh%nbz) then
2151 !      write(msg,'(6a)')ch10,&
2152 !&      ' nkptgw and nbz differs: this is not allowed to compute U in cRPA '
2153 !      ABI_ERROR(msg)
2154 !    endif
2155      write(std_out,*) " ikcalc",ikcalc,size(Sigp%kptgw2bz),nkcalc
2156 
2157      if(mod(ikcalc-1,nprocs)==Wfd%my_rank) then
2158        if(cryst%nsym==1) ik_ibz=Kmesh%tab(ikcalc) ! Index of the irred k-point for GW
2159        if(cryst%nsym>1) ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2160 
2161        call prep_calc_ucrpa(ik_ibz,ikcalc,itypatcor,ib1,ib2,Cryst,qp_ebands,Sigp,Gsph_x,Vcp,&
2162 &       Kmesh,Qmesh,lcor,M1_q_m,Pawtab,Pawang,Paw_pwff,&
2163 &       Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,gwx_ngfft,&
2164 &       ngfftf,Dtset%prtvol,Dtset%pawcross,Dtset%plowan_compute,rhot1_q_m,wanbz,rhot1)
2165      end if
2166    end do
2167    if (dtset%plowan_compute<10)then
2168      call xmpi_sum(rhot1_q_m,Wfd%comm,ierr)
2169      call xmpi_sum(M1_q_m,Wfd%comm,ierr)
2170      M1_q_m=M1_q_m/Kmesh%nbz/Wfd%nsppol
2171      rhot1_q_m=rhot1_q_m/Kmesh%nbz/Wfd%nsppol
2172      ! do pwx=1,sigp%npwx
2173      !   do ibz=1,Qmesh%nibz
2174      !     do ispinor1=1,dtset%nspinor
2175      !       do ispinor2=1,dtset%nspinor
2176      !         do iatom1=1,cryst%nattyp(itypatcor)
2177      !           do im1=1,2*lcor+1
2178      !             do im2=1,2*lcor+1
2179      !               write(67,*)ibz,im1,im2,rhot1_q_m(iatom1,ispinor1,ispinor2,im1,im2,pwx,ibz)
2180      !             enddo
2181      !           enddo
2182      !         enddo
2183      !       enddo
2184      !     enddo
2185      !   enddo
2186      ! enddo
2187    else
2188      call reduce_operwan_realspace(wanbz,rhot1,sigp%npwx,Qmesh%nibz,Wfd%comm,Kmesh%nbz,Wfd%nsppol)
2189      !call cwtime(cpu,wall,gflops,"start") !reduction of rhot1
2190      ! dim=0
2191      ! do pwx=1,sigp%npwx
2192      ! do ibz=1,Qmesh%nibz
2193      !   do spin=1,wanbz%nsppol
2194      !   do ispinor1=1,wanbz%nspinor
2195      !   do ispinor2=1,wanbz%nspinor
2196      !     do iatom1=1,wanbz%natom_wan
2197      !     do iatom2=1,wanbz%natom_wan
2198      !       do pos1=1,size(wanbz%nposition(iatom1)%pos,1)
2199      !       do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
2200      !         do il1=1,wanbz%nbl_atom_wan(iatom1)
2201      !         do il2=1,wanbz%nbl_atom_wan(iatom2)
2202      !           do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
2203      !           do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
2204      ! dim=dim+1
2205      !           enddo!im2
2206      !           enddo!im1
2207      !         enddo!il2
2208      !         enddo!il1
2209      !       enddo!pos2
2210      !       enddo!pos1
2211      !     enddo!iatom2
2212      !     enddo!iatom1
2213      !   enddo!ispinor2
2214      !   enddo!ispinor1
2215      !   enddo!spin
2216      ! enddo!ibz
2217      ! enddo!pwx
2218      ! ABI_MALLOC(buffer,(dim))
2219      ! nnn=0
2220      ! do pwx=1,sigp%npwx
2221      ! do ibz=1,Qmesh%nibz
2222      !     do spin=1,wanbz%nsppol
2223      !     do ispinor1=1,wanbz%nspinor
2224      !     do ispinor2=1,wanbz%nspinor
2225      !       do iatom1=1,wanbz%natom_wan
2226      !       do iatom2=1,wanbz%natom_wan
2227      !         do pos1=1,size(wanbz%nposition(iatom1)%pos,1)
2228      !         do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
2229      !           do il1=1,wanbz%nbl_atom_wan(iatom1)
2230      !           do il2=1,wanbz%nbl_atom_wan(iatom2)
2231      !             do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
2232      !             do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
2233      ! nnn=nnn+1
2234      ! buffer(nnn)=rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)
2235      !             enddo!im2
2236      !             enddo!im1
2237      !           enddo!il2
2238      !           enddo!il1
2239      !         enddo!pos2
2240      !         enddo!pos1
2241      !       enddo!iatom2
2242      !       enddo!iatom1
2243      !     enddo!ispinor2
2244      !     enddo!ispinor1
2245      !   enddo!spin
2246      ! enddo!ibz
2247      ! enddo!pwx
2248      ! call xmpi_barrier(Wfd%comm)
2249      ! call xmpi_sum(buffer,Wfd%comm,ierr)
2250      ! call xmpi_barrier(Wfd%comm)
2251      ! buffer=buffer/Kmesh%nbz/Wfd%nsppol
2252      ! nnn=0
2253      ! do pwx=1,sigp%npwx
2254      ! do ibz=1,Qmesh%nibz
2255      !   do spin=1,wanbz%nsppol
2256      !   do ispinor1=1,wanbz%nspinor
2257      !   do ispinor2=1,wanbz%nspinor
2258      !     do iatom1=1,wanbz%natom_wan
2259      !     do iatom2=1,wanbz%natom_wan
2260      !       do pos1=1,size(wanbz%nposition(iatom1)%pos,1)
2261      !       do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
2262      !         do il1=1,wanbz%nbl_atom_wan(iatom1)
2263      !         do il2=1,wanbz%nbl_atom_wan(iatom2)
2264      !           do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
2265      !           do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
2266      !  nnn=nnn+1
2267      !  rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)=buffer(nnn)
2268      !  !write(67,*)ibz,im1,im2,rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)
2269      !           enddo!im2
2270      !           enddo!im1
2271      !         enddo!il2
2272      !         enddo!il1
2273      !       enddo!pos2
2274      !       enddo!pos1
2275      !     enddo!iatom2
2276      !     enddo!iatom1
2277      !   enddo!ispinor2
2278      !   enddo!ispinor1
2279      !   enddo!spin
2280      ! enddo!ibz
2281      ! enddo!pwx
2282      ! ABI_FREE(buffer)
2283    endif
2284 
2285    !close(67)
2286    ! call xmpi_barrier(Wfd%comm)
2287    ! call cwtime(cpu,wall,gflops,"stop")!reduction of rhot1
2288    ! write(6,*)cpu,wall,gflops
2289 !   if(Cryst%nsym==1) then
2290 !   M1_q_m=M1_q_m/Kmesh%nbz/Wfd%nsppol
2291 !   rhot1_q_m=rhot1_q_m/Kmesh%nbz/Wfd%nsppol
2292 !   endif
2293 !  Calculation of U in cRPA: need to treat with a different cutoff the
2294 !  bare coulomb and the screened coulomb interaction.
2295    ! else !DEBUG
2296    !   write(6,*)"DEBUGGING calc_ucrpa"
2297    !   open(67,file="test.rhot1",status="OLD")
2298    !   rewind(67)
2299    !    do pwx=1,sigp%npwx
2300    !     do ibz=1,Qmesh%nibz
2301    !       do spin=1,wanbz%nsppol
2302    !       do ispinor1=1,wanbz%nspinor
2303    !         do ispinor2=1,wanbz%nspinor
2304    !           do iatom1=1,wanbz%natom_wan
2305    !             do iatom2=1,wanbz%natom_wan
2306    !               do pos1=1,size(wanbz%nposition(iatom1)%pos,1)
2307    !                 do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
2308    !                   do il1=1,wanbz%nbl_atom_wan(iatom1)
2309    !                     do il2=1,wanbz%nbl_atom_wan(iatom2)
2310    !                       do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
2311    !                         do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
2312    !    read(67,*)dummy,dummy,dummy,xx
2313    !    rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)=xx
2314    !                         enddo!im2
2315    !                       enddo!im1
2316    !                     enddo!il2
2317    !                   enddo!il1
2318    !                 enddo!pos2
2319    !               enddo!pos1
2320    !             enddo!iatom2
2321    !           enddo!iatom1
2322    !         enddo!ispinor2
2323    !       enddo!ispinor1
2324    !       enddo!spin
2325    !     enddo!ibz
2326    !   enddo!pwx
2327    !   close(67)
2328    ! endif!DEBUG
2329    call calc_ucrpa(itypatcor,cryst,Kmesh,lcor,M1_q_m,Qmesh,Er%npwe,sigp%npwx,&
2330 &   Cryst%nsym,Sigp%nomegasr,Sigp%minomega_r,Sigp%maxomega_r,ib1,ib2,&
2331 &'Gsum',Cryst%ucvol,Wfd,Er%fname,dtset%plowan_compute,rhot1,wanbz)
2332 
2333 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2334 !Deallocation of wan and rhot1
2335 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2336    if (dtset%plowan_compute >=10) then
2337      do pwx=1,sigp%npwx
2338        do ibz=1,Qmesh%nibz
2339          call destroy_operwan_realspace(wanbz,rhot1(pwx,ibz))
2340        enddo
2341      enddo
2342      ABI_FREE(rhot1)
2343      call destroy_plowannier(wanbz)
2344    endif
2345    ABI_FREE(rhot1_q_m)
2346    ABI_FREE(M1_q_m)
2347 
2348  else
2349    ! MRM: sigmak_todo is an array indicating whether the k-point must be computed for GW correction. This array is set to DO IT (to 1) for any GW calculation
2350    ! and it will only change to NOT DO IT (to 0) in case that a density matrix update calc. is used and we are reading checkpoint files
2351    ABI_MALLOC(sigmak_todo,(Wfd%nkibz))
2352    sigmak_todo(:)=1
2353 
2354    ! ======================================================
2355    ! ==== Initialize arrays for density matrix updates ====
2356    ! ======================================================
2357    if (rdm_update) then
2358      if (Sigp%nsppol/=1) then
2359        ABI_ERROR("1-RDM GW correction only implemented for restricted closed-shell calculations!")
2360        ! Note: all subroutines of 70_gw/m_gwrdm.F90 are implemented assuming Sigp%nsppol==1
2361      end if
2362      ABI_MALLOC(nateigv,(Wfd%mband,Wfd%mband,Wfd%nkibz,Sigp%nsppol))
2363      ABI_MALLOC(nat_occs,(Wfd%mband,Wfd%nkibz))
2364      ABI_MALLOC(xrdm_k_full,(b1gw:b2gw,b1gw:b2gw,Wfd%nkibz))
2365      write(msg,'(a34,2i9)')' Bands used for the GW 1RDM arrays',b1gw,b2gw
2366      call wrtout([std_out, ab_out], msg)
2367      xrdm_k_full=czero; nat_occs=zero; nateigv=czero;
2368      do ik_ibz=1,Wfd%nkibz
2369        do ib=b1gw,b2gw
2370          xrdm_k_full(ib,ib,ik_ibz)=qp_ebands%occ(ib,ik_ibz,1)
2371        enddo
2372        do ib=1,Wfd%mband
2373          nat_occs(ib,ik_ibz)=qp_ebands%occ(ib,ik_ibz,1)   ! Copy initial occ numbers (in principle 2 or 0 from KS-DFT)
2374          nateigv(ib,ib,ik_ibz,1)=cone                  ! Set to identity matrix
2375        end do
2376      enddo
2377      if (readchkprdm) then
2378        gw1rdm_fname=trim(dtfil%fnameabi_chkp_rdm)
2379        call get_chkprdm(Wfd,Kmesh,Sigp,qp_ebands,nat_occs,nateigv,sigmak_todo,my_rank,gw1rdm_fname)
2380      end if
2381      call xmpi_barrier(Wfd%comm)
2382      ! Prepare arrays for the imaginary freq. integration (quadrature) of Sigma_c(iw)
2383      ABI_MALLOC(weights,(Sigp%nomegasi))
2384      call quadrature_sigma_cw(Sigp,Sr,weights)
2385    end if
2386 
2387    ! ===================================
2388    ! ==== Static part (Sigma_x-Vxc) ====
2389    ! ===================================
2390    do ikcalc=1,Sigp%nkptgw
2391      ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2392      if (sigmak_todo(ik_ibz)==1) then ! Do not compute MELS if the k-point was read from the checkpoint file
2393                                       ! this IF only affects GW density matrix update!
2394        ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2395        ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2396        call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,qp_ebands,Sigp,Sr,Gsph_x,Vcp,Kmesh,Qmesh,Ltg_k(ikcalc),&
2397                          Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,&
2398                          gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty)
2399        if (rdm_update) then ! Only recompute exchange update to the 1-RDM if the point read was broken or not precomputed
2400          ABI_MALLOC(pot_k,(ib1:ib2,ib1:ib2))
2401          ABI_MALLOC(rdm_k,(ib1:ib2,ib1:ib2))
2402          rdm_k=czero
2403 !        Compute Sigma_x - Vxc or DELTA Sigma_x - Vxc. (DELTA Sigma_x = Sigma_x - hyb_parameter Vx^exact for hyb Functionals)
2404          pot_k(ib1:ib2,ib1:ib2)=Sr%x_mat(ib1:ib2,ib1:ib2,ik_ibz,1)-KS_me%vxcval(ib1:ib2,ib1:ib2,ik_ibz,1) ! Only restricted calcs
2405          call calc_rdmx(ib1,ib2,ik_ibz,pot_k,rdm_k,qp_ebands)                                                ! Only restricted closed-shell calcs
2406 !        Update the full 1RDM with the exchange corrected one for this k-point
2407          xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz)=xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz)+rdm_k(ib1:ib2,ib1:ib2)
2408 !        Compute NAT ORBS for exchange corrected 1-RDM
2409          do ib=ib1,ib2
2410            rdm_k(ib,ib)=rdm_k(ib,ib)+qp_ebands%occ(ib,ik_ibz,1)           ! Only restricted closed-shell calcs
2411          enddo
2412          call natoccs(ib1,ib2,rdm_k,nateigv,nat_occs,qp_ebands,ik_ibz,0)  ! Only restricted closed-shell calcs
2413          ABI_FREE(pot_k)
2414          ABI_FREE(rdm_k)
2415        end if
2416      else
2417        write(msg,'(a1)')  ' '
2418        call wrtout(std_out,msg)
2419        write(msg,'(a64,i5)')  ' Skipping the calc. of Sigma_x - ( Vxc + hyb*Fock) for k-point: ',ik_ibz
2420        call wrtout(std_out,msg)
2421        write(msg,'(a1)')  ' '
2422        call wrtout(std_out,msg)
2423      end if
2424    end do
2425    ! for the time being, do not remove this barrier!
2426    call xmpi_barrier(Wfd%comm)
2427    call timab(421,2,tsec) ! calc_sigx_me
2428 
2429    ! ==========================================================
2430    ! ==== Correlation part using the coarse gwc_ngfft mesh ====
2431    ! ==========================================================
2432    if (mod10/=SIG_HF) then
2433      do ikcalc=1,Sigp%nkptgw
2434        ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2435        ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2436        ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2437 
2438        ABI_MALLOC(sigcme_k,(nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab))
2439        sigcme_k=czero
2440        if (any(mod10 == [SIG_SEX, SIG_COHSEX])) then
2441          ! Calculate static COHSEX or SEX using the coarse gwc_ngfft mesh.
2442          call cohsex_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Cryst,qp_ebands,Sigp,Sr,Er,Gsph_c,Vcp,Kmesh,Qmesh,&
2443                         Ltg_k(ikcalc),Pawtab,Pawang,Paw_pwff,Psps,Wfd,QP_sym,&
2444                         gwc_ngfft,Dtset%iomode,Dtset%prtvol,sigcme_k)
2445        else
2446           ! Compute correlated part using the coarse gwc_ngfft mesh.
2447           if (x1rdm/=1 .and. sigmak_todo(ik_ibz)==1) then ! Do not compute correlation MELS if the k-point was read from the checkpoint file
2448                                                           ! this IF only affects GW density matrix update
2449             call calc_sigc_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Dtset,Cryst,qp_ebands, &
2450                               Sigp,Sr,Er,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,&
2451                               Ltg_k(ikcalc),PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,&
2452                               gwc_ngfft,ngfftf,nfftf,ks_rhor,use_aerhor,ks_aepaw_rhor,sigcme_k)
2453           else
2454             write(msg,'(a1)')  ' '
2455             call wrtout(std_out, msg)
2456             write(msg,'(a44,i5)')  ' Skipping the calc. of Sigma_c for k-point: ',ik_ibz
2457             call wrtout(std_out, msg)
2458             write(msg,'(a1)')  ' '
2459             call wrtout(std_out, msg)
2460           end if
2461        end if
2462        ! MRM: compute density matrix numerical correction (freq. integration G0 Sigma_c(iw) G0).
2463        if (rdm_update) then
2464          if (sigmak_todo(ik_ibz)==1) then ! Only recompute correlation update to the 1-RDM if the point read was broken or not precomputed
2465            if (x1rdm/=1) then
2466              ABI_MALLOC(rdm_k,(ib1:ib2,ib1:ib2))
2467              rdm_k=czero
2468              call calc_rdmc(ib1,ib2,ik_ibz,Sr,weights,sigcme_k,qp_ebands,rdm_k)    ! Only restricted closed-shell calcs
2469 !            Update the full 1RDM with the GW corrected one for this k-point
2470              rdm_k(ib1:ib2,ib1:ib2)=xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz)+rdm_k(ib1:ib2,ib1:ib2)
2471 !            Compute nat orbs and occ numbers at k-point ik_ibz
2472              call natoccs(ib1,ib2,rdm_k,nateigv,nat_occs,qp_ebands,ik_ibz,1) ! Only restricted closed-shell calcs
2473              ABI_FREE(rdm_k)
2474            endif
2475 !          Print the checkpoint file if required
2476            if(prtchkprdm) then
2477              gw1rdm_fname=trim(dtfil%fnameabo_chkp_rdm)
2478              call print_chkprdm(Wfd,nat_occs,nateigv,ik_ibz,my_rank,gw1rdm_fname)
2479            end if
2480          end if
2481        else
2482          sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:)=sigcme_k
2483        end if
2484        ABI_FREE(sigcme_k)
2485      end do
2486      call xmpi_barrier(Wfd%comm)
2487      if (rdm_update) then
2488        ABI_FREE(xrdm_k_full)
2489        ABI_FREE(weights)
2490      end if
2491    end if
2492 
2493    call xmpi_barrier(Wfd%comm)
2494    ABI_FREE(sigmak_todo)
2495 
2496    !
2497    ! 1) Print WFK and DEN files. 2) Build band corrections and compute new energies.
2498    !
2499    if (rdm_update) then
2500      ABI_MALLOC(gw_rhor,(nfftf,Dtset%nspden))
2501      gw_rhor=zero
2502      !
2503      ! WARNING!
2504      ! MRM: only the master has bands on Wfd_nato_master so it prints everything and computes gw_rhor
2505      !
2506      ! All procs. update the qp_ebands and the Hdr_sigma
2507      call update_hdr_bst(Wfd_nato_master,nat_occs,b1gw,b2gw,qp_ebands,Hdr_sigma,Dtset%ngfft(1:3))
2508      call print_tot_occ(Sr,Kmesh, qp_ebands) ! Compute unit cell (averaged) occ = \sum _k weight_k occ_k
2509      if (my_rank==master) then
2510        call Wfd_nato_master%rotate(Cryst,nateigv,bdm_mask)                          ! Let it use bdm_mask and build NOs
2511        call Wfd_nato_master%mkrho(Cryst,Psps,Kmesh,qp_ebands,ngfftf,nfftf,gw_rhor)     ! Construct the density
2512        if (dtset%prtwf == 1) then
2513          gw1rdm_fname=dtfil%fnameabo_wfk
2514          ! Print WFK file, here qp_ebands contains nat. orb. occs.
2515          call Wfd_nato_master%write_wfk(Hdr_sigma,qp_ebands,gw1rdm_fname,wfknocheck)
2516        end if
2517        if (dtset%prtden == 1) then
2518          gw1rdm_fname=dtfil%fnameabo_den
2519          call fftdatar_write("density",gw1rdm_fname,dtset%iomode,Hdr_sigma,&        ! Print DEN file
2520          Cryst,ngfftf,cplex1,nfftf,dtset%nspden,gw_rhor,mpi_enreg_seq,ebands=qp_ebands)
2521        end if
2522      end if
2523      ierr=0
2524      call xmpi_bcast(gw_rhor(:,:),master,Wfd%comm,ierr)
2525      if(ierr/=0) then
2526        ABI_ERROR("Error distributing the GW density")
2527      endif
2528      ! We no longer need Wfd_nato_master.
2529      call Wfd_nato_master%free()
2530      call xmpi_barrier(Wfd%comm)
2531      ABI_FREE(bdm_mask) ! The master already used bdm_mask
2532      ABI_FREE(nat_occs)     ! Occs were already placed in qp_ebands
2533      call em1results_free(Er) ! We no longer need Er for GW@KS-DFT 1RDM but we may need space on the RAM memory
2534      if (gw1rdm==2 .and. Sigp%nkptgw==Wfd%nkibz) then  ! Compute energies only if all k-points are available
2535        ABI_MALLOC(old_ks_purex,(b1gw:b2gw,Sigp%nkptgw))   ! We need the hole 1-RDM to build Fock[GW.1RDM]!
2536        ABI_MALLOC(new_hartr,(b1gw:b2gw,Sigp%nkptgw))
2537        ABI_MALLOC(gw_rhog,(2,nfftf))
2538        ABI_MALLOC(gw_vhartr,(nfftf))
2539        gw_rhog=zero
2540        gw_vhartr(:)=zero
2541        old_ks_purex(:,:)=czero
2542        new_hartr(:,:)=czero
2543        !
2544        ! A) Compute Evext = int rho(r) vext(r) dr -> simply dot product on the FFT grid
2545        !
2546        den_int=sum(gw_rhor(:,1))*ucvol_local/nfftf               ! Only restricted closed-shell calcs
2547        evext_energy=sum(gw_rhor(:,1)*vpsp(:))*ucvol_local/nfftf  ! Only restricted closed-shell calcs
2548        !
2549        ! B) Coulomb <KS_i|Vh[NO]|KS_j>
2550        !
2551        call fourdp(1,gw_rhog,gw_rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,tim_fourdp5)                  ! FFT to build gw_rhog
2552        ecutf=dtset%ecut
2553        if (psps%usepaw==1) then
2554          ecutf=dtset%pawecutdg
2555          call wrtout(std_out,ch10//' FFT (fine) grid used in GW update:')
2556        end if
2557        call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf)
2558        call hartre(1,gsqcut,dtset%icutcoul,Psps%usepaw,MPI_enreg_seq,nfftf,ngfftf,dtset%nkpt,dtset%rcut,gw_rhog,&
2559        Cryst%rprimd,dtset%vcutgeo,gw_vhartr)                                                            ! Build Vhartree -> gw_vhartr
2560        ABI_MALLOC(tmp_kstab,(2,Wfd%nkibz,Wfd%nsppol))
2561        tmp_kstab=0
2562        do spin=1,Sigp%nsppol
2563          do ikcalc=1,Sigp%nkptgw
2564            ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
2565            tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin)
2566            tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin)
2567          end do
2568        end do
2569        call calc_vhxc_me(Wfd,KS_mflags,GW1RDM_me,Cryst,Dtset,nfftf,ngfftf,&           ! Build matrix elements from gw_vhartr -> GW1RDM_me
2570        & ks_vtrial,gw_vhartr,ks_vxc,Psps,Pawtab,KS_paw_an,Pawang,Pawfgrtab,KS_paw_ij,dijexc_core,&
2571        & gw_rhor,usexcnhat,ks_nhat,ks_nhatgr,nhatgrdim,tmp_kstab,taur=ks_taur)
2572        call xmpi_barrier(Wfd%comm)
2573        ABI_FREE(tmp_kstab)
2574        ABI_FREE(gw_rhor)
2575        ABI_FREE(gw_rhog)
2576        ABI_FREE(gw_vhartr)
2577        ! Save new <i|Hartree[NO]|i> in KS basis for Delta eik
2578        do ikcalc=1,Sigp%nkptgw
2579          ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW
2580          ib1=MINVAL(Sigp%minbnd(ikcalc,:))       ! min and max band indices for GW corrections (for this k-point)
2581          ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2582          do ib=b1gw,b2gw
2583            new_hartr(ib,ikcalc)=GW1RDM_me%vhartree(ib,ib,ik_ibz,1)
2584          end do
2585        end do
2586        !
2587        ! C) Exchange <KS_i|hyb*K^RANGE-SEP?[KS]|KS_j> and save it in old_ks_purex
2588        !
2589        ! Build Vcp_ks and Vcp_full
2590        call setup_vcp(Vcp_ks,Vcp_full,Dtset,Gsph_x,Gsph_c,Cryst,Qmesh,Kmesh,coef_hyb,comm)
2591        call xmpi_barrier(Wfd%comm)
2592        if(coef_hyb>tol8) then
2593          do ikcalc=1,Sigp%nkptgw
2594            ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW
2595            ib1=MINVAL(Sigp%minbnd(ikcalc,:))       ! min and max band indices for GW corrections (for this k-point)
2596            ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2597            call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,ks_ebands,Sigp,Sr,Gsph_x,Vcp_ks,Kmesh,Qmesh,Ltg_k(ikcalc),& ! Notice we need KS occs
2598            & Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,&                                   ! and band energy diffs
2599            & gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty)
2600            ! Build <KS_i|RS?_Hyb?_Sigma_x[KS]|KS_j> matrix (Use Wfd)
2601            call xmpi_barrier(Wfd%comm)
2602            do ib=b1gw,b2gw
2603              old_ks_purex(ib,ikcalc)=Sr%x_mat(ib,ib,ik_ibz,1)              ! Save old alpha*<i|K[KS]|i> from the GS calc. for Delta eik
2604            enddo
2605          end do
2606        endif
2607        call xmpi_barrier(Wfd%comm)
2608        call Vcp_ks%free()
2609        !
2610        ! Use the Wfd with 'new name' (Wfd_nato_all) because it will contain the GW 1RDM nat. orbs. (bands)
2611        !
2612        ABI_COMMENT("From now on, the Wfd bands will contain the nat. orbs. ones")
2613        Wfd_nato_all => Wfd
2614        call Wfd_nato_all%rotate(Cryst,nateigv)                               ! Let rotate build the NOs in Wfd_nato_all (KS->NO)
2615        call xmpi_barrier(Wfd%comm)
2616        !
2617        ! D) Non-local <NO_i|Vnl|NO_i> terms [saved on nl_bks(band,k,spin)]
2618        !
2619        call Wfd_nato_all%get_nl_me(Cryst,Psps,Pawtab,bdm2_mask,nl_bks)
2620        evextnl_energy=zero
2621        do spin=1,Sigp%nsppol
2622          do ikcalc=1,Sigp%nkptgw
2623            do ib=Sr%b1gw,Sr%b2gw
2624              evextnl_energy=evextnl_energy+Kmesh%wt(ikcalc)*qp_ebands%occ(ib,ikcalc,spin)*nl_bks(ib,ikcalc,spin)
2625            end do
2626          end do
2627        end do
2628        call xmpi_barrier(Wfd%comm)
2629        ABI_FREE(nl_bks)
2630        ABI_FREE(bdm2_mask)
2631        !
2632        ! E) Exchange <NO_i|K[NO]|NO_j>
2633        !
2634        ! Get the matrix elements and save them on Sr%x_mat
2635        tol_empty=0.0001 ! Allow lower occ numbers
2636        do ikcalc=1,Sigp%nkptgw
2637          ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW
2638          ib1=MINVAL(Sigp%minbnd(ikcalc,:))       ! min and max band indices for GW corrections (for this k-point)
2639          ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2640          ! Build <NO_i|Sigma_x[NO]|NO_j> matrix
2641          call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,qp_ebands,Sigp,Sr,Gsph_x,Vcp_full,Kmesh,Qmesh,Ltg_k(ikcalc),&
2642                             Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd_nato_all,Wfdf,QP_sym,&
2643                             gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty)
2644        end do
2645        tol_empty=0.01 ! Recover standard value for tolerance on the occ numbers
2646        call xmpi_barrier(Wfd%comm)
2647        call Vcp_full%free()
2648        ex_energy=sigma_get_exene(Sr,Kmesh,qp_ebands)         ! Save the new total exchange energy Ex = Ex[GW.1RDM]
2649        exc_mbb_energy=sigma_get_excene(Sr,Kmesh,qp_ebands)   ! Save the new total exchange-correlation MBB energy Exc = Exc^MBB[GW.1RDM]
2650        !
2651        ! Transform <NO_i|K[NO]|NO_j> -> <KS_i|K[NO]|KS_j>,
2652        !           <KS_i|J[NO]|KS_j> -> <NO_i|J[NO]|NO_j>,
2653        !       and     <KS_i|T|KS_j> -> <NO_i|T|NO_j>
2654        !
2655        call change_matrix(Sigp,Sr,GW1RDM_me,Kmesh,nateigv)
2656        call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data
2657        !
2658        ! Print Delta eik and band (state) energies
2659        !
2660        call print_band_energies(b1gw,b2gw,Sr,Sigp,KS_me,Kmesh,ks_ebands,new_hartr,old_ks_purex)
2661        call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data
2662        !
2663        ! Print the updated total energy and all energy components
2664        !
2665        eh_energy=mels_get_haene(Sr,GW1RDM_me,Kmesh,qp_ebands)
2666        ekin_energy=mels_get_kiene(Sr,GW1RDM_me,Kmesh,qp_ebands)
2667        etot_sd=ekin_energy+evext_energy+evextnl_energy+QP_energies%e_corepsp+QP_energies%e_ewald+eh_energy+ex_energy       ! SD 2-RDM
2668        etot_mbb=ekin_energy+evext_energy+evextnl_energy+QP_energies%e_corepsp+QP_energies%e_ewald+eh_energy+exc_mbb_energy ! MBB 2-RDM
2669        call print_total_energy(ekin_energy,evext_energy,evextnl_energy,QP_energies%e_corepsp,eh_energy,ex_energy,&
2670        &exc_mbb_energy,QP_energies%e_ewald,etot_sd,etot_mbb,den_int)
2671        !
2672        ! Clean GW1RDM_me and allocated arrays
2673        !
2674        call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data
2675        call GW1RDM_me%free() ! Deallocate GW1RD_me
2676        ABI_FREE(old_ks_purex)
2677        ABI_FREE(new_hartr)
2678      else
2679        ABI_FREE(gw_rhor)
2680        ABI_FREE(bdm2_mask)
2681      end if
2682      ABI_FREE(nateigv)
2683    endif
2684 
2685    call xmpi_barrier(Wfd%comm)
2686 
2687   ! MRM: skip the rest for rdm_update=true (density matrix update)
2688   if (.not.rdm_update) then
2689    !  =====================================================
2690    !  ==== Solve Dyson equation storing results in Sr% ====
2691    !  =====================================================
2692    !  * Use perturbative approach or AC to find QP corrections.
2693    !  * If qp-GW, diagonalize also H0+Sigma in the KS basis set to get the
2694    !  new QP amplitudes and energies (Sr%eigvec_qp and Sr%en_qp_diago.
2695    !  TODO AC with spinor not implemented yet.
2696    !  TODO Diagonalization of Sigma+hhartre with AC is wrong.
2697    !
2698    call timab(425,1,tsec) ! solve_dyson
2699    do ikcalc=1,Sigp%nkptgw
2700      ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2701      ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2702      ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2703 
2704      ABI_MALLOC(sigcme_k,(nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab))
2705      sigcme_k=sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:)
2706 
2707      call solve_dyson(ikcalc,ib1,ib2,nomega_sigc,Sigp,Kmesh,sigcme_k,qp_ebands%eig,Sr,Dtset%prtvol,Dtfil,Wfd%comm)
2708      ABI_FREE(sigcme_k)
2709      !
2710      ! Calculate direct gap for each spin and print out final results.
2711      ! We use the valence index of the KS system because we still do not know
2712      ! all the QP corrections. Ideally one should use the QP valence index
2713      do spin=1,Sigp%nsppol
2714        if (Sigp%maxbnd(ikcalc,spin) >= ks_vbik(ik_ibz,spin)+1 .and. &
2715            Sigp%minbnd(ikcalc,spin) <= ks_vbik(ik_ibz,spin)  ) then
2716          ks_iv=ks_vbik(ik_ibz,spin)
2717          Sr%egwgap (ik_ibz,spin)= Sr%egw(ks_iv+1,ik_ibz,spin) -  Sr%egw(ks_iv,ik_ibz,spin)
2718          Sr%degwgap(ik_ibz,spin)= Sr%degw(ks_iv+1,ik_ibz,spin) - Sr%degw(ks_iv,ik_ibz,spin)
2719        else
2720          ! The "gap" cannot be computed
2721          Sr%e0gap(ik_ibz,spin)=zero; Sr%egwgap (ik_ibz,spin)=zero; Sr%degwgap(ik_ibz,spin)=zero
2722        end if
2723      end do
2724 
2725      if (wfd%my_rank == master) call write_sigma_results(ikcalc,ik_ibz,Sigp,Sr,ks_ebands)
2726    end do !ikcalc
2727     write(msg,'(a20)')'SIGMAS TO DO'
2728 
2729    call timab(425,2,tsec) ! solve_dyson
2730    call timab(426,1,tsec) ! finalize
2731 
2732    ! Update the energies in qp_ebands
2733    ! If QPSCGW, use diagonalized eigenvalues otherwise perturbative results.
2734    if (gwcalctyp>=10) then
2735      do ib=1,Sigp%nbnds
2736        qp_ebands%eig(ib,:,:)=Sr%en_qp_diago(ib,:,:)
2737      end do
2738    else
2739      qp_ebands%eig=Sr%egw
2740    end if
2741 
2742    ! ================================================================================
2743    ! ==== This part is done only if all k-points in the IBZ have been calculated ====
2744    ! ================================================================================
2745    if (Sigp%nkptgw==Kmesh%nibz) then
2746 
2747      ! Recalculate new occupations and Fermi level.
2748      call ebands_update_occ(qp_ebands,Dtset%spinmagntarget,prtvol=Dtset%prtvol)
2749      qp_vbik(:,:) = ebands_get_valence_idx(qp_ebands)
2750 
2751      write(msg,'(2a,3x,2(es16.6,a))')ch10,' New Fermi energy : ',qp_ebands%fermie,' Ha ,',qp_ebands%fermie*Ha_eV,' eV'
2752      call wrtout([std_out, ab_out], msg)
2753 
2754      ! === If all k-points and all occupied bands are calculated, output EXX ===
2755      ! FIXME here be careful about the check on  ks_vbik in case of metals
2756      ! if (my_rank==master.and.Sigp%nkptgw==Kmesh%nibz.and.ALL(Sigp%minbnd(:)==1).and.ALL(Sigp%maxbnd(:)>=MAXVAL(nbv(:)))) then
2757      ! if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=MAXVAL(MAXVAL(ks_vbik(:,:),DIM=1))) ) then
2758      if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=ks_vbik) ) then  ! FIXME here the two arrays use a different indexing.
2759 
2760        ex_energy = sigma_get_exene(Sr,Kmesh,qp_ebands)
2761        write(msg,'(a,2(es16.6,a))')' New Exchange energy : ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV'
2762        call wrtout([std_out, ab_out], msg)
2763      end if
2764 
2765      ! Report the QP gaps (Fundamental and direct)
2766      call ebands_report_gap(qp_ebands,header='QP Band Gaps',unit=ab_out)
2767 
2768      ! Band structure interpolation from QP energies computed on the k-mesh.
2769      if (nint(dtset%einterp(1)) /= 0 .and. all(sigp%minbdgw == sigp%minbnd) .and. all(sigp%maxbdgw == sigp%maxbnd)) then
2770        call ebands_interpolate_kpath(qp_ebands, dtset, cryst, [sigp%minbdgw, sigp%maxbdgw], dtfil%filnam_ds(4), comm)
2771      end if
2772    end if ! Sigp%nkptgw==Kmesh%nibz
2773    !
2774    ! Write SCF data in case of self-consistent calculation ===
2775    ! Save Sr%en_qp_diago, Sr%eigvec_qp and m_ks_to_qp in the _QPS file.
2776    ! Note that in the first iteration qp_rhor contains KS rhor, then the mixed rhor.
2777    if (gwcalctyp>=10) then
2778      ! Calculate the new m_ks_to_qp
2779      call updt_m_ks_to_qp(Sigp,Kmesh,nscf,Sr,Sr%m_ks_to_qp)
2780 
2781      if (wfd%my_rank == master) then
2782        call wrqps(Dtfil%fnameabo_qps,Sigp,Cryst,Kmesh,Psps,Pawtab,QP_Pawrhoij,&
2783          Dtset%nspden,nscf,nfftf,ngfftf,Sr,qp_ebands,Sr%m_ks_to_qp,qp_rhor)
2784      end if
2785 
2786      ! Report the MAX variation for each kptgw and spin
2787      call wrtout(ab_out,ch10//' Convergence of QP corrections ')
2788      do spin=1,Sigp%nsppol
2789        write(msg,'(a,i2,a)')' >>>>> For spin ',spin,' <<<<< '
2790        call wrtout(ab_out, msg)
2791        do ikcalc=1,Sigp%nkptgw
2792          ib1 = Sigp%minbnd(ikcalc,spin); ib2 = Sigp%maxbnd(ikcalc,spin)
2793          ik_bz = Sigp%kptgw2bz(ikcalc); ik_ibz = Kmesh%tab(ik_bz)
2794          ii = imax_loc(ABS(Sr%degw(ib1:ib2,ik_ibz,spin)))
2795          max_degw = Sr%degw(ii,ik_ibz,spin)
2796          write(msg,('(a,i3,a,2f8.3,a,i3)'))'.  kptgw no:',ikcalc,'; Maximum DeltaE = (',max_degw*Ha_eV,') for band index:',ii
2797          call wrtout(ab_out, msg)
2798        end do
2799      end do
2800    end if
2801 
2802    ! ============================================
2803    ! ==== Save the GW results in NETCDF file ====
2804    ! ============================================
2805    if (wfd%my_rank == master) then
2806      NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), '_SIGRES.nc'), xmpi_comm_self))
2807      NCF_CHECK(nctk_defnwrite_ivars(ncid, ["sigres_version"], [1]))
2808      NCF_CHECK(cryst%ncwrite(ncid))
2809      NCF_CHECK(ebands_ncwrite(ks_ebands, ncid))
2810      NCF_CHECK(sigma_ncwrite(Sigp, Er, Sr, ncid)) ! WARNING!! If gw1rdm>0 then Er is no longer present!!
2811      ! Add qp_rhor. Note that qp_rhor == ks_rhor if wavefunctions are not updated.
2812      !ncerr = nctk_write_datar("qp_rhor",path,ngfft,cplex,nfft,nspden,&
2813      !                          comm_fft,fftn3_distrib,ffti3_local,datar,action)
2814      NCF_CHECK(nf90_close(ncid))
2815    end if
2816   end if ! MRM skipped if GW density matrix update
2817  end if ! ucrpa
2818 
2819  !----------------------------- END OF THE CALCULATION ------------------------
2820  !
2821  !=====================
2822  !==== Close Files ====
2823  !=====================
2824  if (wfd%my_rank == master) then
2825    close(unt_gw )
2826    close(unt_gwdiag)
2827    close(unt_sig)
2828    close(unt_sgr)
2829    close(unt_sigc)
2830    if (mod10==SIG_GW_AC) close(unt_sgm)
2831  end if
2832  !
2833  !===============================================
2834  !==== Free arrays and local data structures ====
2835  !===============================================
2836  ABI_FREE(ks_vbik)
2837  ABI_FREE(qp_vbik)
2838  ABI_FREE(ph1d)
2839  ABI_FREE(ph1df)
2840  ABI_FREE(qp_rhor)
2841  ABI_FREE(ks_rhor)
2842  ABI_FREE(ks_rhog)
2843  ABI_FREE(qp_taur)
2844  ABI_FREE(ks_taur)
2845  ABI_FREE(ks_vhartr)
2846  ABI_FREE(ks_vtrial)
2847  ABI_FREE(vpsp)
2848  ABI_FREE(ks_vxc)
2849  ABI_FREE(xccc3d)
2850  ABI_FREE(grchempottn)
2851  ABI_FREE(grewtn)
2852  ABI_FREE(grvdw)
2853  if(.not.rdm_update) then
2854    ABI_FREE(sigcme)
2855  end if
2856 
2857  ABI_SFREE(kxc)
2858  ABI_SFREE(qp_vtrial)
2859 
2860  ABI_FREE(ks_nhat)
2861  ABI_FREE(ks_nhatgr)
2862  ABI_FREE(dijexc_core)
2863  ABI_FREE(ks_aepaw_rhor)
2864  call pawfgr_destroy(Pawfgr)
2865 
2866  ! Deallocation for PAW.
2867  if (Dtset%usepaw==1) then
2868    call pawrhoij_free(KS_Pawrhoij)
2869    ABI_FREE(KS_Pawrhoij)
2870    call pawfgrtab_free(Pawfgrtab)
2871    call paw_ij_free(KS_paw_ij)
2872    call paw_an_free(KS_paw_an)
2873    call pawpwff_free(Paw_pwff)
2874    if (gwcalctyp>=10) then
2875      call pawrhoij_free(QP_pawrhoij)
2876      call paw_ij_free(QP_paw_ij)
2877      call paw_an_free(QP_paw_an)
2878    end if
2879    if (Dtset%pawcross==1) then
2880      call paw_pwaves_lmn_free(Paw_onsite)
2881      call wfdf%free()
2882    end if
2883  end if
2884  ABI_FREE(Paw_onsite)
2885  ABI_FREE(Pawfgrtab)
2886  ABI_FREE(Paw_pwff)
2887  ABI_FREE(KS_paw_ij)
2888  ABI_FREE(KS_paw_an)
2889  if (gwcalctyp>=10) then
2890    ABI_FREE(QP_pawrhoij)
2891    ABI_FREE(QP_paw_an)
2892    ABI_FREE(QP_paw_ij)
2893  end if
2894 
2895  call wfd%free()
2896  call destroy_mpi_enreg(MPI_enreg_seq)
2897  call littlegroup_free(Ltg_k)
2898  ABI_FREE(Ltg_k)
2899  call Kmesh%free()
2900  call Qmesh%free()
2901  call Gsph_Max%free()
2902  call Gsph_x%free()
2903  call Gsph_c%free()
2904  call Vcp%free()
2905  call cryst%free()
2906  call sigma_free(Sr)
2907  if(.not.rdm_update) then
2908    call em1results_free(Er)
2909  end if
2910  call ppm_free(PPm)
2911  call Hdr_sigma%free()
2912  call Hdr_wfk%free()
2913  call ebands_free(ks_ebands)
2914  call ebands_free(qp_ebands)
2915  call KS_me%free()
2916  call esymm_free(KS_sym)
2917  ABI_FREE(KS_sym)
2918 
2919  if (Sigp%symsigma==1.and.gwcalctyp>=20) then
2920    call esymm_free(QP_sym)
2921    ABI_FREE(QP_sym)
2922  end if
2923 
2924  call Sigp%free()
2925 
2926  call timab(426,2,tsec) ! finalize
2927  call timab(401,2,tsec)
2928 
2929  DBG_EXIT('COLL')
2930 
2931 end subroutine sigma