TABLE OF CONTENTS


ABINIT/m_io_screening [ Modules ]

[ Top ] [ Modules ]

NAME

  m_io_screening

FUNCTION

  This module contains the definition of the header of the
  _SCR and _SUSC file as well as methods used to read/write/echo.

COPYRIGHT

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

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_io_screening
24 
25  use defs_basis
26  use m_abicore
27 #if defined HAVE_MPI2
28  use mpi
29 #endif
30  use m_xmpi
31  use m_mpiotk
32  use m_nctk
33  use m_errors
34  use m_dtset
35  use, intrinsic :: iso_c_binding
36  use netcdf
37  use m_hdr
38  use m_sort
39 
40  use m_gwdefs,          only : em1params_t, GW_TOLQ
41  use m_fstrings,        only : sjoin, itoa, endswith, replace_ch0
42  use m_copy,            only : alloc_copy
43  use m_io_tools,        only : open_file, file_exists, iomode2str
44  use m_numeric_tools,   only : print_arr, remove_copies, imax_loc
45  use m_bz_mesh,         only : isequalk
46 
47  implicit none
48 
49  private

m_io_screening/hscr_bcast [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 hscr_bcast

FUNCTION

 This subroutine transmit the header structured datatype initialized
 on one processor (or a group of processor), to the other processors.
 It also allocates the needed part of the header.

INPUTS

  master=ID of the master node.
  my_rank=ID of the node that receives the data.
  comm=MPI communicator.

OUTPUT

  (no output)

SIDE EFFECTS

  Hscr<type(hscr_t)>=the SCR header. For the master, it is already
   initialized entirely, while for the other procs, everything has
   to be transmitted.

NOTES

 This routine is called only in the case of MPI version of the code.

SOURCE

931 subroutine hscr_bcast(hscr, master, my_rank, comm)
932 
933 !Arguments ------------------------------------
934  class(hscr_t),intent(inout) :: hscr
935  integer, intent(in) :: master, my_rank, comm
936 
937 !Local variables-------------------------------
938  integer :: ierr
939 
940 ! *************************************************************************
941 
942  DBG_ENTER("COLL")
943  if (xmpi_comm_size(comm) == 1) return ! Nothing to do
944 
945  ! integer
946  call xmpi_bcast(hscr%id,         master,comm,ierr)
947  call xmpi_bcast(hscr%ikxc,       master,comm,ierr)
948  call xmpi_bcast(hscr%inclvkb,    master,comm,ierr)
949  call xmpi_bcast(hscr%headform,   master,comm,ierr)
950  call xmpi_bcast(hscr%fform,      master,comm,ierr)
951  call xmpi_bcast(hscr%gwcalctyp,  master,comm,ierr)
952  call xmpi_bcast(hscr%nI,         master,comm,ierr)
953  call xmpi_bcast(hscr%nJ,         master,comm,ierr)
954  call xmpi_bcast(hscr%nqibz,      master,comm,ierr)
955  call xmpi_bcast(hscr%nqlwl,      master,comm,ierr)
956  call xmpi_bcast(hscr%nomega,     master,comm,ierr)
957  call xmpi_bcast(hscr%nbnds_used, master,comm,ierr)
958  call xmpi_bcast(hscr%npwe,       master,comm,ierr)
959  call xmpi_bcast(hscr%npwwfn_used,master,comm,ierr)
960  call xmpi_bcast(hscr%spmeth,     master,comm,ierr)
961  call xmpi_bcast(hscr%test_type,  master,comm,ierr)
962  call xmpi_bcast(hscr%tordering,  master,comm,ierr)
963 
964  ! Real
965  call xmpi_bcast(hscr%mbpt_sciss, master,comm,ierr)
966  call xmpi_bcast(hscr%spsmear,    master,comm,ierr)
967  call xmpi_bcast(hscr%zcut,       master,comm,ierr)
968 
969  ! arrays
970  call xmpi_bcast(hscr%titles, master,comm,ierr)
971 
972  if (my_rank /= master) call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
973 
974  call xmpi_bcast(hscr%gvec, master,comm,ierr)
975  call xmpi_bcast(hscr%qibz, master,comm,ierr)
976  call xmpi_bcast(hscr%qlwl, master,comm,ierr)
977  call xmpi_bcast(hscr%omega,master,comm,ierr)
978 
979  ! Communicate the Abinit header.
980  call hscr%Hdr%bcast(master, my_rank, comm)
981 
982  ! HSCR_NEW
983  call xmpi_bcast(hscr%awtr, master, comm, ierr)
984  call xmpi_bcast(hscr%icutcoul, master, comm, ierr)
985  call xmpi_bcast(hscr%vcutgeo, master, comm, ierr)
986  call xmpi_bcast(hscr%gwcomp, master, comm, ierr)
987  call xmpi_bcast(hscr%gwgamma, master, comm, ierr)
988  call xmpi_bcast(hscr%gwencomp, master, comm, ierr)
989  call xmpi_bcast(hscr%kind_cdata, master, comm, ierr)
990  ! HSCR_NEW
991 
992  DBG_EXIT("COLL")
993 
994 end subroutine hscr_bcast

m_io_screening/hscr_copy [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 hscr_copy

FUNCTION

 Deep copy of the header of the _SCR or _SUSC file.

INPUTS

SOURCE

1073 subroutine hscr_copy(Hscr_in, Hscr_cp)
1074 
1075 !Arguments ------------------------------------
1076 !scalars
1077  class(hscr_t),intent(in) :: Hscr_in
1078  class(hscr_t),intent(inout) :: Hscr_cp
1079 
1080 ! *************************************************************************
1081 
1082  !@hscr_t
1083  ! Integer values.
1084  Hscr_cp%id          = Hscr_in%id
1085  Hscr_cp%ikxc        = Hscr_in%ikxc
1086  Hscr_cp%inclvkb     = Hscr_in%inclvkb
1087  Hscr_cp%headform    = Hscr_in%headform
1088  Hscr_cp%fform       = Hscr_in%fform
1089  Hscr_cp%gwcalctyp   = Hscr_in%gwcalctyp
1090  Hscr_cp%nI          = Hscr_in%nI
1091  Hscr_cp%nJ          = Hscr_in%nJ
1092  Hscr_cp%nqibz       = Hscr_in%nqibz
1093  Hscr_cp%nqlwl       = Hscr_in%nqlwl
1094  Hscr_cp%nomega      = Hscr_in%nomega
1095  Hscr_cp%nbnds_used  = Hscr_in%nbnds_used
1096  Hscr_cp%npwe        = Hscr_in%npwe
1097  Hscr_cp%npwwfn_used = Hscr_in%npwwfn_used
1098  Hscr_cp%spmeth      = Hscr_in%spmeth
1099  Hscr_cp%test_type   = Hscr_in%test_type
1100  Hscr_cp%tordering   = Hscr_in%tordering
1101 
1102  ! Real variables
1103  Hscr_cp%mbpt_sciss = Hscr_in%mbpt_sciss
1104  Hscr_cp%spsmear  = Hscr_in%spsmear
1105  Hscr_cp%zcut     = Hscr_in%zcut
1106 
1107  ! Copy the abinit Header
1108  call hdr_copy(Hscr_in%Hdr,Hscr_cp%Hdr)
1109 
1110  Hscr_cp%titles(:) = Hscr_in%titles(:)
1111 
1112  ! Copy allocatable arrays.
1113  call alloc_copy(Hscr_in%gvec , Hscr_cp%gvec)
1114  call alloc_copy(Hscr_in%qibz , Hscr_cp%qibz)
1115  call alloc_copy(Hscr_in%qlwl , Hscr_cp%qlwl)
1116  call alloc_copy(Hscr_in%omega, Hscr_cp%omega)
1117 
1118 ! HSCR_NEW
1119  hscr_cp%awtr      =  hscr_in%awtr
1120  hscr_cp%icutcoul  =  hscr_in%icutcoul
1121  hscr_cp%vcutgeo   =  hscr_in%vcutgeo
1122  hscr_cp%gwcomp    =  hscr_in%gwcomp
1123  hscr_cp%gwgamma   =  hscr_in%gwgamma
1124  hscr_cp%gwencomp  =  hscr_in%gwencomp
1125  hscr_cp%kind_cdata  =  hscr_in%kind_cdata
1126 ! HSCR_NEW
1127 
1128 end subroutine hscr_copy

m_io_screening/hscr_free [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 hscr_free

FUNCTION

 Deallocate the components of the header structured datatype

INPUTS

 hdr <type(hdr_type)>=the header

OUTPUT

  (only deallocate)

SOURCE

1043 subroutine hscr_free(hscr)
1044 
1045 !Arguments ------------------------------------
1046  class(hscr_t),intent(inout) :: hscr
1047 
1048 ! *************************************************************************
1049 
1050  ABI_SFREE(hscr%gvec)
1051  ABI_SFREE(hscr%qibz)
1052  ABI_SFREE(hscr%qlwl)
1053  ABI_SFREE(hscr%omega)
1054 
1055  call hscr%Hdr%free()
1056 
1057 end subroutine hscr_free

m_io_screening/hscr_from_file [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

  hscr_from_file

FUNCTION

  Read the header of the (SCR/SUS) file

INPUTS

  path=File name
  comm = MPI communicator.

OUTPUT

  hscr<hscr_t>=The header.
  fform=Kind of the array in the file (0 signals an error)

SOURCE

309 subroutine hscr_from_file(hscr, path, fform, comm)
310 
311 !Arguments ------------------------------------
312 !scalars
313  class(hscr_t),intent(out) :: hscr
314  character(len=*),intent(in) :: path
315  integer,intent(in) :: comm
316  integer,intent(out) :: fform
317 
318 !Local variables-------------------------------
319 !scalars
320  integer,parameter :: rdwr5 = 5, master = 0
321  integer :: unt,my_rank,ierr
322  character(len=500) :: msg
323 
324 ! *************************************************************************
325 
326  my_rank = xmpi_comm_rank(comm)
327 
328  ! Master reads and broadcasts.
329  if (my_rank == master) then
330    if (.not. endswith(path, ".nc")) then
331      ! Fortran-IO
332      if (open_file(path,msg,newunit=unt,form="unformatted", status="old",action="read") /= 0) then
333        ABI_ERROR(msg)
334      end if
335      call hscr_io(hscr,fform,rdwr5,unt,xmpi_comm_self,master,IO_MODE_FORTRAN)
336      close(unt)
337    else
338      ! Netcdf format
339      NCF_CHECK(nctk_open_read(unt, path, xmpi_comm_self))
340      call hscr_io(hscr,fform,rdwr5,unt,xmpi_comm_self,master,IO_MODE_ETSF)
341      NCF_CHECK(nf90_close(unt))
342    end if
343 
344    ABI_CHECK(fform /= 0, sjoin("hscr_io returned fform == 0 while reading:", path))
345  end if
346 
347  ! Broadcast data.
348  if (xmpi_comm_size(comm) > 1) then
349    call hscr%bcast(master, my_rank, comm)
350    call xmpi_bcast(fform,master,comm,ierr)
351  end if
352 
353 end subroutine hscr_from_file

m_io_screening/hscr_io [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

  hscr_io

FUNCTION

 This subroutine deals with the I/O of the hscr_t structured variables (read/write/echo).
 According to the value of rdwr, it reads the header of a file, writes it, or echo the value
 of the structured variable to a file. Note that, when reading, different records of hscr_t
 are allocated here, according to the values of the read variables. Records of hscr_t should be
 deallocated correctly by a call to hdr_free when hscr_t is not used anymore.

INPUTS

  iomode=Option defining the file format of the external file.
  comm=MPI communicator.
  master=rank of the master node in comm, usually 0
  rdwr= if 1, read the hscr_t structured variable from the header of the file,
        if 2, write the header to unformatted file
        if 3, echo part of the header to formatted file (records 1 and 2)
        if 4, echo the header to formatted file
        if 5, read the hscr_t without rewinding (unformatted)
        if 6, write the hscr_t without rewinding (unformatted)
  unt=unit number of the file (unformatted if rdwr=1, 2, 5 or 6 formatted if rdwr=3,4)

OUTPUT

  (see side effects)

SIDE EFFECTS

  The following variables are both input or output :
  fform=kind of the array in the file
   if rdwr=1,5 : will be output ; if the reading fail, return fform=0
   if rdwr=2,3,4,6 : should be input, will be written or echo to file
  hscr_t <type(hscr_t)>=the header structured variable
   if rdwr=1,5 : will be output
   if rdwr=2,3,4,6 : should be input, will be written or echo to file

NOTES

 In all cases, the file is supposed to be open already
 When reading (rdwr=1) or writing (rdwr=2), rewind the file
 When echoing (rdwr=3) does not rewind the file.
 When reading (rdwr=5) or writing (rdwr=6), DOES NOT rewind the file

 In writing mode, the routine is supposed to called by the master node.
 no check is done, it is up to the developer.

SOURCE

402 subroutine hscr_io(hscr, fform, rdwr, unt, comm, master, iomode)
403 
404 !Arguments ------------------------------------
405 !scalars
406  class(hscr_t),target,intent(inout) :: hscr
407  integer,intent(inout) :: fform
408  integer,intent(in) :: rdwr,unt,iomode,comm,master
409 
410 !Local variables-------------------------------
411 !scalars
412  integer :: my_rank,nprocs,ncerr,ncid,varid,ierr !ii
413  character(len=500) :: errmsg
414  character(len=nctk_slen) :: varname !,head_shape,wing_shape
415 !arrays
416  real(dp),allocatable :: real_omega(:,:)
417  real(dp), ABI_CONTIGUOUS pointer :: r2vals(:,:) !,rvals3(:,:,:)
418 ! *************************************************************************
419 
420  DBG_ENTER("COLL")
421  !@hscr_t
422  ABI_UNUSED(master) ! FIXME
423 
424  ! Initialize MPI info for comm ===
425  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
426 
427  if (rdwr==1 .or. rdwr==5) then
428 
429    if (.True.) then
430    ! TODO: only master should read but then I have to skip the header.
431    !if (my_rank == master) then
432      ! Read the abinit header, rewinding of the file (if any) is done here.
433      if (iomode==IO_MODE_FORTRAN) then
434        call hdr_fort_read(hscr%hdr, unt, fform, rewind=(rdwr==1))
435      else if (iomode==IO_MODE_ETSF) then
436        call hdr_ncread(hscr%hdr, unt, fform)
437      end if
438 
439      ! Reset the variables absent in old versions.
440      Hscr%fform=fform
441 
442      if (iomode==IO_MODE_FORTRAN .or. iomode==IO_MODE_MPI) then
443        select case (fform)
444        case (1003, 1004)
445          ! File format for epsilon^-1, espilon, chi0
446          read(unt, err=10, iomsg=errmsg)hscr%titles
447          read(unt, err=10, iomsg=errmsg)&
448            hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp,&
449            hscr%nI, hscr%nJ, hscr%nqibz, hscr%nqlwl, hscr%nomega, hscr%nbnds_used,&
450            hscr%npwe, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering,&
451            hscr%awtr, hscr%icutcoul, hscr%gwgamma, hscr%vcutgeo(1:3)
452 
453          ! Read real scalars
454          read(unt, err=10, iomsg=errmsg)&
455            hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp,hscr%kind_cdata
456 
457          ! Allocate arrays and read them
458          call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
459 
460          read(unt, err=10, iomsg=errmsg)hscr%gvec(:,:)
461          read(unt, err=10, iomsg=errmsg)hscr%qibz(:,:)
462          read(unt, err=10, iomsg=errmsg)hscr%omega(:)
463 
464          ! Read data for q-->0 limit.
465          if (hscr%nqlwl>0) then
466            read(unt, err=10, iomsg=errmsg)hscr%qlwl(:,:)
467          end if
468 
469        case default
470          ABI_BUG(sjoin('Wrong fform read:', itoa(fform)))
471        end select
472 
473      else if (iomode == IO_MODE_ETSF) then
474        ncid = unt
475 
476        select case (fform)
477        case (1003, 1004)
478          ! Get dimensions and allocate arrays.
479          NCF_CHECK(nctk_get_dim(ncid, "number_of_coefficients_dielectric_function", hscr%npwe))
480          NCF_CHECK(nctk_get_dim(ncid, "number_of_qpoints_dielectric_function", hscr%nqibz))
481          NCF_CHECK(nctk_get_dim(ncid, "number_of_frequencies_dielectric_function", hscr%nomega))
482          NCF_CHECK(nctk_get_dim(ncid, "number_of_qpoints_gamma_limit", hscr%nqlwl))
483          NCF_CHECK(nctk_get_dim(ncid, "nI", hscr%ni))
484          NCF_CHECK(nctk_get_dim(ncid, "nJ", hscr%nj))
485          call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
486 
487          varid = nctk_idname(ncid, 'reduced_coordinates_plane_waves_dielectric_function')
488          NCF_CHECK(nf90_get_var(ncid, varid, hscr%gvec, start=[1,1,1]))
489          NCF_CHECK(nf90_get_var(ncid, vid('qpoints_dielectric_function'), hscr%qibz))
490 
491          ABI_MALLOC(real_omega, (2, hscr%nomega))
492          NCF_CHECK(nf90_get_var(ncid, vid('frequencies_dielectric_function'), real_omega))
493          hscr%omega = dcmplx(real_omega(1,:), real_omega(2,:))
494          ABI_FREE(real_omega)
495 
496          ! Read extra data added in new header.
497          NCF_CHECK(nf90_get_var(ncid, vid("vcutgeo"), hscr%vcutgeo))
498          NCF_CHECK(nf90_get_var(ncid, vid("id"), hscr%id))
499          NCF_CHECK(nf90_get_var(ncid, vid("ikxc"), hscr%ikxc))
500          NCF_CHECK(nf90_get_var(ncid, vid("inclvkb"), hscr%inclvkb))
501          NCF_CHECK(nf90_get_var(ncid, vid("headform"), hscr%headform))
502          NCF_CHECK(nf90_get_var(ncid, vid("fform"), hscr%fform))
503          NCF_CHECK(nf90_get_var(ncid, vid("gwcalctyp"), hscr%gwcalctyp))
504          NCF_CHECK(nf90_get_var(ncid, vid("nbands_used"), hscr%nbnds_used))
505          NCF_CHECK(nf90_get_var(ncid, vid("npwwfn_used"), hscr%npwwfn_used))
506          NCF_CHECK(nf90_get_var(ncid, vid("spmeth"), hscr%spmeth))
507          NCF_CHECK(nf90_get_var(ncid, vid("test_type"), hscr%test_type))
508          NCF_CHECK(nf90_get_var(ncid, vid("tordering"), hscr%tordering))
509          NCF_CHECK(nf90_get_var(ncid, vid("awtr"), hscr%awtr))
510          NCF_CHECK(nf90_get_var(ncid, vid("gw_icutcoul"), hscr%icutcoul))
511          NCF_CHECK(nf90_get_var(ncid, vid("gwcomp"), hscr%gwcomp))
512          NCF_CHECK(nf90_get_var(ncid, vid("gwgamma"), hscr%gwgamma))
513          NCF_CHECK(nf90_get_var(ncid, vid("mbpt_sciss"), hscr%mbpt_sciss))
514          NCF_CHECK(nf90_get_var(ncid, vid("spsmear"), hscr%spsmear))
515          NCF_CHECK(nf90_get_var(ncid, vid("zcut"), hscr%zcut))
516          NCF_CHECK(nf90_get_var(ncid, vid("gwencomp"), hscr%gwencomp))
517          NCF_CHECK(nf90_get_var(ncid, vid("kind_cdata"), hscr%kind_cdata))
518          call replace_ch0(hscr%kind_cdata)
519 
520          NCF_CHECK(nf90_get_var(ncid, vid("titles"), hscr%titles))
521          call replace_ch0(hscr%titles(:))
522 
523          ! TODO Read it
524          if (hscr%nqlwl /= 0) then
525            NCF_CHECK(nf90_get_var(ncid, vid("qpoints_gamma_limit"), hscr%qlwl))
526          end if
527 
528        case default
529          ABI_BUG(sjoin('Unsupported fform read:',itoa(fform)))
530        end select
531      else
532        ABI_ERROR(sjoin("Unsupported value of iomode:", iomode2str(iomode)))
533      end if
534 
535    end if ! master
536 
537    !call hscr%bcast(master, my_rank, comm)
538    !call hscr_mpio_skip(mpio_fh,fform,offset)
539 
540  else if (rdwr == 2 .or. rdwr == 6) then
541    ! Writing the header of an unformatted file.
542    ! Always use the latest version.
543 
544    if (iomode==IO_MODE_FORTRAN .or. iomode==IO_MODE_MPI) then
545      ! Write the abinit header.
546      call hscr%hdr%fort_write(unt, fform, ierr)
547      ABI_CHECK(ierr == 0, "hdr_fort_write retured ierr != 0")
548 
549      write(unt, err=10, iomsg=errmsg)hscr%titles
550 
551      ! Write integers
552      write(unt, err=10, iomsg=errmsg)&
553        hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp,&
554        hscr%nI, hscr%nJ, hscr%nqibz, hscr%nqlwl, hscr%nomega, hscr%nbnds_used,&
555        hscr%npwe, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering,&
556        hscr%awtr, hscr%icutcoul, hscr%gwgamma, hscr%vcutgeo(1:3)
557 
558      ! Write real scalars
559      write(unt, err=10, iomsg=errmsg)&
560        hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp,hscr%kind_cdata
561 
562      ! Write arrays
563      write(unt, err=10, iomsg=errmsg)hscr%gvec(:,:)
564      write(unt, err=10, iomsg=errmsg)hscr%qibz(:,:)
565      write(unt, err=10, iomsg=errmsg)hscr%omega(:)
566 
567      ! Add q-points for heads and wings for q-->0.
568      if (hscr%nqlwl > 0) then
569        write(unt, err=10, iomsg=errmsg)hscr%qlwl(:,:)
570      end if
571 
572    else if (iomode == IO_MODE_ETSF) then
573      ncid = unt
574      ! Write the abinit header, rewinding of the file (if any) is done here.
575      NCF_CHECK(hscr%hdr%ncwrite(ncid, fform, nc_define=.True.))
576 
577      ! Define dimensions
578      ! Part 2) of etsf-io specifications
579      ! FIXME: Spin is only used in particular cases, We usuall get the trace of W in spin space
580      ! and I'm not gonna allocate extra memory just to have up up, down down
581      ! Besides number_of_spins should be replaced by `number_of_spins_dielectric_function`
582      ! Should add spin_dependent attribute.
583      ncerr = nctk_def_dims(ncid, [ &
584        nctkdim_t("complex", 2), nctkdim_t("number_of_reduced_dimensions", 3), &
585        nctkdim_t("number_of_frequencies_dielectric_function", hscr%nomega), &
586        nctkdim_t("number_of_qpoints_dielectric_function", hscr%nqibz), &
587        nctkdim_t("number_of_qpoints_gamma_limit", hscr%nqlwl), &
588        nctkdim_t("number_of_spins", hscr%hdr%nsppol), &
589        nctkdim_t("nI", hscr%nI), nctkdim_t("nJ", hscr%nJ), &
590        nctkdim_t("number_of_coefficients_dielectric_function", hscr%npwe)], defmode=.True.)
591      NCF_CHECK(ncerr)
592 
593      ! Part 3) of the specs
594      ! (note that, in the specs, the Gs depend on the q-point but npwe is a scalar
595      ! basis_set is added by the abinit header.
596      ! FIXME: g-vectors are not written properly.
597      ncerr = nctk_def_arrays(ncid, [&
598        ! Standard
599        nctkarr_t('frequencies_dielectric_function', "dp", 'complex, number_of_frequencies_dielectric_function'), &
600        nctkarr_t('qpoints_dielectric_function', "dp", 'number_of_reduced_dimensions, number_of_qpoints_dielectric_function'),&
601        nctkarr_t('qpoints_gamma_limit', "dp", 'number_of_reduced_dimensions, number_of_qpoints_gamma_limit'), &
602        nctkarr_t('reduced_coordinates_plane_waves_dielectric_function', "i", &
603        'number_of_reduced_dimensions, number_of_coefficients_dielectric_function, number_of_qpoints_dielectric_function'), &
604        ! Abinit
605        nctkarr_t('vcutgeo', "dp", 'number_of_reduced_dimensions'), &
606        nctkarr_t('kind_cdata', "char", 'character_string_length'), &
607        nctkarr_t('titles', "char", 'character_string_length, two') &
608      ])
609      NCF_CHECK(ncerr)
610 
611      ! FIXME problem with q-points, heads and wings?
612      ! The order in P in the specs is wrong, q should be the last dimension here I use the "correct" version
613      ! TODO: Remove. Use abifile_t
614      varname = ncname_from_id(hscr%id)
615      ncerr = nctk_def_arrays(ncid, &
616        nctkarr_t(varname, "dp", &
617 &"complex, number_of_coefficients_dielectric_function, number_of_coefficients_dielectric_function,&
618 &number_of_spins, number_of_spins, number_of_frequencies_dielectric_function, number_of_qpoints_dielectric_function"))
619      NCF_CHECK(ncerr)
620 
621      !write(std_out,*)"nqlwl",hscr%nqlwl
622      NCF_CHECK(nctk_set_datamode(ncid))
623      call c_f_pointer(c_loc(hscr%omega(1)), r2vals, shape=[2, size(hscr%omega)])
624      NCF_CHECK(nf90_put_var(ncid, vid('frequencies_dielectric_function'), r2vals))
625      NCF_CHECK(nf90_put_var(ncid, vid('qpoints_dielectric_function'), hscr%qibz))
626      NCF_CHECK(nf90_put_var(ncid, vid('reduced_coordinates_plane_waves_dielectric_function'), hscr%gvec))
627 
628      NCF_CHECK(nf90_put_var(ncid, vid("titles"), hscr%titles))
629      NCF_CHECK(nf90_put_var(ncid, vid("kind_cdata"), hscr%kind_cdata))
630      NCF_CHECK(nf90_put_var(ncid, vid("vcutgeo"), hscr%vcutgeo))
631 
632      ncerr = nctk_defnwrite_ivars(ncid, [character(len=nctk_slen) :: &
633        "id", "ikxc", "inclvkb", "headform", "fform", "gwcalctyp", &
634        "nbands_used", "npwwfn_used", "spmeth", "test_type", "tordering", "awtr", "gw_icutcoul", &
635        "gwcomp", "gwgamma" &
636       ],&
637       [ hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp, &
638        hscr%nbnds_used, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering, hscr%awtr, hscr%icutcoul, &
639        hscr%gwcomp, hscr%gwgamma &
640      ])
641      NCF_CHECK(ncerr)
642 
643      ncerr = nctk_defnwrite_dpvars(ncid, [character(len=nctk_slen) :: &
644         "mbpt_sciss", "spsmear", "zcut", "gwencomp"], &
645         [hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp &
646      ])
647      NCF_CHECK(ncerr)
648 
649      ! Add q-points for heads and wings for q-->0.
650      if (hscr%nqlwl > 0) then
651        !  MG: This part has been commented out as it's not used
652        !  head_shape = "complex, number_of_spins, number_of_spins, number_of_frequencies_dielectric_function"
653        !  head_shape = trim(head_shape)//", number_of_qpoints_gamma_limit"
654 
655        !  wing_shape = "complex, number_of_coefficients_dielectric_function, number_of_spins, number_of_spins"
656        !  wing_shape = trim(wing_shape)//", number_of_frequencies_dielectric_function, number_of_qpoints_gamma_limit"
657 
658        !  ncerr = nctk_def_arrays(ncid, [&
659        !    nctkarr_t("dielectric_function_head", "dp", head_shape),&
660        !    nctkarr_t("dielectric_function_upper_wing", "dp", wing_shape),&
661        !    nctkarr_t("dielectric_function_lower_wing", "dp", wing_shape)], defmode=.True.)
662        !  NCF_CHECK(ncerr)
663 
664        NCF_CHECK(nctk_set_datamode(ncid))
665        NCF_CHECK(nf90_put_var(ncid, vid('qpoints_gamma_limit'), hscr%qlwl))
666      end if
667    else
668      ABI_ERROR(sjoin('Unsupported iomode:',iomode2str(iomode)))
669    end if
670 
671  else
672    ABI_BUG(sjoin("Wrong value for rdwr:", itoa(rdwr)))
673  end if ! read/write/echo
674 
675  DBG_EXIT("COLL")
676 
677  return
678 
679  ! Handle Fortran IO error
680  10 continue
681  ABI_ERROR(errmsg)
682 
683 contains
684  integer function vid(vname)
685    character(len=*),intent(in) :: vname
686    vid = nctk_idname(ncid, vname)
687  end function vid
688 
689 end subroutine hscr_io

m_io_screening/hscr_malloc [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 hscr_malloc

FUNCTION

 Allocate the components of the header structured datatype except for hscr%hdr

SOURCE

1008 subroutine hscr_malloc(hscr, npwe, nqibz, nomega, nqlwl)
1009 
1010 !Arguments ------------------------------------
1011 !scalars
1012  class(hscr_t),intent(inout) :: Hscr
1013  integer,intent(in) :: npwe, nqibz, nomega, nqlwl
1014 
1015 ! *************************************************************************
1016 
1017  !@hscr_t
1018  ABI_MALLOC(hscr%gvec, (3, npwe))
1019  ABI_MALLOC(hscr%qibz, (3, nqibz))
1020  ABI_MALLOC(hscr%qlwl, (3, nqlwl))
1021  ABI_MALLOC(hscr%omega, (nomega))
1022 
1023 end subroutine hscr_malloc

m_io_screening/hscr_merge [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 hscr_merge

FUNCTION

 This subroutine merges diffrent header structured variable (hscr_t)

INPUTS

  Hscr_in(:) <hscr_t)>=List of headers to be merged.

OUTPUT

  Hscr_out<hscr_t>=The output merged header.

SOURCE

1148 subroutine hscr_merge(Hscr_in, Hscr_out)
1149 
1150 !Arguments ------------------------------------
1151 !scalars
1152  type(hscr_t),intent(in) :: Hscr_in(:)
1153  type(hscr_t),intent(out) :: Hscr_out
1154 
1155 !Local variables-------------------------------
1156 !scalars
1157  integer :: nhds,restart,restartpaw,ihd,ii,nqtot,nqneq
1158  logical :: isok
1159  character(len=500) :: msg
1160 !arrays
1161  real(dp),allocatable :: qset(:,:)
1162 ! *************************************************************************
1163 
1164  !@hscr_t
1165  nhds=SIZE(Hscr_in)
1166 
1167  ! Initial copy of the header ===
1168  ! If multiple headers, select the header containing q-->0 so that we copy also heads and wings
1169  ii = imax_loc(Hscr_in(:)%nqlwl)
1170  call hscr_copy(Hscr_in(ii),Hscr_out)
1171  if (nhds==1) return
1172 
1173  ! Check consistency of the abinit Headers.
1174  ! FFT grid might be q-point dependent so we stop only when restart==0
1175  isok=.TRUE.
1176  do ihd=2,nhds
1177    call hdr_check(Hscr_in(1)%fform,Hscr_in(ihd)%fform,Hscr_in(1)%Hdr,Hscr_in(ihd)%Hdr,'COLL',restart,restartpaw)
1178    if (restart==0) then
1179      isok=.FALSE.
1180      write(msg,'(a,i0,a)')' Abinit header no.',ihd,' is not consistent with the first header '
1181      ABI_WARNING(msg)
1182    end if
1183  end do
1184  if (.not.isok) then
1185    ABI_ERROR('Cannot continue, Check headers')
1186  end if
1187 
1188  ! Now check variables related to polarizability|epsilon^{-1}.
1189  ! 1) Tests quantities that must be equal
1190  ii = assert_eq(Hscr_in(:)%ID,       'Headers have different Identifiers')
1191  ii = assert_eq(Hscr_in(:)%ikxc,     'Headers have different ikxc'       )
1192  ii = assert_eq(Hscr_in(:)%headform, 'Headers have different headform'   )
1193  ii = assert_eq(Hscr_in(:)%fform,    'Headers have different fform'      )
1194  ii = assert_eq(Hscr_in(:)%gwcalctyp,'Headers have different gwcalctyp'  )
1195  ii = assert_eq(Hscr_in(:)%nI,       'Headers have different nI'         )
1196  ii = assert_eq(Hscr_in(:)%nJ,       'Headers have different nJ'         )
1197  ii = assert_eq(Hscr_in(:)%nomega,   'Headers have different nomega'     )
1198  ii = assert_eq(Hscr_in(:)%test_type,'Headers have different test_type'  )
1199  ii = assert_eq(Hscr_in(:)%tordering,'Headers have different tordering'  )
1200 
1201  ! This is not mandatory but makes life easier!
1202  ii = assert_eq(Hscr_in(:)%npwe,'Headers have different number of G-vectors'  )
1203 
1204  do ihd=2,nhds
1205    if (ANY(ABS(Hscr_in(ihd)%omega-Hscr_in(1)%omega)>tol6)) then
1206      write(msg,'(a,i0,a)')' Frequencies in the first and the ',ihd,'-th header differ'
1207      ABI_ERROR(msg)
1208    end if
1209    if (ANY(Hscr_in(ihd)%gvec(:,:)-Hscr_in(1)%gvec(:,:)/=0)) then
1210      write(msg,'(a,i0,a)')' Incompatible G-vector list found in the ',ihd,'-th header'
1211      ABI_ERROR(msg)
1212    end if
1213    if (hscr_in(ihd)%kind_cdata /= hscr_in(1)%kind_cdata) then
1214      write(msg,'(3a,i0,2a)')' Files contain data with different precisions.',ch10,&
1215      "In particular the ",ihd,'-th header has precision:',trim(hscr_in(ihd)%kind_cdata)
1216      ABI_ERROR(msg)
1217    end if
1218  end do !ihd
1219 
1220  ! If error is not fatal, just warn ===
1221  if (ANY(Hscr_in(:)%npwwfn_used/=Hscr_in(1)%npwwfn_used)) then
1222    ABI_COMMENT('Files have been produced with a different number of planewaves for the wavefunctions.')
1223  end if
1224  if (ANY(Hscr_in(:)%nbnds_used/=Hscr_in(1)%nbnds_used)) then
1225    ABI_COMMENT('Files have been produced with a different number of bands.')
1226  end if
1227  if (ANY(Hscr_in(:)%spmeth/=Hscr_in(1)%spmeth)) then
1228    ABI_COMMENT('Files have been produced with different algorithms.')
1229  end if
1230  if (ANY(ABS(Hscr_in(:)%mbpt_sciss-Hscr_in(1)%mbpt_sciss)>tol6)) then
1231    ABI_COMMENT('Files have benn produced with different values of mbpt_sciss.')
1232  end if
1233  if (ANY(ABS(Hscr_in(:)%spsmear-Hscr_in(1)%spsmear)>tol6)) then
1234    ABI_COMMENT('Files have been produced with different values of spsmear.')
1235  end if
1236  if (ANY(ABS(Hscr_in(:)%zcut-Hscr_in(1)%zcut)>tol6)) then
1237    ABI_COMMENT('Files have been produced with different values of zcut.')
1238  end if
1239 
1240  ! Now merge the list of q-points.
1241  ! Take the union of the q-points, remove possible duplicated
1242  ! are change the parameters in hscr_out that depends on q-points.
1243  nqtot=SUM(Hscr_in(:)%nqibz)
1244  ABI_MALLOC(qset,(3,nqtot))
1245 
1246  ii=0
1247  do ihd=1,nhds
1248    qset(:,ii+1:ii+Hscr_in(ihd)%nqibz)=Hscr_in(ihd)%qibz(:,:)
1249    ii=ii+Hscr_in(ihd)%nqibz
1250  end do
1251 
1252  call remove_copies(nqtot,qset,nqneq,isequalk)
1253 
1254  if (nqneq /= nqtot) then
1255    write(msg,'(3a,2(i0,a))')&
1256     'COMMENT: Headers contain duplicated q-points ',ch10,&
1257     'Found ',nqneq,' distinct q-points among the total ',nqtot,' points reported in the headers. '
1258    call wrtout(std_out, msg)
1259  end if
1260 
1261  Hscr_out%nqibz = nqneq
1262  ABI_FREE(Hscr_out%qibz)
1263  ABI_MALLOC(Hscr_out%qibz,(3,nqneq))
1264  Hscr_out%qibz(:,:)=qset(:,1:nqneq)
1265  ABI_FREE(qset)
1266 
1267 end subroutine hscr_merge

m_io_screening/hscr_mpio_skip [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

  hscr_mpio_skip

FUNCTION

   Skip the header of the (SCR|SUSC) file in MPI-IO mode. This routine uses local MPI-IO calls hence
   it can be safely called by master node only. Note however that in this case the
   offset has to be communicated to the other nodes.

INPUTS

  mpio_fh=MPI-IO file handler
  fmarker_bsize   = Byte length of Fortran record marker.
  fmarker_mpi_type= MPI type of the Fortran record marker

OUTPUT

  fform=kind of the array in the file
  offset=The offset of the Fortran record located immediately below the Abinit header.

SOURCE

1710 subroutine hscr_mpio_skip(mpio_fh, fform, offset)
1711 
1712 !Arguments ------------------------------------
1713  integer,intent(in) :: mpio_fh
1714  integer,intent(out) :: fform
1715  integer(kind=XMPI_OFFSET_KIND),intent(out) :: offset
1716 
1717 !Local variables-------------------------------
1718 !scalars
1719  integer :: bsize_frm,mpi_type_frm
1720 #ifdef HAVE_MPI_IO
1721  integer :: ierr,isk,nqlwl
1722  !character(len=500) :: msg
1723 !arrays
1724  integer(kind=MPI_OFFSET_KIND) :: fmarker,positloc
1725  integer :: statux(MPI_STATUS_SIZE)
1726 #endif
1727 
1728 ! *************************************************************************
1729 
1730  offset = 0
1731  bsize_frm    = xmpio_bsize_frm    ! Byte size of the Fortran record marker.
1732  mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker.
1733 
1734  call hdr_mpio_skip(mpio_fh,fform,offset)
1735 
1736  call wrtout(std_out, sjoin("in hdr_mpio_skip with fform = ",itoa(fform)))
1737 
1738 #ifdef HAVE_MPI_IO
1739  select case (fform)
1740  case (1003, 1004)
1741    ! Ship the titles
1742    call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1743 
1744    ! read nqlwl from the 2d record.
1745    positloc  = offset + bsize_frm + 9*xmpi_bsize_int
1746    call MPI_FILE_READ_AT(mpio_fh,positloc,nqlwl,1,MPI_INTEGER,statux,ierr)
1747    call wrtout(std_out, sjoin("nqlwl = ",itoa(nqlwl)))
1748 
1749    do isk=1,5
1750      call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1751    end do
1752 
1753    if (nqlwl>0) then  ! skip qlwl
1754      call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1755    end if
1756 
1757  case default
1758    ABI_BUG(sjoin('Wrong fform read:', itoa(fform)))
1759  end select
1760 
1761 #else
1762  ABI_ERROR("hscr_mpio_skip cannot be used when MPI-IO is not enabled")
1763 #endif
1764 
1765 end subroutine hscr_mpio_skip

m_io_screening/hscr_new [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

  hscr_new

FUNCTION

  Initialize the Hscr datatype and most of its content from the em1params_t data type Ep.

INPUTS

  varname=Name of the netcdf variable (used to get fform and ID).
  ikxc=Integer flag definining the type of XC kernel (0 if None i.e RPA)
  test_type=Integer flag defining the type of probing charge (0 for None)
  tordering=The time-ordering of the Response function.
  gvec(3,Ep%npwe)=The G-vectors used.
  Ep<em1params_t>=Parameters defining the calculation of the screening.
  hdr_abinit<hdr_type>=The abinit header.

OUTPUT

  Hscr<type(hscr_t)>=the header, initialized.

SOURCE

823 type(hscr_t) function hscr_new(varname,dtset,ep,hdr_abinit,ikxc,test_type,tordering,titles,ngvec,gvec) result(hscr)
824 
825 !Arguments ------------------------------------
826 !scalars
827  integer,intent(in) :: ikxc,test_type,tordering,ngvec
828  character(len=*),intent(in) :: varname
829  type(dataset_type),intent(in) :: dtset
830  type(em1params_t),intent(in) :: Ep
831  type(hdr_type),intent(in) :: hdr_abinit
832 !arrays
833  integer,intent(in) :: gvec(3,ngvec)
834  character(len=80),intent(in) :: titles(2)
835 
836 !Local variables-------------------------------
837  integer :: id
838  type(abifile_t) :: abifile
839 ! *************************************************************************
840 
841  !@hscr_t
842  ABI_CHECK(ngvec == Ep%npwe, 'ngvec/=Ep%npwe')
843 
844  ! Identifier used to define the type of Response function (e^-1, chi0)
845  id = 0
846  if (varname == "polarizability") id = 1
847  if (varname == "inverse_dielectric_function") id = 4
848  ABI_CHECK(id /= 0, sjoin("Invalid varname: ",varname))
849 
850  ! Get fform from abifile.
851  abifile = abifile_from_varname(varname)
852  if (abifile%fform == 0) then
853     ABI_ERROR(sjoin("Cannot find any abifile object associated to varname:", varname))
854  end if
855 
856  ! Copy the abinit header.
857  call hdr_copy(hdr_abinit, Hscr%Hdr)
858 
859  ! Initialize quantities related to the screening file
860  hscr%id         =id
861  hscr%ikxc       =ikxc
862  hscr%inclvkb    =Ep%inclvkb
863  hscr%headform   =HSCR_LATEST_HEADFORM
864  hscr%fform      =abifile%fform
865  hscr%gwcalctyp  =Ep%gwcalctyp
866  hscr%nI         =Ep%nI
867  hscr%nJ         =Ep%nJ
868  hscr%nqibz      =Ep%nqcalc  ! nqcalc == nqibz except if we split the calculation with nqptdm
869  hscr%nqlwl      =Ep%nqlwl
870  hscr%nomega     =Ep%nomega
871  hscr%nbnds_used =Ep%nbnds
872  hscr%npwe       =Ep%npwe
873  hscr%npwwfn_used=Ep%npwwfn
874  hscr%spmeth     =Ep%spmeth
875  hscr%test_type  =test_type
876  hscr%tordering  =tordering
877  hscr%mbpt_sciss =Ep%mbpt_sciss
878  hscr%spsmear    =Ep%spsmear
879  hscr%zcut       =Ep%zcut
880 
881  hscr%titles(:)=titles(:)
882 
883  call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
884  hscr%gvec(:,:) = gvec(1:3,1:Ep%npwe)
885  hscr%qibz(:,:) = Ep%qcalc
886  hscr%qlwl(:,:) = Ep%qlwl
887  hscr%omega(:) = Ep%omega
888 
889 ! HSCR_NEW
890  hscr%awtr = dtset%awtr
891  hscr%icutcoul = dtset%gw_icutcoul
892  hscr%vcutgeo = dtset%vcutgeo
893  hscr%gwcomp = dtset%gwcomp
894  hscr%gwgamma = dtset%gwgamma
895  hscr%gwencomp = dtset%gwencomp
896  hscr%kind_cdata = "dpc" ! For the time being, data is always written in double-precision
897 ! HSCR_NEW
898 
899 end function hscr_new

m_io_screening/hscr_print [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 hscr_print

FUNCTION

  Prints info on the header of the SCR|SUSC file.

INPUTS

OUTPUT

SOURCE

707 subroutine hscr_print(Hscr, header, unit, prtvol, mode_paral)
708 
709 !Arguments ------------------------------------
710 !scalars
711  class(hscr_t),intent(in) :: hscr
712  integer,intent(in),optional :: prtvol,unit
713  character(len=4),intent(in),optional :: mode_paral
714  character(len=*),intent(in),optional :: header
715 
716 !Local variables-------------------------------
717 !scalars
718  integer :: iomega,iqibz,unt,verbose
719  character(len=4) :: mode
720  character(len=500) :: msg
721 
722 ! *************************************************************************
723 
724  unt=std_out; if (PRESENT(unit      )) unt    =unit
725  verbose=0  ; if (PRESENT(prtvol    )) verbose=prtvol
726  mode='COLL'; if (PRESENT(mode_paral)) mode   =mode_paral
727 
728  if (PRESENT(header)) then
729    msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
730    call wrtout(unt,msg,mode)
731  end if
732 
733  write(msg,'(1x,a)')TRIM(hscr%titles(1))
734  call wrtout(unt,msg,mode)
735  write(msg,'(1x,a)')TRIM(hscr%titles(2))
736  call wrtout(unt,msg,mode)
737  write(msg,'(a,i8)') ' Identifier                ',hscr%ID
738  call wrtout(unt,msg,mode)
739  write(msg,'(a,i8)') ' Kxc kernel                ',hscr%ikxc
740  call wrtout(unt,msg,mode)
741  write(msg,'(a,i8)') ' Treatment of q-->0 limit  ',hscr%inclvkb
742  call wrtout(unt,msg,mode)
743  write(msg,'(a,i8)') ' headform                  ',hscr%headform
744  call wrtout(unt,msg,mode)
745  write(msg,'(a,i8)') ' fform                     ',hscr%fform
746  call wrtout(unt,msg,mode)
747  write(msg,'(a,i8)') ' gwcalctyp                 ',hscr%gwcalctyp
748  call wrtout(unt,msg,mode)
749  write(msg,'(a,2i8)')' Number of components      ',hscr%nI,hscr%nJ
750  call wrtout(unt,msg,mode)
751  write(msg,'(a,i8)') ' Number of q-points        ',hscr%nqibz
752  call wrtout(unt,msg,mode)
753  write(msg,'(a,i8)') ' Number of q-directions    ',hscr%nqlwl
754  call wrtout(unt,msg,mode)
755  write(msg,'(a,i8)') ' Number of frequencies     ',hscr%nomega
756  call wrtout(unt,msg,mode)
757  write(msg,'(a,i8)') ' Number of bands used      ',hscr%nbnds_used
758  call wrtout(unt,msg,mode)
759  write(msg,'(a,i8)') ' Dimension of matrix       ',hscr%npwe
760  call wrtout(unt,msg,mode)
761  write(msg,'(a,i8)') ' Number of planewaves used ',hscr%npwwfn_used
762  call wrtout(unt,msg,mode)
763  write(msg,'(a,i8)') ' Spectral method           ',hscr%spmeth
764  call wrtout(unt,msg,mode)
765  write(msg,'(a,i8)') ' Test_type                 ',hscr%test_type
766  call wrtout(unt,msg,mode)
767  write(msg,'(a,i8)') ' Time-ordering             ',hscr%tordering
768  call wrtout(unt,msg,mode)
769  write(msg,'(a,es16.6)')' Scissor Energy             ',hscr%mbpt_sciss
770  call wrtout(unt,msg,mode)
771  write(msg,'(a,es16.6)')' Spectral smearing          ',hscr%spsmear
772  call wrtout(unt,msg,mode)
773  write(msg,'(a,es16.6)')' Complex Imaginary Shift    ',hscr%zcut
774  call wrtout(unt,msg,mode)
775 
776  if (verbose==0) then
777    call wrtout(unt,' The header contains additional records.',mode)
778  else
779    write(msg,'(2a)')ch10,' q-points [r.l.u.]:'
780    call wrtout(unt,msg,mode)
781    do iqibz=1,hscr%nqibz
782      write(msg,'(i5,3f12.6)')iqibz,hscr%qibz(:,iqibz)
783      call wrtout(unt,msg,mode)
784    end do
785 
786    write(msg,'(2a)')ch10,' Frequencies used [eV]:'
787    call wrtout(unt,msg,mode)
788    do iomega=1,hscr%nomega
789      write(msg,'(i3,2f7.2)')iomega,REAL(hscr%omega(iomega))*Ha_eV,AIMAG(hscr%omega(iomega))*Ha_eV
790      call wrtout(unt,msg,mode)
791    end do
792  end if
793 
794  ! Echo the abinit header.
795  !if (prtvol>0) call hdr_echo(hscr%hdr,fform,rdwr,unit=unt)
796 
797 end subroutine hscr_print

m_io_screening/hscr_t [ Types ]

[ Top ] [ m_io_screening ] [ Types ]

NAME

  hscr_t

FUNCTION

  The structure defining the header of the SCR/SUSC file.
  hscr_t contains the most important dimensions associated to the SCR/SUSC matrix,
  important GW metadata and the Abint header. The SCR/SUS matrices are saved
  with the same format. There are nqibz blocks, each block contains (npwe,npwe,nomega) matrices.
  SCR and SUS files mainly differ for what concerns the treatment of the q-->0 limit.
  The treatment of the non-analytic behaviour is not yet implemented but the main ideas are
  sketched in the NOTES below.

NOTES

  On the treatment of the q-->0 limit

  1) q=Gamma should be the first q-point

  2) This point contains an initial section with data used to treat the q-->0 limit, followed
     by the SCR/SUS matrix evaluated for the small q-point (qlwl). The header should contains
     enough info so that we can skip this section and use the pre-existing routines to read
     the matrices.

  3) The data stored in the q-->0 section depends on the type of matrix stored in the file:

     SUS file: we store the tensor and the wings needed to reconstruct the q-dependence around Gamma.
       This section if followed by the X(G,G') matrix evaluated at qlwl.
       Note that the content of the SUS file does not depend on a possible cutoff in vcoul.

     SCR file: two cases must be considered:

     No cutoff in vcoul:
        The q-->0 section stores the tensor, the wings as well as the inverse of the body B^{-1}.
        This info is used to integrate W(q) around q==Gamma.

     cutoff in vcoul:
        The content of the q-->0 section depends on dimensionality of the system.
        In 2D we need info on chi0 as well as A and a. See http://arxiv.org/pdf/1511.00129v1.pdf
        I think that this kind of calculations are easy to implement if we start from the SUS file.
        Ok, we have to recompute e-1 at each run but the logic is easier to implement.

SOURCE

105  type,public :: hscr_t
106 
107   integer :: id = -1
108     ! Matrix identifier: 1 for chi0, 2 for chi, 3 for epsilon, 4 for espilon^{-1}
109 
110   integer :: ikxc = 0
111     ! Kxc kernel used,
112     ! 0 for None (RPA), >0 for static TDDFT (=ixc), <0 for frequency-dependent TDDFT
113 
114   integer :: inclvkb = 2
115     ! q-->0 treatment, 0 for None, 1-2 for transversal gauge, 3 for longitudinal
116 
117   integer :: headform
118     ! format of the SCR header
119 
120   integer :: fform
121     ! File format
122 
123   integer :: gwcalctyp = 0
124     ! Calculation type (G0W0, G0W, GW ...)
125 
126   integer :: nI = 1, nJ = 1
127     ! Number of spin components (rows,columns) in chi|eps^-1. (1,1) if collinear.
128     ! The internal representation of the matrix is eps(nI*npwe,nJ*npwe)
129 
130   integer :: nqibz = -1
131     ! Number of q-points in the IBZ.
132 
133   integer :: nqlwl = 1
134     ! Number of points for the treatment of the long wavelength limit.
135 
136   integer :: nomega = -1
137     ! Total number of frequencies.
138 
139   integer :: nbnds_used = -1
140     ! Number of bands used during the screening calculation (only for info)
141 
142   integer :: npwe = -1
143     ! Number of G vectors reported on the file.
144 
145   integer :: npwwfn_used = -1
146     ! Number of G vectors for wavefunctions used during the screening calculation (only for info)
147 
148   integer :: spmeth = 0
149     ! Method used to approximate the delta function in the expression for Im Chi_0
150 
151   integer :: test_type
152     ! 1 for TEST-PARTICLE, 2 for TEST-ELECTRON.
153 
154   integer :: tordering = 1
155     ! 1 for Time-Ordered, 2 for Advanced, 3 for Retarded.
156 
157 ! HSCR_NEW
158   integer :: awtr = 1
159   ! Input variable (time-reversal symmetry in RPA expression)
160 
161   integer :: icutcoul = 0
162   ! Input variable (Coulomb singularity treatment)
163 
164   integer :: gwcomp = 0
165   ! Input variable (GW compensation energy technique)
166 
167   integer :: gwgamma = 0
168   ! Input variable Vertex correction
169 ! HSCR_NEW
170 
171   real(dp) :: mbpt_sciss = zero
172     ! Scissor Energy, zero if not used
173 
174   real(dp) :: spsmear = zero
175     ! Smearing of the delta in case of spmeth==2
176 
177   real(dp) :: zcut = -one
178     ! Imaginary shift to avoid the poles along the real axis.
179 
180 ! HSCR_NEW
181   real(dp) :: gwencomp = -one
182    ! Input variable (GW compensation energy technique)
183 
184   character(len=3) :: kind_cdata
185   ! Flag to signal whether the data is in single or double precision ("spc" or "dpc")
186   ! For the time being, we always write/read in double precision.
187   ! This flag could be use to reduce the memory requirements if spc:
188   ! we run calculations in single precision dump the results with the same precision without
189   ! having to allocate extra memory.
190 ! HSCR_NEW
191 
192 !arrays
193 
194 ! HSCR_NEW
195   real(dp) :: vcutgeo(3) = zero
196    ! Input variable (defines coulomb cutoff)
197 ! HSCR_NEW
198 
199   character(len=80) :: titles(2)
200     ! Titles describing the content of the file.
201 
202   integer,allocatable  :: gvec(:,:)
203     ! gvec(3,npwe)
204     ! G vectors in reduced coordinates.
205 
206   real(dp),allocatable :: qibz(:,:)
207     ! qibz(3,nqibz)
208     ! q-points in the IBZ in reduced coordinates.
209 
210   real(dp),allocatable :: qlwl(:,:)
211     ! qlwl(3,nqlwl)
212     ! q-points for the long wave-length limit treatment (r.l.u)
213 
214   complex(dpc),allocatable :: omega(:)
215     ! omega(nomega)
216     ! All frequencies calculated both along the real and the imaginary axis.
217     ! Real frequencies are packed in the first section.
218     ! TODO: Add frequency mesh type?
219 
220   type(hdr_type) :: hdr
221     ! The abinit header.
222 
223   contains
224 
225     procedure :: from_file => hscr_from_file    ! Read the header from file.
226     procedure :: print => hscr_print            ! Print the SCR-related part of the header.
227     procedure :: bcast => hscr_bcast            ! Broadcast the header.
228     procedure :: free => hscr_free              ! Free the header.
229     procedure :: copy => hscr_copy              ! Copy the SCR|SUSC header.
230 
231  end type hscr_t

m_io_screening/ioscr_qmerge [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 ioscr_qmerge

FUNCTION

  Produce new file by merging the q-points stored in other files.
  This routine should be called by a single MPI process.

INPUTS

  nfiles=Number of files to be merged.
  filenames(nfiles)=Paths of files to be merged.
  hscr_files(nfiles)<hscr_t>=Heades of the files to be merged.
  fname_out=Name of the file to be produced.

OUTPUT

  ohscr<hscr_t>=The header of the output file.

SOURCE

1789 subroutine ioscr_qmerge(nfiles, filenames, hscr_files, fname_out, ohscr)
1790 
1791 !Arguments ------------------------------------
1792 !scalars
1793  integer,intent(in) :: nfiles
1794  character(len=*),intent(in) :: fname_out
1795  type(hscr_t),intent(out) :: ohscr
1796 !arrays
1797  character(len=*),intent(in) :: filenames(nfiles)
1798  type(hscr_t),intent(in) :: hscr_files(nfiles)
1799 
1800 !Local variables-------------------------------
1801 !scalars
1802  integer,parameter :: rdwr2=2,master=0
1803  integer :: iqibz,ifound,ifile,iqf,ount,iomode,fform_merge,comm,nomega4m,npwe4m,iqiA,ierr
1804  character(len=500) :: msg
1805  character(len=nctk_slen) :: varname
1806  type(abifile_t) :: abifile
1807 !arrays
1808  integer,allocatable :: merge_table(:,:)
1809  real(dp) :: qdiff(3)
1810  complex(gwpc),allocatable :: epsm1(:,:,:,:)
1811 
1812 ! *************************************************************************
1813 
1814  comm = xmpi_comm_self
1815 
1816  if (file_exists(fname_out)) then
1817    ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
1818  end if
1819 
1820  ! Merge the headers creating the full list of q-points.
1821  call hscr_merge(Hscr_files(1:nfiles), ohscr)
1822  call hscr_print(ohscr, header='Header of the final file', unit=std_out, prtvol=1)
1823 
1824  ! For each q to be merged, save the index of the file where q is stored as well as its sequential index.
1825  ! Useful to do the merge point-by-point thus avoiding the allocation of the entire epsm1 array.
1826  ABI_MALLOC(merge_table,(ohscr%nqibz,2))
1827  do iqibz=1,ohscr%nqibz
1828    ifound=0
1829    fl: do ifile=1,nfiles
1830      do iqf=1,Hscr_files(ifile)%nqibz
1831        qdiff(:)=ohscr%qibz(:,iqibz)-Hscr_files(ifile)%qibz(:,iqf)
1832        if (all(abs(qdiff) < GW_TOLQ)) then
1833          merge_table(iqibz,1)=ifile
1834          merge_table(iqibz,2)=iqf
1835          ifound=ifound+1
1836          write(msg,'(a,3f12.6,2a)')'. q-point:',ohscr%qibz(:,iqibz),' will be taken from ',TRIM(filenames(ifile))
1837          call wrtout(std_out, msg)
1838          EXIT fl
1839        end if
1840      end do
1841    end do fl
1842    ! Check if q-point has been found, multiple q-points not allowed.
1843    ABI_CHECK(ifound == 1, 'ifound/=1')
1844  end do
1845 
1846  iomode = IO_MODE_FORTRAN; if (endswith(fname_out, ".nc")) iomode = IO_MODE_ETSF
1847  if (iomode == IO_MODE_FORTRAN) then
1848    if (open_file(fname_out,msg,newunit=ount,status='new',form='unformatted') /= 0) then
1849      ABI_ERROR(msg)
1850    end if
1851  else
1852    NCF_CHECK(nctk_open_create(ount, fname_out, comm))
1853  end if
1854 
1855  ! Write the header.
1856  fform_merge = hscr_files(1)%fform
1857  abifile = abifile_from_fform(fform_merge)
1858  if (abifile%fform == 0) then
1859    ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
1860  end if
1861  varname = abifile%varname
1862 
1863  if (any(hscr_files(:)%fform /= hscr_files(1)%fform)) then
1864    write(std_out,*)"fforms: ",hscr_files(:)%fform
1865    ABI_ERROR("Files to be merged have different fform. Cannot merge data")
1866  end if
1867 
1868  call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode)
1869 
1870  npwe4m   = ohscr%npwe
1871  nomega4m = ohscr%nomega
1872 
1873  ABI_MALLOC_OR_DIE(epsm1,(npwe4m,npwe4m,nomega4m,1), ierr)
1874 
1875  do iqibz=1,ohscr%nqibz
1876    ifile = merge_table(iqibz,1)
1877    iqiA  = merge_table(iqibz,2)
1878    call read_screening(varname,filenames(ifile),npwe4m,1,nomega4m,epsm1,iomode,comm,iqiA=iqiA)
1879    call write_screening(varname,ount,iomode,npwe4m,nomega4m,iqibz,epsm1)
1880  end do
1881 
1882  ABI_FREE(epsm1)
1883  ABI_FREE(merge_table)
1884 
1885  if (iomode == IO_MODE_FORTRAN) then
1886    close(ount)
1887  else
1888    NCF_CHECK(nf90_close(ount))
1889  end if
1890 
1891  write(msg,'(3a)')ch10,' ==== Files have been merged successfully === ',ch10
1892  call wrtout(std_out, msg)
1893 
1894 end subroutine ioscr_qmerge

m_io_screening/ioscr_qrecover [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 ioscr_qrecover

FUNCTION

  Recover q-points from a corrupted file produced e.g. from an interrupted run
  This routine should be called by a single MPI process.

INPUTS

  path=Corrupted file.
  nqrec=Number of q-points to recover.
  fname_out=Name of the file to be produced.

OUTPUT

  Output is written to file.

SOURCE

1917 subroutine ioscr_qrecover(ipath, nqrec, fname_out)
1918 
1919 !Arguments ------------------------------------
1920 !scalars
1921  integer,intent(in) :: nqrec
1922  character(len=*),intent(in) :: ipath,fname_out
1923 
1924 !Local variables-------------------------------
1925 !scalars
1926  integer,parameter :: rdwr2=2,master=0
1927  integer :: iqiA,nqibzA,nomega_asked,unt,npwe_asked,iomode,comm,fform1,ifform,ierr
1928  character(len=500) :: msg
1929  character(len=nctk_slen) :: varname
1930  type(hscr_t) :: hscr_recov,hscr
1931  type(abifile_t) :: abifile
1932 !arrays
1933  complex(gwpc),allocatable :: epsm1(:,:,:,:)
1934 
1935 ! *************************************************************************
1936 
1937  comm = xmpi_comm_self
1938 
1939  call wrtout(std_out, sjoin(". Recovering q-points in file:", ipath))
1940  call wrtout(std_out, sjoin(". Data written to file:", fname_out))
1941 
1942  if (file_exists(fname_out)) then
1943    ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
1944  end if
1945 
1946  ! Find iomode from file extension and open output file.
1947  if (endswith(fname_out, ".nc")) then
1948    iomode = IO_MODE_ETSF
1949    NCF_CHECK(nctk_open_create(unt, fname_out, comm))
1950  else
1951    iomode = IO_MODE_FORTRAN
1952    if (open_file(fname_out, msg, newunit=unt, status='new', form='unformatted') /= 0) then
1953      ABI_ERROR(msg)
1954    end if
1955  end if
1956 
1957  ! Read header.
1958  call hscr_from_file(hscr, ipath, ifform, comm)
1959  ABI_CHECK(ifform /= 0, sjoin("fform = 0 while reading:", ipath))
1960 
1961  if (nqrec < 1 .or. nqrec > hscr%nqibz) then
1962    ABI_ERROR(sjoin("Wrong input. nqibz on file:", itoa(hscr%nqibz)))
1963  end if
1964 
1965  ! Copy header
1966  call hscr_copy(hscr, hscr_recov)
1967 
1968  ! Change dimensions and arrays associated to nqibz.
1969  hscr_recov%nqibz = nqrec
1970  ABI_FREE(hscr_recov%qibz)
1971  ABI_MALLOC(hscr_recov%qibz, (3,nqrec))
1972  hscr_recov%qibz = hscr%qibz(:,1:nqrec)
1973 
1974  call hscr_print(hscr_recov,header="Header of the new SCR file",unit=std_out,prtvol=1)
1975 
1976  ! Write the header of the recovered file.
1977  fform1 = hscr%fform
1978 
1979  abifile = abifile_from_fform(fform1)
1980  if (abifile%fform == 0) then
1981     ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform1)))
1982  end if
1983  varname = abifile%varname
1984 
1985  call hscr_io(hscr_recov,fform1,rdwr2,unt,comm,master,iomode)
1986 
1987  nqibzA=1; nomega_asked=hscr%nomega; npwe_asked=hscr%npwe
1988 
1989  ABI_MALLOC_OR_DIE(epsm1,(npwe_asked,npwe_asked,nomega_asked,1), ierr)
1990 
1991  do iqiA=1,hscr_recov%nqibz
1992    call read_screening(varname,ipath,npwe_asked,nqibzA,nomega_asked,epsm1,iomode,comm,iqiA=iqiA)
1993    call write_screening(varname,unt,iomode,npwe_asked,nomega_asked,iqiA,epsm1)
1994  end do
1995 
1996  if (iomode == IO_MODE_FORTRAN) close(unt)
1997  if (iomode == IO_MODE_ETSF) then
1998    NCF_CHECK(nf90_close(unt))
1999  end if
2000 
2001  ABI_FREE(epsm1)
2002  call hscr%free()
2003  call hscr_recov%free()
2004 
2005  call wrtout(std_out, "Recovery completed")
2006 
2007 end subroutine ioscr_qrecover

m_io_screening/ioscr_wmerge [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 ioscr_wmerge

FUNCTION

  Produce new file by merging the frequencies stored in other files.
  This routine should be called by a single MPI process.

INPUTS

  nfiles=Number of files to be merged.
  filenames(nfiles)=Paths of files to be merged.
  hscr_files(nfiles)<hscr_t>=Heades of the files to be merged.
  fname_out=Name of the file to be produced.

OUTPUT

  ohscr<hscr_t>=The header of the output file.

SOURCE

2031 subroutine ioscr_wmerge(nfiles, filenames, hscr_file, freqremax, fname_out, ohscr)
2032 
2033 !Arguments ------------------------------------
2034 !scalars
2035  integer,intent(in) :: nfiles
2036  real(dp),intent(in) :: freqremax
2037  character(len=*),intent(in) :: fname_out
2038  type(hscr_t),intent(out) :: ohscr
2039 !arrays
2040  character(len=*),intent(in) :: filenames(nfiles)
2041  type(hscr_t),intent(in) :: hscr_file(nfiles)
2042 
2043 !Local variables-------------------------------
2044 !scalars
2045  integer,parameter :: rdwr2=2,master=0
2046  integer :: ii,iqibz,ifile,ount,iomode,fform_merge,comm,nomega4m
2047  integer :: nfreq_tot,nfreqre,nfreqim,ifrq,npwe4mI,npwe4mJ,ierr
2048  character(len=500) :: msg
2049  logical :: skip
2050  character(len=nctk_slen) :: varname
2051  type(abifile_t) :: abifile
2052 !arrays
2053  integer,allocatable :: freq_indx(:,:),ifile_indx(:),pos_indx(:),i_temp(:),i2_temp(:,:)
2054  real(dp),allocatable :: real_omega(:),imag_omega(:)
2055  complex(gwpc),allocatable :: epsm1(:,:,:,:),epsm1_temp(:,:,:,:)
2056  complex(dpc),allocatable :: omega_storage(:)
2057 
2058 ! *************************************************************************
2059 
2060  comm = xmpi_comm_self
2061 
2062  ! Check that q-point sets are the same
2063  do ifile=1,nfiles
2064    do ii=1,nfiles
2065      if (ii == ifile) CYCLE
2066      if (Hscr_file(ifile)%nqibz /= Hscr_file(ii)%nqibz) then
2067        ABI_ERROR(' One or more files do not have the same number of q-points!')
2068      end if
2069      do iqibz=1,Hscr_file(1)%nqibz
2070        if (ABS(SUM(Hscr_file(ifile)%qibz(:,iqibz)-Hscr_file(ii)%qibz(:,iqibz))) > tol6) then
2071          ABI_ERROR('Q-point set differs between one or more files!')
2072        end if
2073      end do
2074    end do
2075  end do
2076 
2077  ! nfreq_tot here is the total *possible* number of freq.
2078  nfreq_tot=0
2079  do ifile=1,nfiles
2080    nfreq_tot = nfreq_tot + Hscr_file(ifile)%nomega
2081  end do
2082 
2083  ABI_MALLOC(omega_storage, (nfreq_tot))
2084  ABI_MALLOC(freq_indx, (nfreq_tot,nfiles))
2085  ABI_MALLOC(ifile_indx, (nfreq_tot))
2086  omega_storage = CMPLX(-one,-one); freq_indx = 0; ifile_indx = 0
2087 
2088  ! Calculate the total number of real freq and store
2089  nfreqre = 0
2090  do ifile=1,nfiles
2091    do ifrq=1,Hscr_file(ifile)%nomega
2092      skip = .FALSE.
2093      ! Check whether to skip this point
2094      if (AIMAG(Hscr_file(ifile)%omega(ifrq)) > tol16) skip = .TRUE.
2095      if (REAL(Hscr_file(ifile)%omega(ifrq)) > freqremax) skip = .TRUE.
2096      ! Check for repetition or non-monotonic points
2097      if (nfreqre>1) then
2098        do ii=1,nfreqre
2099          if (ABS(REAL(Hscr_file(ifile)%omega(ifrq)) - REAL(omega_storage(ii))) < tol6) skip = .TRUE.
2100        end do
2101      end if
2102      if (skip) CYCLE
2103      nfreqre = nfreqre + 1
2104      ! Store (complex) frequency and index
2105      omega_storage(nfreqre) = Hscr_file(ifile)%omega(ifrq)
2106      ifile_indx(nfreqre) = ifile
2107      freq_indx(nfreqre,ifile) = ifrq
2108      write(std_out,'(a,es16.6,a,i0,2(a,i0))')&
2109        ' Found real frequency: ',REAL(omega_storage(nfreqre))*Ha_eV,' [eV], number: ',nfreqre,&
2110        ', in file: ',ifile,' local index: ',ifrq
2111    end do
2112  end do
2113 
2114  if (nfreqre > 0) then
2115    ! Sort real frequencies
2116    ABI_MALLOC(real_omega,(nfreqre))
2117    ABI_MALLOC(pos_indx,(nfreqre))
2118    ABI_MALLOC(i_temp,(nfreqre))
2119    ABI_MALLOC(i2_temp,(nfreqre,nfiles))
2120    real_omega(1:nfreqre) = REAL(omega_storage(1:nfreqre)) ! Copy real frequencies to temp. sorting array
2121    do ii=1,nfreqre ! Set up indexing array
2122      pos_indx(ii) = ii
2123    end do
2124 
2125    ! Sort frequencies while keeping track of index
2126    call sort_dp(nfreqre,real_omega,pos_indx,tol16)
2127    i_temp(1:nfreqre) = ifile_indx(1:nfreqre)
2128    i2_temp(1:nfreqre,1:nfiles) = freq_indx(1:nfreqre,1:nfiles)
2129 
2130    ! Copy sorted frequencies plus file and frequency index
2131    do ii=1,nfreqre
2132      omega_storage(ii) = CMPLX(real_omega(ii),zero)
2133      ifile_indx(ii) = i_temp(pos_indx(ii))
2134      freq_indx(ii,1:nfiles) = i2_temp(pos_indx(ii),1:nfiles)
2135    end do
2136 
2137    ABI_FREE(real_omega)
2138    ABI_FREE(pos_indx)
2139    ABI_FREE(i_temp)
2140    ABI_FREE(i2_temp)
2141  end if
2142 
2143  ! Check imaginary frequencies and store them
2144  nfreqim = 0
2145  do ifile=1,nfiles
2146    do ifrq=1,Hscr_file(ifile)%nomega
2147      if (REAL(Hscr_file(ifile)%omega(ifrq)) > tol8) CYCLE
2148      if (AIMAG(Hscr_file(ifile)%omega(ifrq)) < tol8) CYCLE
2149      nfreqim = nfreqim + 1
2150      omega_storage(nfreqre+nfreqim) = Hscr_file(ifile)%omega(ifrq)
2151      ifile_indx(nfreqre+nfreqim) = ifile
2152      freq_indx(nfreqre+nfreqim,ifile) = ifrq
2153      write(std_out,'(a,es16.6,a,i0,2(a,i0))')&
2154       ' Found imaginary frequency: ',AIMAG(omega_storage(nfreqre+nfreqim))*Ha_eV,' [eV], number: ',nfreqim,&
2155       ', in file: ',ifile,' local index: ',ifrq
2156    end do
2157  end do
2158 
2159  ! Sort imaginary frequencies
2160  ABI_MALLOC(imag_omega,(nfreqim))
2161  ABI_MALLOC(pos_indx,(nfreqim))
2162  ABI_MALLOC(i_temp,(nfreqim))
2163  ABI_MALLOC(i2_temp,(nfreqim,nfiles))
2164 
2165  ! Copy imaginary frequencies to temp. sorting array
2166  imag_omega(1:nfreqim) = AIMAG(omega_storage(nfreqre+1:nfreqre+nfreqim))
2167  do ii=1,nfreqim ! Set up indexing array
2168    pos_indx(ii) = ii
2169  end do
2170 
2171  ! Sort frequencies while keeping track of index
2172  call sort_dp(nfreqim,imag_omega,pos_indx,tol16)
2173  i_temp(1:nfreqim) = ifile_indx(nfreqre+1:nfreqre+nfreqim)
2174  i2_temp(1:nfreqim,1:nfiles) = freq_indx(nfreqre+1:nfreqre+nfreqim,1:nfiles)
2175 
2176  ! Copy sorted frequencies plus file and frequency index
2177  do ii=1,nfreqim
2178    omega_storage(nfreqre+ii) = CMPLX(zero,imag_omega(ii))
2179    ifile_indx(nfreqre+ii) = i_temp(pos_indx(ii))
2180    freq_indx(nfreqre+ii,1:nfiles) = i2_temp(pos_indx(ii),1:nfiles)
2181  end do
2182 
2183  ABI_FREE(imag_omega)
2184  ABI_FREE(pos_indx)
2185  ABI_FREE(i_temp)
2186  ABI_FREE(i2_temp)
2187 
2188  nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq
2189  write(std_out,'(2a,i0,a)') ch10,' Merging ',nfreq_tot,' frequencies.'
2190  write(std_out,'(2(a,i0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10
2191 
2192  ! Copy old header
2193  call hscr_copy(Hscr_file(1),ohscr)
2194 
2195  ! TODO: hscr_wmerge
2196  ! Then modify entries for new frequency grid
2197  ohscr%nomega = nfreq_tot
2198  ABI_FREE(ohscr%omega)
2199  ABI_MALLOC(ohscr%omega,(nfreq_tot))
2200  ohscr%omega(1:nfreq_tot) = omega_storage(1:nfreq_tot)
2201 
2202  npwe4mI = ohscr%npwe*ohscr%nI
2203  npwe4mJ = ohscr%npwe*ohscr%nJ
2204 
2205  ! Print new header for info
2206  call hscr_print(ohscr,header='Header of the final file',unit=std_out,prtvol=1)
2207 
2208  if (file_exists(fname_out)) then
2209    ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
2210  end if
2211 
2212  if (endswith(fname_out, ".nc")) then
2213    iomode = IO_MODE_ETSF
2214    NCF_CHECK(nctk_open_create(ount, fname_out, comm))
2215  else
2216    iomode = IO_MODE_FORTRAN
2217    if (open_file(fname_out, msg, newunit=ount, status='new',form='unformatted') /= 0) then
2218      ABI_ERROR(msg)
2219    end if
2220  end if
2221 
2222  ! * Write the header.
2223  fform_merge = ohscr%fform
2224 
2225  abifile = abifile_from_fform(fform_merge)
2226  if (abifile%fform == 0) then
2227     ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
2228  end if
2229  varname = abifile%varname
2230 
2231  call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode)
2232 
2233  npwe4mI = ohscr%npwe*ohscr%nI
2234  npwe4mJ = ohscr%npwe*ohscr%nJ
2235  nomega4m = ohscr%nomega
2236  ABI_MALLOC_OR_DIE(epsm1,(npwe4mI,npwe4mJ,nomega4m,1), ierr)
2237 
2238  do iqibz=1,ohscr%nqibz
2239 
2240    do ifile=1,nfiles
2241      ! allocate temporary array
2242      npwe4mI = Hscr_file(ifile)%npwe*Hscr_file(ifile)%nI
2243      npwe4mJ = Hscr_file(ifile)%npwe*Hscr_file(ifile)%nJ
2244      nomega4m = Hscr_file(ifile)%nomega
2245      ABI_MALLOC_OR_DIE(epsm1_temp,(npwe4mI,npwe4mJ,nomega4m,1), ierr)
2246 
2247      ! read screening
2248      call read_screening(varname,filenames(ifile),npwe4mI,1,nomega4m,epsm1_temp,iomode,comm,iqiA=iqibz)
2249 
2250      ! Copy matrices for relevant frequencies
2251      do ifrq=1,nfreq_tot
2252        if (ifile_indx(ifrq)==ifile) then
2253          epsm1(:,:,ifrq,1)=epsm1_temp(:,:,freq_indx(ifrq,ifile),1)
2254        end if
2255      end do
2256 
2257      !Add imaginary part if needed
2258      !if (indx_imfreq_file==ifile) then
2259      !do ifrq=nfreqre+1,ohscr%nomega
2260      !epsm1(:,:,ifrq,1)=epsm1_temp(:,:,freq_indx(ifrq,ifile),1)
2261      !end do
2262      !end if
2263 
2264      ABI_FREE(epsm1_temp)
2265    end do !ifile
2266 
2267    ! Write data.
2268    npwe4mI = ohscr%npwe*ohscr%nI
2269    nomega4m = ohscr%nomega
2270    call write_screening(varname,ount,iomode,npwe4mI,nomega4m,iqibz,epsm1)
2271  end do ! iqibz
2272 
2273  ABI_FREE(epsm1)
2274  ABI_FREE(omega_storage)
2275  ABI_FREE(freq_indx)
2276  ABI_FREE(ifile_indx)
2277 
2278  if (iomode == IO_MODE_FORTRAN) then
2279    close(ount)
2280  else
2281    NCF_CHECK(nf90_close(ount))
2282  end if
2283 
2284  write(msg,'(3a)')ch10,' ==== Files have been merged successfully === ',ch10
2285  call wrtout(std_out, msg)
2286 
2287 end subroutine ioscr_wmerge

m_io_screening/ioscr_wremove [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 ioscr_wremove

FUNCTION

  Produce new file by removing selected frequencies in the initial file `inpath`.
  This routine should be called by a single MPI process.

INPUTS

  inpath=Input file
  ihscr<hscr_t>=Headerf of the input file.
  fname_out=Output file.
  nfreq_tot=Number of frequencies in new file.
  freq_indx(nfreq_tot)=Index of frequency to be kept in input file.

OUTPUT

  ohscr<hscr_t>=The header of the output file.

SOURCE

2312 subroutine ioscr_wremove(inpath, ihscr, fname_out, nfreq_tot, freq_indx, ohscr)
2313 
2314 !Arguments ------------------------------------
2315 !scalars
2316  integer,intent(in) :: nfreq_tot
2317  character(len=*),intent(in) :: inpath,fname_out
2318  type(hscr_t),intent(in) :: ihscr
2319  type(hscr_t),intent(out) :: ohscr
2320 !arrays
2321  integer,intent(in) :: freq_indx(nfreq_tot)
2322 
2323 !Local variables-------------------------------
2324 !scalars
2325  integer,parameter :: rdwr2=2,master=0
2326  integer :: iqibz,fform_merge,comm,nomega4m,ierr
2327  integer :: ifrq,npwe4mI,npwe4mJ,iomode,ount
2328  character(len=500) :: msg
2329  character(len=nctk_slen) :: varname
2330  type(abifile_t) :: abifile
2331 !arrays
2332  complex(gwpc),allocatable :: epsm1(:,:,:),epsm1_temp(:,:,:)
2333 
2334 ! *************************************************************************
2335 
2336  comm = xmpi_comm_self
2337 
2338  ! check ifreq_idx
2339  ABI_CHECK(nfreq_tot > 0, "nfreq_tot <= 0!")
2340  if (all(freq_indx == 0)) ABI_ERROR("all(freq_indx == 0)")
2341 
2342  ! Copy the old header
2343  call hscr_copy(ihscr, ohscr)
2344 
2345  ! Then modify entries for new frequency grid.
2346  ohscr%nomega = nfreq_tot
2347  ABI_FREE(ohscr%omega)
2348  ABI_MALLOC(ohscr%omega,(nfreq_tot))
2349  do ifrq=1,nfreq_tot
2350    ohscr%omega(ifrq) = ihscr%omega(freq_indx(ifrq))
2351  end do
2352 
2353  npwe4mI = ohscr%npwe*ohscr%nI
2354  npwe4mJ = ohscr%npwe*ohscr%nJ
2355 
2356  ! Print new header for info
2357  call hscr_print(ohscr,header='Header of the final file',unit=std_out,prtvol=1)
2358 
2359  ! Open output file.
2360  if (endswith(fname_out, ".nc")) then
2361    iomode = IO_MODE_ETSF
2362    NCF_CHECK(nctk_open_create(ount, fname_out, comm))
2363  else
2364    iomode = IO_MODE_FORTRAN
2365    if (open_file(fname_out, msg, newunit=ount, status='new', form='unformatted') /= 0) then
2366      ABI_ERROR(msg)
2367    end if
2368  end if
2369 
2370  ! Write the header.
2371  fform_merge = ohscr%fform
2372  abifile = abifile_from_fform(fform_merge)
2373  if (abifile%fform == 0) then
2374     ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
2375  end if
2376  varname = abifile%varname
2377 
2378  call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode)
2379 
2380  npwe4mI = ohscr%npwe*ohscr%nI; npwe4mJ = ohscr%npwe*ohscr%nJ
2381  nomega4m = ohscr%nomega
2382 
2383  ABI_MALLOC_OR_DIE(epsm1, (npwe4mI,npwe4mJ,nomega4m), ierr)
2384 
2385  do iqibz=1,ohscr%nqibz
2386    ! allocate temporary array
2387    npwe4mI = ihscr%npwe * ihscr%nI
2388    npwe4mJ = ihscr%npwe * ihscr%nJ
2389    nomega4m = ihscr%nomega
2390    ABI_MALLOC_OR_DIE(epsm1_temp,(npwe4mI,npwe4mJ,nomega4m), ierr)
2391 
2392    ! read full screening matrix for this q-point
2393    call read_screening(varname,inpath,npwe4mI,1,nomega4m,epsm1_temp,iomode,comm,iqiA=iqibz)
2394 
2395    ! Copy relevant frequencies
2396    do ifrq=1,nfreq_tot
2397      epsm1(:,:,ifrq) = epsm1_temp(:,:,freq_indx(ifrq))
2398    end do
2399 
2400    ABI_FREE(epsm1_temp)
2401 
2402    npwe4mI = ohscr%npwe*ohscr%nI; nomega4m = ohscr%nomega
2403    call write_screening(varname,ount,iomode,npwe4mI,nomega4m,iqibz,epsm1)
2404  end do ! iqibz
2405 
2406  ABI_FREE(epsm1)
2407 
2408  if (iomode == IO_MODE_FORTRAN) then
2409    close(ount)
2410  else
2411    NCF_CHECK(nf90_close(ount))
2412  end if
2413 
2414  write(msg,'(3a)')ch10,' ==== Frequencies have been removed successfully === ',ch10
2415  call wrtout(std_out, msg)
2416 
2417 end subroutine ioscr_wremove

m_io_screening/ncname_from_id [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

  ncname_from_id

FUNCTION

  Return the name of the netcdf variable (chi0, espilon, em1...) from the id.

SOURCE

279 character(len=nctk_slen) function ncname_from_id(id) result(varname)
280 
281   integer,intent(in) :: id
282 
283   varname = "None"
284   if (id == 1) varname = chi0_ncname
285   if (id == 3) varname = e_ncname
286   if (id == 4) varname = em1_ncname
287   ABI_CHECK(varname /= "None", "Wrong id")
288 
289 end function ncname_from_id

m_io_screening/read_screening [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 read_screening

FUNCTION

 Read either a screening (\tilde epsilon^{-1}) file in the SCR format or
 the irreducible polarizability (chi0) in the SUSC format.

INPUTS

  varname=Name of the array to read. Used for ETSF-IO files.
  iomode=Integer flag defining the format of the output file. Available options:
    IO_MODE_FORTRAN--> Plain Fortran file
    IO_MODE_ETSF--> ETSF format
  iqiA[optional]=Used if only a particular q-point is required. In this case iqiA define the index
   of the required q-point in the array qibz(3,Hscr%nqibz)
  nqibzA=number of asked q-points (used to dimension the output arrays).
   Equal to Hscr%nqibz if the full matrix is required
  comm=MPI communicator.
  npweA=number of asked planewaves
  nomegaA=number of asked frequencies

OUTPUT

  epsm1(npweA,npweA,nomegaA,nqibzA) = \tilde\epsilon^{-1}(Ng,Ng,Nw,Nq)

NOTES

  * If the epsilon matrix read is bigger than npweA x npweA, it will be truncated;
    if it is smaller, an error will occur
  * If the number of frequencies asked for is smaller than that reported in the file, the matrix
    will be truncated. If nomegaA > Hscr%nomega an error will occur

SOURCE

1400 subroutine read_screening(varname,fname,npweA,nqibzA,nomegaA,epsm1,iomode,comm, &
1401 & iqiA) ! Optional
1402 
1403 !Arguments ------------------------------------
1404 !scalars
1405  integer,intent(in) :: iomode,nomegaA,npweA,nqibzA,comm
1406  integer,optional,intent(in) :: iqiA
1407  character(len=*),intent(in) :: varname,fname
1408 !arrays
1409  complex(gwpc),target,intent(inout) :: epsm1(npweA,npweA,nomegaA,nqibzA)
1410 
1411 !Local variables-------------------------------
1412 !scalars
1413  integer,parameter :: master = 0
1414  integer :: ipwe,fform,iomega,iqibz,unt,rdwr,my_rank,nprocs,my_iomode
1415  integer :: varid,ncerr
1416 #ifdef HAVE_MPI_IO
1417  integer :: test_fform,mpi_err,ierr,sc_mode
1418  integer :: bsize_frm,mpi_type_frm
1419  integer :: mpi_fh,buf_dim !,mat_ggw,mat_ggwq
1420  integer(XMPI_OFFSET_KIND) :: offset,displ_wq !,my_offpad
1421  !complex(dpc) :: ctmp
1422 #endif
1423  character(len=500) :: msg,errmsg
1424  logical :: read_qslice
1425  type(hscr_t) :: Hscr
1426 !arrays
1427 #ifdef HAVE_MPI_IO
1428  integer(MPI_OFFSET_KIND),allocatable :: offset_wq(:,:)
1429 #endif
1430  complex(dpc),allocatable :: bufdc2d(:,:),bufdc3d(:,:,:)
1431  ! pointers passed to netcdf4 routines (complex datatypes are not supported).
1432 #ifdef HAVE_GW_DPC
1433  real(dp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:)
1434 #else
1435  real(sp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:)
1436 #endif
1437  integer :: spins(2),s1,s2
1438 
1439 ! *************************************************************************
1440 
1441  DBG_ENTER("COLL")
1442 
1443  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1444  my_iomode = iomode
1445  if (endswith(fname, ".nc")) my_iomode = IO_MODE_ETSF
1446  !my_iomode = IO_MODE_MPI
1447  !if (my_iomode  == IO_MODE_MPI) my_iomode = IO_MODE_FORTRAN
1448 
1449  rdwr=1
1450  select case (my_iomode)
1451  case (IO_MODE_MPI)
1452 #ifdef HAVE_MPI_IO
1453    bsize_frm    = xmpio_bsize_frm    ! bsize_frm= Byte length of the Fortran record marker.
1454    mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker.
1455    sc_mode = xmpio_collective
1456 
1457    ! Master reads the header via Fortran IO then bcast the data.
1458    call hscr_from_file(hscr, fname, fform, comm)
1459 
1460    ! Open the file with MPI-IO
1461    call MPI_FILE_OPEN(comm, fname, MPI_MODE_RDONLY, xmpio_info ,mpi_fh, mpi_err)
1462    ABI_CHECK_MPI(mpi_err, sjoin("MPI_FILE_OPEN:", fname))
1463 
1464    ! Retrieve the offset of the section immediately below the header.
1465    call hscr_mpio_skip(mpi_fh,test_fform,offset)
1466    ABI_CHECK(test_fform == fform, "mismatch in fform!")
1467 
1468    ! Offsets of the Fortran markers corresponding to the (w,q) slices.
1469    ABI_MALLOC(offset_wq,(HScr%nomega,HScr%nqibz))
1470    displ_wq = offset
1471    do iqibz=1,Hscr%nqibz
1472      do iomega=1,Hscr%nomega
1473        ABI_CHECK(displ_wq > 0, "displ_wq < 0, your SCR|SUSC file is too big for MPI-IO!")
1474        offset_wq(iomega,iqibz) = displ_wq
1475        displ_wq = displ_wq + Hscr%npwe**2 * xmpi_bsize_dpc + Hscr%npwe * 2 * bsize_frm
1476      end do
1477    end do
1478 #else
1479    ABI_ERROR("MPI-IO support not enabled at configure-time")
1480 #endif
1481 
1482  case (IO_MODE_FORTRAN)
1483    ! Plain Fortran IO, all nodes read.
1484    if (open_file(fname,msg,newunit=unt,form="unformatted",status="old",action="read") /= 0) then
1485      ABI_ERROR(msg)
1486    end if
1487    call hscr_io(hscr,fform,rdwr,unt,comm,master,my_iomode)
1488 
1489  case (IO_MODE_ETSF)
1490    NCF_CHECK(nctk_open_read(unt, fname, xmpi_comm_self))
1491    call hscr_io(hscr,fform,rdwr,unt,comm,master,my_iomode)
1492 
1493  case default
1494    ABI_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode)))
1495  end select
1496 
1497  ! Slice or full array?
1498  read_qslice = .False.
1499  if (PRESENT(iqiA)) then
1500    read_qslice = .True.
1501    !call wrtout(std_out, sjoin('. Reading q-slice for iq = ',itoa(iqiA),' from: ', fname))
1502    if (iqiA <= 0 .or. iqiA > Hscr%nqibz) then
1503      ABI_BUG('iqiA out of range')
1504    end if
1505  end if
1506 
1507  ! Do some check
1508  if (Hscr%npwe>npweA) then
1509    write(msg,'(a,i0,2a,i0)')&
1510     'Total number of G-vectors reported on file = ',Hscr%npwe,ch10,&
1511     'Reading a smaller matrix of dimension      = ',npweA
1512    ABI_COMMENT(msg)
1513  end if
1514 
1515  if (npweA>Hscr%npwe) then
1516    write(msg,'(2(a,i0))')' Dimension of matrix = ',Hscr%npwe," requiring a too big matrix = ",npweA
1517    ABI_ERROR(msg)
1518  end if
1519 
1520  ABI_CHECK(nqibzA  <= Hscr%nqibz, 'Requiring too many q-points')
1521  ABI_CHECK(nomegaA <= Hscr%nomega,'Requiring too many frequencies')
1522 
1523  select case (my_iomode)
1524  case (IO_MODE_MPI)
1525 #ifdef HAVE_MPI_IO
1526    if (read_qslice) then
1527       call wrtout(std_out, "calling mpiotk to read_qslice")
1528       buf_dim = (npweA)**2 * nomegaA
1529       offset = offset_wq(1,iqiA)
1530       sc_mode = xmpio_collective
1531 
1532 #ifdef HAVE_GW_DPC
1533      ! Read in-place.
1534      call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1535         buf_dim,epsm1,xmpio_chunk_bsize,sc_mode,comm,ierr)
1536      ABI_CHECK(ierr==0,"Fortran matrix too big")
1537 #else
1538      ! Have to allocate workspace for dpc data.
1539      ! FIXME: Change the file format of the SCR and SUC file so that
1540      ! they are written in single precision if not HAVE_GW_DPC
1541      ABI_MALLOC_OR_DIE(bufdc3d,(npweA,npweA,nomegaA), ierr)
1542 
1543      call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1544         buf_dim,bufdc3d,xmpio_chunk_bsize,sc_mode,comm,ierr)
1545      ABI_CHECK(ierr==0,"Fortran matrix too big")
1546 
1547      epsm1(:,:,:,1) = bufdc3d
1548      ABI_FREE(bufdc3d)
1549 #endif
1550 
1551    else
1552      ! Full matrix (G,G',w,q) is needed.
1553      call wrtout(std_out, "calling mpiotk: Full matrix (G,G',w,q) is needed.")
1554 
1555 #ifdef HAVE_GW_DPC
1556      ! Can read all data at once.
1557      buf_dim = (npweA)**2 * nomegaA * HScr%nqibz
1558      offset = offset_wq(1,1)
1559      sc_mode = xmpio_collective
1560 
1561      call mpiotk_read_fsuba_dpc4D(mpi_fh,offset,&
1562        [HScr%npwe,HScr%npwe,HScr%nomega,HScr%nqibz], [npweA,npweA,nomegaA,HScr%nqibz], [1,1,1,1],&
1563        buf_dim,epsm1,xmpio_chunk_bsize,sc_mode,comm,ierr)
1564      ABI_CHECK(ierr==0,"Fortran record too big")
1565 #else
1566      ! Have to allocate workspace for dpc data.
1567      ABI_MALLOC_OR_DIE(bufdc3d,(npweA,npweA,nomegaA), ierr)
1568      sc_mode = xmpio_collective
1569 
1570      do iqibz=1,Hscr%nqibz
1571        offset = offset_wq(1,iqibz)
1572        buf_dim = (2*npweA)**2 * nomegaA
1573 
1574        call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, &
1575         [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1576          buf_dim,bufdc3d,xmpio_chunk_bsize,sc_mode,comm,ierr)
1577        ABI_CHECK(ierr==0,"Fortran matrix too big")
1578 
1579        epsm1(:,:,:,iqibz) = bufdc3d
1580      end do
1581 
1582      ABI_FREE(bufdc3d)
1583 #endif
1584    end if
1585 
1586    call MPI_FILE_CLOSE(mpi_fh,mpi_err)
1587    ABI_FREE(offset_wq)
1588 #endif
1589 
1590  case (IO_MODE_FORTRAN)
1591    ! Read epsilon^-1 with Fortran IO
1592    ! Allocate a single column to save memory.
1593    ! TODO re-merge the two cases.
1594    ABI_MALLOC(bufdc2d,(Hscr%npwe,1))
1595 
1596    ! Two coding for different case just to keep it readable.
1597    select case (read_qslice)
1598    case (.True.)
1599      ! Read only a slice of the full array (useful if the entire array is huge).
1600      !if (dim_wings==1) STOP 'not implemented'
1601      !TODO this has to be done in a cleaner way.
1602      qread_loop: &
1603 &    do iqibz=1,Hscr%nqibz
1604        if (iqibz==iqiA) then
1605          do iomega=1,nomegaA
1606            do ipwe=1,Hscr%npwe
1607              read(unt, err=10, iomsg=errmsg) bufdc2d(1:Hscr%npwe,1)
1608              if (ipwe<=npweA) epsm1(1:npweA,ipwe,iomega,1)=bufdc2d(1:npweA,1)
1609            end do
1610          end do
1611          EXIT qread_loop ! Got data. Do not need to read file till the end.
1612        else
1613          ! Skip other q-points i.e bufdc2d(1:Hscr%npwe,1:Hscr%npwe)
1614          do iomega=1,Hscr%nomega
1615            do ipwe=1,Hscr%npwe
1616             read(unt, err=10, iomsg=errmsg)
1617            end do
1618          end do
1619        end if ! iqibz==iqiA
1620      end do qread_loop ! iqibz
1621 
1622    case (.False.)
1623      ! Read the entire array.
1624      do iqibz=1,Hscr%nqibz
1625        do iomega=1,nomegaA
1626          do ipwe=1,Hscr%npwe
1627           read(unt, err=10, iomsg=errmsg) bufdc2d(1:Hscr%npwe,1)
1628           if (ipwe<=npweA) epsm1(1:npweA,ipwe,iomega,iqibz)=bufdc2d(1:npweA,1)
1629          end do
1630        end do
1631        ! Skip other frequencies
1632        do iomega=nomegaA+1,Hscr%nomega
1633          do ipwe=1,Hscr%npwe
1634            read(unt, err=10, iomsg=errmsg)
1635          end do
1636        end do
1637      end do !iqibz
1638    end select
1639 
1640    close(unt)
1641 
1642  case (IO_MODE_ETSF)
1643    ! netcdf does not support complex datatypes. Here I use some C-magic to  associate the memory
1644    ! to a Fortran real pointer with the correct type and shape. Note that the data on file is always in double precision.
1645    ! nf90_get_var will automatically convert from double to single if the GW code is in single precision mode.
1646    ! This is the reason why I'm using CPP option in the declaration of real_epsm1.
1647 
1648    ! FIXME: Need to know the type to read
1649    !write(std_out,*)"in read_screening"
1650    varid = nctk_idname(unt, varname)
1651 
1652    ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt]
1653    call c_f_pointer(c_loc(epsm1(1,1,1,1)), real_epsm1, [2,npweA,npweA,1,1,nomegaA,nqibzA])
1654    spins = 1; s1 = spins(1); s2 = spins(2)
1655    if (read_qslice) then
1656      ncerr = nf90_get_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,iqia], count=[2,npweA,npweA,1,1,nomegaA,1])
1657    else
1658      ncerr = nf90_get_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,1], count=[2,npweA,npweA,1,1,nomegaA,nqibzA])
1659      !do iqibz=1,nqibzA
1660      !  !write(*,*)"epsm1: in read ",iqibz,epsm1(1:3,1,1,iqibz)
1661      !end do
1662    end if
1663    NCF_CHECK_MSG(ncerr, sjoin("getting var:", varname))
1664    NCF_CHECK(nf90_close(unt))
1665    !write(std_out,*)"read_screening done"
1666 
1667  case default
1668    ABI_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode)))
1669  end select
1670 
1671  ! Free memory
1672  ABI_SFREE(bufdc2d)
1673  ABI_SFREE(bufdc3d)
1674 
1675  call Hscr%free()
1676 
1677  DBG_EXIT("COLL")
1678 
1679  return
1680 
1681  ! Handle Fortran IO error.
1682 10 continue
1683  ABI_ERROR(errmsg)
1684 
1685 end subroutine read_screening

