TABLE OF CONTENTS


ABINIT/m_screening_driver [ Modules ]

[ Top ] [ Modules ]

NAME

  m_screening_driver

FUNCTION

 Calculate screening and dielectric functions

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

 21 #if defined HAVE_CONFIG_H
 22 #include "config.h"
 23 #endif
 24 
 25 #include "abi_common.h"
 26 
 27 module m_screening_driver
 28 
 29  use defs_basis
 30  use defs_datatypes
 31  use defs_abitypes
 32  use defs_wvltypes
 33  use m_abicore
 34  use m_xmpi
 35  use m_xomp
 36  use m_errors
 37  use m_ab7_mixing
 38  use m_kxc
 39  use m_nctk
 40 #ifdef HAVE_NETCDF
 41  use netcdf
 42 #endif
 43  use libxc_functionals
 44  use m_hdr
 45 
 46  use m_time,          only : timab
 47  use m_io_tools,      only : open_file, file_exists, iomode_from_fname
 48  use m_fstrings,      only : int2char10, sjoin, strcat, itoa, ltoa, itoa
 49  use m_energies,      only : energies_type, energies_init
 50  use m_numeric_tools, only : print_arr, iseven, coeffs_gausslegint
 51  use m_geometry,      only : normv, vdotw, mkrdim, metric
 52  use m_gwdefs,        only : GW_TOLQ0, GW_TOLQ, em1params_free, em1params_t, GW_Q0_DEFAULT
 53  use m_mpinfo,        only : destroy_mpi_enreg, initmpi_seq
 54  use m_crystal,       only : crystal_free, crystal_t, crystal_print
 55  use m_crystal_io,    only : crystal_ncwrite, crystal_from_hdr
 56  use m_ebands,        only : ebands_update_occ, ebands_copy, get_valence_idx, get_occupied, apply_scissor, &
 57                              ebands_free, ebands_has_metal_scheme, ebands_ncwrite, ebands_init
 58  use m_bz_mesh,       only : kmesh_t, kmesh_init, kmesh_free, littlegroup_t, littlegroup_free, littlegroup_init, &
 59                              get_ng0sh, kmesh_print, find_qmesh, get_BZ_item
 60  use m_kg,            only : getph
 61  use m_gsphere,       only : gsph_free, gsphere_t, gsph_init, merge_and_sort_kg, setshells
 62  use m_vcoul,         only : vcoul_t, vcoul_free, vcoul_init
 63  use m_qparticles,    only : rdqps, rdgw, show_QP
 64  use m_screening,     only : make_epsm1_driver, lwl_write, chi_t, chi_free, chi_new
 65  use m_io_screening,  only : hscr_new, hscr_io, write_screening, hscr_free, hscr_t
 66  use m_spectra,       only : spectra_t, spectra_write, spectra_repr, spectra_free, W_EM_LF, W_EM_NLF, W_EELF
 67  use m_fftcore,       only : print_ngfft
 68  use m_fft_mesh,      only : rotate_FFT_mesh, cigfft, get_gftt, setmesh
 69  use m_fft,           only : fourdp
 70  use m_wfd,           only : wfd_init, wfd_free,  wfd_nullify, wfd_print, wfd_t, wfd_rotate, wfd_test_ortho,&
 71                              wfd_read_wfk, wfd_test_ortho, wfd_copy, wfd_change_ngfft, wfd_mkrho, test_charge
 72  use m_wfk,           only : wfk_read_eigenvalues
 73  use m_io_kss,        only : make_gvec_kss
 74  use m_chi0tk,        only : output_chi0sumrule
 75  use m_pawang,        only : pawang_type
 76  use m_pawrad,        only : pawrad_type
 77  use m_pawtab,        only : pawtab_type, pawtab_print, pawtab_get_lsize
 78  use m_paw_an,        only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
 79  use m_paw_ij,        only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
 80  use m_pawfgrtab,     only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
 81  use m_pawrhoij,      only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy,&
 82 &                            pawrhoij_free, symrhoij, pawrhoij_get_nspden
 83  use m_pawdij,        only : pawdij, symdij_all
 84  use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
 85  use m_pawpwij,       only : pawpwff_t, pawpwff_init, pawpwff_free
 86  use m_pawfgr,        only : pawfgr_type, pawfgr_init, pawfgr_destroy
 87  use m_paw_sphharm,   only : setsym_ylm
 88  use m_paw_onsite,    only : pawnabla_init
 89  use m_paw_nhat,      only : nhatgrid,pawmknhat
 90  use m_paw_denpot,    only : pawdenpot
 91  use m_paw_init,      only : pawinit,paw_gencond
 92  use m_paw_tools,     only : chkpawovlp,pawprt
 93  use m_chi0,          only : cchi0, cchi0q0, chi0q0_intraband
 94  use m_setvtr,        only : setvtr
 95  use m_mkrho,         only : prtrhomxmn
 96  use m_pspini,        only : pspini
 97 
 98  use m_paw_correlations, only : pawpuxinit
 99 
100  implicit none
101 
102  private

m_screening_driver/calc_rpa_functional [ Functions ]

[ Top ] [ m_screening_driver ] [ Functions ]

NAME

 calc_rpa_functional

FUNCTION

  Routine used to calculate the RPA approximation to the correlation energy
  from the irreducible polarizability.

INPUTS

  iq=index of the q-point in the array Qmesh%ibz where epsilon^-1 has to be calculated
  Ep<em1params_t>=Structure with parameters and dimensions related to the inverse dielectric matrix.
  Pvc<vcoul_t>=Structure gathering data on the Coulombian interaction
  Qmesh<kmesh_t>=Data type with information on the q-sampling
  Dtfil<Datafiles_type)>=variables related to files
  spaceComm=MPI communicator.

OUTPUT

PARENTS

      screening

CHILDREN

      coeffs_gausslegint,wrtout,xginv,xheev,xmpi_sum_master

SOURCE

2615 subroutine calc_rpa_functional(gwrpacorr,iqcalc,iq,Ep,Pvc,Qmesh,Dtfil,gmet,chi0,spaceComm,ec_rpa)
2616 
2617  use m_hide_lapack, only : xginv, xheev
2618 
2619 !This section has been created automatically by the script Abilint (TD).
2620 !Do not modify the following lines by hand.
2621 #undef ABI_FUNC
2622 #define ABI_FUNC 'calc_rpa_functional'
2623 !End of the abilint section
2624 
2625  implicit none
2626 
2627 !Arguments ------------------------------------
2628 !scalars
2629  integer,intent(in) :: iqcalc,iq,gwrpacorr,spaceComm
2630  type(kmesh_t),intent(in) :: Qmesh
2631  type(vcoul_t),intent(in) :: Pvc
2632  type(Datafiles_type),intent(in) :: Dtfil
2633  type(em1params_t),intent(in) :: Ep
2634 !arrays
2635  real(dp),intent(in) :: gmet(3,3)
2636  real(dp),intent(inout) :: ec_rpa(gwrpacorr)
2637  complex(gwpc),intent(inout) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega)
2638 
2639 !Local variables-------------------------------
2640 !scalars
2641  integer :: ig1,ig2,ilambda,io,master,rank,nprocs,unt,ierr
2642  real(dp) :: ecorr
2643  real(dp) :: lambda
2644  logical :: qeq0
2645  character(len=500) :: msg
2646 !arrays
2647  real(dp),allocatable :: z(:),zl(:),zlw(:),zw(:)
2648  complex(gwpc),allocatable :: chi0_diag(:),chitmp(:,:)
2649  real(gwpc),allocatable :: eig(:)
2650 
2651 ! *************************************************************************
2652 
2653  DBG_ENTER("COLL")
2654 
2655 ! initialize MPI data
2656  master=0
2657  rank   = xmpi_comm_rank(spaceComm)
2658  nprocs = xmpi_comm_size(spaceComm)
2659 
2660  !if (rank==master) then ! presently only master has chi0 in screening
2661 
2662  ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff
2663  qeq0=(normv(Qmesh%ibz(:,iq),gmet,'G')<GW_TOLQ0)
2664 
2665  ! Calculate Gauss-Legendre quadrature knots and weights for the omega integration
2666  ABI_ALLOCATE(zw,(Ep%nomegaei))
2667  ABI_ALLOCATE(z,(Ep%nomegaei))
2668  call coeffs_gausslegint(zero,one,z,zw,Ep%nomegaei)
2669 
2670  ! Calculate Gauss-Legendre quadrature knots and weights for the lambda integration
2671  ABI_ALLOCATE(zlw,(gwrpacorr))
2672  ABI_ALLOCATE(zl,(gwrpacorr))
2673  call coeffs_gausslegint(zero,one,zl,zlw,gwrpacorr)
2674 
2675 
2676  ABI_ALLOCATE(chi0_diag,(Ep%npwe))
2677  ABI_STAT_ALLOCATE(chitmp,(Ep%npwe,Ep%npwe), ierr)
2678  ABI_CHECK(ierr==0, "out-of-memory in chitmp")
2679 
2680  do io=2,Ep%nomega
2681 
2682    if(gwrpacorr==1) then ! exact integration over the coupling constant
2683 
2684      if(modulo(io-2,nprocs)/=rank) cycle ! distributing the workload
2685 
2686      do ig2=1,Ep%npwe
2687        do ig1=1,Ep%npwe
2688          chitmp(ig1,ig2) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig2,iq) * chi0(ig1,ig2,io)
2689        end do !ig1
2690      end do !ig2
2691      ABI_ALLOCATE(eig,(Ep%npwe))
2692      call xheev('V','U',Ep%npwe,chitmp,eig)
2693 
2694      do ig1=1,Ep%npwe
2695        ec_rpa(:) = ec_rpa(:) &
2696 &         - zw(io-1) / ( z(io-1) * z(io-1) ) &
2697 &              * Qmesh%wt(iq) * (-log( 1.0_dp-eig(ig1) )  - eig(ig1) ) / (2.0_dp * pi )
2698      end do
2699      ABI_DEALLOCATE(eig)
2700 
2701    else ! numerical integration over the coupling constant
2702 
2703       if(modulo( (ilambda-1)+gwrpacorr*(io-2),nprocs)/=rank) cycle ! distributing the workload
2704 
2705      do ilambda=1,gwrpacorr
2706        lambda=zl(ilambda)
2707        do ig1=1,Ep%npwe
2708          chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq)**2 * chi0(ig1,ig1,io)
2709        end do
2710 
2711        do ig2=1,Ep%npwe
2712          do ig1=1,Ep%npwe
2713            chitmp(ig1,ig2) = - lambda * Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chi0(ig1,ig2,io)
2714          end do !ig1
2715          chitmp(ig2,ig2) = chitmp(ig2,ig2) + 1.0_dp
2716        end do !ig2
2717        call xginv(chitmp(:,:),Ep%npwe)
2718        chitmp(:,:) = matmul( chi0(:,:,io) , chitmp(:,:) )
2719 
2720        do ig1=1,Ep%npwe
2721          chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chitmp(ig1,ig1) - chi0_diag(ig1)
2722        end do
2723 
2724        do ig1=1,Ep%npwe
2725          ec_rpa(ilambda) = ec_rpa(ilambda) &
2726 &           - zw(io-1) / ( z(io-1) * z(io-1) ) * Qmesh%wt(iq) * real(  chi0_diag(ig1) ) / (2.0_dp * pi )
2727        end do
2728 
2729      end do ! ilambda
2730 
2731    end if ! exact or numerical integration over the coupling constant
2732 
2733  end do ! io
2734 
2735 
2736  ! Output the correlation energy when the last q-point to be calculated is reached
2737  ! This would allow for a manual parallelization over q-points
2738  if(iqcalc==Ep%nqcalc) then
2739 
2740    call xmpi_sum_master(ec_rpa,master,spaceComm,ierr)
2741 
2742    if(rank==master) then
2743      ecorr = sum( zlw(:)*ec_rpa(:) )
2744      if (open_file(dtfil%fnameabo_rpa, msg, newunit=unt) /=0) then
2745        MSG_ERROR(msg)
2746      end if
2747      write(unt,'(a,(2x,f14.8))') '#RPA',ecorr
2748      write(msg,'(2a,(2x,f14.8))') ch10,' RPA energy [Ha] :',ecorr
2749      call wrtout(std_out,msg,'COLL')
2750      call wrtout(ab_out,msg,'COLL')
2751      if(gwrpacorr>1) then
2752        do ilambda=1,gwrpacorr
2753          write(unt,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda)
2754          write(msg,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda)
2755          call wrtout(std_out,msg,'COLL')
2756          call wrtout(ab_out,msg,'COLL')
2757        end do
2758      end if
2759      close(unt)
2760    end if
2761 
2762  end if
2763 
2764  ABI_DEALLOCATE(chi0_diag)
2765  ABI_DEALLOCATE(chitmp)
2766  ABI_DEALLOCATE(zl)
2767  ABI_DEALLOCATE(zlw)
2768  ABI_DEALLOCATE(z)
2769  ABI_DEALLOCATE(zw)
2770 
2771  DBG_EXIT("COLL")
2772 
2773 end subroutine calc_rpa_functional

m_screening_driver/chi0_bksmask [ Functions ]

[ Top ] [ m_screening_driver ] [ Functions ]

NAME

 chi0_bksmask

FUNCTION

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

INPUTS

 Dtset<type(dataset_type)>=all input variables for this dataset
 Ep<em1params_t>=Parameters for the screening calculation.
 Kmesh <kmesh_t>=Structure describing the k-point sampling.
 nbvw = Max. number of fully/partially occupied states over spin
 nbcw = Max. number of unoccupied states considering the spin
 nprocs=Total number of MPI processors
 my_rank=Rank of this this processor.

