TABLE OF CONTENTS


ABINIT/mrgscr [ Programs ]

[ Top ] [ Programs ]

NAME

 mrgscr

FUNCTION

 This code reads partial (SCR|SUSC) files for different q points creating a single file that
 can be used to perform a sigma calculation.

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (RS, MG, MS)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

NOTES

 If the number of SCR files to be merged is equal to 1, the program checks
 the integrity of the file reporting the list of missing q-points.
 Note that the list of required q-points depends on the k-mesh
 used during the calculation of the KSS file. We assume indeed that the same k-mesh
 is used during the calculation of the matrix elements of sigma.

INPUTS

  (Main program)

OUTPUT

  Only checking and writing

PARENTS

CHILDREN

      abi_io_redirect,abimem_init,abinit_doctor,cqratio,crystal_free
      crystal_from_hdr,cspint,decompose_epsm1,destroy_mpi_enreg
      em1results_free,em1results_print,find_qmesh,flush_unit,fourdp
      get_ppm_eigenvalues,getem1_from_ppm_one_ggp,getng,gsph_free,gsph_init
      hdr_free,herald,hscr_free,hscr_from_file,hscr_print,init_er_from_file
      initmpi_seq,int2char4,ioscr_qmerge,ioscr_qrecover,ioscr_wmerge
      ioscr_wremove,kmesh_free,kmesh_init,kmesh_print,metric,mkdump_er
      pawrhoij_free,ppm_free,ppm_init,prompt,read_rhor,read_screening
      remove_phase,setup_ppmodel,test_charge,timein,vcoul_free,vcoul_init
      wrtout,xmpi_init

