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

PARENTS

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 MODULE m_io_screening
27 
28  use defs_basis
29  use defs_abitypes
30  use m_abicore
31 #if defined HAVE_MPI2
32  use mpi
33 #endif
34  use m_xmpi
35  use m_mpiotk
36  use m_nctk
37  use m_errors
38  use iso_c_binding
39 #ifdef HAVE_NETCDF
40  use netcdf
41 #endif
42  use m_hdr
43  use m_sort
44 
45  use m_gwdefs,          only : em1params_t, GW_TOLQ
46  use m_fstrings,        only : sjoin, itoa, endswith, replace_ch0
47  use m_copy,            only : alloc_copy
48  use m_io_tools,        only : open_file, file_exists, iomode2str
49  use m_numeric_tools,   only : print_arr, remove_copies, imax_loc
50  use m_bz_mesh,         only : isequalk
51 
52  implicit none
53 
54  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.

PARENTS

      m_io_screening,setup_bse

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

1016 subroutine hscr_bcast(hscr,master,my_rank,comm)
1017 
1018 
1019 !This section has been created automatically by the script Abilint (TD).
1020 !Do not modify the following lines by hand.
1021 #undef ABI_FUNC
1022 #define ABI_FUNC 'hscr_bcast'
1023 !End of the abilint section
1024 
1025  implicit none
1026 
1027 !Arguments ------------------------------------
1028  integer, intent(in) :: master,my_rank,comm
1029  type(hscr_t),intent(inout) :: hscr
1030 
1031 !Local variables-------------------------------
1032  integer :: ierr
1033 
1034 ! *************************************************************************
1035 
1036  DBG_ENTER("COLL")
1037  if (xmpi_comm_size(comm) == 1) return ! Nothing to do
1038 
1039  !@hscr_t
1040 ! integer
1041  call xmpi_bcast(hscr%id,         master,comm,ierr)
1042  call xmpi_bcast(hscr%ikxc,       master,comm,ierr)
1043  call xmpi_bcast(hscr%inclvkb,    master,comm,ierr)
1044  call xmpi_bcast(hscr%headform,   master,comm,ierr)
1045  call xmpi_bcast(hscr%fform,      master,comm,ierr)
1046  call xmpi_bcast(hscr%gwcalctyp,  master,comm,ierr)
1047  call xmpi_bcast(hscr%nI,         master,comm,ierr)
1048  call xmpi_bcast(hscr%nJ,         master,comm,ierr)
1049  call xmpi_bcast(hscr%nqibz,      master,comm,ierr)
1050  call xmpi_bcast(hscr%nqlwl,      master,comm,ierr)
1051  call xmpi_bcast(hscr%nomega,     master,comm,ierr)
1052  call xmpi_bcast(hscr%nbnds_used, master,comm,ierr)
1053  call xmpi_bcast(hscr%npwe,       master,comm,ierr)
1054  call xmpi_bcast(hscr%npwwfn_used,master,comm,ierr)
1055  call xmpi_bcast(hscr%spmeth,     master,comm,ierr)
1056  call xmpi_bcast(hscr%test_type,  master,comm,ierr)
1057  call xmpi_bcast(hscr%tordering,  master,comm,ierr)
1058 
1059  ! Real
1060  call xmpi_bcast(hscr%mbpt_sciss, master,comm,ierr)
1061  call xmpi_bcast(hscr%spsmear,    master,comm,ierr)
1062  call xmpi_bcast(hscr%zcut,       master,comm,ierr)
1063 
1064  ! arrays
1065  call xmpi_bcast(hscr%titles, master,comm,ierr)
1066 
1067  if (my_rank /= master) then
1068    call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
1069  end if
1070 
1071  call xmpi_bcast(hscr%gvec, master,comm,ierr)
1072  call xmpi_bcast(hscr%qibz, master,comm,ierr)
1073  call xmpi_bcast(hscr%qlwl, master,comm,ierr)
1074  call xmpi_bcast(hscr%omega,master,comm,ierr)
1075 
1076  ! Communicate the Abinit header.
1077  call hdr_bcast(hscr%Hdr,master,my_rank,comm)
1078 
1079 ! HSCR_NEW
1080  call xmpi_bcast(hscr%awtr, master, comm, ierr)
1081  call xmpi_bcast(hscr%icutcoul, master, comm, ierr)
1082  call xmpi_bcast(hscr%vcutgeo, master, comm, ierr)
1083  call xmpi_bcast(hscr%gwcomp, master, comm, ierr)
1084  call xmpi_bcast(hscr%gwgamma, master, comm, ierr)
1085  call xmpi_bcast(hscr%gwencomp, master, comm, ierr)
1086  call xmpi_bcast(hscr%kind_cdata, master, comm, ierr)
1087 ! HSCR_NEW
1088 
1089  DBG_EXIT("COLL")
1090 
1091 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

PARENTS

      m_io_screening,m_screening

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

1220 subroutine hscr_copy(Hscr_in,Hscr_cp)
1221 
1222 
1223 !This section has been created automatically by the script Abilint (TD).
1224 !Do not modify the following lines by hand.
1225 #undef ABI_FUNC
1226 #define ABI_FUNC 'hscr_copy'
1227 !End of the abilint section
1228 
1229  implicit none
1230 
1231 !Arguments ------------------------------------
1232 !scalars
1233  type(hscr_t),intent(in) :: Hscr_in
1234  type(hscr_t),intent(inout) :: Hscr_cp
1235 
1236 !Local variables-------------------------------
1237 !scalars
1238  !character(len=500) :: msg
1239 ! *************************************************************************
1240 
1241  !@hscr_t
1242  ! Integer values.
1243  Hscr_cp%id          = Hscr_in%id
1244  Hscr_cp%ikxc        = Hscr_in%ikxc
1245  Hscr_cp%inclvkb     = Hscr_in%inclvkb
1246  Hscr_cp%headform    = Hscr_in%headform
1247  Hscr_cp%fform       = Hscr_in%fform
1248  Hscr_cp%gwcalctyp   = Hscr_in%gwcalctyp
1249  Hscr_cp%nI          = Hscr_in%nI
1250  Hscr_cp%nJ          = Hscr_in%nJ
1251  Hscr_cp%nqibz       = Hscr_in%nqibz
1252  Hscr_cp%nqlwl       = Hscr_in%nqlwl
1253  Hscr_cp%nomega      = Hscr_in%nomega
1254  Hscr_cp%nbnds_used  = Hscr_in%nbnds_used
1255  Hscr_cp%npwe        = Hscr_in%npwe
1256  Hscr_cp%npwwfn_used = Hscr_in%npwwfn_used
1257  Hscr_cp%spmeth      = Hscr_in%spmeth
1258  Hscr_cp%test_type   = Hscr_in%test_type
1259  Hscr_cp%tordering   = Hscr_in%tordering
1260 
1261  ! Real variables
1262  Hscr_cp%mbpt_sciss = Hscr_in%mbpt_sciss
1263  Hscr_cp%spsmear  = Hscr_in%spsmear
1264  Hscr_cp%zcut     = Hscr_in%zcut
1265 
1266  ! Copy the abinit Header
1267  call hdr_copy(Hscr_in%Hdr,Hscr_cp%Hdr)
1268 
1269  Hscr_cp%titles(:) = Hscr_in%titles(:)
1270 
1271  ! Copy allocatable arrays.
1272  call alloc_copy(Hscr_in%gvec , Hscr_cp%gvec)
1273  call alloc_copy(Hscr_in%qibz , Hscr_cp%qibz)
1274  call alloc_copy(Hscr_in%qlwl , Hscr_cp%qlwl)
1275  call alloc_copy(Hscr_in%omega, Hscr_cp%omega)
1276 
1277 ! HSCR_NEW
1278  hscr_cp%awtr      =  hscr_in%awtr
1279  hscr_cp%icutcoul  =  hscr_in%icutcoul
1280  hscr_cp%vcutgeo   =  hscr_in%vcutgeo
1281  hscr_cp%gwcomp    =  hscr_in%gwcomp
1282  hscr_cp%gwgamma   =  hscr_in%gwgamma
1283  hscr_cp%gwencomp  =  hscr_in%gwencomp
1284  hscr_cp%kind_cdata  =  hscr_in%kind_cdata
1285 ! HSCR_NEW
1286 
1287 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)

PARENTS

      m_io_screening,m_screen,m_screening,mrgscr,screening,setup_bse

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

1161 subroutine hscr_free(hscr)
1162 
1163 
1164 !This section has been created automatically by the script Abilint (TD).
1165 !Do not modify the following lines by hand.
1166 #undef ABI_FUNC
1167 #define ABI_FUNC 'hscr_free'
1168 !End of the abilint section
1169 
1170  implicit none
1171 
1172 !Arguments ------------------------------------
1173 !scalars
1174  type(hscr_t),intent(inout) :: hscr
1175 
1176 ! *************************************************************************
1177 
1178  !@hscr_t
1179  DBG_ENTER("COLL")
1180 
1181  if (allocated(hscr%gvec)) then
1182    ABI_FREE(hscr%gvec)
1183  end if
1184  if (allocated(hscr%qibz)) then
1185    ABI_FREE(hscr%qibz)
1186  end if
1187  if (allocated(hscr%qlwl)) then
1188    ABI_FREE(hscr%qlwl)
1189  end if
1190  if (allocated(hscr%omega)) then
1191    ABI_FREE(hscr%omega)
1192  end if
1193 
1194  call hdr_free(hscr%Hdr)
1195 
1196  DBG_EXIT("COLL")
1197 
1198 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)

PARENTS

      m_io_screening,m_screen,m_screening,mrgscr,setup_bse

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