OUTPUT

 bks_mask(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state.
 keep_ur(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space.
 ierr=Exit status.

PARENTS

      screening

CHILDREN

      xmpi_split_work2_i4b

SOURCE

2237 subroutine chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr)
2238 
2239 
2240 !This section has been created automatically by the script Abilint (TD).
2241 !Do not modify the following lines by hand.
2242 #undef ABI_FUNC
2243 #define ABI_FUNC 'chi0_bksmask'
2244 !End of the abilint section
2245 
2246  implicit none
2247 
2248 !Arguments ------------------------------------
2249 !scalars
2250  integer,intent(in) :: my_rank,nprocs,nbvw,nbcw
2251  integer,intent(out) :: ierr
2252  type(Dataset_type),intent(in) :: Dtset
2253  type(em1params_t),intent(in) :: Ep
2254  type(kmesh_t),intent(in) :: Kmesh
2255 !arrays
2256  logical,intent(out) :: bks_mask(Ep%nbnds,Kmesh%nibz,Dtset%nsppol)
2257  logical,intent(out) :: keep_ur(Ep%nbnds,Kmesh%nibz,Dtset%nsppol)
2258 
2259 !Local variables-------------------------------
2260 !scalars
2261  integer :: my_nspins,my_maxb,my_minb,isp,spin,distr_err,nsppol,band,rank_spin,ib
2262  character(len=500) :: msg
2263  logical :: store_ur
2264 !arrays
2265  integer :: my_spins(Dtset%nsppol),nprocs_spin(Dtset%nsppol)
2266  integer,allocatable :: istart(:),istop(:)
2267 
2268 ! *************************************************************************
2269 
2270  ierr=0; nsppol=Dtset%nsppol
2271 
2272  my_nspins=Dtset%nsppol; my_spins= [(isp,isp=1,nsppol)]
2273 
2274  ! List of spins for each node, number of processors per each spin
2275  ! and the MPI rank in the "spin" communicator.
2276  nprocs_spin = nprocs; rank_spin = my_rank
2277 
2278  if (nsppol==2.and.nprocs>1) then
2279    ! Distribute spins (optimal distribution if nprocs is even)
2280    nprocs_spin(1) = nprocs/2
2281    nprocs_spin(2) = nprocs - nprocs/2
2282    my_nspins=1; my_spins(1)=1
2283    if (my_rank+1>nprocs/2) then
2284      my_spins(1)=2
2285      rank_spin = my_rank - nprocs/2
2286    end if
2287  end if
2288 
2289  store_ur = (MODULO(Dtset%gwmem,10)==1)
2290  bks_mask=.FALSE.; keep_ur=.FALSE.
2291 
2292  select case (Dtset%gwpara)
2293  case (1)
2294    ! Parallelization over transitions **without** memory distributions (Except for the spin).
2295    my_minb=1; my_maxb=Ep%nbnds
2296    do isp=1,my_nspins
2297      spin = my_spins(isp)
2298      bks_mask(my_minb:my_maxb,:,spin) = .TRUE.
2299      if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
2300    end do
2301 
2302  case (2)
2303    ! Distribute bands and spin.
2304    do isp=1,my_nspins
2305      spin = my_spins(isp)
2306 
2307      if (nprocs_spin(spin) <= nbcw) then
2308        ! Distribute nbcw empty bands among nprocs_spin (block of bands without replicas).
2309        ! Bands are distributed in contiguous blocks because
2310        ! this distribution is well suited for the Hilber transform
2311        ! since each node will allocate only a smaller frequency interval
2312        ! for the spectral function whose size scales with the number of MPI nodes.
2313        ! Note it is now meaningless to distinguish gwcomp=0 or 1 since the workload is well balanced later on
2314        ABI_MALLOC(istart,(nprocs_spin(spin)))
2315        ABI_MALLOC(istop,(nprocs_spin(spin)))
2316 
2317        call xmpi_split_work2_i4b(nbcw,nprocs_spin(spin),istart,istop,msg,distr_err)
2318 
2319        if (distr_err==2) then
2320          ! Too many processors.
2321          !MSG_WARNING(msg)
2322          ierr=1
2323        end if
2324 
2325        my_minb = nbvw + istart(rank_spin+1)
2326        my_maxb = nbvw + istop (rank_spin+1)
2327 
2328        ABI_FREE(istart)
2329        ABI_FREE(istop)
2330 
2331        if (my_maxb-my_minb+1<=0) then
2332          write(msg,'(3a,2(i0,a),2a)')&
2333 &         'One or more processors has zero number of bands ',ch10,&
2334 &         'my_minb= ',my_minb,' my_maxb= ',my_maxb,ch10,&
2335 &         'This is a waste, decrease the number of processors.'
2336          MSG_ERROR(msg)
2337        end if
2338 
2339        bks_mask(my_minb:my_maxb,:,spin)=.TRUE.
2340        if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
2341 
2342      else
2343        ! New version (alternate bands with replicas if nprocs > nbcw)
2344        ! FIXME: Fix segmentation fault with Hilbert transform.
2345        do ib=1,nbcw
2346          if (xmpi_distrib_with_replicas(ib,nbcw,rank_spin,nprocs_spin(spin))) then
2347            band = ib + nbvw
2348            bks_mask(band,:,spin)=.TRUE.
2349            if (store_ur) keep_ur(band,:,spin)=.TRUE.
2350          end if
2351        end do
2352      end if
2353 
2354      ! This is needed to have all the occupied states on each node.
2355      bks_mask(1:nbvw,:,spin) = .TRUE.
2356      if (store_ur) keep_ur(1:nbvw,:,spin)=.TRUE.
2357    end do ! isp
2358 
2359  case default
2360    ierr = 1
2361    MSG_WARNING("Wrong value for gwpara")
2362  end select
2363 
2364 end subroutine chi0_bksmask

m_screening_driver/random_stopping_power [ Functions ]

[ Top ] [ m_screening_driver ] [ Functions ]

NAME

 random_stopping_power

FUNCTION

 Calculate the electronic random stopping power

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      screening

CHILDREN

      get_bz_item,spline,splint,wrtout

SOURCE

2388 subroutine random_stopping_power(iqibz,npvel,pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower)
2389 
2390  use m_splines
2391 
2392 !This section has been created automatically by the script Abilint (TD).
2393 !Do not modify the following lines by hand.
2394 #undef ABI_FUNC
2395 #define ABI_FUNC 'random_stopping_power'
2396 !End of the abilint section
2397 
2398  implicit  none
2399 
2400 !Arguments ------------------------------------
2401 !scalars
2402  integer,intent(in)                    :: iqibz,npvel
2403  real(dp),intent(in)                   :: pvelmax(3)
2404 
2405  type(em1params_t),intent(in) :: Ep
2406  type(gsphere_t),intent(in)            :: Gsph_epsG0
2407  type(kmesh_t),intent(in)              :: Qmesh
2408  type(vcoul_t),intent(in)              :: Vcp
2409  type(crystal_t),intent(in)            :: Cryst
2410  type(Datafiles_type),intent(in)       :: Dtfil
2411 
2412  complex(gwpc),intent(in)              :: epsm1(Ep%npwe,Ep%npwe,Ep%nomega)
2413 
2414  real(dp),intent(inout)                :: rspower(npvel)
2415 
2416 !Local variables ------------------------------
2417  integer :: ipvel,ig
2418  integer :: iq_bz,iq_ibz,isym_q,itim_q
2419  integer :: iomega,iomegap,nomega_re
2420  integer :: unt_rsp
2421  integer,allocatable :: iomega_re(:)
2422 
2423  real(dp),parameter :: zp=1.0_dp              ! Hard-coded charge of the impinging particle
2424  real(dp) :: omega_p
2425  real(dp) :: im_epsm1_int(1)
2426  real(dp) :: qbz(3),qpgcart(3),qpgred(3)
2427  real(dp) :: pvel(3,npvel),pvel_norm(npvel)
2428  real(dp) :: ypp_i(Ep%nomega)
2429  real(dp) :: vcoul(Ep%npwe)
2430  real(dp),allocatable :: im_epsm1_diag_qbz(:,:),tmp_data(:)
2431  real(dp),allocatable :: omega_re(:)
2432 
2433  character(len=500)     :: msg
2434  character(len=fnlen+4) :: fname
2435 
2436 !************************************************************************
2437 
2438  !
2439  ! First set up the velocities array from the input variables npvel and pvelmax(3)
2440  ! Remember pvelmax is in Cartesian coordinates and so is pvel
2441  do ipvel=1,npvel
2442    pvel(:,ipvel)    = REAL(ipvel,dp) / REAL(npvel,dp) * pvelmax(:)
2443    pvel_norm(ipvel) = SQRT( SUM( pvel(:,ipvel)**2 ) )
2444  enddo
2445  !
2446  ! Select the purely real frequency in Ep%omega
2447  nomega_re=0
2448  do iomega=1,Ep%nomega
2449    if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then
2450      nomega_re=nomega_re+1
2451    endif
2452  enddo
2453  ABI_ALLOCATE(omega_re,(nomega_re))
2454  ABI_ALLOCATE(iomega_re,(nomega_re))
2455  ABI_ALLOCATE(im_epsm1_diag_qbz,(Ep%npwe,Ep%nomega))
2456  ABI_ALLOCATE(tmp_data,(Ep%nomega))
2457 
2458  iomegap=0
2459  do iomega=1,Ep%nomega
2460    if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then
2461      iomegap=iomegap+1
2462      iomega_re(iomegap)=iomega
2463      omega_re(iomegap)=REAL(Ep%omega(iomega),dp)
2464    endif
2465  enddo
2466 
2467  !
2468  ! Loop over all the q-points in the full Brillouin zone and select only the
2469  ! ones that corresponds to the correct q-point in the irreducible wedge we are
2470  ! currently treating (index iqibz)
2471  do iq_bz=1,Qmesh%nbz
2472 
2473    ! Perform the check and obtain the symmetry information
2474    call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q)
2475    if( iqibz /= iq_ibz ) cycle
2476 
2477    ! Apply the symmetry operation to the diagonal of epsm1
2478    do iomega=1,nomega_re
2479      do ig=1,Ep%npwe
2480        im_epsm1_diag_qbz(Gsph_epsG0%rottb(ig,itim_q,isym_q),iomega)= AIMAG( epsm1(ig,ig,iomega_re(iomega)) )
2481      enddo
2482    enddo
2483    ! Apply the symmetry operation to the Coulomb interaction
2484    do ig=1,Ep%npwe
2485      vcoul(Gsph_epsG0%rottb(ig,itim_q,isym_q))=Vcp%vc_sqrt(ig,iqibz)**2
2486    enddo
2487 
2488    !
2489    ! Sum over G vectors
2490    do ig=1,Ep%npwe
2491      !
2492      ! Loop over velocities
2493      do ipvel=1,npvel
2494 
2495        qpgred(:) = qbz(:) + Gsph_epsG0%gvec(:,ig)
2496        ! Transform q + G from reduced to cartesian with the symmetry operation
2497        qpgcart(:) = two_pi * Cryst%gprimd(:,1) * qpgred(1) &
2498 &                 + two_pi * Cryst%gprimd(:,2) * qpgred(2) &
2499 &                 + two_pi * Cryst%gprimd(:,3) * qpgred(3)
2500 
2501        ! omega_p = ( q + G ) . v
2502        omega_p =  DOT_PRODUCT( qpgcart(:) , pvel(:,ipvel) )
2503 
2504        ! Check that the calculated frequency omega_p is within the omega
2505        ! range of epsm1 and thus that the interpolation will go fine
2506        if ( ABS(omega_p) > omega_re(nomega_re) ) then
2507          write(msg,'(a,e16.4,2a,e16.4)') ' freqremax is currently ',omega_re(nomega_re),ch10,&
2508 &                                        ' increase it to at least ',omega_p
2509          MSG_WARNING(msg)
2510        endif
2511 
2512        ! Perform the spline interpolation to obtain epsm1 at the desired omega = omega_p
2513        tmp_data = im_epsm1_diag_qbz(ig,:)
2514 
2515        call spline( omega_re, tmp_data, nomega_re, 1.0e+32_dp, 1.0e+32_dp, ypp_i)
2516        call splint( nomega_re, omega_re, tmp_data, ypp_i, 1, (/ ABS(omega_p) /),  im_epsm1_int )
2517 
2518        !
2519        ! Apply the odd parity of Im epsm1 in  omega to recover the causal response function
2520        if (omega_p<zero) im_epsm1_int(1)=-im_epsm1_int(1)
2521 
2522        !
2523        ! Calculate 4 * pi / |q+G|**2 * omega_p * Im{ epsm1_GG(q,omega_p) }
2524        !
2525        im_epsm1_int(1) = omega_p * vcoul(ig) * im_epsm1_int(1)
2526 
2527        ! Accumulate the final result without the prefactor
2528        ! (It will be included at the very end)
2529        rspower(ipvel) = rspower(ipvel) + im_epsm1_int(1)
2530 
2531      end do ! end of velocity loop
2532    end do ! end G-loop
2533 
2534  enddo ! end of q loop in the full BZ
2535 
2536  ! If it is the last q, write down the result in the main output file and in a
2537  ! separate file _RSP (for Random Stopping Power)
2538  if (iqibz == Qmesh%nibz ) then
2539 
2540    ! Multiply by the prefactors
2541    ! Note that this expression differs from Eq. (3.11) in Campillo PRB 58, 10307 (1998) [[cite:Campillo1998]].
2542    ! A factor one half is missing in the paper.
2543    rspower(:) = - zp**2 / ( Cryst%ucvol * Qmesh%nbz * pvel_norm(:) ) * rspower(:)
2544 
2545    write(msg,'(2a)')         ch10,' ==== Random stopping power along Cartesian direction  === '
2546    call wrtout(std_out,msg,'COLL')
2547    call wrtout(ab_out,msg,'COLL')
2548    write(msg,'(a,3(f12.4,2x),a)') ' ====  ',pvelmax(:),'===='
2549    call wrtout(std_out,msg,'COLL')
2550    call wrtout(ab_out,msg,'COLL')
2551    write(msg,'(a)')               '#  |v| (a.u.) , RSP (a.u.) '
2552    call wrtout(std_out,msg,'COLL')
2553    call wrtout(ab_out,msg,'COLL')
2554    do ipvel=1,npvel
2555      write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel)
2556      call wrtout(std_out,msg,'COLL')
2557      call wrtout(ab_out,msg,'COLL')
2558    enddo
2559    write(msg,'(2a)')              ' ========================================================= ',ch10
2560    call wrtout(std_out,msg,'COLL')
2561    call wrtout(ab_out,msg,'COLL')
2562 
2563    fname=TRIM(Dtfil%filnam_ds(4))//'_RSP'
2564    if (open_file(fname,msg,newunit=unt_rsp,status='unknown',form='formatted') /= 0) then
2565      MSG_ERROR(msg)
2566    end if
2567 
2568    write(msg,'(a)')               '# ==== Random stopping power along Cartesian direction  === '
2569    call wrtout(unt_rsp,msg,'COLL')
2570    write(msg,'(a,3(f12.4,2x))')   '# ====  ',pvelmax(:)
2571    call wrtout(unt_rsp,msg,'COLL')
2572    write(msg,'(a)')               '#  |v| (a.u.) , RSP (a.u.) '
2573    call wrtout(unt_rsp,msg,'COLL')
2574    do ipvel=1,npvel
2575      write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel)
2576      call wrtout(unt_rsp,msg,'COLL')
2577    enddo
2578    close(unt_rsp)
2579  end if
2580 
2581  ABI_DEALLOCATE(omega_re)
2582  ABI_DEALLOCATE(iomega_re)
2583  ABI_DEALLOCATE(im_epsm1_diag_qbz)
2584  ABI_DEALLOCATE(tmp_data)
2585 
2586 end subroutine random_stopping_power

m_screening_driver/screening [ Functions ]

[ Top ] [ m_screening_driver ] [ Functions ]

NAME

 screening

FUNCTION

 Calculate screening and dielectric functions

INPUTS

 acell(3)=length scales of primitive translations (bohr)
 codvsn=code version
 Dtfil<datafiles_type)>=variables related to file names and unit numbers.
 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
 Psps <type(pseudopotential_type)>=variables related to pseudopotentials
  Before entering the first time in screening, 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 output file.
 The symmetrical inverse dielectric matrix is stored in the _SCR file

SIDE EFFECTS

  Dtset<type(dataset_type)>=all input variables for this dataset

NOTES

 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.

PARENTS

      driver

CHILDREN

      apply_scissor,calc_rpa_functional,cchi0,cchi0q0,chi0_bksmask
      chi0q0_intraband,chi_free,chkpawovlp,coeffs_gausslegint,crystal_free
      destroy_mpi_enreg,ebands_copy,ebands_free,ebands_update_occ
      em1params_free,energies_init,fourdp,get_gftt,getph,gsph_free,hdr_free
      hscr_free,hscr_io,init_distribfft_seq,initmpi_seq,kmesh_free,kxc_ada
      kxc_driver,littlegroup_free,lwl_write,make_epsm1_driver,metric,mkrdim
      nhatgrid,output_chi0sumrule,paw_an_free,paw_an_init,paw_an_nullify
      paw_gencond,paw_ij_free,paw_ij_init,paw_ij_nullify,paw_pwaves_lmn_free
      paw_pwaves_lmn_init,pawdenpot,pawdij,pawfgr_destroy,pawfgr_init
      pawfgrtab_free,pawfgrtab_init,pawinit,pawmknhat,pawnabla_init,pawprt
      pawpuxinit,pawpwff_free,pawpwff_init,pawrhoij_alloc,pawrhoij_copy
      pawrhoij_free,pawtab_get_lsize,pawtab_print,print_arr,print_ngfft
      prtrhomxmn,pspini,random_stopping_power,rdgw,rdqps,rotate_fft_mesh
      setsym_ylm,setup_screening,setvtr,spectra_free,spectra_repr
      spectra_write,symdij,symdij_all,test_charge,timab,vcoul_free
      wfd_change_ngfft,wfd_copy,wfd_free,wfd_init,wfd_mkrho,wfd_print
      wfd_read_wfk,wfd_rotate,wfd_test_ortho,write_screening,wrtout
      xmpi_bcast