SOURCE

  46 #if defined HAVE_CONFIG_H
  47 #include "config.h"
  48 #endif
  49 
  50 #include "abi_common.h"
  51 
  52 program mrgscr
  53 
  54  use defs_basis
  55  use defs_abitypes
  56  use m_xmpi
  57  use m_abicore
  58  use m_build_info
  59  use m_errors
  60  use m_nctk
  61 #ifdef HAVE_NETCDF
  62  use netcdf
  63 #endif
  64  use m_hdr
  65  use m_crystal
  66  use m_pawrhoij
  67 
  68  use m_specialmsg,          only : specialmsg_getcount, herald
  69  use m_time,                only : timein
  70  use m_gwdefs,              only : GW_TOLQ, GW_TOLQ0, GW_Q0_DEFAULT
  71  use m_io_tools,            only : prompt, file_exists, flush_unit, open_file
  72  use m_fstrings,            only : int2char4, endswith, itoa, sjoin, strcat
  73  use m_fft_mesh,            only : g2ifft
  74  use m_fftcore,             only : get_cache_kb, getng
  75  use m_fft,                 only : fourdp
  76  use m_numeric_tools,       only : iseven, cspint
  77  use m_mpinfo,              only : destroy_mpi_enreg, initmpi_seq
  78  use m_geometry,            only : normv, metric
  79  use m_crystal_io,          only : crystal_from_hdr
  80  use m_gsphere,             only : gsph_init, gsph_free, gsphere_t
  81  use m_bz_mesh,             only : kmesh_t, find_qmesh, kmesh_init, kmesh_print, kmesh_free
  82  use m_vcoul,               only : vcoul_t, vcoul_init, vcoul_free
  83  use m_ioarr,               only : read_rhor
  84  use m_io_screening,        only : hscr_io, hscr_print, hscr_copy, read_screening, &
  85                                    write_screening, hscr_free, hscr_t, hscr_from_file,&
  86                                    em1_ncname, chi0_ncname, &
  87                                    ioscr_qmerge, ioscr_qrecover, ioscr_wmerge, ioscr_wremove
  88  use m_ppmodel,             only : ppm_init, ppm_free, setup_ppmodel, getem1_from_PPm_one_ggp, &
  89 &                                  get_PPm_eigenvalues, ppmodel_t, cqratio
  90  use m_model_screening,     only : init_peaks_from_grid,im_screening,re_screening,remove_phase, &
  91 &                                  init_peaks_even_dist,init_single_peak,sequential_fitting, &
  92 &                                  re_and_im_screening_with_phase
  93  use m_levenberg_marquardt, only : lmdif1, lm_fit_print_info
  94  use m_screening,           only : mkdump_er, em1results_free, em1results_print,decompose_epsm1,&
  95 &                                  init_er_from_file, Epsilonm1_results
  96  use m_wfd,                 only : test_charge
  97 
  98 !This section has been created automatically by the script Abilint (TD).
  99 !Do not modify the following lines by hand.
 100 #undef ABI_FUNC
 101 #define ABI_FUNC 'mrgscr'
 102 !End of the abilint section
 103 
 104  implicit none
 105 
 106 !Local variables-------------------------------
 107 !scalars
 108  integer,parameter :: MAX_NUMFILES=100,master=0,paral_kgb0=0,rdwr2=2,prtvol=0,cplex1=1
 109  integer :: iomode,fform1,ifile,ierr,ii,ios,iqibz,iqf,nfiles,timrev
 110  integer :: unt_dump,idx,ig1,ig2,iomega,ppmodel,npwe_asked,mqmem,io,unt_dump2
 111  integer :: id_required,ikxc,approx_type,option_test,dim_kxcg,usexcnhat,usefinegrid
 112  integer :: mgfft,nqlwl,nfft,igmax,comm,nq_selected,kptopt
 113  integer :: choice,nfreq_tot,nfreqre,nfreqim,nfreqc,ifrq,imax
 114  integer :: ig1_start,ig1_end,ig2_start,ig2_end,gmgp_idx,orig_npwe
 115  real(dp) :: ucvol,boxcutmin,ecut,drude_plsmf,compch_fft,compch_sph
 116  real(dp) :: nelectron_exp,freqremax,eps_diff,eps_norm,eps_ppm_norm
 117  real(dp) :: value1,value2,factor,GN_drude_plsmf
 118  real(dp) :: tcpu,tcpui,twall,twalli
 119  real(gwp) :: phase
 120  logical :: is_sus,is_scr,same_freqs,calc_epsilon,only_diag
 121  character(len=1) :: ans
 122  character(len=10) :: tagq
 123  character(len=24) :: codename
 124  character(len=500) :: msg
 125  character(len=nctk_slen) :: varname
 126  character(len=fnlen) :: fname_out,fname,fname_dump,fname_rho,prefix,fname_eigen,fname_dump2
 127  type(hdr_type) :: hdr_rhor
 128  type(abifile_t) :: abifile
 129  type(hscr_t),pointer :: Hscr0
 130  type(hscr_t),target :: Hscr_merge
 131  type(MPI_type) :: MPI_enreg
 132  type(kmesh_t) :: Kmesh,Qmesh
 133  type(crystal_t) :: Cryst
 134  type(gsphere_t)  :: Gsphere
 135  type(ppmodel_t) :: PPm
 136  type(Epsilonm1_results) :: Er
 137  type(vcoul_t),target :: Vcp
 138  type(Dataset_type) :: Dtset
 139 !arrays
 140  integer :: ngfft(18)
 141  integer,allocatable :: foundq(:),freq_indx(:,:)
 142  real(dp),parameter :: k0(3) = [zero,zero,zero]
 143  real(dp) :: gmet(3,3),gprimd(3,3),qdiff(3),rmet(3,3),qtmp(3),tsec(2)
 144  real(dp),allocatable :: qlwl(:,:),real_omega(:),rhor(:,:),rhog(:,:),nhat(:,:)
 145  real(dp),allocatable :: work(:),ftab(:),ysp(:,:),eint(:),qratio(:,:)
 146  complex(gwpc),pointer :: vc_sqrt(:)
 147  complex(gwpc),allocatable :: epsm1(:,:,:,:),kxcg(:,:)
 148  complex(dpc),allocatable :: omega(:),em1_ppm(:),epsm1_eigen(:,:),ppm_eigen(:,:),rhoggp(:,:)
 149  character(len=fnlen),allocatable :: filenames(:)
 150  type(pawrhoij_type),allocatable :: pawrhoij(:)
 151  !type(hscr_t),target,allocatable :: Hscr_file(:) ! Cannot use allocatable as pathscale5 miscompiles the code
 152  type(hscr_t),target :: Hscr_file(MAX_NUMFILES)
 153 
 154 ! *************************************************************************
 155 
 156  ! Change communicator for I/O (mandatory!)
 157  call abi_io_redirect(new_io_comm=xmpi_world)
 158 
 159  ! Initialize MPI
 160  call xmpi_init()
 161 
 162 !Initialize memory profiling if it is activated
 163 !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
 164 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
 165 #ifdef HAVE_MEM_PROFILING
 166  call abimem_init(0)
 167 #endif
 168 
 169  call timein(tcpui,twalli)
 170 
 171  ! Default for sequential use
 172  call initmpi_seq(MPI_enreg); comm = MPI_enreg%comm_world
 173 
 174  is_sus=.FALSE.; is_scr=.FALSE.
 175 
 176 !=== Write greetings, and read the number of files ===
 177  codename='MRGSCR'//REPEAT(' ',18)
 178  call herald(codename,abinit_version,std_out)
 179 
 180  call prompt(' Enter the number of files to merge: ',nfiles)
 181  ABI_CHECK(nfiles>0,'nfiles must be >0')
 182  if (nfiles>MAX_NUMFILES) then
 183    write(msg,"(3a,i0,2a)")&
 184 &   "Do to a bug in the PATHSCALE-v5 compiler, ",ch10,&
 185 &   "the maximum number of files that can be merged is limited to: ",MAX_NUMFILES,ch10,&
 186 &   "Increase the value of MAX_NUMFILES in mrgscr.F90 and recompile"
 187    MSG_ERROR(msg)
 188  end if
 189 
 190  ABI_MALLOC(filenames,(nfiles))
 191  !ABI_DT_MALLOC(Hscr_file,(nfiles))
 192 
 193  if (nfiles == 1) then
 194    call prompt(' Enter the name of the file to be analyzed: ',filenames(1))
 195    write(msg,'(7a)')ch10,&
 196 &   ' Running single-file mode:',ch10,&
 197 &   ' Checking the integrity of file: ',TRIM(filenames(1)),ch10,&
 198 &   ' reporting the list of q-points that are missing. '
 199    call wrtout(std_out,msg,'COLL')
 200 
 201    if (nctk_try_fort_or_ncfile(filenames(1), msg) /= 0) then
 202      MSG_ERROR(msg)
 203    end if
 204 
 205  else if (nfiles > 1) then
 206    ! Read name of files to be merged and check for existence.
 207    call prompt(' Enter the prefix for the final output file: ',fname_out)
 208 
 209    do ifile=1,nfiles
 210      write(msg,'(a,i4)')' Enter the name for the partial screening file no.',ifile
 211      call prompt(msg,filenames(ifile))
 212 
 213      if (nctk_try_fort_or_ncfile(filenames(ifile), msg) /= 0) then
 214        MSG_ERROR(msg)
 215      end if
 216 
 217    end do
 218  end if
 219 
 220  ! Read the header of each file.
 221  do ifile=1,nfiles
 222    iomode = IO_MODE_FORTRAN; if (endswith(filenames(ifile), ".nc")) iomode = IO_MODE_ETSF
 223 
 224    call hscr_from_file(Hscr_file(ifile), filenames(ifile), fform1, comm)
 225    ABI_CHECK(fform1 /= 0, sjoin("fform == 0 in", filenames(ifile)))
 226 
 227    abifile = abifile_from_fform(fform1)
 228    if (abifile%fform == 0) then
 229      MSG_ERROR(sjoin("Cannot find any abifile object associated to fform1:", itoa(fform1)))
 230    end if
 231    if (abifile%class /= "polariz" .and. abifile%class /= "epsm1") then
 232      MSG_ERROR(sjoin('Error while reading header, fform= ',itoa(fform1)))
 233    end if
 234    is_scr = abifile%class == "epsm1"
 235    is_sus = abifile%class == "polariz"
 236 
 237    call hscr_print(Hscr_file(ifile),unit=std_out,prtvol=1)
 238 
 239    if (ifile==1) then
 240      call metric(gmet,gprimd,-1,rmet,Hscr_file(ifile)%Hdr%rprimd,ucvol)
 241    end if
 242  end do !ifile
 243 
 244  if (nfiles>1) then
 245    ! Put the correct ending on the output file
 246    if (is_scr) fname_out=TRIM(fname_out)//'_SCR'
 247    if (is_sus) fname_out=TRIM(fname_out)//'_SUS'
 248  end if
 249 
 250  ! Produce output file in netcdf format we are merging netcdf files.
 251  if (iomode == IO_MODE_ETSF .and. .not. endswith(fname_out, ".nc")) fname_out = nctk_ncify(fname_out)
 252 
 253 !============================
 254 !=== Merge multiple files ===
 255 !============================
 256  if (nfiles>1) then
 257 
 258    ! Check what kind of merging is to be performed
 259    write(std_out,'(2(a))') ch10,' Do you want to merge q-points        (= 1) ?'
 260    write(std_out,'(a)')         '  or do you want to merge frequencies (= 2) ?'
 261    read(std_in,*)choice
 262 
 263    select case(choice)
 264    case(1)
 265      write(std_out,'(3a)') ch10,' 1 => merging q-points',ch10
 266      call ioscr_qmerge(nfiles, filenames, hscr_file, fname_out, hscr_merge)
 267 
 268    case (2)
 269      ! Merge frequencies
 270      write(std_out,'(3a)') ch10,' 2 => merging frequency grids',ch10
 271      !MSG_WARNING("Advanced user option, consistency in fform etc. will not be checked.")
 272 
 273      write(std_out,'(2a)') ch10,' Enter freqremax [eV] for the merged file (Enter 0 to use all freq. found):'
 274      read(std_in,*)freqremax
 275 
 276      freqremax = freqremax/Ha_eV; if (freqremax<tol16) freqremax = HUGE(freqremax)
 277      call ioscr_wmerge(nfiles, filenames, hscr_file, freqremax, fname_out, hscr_merge)
 278 
 279    case default
 280      MSG_ERROR("Invalid choice!")
 281    end select
 282 
 283  end if ! nfiles>1
 284 
 285 !=== Now check if the list of q-points is complete ===
 286 !* Here we assume that the k-mesh reported in the header is the same as that used during the sigma calculation.
 287  write(msg,'(3a)') ch10,' Checking if the list of q-points is complete. ',ch10
 288  call wrtout(std_out,msg,'COLL')
 289 
 290  !call hscr_check_qpoints(hscr0)
 291 
 292  Hscr0 => Hscr_file(1)
 293  fname =filenames(1)
 294  if (nfiles>1) then
 295    Hscr0 => Hscr_merge
 296    fname =fname_out
 297  end if
 298 
 299  timrev=2 ! This should be read from kptopt
 300  call crystal_from_hdr(Cryst,HScr0%Hdr,timrev,remove_inv=.FALSE.)
 301 
 302  kptopt=1
 303  call kmesh_init(Kmesh,Cryst,HScr0%Hdr%nkpt,Hscr0%Hdr%kptns,kptopt)
 304  call kmesh_print(Kmesh,"K-mesh for the wavefunctions",prtvol=prtvol)
 305 
 306  call find_qmesh(Qmesh,Cryst,Kmesh)
 307  call kmesh_print(Qmesh,"Q-mesh for the screening function",prtvol=prtvol)
 308 
 309  ABI_MALLOC(foundq,(Qmesh%nibz))
 310  foundq(:)=0
 311  do iqibz=1,Qmesh%nibz
 312    do iqf=1,Hscr0%nqibz
 313      qdiff(:)=Qmesh%ibz(:,iqibz)-Hscr0%qibz(:,iqf)
 314      if (normv(qdiff,gmet,'G')<GW_TOLQ) foundq(iqibz)=foundq(iqibz)+1
 315    end do
 316  end do
 317 
 318  if (ANY(foundq==0)) then
 319    write(msg,'(6a)')ch10,&
 320 &   ' File ',TRIM(fname),' is not complete ',ch10,&
 321 &   ' The following q-points are missing:'
 322    call wrtout(std_out,msg,'COLL')
 323    ii=0
 324    do iqibz=1,Qmesh%nibz
 325      if (foundq(iqibz)==0) then
 326        ii=ii+1
 327        write(msg,'(i3,a,3f12.6)')ii,') ',Qmesh%ibz(:,iqibz)
 328        call wrtout(std_out,msg,'COLL')
 329      end if
 330    end do
 331  end if
 332 
 333  if (ANY(foundq>1)) then
 334    write(msg,'(6a)')ch10,&
 335 &   ' File ',TRIM(fname),' is overcomplete ',ch10,&
 336 &   ' The following q-points are present more than once:'
 337    call wrtout(std_out,msg,'COLL')
 338    ii=0
 339    do iqibz=1,Qmesh%nibz
 340      if (foundq(iqibz)>1) then
 341        ii=ii+1
 342        write(msg,'(i3,a,3f12.6)')ii,') ',Qmesh%ibz(:,iqibz)
 343        call wrtout(std_out,msg,'COLL')
 344      end if
 345    end do
 346  end if
 347 
 348  if (ALL(foundq==1)) then
 349    write(msg,'(5a)')ch10,&
 350 &   '.File ',TRIM(fname),' contains a complete list of q-points ',ch10
 351    call wrtout(std_out,msg,'COLL')
 352  end if
 353 
 354 !=====================
 355 !=== Recovery mode ===
 356 !=====================
 357  if (nfiles==1) then
 358 
 359    write(std_out,'(2(a))') ch10,' Do you want to recover a subset of q-points    (= 1) ?'
 360    write(std_out,'(a)')         '  or extract the contents of the file           (= 2) ?'
 361    write(std_out,'(a)')         '  or create dielectric function (SCR file)'
 362    write(std_out,'(a)')         '    and/or extract plasmon-pole parameters      (= 3) ?'
 363    write(std_out,'(a)')         '  or remove real frequencies                    (= 4) ?'
 364    write(std_out,'(a)')         '  or remove imaginary frequencies               (= 5) ?'
 365    write(std_out,'(a)')         '  or calculate a model screening                (= 6) ?'
 366    !write(std_out,'(a)')         '  or interpolate a new real freq. grid          (= 7) ?'
 367    !write(std_out,'(a)')         '  or interpolate the screening in k-space       (= 8) ?'
 368    !write(std_out,'(a)')         '  or convert a netcd file to Fortran            (= 9) ?'
 369    read(std_in,*)choice
 370 
 371    select case(choice)
 372    case(1)
 373        ! Recover subset of q-points --------------------------------------------------
 374      write(std_out,'(a)') ' 1 => Recovering subset of q-points'
 375      call prompt(' Enter the number of q-points to be extracted: ',nq_selected)
 376      call prompt(' Enter the name of the final output file: ',fname_out)
 377 
 378      if (endswith(filenames(1), ".nc") .and. .not. endswith(fname_out, ".nc")) then
 379        fname_out = nctk_ncify(fname_out)
 380        call wrtout(std_out,"- Added .nc extension to output file as input data is in netcdf format.")
 381      end if
 382 
 383      call ioscr_qrecover(filenames(1), nq_selected, fname_out)
 384 
 385    case(2)
 386        ! Analyse file ----------------------------------------------------------------
 387      ABI_CHECK(iomode==IO_MODE_FORTRAN, "netcdf output not coded")
 388 
 389      write(std_out,'(a)') ' 2 => Extraction of file contents'
 390 
 391      ! Initialize the G-sphere.
 392      call gsph_init(Gsphere,Cryst,Hscr0%npwe,gvec=Hscr0%gvec)
 393 
 394      ABI_STAT_MALLOC(epsm1,(Hscr0%npwe,Hscr0%npwe,Hscr0%nomega,1), ierr)
 395      ABI_CHECK(ierr==0, 'out of memory in epsm1')
 396 
 397      ! Give option to output epsilon instead of chi0
 398      calc_epsilon = .FALSE.
 399      if (is_sus) then
 400        write(std_out,'(2a)') ch10,&
 401 &       ' You have provided a chi_0 file for analysis. Would you like to output'
 402        write(std_out,'(2a)',advance='no') ' the dielectric function epsilon_GG'' ',&
 403 &       '= delta_GG'' - v_G*chi0_GG''[Y/N] ? '
 404        read(std_in,*)ans
 405 
 406        if (ans=='Y'.or.ans=='y') then
 407          ! Initialise Coulomb terms
 408          if (Er%Hscr%nqlwl==0) then
 409            nqlwl=1
 410            ABI_MALLOC(qlwl,(3,nqlwl))
 411            qlwl(:,1)= GW_Q0_DEFAULT
 412          else
 413            nqlwl=Er%Hscr%nqlwl
 414            ABI_MALLOC(qlwl,(3,nqlwl))
 415            qlwl(:,:)=Er%Hscr%qlwl(:,1:nqlwl)
 416          end if
 417 
 418          Dtset%icutcoul=3; Dtset%rcut=zero
 419          Dtset%vcutgeo=(/zero,zero,zero/);
 420          Dtset%boxcenter=(/zero,zero,zero/)
 421 
 422          write(std_out,'(2a)',advance='no') ch10,' Was a Coulomb cutoff technique used [Y/N] ? '
 423          read(std_in,*)ans
 424          if (ans=='Y'.or.ans=='y') then
 425            write(std_out,'(2a)',advance='no') ' Enter icutcoul: '
 426            read(std_in,*)Dtset%icutcoul
 427            write(std_out,'(2a)',advance='no') ' Enter vcutgeo: '
 428            read(std_in,*)Dtset%vcutgeo
 429            write(std_out,'(2a)',advance='no') ' Enter boxcenter: '
 430            read(std_in,*)Dtset%boxcenter
 431          end if
 432          dtset%ecutsigx = -one
 433 
 434          call vcoul_init(Vcp,Gsphere,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%icutcoul,&
 435 &         Dtset%vcutgeo,Dtset%ecutsigx,Hscr0%npwe,nqlwl,qlwl,ngfft,comm)
 436          ABI_FREE(qlwl)
 437 
 438          calc_epsilon = .TRUE.
 439        end if
 440      end if
 441 
 442      ig1 = 0; ig2 = 0
 443      write(std_out,'(2(a),I0,a)',advance='NO') ch10,' Enter the starting index for G (1 - ',Hscr0%npwe,' ): '
 444      read(std_in,*)ig1_start
 445      if (ig1_start<1.OR.ig1_start>Hscr0%npwe) then
 446        MSG_ERROR(' Starting index out of bounds')
 447      end if
 448      write(std_out,'(a,I0,a,I0,a)',advance='NO')    ' Enter the ending index for G ( ',ig1_start,' - ',Hscr0%npwe,' ): '
 449      read(std_in,*)ig1_end
 450      if (ig1_end<ig1_start.OR.ig1_end>Hscr0%npwe) then
 451        MSG_ERROR(' Ending index out of bounds')
 452      end if
 453      write(std_out,'(a,I0,a)',advance='NO')         ' Enter the starting index for G'' (1 - ',Hscr0%npwe,' ): '
 454      read(std_in,*)ig2_start
 455      if (ig2_start<1.OR.ig2_start>Hscr0%npwe) then
 456        MSG_ERROR(' Starting index out of bounds')
 457      end if
 458      write(std_out,'(a,I0,a,I0,a)',advance='NO')    ' Enter the ending index for G'' ( ',ig2_start,' - ',Hscr0%npwe,' ): '
 459      read(std_in,*)ig2_end
 460      if (ig2_end<ig2_start.OR.ig2_end>Hscr0%npwe) then
 461        MSG_ERROR(' Ending index out of bounds')
 462      end if
 463 
 464      only_diag = .FALSE.
 465      write(std_out,'(a)',advance='no') ' Would you like to output only the diagonal [Y/N] ? '
 466      read(std_in,*)ans
 467      if (ans=='Y'.or.ans=='y') only_diag = .TRUE.
 468 
 469      do iqibz=1,Hscr0%nqibz
 470        ! In the long wavelength limit we set q==0, because we still can use symmetries for the Body.
 471        qtmp(:)=Hscr0%qibz(:,iqibz); if (normv(qtmp,Cryst%gmet,'G')<GW_TOLQ0) qtmp(:)=zero
 472 
 473        ! FIXME
 474        varname = "none"
 475        call read_screening(varname,fname,Hscr0%npwe,1,Hscr0%nomega,epsm1,iomode,comm,iqiA=iqibz)
 476 
 477        if (calc_epsilon) then ! Calculate epsilon
 478          do iomega=1,Hscr0%nomega
 479            if (iqibz==1) then
 480              if (nqlwl>1) then
 481                MSG_ERROR('nqlwl>1 not coded yet!')
 482              end if
 483              vc_sqrt => Vcp%vcqlwl_sqrt(:,iqibz)  ! Use Coulomb term for q-->0
 484            else
 485              vc_sqrt => Vcp%vc_sqrt(:,iqibz)
 486            end if
 487            do ig2=ig2_start,ig2_end
 488              do ig1=ig1_start,ig1_end
 489                epsm1(ig1,ig2,iomega,1) = -(vc_sqrt(ig1)**2)*epsm1(ig1,ig2,iomega,1)
 490              end do ! ig1
 491              epsm1(ig2,ig2,iomega,1) = one + epsm1(ig2,ig2,iomega,1)
 492            end do ! ig2
 493          end do ! iomega
 494        end if ! Do we calculate epsilon
 495 
 496        ! Find out the total number of frequencies along real/imaginary axes
 497        ! and possibly in the z-plane
 498        nfreqre=0; nfreqim=0; nfreqc=0;
 499        do iomega=1,Hscr0%nomega
 500          if (ABS(REAL(Hscr0%omega(iomega)))<tol8.AND.&
 501 &         ABS(AIMAG(Hscr0%omega(iomega)))<tol8) nfreqre = nfreqre + 1
 502          if (ABS(REAL(Hscr0%omega(iomega)))>tol8.AND.&
 503 &         ABS(AIMAG(Hscr0%omega(iomega)))<tol8) nfreqre = nfreqre + 1
 504          if (ABS(REAL(Hscr0%omega(iomega)))<tol8.AND.&
 505 &         ABS(AIMAG(Hscr0%omega(iomega)))>tol8) nfreqim = nfreqim + 1
 506        end do
 507        if (Hscr0%nomega-nfreqre-nfreqim/=0) then
 508          write(std_out,'(/,a)') ' WARNING: There are frequencies in the full complex plane.'
 509          write(std_out,'(a)')   '          The _SCR or _SUS file might not be suitable'
 510          write(std_out,'(a,/)') '          for self-energy calculations.'
 511          nfreqc = Hscr0%nomega-nfreqre-nfreqim
 512        end if
 513        write(std_out,'(2a,I0,a)') ch10,' Found ',Hscr0%nomega,' frequencies.'
 514        write(std_out,'(2(a,I0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10
 515        if (nfreqc>0) then
 516          write(std_out,'(a,I0)') ' There is a grid in the complex plane with ',nfreqc
 517          write(std_out,'(2a)')   '  extra frequencies in the list.',ch10
 518        end if
 519 
 520        ! Get Q index for name
 521        call int2char4(iqibz,tagq)
 522        ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!')
 523 
 524        if (nfreqre>0) then ! Output real frequency axis
 525          if (calc_epsilon) then
 526            fname_dump=TRIM(fname)//'_EPS_Q'//TRIM(tagq)
 527          else
 528            fname_dump=TRIM(fname)//'_Q'//TRIM(tagq)
 529          end if
 530 
 531          if (open_file(fname_dump, msg, newunit=unt_dump, status='replace', form='formatted') /= 0) then
 532            MSG_ERROR(msg)
 533          end if
 534 
 535          do ig1=ig1_start,ig1_end
 536            do ig2=ig2_start,ig2_end
 537              if (only_diag.AND.ig1/=ig2) CYCLE
 538              write(unt_dump,'(2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,/,a,/)')&
 539 &             '# ig1= ',ig1,'    ig2= ',ig2,&
 540 &             '# q = ',Hscr0%qibz(:,iqibz),&
 541 &             '# G = ',Hscr0%gvec(:,ig1),'  G''= ',Hscr0%gvec(:,ig2),&
 542 &             '#   omega [eV]           Re             Im '
 543              do iomega=1,nfreqre
 544                write(unt_dump,'(f8.2,4x,2es16.8)') REAL(Hscr0%omega(iomega))*Ha_eV,&
 545 &               REAL(epsm1(ig1,ig2,iomega,1)),AIMAG(epsm1(ig1,ig2,iomega,1))
 546              end do
 547              write(unt_dump,*)
 548              write(unt_dump,*)
 549            end do !ig2
 550          end do !ig1
 551 !            end if
 552          close(unt_dump)
 553        end if ! Output real frequency axis
 554 
 555        if (nfreqim>0) then ! output imaginary frequency axis
 556          if (calc_epsilon) then
 557            fname_dump=TRIM(fname)//'_EPS_Imfrq_Q'//TRIM(tagq)
 558          else
 559            fname_dump=TRIM(fname)//'_Imfrq_Q'//TRIM(tagq)
 560          end if
 561          if (open_file(fname_dump,msg,newunit=unt_dump,status='replace',form='formatted') /= 0) then
 562            MSG_ERROR(msg)
 563          end if
 564          do ig1=ig1_start,ig1_end
 565            do ig2=ig2_start,ig2_end
 566              if (only_diag.AND.ig1/=ig2) CYCLE
 567              write(unt_dump,'(a,i4,2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,/,a,/)')&
 568 &             '# index= ',idx,'    ig1= ',ig1,'    ig2= ',ig2,&
 569 &             '# q = ',Hscr0%qibz(:,iqibz),&
 570 &             '# G = ',Hscr0%gvec(:,ig1),'  G''= ',Hscr0%gvec(:,ig2),&
 571 &             '#   omega [eV]           Re             Im '
 572              do iomega=nfreqre+1,nfreqre+nfreqim
 573                write(unt_dump,'(f8.2,4x,2es16.8)') AIMAG(Hscr0%omega(iomega))*Ha_eV,epsm1(ig1,ig2,iomega,1)
 574              end do
 575              write(unt_dump,*)
 576              write(unt_dump,*)
 577            end do !ig2
 578          end do !ig1
 579 !            end if
 580          close(unt_dump)
 581        end if ! Check for imaginary frequencies
 582 
 583        ! Check for complex plane values
 584        if (nfreqc>0) then
 585          if (calc_epsilon) then
 586            fname_dump=TRIM(fname)//'_EPS_ZPLANE_Q'//TRIM(tagq)
 587          else
 588            fname_dump=TRIM(fname)//'_ZPLANE_Q'//TRIM(tagq)
 589          end if
 590 
 591          if (open_file(fname_dump,msg,newunit=unt_dump,status='replace',form='formatted') /= 0) then
 592            MSG_ERROR(msg)
 593          end if
 594 
 595          do ig1=ig1_start,ig1_end
 596            do ig2=ig2_start,ig2_end
 597              if (only_diag.AND.ig1/=ig2) CYCLE
 598              write(unt_dump,'(a,i4,2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,/,a,/)')&
 599 &             '# index= ',idx,'    ig1= ',ig1,'    ig2= ',ig2,&
 600 &             '# q = ',Hscr0%qibz(:,iqibz),&
 601 &             '# G = ',Hscr0%gvec(:,ig1),'  G''= ',Hscr0%gvec(:,ig2),&
 602 &             '#   omega [eV]           Re             Im '
 603              do iomega=1,nfreqre
 604                write(unt_dump,'(2(f8.2),4x,2es16.8)') REAL(Hscr0%omega(iomega))*Ha_eV,&
 605 &               AIMAG(Hscr0%omega(iomega))*Ha_eV,epsm1(ig1,ig2,iomega,1)
 606              end do
 607              write(unt_dump,*)
 608              do ios=1,nfreqim
 609                do iomega=1,nfreqre
 610                  if (iomega==1) then
 611                    io = nfreqre + ios
 612                  else
 613                    io = nfreqre + nfreqim + (ios-1)*(nfreqre-1) + (iomega-1)
 614                  end if
 615                  write(unt_dump,'(2(f8.2),4x,2es16.8)') REAL(Hscr0%omega(io))*Ha_eV,&
 616 &                 AIMAG(Hscr0%omega(io))*Ha_eV,epsm1(ig1,ig2,io,1)
 617                end do
 618                write(unt_dump,*)
 619              end do
 620              write(unt_dump,*)
 621              write(unt_dump,*)
 622            end do !ig2
 623          end do !ig1
 624          close(unt_dump)
 625        end if ! Check for complex plane freqs
 626 
 627      end do !iqibz
 628 
 629      ABI_FREE(epsm1)
 630      call gsph_free(Gsphere)
 631 
 632    case(3)
 633        ! Extract dielectric function and plasmon-pole stuff --------------------------
 634      ABI_CHECK(iomode==IO_MODE_FORTRAN, "netcdf output not coded")
 635      write(std_out,'(a)') ' 3 => Calculation of dielectric function and plasmon-pole model'
 636 
 637      npwe_asked=Hscr0%npwe; mqmem=Hscr0%nqibz
 638      call init_Er_from_file(Er,fname,mqmem,npwe_asked,comm)
 639 
 640 !      === Initialize the G-sphere ===
 641      call gsph_init(Gsphere,Cryst,Hscr0%npwe,gvec=Hscr0%gvec)
 642 
 643      boxcutmin=two; igmax=Gsphere%shlim(Gsphere%nsh)
 644      ecut=Er%Hscr%Hdr%ecutdg
 645 
 646      call getng(boxcutmin,ecut,Gsphere%gmet,k0,MPI_enreg%me_fft,&
 647 &     mgfft,nfft,ngfft,MPI_enreg%nproc_fft,Cryst%nsym,paral_kgb0,Cryst%symrel)
 648 
 649 !      I am using standard valued, it would be better to call indefo
 650 !      ngfft(1:3)=Er%Hscr%Hdr%ngfft(1:3)
 651      ngfft(7)=112
 652      ngfft(8)=get_cache_kb()
 653      nfft = PRODUCT(ngfft(1:3))
 654 
 655      Dtset%icutcoul=3; Dtset%rcut=zero
 656      Dtset%vcutgeo=(/zero,zero,zero/); Dtset%boxcenter=(/zero,zero,zero/)
 657      Dtset%ecutsigx = -1
 658 
 659      if (Er%Hscr%nqlwl==0) then
 660        nqlwl=1
 661        ABI_MALLOC(qlwl,(3,nqlwl))
 662        qlwl(:,1)= GW_Q0_DEFAULT
 663      else
 664        nqlwl=Er%Hscr%nqlwl
 665        ABI_MALLOC(qlwl,(3,nqlwl))
 666        qlwl(:,:)=Er%Hscr%qlwl(:,1:nqlwl)
 667      end if
 668 
 669      call vcoul_init(Vcp,Gsphere,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Hscr0%npwe,nqlwl,&
 670 &     qlwl,ngfft,comm)
 671      ABI_FREE(qlwl)
 672 
 673      ! Get the density from an external file ===
 674      ! If meshes are not the same, do an FFT interpolation to have rhor on ngfft.
 675      call prompt(' Enter name for external DEN (or PAWDEN) file: ',fname_rho)
 676 
 677      ABI_MALLOC(rhor,(nfft,Hscr0%Hdr%nspden))
 678      ABI_DT_MALLOC(pawrhoij,(Hscr0%Hdr%natom*Hscr0%Hdr%usepaw))
 679 
 680      call read_rhor(fname_rho, cplex1, nfft, Hscr0%Hdr%nspden, ngfft, 1, MPI_enreg, rhor, hdr_rhor, pawrhoij, comm)
 681 
 682      call hdr_free(hdr_rhor)
 683      call pawrhoij_free(pawrhoij)
 684      ABI_DT_FREE(pawrhoij)
 685 
 686      ABI_MALLOC(rhog,(2,nfft))
 687      call fourdp(1,rhog,rhor(:,1),-1,MPI_enreg,nfft,ngfft,paral_kgb0,0)
 688 
 689      ABI_MALLOC(nhat,(nfft,Hscr0%Hdr%nspden*Hscr0%Hdr%usepaw))
 690      compch_sph=greatest_real; compch_fft=greatest_real
 691      usexcnhat=0; usefinegrid=0
 692 
 693      nelectron_exp = Hscr0%Hdr%nelect
 694 
 695      call test_charge(nfft,nelectron_exp,Hscr0%Hdr%nspden,rhor,Cryst%ucvol,&
 696 &     Hscr0%Hdr%usepaw,usexcnhat,usefinegrid,compch_sph,compch_fft,drude_plsmf)
 697      GN_drude_plsmf = drude_plsmf
 698 
 699      ! Read and in case make Epsilon^{-1} according the the options specified
 700      id_required=4; ikxc=0; approx_type=0; option_test=0; dim_kxcg=0
 701      ABI_MALLOC(kxcg,(nfft,dim_kxcg))
 702 
 703      call prompt(' Enter prefix for output files: ',prefix)
 704      fname_dump=TRIM(prefix)//'_SCR'
 705 
 706      orig_npwe = Er%npwe
 707      write(std_out,'(2a,I0)') ch10,' Number of plane waves is: ',Er%npwe
 708      write(std_out,'(a)',advance='no') ' Would you like to change it [Y/N] ?'
 709      read(std_in,*) ans
 710      if (ans=='Y'.or.ans=='y') then
 711        write(std_out,'(a)',advance='no') ' Enter new no. of plane waves (0 means use old value): '
 712        read(std_in,*) ii
 713        if (ii>0.or.ii<=Er%npwe) Er%npwe = ii
 714        if (ii<0.or.ii>Er%npwe) then
 715          MSG_ERROR(' Wrong value for no. of plane waves!')
 716        end if
 717      end if
 718 
 719      if (is_scr) Er%mqmem=1
 720      if (is_sus) Er%mqmem=0
 721      call mkdump_Er(Er,Vcp,Er%npwe,Gsphere%gvec,dim_kxcg,kxcg,id_required,approx_type,ikxc,option_test,&
 722 &     fname_dump,iomode,nfft,ngfft,comm)
 723      Er%mqmem=1
 724 
 725      call em1results_print(Er)
 726 
 727      write(std_out,'(2a)',advance='no') ch10,&
 728      ' Would you like to calculate the eigenvalues of eps^{-1}_GG''(omega) [Y/N] ? '
 729      read(std_in,*) ans
 730      if (ans=='Y'.or.ans=='y') then
 731 
 732        ABI_MALLOC(epsm1_eigen,(Er%npwe,Er%nomega))
 733        imax = 10
 734        if (Er%npwe < imax) imax = Er%npwe
 735        do iqibz=1,Er%nqibz
 736          call int2char4(iqibz,tagq)
 737          ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!')
 738          fname_eigen=TRIM(prefix)//'_EM1_EIG_Q'//TRIM(tagq)
 739          if (open_file(fname_eigen,msg,newunit=unt_dump,status='replace',form='formatted') /= 0) then
 740            MSG_ERROR(msg)
 741          end if
 742          call decompose_epsm1(Er,iqibz,epsm1_eigen)
 743          write(unt_dump,'(a)')       '# First (max 10) eigenvalues of eps^{-1}(omega)'
 744          write(unt_dump,'(a,3f12.6)')'# q = ',Hscr0%qibz(:,iqibz)
 745          write(unt_dump,'(a)')       '# REAL omega [eV]  REAL(eigen(esp^-1(1,w)))  AIMAG(eigen(esp^-1(1,w))  ...'
 746          do iomega=1,Er%nomega_r
 747            write(unt_dump,'(21(es16.8))')REAL(Er%omega(iomega))*Ha_eV,&
 748 &           (REAL(epsm1_eigen(ii,iomega)),ii=1,imax),(AIMAG(epsm1_eigen(ii,iomega)),ii=1,imax)
 749          end do
 750          close(unt_dump)
 751        end do
 752        ABI_FREE(epsm1_eigen)
 753        ABI_FREE(kxcg)
 754 
 755      end if ! Calculate eigenvalues
 756 
 757      ! Analyze the PPmodel.
 758      write(std_out,'(2a)') ch10,' Would you like to analyse plasmon-pole models [Y/N] ? '
 759      read(std_in,*)ans
 760 
 761      if (ans=='Y'.or.ans=='y') then
 762 
 763        write(std_out,'(2a,f6.2,a)') ch10,' Plasma frequency for GN PPM is: ',GN_drude_plsmf*Ha_eV, ' eV'
 764        write(std_out,'(a)',advance='no') ' Would you like to change it [Y/N] ?'
 765        read(std_in,*) ans
 766        if (ans=='Y'.or.ans=='y') then
 767          write(std_out,'(2a)',advance='no') ch10,' Enter plasma frequency [eV]: '
 768          read(std_in,*) GN_drude_plsmf
 769          GN_drude_plsmf = GN_drude_plsmf/Ha_eV
 770        end if
 771 
 772        write(std_out,'(2a)') ch10,' Would you like to calculate the plasmon-pole model'
 773        write(std_out,'(a)',advance='no')       '        eigenvalues of eps^{-1}_GG''(omega) [Y/N] ? '
 774        read(std_in,*) ans
 775 
 776        if (ans=='Y'.or.ans=='y') then
 777 
 778          ABI_MALLOC(ppm_eigen,(PPm%npwc,Er%nomega))
 779          imax = 10; if (Er%npwe < imax) imax = Er%npwe
 780          do iqibz=1,Er%nqibz
 781            do ppmodel=1,2
 782 
 783              call int2char4(iqibz,tagq)
 784              ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!')
 785              if (ppmodel==1) fname_dump=TRIM(prefix)//'_PPM_GN_EM1_EIG_Q'//TRIM(tagq)
 786              if (ppmodel==2) fname_dump=TRIM(prefix)//'_PPM_HL_EM1_EIG_Q'//TRIM(tagq)
 787              if (ppmodel==3) fname_dump=TRIM(prefix)//'_PPM_vdLH_EM1_EIG_Q'//TRIM(tagq)
 788              if (ppmodel==4) fname_dump=TRIM(prefix)//'_PPM_EF_EM1_EIG_Q'//TRIM(tagq)
 789 
 790              if (open_file(fname_eigen,msg,newunit=unt_dump,status='new',form='formatted') /= 0) then
 791                MSG_ERROR(msg)
 792              end if
 793 
 794              call ppm_free(PPm)
 795              if (ppmodel==1) then
 796                call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,GN_drude_plsmf,Dtset%gw_invalid_freq)
 797              else
 798                call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,drude_plsmf,Dtset%gw_invalid_freq)
 799              end if
 800              call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfft,Gsphere%gvec,ngfft,rhor(:,1),iqibz)
 801 
 802              call get_PPm_eigenvalues(PPm,iqibz,Er%Hscr%zcut,Er%nomega,Er%omega,Vcp,ppm_eigen)
 803 
 804              write(unt_dump,'(a)')       '# First (max 10) eigenvalues of eps^{-1}(omega) from Plasmon-pole model'
 805              write(unt_dump,'(a,3f12.6)')'# q = ',Hscr0%qibz(:,iqibz)
 806              select case(ppmodel)
 807              case(1)
 808                write(unt_dump,'(a)')     '# ppmodel = 1 : Godby - Needs'
 809              case(2)
 810                write(unt_dump,'(a)')     '# ppmodel = 2 : Hybertsen - Louie'
 811              case(3)
 812                write(unt_dump,'(a)')     '# ppmodel = 3 : von der Linden - Horsch'
 813              case(4)
 814                write(unt_dump,'(a)')     '# ppmodel = 4 : Engel - Farid'
 815              end select
 816              write(unt_dump,'(a)')       '# REAL omega [eV]  REAL(eigen(ppm_eps^-1(1,w)))  AIMAG(eigen(ppm_eps^-1(1,w))  ...'
 817              do iomega=1,Er%nomega_r
 818                write(unt_dump,'(21(es16.8))')REAL(Er%omega(iomega))*Ha_eV,&
 819 &               (REAL(ppm_eigen(ii,iomega)),ii=1,imax),(AIMAG(ppm_eigen(ii,iomega)),ii=1,imax)
 820              end do
 821              close(unt_dump)
 822 
 823            end do !ppmodel
 824 
 825          end do ! iqibz
 826          ABI_FREE(ppm_eigen)
 827 
 828        end if ! Calculate PPM eigenvalues
 829 
 830 !        Optionally output eps^{-1}_GG''(w) for a given set of GG' and gridpoints
 831        write(std_out,'(2a)',advance='no') ch10,' Would you like to extract eps^{-1}_GG''(omega) for the PPM [Y/N] ?'
 832        read(std_in,*) ans
 833 
 834        if (ans=='Y'.or.ans=='y') then
 835 !          Reconstruct e^{-1}_GG'(w) according to PPmodel for statistical analysis.
 836          write(std_out,'(a)') ' Enter the number of frequency points in the'
 837          write(std_out,'(a)') '  interval 0 - freqremax (0 means same as input file ): '
 838          read(std_in,*) nfreqre
 839          if (nfreqre==0) then
 840            nfreqre   = Er%nomega_r
 841            nfreqim   = Er%nomega_i
 842            nfreq_tot = Er%nomega
 843            freqremax = REAL(Er%omega(Er%nomega_r))
 844            ABI_MALLOC(omega,(nfreq_tot))
 845            omega(:) = Er%omega(:)
 846            same_freqs = .TRUE.
 847          else
 848            write(std_out,'(a)') ' Enter the value of freqremax (in eV): '
 849            read(std_in,*) freqremax
 850            nfreqim   = Er%nomega_i
 851            nfreq_tot = nfreqre+nfreqim
 852            ABI_MALLOC(omega,(nfreqre+Er%nomega_i))
 853            do iomega=1,nfreqre
 854              omega(iomega) =  CMPLX((freqremax/REAL((nfreqre-1)))*(iomega-1),zero)
 855            end do
 856            omega(nfreqre+1:nfreq_tot) = Er%omega(Er%nomega_r+1:Er%nomega)
 857            same_freqs = .FALSE.
 858          end if ! frequencies
 859 
 860          do iqibz=1,Er%nqibz
 861 
 862            qtmp(:)=Er%qibz(:,iqibz); if (normv(qtmp,Cryst%gmet,'G')<GW_TOLQ0) qtmp(:)=zero
 863            call int2char4(iqibz,tagq)
 864            ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!')
 865 
 866 !            At this time only the Godby-Needs and Hybertsen-Louie models
 867 !            TODO: Check the results from the others
 868            do ppmodel=1,2
 869 
 870              call ppm_free(PPm)
 871              if (ppmodel==1) then
 872                call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,GN_drude_plsmf,Dtset%gw_invalid_freq)
 873              else
 874                call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,drude_plsmf,Dtset%gw_invalid_freq)
 875              end if
 876              call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,&
 877 &             nfft,Gsphere%gvec,ngfft,rhor(:,1),iqibz)
 878 
 879 !              Prepare file for data on real omega axis
 880              if (ppmodel==1) fname_dump=TRIM(prefix)//'_PPM_w_GN_Q'//TRIM(tagq)
 881              if (ppmodel==2) fname_dump=TRIM(prefix)//'_PPM_w_HL_Q'//TRIM(tagq)
 882              if (ppmodel==3) fname_dump=TRIM(prefix)//'_PPM_w_vdLH_Q'//TRIM(tagq)
 883              if (ppmodel==4) fname_dump=TRIM(prefix)//'_PPM_w_EF_Q'//TRIM(tagq)
 884 
 885              if (open_file(fname_dump,msg,newunit=unt_dump,status='replace',form='formatted') /= 0) then
 886                MSG_ERROR(msg)
 887              end if
 888 
 889 !              Prepare file for data on imaginary omega axis
 890              if (ppmodel==1) fname_dump2=TRIM(prefix)//'_PPM_iw_GN_Q'//TRIM(tagq)
 891              if (ppmodel==2) fname_dump2=TRIM(prefix)//'_PPM_iw_HL_Q'//TRIM(tagq)
 892              if (ppmodel==3) fname_dump2=TRIM(prefix)//'_PPM_iw_vdLH_Q'//TRIM(tagq)
 893              if (ppmodel==4) fname_dump2=TRIM(prefix)//'_PPM_iw_EF_Q'//TRIM(tagq)
 894 
 895              if (open_file(fname_dump2,msg,newunit=unt_dump2,status='replace',form='formatted') /= 0) then
 896                MSG_ERROR(msg)
 897              end if
 898 
 899              ABI_MALLOC(em1_ppm,(nfreq_tot))
 900 
 901              ig1 = 0; ig2 = 0
 902              write(std_out,'(3a,I0,a,I0)') ch10,' Enter indices for G and G''.',&
 903 &             'Entering 0 exits the loop. iqibz = ',iqibz,' ppmodel = ',ppmodel
 904 
 905              do
 906                write(std_out,'(2(a),I0,a)',advance='NO') ch10,' Enter index for G (1 - ',Er%npwe,' ): '
 907                read(std_in,*)ig1
 908                if (ig1==0) EXIT
 909                if (ig1<0.OR.ig1>Er%npwe) MSG_ERROR(' index out of bounds')
 910                write(std_out,'(2(a),I0,a)',advance='NO') ch10,' Enter index for G'' (1 - ',Er%npwe,' ): '
 911                read(std_in,*)ig2
 912                if (ig2==0) EXIT
 913                if (ig2<0.OR.ig2>Er%npwe) MSG_ERROR(' index out of bounds')
 914 
 915 !                Generate the PPM representation of epsilon^-1
 916                call getem1_from_PPm_one_ggp(PPm,iqibz,Er%Hscr%zcut,nfreq_tot,omega,Vcp,em1_ppm,ig1,ig2)
 917 
 918                write(unt_dump,'(a,I1)') '# epsilon^-1_GG''(omega) from ppmodel = ',ppmodel
 919                write(unt_dump,'(2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,&
 920 &               /,a,3F9.4,a,3F9.4,a,/a,f9.4,a,f9.4,a,/,a,/)')&
 921 &               '# ig1= ',ig1,'    ig2= ',ig2,&
 922 &               '# q = ',Er%qibz(:,iqibz),&
 923 &               '# G = ',Er%gvec(:,ig1),'  G''= ',Er%gvec(:,ig2),&
 924 &               '# G = (',MATMUL(two_pi*Cryst%gmet,Er%gvec(:,ig1)),&
 925 &               ')  G''= (',MATMUL(two_pi*Cryst%gmet,Er%gvec(:,ig2)),')',&
 926 &               '# 1/2|G|^2 =',half*normv(Er%gvec(:,ig1),Cryst%gmet,'G')**2,&
 927 &               ' Ha 1/2|G''|^2 =',half*normv(Er%gvec(:,ig2),Cryst%gmet,'G')**2,' Ha',&
 928 &               '#   omega [eV]           Re             Im '
 929                write(unt_dump2,'(a,I1)') '# epsilon^-1_GG''(iomega) from ppmodel = ',ppmodel
 930                write(unt_dump2,'(2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,&
 931 &               /,a,3F9.4,a,3F9.4,a,/a,f9.4,a,f9.4,a,/,a,/)')&
 932 &               '# ig1= ',ig1,'    ig2= ',ig2,&
 933 &               '# q = ',Er%qibz(:,iqibz),&
 934 &               '# G = ',Er%gvec(:,ig1),'  G''= ',Er%gvec(:,ig2),&
 935 &               '# G = (',MATMUL(two_pi*Cryst%gmet,Er%gvec(:,ig1)),&
 936 &               ')  G''= (',MATMUL(two_pi*Cryst%gmet,Er%gvec(:,ig2)),')',&
 937 &               '# 1/2|G|^2 =',half*normv(Er%gvec(:,ig1),Cryst%gmet,'G')**2,&
 938 &               ' Ha 1/2|G''|^2 =',half*normv(Er%gvec(:,ig2),Cryst%gmet,'G')**2,' Ha',&
 939 &               '#   iomega [eV]           Re             Im '
 940 
 941                do iomega=1,nfreqre
 942                  if (same_freqs) then
 943                    write(unt_dump,'(f8.2,4x,4es16.8)') REAL(omega(iomega))*Ha_eV,em1_ppm(iomega),&
 944 &                   Er%epsm1(ig1,ig2,iomega,iqibz)
 945                  else
 946                    write(unt_dump,'(f8.2,4x,2es16.8)') REAL(omega(iomega))*Ha_eV,em1_ppm(iomega)
 947                  end if
 948                end do
 949 !                First output the iomega = 0 point
 950                write(unt_dump2,'(f8.2,4x,4es16.8)') AIMAG(omega(1))*Ha_eV,em1_ppm(1),&
 951 &               Er%epsm1(ig1,ig2,1,iqibz)
 952 !                Then the rest
 953                do iomega=nfreqre+1,nfreq_tot
 954                  write(unt_dump2,'(f8.2,4x,4es16.8)') AIMAG(omega(iomega))*Ha_eV,em1_ppm(iomega),&
 955 &                 Er%epsm1(ig1,ig2,iomega,iqibz)
 956                end do
 957                write(unt_dump,*)
 958                write(unt_dump,*)
 959                write(unt_dump2,*)
 960                write(unt_dump2,*)
 961              end do ! Empty
 962              ABI_FREE(em1_ppm)
 963              close(unt_dump)
 964              close(unt_dump2)
 965 
 966            end do ! ppmodel
 967 
 968          end do ! iqibz
 969          ABI_FREE(omega)
 970        end if ! Output epsilon for PPM
 971 
 972 !        Optionally statistics for all PPMs
 973        write(std_out,'(2a)',advance='no') ch10,' Would you like to output statistics for all PPMs [Y/N] ?'
 974        read(std_in,*) ans
 975 
 976        if (ans=='Y'.or.ans=='y') then
 977          nfreqre   = Er%nomega_r
 978          nfreq_tot = Er%nomega
 979          freqremax = REAL(Er%omega(Er%nomega_r))
 980          ABI_MALLOC(real_omega,(nfreqre))
 981          real_omega(:) = REAL(Er%omega(1:nfreqre))
 982 
 983          do iqibz=1,Er%nqibz
 984            do ppmodel=1,2
 985 
 986              qtmp(:)=Er%qibz(:,iqibz)
 987              if (normv(qtmp,Cryst%gmet,'G')<GW_TOLQ0) qtmp(:)=zero
 988 
 989              call ppm_free(PPm)
 990              if (ppmodel==1) then
 991                call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,GN_drude_plsmf,Dtset%gw_invalid_freq)
 992              else
 993                call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,drude_plsmf,Dtset%gw_invalid_freq)
 994              end if
 995              call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,&
 996              nfft,Gsphere%gvec,ngfft,rhor(:,1),iqibz)
 997 
 998              ! Prepare ratios and density for the f-sum rule
 999              ABI_STAT_MALLOC(qratio,(orig_npwe,orig_npwe), ierr)