315 subroutine hscr_from_file(hscr,path,fform,comm)
316 
317 
318 !This section has been created automatically by the script Abilint (TD).
319 !Do not modify the following lines by hand.
320 #undef ABI_FUNC
321 #define ABI_FUNC 'hscr_from_file'
322 !End of the abilint section
323 
324  implicit none
325 
326 !Arguments ------------------------------------
327 !scalars
328  character(len=*),intent(in) :: path
329  integer,intent(in) :: comm
330  integer,intent(out) :: fform
331  type(hscr_t),intent(out) :: hscr
332 
333 !Local variables-------------------------------
334 !scalars
335  integer,parameter :: rdwr5=5,master=0
336  integer :: unt,my_rank,ierr
337  character(len=500) :: msg
338 
339 ! *************************************************************************
340 
341  my_rank = xmpi_comm_rank(comm)
342 
343  ! Master reads and broadcasts.
344  if (my_rank == master) then
345    if (.not. endswith(path, ".nc")) then
346      ! Fortran-IO
347      if (open_file(path,msg,newunit=unt,form="unformatted", status="old",action="read") /= 0) then
348        MSG_ERROR(msg)
349      end if
350      call hscr_io(hscr,fform,rdwr5,unt,xmpi_comm_self,master,IO_MODE_FORTRAN)
351      close(unt)
352    else
353      ! Netcdf format
354 #ifdef HAVE_NETCDF
355      NCF_CHECK(nctk_open_read(unt, path, xmpi_comm_self))
356      call hscr_io(hscr,fform,rdwr5,unt,xmpi_comm_self,master,IO_MODE_ETSF)
357      NCF_CHECK(nf90_close(unt))
358 #else
359      NETCDF_NOTENABLED_ERROR()
360 #endif
361    end if
362 
363    ABI_CHECK(fform /= 0, sjoin("hscr_io returned fform == 0 while reading:", path))
364  end if
365 
366  ! Broadcast data.
367  if (xmpi_comm_size(comm) > 1) then
368    call hscr_bcast(hscr,master,my_rank,comm)
369    call xmpi_bcast(fform,master,comm,ierr)
370  end if
371 
372 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.

PARENTS

      m_io_screening,m_screening,screening

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

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

PARENTS

      m_io_screening

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

1111 subroutine hscr_malloc(hscr, npwe, nqibz, nomega, nqlwl)
1112 
1113 
1114 !This section has been created automatically by the script Abilint (TD).
1115 !Do not modify the following lines by hand.
1116 #undef ABI_FUNC
1117 #define ABI_FUNC 'hscr_malloc'
1118 !End of the abilint section
1119 
1120  implicit none
1121 
1122 !Arguments ------------------------------------
1123 !scalars
1124  integer,intent(in) :: npwe, nqibz, nomega, nqlwl
1125  type(hscr_t),intent(inout) :: Hscr
1126 
1127 ! *************************************************************************
1128 
1129  !@hscr_t
1130  ABI_MALLOC(hscr%gvec, (3, npwe))
1131  ABI_MALLOC(hscr%qibz, (3, nqibz))
1132  ABI_MALLOC(hscr%qlwl, (3, nqlwl))
1133  ABI_MALLOC(hscr%omega, (nomega))
1134 
1135 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.

PARENTS

      m_io_screening

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

1313 subroutine hscr_merge(Hscr_in,Hscr_out)
1314 
1315 
1316 !This section has been created automatically by the script Abilint (TD).
1317 !Do not modify the following lines by hand.
1318 #undef ABI_FUNC
1319 #define ABI_FUNC 'hscr_merge'
1320 !End of the abilint section
1321 
1322  implicit none
1323 
1324 !Arguments ------------------------------------
1325 !scalars
1326  type(hscr_t),intent(in) :: Hscr_in(:)
1327  type(hscr_t),intent(out) :: Hscr_out
1328 
1329 !Local variables-------------------------------
1330 !scalars
1331  integer :: nhds,restart,restartpaw,ihd,ii,nqtot,nqneq
1332  logical :: isok
1333  character(len=500) :: msg
1334 !arrays
1335  real(dp),allocatable :: qset(:,:)
1336 ! *************************************************************************
1337 
1338  !@hscr_t
1339  nhds=SIZE(Hscr_in)
1340 
1341  ! Initial copy of the header ===
1342  ! If multiple headers, select the header containing q-->0 so that we copy also heads and wings
1343  ii = imax_loc(Hscr_in(:)%nqlwl)
1344  call hscr_copy(Hscr_in(ii),Hscr_out)
1345  if (nhds==1) return
1346 
1347  ! Check consistency of the abinit Headers.
1348  ! FFT grid might be q-point dependent so we stop only when restart==0
1349  isok=.TRUE.
1350  do ihd=2,nhds
1351    call hdr_check(Hscr_in(1)%fform,Hscr_in(ihd)%fform,Hscr_in(1)%Hdr,Hscr_in(ihd)%Hdr,'COLL',restart,restartpaw)
1352    if (restart==0) then
1353      isok=.FALSE.
1354      write(msg,'(a,i0,a)')' Abinit header no.',ihd,' is not consistent with the first header '
1355      MSG_WARNING(msg)
1356    end if
1357  end do
1358  if (.not.isok) then
1359    MSG_ERROR('Cannot continue, Check headers')
1360  end if
1361 
1362  ! Now check variables related to polarizability|epsilon^{-1}.
1363  ! 1) Tests quantities that must be equal
1364  ii = assert_eq(Hscr_in(:)%ID,       'Headers have different Identifiers')
1365  ii = assert_eq(Hscr_in(:)%ikxc,     'Headers have different ikxc'       )
1366  ii = assert_eq(Hscr_in(:)%headform, 'Headers have different headform'   )
1367  ii = assert_eq(Hscr_in(:)%fform,    'Headers have different fform'      )
1368  ii = assert_eq(Hscr_in(:)%gwcalctyp,'Headers have different gwcalctyp'  )
1369  ii = assert_eq(Hscr_in(:)%nI,       'Headers have different nI'         )
1370  ii = assert_eq(Hscr_in(:)%nJ,       'Headers have different nJ'         )
1371  ii = assert_eq(Hscr_in(:)%nomega,   'Headers have different nomega'     )
1372  ii = assert_eq(Hscr_in(:)%test_type,'Headers have different test_type'  )
1373  ii = assert_eq(Hscr_in(:)%tordering,'Headers have different tordering'  )
1374 
1375  ! This is not mandatory but makes life easier!
1376  ii = assert_eq(Hscr_in(:)%npwe,'Headers have different number of G-vectors'  )
1377 
1378  do ihd=2,nhds
1379    if (ANY(ABS(Hscr_in(ihd)%omega-Hscr_in(1)%omega)>tol6)) then
1380      write(msg,'(a,i0,a)')' Frequencies in the first and the ',ihd,'-th header differ'
1381      MSG_ERROR(msg)
1382    end if
1383    if (ANY(Hscr_in(ihd)%gvec(:,:)-Hscr_in(1)%gvec(:,:)/=0)) then
1384      write(msg,'(a,i0,a)')' Incompatible G-vector list found in the ',ihd,'-th header'
1385      MSG_ERROR(msg)
1386    end if
1387    if (hscr_in(ihd)%kind_cdata /= hscr_in(1)%kind_cdata) then
1388      write(msg,'(3a,i0,2a)')' Files contain data with different precisions.',ch10,&
1389      "In particular the ",ihd,'-th header has precision:',trim(hscr_in(ihd)%kind_cdata)
1390      MSG_ERROR(msg)
1391    end if
1392  end do !ihd
1393 
1394  ! If error is not fatal, just warn ===
1395  if (ANY(Hscr_in(:)%npwwfn_used/=Hscr_in(1)%npwwfn_used)) then
1396    MSG_COMMENT('Files have been produced with a different number of planewaves for the wavefunctions.')
1397  end if
1398  if (ANY(Hscr_in(:)%nbnds_used/=Hscr_in(1)%nbnds_used)) then
1399    MSG_COMMENT('Files have been produced with a different number of bands.')
1400  end if
1401  if (ANY(Hscr_in(:)%spmeth/=Hscr_in(1)%spmeth)) then
1402    MSG_COMMENT('Files have been produced with different algorithms.')
1403  end if
1404  if (ANY(ABS(Hscr_in(:)%mbpt_sciss-Hscr_in(1)%mbpt_sciss)>tol6)) then
1405    MSG_COMMENT('Files have benn produced with different values of mbpt_sciss.')
1406  end if
1407  if (ANY(ABS(Hscr_in(:)%spsmear-Hscr_in(1)%spsmear)>tol6)) then
1408    MSG_COMMENT('Files have been produced with different values of spsmear.')
1409  end if
1410  if (ANY(ABS(Hscr_in(:)%zcut-Hscr_in(1)%zcut)>tol6)) then
1411    MSG_COMMENT('Files have been produced with different values of zcut.')
1412  end if
1413 
1414  ! Now merge the list of q-points.
1415  ! Take the union of the q-points, remove possible duplicated
1416  ! are change the parameters in hscr_out that depends on q-points.
1417  nqtot=SUM(Hscr_in(:)%nqibz)
1418  ABI_MALLOC(qset,(3,nqtot))
1419 
1420  ii=0
1421  do ihd=1,nhds
1422    qset(:,ii+1:ii+Hscr_in(ihd)%nqibz)=Hscr_in(ihd)%qibz(:,:)
1423    ii=ii+Hscr_in(ihd)%nqibz
1424  end do
1425 
1426  call remove_copies(nqtot,qset,nqneq,isequalk)
1427 
1428  if (nqneq /= nqtot) then
1429    write(msg,'(3a,2(i0,a))')&
1430     'COMMENT: Headers contain duplicated q-points ',ch10,&
1431     'Found ',nqneq,' distinct q-points among the total ',nqtot,' points reported in the headers. '
1432    call wrtout(std_out, msg)
1433  end if
1434 
1435  Hscr_out%nqibz = nqneq
1436  ABI_FREE(Hscr_out%qibz)
1437  ABI_MALLOC(Hscr_out%qibz,(3,nqneq))
1438  Hscr_out%qibz(:,:)=qset(:,1:nqneq)
1439  ABI_FREE(qset)
1440 
1441 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.