SOURCE

 181 subroutine screening(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim)
 182 
 183 
 184 !This section has been created automatically by the script Abilint (TD).
 185 !Do not modify the following lines by hand.
 186 #undef ABI_FUNC
 187 #define ABI_FUNC 'screening'
 188 !End of the abilint section
 189 
 190  implicit none
 191 
 192 !Arguments ------------------------------------
 193 !scalars
 194  character(len=6),intent(in) :: codvsn
 195  type(Datafiles_type),intent(in) :: Dtfil
 196  type(Dataset_type),intent(inout) :: Dtset
 197  type(Pawang_type),intent(inout) :: Pawang
 198  type(Pseudopotential_type),intent(inout) :: Psps
 199 !arrays
 200  real(dp),intent(in) :: acell(3),rprim(3,3)
 201  type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Dtset%usepaw)
 202  type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Dtset%usepaw)
 203 
 204 !Local variables ------------------------------
 205  character(len=4) :: ctype='RPA '
 206 !scalars
 207  integer,parameter :: tim_fourdp4=4,NOMEGA_PRINTED=15,master=0
 208  integer :: spin,ik_ibz,my_nbks
 209  integer :: choice,cplex,dim_kxcg,dim_wing,ount,omp_ncpus
 210  integer :: fform_chi0,fform_em1,gnt_option,iat,ider,idir,ierr,band
 211  integer :: ifft,ii,ikbz,ikxc,initialized,iomega,ios,ipert
 212  integer :: iqibz,iqcalc,is_qeq0,isym,izero,ifirst,ilast
 213  integer :: label,mgfftf,mgfftgw
 214  integer :: nt_per_proc,work_size
 215  integer :: moved_atm_inside,moved_rhor
 216  integer :: nbcw,nbsc,nbvw,nkxc,nkxc1,n3xccc,optene,istep
 217  integer :: nfftf,nfftf_tot,nfftgw,nfftgw_tot,ngrvdw,nhatgrdim,nprocs,nspden_rhoij
 218  integer :: nscf,nzlmopt,mband
 219  integer :: optcut,optgr0,optgr1,optgr2,option,approx_type,option_test,optgrad
 220  integer :: optrad,optrhoij,psp_gencond,my_rank
 221  integer :: rhoxsp_method,comm,test_type,tordering,unt_em1,unt_susc,usexcnhat
 222  real(dp) :: compch_fft,compch_sph,domegareal,e0,ecore,ecut_eff,ecutdg_eff
 223  real(dp) :: gsqcutc_eff,gsqcutf_eff,gsqcut_shp,omegaplasma,ucvol,vxcavg,gw_gsq,r_s
 224  real(dp) :: alpha,rhoav,opt_ecut,factor
 225  real(dp):: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem
 226  logical :: found,iscompatibleFFT,use_tr,is_first_qcalc
 227  logical :: add_chi0_intraband,update_energies,call_pawinit
 228  character(len=10) :: string
 229  character(len=500) :: msg
 230  character(len=80) :: bar
 231  type(ebands_t) :: KS_BSt,QP_BSt
 232  type(kmesh_t) :: Kmesh,Qmesh
 233  type(vcoul_t) :: Vcp
 234  type(crystal_t) :: Cryst
 235  type(em1params_t) :: Ep
 236  type(Energies_type) :: KS_energies
 237  type(gsphere_t) :: Gsph_epsG0,Gsph_wfn
 238  type(Hdr_type) :: Hdr_wfk,Hdr_local
 239  type(MPI_type) :: MPI_enreg_seq
 240  type(Pawfgr_type) :: Pawfgr
 241  type(hscr_t) :: Hem1,Hchi0
 242  type(wfd_t) :: Wfd,Wfdf
 243  type(spectra_t) :: spectra
 244  type(chi_t) :: chihw
 245  type(wvl_data) :: wvl_dummy
 246 !arrays
 247  integer :: ibocc(Dtset%nsppol),ngfft_gw(18),ngfftc(18),ngfftf(18)
 248  integer,allocatable :: irottb(:,:),ktabr(:,:),ktabrf(:,:),l_size_atm(:)
 249  integer,allocatable :: ks_vbik(:,:),ks_occ_idx(:,:),qp_vbik(:,:),nband(:,:)
 250  integer,allocatable :: nq_spl(:),nlmn_atm(:),gw_gfft(:,:)
 251  real(dp) :: gmet(3,3),gprimd(3,3),k0(3),qtmp(3),rmet(3,3),rprimd(3,3),tsec(2),strsxc(6)
 252  real(dp),allocatable :: igwene(:,:,:),chi0_sumrule(:),ec_rpa(:),rspower(:)
 253  real(dp),allocatable :: nhat(:,:),nhatgr(:,:,:),ph1d(:,:),ph1df(:,:)
 254  real(dp),allocatable :: rhog(:,:),rhor(:,:),rhor_p(:,:),rhor_kernel(:,:),taug(:,:),taur(:,:)
 255  real(dp),allocatable :: z(:),zw(:),grchempottn(:,:),grewtn(:,:),grvdw(:,:),kxc(:,:),qmax(:)
 256  real(dp),allocatable :: ks_vhartr(:),vpsp(:),ks_vtrial(:,:),ks_vxc(:,:),xccc3d(:)
 257  complex(gwpc),allocatable :: arr_99(:,:),kxcg(:,:),fxc_ADA(:,:,:)
 258  complex(dpc),allocatable :: m_lda_to_qp(:,:,:,:)
 259  complex(dpc),allocatable :: chi0_lwing(:,:,:),chi0_uwing(:,:,:),chi0_head(:,:,:)
 260  complex(dpc),allocatable :: chi0intra_lwing(:,:,:),chi0intra_uwing(:,:,:),chi0intra_head(:,:,:)
 261  complex(gwpc),allocatable,target :: chi0(:,:,:),chi0intra(:,:,:)
 262  complex(gwpc),ABI_CONTIGUOUS pointer :: epsm1(:,:,:)
 263  logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
 264  character(len=80) :: title(2)
 265  character(len=fnlen) :: gw_fname,wfk_fname,lwl_fname
 266  type(littlegroup_t),pointer :: Ltg_q(:)
 267  type(Paw_an_type),allocatable :: Paw_an(:)
 268  type(Paw_ij_type),allocatable :: Paw_ij(:)
 269  type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
 270  type(Pawrhoij_type),allocatable :: Pawrhoij(:),prev_Pawrhoij(:)
 271  type(pawpwff_t),allocatable :: Paw_pwff(:)
 272  type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
 273 
 274 !************************************************************************
 275 
 276  DBG_ENTER("COLL")
 277 
 278  call timab(301,1,tsec) ! overall time
 279  call timab(302,1,tsec) ! screening(init
 280 
 281  write(msg,'(6a)')&
 282 & ' SCREENING: Calculation of the susceptibility and dielectric matrices ',ch10,ch10,&
 283 & ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,&
 284 & ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.'
 285  call wrtout(ab_out, msg,'COLL')
 286  call wrtout(std_out,msg,'COLL')
 287 
 288  if(dtset%ucrpa>0) then
 289    write(msg,'(6a)')ch10,&
 290 &   ' cRPA Calculation: The calculation of the polarisability is constrained (ucrpa/=0)',ch10
 291    call wrtout(ab_out, msg,'COLL')
 292    call wrtout(std_out,msg,'COLL')
 293  end if
 294 #if defined HAVE_GW_DPC
 295  if (gwpc/=8) then
 296    write(msg,'(6a)')ch10,&
 297 &   ' Number of bytes for double precision complex /=8 ',ch10,&
 298 &   ' Cannot continue due to kind mismatch in BLAS library ',ch10,&
 299 &   ' Some BLAS interfaces are not generated by abilint '
 300    MSG_ERROR(msg)
 301  end if
 302  write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10
 303 #else
 304  write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10
 305 #endif
 306  call wrtout(std_out,msg,'COLL')
 307  call wrtout(ab_out, msg,'COLL')
 308 
 309  ! === Initialize MPI variables, and parallelization level ===
 310  ! gwpara: 0--> sequential run, 1--> parallelism over k-points, 2--> parallelism over bands.
 311  ! gwpara==2, each node has both fully and partially occupied states while conduction bands are divided
 312  comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs  = xmpi_comm_size(comm)
 313 
 314  if (my_rank == master) then
 315    wfk_fname = dtfil%fnamewffk
 316    if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then
 317      MSG_ERROR(msg)
 318    end if
 319  end if
 320  call xmpi_bcast(wfk_fname, master, comm, ierr)
 321 
 322  ! Some variables need to be initialized/nullify at start
 323  call energies_init(KS_energies)
 324  usexcnhat=0
 325 
 326  call mkrdim(acell,rprim,rprimd)
 327  call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
 328 
 329 !=== Define FFT grid(s) sizes ===
 330 ! Be careful! This mesh is only used for densities and potentials. It is NOT the (usually coarser)
 331 ! GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90.
 332 ! See also NOTES in the comments at the beginning of this file.
 333 ! NOTE: The mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn.
 334 
 335  k0(:)=zero
 336  call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,&
 337 & gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0)
 338 
 339  call print_ngfft(ngfftf,'Dense FFT mesh used for densities and potentials')
 340  nfftf_tot=PRODUCT(ngfftf(1:3))
 341 
 342  ! We can intialize MPI_enreg and fft distrib here, now ngfft are known
 343  call initmpi_seq(MPI_enreg_seq) ! Fake MPI_type for the sequential part.
 344  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all')
 345  call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all')
 346 
 347 !=============================================
 348 !==== Open and read pseudopotential files ====
 349 !=============================================
 350  call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm)
 351 
 352 !=== Initialize dimensions and basic objects ===
 353  call setup_screening(codvsn,acell,rprim,ngfftf,wfk_fname,dtfil,Dtset,Psps,Pawtab,&
 354 & ngfft_gw,Hdr_wfk,Hdr_local,Cryst,Kmesh,Qmesh,KS_BSt,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm)
 355 
 356  call timab(302,2,tsec) ! screening(init)
 357  call print_ngfft(ngfft_gw,'FFT mesh used for oscillator strengths')
 358 
 359  nfftgw_tot=PRODUCT(ngfft_gw(1:3))
 360  mgfftgw   =MAXVAL (ngfft_gw(1:3))
 361  nfftgw    =nfftgw_tot ! no FFT //
 362 
 363 !TRYING TO RECREATE AN "ABINIT ENVIRONMENT"
 364  KS_energies%e_corepsp=ecore/Cryst%ucvol
 365 
 366 !==========================
 367 !=== PAW initialization ===
 368 !==========================
 369  if (Dtset%usepaw==1) then
 370    call timab(315,1,tsec) ! screening(pawin
 371 
 372    call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,Cryst%xred)
 373 
 374    ABI_DT_MALLOC(Pawrhoij,(Cryst%natom))
 375    nspden_rhoij=pawrhoij_get_nspden(Dtset%nspden,Dtset%nspinor,Dtset%pawspnorb)
 376    call pawrhoij_alloc(Pawrhoij,Dtset%pawcpxocc,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,&
 377 &   Cryst%typat,pawtab=Pawtab)
 378 
 379    ! Initialize values for several basic arrays stored in Pawinit
 380    gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
 381 
 382    ! Test if we have to call pawinit
 383    call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
 384 
 385    if (psp_gencond==1.or.call_pawinit) then
 386      gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
 387      call pawinit(gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,&
 388 &     Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,&
 389 &     Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%xclevel,Dtset%usepotzero)
 390 
 391      ! Update internal values
 392      call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
 393    else
 394      if (Pawtab(1)%has_kij  ==1) Pawtab(1:Cryst%ntypat)%has_kij  =2
 395      if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2
 396    end if
 397    Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore)
 398 
 399 !  Initialize optional flags in Pawtab to zero
 400 !  (Cannot be done in Pawinit since the routine is called only if some pars. are changed)
 401    Pawtab(:)%has_nabla = 0
 402    Pawtab(:)%usepawu   = 0
 403    Pawtab(:)%useexexch = 0
 404    Pawtab(:)%exchmix   =zero
 405 !
 406 !  * Evaluate <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for the long wavelength limit.
 407 !  TODO solve problem with memory leak and clean this part as well as the associated flag
 408    call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab)
 409 
 410    call setsym_ylm(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,rprimd,Cryst%symrec,Pawang%zarot)
 411 
 412 !  * Initialize and compute data for LDA+U.
 413 !  paw_dmft%use_dmft=dtset%usedmft
 414    if (Dtset%usepawu>0.or.Dtset%useexexch>0) then
 415      call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,&
 416 &     Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,Cryst%ntypat,Pawang,Dtset%pawprtvol,&
 417 &     Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa)
 418    end if
 419 
 420    if (my_rank == master) call pawtab_print(Pawtab)
 421 
 422    ! Get Pawrhoij from the header of the WFK file.
 423    call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij)
 424 
 425    ! Re-symmetrize symrhoij.
 426 !  this call leads to a SIGFAULT, likely some pointer is not initialized correctly
 427    choice=1; optrhoij=1; ipert=0; idir=0
 428 !  call symrhoij(Pawrhoij,Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,&
 429 !  &             Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,Pawtab,&
 430 !  &             Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
 431 !
 432    ! Evaluate form factors for the radial part of phi.phj-tphi.tphj ===
 433    ! rhoxsp_method=1 ! Arnaud-Alouani
 434    ! rhoxsp_method=2 ! Shiskin-Kresse
 435    rhoxsp_method=2
 436 
 437    ! At least for ucrpa, the Arnaud Alouani is always a better choice but needs a larger cutoff
 438    if(dtset%ucrpa>0) rhoxsp_method=1
 439    if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc
 440 
 441    ABI_MALLOC(gw_gfft,(3,nfftgw_tot))
 442    call get_gftt(ngfft_gw,(/zero,zero,zero/),gmet,gw_gsq,gw_gfft)
 443    ABI_FREE(gw_gfft)
 444 
 445 !  Set up q grids, make qmax 20% larger than largest expected:
 446    ABI_MALLOC(nq_spl,(Psps%ntypat))
 447    ABI_MALLOC(qmax,(Psps%ntypat))
 448    nq_spl = Psps%mqgrid_ff
 449    qmax = SQRT(gw_gsq)*1.2d0  !qmax = Psps%qgrid_ff(Psps%mqgrid_ff)
 450    ABI_DT_MALLOC(Paw_pwff,(Psps%ntypat))
 451 
 452    call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps)
 453    ABI_FREE(nq_spl)
 454    ABI_FREE(qmax)
 455 
 456    ! Variables/arrays related to the fine FFT grid
 457    ABI_MALLOC(nhat,(nfftf,Dtset%nspden))
 458    nhat=zero; cplex=1
 459    ABI_DT_MALLOC(Pawfgrtab,(Cryst%natom))
 460    call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
 461    call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat)
 462    ABI_FREE(l_size_atm)
 463    compch_fft=greatest_real
 464    usexcnhat=MAXVAL(Pawtab(:)%usexcnhat)
 465 !  * 0 --> Vloc in atomic data is Vbare    (Blochl s formulation)
 466 !  * 1 --> Vloc in atomic data is VH(tnzc) (Kresse s formulation)
 467    write(msg,'(a,i3)')' screening : using usexcnhat = ',usexcnhat
 468    call wrtout(std_out,msg,'COLL')
 469 !
 470 !  Identify parts of the rectangular grid where the density has to be calculated.
 471    optcut=0;optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm
 472    if (Dtset%pawcross==1) optrad=1
 473    if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm
 474 
 475    call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,&
 476 &   optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
 477 
 478    call timab(315,2,tsec) ! screening(pawin
 479  else
 480    ! allocate empty structure for the sake of -fcheck=all...
 481    ABI_DT_MALLOC(Paw_pwff,(0))
 482    ABI_DT_MALLOC(Pawrhoij,(0))
 483    ABI_DT_MALLOC(Pawfgrtab,(0))
 484  end if ! End of PAW initialization.
 485 
 486  ! Consistency check and additional stuff done only for GW with PAW.
 487  ABI_DT_MALLOC(Paw_onsite,(Cryst%natom))
 488  if (Dtset%usepaw==1) then
 489    if (Dtset%ecutwfn < Dtset%ecut) then
 490      write(msg,"(5a)")&
 491 &     "WARNING - ",ch10,&
 492 &     "  It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,&
 493 &     "  an excessive truncation of the planewave basis set can lead to unphysical results."
 494      call wrtout(ab_out,msg,'COLL')
 495    end if
 496 
 497    ABI_CHECK(Dtset%useexexch==0,"LEXX not yet implemented in GW")
 498    ABI_CHECK(Dtset%usedmft==0,"DMFT + GW not available")
 499 
 500    if (Dtset%pawcross==1) then
 501      optgrad=1
 502      call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,Cryst%rprimd,&
 503 &     Cryst%xcart,Pawtab,Pawrad,Pawfgrtab,optgrad)
 504    end if
 505  end if
 506 
 507 !Allocate these arrays anyway, since they are passed to subroutines.
 508  if (.not.allocated(nhat))  then
 509    ABI_MALLOC(nhat,(nfftf,0))
 510  end if
 511 
 512  call timab(316,1,tsec) ! screening(wfs
 513 
 514 !=====================================================
 515 !=== Prepare the distribution of the wavefunctions ===
 516 !=====================================================
 517 ! valence and partially occupied are replicate on each node  while conduction bands are MPI distributed.
 518 ! This method is mandatory if gwpara==2 and/or we are using awtr==1 or the spectral method.
 519 ! If awtr==1, we evaluate chi0 taking advantage of time-reversal (speed-up~2)
 520 ! Useful indeces:
 521 !       nbvw = Max. number of fully/partially occupied states over spin
 522 !       nbcw = Max. number of unoccupied states considering the spin
 523 !TODO:
 524 !  Here for semiconducting systems we have to be sure that each processor has all the
 525 !  states considered in the SCGW, moreover nbsc<nbvw
 526 !  in case of SCGW vale and conduction has to be recalculated to avoid errors
 527 !  if a metal becomes semiconductor or viceversa.
 528 !  Ideally nbvw should include only the states v such that the transition
 529 !  c-->v is taken into account in cchi0 (see GW_TOLDOCC). In the present implementation
 530 
 531  ABI_MALLOC(ks_occ_idx,(KS_BSt%nkpt,KS_BSt%nsppol))
 532  ABI_MALLOC(ks_vbik   ,(KS_BSt%nkpt,KS_BSt%nsppol))
 533  ABI_MALLOC(qp_vbik   ,(KS_BSt%nkpt,KS_BSt%nsppol))
 534 
 535  call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=0)
 536  ks_occ_idx = get_occupied(KS_BSt,tol8) ! tol8 to be consistent when the density
 537  ks_vbik    = get_valence_idx(KS_BSt)
 538 
 539  ibocc(:)=MAXVAL(ks_occ_idx(:,:),DIM=1) ! Max occupied band index for each spin.
 540  ABI_FREE(ks_occ_idx)
 541 
 542  use_tr=.FALSE.; nbvw=0
 543  if (Dtset%gwpara==2.or.Ep%awtr==1.or.Dtset%spmeth>0) then
 544    use_tr=.TRUE.
 545    nbvw=MAXVAL(ibocc)
 546    nbcw=Ep%nbnds-nbvw
 547    write(msg,'(4a,i0,2a,i0,2a,i0,a)')ch10,&
 548 &   '- screening: taking advantage of time-reversal symmetry ',ch10,&
 549 &   '- Maximum band index for partially occupied states nbvw = ',nbvw,ch10,&
 550 &   '- Remaining bands to be divided among processors   nbcw = ',nbcw,ch10,&
 551 &   '- Number of bands treated by each node ~',nbcw/nprocs,ch10
 552    call wrtout(ab_out,msg,'COLL')
 553    if (Cryst%timrev/=2) then
 554      MSG_ERROR('Time-reversal cannot be used since cryst%timrev/=2')
 555    end if
 556  end if
 557 
 558  mband=Ep%nbnds
 559  ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol))
 560  nband=mband
 561  ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol))
 562  ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol))
 563  bks_mask=.FALSE.; keep_ur=.FALSE.
 564 
 565  ! autoparal section
 566  if (dtset%max_ncpus /=0) then
 567    ount = ab_out
 568    ! Temporary table needed to estimate memory
 569    ABI_MALLOC(nlmn_atm,(Cryst%natom))
 570    if (Dtset%usepaw==1) then
 571      do iat=1,Cryst%natom
 572        nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size
 573      end do
 574    end if
 575 
 576    write(ount,'(a)')"--- !Autoparal"
 577    write(ount,"(a)")"# Autoparal section for Screening runs"
 578    write(ount,"(a)")   "info:"
 579    write(ount,"(a,i0)")"    autoparal: ",dtset%autoparal
 580    write(ount,"(a,i0)")"    max_ncpus: ",dtset%max_ncpus
 581    write(ount,"(a,i0)")"    gwpara: ",dtset%gwpara
 582    write(ount,"(a,i0)")"    nkpt: ",dtset%nkpt
 583    write(ount,"(a,i0)")"    nsppol: ",dtset%nsppol
 584    write(ount,"(a,i0)")"    nspinor: ",dtset%nspinor
 585    write(ount,"(a,i0)")"    nbnds: ",Ep%nbnds
 586 
 587    work_size = nbvw * nbcw * Kmesh%nibz**2 * Dtset%nsppol
 588 
 589    ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI.
 590    nonscal_mem = (two*gwpc*Ep%npwe**2*(Ep%nomega*b2Mb)) * 1.1_dp
 591 
 592    ! List of configurations.
 593    ! Assuming an OpenMP implementation with perfect speedup!
 594    write(ount,"(a)")"configurations:"
 595 
 596    do ii=1,dtset%max_ncpus
 597      nt_per_proc = 0
 598      eff = HUGE(one)
 599      max_wfsmem_mb = zero
 600 
 601      do my_rank=0,ii-1
 602        call chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,ii,bks_mask,keep_ur,ierr)
 603        if (ierr /= 0) exit
 604        nt_per_proc = MAX(nt_per_proc, COUNT(bks_mask(1:nbvw,:,:)) * COUNT(bks_mask(nbvw+1:,:,:)))
 605        eff = MIN(eff, (one * work_size) / (ii * nt_per_proc))
 606 
 607        ! Memory needed for Fourier components ug.
 608        my_nbks = COUNT(bks_mask)
 609        ug_mem = two*gwpc*Dtset%nspinor*Ep%npwwfn*my_nbks*b2Mb
 610 
 611        ! Memory needed for real space ur.
 612        ur_mem = two*gwpc*Dtset%nspinor*nfftgw*COUNT(keep_ur)*b2Mb
 613 
 614        ! Memory needed for PAW projections Cprj
 615        cprj_mem = zero
 616        if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb
 617 
 618        max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem)
 619      end do
 620      if (ierr /= 0) cycle
 621 
 622      ! Add the non-scalable part and increase by 10% to account for other datastructures.
 623      mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp
 624 
 625      do omp_ncpus=1,xomp_get_max_threads()
 626        write(ount,"(a,i0)")"    - tot_ncpus: ",ii * omp_ncpus
 627        write(ount,"(a,i0)")"      mpi_ncpus: ",ii
 628        write(ount,"(a,i0)")"      omp_ncpus: ",omp_ncpus
 629        write(ount,"(a,f12.9)")"      efficiency: ",eff
 630        write(ount,"(a,f12.2)")"      mem_per_cpu: ",mempercpu_mb
 631      end do
 632    end do
 633    write(ount,'(a)')"..."
 634 
 635    ABI_FREE(nlmn_atm)
 636    MSG_ERROR_NODUMP("aborting now")
 637  else
 638    call chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr)
 639  end if
 640 
 641 ! Initialize the Wf_info object (allocate %ug and %ur if required).
 642  opt_ecut=Dtset%ecutwfn
 643 !opt_ecut=zero
 644 
 645 !call gsph_init(Gsph_wfn,Cryst,Ep%npwvec,gvec=gvec_kss)
 646 
 647  call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,Ep%npwwfn,mband,nband,Ep%nkibz,Dtset%nsppol,bks_mask,&
 648 & Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_gw,&
 649 & Gsph_wfn%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut)
 650 
 651  if (Dtset%pawcross==1) then
 652    call wfd_init(Wfdf,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,Ep%npwwfn,mband,nband,Ep%nkibz,Dtset%nsppol,bks_mask,&
 653 &   Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_gw,&
 654 &   Gsph_wfn%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut)
 655  end if
 656 
 657  ABI_FREE(bks_mask)
 658  ABI_FREE(nband)
 659  ABI_FREE(keep_ur)
 660 
 661  call wfd_print(Wfd,mode_paral='PERS')
 662 !FIXME: Rewrite the treatment of use_tr branches in cchi0 ...
 663 !Use a different nbvw for each spin.
 664 !Now use_tr means that one can use time-reversal symmetry.
 665 
 666 !==================================================
 667 !==== Read KS band structure from the KSS file ====
 668 !==================================================
 669  call wfd_read_wfk(Wfd,wfk_fname,iomode_from_fname(wfk_fname))
 670 
 671  if (Dtset%pawcross==1) then
 672    call wfd_copy(Wfd,Wfdf)
 673    call wfd_change_ngfft(Wfdf,Cryst,Psps,ngfftf)
 674  end if
 675 
 676  ! This test has been disabled (too expensive!)
 677  if (.False.) call wfd_test_ortho(Wfd,Cryst,Pawtab,unit=ab_out,mode_paral="COLL")
 678 
 679  call timab(316,2,tsec) ! screening(wfs
 680  call timab(319,1,tsec) ! screening(1)
 681 
 682  if (Cryst%nsym/=Dtset%nsym .and. Dtset%usepaw==1) then
 683    MSG_ERROR('Cryst%nsym/=Dtset%nsym, check pawinit and symrhoij')
 684  end if
 685 
 686 ! Get the FFT index of $ (R^{-1}(r-\tau)) $
 687 ! S= $\transpose R^{-1}$ and k_BZ = S k_IBZ
 688 ! irottb is the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Sk.
 689  ABI_MALLOC(irottb,(nfftgw,Cryst%nsym))
 690  call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfft_gw,irottb,iscompatibleFFT)
 691 
 692  ABI_MALLOC(ktabr,(nfftgw,Kmesh%nbz))
 693  do ikbz=1,Kmesh%nbz
 694    isym=Kmesh%tabo(ikbz)
 695    do ifft=1,nfftgw
 696      ktabr(ifft,ikbz)=irottb(ifft,isym)
 697    end do
 698  end do
 699  ABI_FREE(irottb)
 700 
 701  if (Dtset%usepaw==1 .and. Dtset%pawcross==1) then
 702    ABI_MALLOC(irottb,(nfftf,Cryst%nsym))
 703    call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfftf,irottb,iscompatibleFFT)
 704 
 705    ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz))
 706    do ikbz=1,Kmesh%nbz
 707      isym=Kmesh%tabo(ikbz)
 708      do ifft=1,nfftf
 709        ktabrf(ifft,ikbz)=irottb(ifft,isym)
 710      end do
 711    end do
 712    ABI_FREE(irottb)
 713  else
 714    ABI_MALLOC(ktabrf,(0,0))
 715  end if
 716 
 717 !=== Compute structure factor phases and large sphere cut-off ===
 718 !WARNING cannot use Dtset%mgfft, this has to be checked better
 719 !mgfft=MAXVAL(ngfftc(:))
 720 !allocate(ph1d(2,3*(2*mgfft+1)*Cryst%natom),ph1df(2,3*(2*mgfftf+1)*Cryst%natom))
 721  write(std_out,*)' CHECK ',Dtset%mgfftdg,mgfftf
 722  !if (Dtset%mgfftdg/=mgfftf) write(std_out,*)"WARNING Dtset%mgfftf /= mgfftf"
 723  ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom))
 724  ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom))
 725  call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred)
 726 
 727  if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then
 728    call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred)
 729  else
 730    ph1df(:,:)=ph1d(:,:)
 731  end if
 732 
 733 ! Initialize QP_BSt using KS bands
 734 ! In case of SCGW, update QP_BSt using the QPS file.
 735  call ebands_copy(KS_BSt,QP_BSt)
 736 
 737  call timab(319,2,tsec) ! screening(1)
 738 
 739 !============================
 740 !==== Self-consistent GW ====
 741 !============================
 742  if (Ep%gwcalctyp>=10) then
 743    call timab(304,1,tsec) ! KS => QP; [wfrg]
 744 
 745    ! Initialize with KS eigenvalues and eigenfunctions.
 746    ABI_MALLOC(m_lda_to_qp,(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol))
 747    m_lda_to_qp = czero
 748    do spin=1,Wfd%nsppol
 749      do ik_ibz=1,Wfd%nkibz
 750        do band=1,Wfd%nband(ik_ibz,spin)
 751          m_lda_to_qp(band,band,:,:) = cone
 752        end do
 753      end do
 754    end do
 755 
 756    ! Read unitary transformation and QP energies.
 757    ! TODO switch on the renormalization of n in screening, QPS should report bdgw
 758    ABI_MALLOC(rhor_p,(nfftf,Dtset%nspden))
 759    ABI_DT_MALLOC(prev_Pawrhoij,(Cryst%natom*Psps%usepaw))
 760 
 761    call rdqps(QP_BSt,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,&
 762 &   nfftf,ngfftf,Cryst%ucvol,Dtset%paral_kgb,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_lda_to_qp,rhor_p,prev_Pawrhoij)
 763 
 764    ABI_FREE(rhor_p)
 765    ABI_DT_FREE(prev_Pawrhoij)
 766 
 767    ! FIXME this is to preserve the old implementation for the head and the wings in ccchi0q0
 768    ! But has to be rationalized
 769    KS_BSt%eig=QP_BSt%eig
 770 
 771    ! Calculate new occ. factors and fermi level.
 772    call ebands_update_occ(QP_BSt,Dtset%spinmagntarget)
 773    qp_vbik(:,:) = get_valence_idx(QP_BSt)
 774 
 775    ! === Update only the wfg treated with GW ===
 776    ! For PAW update and re-symmetrize cprj in the full BZ, TODO add rotation in spinor space
 777    if (nscf/=0) call wfd_rotate(Wfd,Cryst,m_lda_to_qp)
 778 
 779    ABI_FREE(m_lda_to_qp)
 780    call timab(304,2,tsec)
 781  end if ! gwcalctyp>=10
 782 
 783  call timab(305,1,tsec) ! screening(densit
 784 !
 785 !=== In case update the eigenvalues ===
 786 !* Either use a scissor operator or an external GW file.
 787  gw_fname = "__in.gw__"
 788  update_energies = file_exists(gw_fname)
 789 
 790  if (ABS(Ep%mbpt_sciss)>tol6) then
 791    write(msg,'(5a,f7.3,a)')&
 792 &   ' screening : performing a first self-consistency',ch10,&
 793 &   ' update of the energies in W by a scissor operator',ch10,&
 794 &   ' applying a scissor operator of [eV] : ',Ep%mbpt_sciss*Ha_eV,ch10
 795    call wrtout(std_out,msg,'COLL')
 796    call wrtout(ab_out,msg,'COLL')
 797    call apply_scissor(QP_BSt,Ep%mbpt_sciss)
 798  else if (update_energies) then
 799    write(msg,'(4a)')&
 800 &   ' screening : performing a first self-consistency',ch10,&
 801 &   ' update of the energies in W by a previous GW calculation via GW file: ',TRIM(gw_fname)
 802    call wrtout(std_out,msg,'COLL')
 803    call wrtout(ab_out,msg,'COLL')
 804    ABI_MALLOC(igwene,(QP_Bst%mband,QP_Bst%nkpt,QP_Bst%nsppol))
 805    call rdgw(QP_BSt,gw_fname,igwene,extrapolate=.FALSE.)
 806 !  call rdgw(QP_BSt,gw_fname,igwene,extrapolate=.TRUE.)
 807    ABI_FREE(igwene)
 808    call ebands_update_occ(QP_BSt,Dtset%spinmagntarget)
 809  end if
 810 
 811 !========================
 812 !=== COMPUTE DENSITY ====
 813 !========================
 814 !* Evaluate PW part (complete charge in case of NC pseudos)
 815 !TODO this part has to be rewritten. If I decrease the tol on the occupations
 816 !I have to code some MPI stuff also if use_tr==.TRUE.
 817 
 818  ABI_MALLOC(rhor,(nfftf,Dtset%nspden))
 819  ABI_MALLOC(taur,(nfftf,Dtset%nspden*Dtset%usekden))
 820 
 821  call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,rhor)
 822  if (Dtset%usekden==1) then
 823    call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,taur,optcalc=1)
 824  end if
 825 
 826  call timab(305,2,tsec) ! screening(densit
 827 
 828  nhatgrdim = 0
 829  if (Dtset%usepaw==1) then ! Additional computation for PAW.
 830    call timab(320,1,tsec) ! screening(paw
 831 
 832 !  Add the compensation charge to the PW density.
 833    nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc
 834    cplex=1; ider=2*nhatgrdim; izero=0
 835    if (nhatgrdim>0)  then
 836      ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,3))
 837    else
 838      ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0))
 839    end if
 840    call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,Cryst%gprimd,&
 841 &   Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
 842 &   Pawfgrtab,nhatgr,nhat,Pawrhoij,Pawrhoij,Pawtab,k0,Cryst%rprimd,Cryst%ucvol,dtset%usewvl,Cryst%xred)
 843 
 844 !  === Evaluate onsite energies, potentials, densities ===
 845 !  * Initialize variables/arrays related to the PAW spheres.
 846 !  * Initialize also lmselect (index of non-zero LM-moments of densities).
 847    cplex=1
 848    ABI_DT_MALLOC(Paw_ij,(Cryst%natom))
 849    call paw_ij_nullify(Paw_ij)
 850    call paw_ij_init(Paw_ij,cplex,Dtset%nspinor,Wfd%nsppol,&
 851 &   Wfd%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
 852 &   has_dij=1,has_dijhartree=1,has_exexch_pot=1,has_pawu_occ=1)
 853 
 854    nkxc1=0
 855    ABI_DT_MALLOC(Paw_an,(Cryst%natom))
 856    call paw_an_nullify(Paw_an)
 857    call paw_an_init(Paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
 858 &   cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0)
 859 
 860    nzlmopt=-1; option=0; compch_sph=greatest_real
 861    call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert,Dtset%ixc,&
 862 &   Cryst%natom,Cryst%natom,Dtset%nspden,Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,Paw_an,Paw_an,&
 863 &   Paw_ij,Pawang,Dtset%pawprtvol,Pawrad,Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,&
 864 &   Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp)
 865    call timab(320,2,tsec) ! screening(paw
 866  else
 867    ABI_DT_MALLOC(Paw_ij,(0))
 868    ABI_DT_MALLOC(Paw_an,(0))
 869  end if ! usepaw
 870 
 871  call timab(321,1,tsec) ! screening(2)
 872 
 873  !JB : Should be remove : cf. l 839
 874  if (.not.allocated(nhatgr))  then
 875    ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0))
 876  end if
 877 
 878  call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,rhor,ucvol,&
 879 & Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,omegaplasma)
 880 
 881 !For PAW, add the compensation charge the FFT mesh, then get rho(G).
 882  if (Dtset%usepaw==1) rhor(:,:)=rhor(:,:)+nhat(:,:)
 883 
 884  call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,rhor,ucvol=ucvol)
 885  if(Dtset%usekden==1)then
 886    call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,taur,ucvol=ucvol,optrhor=1)
 887  end if
 888 
 889  if (dtset%gwgamma>0) then
 890    ABI_MALLOC(rhor_kernel,(nfftf,Dtset%nspden))
 891  end if
 892 
 893  ABI_MALLOC(rhog,(2,nfftf))
 894  ABI_MALLOC(taug,(2,nfftf*Dtset%usekden))
 895  call fourdp(1,rhog,rhor(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp4)
 896  if(Dtset%usekden==1)then
 897    call fourdp(1,taug,taur(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp4)
 898  end if
 899 
 900 !The following steps have been gathered in the setvtr routine:
 901 !- get Ewald energy and Ewald forces
 902 !- compute local ionic pseudopotential vpsp
 903 !- eventually compute 3D core electron density xccc3d
 904 !- eventually compute vxc and vhartr
 905 !- set up ks_vtrial
 906 !**************************************************************
 907 !**** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
 908 !**************************************************************
 909 
 910  ngrvdw=0
 911  ABI_MALLOC(grvdw,(3,ngrvdw))
 912  ABI_MALLOC(grchempottn,(3,Cryst%natom))
 913  ABI_MALLOC(grewtn,(3,Cryst%natom))
 914  nkxc=0
 915  if (Dtset%nspden==1) nkxc=2
 916  if (Dtset%nspden>=2) nkxc=3 ! check GGA and spinor that is messy !!!
 917  ! If MGGA, fxc and kxc are not available and we dont need them for the screening part (for now ...)
 918  if (Dtset%ixc<0 .and. libxc_functionals_ismgga()) nkxc=0
 919  if (nkxc/=0)  then
 920    ABI_MALLOC(kxc,(nfftf,nkxc))
 921  end if
 922 
 923  n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf
 924  ABI_MALLOC(xccc3d,(n3xccc))
 925  ABI_MALLOC(ks_vhartr,(nfftf))
 926  ABI_MALLOC(ks_vtrial,(nfftf,Dtset%nspden))
 927  ABI_MALLOC(vpsp,(nfftf))
 928  ABI_MALLOC(ks_vxc,(nfftf,Dtset%nspden))
 929 
 930  optene=4; moved_atm_inside=0; moved_rhor=0; initialized=1; istep=1
 931  call setvtr(Cryst%atindx1,Dtset,KS_energies,Cryst%gmet,Cryst%gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,istep,kxc,mgfftf,&
 932 & moved_atm_inside,moved_rhor,MPI_enreg_seq, &
 933 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,Cryst%ntypat,&
 934 & Psps%n1xccc,n3xccc,optene,pawrad,Pawtab,ph1df,Psps,rhog,rhor,Cryst%rmet,Cryst%rprimd,strsxc,Cryst%ucvol,usexcnhat,&
 935 & ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,wvl_dummy,xccc3d,Cryst%xred,taug=taug,taur=taur)
 936 
 937  if (nkxc/=0)  then
 938    ABI_FREE(kxc)
 939  end if
 940  ABI_FREE(grchempottn)
 941  ABI_FREE(grewtn)
 942  ABI_FREE(grvdw)
 943  ABI_FREE(xccc3d)
 944 
 945 !============================
 946 !==== Compute KS PAW Dij ====
 947 !============================
 948  if (Dtset%usepaw==1) then
 949    call timab(561,1,tsec)
 950 
 951 !  Calculate unsymmetrized Dij.
 952    cplex=1; ipert=0; idir=0
 953    call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert,&
 954 &   Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
 955 &   Dtset%nspden,Cryst%ntypat,Paw_an,Paw_ij,Pawang,Pawfgrtab,Dtset%pawprtvol,&
 956 &   Pawrad,Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,k0,Dtset%spnorbscl,&
 957 &   Cryst%ucvol,dtset%charge,ks_vtrial,ks_vxc,Cryst%xred,&
 958 &   nucdipmom=Dtset%nucdipmom)
 959 
 960 !  Symmetrize KS Dij
 961 #if 0
 962    call symdij(Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,Cryst%natom,&
 963 &   Cryst%nsym,Cryst%ntypat,0,Paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,&
 964 &   Cryst%symrec)
 965 #else
 966    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,Cryst%natom,&
 967 &   Cryst%nsym,Cryst%ntypat,Paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,&
 968 &   Cryst%symrec)
 969 #endif
 970 !
 971 !  Output of the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
 972    call pawprt(Dtset,Cryst%natom,Paw_ij,Pawrhoij,Pawtab)
 973    call timab(561,2,tsec)
 974  end if
 975 !
 976 !=== Calculate the frequency mesh ===
 977 !* First omega is always zero without broadening.
 978 !FIXME what about metals? I think we should add eta, this means we need to know if the system is metallic, for example using occopt
 979 !MS Modified to account for non-zero starting frequency (19-11-2010)
 980 !MS Modified for tangent grid (07-01-2011)
 981  ABI_MALLOC(Ep%omega,(Ep%nomega))
 982  Ep%omega(1)=CMPLX(Ep%omegaermin,zero,kind=dpc)
 983 
 984  if (Ep%nomegaer>1) then ! Avoid division by zero.
 985    if (Dtset%gw_frqre_tangrid==0.and.Dtset%gw_frqre_inzgrid==0) then
 986      domegareal=(Ep%omegaermax-Ep%omegaermin)/(Ep%nomegaer-1)
 987      do iomega=2,Ep%nomegaer
 988        Ep%omega(iomega)=CMPLX(Ep%omegaermin+(iomega-1)*domegareal,zero,kind=dpc)
 989      end do
 990    else if (Dtset%gw_frqre_tangrid==1.and.Dtset%gw_frqre_inzgrid==0) then ! We have tangent transformed grid
 991      MSG_WARNING('EXPERIMENTAL - Using tangent transform grid for contour deformation.')
 992      Ep%omegaermax = Dtset%cd_max_freq
 993      Ep%omegaermin = zero
 994      ifirst=1; ilast=Ep%nomegaer
 995      if (Dtset%cd_subset_freq(1)/=0) then ! Only a subset of frequencies is being calculated
 996        ifirst=Dtset%cd_subset_freq(1); ilast=Dtset%cd_subset_freq(2)
 997      end if
 998      factor = Dtset%cd_halfway_freq/TAN(pi*quarter)
 999 !    *Important*-here nfreqre is used because the step is set by the original grid