m_io_screening/write_screening [ Functions ]

[ Top ] [ m_io_screening ] [ Functions ]

NAME

 write_screening

FUNCTION

 For a single q-point, write either \tilde epsilon^{-1} on the _SCR file
 or chi0 on the _SUSC file. The file is supposed to have been open in the calling routine.

INPUTS

  varname=The name of the array to write (used if etsf-io format).
  unt=The unit number of the file to be written (supposed to be already open)
  iomode=Integer flag defining the format of the output file. Available options:
    IO_MODE_FORTRAN--> Plain Fortran file
    IO_MODE_ETSF--> ETSF format
  npwe=Number of plane waves in epsm1.
  nomega=Number of frequencies
  iqibz=Index of the q-points in the IBZ.
  epsm1(npwe,npwe,nomega)=The matrix to be written, for different frequencies, and a single q-point.

NOTES

  On some architecture, the code crashes when trying to write or read a record containing the
  entire (G1,G2) matrix thus we use smaller records containing the columns of the two-point function.

OUTPUT

  (only writing on file)

SOURCE

1300 subroutine write_screening(varname, unt, iomode, npwe, nomega, iqibz, epsm1)
1301 
1302 !Arguments ------------------------------------
1303 !scalars
1304  character(len=*),intent(in) :: varname
1305  integer,intent(in) :: nomega,npwe,iqibz,unt,iomode
1306 !arrays
1307  complex(gwpc),target,intent(in) :: epsm1(npwe,npwe,nomega)
1308 
1309 !Local variables-------------------------------
1310 !scalars
1311  integer :: ipwe,iomega,spins(2),s1,s2
1312  character(len=500) :: errmsg
1313 !arrays
1314  complex(dpc),allocatable :: epsm1d(:,:)
1315  integer :: varid,ncerr
1316 #ifdef HAVE_GW_DPC
1317  real(dp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:)
1318 #else
1319  real(sp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:)
1320 #endif
1321 
1322 ! *************************************************************************
1323 
1324  DBG_ENTER("COLL")
1325 
1326  select case (iomode)
1327  case (IO_MODE_FORTRAN, IO_MODE_MPI)
1328    ! Write a record for each omega, Always use double precision.
1329    ABI_MALLOC(epsm1d,(npwe,1))
1330 
1331    do iomega=1,nomega
1332      do ipwe=1,npwe
1333        epsm1d(:,1) = epsm1(:,ipwe,iomega) !spc ==> dpc
1334        write(unt, err=10, iomsg=errmsg)epsm1d(1:npwe,1)
1335      end do
1336    end do
1337    ABI_FREE(epsm1d)
1338 
1339  case (IO_MODE_ETSF)
1340    ! netcdf does not support complex datatypes. Here I use some C-magic to  associate the memory
1341    ! to a Fortran real pointer with the correct type and shape. Note that the data on file is always in double precision.
1342    ! but this is ok since: if the type of data differs from the netCDF variable type, type conversion will occur
1343    ! inside nf90_put_var
1344    varid = nctk_idname(unt, varname)
1345    call c_f_pointer(c_loc(epsm1(1,1,1)), real_epsm1, [2, npwe, npwe, 1, 1, nomega, 1])
1346    ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt]
1347    spins = 1; s1 = spins(1); s2 = spins(2)
1348    ncerr = nf90_put_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,iqibz], count=[2,npwe,npwe,1,1,nomega,1])
1349    NCF_CHECK_MSG(ncerr, sjoin("putting var:", varname))
1350 
1351  case default
1352    ABI_ERROR(sjoin("Wrong iomode:", iomode2str(iomode)))
1353  end select
1354 
1355  DBG_EXIT("COLL")
1356 
1357  return
1358 
1359  ! Handle IO error
1360  10 continue
1361  ABI_ERROR(errmsg)
1362 
1363 end subroutine write_screening