PARENTS

      m_io_screening

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

1934 subroutine hscr_mpio_skip(mpio_fh,fform,offset)
1935 
1936 
1937 !This section has been created automatically by the script Abilint (TD).
1938 !Do not modify the following lines by hand.
1939 #undef ABI_FUNC
1940 #define ABI_FUNC 'hscr_mpio_skip'
1941 !End of the abilint section
1942 
1943  implicit none
1944 
1945 !Arguments ------------------------------------
1946  integer,intent(in) :: mpio_fh
1947  integer,intent(out) :: fform
1948  integer(kind=XMPI_OFFSET_KIND),intent(out) :: offset
1949 
1950 !Local variables-------------------------------
1951 !scalars
1952  integer :: bsize_frm,mpi_type_frm
1953 #ifdef HAVE_MPI_IO
1954  integer :: ierr,isk,nqlwl
1955  !character(len=500) :: msg
1956 !arrays
1957  integer(kind=MPI_OFFSET_KIND) :: fmarker,positloc
1958  integer :: statux(MPI_STATUS_SIZE)
1959 #endif
1960 
1961 ! *************************************************************************
1962 
1963  offset = 0
1964  bsize_frm    = xmpio_bsize_frm    ! Byte size of the Fortran record marker.
1965  mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker.
1966 
1967  call hdr_mpio_skip(mpio_fh,fform,offset)
1968 
1969  call wrtout(std_out, sjoin("in hdr_mpio_skip with fform = ",itoa(fform)))
1970 
1971 #ifdef HAVE_MPI_IO
1972  select case (fform)
1973  case (1003, 1004)
1974    ! Ship the titles
1975    call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1976 
1977    ! read nqlwl from the 2d record.
1978    positloc  = offset + bsize_frm + 9*xmpi_bsize_int
1979    call MPI_FILE_READ_AT(mpio_fh,positloc,nqlwl,1,MPI_INTEGER,statux,ierr)
1980    call wrtout(std_out, sjoin("nqlwl = ",itoa(nqlwl)))
1981 
1982    do isk=1,5
1983      call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1984    end do
1985 
1986    if (nqlwl>0) then  ! skip qlwl
1987      call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1988    end if
1989 
1990  case default
1991    MSG_BUG(sjoin('Wrong fform read:', itoa(fform)))
1992  end select
1993 
1994 #else
1995  MSG_ERROR("hscr_mpio_skip cannot be used when MPI-IO is not enabled")
1996 #endif
1997 
1998 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.

PARENTS

      screening

CHILDREN

SOURCE

893 type(hscr_t) function hscr_new(varname,dtset,ep,hdr_abinit,ikxc,test_type,tordering,titles,ngvec,gvec) result(hscr)
894 
895 
896 !This section has been created automatically by the script Abilint (TD).
897 !Do not modify the following lines by hand.
898 #undef ABI_FUNC
899 #define ABI_FUNC 'hscr_new'
900 !End of the abilint section
901 
902  implicit none
903 
904 !Arguments ------------------------------------
905 !scalars
906  integer,intent(in) :: ikxc,test_type,tordering,ngvec
907  character(len=*),intent(in) :: varname
908  type(dataset_type),intent(in) :: dtset
909  type(em1params_t),intent(in) :: Ep
910  type(hdr_type),intent(in) :: hdr_abinit
911 !arrays
912  integer,intent(in) :: gvec(3,ngvec)
913  character(len=80),intent(in) :: titles(2)
914 
915 !Local variables-------------------------------
916  integer :: id
917  type(abifile_t) :: abifile
918 ! *************************************************************************
919 
920  !@hscr_t
921  ABI_CHECK(ngvec==Ep%npwe,'ngvec/=Ep%npwe')
922 
923  ! ID=Identifier used to define the type of Response function (e^-1, chi0)
924  id = 0
925  if (varname == "polarizability") id = 1
926  if (varname == "inverse_dielectric_function") id = 4
927  ABI_CHECK(id /= 0, sjoin("Invalid varname: ",varname))
928 
929  ! Get fform from abifile.
930  abifile = abifile_from_varname(varname)
931  if (abifile%fform == 0) then
932     MSG_ERROR(sjoin("Cannot find any abifile object associated to varname:", varname))
933  end if
934 
935  ! Copy the abinit header.
936  call hdr_copy(hdr_abinit,Hscr%Hdr)
937 
938  ! Initialize quantities related to the screening file
939  hscr%id         =id
940  hscr%ikxc       =ikxc
941  hscr%inclvkb    =Ep%inclvkb
942  hscr%headform   =HSCR_LATEST_HEADFORM
943  hscr%fform      =abifile%fform
944  hscr%gwcalctyp  =Ep%gwcalctyp
945  hscr%nI         =Ep%nI
946  hscr%nJ         =Ep%nJ
947  hscr%nqibz      =Ep%nqcalc  ! nqcalc == nqibz except if we split the calculation with nqptdm
948  hscr%nqlwl      =Ep%nqlwl
949  hscr%nomega     =Ep%nomega
950  hscr%nbnds_used =Ep%nbnds
951  hscr%npwe       =Ep%npwe
952  hscr%npwwfn_used=Ep%npwwfn
953  hscr%spmeth     =Ep%spmeth
954  hscr%test_type  =test_type
955  hscr%tordering  =tordering
956  hscr%mbpt_sciss =Ep%mbpt_sciss
957  hscr%spsmear    =Ep%spsmear
958  hscr%zcut       =Ep%zcut
959 
960  hscr%titles(:)=titles(:)
961 
962  call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
963  hscr%gvec(:,:) = gvec(1:3,1:Ep%npwe)
964  hscr%qibz(:,:) = Ep%qcalc
965  hscr%qlwl(:,:) = Ep%qlwl
966  hscr%omega(:) = Ep%omega
967 
968 ! HSCR_NEW
969  hscr%awtr = dtset%awtr
970  hscr%icutcoul = dtset%icutcoul
971  hscr%vcutgeo = dtset%vcutgeo
972  hscr%gwcomp = dtset%gwcomp
973  hscr%gwgamma = dtset%gwgamma
974  hscr%gwencomp = dtset%gwencomp
975  hscr%kind_cdata = "dpc" ! For the time being, data is always written in double-precision
976 ! HSCR_NEW
977 
978 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

PARENTS

      m_io_screening,m_screen,m_screening,mrgscr,setup_bse

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

759 subroutine hscr_print(Hscr,header,unit,prtvol,mode_paral)
760 
761 
762 !This section has been created automatically by the script Abilint (TD).
763 !Do not modify the following lines by hand.
764 #undef ABI_FUNC
765 #define ABI_FUNC 'hscr_print'
766 !End of the abilint section
767 
768  implicit none
769 
770 !Arguments ------------------------------------
771 !scalars
772  integer,intent(in),optional :: prtvol,unit
773  character(len=4),intent(in),optional :: mode_paral
774  character(len=*),intent(in),optional :: header
775  type(hscr_t),intent(in) :: hscr
776 
777 !Local variables-------------------------------
778 !scalars
779  integer :: iomega,iqibz,unt,verbose
780  character(len=4) :: mode
781  character(len=500) :: msg
782 
783 ! *************************************************************************
784 
785  unt=std_out; if (PRESENT(unit      )) unt    =unit
786  verbose=0  ; if (PRESENT(prtvol    )) verbose=prtvol
787  mode='COLL'; if (PRESENT(mode_paral)) mode   =mode_paral
788 
789  if (PRESENT(header)) then
790    msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
791    call wrtout(unt,msg,mode)
792  end if
793 
794  write(msg,'(1x,a)')TRIM(hscr%titles(1))
795  call wrtout(unt,msg,mode)
796  write(msg,'(1x,a)')TRIM(hscr%titles(2))
797  call wrtout(unt,msg,mode)
798  write(msg,'(a,i8)') ' Identifier                ',hscr%ID
799  call wrtout(unt,msg,mode)
800  write(msg,'(a,i8)') ' Kxc kernel                ',hscr%ikxc
801  call wrtout(unt,msg,mode)
802  write(msg,'(a,i8)') ' Treatment of q-->0 limit  ',hscr%inclvkb
803  call wrtout(unt,msg,mode)
804  write(msg,'(a,i8)') ' headform                  ',hscr%headform
805  call wrtout(unt,msg,mode)
806  write(msg,'(a,i8)') ' fform                     ',hscr%fform
807  call wrtout(unt,msg,mode)
808  write(msg,'(a,i8)') ' gwcalctyp                 ',hscr%gwcalctyp
809  call wrtout(unt,msg,mode)
810  write(msg,'(a,2i8)')' Number of components      ',hscr%nI,hscr%nJ
811  call wrtout(unt,msg,mode)
812  write(msg,'(a,i8)') ' Number of q-points        ',hscr%nqibz
813  call wrtout(unt,msg,mode)
814  write(msg,'(a,i8)') ' Number of q-directions    ',hscr%nqlwl
815  call wrtout(unt,msg,mode)
816  write(msg,'(a,i8)') ' Number of frequencies     ',hscr%nomega
817  call wrtout(unt,msg,mode)
818  write(msg,'(a,i8)') ' Number of bands used      ',hscr%nbnds_used
819  call wrtout(unt,msg,mode)
820  write(msg,'(a,i8)') ' Dimension of matrix       ',hscr%npwe
821  call wrtout(unt,msg,mode)
822  write(msg,'(a,i8)') ' Number of planewaves used ',hscr%npwwfn_used
823  call wrtout(unt,msg,mode)
824  write(msg,'(a,i8)') ' Spectral method           ',hscr%spmeth
825  call wrtout(unt,msg,mode)
826  write(msg,'(a,i8)') ' Test_type                 ',hscr%test_type
827  call wrtout(unt,msg,mode)
828  write(msg,'(a,i8)') ' Time-ordering             ',hscr%tordering
829  call wrtout(unt,msg,mode)
830  write(msg,'(a,es16.6)')' Scissor Energy             ',hscr%mbpt_sciss
831  call wrtout(unt,msg,mode)
832  write(msg,'(a,es16.6)')' Spectral smearing          ',hscr%spsmear
833  call wrtout(unt,msg,mode)
834  write(msg,'(a,es16.6)')' Complex Imaginary Shift    ',hscr%zcut
835  call wrtout(unt,msg,mode)
836 
837  if (verbose==0) then
838    call wrtout(unt,' The header contains additional records.',mode)
839  else
840    write(msg,'(2a)')ch10,' q-points [r.l.u.]:'
841    call wrtout(unt,msg,mode)
842    do iqibz=1,hscr%nqibz
843      write(msg,'(i5,3f12.6)')iqibz,hscr%qibz(:,iqibz)
844      call wrtout(unt,msg,mode)
845    end do
846 
847    write(msg,'(2a)')ch10,' Frequencies used [eV]:'
848    call wrtout(unt,msg,mode)
849    do iomega=1,hscr%nomega
850      write(msg,'(i3,2f7.2)')iomega,REAL(hscr%omega(iomega))*Ha_eV,AIMAG(hscr%omega(iomega))*Ha_eV
851      call wrtout(unt,msg,mode)
852    end do
853  end if
854 
855 ! HSCR_NEW
856 ! HSCR_NEW
857 
858  ! Echo the abinit header.
859  !if (prtvol>0) call hdr_echo(hscr%hdr,fform,rdwr,unit=unt)
860 
861 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

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