1000      domegareal=(ATAN(Ep%omegaermax/factor)*two*piinv)/(Dtset%nfreqre-1) ! Stepsize in transformed variable
1001      do iomega=1,Ep%nomegaer
1002        Ep%omega(iomega)=CMPLX(factor*TAN((iomega+ifirst-2)*domegareal*pi*half),zero,kind=dpc)
1003      end do
1004      Ep%omegaermin = REAL(Ep%omega(1))
1005      Ep%omegaermax = REAL(Ep%omega(Ep%nomegaer))
1006    else if (Dtset%gw_frqre_tangrid==0.and.Dtset%gw_frqre_inzgrid==1) then
1007      e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma
1008      domegareal=one/(Ep%nomegaer)
1009      do iomega=1,Ep%nomegaer
1010        factor = (iomega-1)*domegareal
1011        Ep%omega(iomega)=CMPLX(e0*factor/(one-factor),zero,kind=dpc)
1012      end do
1013      Ep%omegaermin = REAL(Ep%omega(1))
1014      Ep%omegaermax = REAL(Ep%omega(Ep%nomegaer))
1015    else
1016      MSG_ERROR('Error in specification of real frequency grid')
1017    end if
1018  end if
1019 
1020  if (Ep%plasmon_pole_model.and.Ep%nomega==2) then
1021    e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma
1022    Ep%omega(2)=CMPLX(zero,e0,kind=dpc)
1023  end if
1024 !
1025 !=== For AC, use Gauss-Legendre quadrature method ===
1026 !* Replace $ \int_0^\infty dx f(x) $ with $ \int_0^1 dz f(1/z - 1)/z^2 $.
1027 !* Note that the grid is not log as required by CD, so we cannot use the same SCR file.
1028  if (Ep%analytic_continuation) then
1029    ABI_MALLOC(z,(Ep%nomegaei))
1030    ABI_MALLOC(zw,(Ep%nomegaei))
1031    call coeffs_gausslegint(zero,one,z,zw,Ep%nomegaei)
1032    do iomega=1,Ep%nomegaei
1033      Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,one/z(iomega)-one,kind=dpc)
1034    end do
1035    ABI_FREE(z)
1036    ABI_FREE(zw)
1037  else if (Ep%contour_deformation.and.(Dtset%cd_customnimfrqs/=0)) then
1038    Ep%omega(Ep%nomegaer+1)=CMPLX(zero,Dtset%cd_imfrqs(1))
1039    do iomega=2,Ep%nomegaei
1040      if (Dtset%cd_imfrqs(iomega)<=Dtset%cd_imfrqs(iomega-1)) then
1041        MSG_ERROR(' Specified imaginary frequencies need to be strictly increasing!')
1042      end if
1043      Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,Dtset%cd_imfrqs(iomega))
1044    end do
1045  else if (Ep%contour_deformation.and.(Dtset%gw_frqim_inzgrid/=0)) then
1046    e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma
1047    domegareal=one/(Ep%nomegaei+1)
1048    do iomega=1,Ep%nomegaei
1049      factor = iomega*domegareal
1050      Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,e0*factor/(one-factor),kind=dpc)
1051    end do
1052  else if (Ep%contour_deformation.and.(Ep%nomegaei/=0)) then
1053    e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma
1054    do iomega=1,Ep%nomegaei
1055      Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,e0/(Dtset%freqim_alpha-two)&
1056 &     *(EXP(two/(Ep%nomegaei+1)*LOG(Dtset%freqim_alpha-one)*iomega)-one),kind=dpc)
1057    end do
1058  end if
1059 
1060  if (Dtset%cd_full_grid/=0) then ! Full grid will be calculated
1061 !  Grid values are added after the last imaginary freq.
1062    do ios=1,Ep%nomegaei
1063      do iomega=2,Ep%nomegaer
1064        Ep%omega(Ep%nomegaer+Ep%nomegaei+(ios-1)*(Ep%nomegaer-1)+(iomega-1)) = &
1065 &       CMPLX(REAL(Ep%omega(iomega)),AIMAG(Ep%omega(Ep%nomegaer+ios)))
1066      end do
1067    end do
1068  end if
1069 
1070  ! Report frequency mesh for chi0.
1071  write(msg,'(2a)')ch10,' calculating chi0 at frequencies [eV] :'
1072  call wrtout(std_out,msg,'COLL'); call wrtout(ab_out,msg,'COLL')
1073  do iomega=1,Ep%nomega
1074    write(msg,'(i3,2es16.6)')iomega,Ep%omega(iomega)*Ha_eV
1075    call wrtout(std_out,msg,'COLL'); call wrtout(ab_out,msg,'COLL')
1076  end do
1077 
1078  ! Allocate chi0, wings and array for chi0_sumrule check.
1079  ABI_MALLOC(chi0_sumrule,(Ep%npwe))
1080 
1081  write(msg,'(a,f12.1,a)')' Memory required for chi0 matrix= ',two*gwpc*Ep%npwe**2*Ep%nI*Ep%nJ*Ep%nomega*b2Mb," [Mb]."
1082  call wrtout(std_out,msg,'COLL')
1083  ABI_STAT_MALLOC(chi0,(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega), ierr)
1084  ABI_CHECK(ierr==0, "Out of memory in chi0")
1085 !
1086 !============================== END OF THE INITIALIZATION PART ===========================
1087 !
1088 !======================================================================
1089 !==== Loop over q-points. Calculate \epsilon^{-1} and save on disc ====
1090 !======================================================================
1091  call timab(321,2,tsec) ! screening(2)
1092 
1093  iqcalc = 0
1094  do iqibz=1,Qmesh%nibz
1095    call timab(306,1,tsec)
1096    is_first_qcalc=(iqibz==1)
1097 
1098    ! Selective q-point calculation.
1099    found=.FALSE.; label=iqibz
1100    if (Ep%nqcalc/=Ep%nqibz) then
1101      do ii=1,Ep%nqcalc
1102        qtmp(:)=Qmesh%ibz(:,iqibz)-Ep%qcalc(:,ii)
1103        found=(normv(qtmp,gmet,'G')<GW_TOLQ)
1104        if (found) then
1105          label=ii; EXIT !ii
1106        end if
1107      end do
1108      if (.not.found) CYCLE !iqibz
1109      qtmp(:)=Ep%qcalc(:,1)-Qmesh%ibz(:,iqibz)
1110      is_first_qcalc=(normv(qtmp,gmet,'G')<GW_TOLQ)
1111    end if
1112    iqcalc = iqcalc + 1
1113 
1114    bar=REPEAT('-',80)
1115    write(msg,'(4a,1x,a,i2,a,f9.6,2(",",f9.6),3a)')ch10,ch10,bar,ch10,&
1116 &   ' q-point number ',label,'        q = (',(Qmesh%ibz(ii,iqibz),ii=1,3),') [r.l.u.]',ch10,bar
1117    call wrtout(std_out,msg,'COLL')
1118    call wrtout(ab_out,msg,'COLL')
1119    is_qeq0 = 0; if (normv(Qmesh%ibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1120 
1121    call timab(306,2,tsec)
1122 
1123    if (is_qeq0==1) then
1124      ! Special treatment of the long wavelength limit.
1125      call timab(307,1,tsec)
1126 
1127      ABI_MALLOC(chi0_lwing,(Ep%npwe*Ep%nI,Ep%nomega,3))
1128      ABI_MALLOC(chi0_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,3))
1129      ABI_MALLOC(chi0_head,(3,3,Ep%nomega))
1130 
1131      call cchi0q0(use_tr,Dtset,Cryst,Ep,Psps,Kmesh,QP_BSt,KS_BSt,Gsph_epsG0,&
1132 &     Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,Pawfgrtab,Paw_onsite,ktabr,ktabrf,nbvw,ngfft_gw,nfftgw,&
1133 &     ngfftf,nfftf_tot,chi0,chi0_head,chi0_lwing,chi0_uwing,Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf)
1134 
1135      chihw = chi_new(ep%npwe, ep%nomega)
1136      chihw%head = chi0_head
1137      chihw%lwing = chi0_lwing
1138      chihw%uwing = chi0_uwing
1139 
1140      ! Add the intraband term if required and metallic occupation scheme is used.
1141      add_chi0_intraband=.FALSE. !add_chi0_intraband=.TRUE.
1142      if (add_chi0_intraband .and. ebands_has_metal_scheme(QP_BSt)) then
1143 
1144        ABI_STAT_MALLOC(chi0intra,(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega), ierr)
1145        ABI_CHECK(ierr==0, "Out of memory in chi0intra")
1146 
1147        ABI_MALLOC(chi0intra_lwing,(Ep%npwe*Ep%nI,Ep%nomega,3))
1148        ABI_MALLOC(chi0intra_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,3))
1149        ABI_MALLOC(chi0intra_head,(3,3,Ep%nomega))
1150 
1151        call chi0q0_intraband(Wfd,Cryst,Ep,Psps,QP_BSt,Gsph_epsG0,Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,use_tr,Dtset%usepawu,&
1152        ngfft_gw,chi0intra,chi0intra_head,chi0intra_lwing,chi0intra_uwing)
1153 
1154        call wrtout(std_out,"Head of chi0 and chi0_intra","COLL")
1155        do iomega=1,Ep%nomega
1156          write(std_out,*)Ep%omega(iomega)*Ha_eV,REAL(chi0(1,1,iomega)),REAL(chi0intra(1,1,iomega))
1157          write(std_out,*)Ep%omega(iomega)*Ha_eV,AIMAG(chi0(1,1,iomega)),AIMAG(chi0intra(1,1,iomega))
1158        end do
1159 
1160        chi0       = chi0       + chi0intra
1161        chi0_head  = chi0_head  + chi0intra_head
1162        chi0_lwing = chi0_lwing + chi0intra_lwing
1163        chi0_uwing = chi0_uwing + chi0intra_uwing
1164 
1165        ABI_FREE(chi0intra)
1166        ABI_FREE(chi0intra_lwing)
1167        ABI_FREE(chi0intra_uwing)
1168        ABI_FREE(chi0intra_head)
1169      end if
1170 
1171      if (.False.) then
1172        lwl_fname = strcat(dtfil%filnam_ds(4), "_LWL")
1173        call lwl_write(lwl_fname,cryst,vcp,ep%npwe,ep%nomega,gsph_epsg0%gvec,chi0,chi0_head,chi0_lwing,chi0_uwing,comm)
1174      end if
1175 
1176      call timab(307,2,tsec)
1177 
1178    else
1179      ! Calculate cchi0 for q/=0.
1180      call timab(308,1,tsec)
1181 
1182      call cchi0(use_tr,Dtset,Cryst,Qmesh%ibz(:,iqibz),Ep,Psps,Kmesh,QP_BSt,Gsph_epsG0,&
1183 &     Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,nbvw,ngfft_gw,nfftgw,ngfftf,nfftf_tot,chi0,ktabr,ktabrf,&
1184 &     Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf)
1185 
1186      call timab(308,2,tsec)
1187    end if
1188 
1189    ! Print chi0(q,G,Gp,omega), then calculate epsilon and epsilon^-1 for this q-point.
1190    ! Only master works but this part could be parallelized over frequencies.
1191    call timab(309,1,tsec)
1192 
1193    do iomega=1,MIN(Ep%nomega,NOMEGA_PRINTED)
1194      write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]'
1195      call wrtout(std_out,msg,'COLL')
1196      call wrtout(ab_out,msg,'COLL')
1197      write(msg,'(1x,a,i3,a,i4,a)')' chi0(q =',iqibz, ', omega =',iomega,', G,G'')'
1198      if (Ep%nqcalc/=Ep%nqibz) write(msg,'(a,i3,a,i4,a)')'  chi0(q=',iqcalc,', omega=',iomega,', G,G'')'
1199      call wrtout(std_out,msg,'COLL')
1200 !    arr99 is needed to avoid the update of all the tests. Now chi0 is divided by ucvol inside (cchi0|cchi0q0).
1201 !    TODO should be removed but GW tests have to be updated.
1202      ii = MIN(9,Ep%npwe)
1203      ABI_MALLOC(arr_99,(ii,ii))
1204      arr_99 = chi0(1:ii,1:ii,iomega)*ucvol
1205      call print_arr(arr_99,max_r=2,unit=ab_out)
1206      call print_arr(arr_99,unit=std_out)
1207      ABI_FREE(arr_99)
1208      !call print_arr(chi0(:,:,iomega),max_r=2,unit=ab_out)
1209      !call print_arr(chi0(:,:,iomega),unit=std_out)
1210    end do
1211 
1212    if (Ep%nomega>NOMEGA_PRINTED) then
1213      write(msg,'(a,i3,a)')' No. of calculated frequencies > ',NOMEGA_PRINTED,', stop printing '
1214      call wrtout(std_out,msg,'COLL')
1215      call wrtout(ab_out,msg,'COLL')
1216    end if
1217 
1218    ! Write chi0 to _SUSC file
1219    ! Master creates and write the header if this is the first q-point calculated.
1220    if (Dtset%prtsuscep>0 .and. my_rank==master) then
1221      title(1)='CHI0 file: chi0'
1222      title(2)=' '
1223      if (is_qeq0==1) then
1224        string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string)
1225        title(1)=title(1)(1:21)//', calculated using inclvkb = '//string
1226      end if
1227 
1228      ! Open file and write header for polarizability files.
1229      if (is_first_qcalc) then
1230        ikxc=0; test_type=0; tordering=1
1231        hchi0 = hscr_new("polarizability",dtset,ep,hdr_local,ikxc,test_type,tordering,title,Ep%npwe,Gsph_epsG0%gvec)
1232        if (dtset%iomode == IO_MODE_ETSF) then
1233 #ifdef HAVE_NETCDF
1234          NCF_CHECK(nctk_open_create(unt_susc, nctk_ncify(dtfil%fnameabo_sus), xmpi_comm_self))
1235          NCF_CHECK(crystal_ncwrite(cryst, unt_susc))
1236          NCF_CHECK(ebands_ncwrite(QP_BSt, unt_susc))
1237 #endif
1238        else
1239          unt_susc=Dtfil%unchi0
1240          if (open_file(dtfil%fnameabo_sus,msg,unit=unt_susc,status='unknown',form='unformatted') /= 0) then
1241            MSG_ERROR(msg)
1242          end if
1243        end if
1244        fform_chi0 = hchi0%fform
1245        call hscr_io(hchi0,fform_chi0,2,unt_susc,xmpi_comm_self,0,Dtset%iomode)
1246        call hscr_free(Hchi0)
1247      end if
1248      call write_screening("polarizability",unt_susc,Dtset%iomode,Ep%npwe,Ep%nomega,iqcalc,chi0)
1249    end if
1250 
1251 !  Calculate the RPA functional correlation energy if the polarizability on a
1252 !  Gauss-Legendre mesh along imaginary axis is available
1253    if (Ep%analytic_continuation .and. Dtset%gwrpacorr>0 ) then
1254      if (is_first_qcalc) then
1255        ABI_MALLOC(ec_rpa,(Dtset%gwrpacorr))
1256        ec_rpa(:)=zero
1257      end if
1258      call calc_rpa_functional(Dtset%gwrpacorr,label,iqibz,Ep,Vcp,Qmesh,Dtfil,gmet,chi0,comm,ec_rpa)
1259      if (label==Ep%nqcalc) then
1260        ABI_FREE(ec_rpa)
1261      end if
1262    end if
1263 
1264 !  ==========================================================
1265 !  === Calculate RPA \tilde\epsilon^{-1} overwriting chi0 ===
1266 !  ==========================================================
1267    approx_type=0 ! RPA
1268    option_test=0 ! TESTPARTICLE
1269    dim_wing=0; if (is_qeq0==1) dim_wing=3
1270 
1271    if (dim_wing==0) then
1272      dim_wing=1
1273      if (.not.allocated(chi0_lwing))  then
1274        ABI_MALLOC(chi0_lwing,(Ep%npwe*Ep%nI,Ep%nomega,dim_wing))
1275      end if
1276      if (.not.allocated(chi0_uwing))  then
1277        ABI_MALLOC(chi0_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,dim_wing))
1278      end if
1279      if (.not.allocated(chi0_head ))  then
1280        ABI_MALLOC(chi0_head,(dim_wing,dim_wing,Ep%nomega))
1281      end if
1282      dim_wing=0
1283    end if
1284 
1285 #if 0
1286 !  Using the random q for the optical limit is one of the reasons
1287 !  why sigma breaks the initial energy degeneracies.
1288    chi0_lwing=czero
1289    chi0_uwing=czero
1290    chi0_head=czero
1291 #endif
1292 
1293    ! Setup flags for the computation of em1
1294    ! If the vertex is being included for the spectrum, calculate the kernel now and pass it on
1295    if (dtset%gwgamma>0) rhor_kernel = rhor
1296 
1297    select case (dtset%gwgamma)
1298    case (0)
1299      approx_type=0; option_test=0; dim_kxcg=0
1300      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1301 
1302    case (1, 2)
1303      ! ALDA TDDFT kernel vertex
1304      ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available")
1305      MSG_WARNING('EXPERIMENTAL: Kernel is being added to screening, the SCR file will be non-standard!!')
1306      ikxc=7; approx_type=1; dim_kxcg=1
1307      if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1308      if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1309      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1310      call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,rhor_kernel,&
1311      Ep%npwe,dim_kxcg,kxcg,Gsph_epsG0%gvec,xmpi_comm_self)
1312 
1313    case (3, 4)
1314      ! ADA non-local kernel vertex
1315      ABI_CHECK(Wfd%usepaw==0,"ADA vertex + PAW not available")
1316      ABI_CHECK(Wfd%nsppol==1,"ADA vertex for GWGamma not available yet for spin-polarised cases")
1317      MSG_WARNING('EXPERIMENTAL: Kernel is being added to screening, the SCR file will be non-standard!!')
1318      ikxc=7; approx_type=2
1319      if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1320      if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1321      ABI_MALLOC(fxc_ADA,(Ep%npwe,Ep%npwe,Ep%nqibz))
1322 !    Use userrd to set kappa
1323      if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp
1324 !    Set correct value of kappa (should be scaled with alpha*r_s where)
1325 !    r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3)
1326      rhoav = (omegaplasma*omegaplasma)/four_pi
1327      r_s = (three/(four_pi*rhoav))**third
1328      alpha = (four*ninth*piinv)**third
1329      Dtset%userrd = Dtset%userrd*alpha*r_s
1330 
1331      call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf,Wfd%nspden,rhor_kernel,Ep%npwe,Ep%nqibz,Ep%qibz,&
1332      fxc_ADA,Gsph_epsG0%gvec,xmpi_comm_self,kappa_init=Dtset%userrd)
1333 
1334      dim_kxcg=0
1335      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1336 
1337 !  @WC: bootstrap --
1338    case (-3, -4, -5, -6, -7, -8)
1339      ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available")
1340      if (Dtset%gwgamma>-5) then
1341        MSG_WARNING('EXPERIMENTAL: Bootstrap kernel is being added to screening')
1342        approx_type=4
1343      else if (Dtset%gwgamma>-7) then
1344        MSG_WARNING('EXPERIMENTAL: Bootstrap kernel (head-only) is being added to screening')
1345        approx_type=5
1346      else
1347        MSG_WARNING('EXPERIMENTAL: Bootstrap kernel (RPA-type, head-only) is being added to screening')
1348        approx_type=6
1349      end if
1350      dim_kxcg=0
1351      option_test=MOD(Dtset%gwgamma,2)
1352      ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma
1353      ! 0 -> TESTPARTICLE, vertex in chi0 only
1354      ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1355 !--@WC
1356 
1357    case default
1358      MSG_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma)))
1359    end select
1360 
1361    if (approx_type<2) then
1362      call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,&
1363      approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,&
1364      chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm)
1365    else if (approx_type<3) then !@WC: ADA
1366      call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,&
1367      approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,&
1368      chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm,fxc_ADA=fxc_ADA(:,:,iqibz))
1369    else if (approx_type<7) then !@WC: bootstrap
1370      call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,&
1371      approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,&
1372      chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm)
1373    else
1374      MSG_ERROR(sjoin("Wrong approx_type:", itoa(approx_type)))
1375    end if
1376 
1377    ABI_FREE(chi0_lwing)
1378    ABI_FREE(chi0_uwing)
1379    ABI_FREE(chi0_head)
1380 
1381    if (my_rank==master .and. is_qeq0==1) then
1382      call spectra_repr(spectra,msg)
1383      call wrtout(std_out,msg,'COLL')
1384      call wrtout(ab_out,msg,'COLL')
1385 
1386      if (Ep%nomegaer>2) then
1387        call spectra_write(spectra,W_EELF  ,Dtfil%fnameabo_eelf)
1388        call spectra_write(spectra,W_EM_LF ,Dtfil%fnameabo_em1_lf)
1389        call spectra_write(spectra,W_EM_NLF,Dtfil%fnameabo_em1_nlf)
1390      end if
1391    end if ! master and is_qeq0==1
1392 
1393    if (is_qeq0==1) call chi_free(chihw)
1394 
1395    call spectra_free(spectra)
1396    if (allocated(kxcg)) then
1397      ABI_FREE(kxcg)
1398    end if
1399    if (allocated(fxc_ADA)) then
1400      ABI_FREE(fxc_ADA)
1401    end if
1402    !
1403    ! Output the sum rule evaluation.
1404    ! Vcp%vc_sqrt(:,iqibz) Contains vc^{1/2}(q,G), complex-valued due to a possible cutoff
1405    epsm1 => chi0
1406    call output_chi0sumrule((is_qeq0==1),iqibz,Ep%npwe,omegaplasma,chi0_sumrule,epsm1(:,:,1),Vcp%vc_sqrt(:,iqibz))
1407 
1408    ! If input variable npvel is larger than 0, trigger the Random Stopping Power calculation
1409    ! Only the masternode is used
1410    if (my_rank==master .and. Dtset%npvel>0) then
1411      if (is_first_qcalc) then
1412        ABI_MALLOC(rspower,(Dtset%npvel))
1413        rspower(:)=zero
1414      end if
1415      call random_stopping_power(iqibz,Dtset%npvel,Dtset%pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower)
1416      if (label==Ep%nqcalc) then
1417        ABI_FREE(rspower)
1418      end if
1419    end if
1420 
1421    ! Write heads and wings to main output file.
1422    if (is_qeq0==1) then
1423      write(msg,'(1x,2a)')' Heads and wings of the symmetrical epsilon^-1(G,G'') ',ch10
1424      call wrtout(ab_out,msg,'COLL')
1425      do iomega=1,Ep%nomega
1426        write(msg,'(2x,a,i4,a,2f9.4,a)')' Upper and lower wings at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]'
1427        call wrtout(ab_out,msg,'COLL')
1428        call print_arr(epsm1(1,:,iomega),max_r=9,unit=ab_out)
1429        call print_arr(epsm1(:,1,iomega),max_r=9,unit=ab_out)
1430        call wrtout(ab_out,ch10,'COLL')
1431      end do
1432    end if
1433 
1434    call timab(309,2,tsec)
1435    call timab(310,1,tsec) ! wrscr
1436 
1437    if (my_rank==master) then
1438      ! === Write the symmetrical epsilon^-1 on file ===
1439      title(1)='SCR file: epsilon^-1'
1440      if (is_qeq0==1) then
1441        string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string)
1442        title(1)=title(1)(1:21)//', calculated using inclvkb = '//string
1443      end if
1444      title(2)='TESTPARTICLE'
1445      ctype='RPA'
1446      title(2)(14:17)=ctype !this has to be modified
1447 
1448      if (is_first_qcalc) then
1449        ! === Open file and write the header for the SCR file ===
1450        ! * Here we write the RPA approximation for \tilde\epsilon^{-1}
1451        ikxc=0; test_type=0; tordering=1
1452        hem1 = hscr_new("inverse_dielectric_function",dtset,ep,hdr_local,ikxc,test_type,tordering,title,&
1453 &       Ep%npwe,Gsph_epsG0%gvec)
1454        fform_em1 = hem1%fform
1455        if (dtset%iomode == IO_MODE_ETSF) then
1456 #ifdef HAVE_NETCDF
1457          NCF_CHECK(nctk_open_create(unt_em1, nctk_ncify(dtfil%fnameabo_scr), xmpi_comm_self))
1458          NCF_CHECK(crystal_ncwrite(cryst, unt_em1))
1459          NCF_CHECK(ebands_ncwrite(QP_BSt, unt_em1))
1460 #endif
1461        else
1462          unt_em1=Dtfil%unscr
1463          if (open_file(dtfil%fnameabo_scr,msg,unit=unt_em1,status='unknown',form='unformatted') /= 0) then
1464            MSG_ERROR(msg)
1465          end if
1466        end if
1467        call hscr_io(hem1,fform_em1,2,unt_em1,xmpi_comm_self,0,Dtset%iomode)
1468        call hscr_free(Hem1)
1469      end if
1470 
1471      call write_screening("inverse_dielectric_function",unt_em1,Dtset%iomode,Ep%npwe,Ep%nomega,iqcalc,epsm1)
1472    end if ! my_rank==master
1473 
1474    call timab(310,2,tsec)
1475  end do ! Loop over q-points
1476 
1477  !----------------------------- END OF THE CALCULATION ------------------------
1478 
1479  ! Close Files.
1480  if (my_rank==master) then
1481    if (dtset%iomode == IO_MODE_ETSF) then
1482 #ifdef HAVE_NETCDF
1483      NCF_CHECK(nf90_close(unt_em1))
1484      if (dtset%prtsuscep>0) then
1485        NCF_CHECK(nf90_close(unt_susc))
1486      end if
1487 #endif
1488    else
1489      close(unt_em1)
1490      if (dtset%prtsuscep>0) close(unt_susc)
1491    end if
1492  end if
1493 !
1494 !=====================
1495 !==== Free memory ====
1496 !=====================
1497  ABI_FREE(chi0_sumrule)
1498  ABI_FREE(chi0)
1499  if (allocated(rhor_kernel)) then
1500    ABI_FREE(rhor_kernel)
1501  end if
1502 
1503  ABI_FREE(rhor)
1504  ABI_FREE(rhog)
1505  ABI_FREE(ks_vbik)
1506  ABI_FREE(qp_vbik)
1507  ABI_FREE(ktabr)
1508  ABI_FREE(taur)
1509  ABI_FREE(taug)
1510  ABI_FREE(ks_vhartr)
1511  ABI_FREE(ks_vtrial)
1512  ABI_FREE(vpsp)
1513  ABI_FREE(ks_vxc)
1514  ABI_FREE(ph1d)
1515  ABI_FREE(ph1df)
1516  ABI_FREE(nhatgr)
1517  ABI_FREE(nhat)
1518  call pawfgr_destroy(Pawfgr)
1519 
1520  if (Dtset%usepaw==1) then ! Optional deallocation for PAW.
1521    call pawrhoij_free(Pawrhoij)
1522    call pawfgrtab_free(Pawfgrtab)
1523    call paw_ij_free(Paw_ij)
1524    call paw_an_free(Paw_an)
1525    call pawpwff_free(Paw_pwff)
1526    if (Dtset%pawcross==1) then
1527      call paw_pwaves_lmn_free(Paw_onsite)
1528      Wfdf%bks_comm = xmpi_comm_null
1529      call wfd_free(Wfdf)
1530    end if
1531  end if
1532 
1533  ABI_DT_FREE(Pawfgrtab)
1534  ABI_DT_FREE(Paw_pwff)
1535  ABI_DT_FREE(Pawrhoij)
1536  ABI_DT_FREE(Paw_ij)
1537  ABI_DT_FREE(Paw_an)
1538  ABI_FREE(ktabrf)
1539  ABI_DT_FREE(Paw_onsite)
1540 
1541  call wfd_free(Wfd)
1542  call kmesh_free(Kmesh)
1543  call kmesh_free(Qmesh)
1544  call crystal_free(Cryst)
1545  call gsph_free(Gsph_epsG0)
1546  call gsph_free(Gsph_wfn)
1547  call vcoul_free(Vcp)
1548  call em1params_free(Ep)
1549  call hdr_free(Hdr_wfk)
1550  call hdr_free(Hdr_local)
1551  call ebands_free(KS_BSt)
1552  call ebands_free(QP_BSt)
1553  call littlegroup_free(Ltg_q)
1554  call destroy_mpi_enreg(MPI_enreg_seq)
1555  ABI_DT_FREE(Ltg_q)
1556 
1557  call timab(301,2,tsec)
1558 
1559  DBG_EXIT("COLL")
1560 
1561 end subroutine screening