1000              ABI_CHECK(ierr==0, "out-of-memory in qratio")
1001 
1002              ABI_STAT_MALLOC(rhoggp,(Er%npwe,Er%npwe), ierr)
1003              ABI_CHECK(ierr==0, "out-of-memory in rhoggp")
1004 
1005              call cqratio(orig_npwe,Gsphere%gvec,qtmp,Cryst%gmet,Cryst%gprimd,qratio)
1006              ! Arrange n(G-G')->n(G,G')
1007              ierr=0
1008              do ig1=1,Er%npwe
1009                do ig2=1,Er%npwe
1010                  gmgp_idx = g2ifft(Gsphere%gvec(:,ig1)-Gsphere%gvec(:,ig2),ngfft)
1011                  if (gmgp_idx/=0) then
1012                    rhoggp(ig1,ig2)=CMPLX(rhog(1,gmgp_idx),rhog(2,gmgp_idx))
1013                  else
1014                    ierr=ierr+1
1015                    rhoggp(ig1,ig2)=czero
1016                  end if
1017                end do
1018              end do
1019              if (ierr/=0) then
1020                write(std_out,'(a,i0,a)')' Found ',ierr,' G1-G2 vectors falling outside the FFT box. '
1021              end if
1022 
1023              ! Prepare files
1024              call int2char4(iqibz,tagq)
1025              ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!')
1026              if (ppmodel==1) fname_dump=TRIM(prefix)//'_norms_GN_Q'//TRIM(tagq)
1027              if (ppmodel==2) fname_dump=TRIM(prefix)//'_norms_HL_Q'//TRIM(tagq)
1028              if (open_file(fname_dump,msg, newunit=unt_dump, status='replace',form='formatted') /= 0) then
1029                MSG_ERROR(msg)
1030              end if
1031              write(unt_dump,'(a)') '# Various norms integrated through spline interpolation'
1032              write(unt_dump,'(a)') '# over all frequencies in the input file,'
1033              write(unt_dump,'(a)') '# for all G and G'' vectors.'
1034              write(unt_dump,'(a,I0)') '#               ppmodel: ',ppmodel
1035              write(unt_dump,'(a,I0)') '# Number of frequencies: ',nfreqre
1036              write(unt_dump,'(a,f12.6)') '# Maximum frequency    : ',freqremax
1037              write(unt_dump,'(a)') '# Columns:'
1038              write(unt_dump,'(2a)') '#  ig1      ig2   |eps-eps_PPM|/|eps|',&
1039 &             '   |eps-eps_PPM|    |eps|     |eps_PPM|            G                  G'''
1040              if (ppmodel==1) fname_dump2=TRIM(prefix)//'_f_sumrule_GN_Q'//TRIM(tagq)
1041              if (ppmodel==2) fname_dump2=TRIM(prefix)//'_f_sumrule_HL_Q'//TRIM(tagq)
1042 
1043              if (open_file(fname_dump2,msg,newunit=unt_dump2,status='replace',form='formatted') /= 0) then
1044                MSG_ERROR(msg)
1045              end if
1046 
1047              write(unt_dump2,'(a)') '# The fulfillment of the f-sum rule: I(epsilon) ='
1048              write(unt_dump2,'(a)') '#   int_0^{inf}{omega*Im[epsilon_G,G''(omega)]}/C_qGG'''
1049              write(unt_dump2,'(a)') '# C_qGG'' = '
1050              write(unt_dump2,'(a)') '#   -Pi/2*omega_p^2*(q+G)*(q+G'')/|q+G|^2*n(G-G'')/n(0)'
1051              write(unt_dump2,'(a)') '# for all G and G'' vectors.'
1052              write(unt_dump2,'(a,I0)') '#               ppmodel: ',ppmodel
1053              write(unt_dump2,'(a,I0)') '# Number of frequencies: ',nfreqre
1054              write(unt_dump2,'(a,f12.6)') '# Maximum frequency    : ',freqremax
1055              write(unt_dump2,'(a)') '# Columns:'
1056              write(unt_dump2,'(3a)') '#  ig1      ig2   I(epsilon)',&
1057 &             '   I(eps_PPM)   Re[n(G-G'')]    Im[n(G-G'')]    qratio      I1*C_qGG''',&
1058 &             ' Re[Omegatwsq] Im[Omegatwsq]   Re[omegatw]   Im[omegatw]    |G|    1/2|G|^2'
1059 
1060              ABI_MALLOC(em1_ppm,(nfreq_tot))
1061              ABI_MALLOC(ftab,(nfreqre))
1062              ABI_MALLOC(ysp,(3,nfreqre))
1063              ABI_MALLOC(work,(nfreqre))
1064              ABI_MALLOC(eint,(nfreqre))
1065 
1066              do ig1=1,Er%npwe
1067                write(std_out,'(2(a,I0))') ' ig1= ',ig1, ' of ',Er%npwe
1068                do ig2=1,Er%npwe
1069                  !ig2 = ig1
1070                  call getem1_from_PPm_one_ggp(PPm,iqibz,Er%Hscr%zcut,nfreq_tot,Er%omega,Vcp,&
1071 &                 em1_ppm,ig1,ig2)
1072 
1073                  ! Calculate norms in real
1074                  eps_diff=0; eps_norm=0; eps_ppm_norm=0
1075                  ftab(1:nfreqre) = ABS(Er%epsm1(ig1,ig2,1:nfreqre,iqibz)-em1_ppm(1:nfreqre))
1076                  call cspint(ftab,real_omega,nfreqre,real_omega(1),&
1077 &                 real_omega(nfreqre),ysp,eint,work,eps_diff)
1078                  ftab(1:nfreqre) = ABS(Er%epsm1(ig1,ig2,1:nfreqre,iqibz))
1079                  call cspint(ftab,real_omega,nfreqre,real_omega(1),&
1080 &                 real_omega(nfreqre),ysp,eint,work,eps_norm)
1081                  ftab(1:nfreqre) = ABS(em1_ppm(1:nfreqre))
1082                  call cspint(ftab,real_omega,nfreqre,real_omega(1),&
1083 &                 real_omega(nfreqre),ysp,eint,work,eps_ppm_norm)
1084                  write(unt_dump,'(2i6,f12.4,3es14.4,6i4)') ig1,ig2,eps_diff/eps_norm,eps_diff,&
1085 &                 eps_norm,eps_ppm_norm,Er%gvec(:,ig1),Er%gvec(:,ig2)
1086 
1087                  ! Evaluate the f-sum rule
1088                  if (ig1==ig2) then
1089                    ftab(1:nfreqre) = real_omega(1:nfreqre)*AIMAG(Er%epsm1(ig1,ig2,1:nfreqre,iqibz))
1090                  else ! Dephase first - HERE epsm1 is changed!
1091                    call remove_phase(epsm1(ig1,ig2,:,1),Hscr_file(1)%nomega,phase)
1092                    ftab(1:nfreqre) = real_omega(1:nfreqre)*AIMAG(Er%epsm1(ig1,ig2,1:nfreqre,iqibz))
1093                  end if
1094 
1095                  call cspint(ftab,real_omega,nfreqre,real_omega(1),&
1096 &                 real_omega(nfreqre),ysp,eint,work,eps_diff)
1097 
1098                  if (ig1==ig2) then
1099                    factor = -two*pi*pi*REAL(rhoggp(ig1,ig2))*qratio(ig1,ig2)
1100                  else
1101                    rhoggp(ig1,ig2) = CMPLX(COS(phase),-SIN(phase))*rhoggp(ig1,ig2)
1102                    factor = -two*pi*pi*REAL(rhoggp(ig1,ig2))*qratio(ig1,ig2)
1103                  end if
1104 
1105                  if (ABS(qratio(ig1,ig2))>zero) then
1106                    value1 = eps_diff/factor
1107                    if (ppmodel==1) then
1108                      value2 = -pi*half*(REAL(PPm%bigomegatwsq(iqibz)%vals(ig1,ig2))&
1109 &                     /(REAL(PPm%omegatw(iqibz)%vals(ig1,ig2))))&
1110 &                     /factor*(2*sqrt(pi*rhoggp(1,1)))
1111                    else
1112                      value2 = -pi*half*(SQRT(REAL(PPm%bigomegatwsq(iqibz)%vals(ig1,ig2))))&
1113 &                     /factor*(2*sqrt(pi*rhoggp(1,1)))
1114                    end if
1115                  else
1116                    value1 = zero
1117                    value2 = zero
1118                  end if
1119 
1120                  write(unt_dump2,'(2i6,12es14.4)') ig1,ig2,value1,value2,&
1121 &                 REAL(rhoggp(ig1,ig2)),AIMAG(rhoggp(ig1,ig2)),qratio(ig1,ig2),&
1122 &                 eps_diff,REAL(PPm%bigomegatwsq(iqibz)%vals(ig1,ig2)),&
1123 &                 AIMAG(PPm%bigomegatwsq(iqibz)%vals(ig1,ig2)),&
1124 &                 REAL(PPm%omegatw(iqibz)%vals(ig1,ig2)),&
1125 &                 AIMAG(PPm%omegatw(iqibz)%vals(ig1,ig2)),&
1126 &                 normv(Er%gvec(:,ig1),Cryst%gmet,'G'),&
1127 &                 half*normv(Er%gvec(:,ig1),Cryst%gmet,'G')**2
1128 
1129                end do !ig2
1130              end do !ig1
1131 
1132              ABI_FREE(em1_ppm)
1133              ABI_FREE(ftab)
1134              ABI_FREE(ysp)
1135              ABI_FREE(work)
1136              ABI_FREE(eint)
1137              ABI_FREE(qratio)
1138              ABI_FREE(rhoggp)
1139              close(unt_dump); close(unt_dump2)
1140 
1141            end do ! ppmodel
1142          end do ! iqibz
1143 
1144          ABI_FREE(real_omega)
1145        end if ! Output statistics
1146 
1147        call ppm_free(PPm)
1148      end if ! If ppmodel>0
1149 
1150      ABI_FREE(rhor)
1151      ABI_FREE(rhog)
1152      ABI_FREE(nhat)
1153 
1154      call vcoul_free(Vcp)
1155      call em1results_free(Er)
1156      call gsph_free(Gsphere)
1157 
1158    case(4)
1159      ! Remove real frequencies ----------------------------------------------------------
1160      write(std_out,'(2(a))') ch10,' Do you want to remove every other real frequency  (= 1) ?'
1161      write(std_out,'(a)')         '  or specify for each real frequency individually  (= 2) ?'
1162      write(std_out,'(a)')         '  or remove ALL real frequencies                   (= 3) ?'
1163      read(std_in,*)choice
1164 
1165        ! Calculate the total number of real freq
1166      nfreqre = 0; nfreqim = 0
1167      do ifrq=1,Hscr_file(1)%nomega
1168          ! If frequency is not imaginary, count.
1169        if (AIMAG(Hscr_file(1)%omega(ifrq)) < tol8) nfreqre = nfreqre + 1
1170        if (REAL(Hscr_file(1)%omega(ifrq)) < tol8 .and. AIMAG(Hscr_file(1)%omega(ifrq))>tol8)  nfreqim = nfreqim + 1
1171      end do
1172 
1173      nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq
1174      write(std_out,'(2a,I0,a)') ch10,' Found ',nfreq_tot,' frequencies.'
1175      write(std_out,'(2(a,I0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10
1176 
1177      ! Array with the index of frequencies to be kept.
1178      ABI_MALLOC(freq_indx,(nfreq_tot,1))
1179      freq_indx = 0
1180 
1181      select case(choice)
1182      case(1)
1183        ! Remove every other frequency
1184        write(std_out,'(2(a))') ch10,' Removing every other real frequency, i.e. every even one.'
1185        write(std_out,'(a)')         ' If the total number of frequencies is odd, the first and last one will be kept.'
1186        write(std_out,'(a)')         ' If the total number is even, the first one will still be in the final set.'
1187 
1188        ! Test for no real frequencies
1189        ABI_CHECK(nfreqre /= 0, "No real frequencies in file!")
1190 
1191        ii=nfreqre; nfreqre = 0
1192        do ifrq=1,ii
1193          if (.not. iseven(ifrq)) then
1194            nfreqre = nfreqre + 1; freq_indx(nfreqre,1) = ifrq
1195          end if
1196        end do
1197        write(std_out,'(2a,I0,a)') ch10,' ',nfreqre,' real frequencies will be kept.'
1198 
1199      case(2)
1200          ! Specify freq. individually
1201        ii = nfreqre; nfreqre = 0
1202        do ifrq=1,ii
1203          write(std_out,'(a,f12.6,a)') ' Would you like to keep freq. at: ',REAL(Hscr_file(1)%omega(ifrq))*Ha_eV,' eV? [y/n]'
1204          read(std_in,*) ans
1205          if (ans=='Y'.or.ans=='y') then
1206            nfreqre = nfreqre + 1; freq_indx(nfreqre,1) = ifrq
1207          end if
1208        end do
1209        write(std_out,'(2a,I0,a)') ch10,' ',nfreqre,' real frequencies will be kept.'
1210 
1211      case(3)
1212        ! Remove all real freq.
1213        nfreqre = 0
1214 
1215      case default
1216        MSG_ERROR("Invalid choice!")
1217      end select
1218 
1219      ! Add imaginary frequencies if any
1220      if (nfreqim > 0) then
1221        nfreqim = 0
1222        do ifrq=1,Hscr_file(1)%nomega
1223          if (AIMAG(Hscr_file(1)%omega(ifrq)) > tol8) then
1224            nfreqim = nfreqim + 1; freq_indx(nfreqre+nfreqim,1) = ifrq
1225          end if
1226        end do
1227      end if
1228 
1229      nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq
1230      write(std_out,'(3(a,i0),a)')' Finally, we have ',nfreq_tot,' frequencies. ',nfreqre,' real, and ',nfreqim,' imaginary.'
1231 
1232      call prompt(' Enter the full name of the final output file: ',fname_out)
1233 
1234      if (endswith(filenames(1), ".nc") .and. .not. endswith(fname_out, ".nc")) then
1235        fname_out = nctk_ncify(fname_out)
1236        call wrtout(std_out,"- Added .nc extension to output file as input data is in netcdf format.")
1237      end if
1238 
1239      call ioscr_wremove(filenames(1), hscr_file(1), fname_out, nfreq_tot, freq_indx, hscr_merge)
1240      ABI_FREE(freq_indx)
1241 
1242    case(5)
1243      ! Remove imaginary frequencies ------------------------------------------------
1244 
1245      ! Calculate the total number of freq
1246      nfreqre = 0; nfreqim = 0
1247      do ifrq=1,Hscr_file(1)%nomega
1248          ! If frequency is not imaginary, count.
1249        if (AIMAG(Hscr_file(1)%omega(ifrq))<tol8) nfreqre = nfreqre + 1
1250        if (REAL(Hscr_file(1)%omega(ifrq))<tol8.and.AIMAG(Hscr_file(1)%omega(ifrq))>tol8)  nfreqim = nfreqim + 1
1251      end do ! ifrq
1252 
1253        ! Test for no real frequencies
1254      if (nfreqim==0) then
1255        MSG_ERROR("No imaginary frequencies in file!")
1256      end if
1257 
1258      nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq
1259      write(std_out,'(2a,I0,a)') ch10,' Found ',nfreq_tot,' frequencies.'
1260      write(std_out,'(2(a,I0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10
1261 
1262      ABI_MALLOC(freq_indx,(nfreq_tot,1))
1263      freq_indx = 0
1264 
1265      ! Specify freq. individually
1266      ii=nfreq_tot; nfreqim = 0
1267      do ifrq=nfreqre+1,ii
1268        write(std_out,'(a,f12.6,a)')&
1269 &       ' Would you like to keep imaginary freq. at: ',AIMAG(Hscr_file(1)%omega(ifrq))*Ha_eV,' eV? [y/n]'
1270        read(std_in,*) ans
1271        if (ans=='Y'.or.ans=='y') then
1272          nfreqim = nfreqim + 1; freq_indx(nfreqre+nfreqim,1) = ifrq
1273        end if
1274      end do ! ifrq
1275      write(std_out,'(2a,I0,a)') ch10,' ',nfreqim,' imaginary frequencies will be kept.'
1276 
1277      ! Add real frequencies if any
1278      if (nfreqre > 0) then
1279        nfreqre = 0
1280        do ifrq=1,Hscr_file(1)%nomega
1281          if (AIMAG(Hscr_file(1)%omega(ifrq)) < tol8) then
1282            nfreqre = nfreqre + 1; freq_indx(nfreqre,1) = ifrq
1283          end if
1284        end do
1285      end if
1286 
1287      nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq
1288      write(std_out,'(2a,I0,a)') ch10,' Finally, we have ',nfreq_tot,' frequencies.'
1289      write(std_out,'(2(a,I0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10
1290 
1291      call prompt(' Enter the full name of the final output file: ',fname_out)
1292 
1293      if (endswith(filenames(1), ".nc") .and. .not. endswith(fname_out, ".nc")) then
1294        fname_out = nctk_ncify(fname_out)
1295        call wrtout(std_out,"- Added .nc extension to output file as input data is in netcdf format.")
1296      end if
1297 
1298      call ioscr_wremove(filenames(1), hscr_file(1), fname_out, nfreq_tot, freq_indx, hscr_merge)
1299 
1300      ABI_FREE(freq_indx)
1301 
1302    case(6)
1303      ! Model screening -------------------------------------------------------------
1304      MSG_ERROR("Model screening has been removed")
1305 
1306    case (7)
1307      ! Interpolation for real frequency -------------------------------------
1308      MSG_ERROR("Still under development")
1309 
1310    !case(9)
1311    !  TODO
1312    !  ! netcdf --> Fortran converter -------------------------------------------------------------
1313    !  call prompt(' Enter the name of the final output Fortran file: ',fname_out)
1314    !  ! fname_out extension should be consistent with filenames(1)
1315    !  call ioscr_nc2fort(filenames(1), fname_out)
1316 
1317    case default
1318        ! Bail if choice is wrong
1319      write(std_out,*) ' Invalid choice! Exiting...'
1320      goto 100
1321    end select
1322 
1323  end if ! Single file mode
1324 
1325  call timein(tcpu,twall)
1326 
1327  tsec(1)=tcpu-tcpui
1328  tsec(2)=twall-twalli
1329 
1330  write(std_out, '(a,a,a,f13.1,a,f13.1)' ) &
1331 & '-',ch10,'- Proc.   0 individual time (sec): cpu=',tsec(1),'  wall=',tsec(2)
1332 
1333 !=====================
1334 !==== Free memory ====
1335 !=====================
1336  ABI_FREE(filenames)
1337 
1338  if (allocated(kxcg)) then
1339    ABI_FREE(kxcg)
1340  end if
1341  if (allocated(foundq)) then
1342    ABI_FREE(foundq)
1343  end if
1344 
1345  call crystal_free(Cryst)
1346  call kmesh_free(Kmesh)
1347  call kmesh_free(Qmesh)
1348  call destroy_mpi_enreg(MPI_enreg)
1349 
1350  nullify(Hscr0)
1351  !if (nfiles>1) call hscr_free(Hscr_merge)
1352  call hscr_free(Hscr_merge)
1353 
1354  do ifile=1,nfiles
1355    call hscr_free(Hscr_file(ifile))
1356  end do
1357  !ABI_DT_FREE(Hscr_file)
1358 
1359  call flush_unit(std_out)
1360 
1361  call abinit_doctor("__mrgscr")
1362 
1363  100 call xmpi_end()
1364 
1365  end program mrgscr