PARENTS

      mrgscr

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

2028 subroutine ioscr_qmerge(nfiles, filenames, hscr_files, fname_out, ohscr)
2029 
2030 
2031 !This section has been created automatically by the script Abilint (TD).
2032 !Do not modify the following lines by hand.
2033 #undef ABI_FUNC
2034 #define ABI_FUNC 'ioscr_qmerge'
2035 !End of the abilint section
2036 
2037  implicit none
2038 
2039 !Arguments ------------------------------------
2040 !scalars
2041  integer,intent(in) :: nfiles
2042  character(len=*),intent(in) :: fname_out
2043  type(hscr_t),intent(out) :: ohscr
2044 !arrays
2045  character(len=*),intent(in) :: filenames(nfiles)
2046  type(hscr_t),intent(in) :: hscr_files(nfiles)
2047 
2048 !Local variables-------------------------------
2049 !scalars
2050  integer,parameter :: rdwr2=2,master=0
2051  integer :: iqibz,ifound,ifile,iqf,ount,iomode,fform_merge,comm,nomega4m,npwe4m,iqiA,ierr
2052  character(len=500) :: msg
2053  character(len=nctk_slen) :: varname
2054  type(abifile_t) :: abifile
2055 !arrays
2056  integer,allocatable :: merge_table(:,:)
2057  real(dp) :: qdiff(3)
2058  complex(gwpc),allocatable :: epsm1(:,:,:,:)
2059 
2060 ! *************************************************************************
2061 
2062  comm = xmpi_comm_self
2063 
2064  if (file_exists(fname_out)) then
2065    MSG_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
2066  end if
2067 
2068  ! Merge the headers creating the full list of q-points.
2069  call hscr_merge(Hscr_files(1:nfiles), ohscr)
2070  call hscr_print(ohscr, header='Header of the final file', unit=std_out, prtvol=1)
2071 
2072  ! For each q to be merged, save the index of the file where q is stored as well as its sequential index.
2073  ! Useful to do the merge point-by-point thus avoiding the allocation of the entire epsm1 array.
2074  ABI_MALLOC(merge_table,(ohscr%nqibz,2))
2075  do iqibz=1,ohscr%nqibz
2076    ifound=0
2077    fl: do ifile=1,nfiles
2078      do iqf=1,Hscr_files(ifile)%nqibz
2079        qdiff(:)=ohscr%qibz(:,iqibz)-Hscr_files(ifile)%qibz(:,iqf)
2080        if (all(abs(qdiff) < GW_TOLQ)) then
2081          merge_table(iqibz,1)=ifile
2082          merge_table(iqibz,2)=iqf
2083          ifound=ifound+1
2084          write(msg,'(a,3f12.6,2a)')'. q-point:',ohscr%qibz(:,iqibz),' will be taken from ',TRIM(filenames(ifile))
2085          call wrtout(std_out,msg,'COLL')
2086          EXIT fl
2087        end if
2088      end do
2089    end do fl
2090    ! Check if q-point has been found, multiple q-points not allowed.
2091    ABI_CHECK(ifound == 1, 'ifound/=1')
2092  end do
2093 
2094  iomode = IO_MODE_FORTRAN; if (endswith(fname_out, ".nc")) iomode = IO_MODE_ETSF
2095  if (iomode == IO_MODE_FORTRAN) then
2096    if (open_file(fname_out,msg,newunit=ount,status='new',form='unformatted') /= 0) then
2097      MSG_ERROR(msg)
2098    end if
2099  else
2100 #ifdef HAVE_NETCDF
2101    NCF_CHECK(nctk_open_create(ount, fname_out, comm))
2102 #endif
2103  end if
2104 
2105  ! Write the header.
2106  fform_merge = hscr_files(1)%fform
2107  abifile = abifile_from_fform(fform_merge)
2108  if (abifile%fform == 0) then
2109    MSG_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
2110  end if
2111  varname = abifile%varname
2112 
2113  if (any(hscr_files(:)%fform /= hscr_files(1)%fform)) then
2114    write(std_out,*)"fforms: ",hscr_files(:)%fform
2115    MSG_ERROR("Files to be merged have different fform. Cannot merge data")
2116  end if
2117 
2118  call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode)
2119 
2120  npwe4m   = ohscr%npwe
2121  nomega4m = ohscr%nomega
2122 
2123  ABI_STAT_MALLOC(epsm1,(npwe4m,npwe4m,nomega4m,1), ierr)
2124  ABI_CHECK(ierr==0, 'out of memory in epsm1')
2125 
2126  do iqibz=1,ohscr%nqibz
2127    ifile = merge_table(iqibz,1)
2128    iqiA  = merge_table(iqibz,2)
2129    call read_screening(varname,filenames(ifile),npwe4m,1,nomega4m,epsm1,iomode,comm,iqiA=iqiA)
2130    call write_screening(varname,ount,iomode,npwe4m,nomega4m,iqibz,epsm1)
2131  end do
2132 
2133  ABI_FREE(epsm1)
2134  ABI_FREE(merge_table)
2135 
2136  if (iomode == IO_MODE_FORTRAN) then
2137    close(ount)
2138  else
2139 #ifdef HAVE_NETCDF
2140    NCF_CHECK(nf90_close(ount))
2141 #endif
2142  end if
2143 
2144  write(msg,'(3a)')ch10,' ==== Files have been merged successfully === ',ch10
2145  call wrtout(std_out,msg,'COLL')
2146 
2147 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.

PARENTS

      mrgscr

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

2176 subroutine ioscr_qrecover(ipath, nqrec, fname_out)
2177 
2178 
2179 !This section has been created automatically by the script Abilint (TD).
2180 !Do not modify the following lines by hand.
2181 #undef ABI_FUNC
2182 #define ABI_FUNC 'ioscr_qrecover'
2183 !End of the abilint section
2184 
2185  implicit none
2186 
2187 !Arguments ------------------------------------
2188 !scalars
2189  integer,intent(in) :: nqrec
2190  character(len=*),intent(in) :: ipath,fname_out
2191 
2192 !Local variables-------------------------------
2193 !scalars
2194  integer,parameter :: rdwr2=2,master=0
2195  integer :: iqiA,nqibzA,nomega_asked,unt,npwe_asked,iomode,comm,fform1,ifform,ierr
2196  character(len=500) :: msg
2197  character(len=nctk_slen) :: varname
2198  type(hscr_t) :: hscr_recov,hscr
2199  type(abifile_t) :: abifile
2200 !arrays
2201  complex(gwpc),allocatable :: epsm1(:,:,:,:)
2202 
2203 ! *************************************************************************
2204 
2205  comm = xmpi_comm_self
2206 
2207  call wrtout(std_out, sjoin(". Recovering q-points in file:", ipath))
2208  call wrtout(std_out, sjoin(". Data written to file:", fname_out))
2209 
2210  if (file_exists(fname_out)) then
2211    MSG_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
2212  end if
2213 
2214  ! Find iomode from file extension and open output file.
2215  if (endswith(fname_out, ".nc")) then
2216    iomode = IO_MODE_ETSF
2217 #ifdef HAVE_NETCDF
2218    NCF_CHECK(nctk_open_create(unt, fname_out, comm))
2219 #else
2220    MSG_ERROR("Netcdf support is not available")
2221 #endif
2222  else
2223    iomode = IO_MODE_FORTRAN
2224    if (open_file(fname_out, msg, newunit=unt, status='new', form='unformatted') /= 0) then
2225      MSG_ERROR(msg)
2226    end if
2227  end if
2228 
2229  ! Read header.
2230  call hscr_from_file(hscr, ipath, ifform, comm)
2231  ABI_CHECK(ifform /= 0, sjoin("fform = 0 while reading:", ipath))
2232 
2233  if (nqrec < 1 .or. nqrec > hscr%nqibz) then
2234    MSG_ERROR(sjoin("Wrong input. nqibz on file:", itoa(hscr%nqibz)))
2235  end if
2236 
2237  ! Copy header
2238  call hscr_copy(hscr, hscr_recov)
2239 
2240  ! Change dimensions and arrays associated to nqibz.
2241  hscr_recov%nqibz = nqrec
2242  ABI_FREE(hscr_recov%qibz)
2243  ABI_MALLOC(hscr_recov%qibz, (3,nqrec))
2244  hscr_recov%qibz = hscr%qibz(:,1:nqrec)
2245 
2246  call hscr_print(hscr_recov,header="Header of the new SCR file",unit=std_out,prtvol=1)
2247 
2248  ! Write the header of the recovered file.
2249  fform1 = hscr%fform
2250 
2251  abifile = abifile_from_fform(fform1)
2252  if (abifile%fform == 0) then
2253     MSG_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform1)))
2254  end if
2255  varname = abifile%varname
2256 
2257  call hscr_io(hscr_recov,fform1,rdwr2,unt,comm,master,iomode)
2258 
2259  nqibzA=1; nomega_asked=hscr%nomega; npwe_asked=hscr%npwe
2260 
2261  ABI_STAT_MALLOC(epsm1,(npwe_asked,npwe_asked,nomega_asked,1), ierr)
2262  ABI_CHECK(ierr==0, 'out of memory in epsm1')
2263 
2264  do iqiA=1,hscr_recov%nqibz
2265    call read_screening(varname,ipath,npwe_asked,nqibzA,nomega_asked,epsm1,iomode,comm,iqiA=iqiA)
2266    call write_screening(varname,unt,iomode,npwe_asked,nomega_asked,iqiA,epsm1)
2267  end do
2268 
2269  if (iomode == IO_MODE_FORTRAN) close(unt)
2270 #ifdef HAVE_NETCDF
2271  if (iomode == IO_MODE_ETSF) then
2272    NCF_CHECK(nf90_close(unt))
2273  end if
2274 #endif
2275 
2276  ABI_FREE(epsm1)
2277  call hscr_free(hscr)
2278  call hscr_free(hscr_recov)
2279 
2280  call wrtout(std_out,"Recovery completed",'COLL')
2281 
2282 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.