m_screening_driver/setup_screening [ Functions ]

[ Top ] [ m_screening_driver ] [ Functions ]

NAME

 setup_screening

FUNCTION

  Initialize the Ep% data type containing the parameters used during the screening calculation.
  as well as basic objects describing the BZ sampling .... TODO list to be completed

INPUTS

 wfk_fname=Name of the input WFK file.
 acell(3)=length scales of primitive translations (Bohr).
 rprim(3,3)=dimensionless real space primitive translations.
 ngfftf(18)=Contain all needed information about the 3D FFT for densities and potentials.
 dtfil <type(datafiles_type)>=variables related to files

OUTPUT

 ngfft_gw(18)=Contain all needed information about the 3D FFT for the oscillator strengths.
  See ~abinit/doc/variables/vargs.htm#ngfft
 Ltg_q(:)<littlegroup_t>,=
 Ep<em1params_t>=Parameters for the screening calculation.
  Most part of it is Initialized and checked.
 Hdr_wfk type(Hdr_type)=Header of the KSS file.
 Cryst<crystal_t>=Definition of the unit cell and its symmetries.
 Kmesh<kmesh_t>=Structure defining the k-point sampling (wavefunctions).
 Qmesh<kmesh_t>=Structure defining the q-point sampling (screening)
 Gsph_wfn<gsphere_t>=Structure defining the G-sphere for the wavefunctions (not k-dependent).
 Gsph_epsG0<gsphere_t>=The G-sphere for the screening, enlarged to take into account for umklapps.
 Psps <Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the KSS file
 Vcp <type vcoul_t> datatype gathering information on the coulombian cutoff technique
 comm=MPI communicator.