PARENTS

      mrgscr

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

2312 subroutine ioscr_wmerge(nfiles, filenames, hscr_file, freqremax, fname_out, ohscr)
2313 
2314 
2315 !This section has been created automatically by the script Abilint (TD).
2316 !Do not modify the following lines by hand.
2317 #undef ABI_FUNC
2318 #define ABI_FUNC 'ioscr_wmerge'
2319 !End of the abilint section
2320 
2321  implicit none
2322 
2323 !Arguments ------------------------------------
2324 !scalars
2325  integer,intent(in) :: nfiles
2326  real(dp),intent(in) :: freqremax
2327  character(len=*),intent(in) :: fname_out
2328  type(hscr_t),intent(out) :: ohscr
2329 !arrays
2330  character(len=*),intent(in) :: filenames(nfiles)
2331  type(hscr_t),intent(in) :: hscr_file(nfiles)
2332 
2333 !Local variables-------------------------------
2334 !scalars
2335  integer,parameter :: rdwr2=2,master=0
2336  integer :: ii,iqibz,ifile,ount,iomode,fform_merge,comm,nomega4m
2337  integer :: nfreq_tot,nfreqre,nfreqim,ifrq,npwe4mI,npwe4mJ,ierr
2338  character(len=500) :: msg
2339  logical :: skip
2340  character(len=nctk_slen) :: varname
2341  type(abifile_t) :: abifile
2342 !arrays
2343  integer,allocatable :: freq_indx(:,:),ifile_indx(:),pos_indx(:),i_temp(:),i2_temp(:,:)
2344  real(dp),allocatable :: real_omega(:),imag_omega(:)
2345  complex(gwpc),allocatable :: epsm1(:,:,:,:),epsm1_temp(:,:,:,:)
2346  complex(dpc),allocatable :: omega_storage(:)
2347 
2348 ! *************************************************************************
2349 
2350  comm = xmpi_comm_self
2351 
2352  ! Check that q-point sets are the same
2353  do ifile=1,nfiles
2354    do ii=1,nfiles
2355      if (ii == ifile) CYCLE
2356      if (Hscr_file(ifile)%nqibz /= Hscr_file(ii)%nqibz) then
2357        MSG_ERROR(' One or more files do not have the same number of q-points!')
2358      end if
2359      do iqibz=1,Hscr_file(1)%nqibz
2360        if (ABS(SUM(Hscr_file(ifile)%qibz(:,iqibz)-Hscr_file(ii)%qibz(:,iqibz))) > tol6) then
2361          MSG_ERROR('Q-point set differs between one or more files!')
2362        end if
2363      end do
2364    end do
2365  end do
2366 
2367  ! nfreq_tot here is the total *possible* number of freq.
2368  nfreq_tot=0
2369  do ifile=1,nfiles
2370    nfreq_tot = nfreq_tot + Hscr_file(ifile)%nomega
2371  end do
2372 
2373  ABI_MALLOC(omega_storage, (nfreq_tot))
2374  ABI_MALLOC(freq_indx, (nfreq_tot,nfiles))
2375  ABI_MALLOC(ifile_indx, (nfreq_tot))
2376  omega_storage = CMPLX(-one,-one); freq_indx = 0; ifile_indx = 0
2377 
2378  ! Calculate the total number of real freq and store
2379  nfreqre = 0
2380  do ifile=1,nfiles
2381    do ifrq=1,Hscr_file(ifile)%nomega
2382      skip = .FALSE.
2383      ! Check whether to skip this point
2384      if (AIMAG(Hscr_file(ifile)%omega(ifrq)) > tol16) skip = .TRUE.
2385      if (REAL(Hscr_file(ifile)%omega(ifrq)) > freqremax) skip = .TRUE.
2386      ! Check for repetition or non-monotonic points
2387      if (nfreqre>1) then
2388        do ii=1,nfreqre
2389          if (ABS(REAL(Hscr_file(ifile)%omega(ifrq)) - REAL(omega_storage(ii))) < tol6) skip = .TRUE.
2390        end do
2391      end if
2392      if (skip) CYCLE
2393      nfreqre = nfreqre + 1
2394      ! Store (complex) frequency and index
2395      omega_storage(nfreqre) = Hscr_file(ifile)%omega(ifrq)
2396      ifile_indx(nfreqre) = ifile
2397      freq_indx(nfreqre,ifile) = ifrq
2398      write(std_out,'(a,es16.6,a,i0,2(a,i0))')&
2399        ' Found real frequency: ',REAL(omega_storage(nfreqre))*Ha_eV,' [eV], number: ',nfreqre,&
2400        ', in file: ',ifile,' local index: ',ifrq
2401    end do
2402  end do
2403 
2404  if (nfreqre > 0) then
2405    ! Sort real frequencies
2406    ABI_MALLOC(real_omega,(nfreqre))
2407    ABI_MALLOC(pos_indx,(nfreqre))
2408    ABI_MALLOC(i_temp,(nfreqre))
2409    ABI_MALLOC(i2_temp,(nfreqre,nfiles))
2410    real_omega(1:nfreqre) = REAL(omega_storage(1:nfreqre)) ! Copy real frequencies to temp. sorting array
2411    do ii=1,nfreqre ! Set up indexing array
2412      pos_indx(ii) = ii
2413    end do
2414 
2415    ! Sort frequencies while keeping track of index
2416    call sort_dp(nfreqre,real_omega,pos_indx,tol16)
2417    i_temp(1:nfreqre) = ifile_indx(1:nfreqre)
2418    i2_temp(1:nfreqre,1:nfiles) = freq_indx(1:nfreqre,1:nfiles)
2419 
2420    ! Copy sorted frequencies plus file and frequency index
2421    do ii=1,nfreqre
2422      omega_storage(ii) = CMPLX(real_omega(ii),zero)
2423      ifile_indx(ii) = i_temp(pos_indx(ii))
2424      freq_indx(ii,1:nfiles) = i2_temp(pos_indx(ii),1:nfiles)
2425    end do
2426 
2427    ABI_FREE(real_omega)
2428    ABI_FREE(pos_indx)
2429    ABI_FREE(i_temp)
2430    ABI_FREE(i2_temp)
2431  end if
2432 
2433  ! Check imaginary frequencies and store them
2434  nfreqim = 0
2435  do ifile=1,nfiles
2436    do ifrq=1,Hscr_file(ifile)%nomega
2437      if (REAL(Hscr_file(ifile)%omega(ifrq)) > tol8) CYCLE
2438      if (AIMAG(Hscr_file(ifile)%omega(ifrq)) < tol8) CYCLE
2439      nfreqim = nfreqim + 1
2440      omega_storage(nfreqre+nfreqim) = Hscr_file(ifile)%omega(ifrq)
2441      ifile_indx(nfreqre+nfreqim) = ifile
2442      freq_indx(nfreqre+nfreqim,ifile) = ifrq
2443      write(std_out,'(a,es16.6,a,i0,2(a,i0))')&
2444       ' Found imaginary frequency: ',AIMAG(omega_storage(nfreqre+nfreqim))*Ha_eV,' [eV], number: ',nfreqim,&
2445       ', in file: ',ifile,' local index: ',ifrq
2446    end do
2447  end do
2448 
2449  ! Sort imaginary frequencies
2450  ABI_MALLOC(imag_omega,(nfreqim))
2451  ABI_MALLOC(pos_indx,(nfreqim))
2452  ABI_MALLOC(i_temp,(nfreqim))
2453  ABI_MALLOC(i2_temp,(nfreqim,nfiles))
2454 
2455  ! Copy imaginary frequencies to temp. sorting array
2456  imag_omega(1:nfreqim) = AIMAG(omega_storage(nfreqre+1:nfreqre+nfreqim))
2457  do ii=1,nfreqim ! Set up indexing array
2458    pos_indx(ii) = ii
2459  end do
2460 
2461  ! Sort frequencies while keeping track of index
2462  call sort_dp(nfreqim,imag_omega,pos_indx,tol16)
2463  i_temp(1:nfreqim) = ifile_indx(nfreqre+1:nfreqre+nfreqim)
2464  i2_temp(1:nfreqim,1:nfiles) = freq_indx(nfreqre+1:nfreqre+nfreqim,1:nfiles)
2465 
2466  ! Copy sorted frequencies plus file and frequency index
2467  do ii=1,nfreqim
2468    omega_storage(nfreqre+ii) = CMPLX(zero,imag_omega(ii))
2469    ifile_indx(nfreqre+ii) = i_temp(pos_indx(ii))
2470    freq_indx(nfreqre+ii,1:nfiles) = i2_temp(pos_indx(ii),1:nfiles)
2471  end do
2472 
2473  ABI_FREE(imag_omega)
2474  ABI_FREE(pos_indx)
2475  ABI_FREE(i_temp)
2476  ABI_FREE(i2_temp)
2477 
2478  nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq
2479  write(std_out,'(2a,i0,a)') ch10,' Merging ',nfreq_tot,' frequencies.'
2480  write(std_out,'(2(a,i0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10
2481 
2482  ! Copy old header
2483  call hscr_copy(Hscr_file(1),ohscr)
2484 
2485  ! TODO: hscr_wmerge
2486  ! Then modify entries for new frequency grid
2487  ohscr%nomega = nfreq_tot
2488  ABI_FREE(ohscr%omega)
2489  ABI_MALLOC(ohscr%omega,(nfreq_tot))
2490  ohscr%omega(1:nfreq_tot) = omega_storage(1:nfreq_tot)
2491 
2492  npwe4mI = ohscr%npwe*ohscr%nI
2493  npwe4mJ = ohscr%npwe*ohscr%nJ
2494 
2495  ! Print new header for info
2496  call hscr_print(ohscr,header='Header of the final file',unit=std_out,prtvol=1)
2497 
2498  if (file_exists(fname_out)) then
2499    MSG_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
2500  end if
2501 
2502  if (endswith(fname_out, ".nc")) then
2503    iomode = IO_MODE_ETSF
2504 #ifdef HAVE_NETCDF
2505    NCF_CHECK(nctk_open_create(ount, fname_out, comm))
2506 #else
2507    MSG_ERROR("netcdf support is not activated")
2508 #endif
2509  else
2510    iomode = IO_MODE_FORTRAN
2511    if (open_file(fname_out, msg, newunit=ount, status='new',form='unformatted') /= 0) then
2512      MSG_ERROR(msg)
2513    end if
2514  end if
2515 
2516  ! * Write the header.
2517  fform_merge = ohscr%fform
2518 
2519  abifile = abifile_from_fform(fform_merge)
2520  if (abifile%fform == 0) then
2521     MSG_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
2522  end if
2523  varname = abifile%varname
2524 
2525  call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode)
2526 
2527  npwe4mI = ohscr%npwe*ohscr%nI
2528  npwe4mJ = ohscr%npwe*ohscr%nJ
2529  nomega4m = ohscr%nomega
2530  ABI_STAT_MALLOC(epsm1,(npwe4mI,npwe4mJ,nomega4m,1), ierr)
2531  ABI_CHECK(ierr==0, 'out of memory in epsm1')
2532 
2533  do iqibz=1,ohscr%nqibz
2534 
2535    do ifile=1,nfiles
2536      ! allocate temporary array
2537      npwe4mI = Hscr_file(ifile)%npwe*Hscr_file(ifile)%nI
2538      npwe4mJ = Hscr_file(ifile)%npwe*Hscr_file(ifile)%nJ
2539      nomega4m = Hscr_file(ifile)%nomega
2540      ABI_STAT_MALLOC(epsm1_temp,(npwe4mI,npwe4mJ,nomega4m,1), ierr)
2541      ABI_CHECK(ierr==0, 'out of memory in epsm1_temp')
2542 
2543      ! read screening
2544      call read_screening(varname,filenames(ifile),npwe4mI,1,nomega4m,epsm1_temp,iomode,comm,iqiA=iqibz)
2545 
2546      ! Copy matrices for relevant frequencies
2547      do ifrq=1,nfreq_tot
2548        if (ifile_indx(ifrq)==ifile) then
2549          epsm1(:,:,ifrq,1)=epsm1_temp(:,:,freq_indx(ifrq,ifile),1)
2550        end if
2551      end do
2552 
2553      !Add imaginary part if needed
2554      !if (indx_imfreq_file==ifile) then
2555      !do ifrq=nfreqre+1,ohscr%nomega
2556      !epsm1(:,:,ifrq,1)=epsm1_temp(:,:,freq_indx(ifrq,ifile),1)
2557      !end do
2558      !end if
2559 
2560      ABI_FREE(epsm1_temp)
2561    end do !ifile
2562 
2563    ! Write data.
2564    npwe4mI = ohscr%npwe*ohscr%nI
2565    nomega4m = ohscr%nomega
2566    call write_screening(varname,ount,iomode,npwe4mI,nomega4m,iqibz,epsm1)
2567  end do ! iqibz
2568 
2569  ABI_FREE(epsm1)
2570  ABI_FREE(omega_storage)
2571  ABI_FREE(freq_indx)
2572  ABI_FREE(ifile_indx)
2573 
2574  if (iomode == IO_MODE_FORTRAN) then
2575    close(ount)
2576  else
2577 #ifdef HAVE_NETCDF
2578    NCF_CHECK(nf90_close(ount))
2579 #endif
2580  end if
2581 
2582  write(msg,'(3a)')ch10,' ==== Files have been merged successfully === ',ch10
2583  call wrtout(std_out,msg,'COLL')
2584 
2585 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.

PARENTS

      mrgscr

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

2616 subroutine ioscr_wremove(inpath, ihscr, fname_out, nfreq_tot, freq_indx, ohscr)
2617 
2618 
2619 !This section has been created automatically by the script Abilint (TD).
2620 !Do not modify the following lines by hand.
2621 #undef ABI_FUNC
2622 #define ABI_FUNC 'ioscr_wremove'
2623 !End of the abilint section
2624 
2625  implicit none
2626 
2627 !Arguments ------------------------------------
2628 !scalars
2629  integer,intent(in) :: nfreq_tot
2630  character(len=*),intent(in) :: inpath,fname_out
2631  type(hscr_t),intent(in) :: ihscr
2632  type(hscr_t),intent(out) :: ohscr
2633 !arrays
2634  integer,intent(in) :: freq_indx(nfreq_tot)
2635 
2636 !Local variables-------------------------------
2637 !scalars
2638  integer,parameter :: rdwr2=2,master=0
2639  integer :: iqibz,fform_merge,comm,nomega4m,ierr
2640  integer :: ifrq,npwe4mI,npwe4mJ,iomode,ount
2641  character(len=500) :: msg
2642  character(len=nctk_slen) :: varname
2643  type(abifile_t) :: abifile
2644 !arrays
2645  complex(gwpc),allocatable :: epsm1(:,:,:),epsm1_temp(:,:,:)
2646 
2647 ! *************************************************************************
2648 
2649  comm = xmpi_comm_self
2650 
2651  ! check ifreq_idx
2652  ABI_CHECK(nfreq_tot > 0, "nfreq_tot <= 0!")
2653  if (all(freq_indx == 0)) MSG_ERROR("all(freq_indx == 0)")
2654 
2655  ! Copy the old header
2656  call hscr_copy(ihscr, ohscr)
2657 
2658  ! Then modify entries for new frequency grid.
2659  ohscr%nomega = nfreq_tot
2660  ABI_FREE(ohscr%omega)
2661  ABI_MALLOC(ohscr%omega,(nfreq_tot))
2662  do ifrq=1,nfreq_tot
2663    ohscr%omega(ifrq) = ihscr%omega(freq_indx(ifrq))
2664  end do
2665 
2666  npwe4mI = ohscr%npwe*ohscr%nI
2667  npwe4mJ = ohscr%npwe*ohscr%nJ
2668 
2669  ! Print new header for info
2670  call hscr_print(ohscr,header='Header of the final file',unit=std_out,prtvol=1)
2671 
2672  ! Open output file.
2673  if (endswith(fname_out, ".nc")) then
2674    iomode = IO_MODE_ETSF
2675 #ifdef HAVE_NETCDF
2676    NCF_CHECK(nctk_open_create(ount, fname_out, comm))
2677 #else
2678    MSG_ERROR("Netcdf support is not available")
2679 #endif
2680  else
2681    iomode = IO_MODE_FORTRAN
2682    if (open_file(fname_out, msg, newunit=ount, status='new', form='unformatted') /= 0) then
2683      MSG_ERROR(msg)
2684    end if
2685  end if
2686 
2687  ! Write the header.
2688  fform_merge = ohscr%fform
2689  abifile = abifile_from_fform(fform_merge)
2690  if (abifile%fform == 0) then
2691     MSG_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
2692  end if
2693  varname = abifile%varname
2694 
2695  call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode)
2696 
2697  npwe4mI = ohscr%npwe*ohscr%nI; npwe4mJ = ohscr%npwe*ohscr%nJ
2698  nomega4m = ohscr%nomega
2699 
2700  ABI_STAT_MALLOC(epsm1, (npwe4mI,npwe4mJ,nomega4m), ierr)
2701  ABI_CHECK(ierr==0, 'out of memory in epsm1')
2702 
2703  do iqibz=1,ohscr%nqibz
2704    ! allocate temporary array
2705    npwe4mI = ihscr%npwe * ihscr%nI
2706    npwe4mJ = ihscr%npwe * ihscr%nJ
2707    nomega4m = ihscr%nomega
2708    ABI_STAT_MALLOC(epsm1_temp,(npwe4mI,npwe4mJ,nomega4m), ierr)
2709    ABI_CHECK(ierr==0, 'out of memory in epsm1_temp')
2710 
2711    ! read full screening matrix for this q-point
2712    call read_screening(varname,inpath,npwe4mI,1,nomega4m,epsm1_temp,iomode,comm,iqiA=iqibz)
2713 
2714    ! Copy relevant frequencies
2715    do ifrq=1,nfreq_tot
2716      epsm1(:,:,ifrq) = epsm1_temp(:,:,freq_indx(ifrq))
2717    end do
2718 
2719    ABI_FREE(epsm1_temp)
2720 
2721    npwe4mI = ohscr%npwe*ohscr%nI; nomega4m = ohscr%nomega
2722    call write_screening(varname,ount,iomode,npwe4mI,nomega4m,iqibz,epsm1)
2723  end do ! iqibz
2724 
2725  ABI_FREE(epsm1)
2726 
2727  if (iomode == IO_MODE_FORTRAN) then
2728    close(ount)
2729  else
2730 #ifdef HAVE_NETCDF
2731    NCF_CHECK(nf90_close(ount))
2732 #endif
2733  end if
2734 
2735  write(msg,'(3a)')ch10,' ==== Frequencies have been removed successfully === ',ch10
2736  call wrtout(std_out,msg,'COLL')
2737 
2738 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