SIDE EFFECTS

 Dtset<Dataset_type>=All input variables for this dataset.
  %ecutwfn, %npwwfn,
  %ecuteps, %npweps
   might be redefined in setshells in order to close the shell.

PARENTS

      screening

CHILDREN

      xmpi_split_work2_i4b

SOURCE

1609 subroutine setup_screening(codvsn,acell,rprim,ngfftf,wfk_fname,dtfil,Dtset,Psps,Pawtab,&
1610 & ngfft_gw,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,KS_BSt,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm)
1611 
1612 
1613 !This section has been created automatically by the script Abilint (TD).
1614 !Do not modify the following lines by hand.
1615 #undef ABI_FUNC
1616 #define ABI_FUNC 'setup_screening'
1617 !End of the abilint section
1618 
1619  implicit none
1620 
1621 !Arguments ------------------------------------
1622 !scalars
1623  integer,intent(in) :: comm
1624  character(len=6),intent(in) :: codvsn
1625  character(len=fnlen),intent(in) :: wfk_fname
1626  type(Dataset_type),intent(inout) :: Dtset !INOUT is due to setshells
1627  type(datafiles_type),intent(in) :: dtfil
1628  type(Pseudopotential_type),intent(in) :: Psps
1629  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw)
1630  type(em1params_t),intent(out) :: Ep
1631  type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out
1632  type(ebands_t),intent(out) :: KS_BSt
1633  type(kmesh_t),intent(out) :: Kmesh,Qmesh
1634  type(crystal_t),intent(out) :: Cryst
1635  type(gsphere_t),intent(out) :: Gsph_epsG0,Gsph_wfn
1636  type(vcoul_t),intent(out) :: Vcp
1637 !arrays
1638  integer,intent(in) :: ngfftf(18)
1639  integer,intent(out) :: ngfft_gw(18)
1640  real(dp),intent(in) :: acell(3),rprim(3,3)
1641  type(littlegroup_t),pointer :: Ltg_q(:)
1642 
1643 !Local variables-------------------------------
1644 !scalars
1645  integer,parameter :: NOMEGAGAUSS=30,NOMEGAREAL=201,pertcase0=0,master=0
1646  integer :: bantot,ib,ibtot,ikibz,iq,iqp,isppol,ig,ng,ierr
1647  integer :: jj,mod10,mband,ng_kss,iqbz,isym,iq_ibz,itim
1648  integer :: timrev,use_umklp,ncerr
1649  integer :: npwepG0,nshepspG0,method,enforce_sym,nfftgw_tot !,spin,band,ik_ibz,
1650  integer :: istart,iend,test_npwkss,my_rank,nprocs !ii
1651  real(dp),parameter :: OMEGAERMAX=100.0/Ha_eV
1652  real(dp) :: ecutepspG0,ucvol,domegareal
1653  logical :: remove_inv,ltest,found,is_static,has_q0
1654  character(len=500) :: msg
1655  type(wvl_internal_type) :: wvl
1656 !arrays
1657  integer :: ng0sh_opt(3)
1658  integer,allocatable :: npwarr(:)
1659  integer,pointer :: gvec_kss(:,:)
1660  integer,pointer :: test_gvec_kss(:,:)
1661  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),qtmp(3),sq(3),qbz(3)
1662  real(dp),pointer :: energies_p(:,:,:)
1663  real(dp),allocatable :: doccde(:),eigen(:),occfact(:)
1664  type(Pawrhoij_type),allocatable :: Pawrhoij(:)
1665 
1666 ! *************************************************************************
1667 
1668  DBG_ENTER('COLL')
1669 
1670  ! Check for calculations that are not implemented
1671  ltest = ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol) == Dtset%nband(1))
1672  ABI_CHECK(ltest, 'dtset%nband(:) must be constant in the GW code.')
1673 
1674  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1675 
1676  call mkrdim(acell,rprim,rprimd)
1677  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1678 
1679  ! Set up basic parameters of the calculation
1680  Ep%gwcalctyp            =Dtset%gwcalctyp
1681  Ep%plasmon_pole_model   =.TRUE.
1682  Ep%analytic_continuation=.FALSE.
1683  Ep%contour_deformation  =.FALSE.
1684 
1685  mod10=MOD(Ep%gwcalctyp,10)
1686  if (mod10/=0.and.mod10/=8) Ep%plasmon_pole_model   =.FALSE.
1687  if (mod10==1)              Ep%analytic_continuation=.TRUE.
1688  if (mod10==2.or.mod10==9)  Ep%contour_deformation  =.TRUE.
1689  is_static=(mod10==5.or.mod10==6.or.mod10==7)
1690 
1691  Ep%nbnds  =Dtset%nband(1)
1692  Ep%symchi =Dtset%symchi
1693  Ep%inclvkb=Dtset%inclvkb; if (Dtset%usepaw/=0) Ep%inclvkb=0
1694  Ep%zcut   =Dtset%zcut
1695 
1696  write(msg,'(2a,i4,2a,f10.6,a)')ch10,&
1697 &  ' GW calculation type              = ',Ep%gwcalctyp,ch10,&
1698 &  ' zcut to avoid poles in chi0 [eV] = ',Ep%zcut*Ha_eV,ch10
1699  call wrtout(std_out,msg,'COLL')
1700 
1701  Ep%awtr  =Dtset%awtr
1702  Ep%npwe  =Dtset%npweps
1703  Ep%npwwfn=Dtset%npwwfn
1704  Ep%npwvec=MAX(Ep%npwe,Ep%npwwfn)
1705 
1706  timrev = 2 ! This information is not reported in the header
1707             ! 1 --> do not use time-reversal symmetry
1708             ! 2 --> take advantage of time-reversal symmetry
1709  if (any(dtset%kptopt == [3, 4])) timrev = 1
1710 
1711  if (timrev==1.and.Dtset%awtr/=0) then
1712    MSG_ERROR("awtr/=0 cannot be used when time-reversal symmetry doesn't hold")
1713  end if
1714 
1715  ! Read parameters from WFK and verifify them.
1716  call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm)
1717  mband = MAXVAL(Hdr_wfk%nband)
1718  call hdr_vs_dtset(Hdr_wfk,Dtset)
1719  remove_inv=.FALSE.
1720 
1721  test_npwkss = 0
1722  call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,&
1723 &   gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr)
1724  ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss")
1725 
1726  ABI_MALLOC(gvec_kss,(3,test_npwkss))
1727  gvec_kss = test_gvec_kss
1728  ng_kss = test_npwkss
1729 
1730  if (Ep%npwvec>ng_kss) then
1731    Ep%npwvec=ng_kss
1732    if (Ep%npwwfn> ng_kss) Ep%npwwfn=ng_kss
1733    if (Ep%npwe  > ng_kss) Ep%npwe  =ng_kss
1734    write(msg,'(3a,3(a,i6,a))')ch10,&
1735 &    ' Number of G-vectors found less then required. Calculation will proceed with ',ch10,&
1736 &    '  npwvec = ',Ep%npwvec,ch10,&
1737 &    '  npweps = ',Ep%npwe  ,ch10,&
1738 &    '  npwwfn = ',Ep%npwwfn,ch10
1739    MSG_WARNING(msg)
1740  end if
1741 
1742  ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2))
1743  ierr = 0
1744  do ig=1,ng
1745    if (ANY(gvec_kss(:,ig)/=test_gvec_kss(:,ig))) then
1746      ierr=ierr+1
1747      write(std_out,*)" gvec_kss ",ig,"/",ng,gvec_kss(:,ig),test_gvec_kss(:,ig)
1748    end if
1749  end do
1750  ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss")
1751  ABI_FREE(test_gvec_kss)
1752 
1753  ! Get important dimension from Hdr_wfk
1754  ! Check also the consistency btw Hdr_wfk and Dtset.
1755  Ep%nsppol=Hdr_wfk%nsppol
1756  Ep%nkibz =Hdr_wfk%nkpt
1757 
1758  if (Ep%nbnds>mband) then
1759    Ep%nbnds=mband
1760    Dtset%nband(:)=mband
1761    write(msg,'(4a,i4,a)')ch10,&
1762 &    ' Number of bands found less then required. ',ch10,&
1763 &    ' Calculation will proceed with nbnds = ',mband,ch10
1764    MSG_WARNING(msg)
1765  end if
1766 
1767  call crystal_from_hdr(Cryst,Hdr_wfk,timrev,remove_inv)
1768  call crystal_print(Cryst,mode_paral='COLL')
1769 
1770  ! === Create basic data types for the calculation ===
1771  ! * Kmesh defines the k-point sampling for the wavefunctions.
1772  ! * Qmesh defines the q-point sampling for chi0, all possible differences k1-k2 reduced to the IBZ.
1773  ! TODO Kmesh%bz should be [-half,half[ but this modification will be painful!
1774 
1775  call kmesh_init(Kmesh,Cryst,Ep%nkibz,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.FALSE.)
1776  !call kmesh_init(Kmesh,Cryst,Ep%nkibz,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.TRUE.)
1777  ! Some required information are not filled up inside kmesh_init
1778  ! So doing it here, even though it is not clean
1779  Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:)
1780  Kmesh%nshift        =Dtset%nshiftk
1781  ABI_ALLOCATE(Kmesh%shift,(3,Kmesh%nshift))
1782  Kmesh%shift(:,:)    =Dtset%shiftk(:,1:Dtset%nshiftk)
1783 
1784  call kmesh_print(Kmesh,"K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL")
1785  call kmesh_print(Kmesh,"K-mesh for the wavefunctions",ab_out, 0,           "COLL")
1786 
1787  ! === Find Q-mesh, and do setup for long wavelength limit ===
1788  ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed,
1789  !   epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed)
1790  call find_qmesh(Qmesh,Cryst,Kmesh)
1791 
1792  call kmesh_print(Qmesh,"Q-mesh for the screening function",std_out,Dtset%prtvol,"COLL")
1793  call kmesh_print(Qmesh,"Q-mesh for the screening function",ab_out ,0           ,"COLL")
1794 
1795  do iqbz=1,Qmesh%nbz
1796    call get_BZ_item(Qmesh,iqbz,qbz,iq_ibz,isym,itim)
1797    sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz))
1798    if (ANY(ABS(qbz-sq )>1.0d-4)) then
1799      write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')&
1800 &      ' qpoint ',qbz,' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,&
1801 &      ' through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,&
1802 &      ' however a non zero umklapp G_o vector is required and this is not yet allowed'
1803      MSG_ERROR(msg)
1804    end if
1805  end do
1806 
1807  ! Write the list of qpoints for the screening in netcdf format and exit.
1808  ! This file is used by abipy to generate multiple input files.
1809  if (Dtset%nqptdm == -1) then
1810    if (my_rank==master) then
1811 #ifdef HAVE_NETCDF
1812       ncerr = nctk_write_ibz(strcat(dtfil%filnam_ds(4), "_qptdms.nc"), qmesh%ibz, qmesh%wt)
1813       NCF_CHECK(ncerr)
1814 #endif
1815    end if
1816    MSG_ERROR_NODUMP("Aborting now")
1817  end if
1818 
1819  if (Dtset%gw_nqlwl==0) then
1820    Ep%nqlwl=1
1821    ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl))
1822    Ep%qlwl(:,1)=GW_Q0_DEFAULT ! Use default shift 0.000010, 0.000020, 0.000030
1823  else
1824    Ep%nqlwl=Dtset%gw_nqlwl
1825    ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl))
1826    Ep%qlwl(:,:)=Dtset%gw_qlwl(:,1:Ep%nqlwl)
1827    ABI_CHECK(Ep%nqlwl==1,"nqlwl/=1 not coded")
1828  end if
1829  !write(std_out,*)" Using qlwl = ",Ep%qlwl
1830 
1831  ! Find optimal value for G-sphere enlargment due to oscillator matrix elements
1832  call get_ng0sh(Kmesh%nbz,Kmesh%bz,Qmesh%nibz,Qmesh%ibz,Kmesh%nbz,Kmesh%bz,GW_TOLQ0,ng0sh_opt)
1833  call wrtout(std_out,sjoin(' Optimal value for ng0sh:',ltoa(ng0sh_opt)),"COLL")
1834 
1835  Ep%mG0(:)=ng0sh_opt(:) !Ep%mG0(:) = [3, 3, 3]
1836 
1837  ! === In case of symmetrization, find the little group of the q"s ===
1838  ! * For the long-wavelength limit we consider a small but finite q. However the oscillators are
1839  !  evaluated setting q==0. Thus it is possible to take advantage of symmetries also when q --> 0.
1840  ! * Here we calculate the enlargement of the G-sphere, npwepG0, needed to account for umklapps.
1841  ! TODO Switch on use_umklp, write all this stuff to ab_out
1842 
1843  Ep%npwepG0=Ep%npwe
1844  ABI_DT_MALLOC(Ltg_q,(Qmesh%nibz))
1845 
1846  do iq=1,Qmesh%nibz
1847    qtmp(:)=Qmesh%ibz(:,iq); if (normv(qtmp,gmet,'G')<GW_TOLQ0) qtmp(:)=zero; use_umklp=0
1848    call littlegroup_init(qtmp,Kmesh,Cryst,use_umklp,Ltg_q(iq),Ep%npwe,gvec=gvec_kss)
1849  end do
1850 
1851  ecutepspG0 = Dtset%ecuteps
1852  ABI_CHECK(ecutepspG0 > zero, "ecuteps must be > 0")
1853  if (Ep%symchi/=0) then
1854    ecutepspG0=MAXVAL(Ltg_q(:)%max_kin_gmG0)+tol6; npwepG0=0; nshepspG0=0
1855    write(std_out,*)" Due to umklapp processes : ecutepspg0= ",ecutepspG0
1856    call setshells(ecutepspG0,npwepG0,nshepspG0,Cryst%nsym,gmet,gprimd,Cryst%symrel,'eps_pG0',Cryst%ucvol)
1857    Ep%npwepG0=npwepG0
1858  end if
1859 
1860  if (Ep%npwepG0>Ep%npwvec) then
1861    write(msg,'(3a,i5,a,i5)')&
1862 &    ' npwepG0 > npwvec, decrease npweps or increase npwwfn. ',ch10,&
1863 &    ' npwepG0 = ',Ep%npwepG0,' npwvec = ',Ep%npwvec
1864    MSG_ERROR(msg)
1865  end if
1866 
1867  ! === Create structure describing the G-sphere used for chi0/espilon and Wfns ===
1868  ! * The cutoff is >= ecuteps to allow for umklapp
1869 #if 0
1870  call gsph_init(Gsph_wfn,Cryst,0,ecut=Dtset%ecutwfn)
1871  Dtset%npwwfn = Gsph_wfn%ng
1872  Ep%npwwfn = Gsph_wfn%ng
1873  ierr = 0
1874  do ig=1,MIN(Gsph_wfn%ng, ng_kss)
1875    if ( ANY(Gsph_wfn%gvec(:,ig) /= gvec_kss(:,ig)) ) then
1876      write(std_out,*)ig, Gsph_wfn%gvec(:,ig), gvec_kss(:,ig)
1877    end if
1878  end do
1879  ABI_CHECK(ierr==0,"Wrong gvec_wfn")
1880 #else
1881  call gsph_init(Gsph_wfn,Cryst,Ep%npwvec,gvec=gvec_kss)
1882 #endif
1883 
1884 #if 0
1885  call gsph_init(Gsph_epsG0,Cryst,0,ecut=ecutepspG0)
1886  Ep%npwepG0 = Gsph_epsG0%ng
1887  ierr = 0
1888  do ig=1,MIN(Gsph_epsG0%ng, ng_kss)
1889    if ( ANY(Gsph_epsG0%gvec(:,ig) /= gvec_kss(:,ig)) ) then
1890      write(std_out,*)ig, Gsph_epsG0%gvec(:,ig), gvec_kss(:,ig)
1891    end if
1892  end do
1893  ABI_CHECK(ierr==0,"Wrong gvec_epsG0")
1894 #else
1895  call gsph_init(Gsph_epsG0,Cryst,Ep%npwepG0,gvec=gvec_kss)
1896 #endif
1897  !
1898  ! =======================================================================
1899  ! ==== Setup of the FFT mesh used for the oscillator matrix elements ====
1900  ! =======================================================================
1901  ! * ngfft_gw(7:18) is the same as Dtset%ngfft(7:18), initialized before entering setup_screening.
1902  !   Here we just redefine ngfft_gw(1:6) according to the following options:
1903  !
1904  !    method==0 ==> FFT grid read from __fft.in__ (only for debugging purpose)
1905  !    method==1 ==> normal FFT grid
1906  !    method==2 ==> slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
1907  !    method==3 ==> doubled FFT grid, to treat exactly the convolution defining the density,
1908  !      Useful in sigma if ppmodel=[2,3,4] since rho(G-Gp) or to calculate matrix elements of v_Hxc.
1909  !
1910  !    enforce_sym==1 ==> enforce a direct space FFT mesh compatible with all symmetries operation
1911  !    enforce_sym==0 ==> Find the smallest FFT grid compatibile with the library, do not care about symmetries
1912  !
1913  ngfft_gw(1:18)=Dtset%ngfft(1:18); method=2
1914  if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0
1915  if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1
1916  if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2
1917  if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3
1918  enforce_sym=MOD(Dtset%fftgw,10)
1919 
1920  ! Use npwepG0 to account for umklapps.
1921  call setmesh(gmet,gvec_kss,ngfft_gw,Ep%npwvec,Ep%npwepG0,Ep%npwwfn,nfftgw_tot,method,Ep%mG0,Cryst,enforce_sym)
1922  !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Ep%mG0,enforce_sym,ngfft_gw,nfftgw_tot)
1923 
1924  ABI_FREE(gvec_kss)
1925 
1926  ! FIXME this wont work if nqptdm/=0
1927  call vcoul_init(Vcp,Gsph_epsG0,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%icutcoul,Dtset%vcutgeo,Dtset%ecuteps,Ep%npwe,Ep%nqlwl,&
1928 &  Ep%qlwl,ngfftf,comm)
1929 
1930 #if 0
1931  ! Using the random q for the optical limit is one of the reasons
1932  ! why sigma breaks the initial energy degeneracies.
1933  Vcp%i_sz=zero
1934  Vcp%vc_sqrt(1,:)=czero
1935  Vcp%vcqlwl_sqrt(1,:)=czero
1936 #endif
1937 
1938  ! Value of scissor energy
1939  Ep%mbpt_sciss=Dtset%mbpt_sciss
1940 
1941  ! Define the frequency mesh for epsilon according to the method used.
1942  Ep%nomegaei=1
1943  Ep%nomegaer=1; if (is_static) Ep%nomegaer=0
1944  Ep%nomegaec=0
1945  Ep%omegaermax=zero
1946 
1947  ! For ppmodels 2,3,4, only omega=0 is needed.
1948  if (Ep%plasmon_pole_model.and.Dtset%nfreqre==1.and.Dtset%nfreqim==0) then
1949    Ep%nomegaer=1; Ep%nomegaei=0
1950    write(msg,'(7a)')ch10,&
1951 &    ' The inverse dielectric matrix will be calculated on zero frequency only',ch10,&
1952 &    ' please note that the calculated epsilon^-1 cannot be used ',ch10,&
1953 &    ' to calculate QP corrections using plasmonpole model 1',ch10
1954    call wrtout(std_out,msg,'COLL')
1955    call wrtout(ab_out,msg,'COLL')
1956  end if
1957 
1958  ! Max number of omega along the imaginary axis
1959  if (Ep%analytic_continuation.or.Ep%contour_deformation) then
1960    Ep%nomegaei=Dtset%nfreqim
1961    if (Dtset%gw_frqim_inzgrid==1) then
1962      MSG_WARNING('iomega = z/1-z transfom grid will be used for imaginary frequency grid')
1963    end if
1964    if (Dtset%cd_customnimfrqs/=0) then
1965      MSG_WARNING('Custom imaginary grid specified. Assuming experienced user.')
1966      Ep%nomegaei=Dtset%cd_customnimfrqs
1967    end if
1968    if (Ep%nomegaei==-1) then
1969      Ep%nomegaei=NOMEGAGAUSS
1970      MSG_WARNING(sjoin('Number of imaginary frequencies set to default= ',itoa(NOMEGAGAUSS)))
1971    end if
1972    if (Ep%nomegaei==0) then
1973      MSG_WARNING(' nfreqim = 0! Assuming experienced user merging several frequency calculations.')
1974    end if
1975  end if
1976 
1977  ! Range and total number of real frequencies.
1978  Ep%omegaermin = zero
1979  if (Ep%contour_deformation) then
1980    Ep%nomegaer=Dtset%nfreqre; Ep%omegaermin=Dtset%freqremin; Ep%omegaermax=Dtset%freqremax
1981    if (Dtset%gw_frqre_tangrid==1) then
1982      Ep%omegaermax=Dtset%cd_max_freq
1983      MSG_WARNING('Tangent transfom grid will be used for real frequency grid')
1984    end if
1985    if (Dtset%gw_frqre_tangrid==1) then
1986      MSG_WARNING('Tangent transfom grid will be used for real frequency grid')
1987    end if
1988    if (Ep%nomegaer==-1) then
1989      Ep%nomegaer=NOMEGAREAL
1990      MSG_WARNING(sjoin('Number of real frequencies set to default= ',itoa(NOMEGAREAL)))
1991    end if
1992    if (Ep%nomegaer==0) then
1993      MSG_WARNING('nfreqre = 0 ! Assuming experienced user.')
1994    end if
1995    if (ABS(Ep%omegaermin)<TOL16) then
1996      Ep%omegaermin=zero
1997      write(msg,'(a,f8.4)')' Min real frequency set to default [Ha] = ',Ep%omegaermin
1998      MSG_WARNING(msg)
1999    end if
2000    if (Ep%omegaermin>Ep%omegaermax) then
2001      MSG_ERROR('freqremin > freqremax !')
2002    end if
2003    if (Ep%omegaermax<TOL16) then
2004      Ep%omegaermax=OMEGAERMAX
2005      write(msg,'(a,f8.4)')' Max real frequency set to default [Ha] = ',OMEGAERMAX
2006      MSG_WARNING(msg)
2007    end if
2008    ! Check if a subset of the frequencies is to be used
2009    if (Dtset%cd_subset_freq(1)/=0) then
2010      istart = Dtset%cd_subset_freq(1)
2011      iend   = Dtset%cd_subset_freq(2)
2012      if (istart>iend.or.istart<0.or.iend<0) then
2013        MSG_ERROR(' check indices of cd_subset_freq!')
2014      end if
2015      write(msg,'(2(a,i0))')' Using cd_subset_freq to only do freq. from ',istart,' to ',iend
2016      MSG_WARNING(msg)
2017      ! Reset the numbers
2018      if (Dtset%gw_frqre_tangrid/=1) then ! Normal equidistant grid
2019        Ep%nomegaer = iend-istart+1
2020        domegareal=(Ep%omegaermax-Ep%omegaermin)/(Ep%nomegaer-1)
2021        Ep%omegaermin = Ep%omegaermin+(istart-1)*domegareal
2022        Ep%omegaermax = Ep%omegaermin+(iend-1)*domegareal
2023      else
2024        Ep%nomegaer = iend-istart+1
2025      end if
2026    end if
2027  end if
2028 
2029  ! Check full grid calculations
2030  if (Dtset%cd_full_grid/=0) then
2031    MSG_WARNING("FULL GRID IN COMPLEX PLANE CALCULATED. YOU MIGHT NOT BE ABLE TO USE SCREENING FILES!")
2032    if (Dtset%cd_subset_freq(1)/=0) then
2033      MSG_ERROR('cd_subset_freq cannot be used with cd_full_grid!')
2034    end if
2035    Ep%nomegaec = Ep%nomegaei*(Ep%nomegaer-1)
2036  end if
2037 
2038  Ep%nomega=Ep%nomegaer+Ep%nomegaei+Ep%nomegaec ! Total number of frequencies.
2039 
2040  ! ==== Setup of the spectral method ====
2041  Ep%spmeth  =Dtset%spmeth; Ep%nomegasf=Dtset%nomegasf; Ep%spsmear =Dtset%spbroad
2042 
2043  if (Ep%spmeth/=0) then
2044    write(msg,'(2a,i3,2a,i8)')ch10,&
2045 &    ' setup_screening: using spectral method: ',Ep%spmeth,ch10,&
2046 &    ' Number of frequencies for imaginary part: ',Ep%nomegasf
2047    call wrtout(std_out,msg,'COLL')
2048    if (Ep%spmeth==2) then
2049      write(msg,'(a,f8.5,a)')' Gaussian broadening = ',Ep%spsmear*Ha_eV,' [eV]'
2050      call wrtout(std_out,msg,'COLL')
2051    end if
2052  end if
2053 
2054  Ep%nI=1; Ep%nJ=1
2055  if (Dtset%nspinor==2) then
2056    !if (Dtset%usepaw==1.and.Dtset%pawspnorb>0) then
2057    !  Ep%nI=1; Ep%nJ=4
2058    !end if
2059    ! For spin-spin interaction
2060    ! Ep%nI=4; Ep%nJ=4
2061    ABI_CHECK(Ep%npwepG0 == Ep%npwe, "npwepG0 must be == npwe if spinor==2")
2062    !ABI_CHECK(Ep%symchi == 0, "symchi/=0 and nspinor=2 not available")
2063  end if
2064 
2065  ! === Enable the calculations of chi0 on user-specified q-points ===
2066  Ep%nqibz=Qmesh%nibz
2067  ABI_MALLOC(Ep%qibz,(3,Ep%nqibz))
2068  Ep%qibz(:,:)=Qmesh%ibz(:,:)
2069 
2070  Ep%nqcalc=Ep%nqibz
2071  if (Dtset%nqptdm>0) Ep%nqcalc=Dtset%nqptdm
2072 
2073  ABI_MALLOC(Ep%qcalc,(3,Ep%nqcalc))
2074  if (Ep%nqcalc/=Ep%nqibz) then
2075    write(msg,'(6a)')ch10,&
2076 &    ' Dielectric matrix will be calculated only for some ',ch10,&
2077 &    ' selected q points provided by the user through the input variables ',ch10,&
2078 &    ' nqptdm and qptdm'
2079    call wrtout(std_out,msg,'COLL')
2080    call wrtout(ab_out,msg,'COLL')
2081    ltest= Ep%nqcalc <= Qmesh%nibz
2082    ABI_CHECK(ltest, 'nqptdm should not exceed the number of q points in the IBZ')
2083    Ep%qcalc(:,:)=Dtset%qptdm(:,1:Ep%nqcalc)
2084    ! Check whether the q-points provided are correct.
2085    do iq=1,Ep%nqcalc
2086      found=.FALSE.
2087      do iqp=1,Qmesh%nibz
2088        qtmp(:)=Ep%qcalc(:,iq)-Qmesh%ibz(:,iqp)
2089        found=(normv(qtmp,gmet,'G')<GW_TOLQ)
2090        if (found) EXIT
2091      end do
2092      ABI_CHECK(found, 'One or more points specified by Dtset%qptdm do not satisfy q=k1-k2')
2093    end do
2094  else
2095    Ep%qcalc(:,:)=Ep%qibz(:,:)
2096  end if
2097 
2098  ! To write the SCR header correctly, with heads and wings, we have
2099  ! to make sure that q==0, if present, is the first q-point in the list.
2100  !has_q0=(ANY(normv(Ep%qcalc(:,:),gmet,'G')<GW_TOLQ0)) !commented to avoid problems with sunstudio12
2101  has_q0=.FALSE.
2102  do iq=1,Ep%nqcalc
2103    if (normv(Ep%qcalc(:,iq),gmet,'G')<GW_TOLQ0) then
2104      has_q0=.TRUE.; EXIT
2105    end if
2106  end do
2107 
2108  if (has_q0.and.normv(Ep%qcalc(:,1),gmet,'G')>=GW_TOLQ0) then
2109    write(msg,'(5a)')&
2110 &    'The list of q-points to be calculated contains the Gamma point, ',ch10,&
2111 &    'however Gamma is not the first point in the list. ' ,ch10,&
2112 &    'Please, change your input file accordingly. '
2113    MSG_ERROR(msg)
2114  end if
2115 
2116  ! === Initialize the band structure datatype ===
2117  ! * Copy KSS energies and occupations up to Ep%nbnds==Dtset%nband(:)
2118  ! TODO Recheck symmorphy and inversion
2119  bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol))
2120 
2121  ABI_MALLOC(doccde,(bantot))
2122  ABI_MALLOC(eigen,(bantot))
2123  ABI_MALLOC(occfact,(bantot))
2124  doccde(:)=zero; eigen(:)=zero; occfact(:)=zero
2125 
2126  jj=0; ibtot=0
2127  do isppol=1,Dtset%nsppol
2128    do ikibz=1,Dtset%nkpt
2129      do ib=1,Hdr_wfk%nband(ikibz+Dtset%nkpt*(isppol-1))
2130        ibtot=ibtot+1
2131        if (ib<=Ep%nbnds) then
2132          jj=jj+1
2133          occfact(jj)=Hdr_wfk%occ(ibtot)
2134          eigen  (jj)=energies_p(ib,ikibz,isppol)
2135        end if
2136      end do
2137    end do
2138  end do
2139  ABI_FREE(energies_p)
2140 
2141  ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of
2142  ! the symmetry operations in the old GW code (symmorphy and inversion)
2143  ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6))
2144  ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt')
2145 
2146  ABI_MALLOC(npwarr,(Hdr_wfk%nkpt))
2147  npwarr(:)=Ep%npwwfn
2148 
2149  call ebands_init(bantot,KS_BSt,Dtset%nelect,doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,&
2150 & Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,&
2151 & dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
2152 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
2153 
2154  ! TODO modify outkss in order to calculate the eigenvalues also if NSCF calculation.
2155  ! this fails simply because in case of NSCF occ  are zero
2156  !ltest=(ALL(ABS(occfact-KS_BSt%occ)<1.d-2))
2157  !call assert(ltest,'difference in occfact')
2158  !write(std_out,*)MAXVAL(ABS(occfact(:)-KS_BSt%occ(:)))
2159 
2160  !TODO call ebands_update_occ here
2161  !$call ebands_update_occ(KS_BSt,spinmagntarget,stmbias,Dtset%prtvol)
2162 
2163  ABI_FREE(doccde)
2164  ABI_FREE(eigen)
2165  ABI_FREE(npwarr)
2166 
2167  ! Initialize abinit header for the screening part
2168  call hdr_init(KS_BSt,codvsn,Dtset,Hdr_out,Pawtab,pertcase0,Psps,wvl)
2169 
2170  ! Get Pawrhoij from the header.
2171  ABI_DT_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw))
2172  if (Dtset%usepaw==1) then
2173    call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab)
2174    call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij)
2175  end if
2176  call hdr_update(Hdr_out,bantot,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
2177 
2178  ABI_FREE(occfact)
2179  call pawrhoij_free(Pawrhoij)
2180  ABI_DT_FREE(Pawrhoij)
2181 
2182  ! ==== Setup of extrapolar technique ====
2183  Ep%gwcomp = Dtset%gwcomp; Ep%gwencomp = Dtset%gwencomp
2184  if (Ep%gwcomp == 1) then
2185    write(msg,'(a,f8.2,a)')' Using the completeness correction with gwencomp ',Ep%gwencomp*Ha_eV,' [eV] '
2186    call wrtout(std_out,msg,'COLL')
2187  end if
2188 
2189  ! Final compatibility tests
2190  if (ANY(KS_BSt%istwfk /= 1)) then
2191    MSG_WARNING('istwfk/=1 is still under development')
2192  end if
2193 
2194  ltest = (KS_BSt%mband == Ep%nbnds .and. ALL(KS_BSt%nband == Ep%nbnds))
2195  ABI_CHECK(ltest,'BUG in definition of KS_BSt%nband')
2196 
2197  if (Ep%gwcomp==1 .and. Ep%spmeth>0) then
2198    MSG_ERROR("Hilbert transform and extrapolar method are not compatible")
2199  end if
2200 
2201  DBG_EXIT('COLL')
2202 
2203 end subroutine setup_screening