272 character(len=nctk_slen) function ncname_from_id(id) result(varname)
273 
274 
275 !This section has been created automatically by the script Abilint (TD).
276 !Do not modify the following lines by hand.
277 #undef ABI_FUNC
278 #define ABI_FUNC 'ncname_from_id'
279 !End of the abilint section
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

PARENTS

      calc_ucrpa,m_io_screening,m_screen,m_screening,mrgscr

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

1597 subroutine read_screening(varname,fname,npweA,nqibzA,nomegaA,epsm1,iomode,comm,&
1598 & iqiA) ! Optional
1599 
1600 
1601 !This section has been created automatically by the script Abilint (TD).
1602 !Do not modify the following lines by hand.
1603 #undef ABI_FUNC
1604 #define ABI_FUNC 'read_screening'
1605 !End of the abilint section
1606 
1607  implicit none
1608 
1609 !Arguments ------------------------------------
1610 !scalars
1611  integer,intent(in) :: iomode,nomegaA,npweA,nqibzA,comm
1612  integer,optional,intent(in) :: iqiA
1613  character(len=*),intent(in) :: varname,fname
1614 !arrays
1615  complex(gwpc),target,intent(inout) :: epsm1(npweA,npweA,nomegaA,nqibzA)
1616 
1617 !Local variables-------------------------------
1618 !scalars
1619  integer,parameter :: master=0
1620  integer :: ipwe,fform,iomega,iqibz,unt,rdwr,my_rank,nprocs,my_iomode
1621 #ifdef HAVE_NETCDF
1622  integer :: varid,ncerr
1623 #endif
1624 #ifdef HAVE_MPI_IO
1625  integer :: test_fform,mpi_err,ierr,sc_mode
1626  integer :: bsize_frm,mpi_type_frm
1627  integer :: mpi_fh,buf_dim !,mat_ggw,mat_ggwq
1628  integer(XMPI_OFFSET_KIND) :: offset,displ_wq !,my_offpad
1629  !complex(dpc) :: ctmp
1630 #endif
1631  character(len=500) :: msg,errmsg
1632  logical :: read_qslice
1633  type(hscr_t) :: Hscr
1634 !arrays
1635 #ifdef HAVE_MPI_IO
1636  integer(MPI_OFFSET_KIND),allocatable :: offset_wq(:,:)
1637 #endif
1638  complex(dpc),allocatable :: bufdc2d(:,:),bufdc3d(:,:,:)
1639  ! pointers passed to netcdf4 routines (complex datatypes are not supported).
1640 #ifdef HAVE_GW_DPC
1641  real(dp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:)
1642 #else
1643  real(sp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:)
1644 #endif
1645  integer :: spins(2),s1,s2
1646 
1647 ! *************************************************************************
1648 
1649  DBG_ENTER("COLL")
1650 
1651  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1652  my_iomode = iomode
1653  if (endswith(fname, ".nc")) my_iomode = IO_MODE_ETSF
1654  !my_iomode = IO_MODE_MPI
1655  !if (my_iomode  == IO_MODE_MPI) my_iomode = IO_MODE_FORTRAN
1656 
1657  rdwr=1
1658  select case (my_iomode)
1659  case (IO_MODE_MPI)
1660 #ifdef HAVE_MPI_IO
1661    bsize_frm    = xmpio_bsize_frm    ! bsize_frm= Byte length of the Fortran record marker.
1662    mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker.
1663    sc_mode = xmpio_collective
1664 
1665    ! Master reads the header via Fortran IO then bcast the data.
1666    call hscr_from_file(hscr, fname, fform, comm)
1667 
1668    ! Open the file with MPI-IO
1669    call MPI_FILE_OPEN(comm, fname, MPI_MODE_RDONLY, xmpio_info ,mpi_fh, mpi_err)
1670    ABI_CHECK_MPI(mpi_err, sjoin("MPI_FILE_OPEN:", fname))
1671 
1672    ! Retrieve the offset of the section immediately below the header.
1673    call hscr_mpio_skip(mpi_fh,test_fform,offset)
1674    ABI_CHECK(test_fform == fform, "mismatch in fform!")
1675 
1676    ! Offsets of the Fortran markers corresponding to the (w,q) slices.
1677    ABI_MALLOC(offset_wq,(HScr%nomega,HScr%nqibz))
1678    displ_wq = offset
1679    do iqibz=1,Hscr%nqibz
1680      do iomega=1,Hscr%nomega
1681        ABI_CHECK(displ_wq > 0, "displ_wq < 0, your SCR|SUSC file is too big for MPI-IO!")
1682        offset_wq(iomega,iqibz) = displ_wq
1683        displ_wq = displ_wq + Hscr%npwe**2 * xmpi_bsize_dpc + Hscr%npwe * 2 * bsize_frm
1684      end do
1685    end do
1686 #else
1687    MSG_ERROR("MPI-IO support not enabled at configure-time")
1688 #endif
1689 
1690  case (IO_MODE_FORTRAN)
1691    ! Plain Fortran IO, all nodes read.
1692    if (open_file(fname,msg,newunit=unt,form="unformatted",status="old",action="read") /= 0) then
1693      MSG_ERROR(msg)
1694    end if
1695    call hscr_io(hscr,fform,rdwr,unt,comm,master,my_iomode)
1696 
1697  case (IO_MODE_ETSF)
1698 #ifdef HAVE_NETCDF
1699    NCF_CHECK(nctk_open_read(unt, fname, xmpi_comm_self))
1700    call hscr_io(hscr,fform,rdwr,unt,comm,master,my_iomode)
1701 #endif
1702 
1703  case default
1704    MSG_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode)))
1705  end select
1706 
1707  ! Slice or full array?
1708  read_qslice = .False.
1709  if (PRESENT(iqiA)) then
1710    read_qslice = .True.
1711    !call wrtout(std_out, sjoin('. Reading q-slice for iq = ',itoa(iqiA),' from: ', fname))
1712    if (iqiA <= 0 .or. iqiA > Hscr%nqibz) then
1713      MSG_BUG('iqiA out of range')
1714    end if
1715  end if
1716 
1717  ! Do some check
1718  if (Hscr%npwe>npweA) then
1719    write(msg,'(a,i0,2a,i0)')&
1720 &    'Total number of G-vectors reported on file = ',Hscr%npwe,ch10,&
1721 &    'Reading a smaller matrix of dimension      = ',npweA
1722    MSG_COMMENT(msg)
1723  end if
1724 
1725  if (npweA>Hscr%npwe) then
1726    write(msg,'(2(a,i0))')' Dimension of matrix = ',Hscr%npwe," requiring a too big matrix = ",npweA
1727    MSG_ERROR(msg)
1728  end if
1729 
1730  ABI_CHECK(nqibzA  <= Hscr%nqibz, 'Requiring too many q-points')
1731  ABI_CHECK(nomegaA <= Hscr%nomega,'Requiring too many frequencies')
1732 
1733  select case (my_iomode)
1734  case (IO_MODE_MPI)
1735 #ifdef HAVE_MPI_IO
1736    if (read_qslice) then
1737       call wrtout(std_out, "calling mpiotk to read_qslice", "COLL")
1738       buf_dim = (npweA)**2 * nomegaA
1739       offset = offset_wq(1,iqiA)
1740       sc_mode = xmpio_collective
1741 
1742 #ifdef HAVE_GW_DPC
1743      ! Read in-place.
1744      call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1745         buf_dim,epsm1,xmpio_chunk_bsize,sc_mode,comm,ierr)
1746      ABI_CHECK(ierr==0,"Fortran matrix too big")
1747 #else
1748      ! Have to allocate workspace for dpc data.
1749      ! FIXME: Change the file format of the SCR and SUC file so that
1750      ! they are written in single precision if not HAVE_GW_DPC
1751      ABI_STAT_MALLOC(bufdc3d,(npweA,npweA,nomegaA), ierr)
1752      ABI_CHECK(ierr==0, "out of memory bufdc3d")
1753 
1754      call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1755         buf_dim,bufdc3d,xmpio_chunk_bsize,sc_mode,comm,ierr)
1756      ABI_CHECK(ierr==0,"Fortran matrix too big")
1757 
1758      epsm1(:,:,:,1) = bufdc3d
1759      ABI_FREE(bufdc3d)
1760 #endif
1761 
1762    else
1763      ! Full matrix (G,G',w,q) is needed.
1764      call wrtout(std_out, "calling mpiotk: Full matrix (G,G',w,q) is needed.")
1765 
1766 #ifdef HAVE_GW_DPC
1767      ! Can read all data at once.
1768      buf_dim = (npweA)**2 * nomegaA * HScr%nqibz
1769      offset = offset_wq(1,1)
1770      sc_mode = xmpio_collective
1771 
1772      call mpiotk_read_fsuba_dpc4D(mpi_fh,offset,&
1773        [HScr%npwe,HScr%npwe,HScr%nomega,HScr%nqibz], [npweA,npweA,nomegaA,HScr%nqibz], [1,1,1,1],&
1774        buf_dim,epsm1,xmpio_chunk_bsize,sc_mode,comm,ierr)
1775      ABI_CHECK(ierr==0,"Fortran record too big")
1776 #else
1777      ! Have to allocate workspace for dpc data.
1778      ABI_STAT_MALLOC(bufdc3d,(npweA,npweA,nomegaA), ierr)
1779      ABI_CHECK(ierr==0, 'out of memory in bufdc3d, change the code to loop over frequencies!')
1780      sc_mode = xmpio_collective
1781 
1782      do iqibz=1,Hscr%nqibz
1783        offset = offset_wq(1,iqibz)
1784        buf_dim = (2*npweA)**2 * nomegaA
1785 
1786        call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, &
1787         [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1788          buf_dim,bufdc3d,xmpio_chunk_bsize,sc_mode,comm,ierr)
1789        ABI_CHECK(ierr==0,"Fortran matrix too big")
1790 
1791        epsm1(:,:,:,iqibz) = bufdc3d
1792      end do
1793 
1794      ABI_FREE(bufdc3d)
1795 #endif
1796    end if
1797 
1798    call MPI_FILE_CLOSE(mpi_fh,mpi_err)
1799    ABI_FREE(offset_wq)
1800 #endif
1801 
1802  case (IO_MODE_FORTRAN)
1803    ! Read epsilon^-1 with Fortran IO
1804    ! Allocate a single column to save memory.
1805    ! TODO re-merge the two cases.
1806    ABI_MALLOC(bufdc2d,(Hscr%npwe,1))
1807 
1808    ! Two coding for different case just to keep it readable.
1809    select case (read_qslice)
1810    case (.True.)
1811      ! Read only a slice of the full array (useful if the entire array is huge).
1812      !if (dim_wings==1) STOP 'not implemented'
1813      !TODO this has to be done in a cleaner way.
1814      qread_loop: &
1815 &    do iqibz=1,Hscr%nqibz
1816        if (iqibz==iqiA) then
1817          do iomega=1,nomegaA
1818            do ipwe=1,Hscr%npwe
1819              read(unt, err=10, iomsg=errmsg) bufdc2d(1:Hscr%npwe,1)
1820              if (ipwe<=npweA) epsm1(1:npweA,ipwe,iomega,1)=bufdc2d(1:npweA,1)
1821            end do
1822          end do
1823          EXIT qread_loop ! Got data. Do not need to read file till the end.
1824        else
1825          ! Skip other q-points i.e bufdc2d(1:Hscr%npwe,1:Hscr%npwe)
1826          do iomega=1,Hscr%nomega
1827            do ipwe=1,Hscr%npwe
1828             read(unt, err=10, iomsg=errmsg)
1829            end do
1830          end do
1831        end if ! iqibz==iqiA
1832      end do qread_loop ! iqibz
1833 
1834    case (.False.)
1835      ! Read the entire array.
1836      do iqibz=1,Hscr%nqibz
1837        do iomega=1,nomegaA
1838          do ipwe=1,Hscr%npwe
1839           read(unt, err=10, iomsg=errmsg) bufdc2d(1:Hscr%npwe,1)
1840           if (ipwe<=npweA) epsm1(1:npweA,ipwe,iomega,iqibz)=bufdc2d(1:npweA,1)
1841          end do
1842        end do
1843        ! Skip other frequencies
1844        do iomega=nomegaA+1,Hscr%nomega
1845          do ipwe=1,Hscr%npwe
1846            read(unt, err=10, iomsg=errmsg)
1847          end do
1848        end do
1849      end do !iqibz
1850    end select
1851 
1852    close(unt)
1853 
1854  case (IO_MODE_ETSF)
1855 #ifdef HAVE_NETCDF
1856    ! netcdf does not support complex datatypes. Here I use some C-magic to  associate the memory
1857    ! to a Fortran real pointer with the correct type and shape. Note that the data on file is always in double precision.
1858    ! nf90_get_var will automatically convert from double to single if the GW code is in single precision mode.
1859    ! This is the reason why I'm using CPP option in the declaration of real_epsm1.
1860 
1861    ! FIXME: Need to know the type to read
1862    !write(std_out,*)"in read_screening"
1863    varid = nctk_idname(unt, varname)
1864 
1865    ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt]
1866    call c_f_pointer(c_loc(epsm1(1,1,1,1)), real_epsm1, [2,npweA,npweA,1,1,nomegaA,nqibzA])
1867    spins = 1; s1 = spins(1); s2 = spins(2)
1868    if (read_qslice) then
1869      ncerr = nf90_get_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,iqia], count=[2,npweA,npweA,1,1,nomegaA,1])
1870    else
1871      ncerr = nf90_get_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,1], count=[2,npweA,npweA,1,1,nomegaA,nqibzA])
1872      !do iqibz=1,nqibzA
1873      !  !write(*,*)"epsm1: in read ",iqibz,epsm1(1:3,1,1,iqibz)
1874      !end do
1875    end if
1876    NCF_CHECK_MSG(ncerr, sjoin("getting var:", varname))
1877    NCF_CHECK(nf90_close(unt))
1878    !write(std_out,*)"read_screening done"
1879 #endif
1880 
1881  case default
1882    MSG_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode)))
1883  end select
1884 
1885  ! Free memory
1886  if (allocated(bufdc2d)) then
1887    ABI_FREE(bufdc2d)
1888  end if
1889  if (allocated(bufdc3d)) then
1890    ABI_FREE(bufdc3d)
1891  end if
1892 
1893  call hscr_free(Hscr)
1894 
1895  DBG_EXIT("COLL")
1896 
1897  return
1898 
1899  ! Handle Fortran IO error.
1900 10 continue
1901  MSG_ERROR(errmsg)
1902 
1903 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)

PARENTS

      m_io_screening,m_screening,screening

CHILDREN

      hscr_copy,hscr_io,hscr_print,read_screening,write_screening,wrtout

SOURCE

1480 subroutine write_screening(varname,unt,iomode,npwe,nomega,iqibz,epsm1)
1481 
1482 
1483 !This section has been created automatically by the script Abilint (TD).
1484 !Do not modify the following lines by hand.
1485 #undef ABI_FUNC
1486 #define ABI_FUNC 'write_screening'
1487 !End of the abilint section
1488 
1489  implicit none
1490 
1491 !Arguments ------------------------------------
1492 !scalars
1493  character(len=*),intent(in) :: varname
1494  integer,intent(in) :: nomega,npwe,iqibz,unt,iomode
1495 !arrays
1496  complex(gwpc),target,intent(in) :: epsm1(npwe,npwe,nomega)
1497 
1498 !Local variables-------------------------------
1499 !scalars
1500  integer :: ipwe,iomega,spins(2),s1,s2
1501  character(len=500) :: errmsg
1502 !arrays
1503  complex(dpc),allocatable :: epsm1d(:,:)
1504 #ifdef HAVE_NETCDF
1505  integer :: varid,ncerr
1506 #endif
1507 #ifdef HAVE_GW_DPC
1508  real(dp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:)
1509 #else
1510  real(sp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:)
1511 #endif
1512 
1513 ! *************************************************************************
1514 
1515  DBG_ENTER("COLL")
1516 
1517  select case (iomode)
1518  case (IO_MODE_FORTRAN, IO_MODE_MPI)
1519    ! Write a record for each omega, Always use double precision.
1520    ABI_MALLOC(epsm1d,(npwe,1))
1521 
1522    do iomega=1,nomega
1523      do ipwe=1,npwe
1524        epsm1d(:,1) = epsm1(:,ipwe,iomega) !spc ==> dpc
1525        write(unt, err=10, iomsg=errmsg)epsm1d(1:npwe,1)
1526      end do
1527    end do
1528    ABI_FREE(epsm1d)
1529 
1530 #ifdef HAVE_NETCDF
1531  case (IO_MODE_ETSF)
1532    ! netcdf does not support complex datatypes. Here I use some C-magic to  associate the memory
1533    ! to a Fortran real pointer with the correct type and shape. Note that the data on file is always in double precision.
1534    varid = nctk_idname(unt, varname)
1535    call c_f_pointer(c_loc(epsm1(1,1,1)), real_epsm1, [2,npwe,npwe,1,1,nomega,1])
1536    ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt]
1537    spins = 1; s1 = spins(1); s2 = spins(2)
1538    ncerr = nf90_put_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,iqibz], count=[2,npwe,npwe,1,1,nomega,1])
1539    NCF_CHECK_MSG(ncerr, sjoin("putting var:", varname))
1540 #endif
1541 
1542  case default
1543    MSG_ERROR(sjoin("Wrong iomode:", iomode2str(iomode)))
1544  end select
1545 
1546  DBG_EXIT("COLL")
1547 
1548  return
1549 
1550  ! Handle IO error
1551  10 continue
1552  MSG_ERROR(errmsg)
1553 
1554 end subroutine write_screening