TABLE OF CONTENTS
- ABINIT/m_io_screening
- m_io_screening/hscr_bcast
- m_io_screening/hscr_copy
- m_io_screening/hscr_free
- m_io_screening/hscr_from_file
- m_io_screening/hscr_io
- m_io_screening/hscr_malloc
- m_io_screening/hscr_merge
- m_io_screening/hscr_mpio_skip
- m_io_screening/hscr_new
- m_io_screening/hscr_print
- m_io_screening/hscr_t
- m_io_screening/ioscr_qmerge
- m_io_screening/ioscr_qrecover
- m_io_screening/ioscr_wmerge
- m_io_screening/ioscr_wremove
- m_io_screening/ncname_from_id
- m_io_screening/read_screening
- m_io_screening/write_screening
ABINIT/m_io_screening [ Modules ]
NAME
m_io_screening
FUNCTION
This module contains the definition of the header of the _SCR and _SUSC file as well as methods used to read/write/echo.
COPYRIGHT
Copyright (C) 2008-2022 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_io_screening 24 25 use defs_basis 26 use m_abicore 27 #if defined HAVE_MPI2 28 use mpi 29 #endif 30 use m_xmpi 31 use m_mpiotk 32 use m_nctk 33 use m_errors 34 use m_dtset 35 use, intrinsic :: iso_c_binding 36 use netcdf 37 use m_hdr 38 use m_sort 39 40 use m_gwdefs, only : em1params_t, GW_TOLQ 41 use m_fstrings, only : sjoin, itoa, endswith, replace_ch0 42 use m_copy, only : alloc_copy 43 use m_io_tools, only : open_file, file_exists, iomode2str 44 use m_numeric_tools, only : print_arr, remove_copies, imax_loc 45 use m_bz_mesh, only : isequalk 46 47 implicit none 48 49 private
m_io_screening/hscr_bcast [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_bcast
FUNCTION
This subroutine transmit the header structured datatype initialized on one processor (or a group of processor), to the other processors. It also allocates the needed part of the header.
INPUTS
master=ID of the master node. my_rank=ID of the node that receives the data. comm=MPI communicator.
OUTPUT
(no output)
SIDE EFFECTS
Hscr<type(hscr_t)>=the SCR header. For the master, it is already initialized entirely, while for the other procs, everything has to be transmitted.
NOTES
This routine is called only in the case of MPI version of the code.
SOURCE
931 subroutine hscr_bcast(hscr, master, my_rank, comm) 932 933 !Arguments ------------------------------------ 934 class(hscr_t),intent(inout) :: hscr 935 integer, intent(in) :: master, my_rank, comm 936 937 !Local variables------------------------------- 938 integer :: ierr 939 940 ! ************************************************************************* 941 942 DBG_ENTER("COLL") 943 if (xmpi_comm_size(comm) == 1) return ! Nothing to do 944 945 ! integer 946 call xmpi_bcast(hscr%id, master,comm,ierr) 947 call xmpi_bcast(hscr%ikxc, master,comm,ierr) 948 call xmpi_bcast(hscr%inclvkb, master,comm,ierr) 949 call xmpi_bcast(hscr%headform, master,comm,ierr) 950 call xmpi_bcast(hscr%fform, master,comm,ierr) 951 call xmpi_bcast(hscr%gwcalctyp, master,comm,ierr) 952 call xmpi_bcast(hscr%nI, master,comm,ierr) 953 call xmpi_bcast(hscr%nJ, master,comm,ierr) 954 call xmpi_bcast(hscr%nqibz, master,comm,ierr) 955 call xmpi_bcast(hscr%nqlwl, master,comm,ierr) 956 call xmpi_bcast(hscr%nomega, master,comm,ierr) 957 call xmpi_bcast(hscr%nbnds_used, master,comm,ierr) 958 call xmpi_bcast(hscr%npwe, master,comm,ierr) 959 call xmpi_bcast(hscr%npwwfn_used,master,comm,ierr) 960 call xmpi_bcast(hscr%spmeth, master,comm,ierr) 961 call xmpi_bcast(hscr%test_type, master,comm,ierr) 962 call xmpi_bcast(hscr%tordering, master,comm,ierr) 963 964 ! Real 965 call xmpi_bcast(hscr%mbpt_sciss, master,comm,ierr) 966 call xmpi_bcast(hscr%spsmear, master,comm,ierr) 967 call xmpi_bcast(hscr%zcut, master,comm,ierr) 968 969 ! arrays 970 call xmpi_bcast(hscr%titles, master,comm,ierr) 971 972 if (my_rank /= master) call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl) 973 974 call xmpi_bcast(hscr%gvec, master,comm,ierr) 975 call xmpi_bcast(hscr%qibz, master,comm,ierr) 976 call xmpi_bcast(hscr%qlwl, master,comm,ierr) 977 call xmpi_bcast(hscr%omega,master,comm,ierr) 978 979 ! Communicate the Abinit header. 980 call hscr%Hdr%bcast(master, my_rank, comm) 981 982 ! HSCR_NEW 983 call xmpi_bcast(hscr%awtr, master, comm, ierr) 984 call xmpi_bcast(hscr%icutcoul, master, comm, ierr) 985 call xmpi_bcast(hscr%vcutgeo, master, comm, ierr) 986 call xmpi_bcast(hscr%gwcomp, master, comm, ierr) 987 call xmpi_bcast(hscr%gwgamma, master, comm, ierr) 988 call xmpi_bcast(hscr%gwencomp, master, comm, ierr) 989 call xmpi_bcast(hscr%kind_cdata, master, comm, ierr) 990 ! HSCR_NEW 991 992 DBG_EXIT("COLL") 993 994 end subroutine hscr_bcast
m_io_screening/hscr_copy [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_copy
FUNCTION
Deep copy of the header of the _SCR or _SUSC file.
INPUTS
SOURCE
1073 subroutine hscr_copy(Hscr_in, Hscr_cp) 1074 1075 !Arguments ------------------------------------ 1076 !scalars 1077 class(hscr_t),intent(in) :: Hscr_in 1078 class(hscr_t),intent(inout) :: Hscr_cp 1079 1080 ! ************************************************************************* 1081 1082 !@hscr_t 1083 ! Integer values. 1084 Hscr_cp%id = Hscr_in%id 1085 Hscr_cp%ikxc = Hscr_in%ikxc 1086 Hscr_cp%inclvkb = Hscr_in%inclvkb 1087 Hscr_cp%headform = Hscr_in%headform 1088 Hscr_cp%fform = Hscr_in%fform 1089 Hscr_cp%gwcalctyp = Hscr_in%gwcalctyp 1090 Hscr_cp%nI = Hscr_in%nI 1091 Hscr_cp%nJ = Hscr_in%nJ 1092 Hscr_cp%nqibz = Hscr_in%nqibz 1093 Hscr_cp%nqlwl = Hscr_in%nqlwl 1094 Hscr_cp%nomega = Hscr_in%nomega 1095 Hscr_cp%nbnds_used = Hscr_in%nbnds_used 1096 Hscr_cp%npwe = Hscr_in%npwe 1097 Hscr_cp%npwwfn_used = Hscr_in%npwwfn_used 1098 Hscr_cp%spmeth = Hscr_in%spmeth 1099 Hscr_cp%test_type = Hscr_in%test_type 1100 Hscr_cp%tordering = Hscr_in%tordering 1101 1102 ! Real variables 1103 Hscr_cp%mbpt_sciss = Hscr_in%mbpt_sciss 1104 Hscr_cp%spsmear = Hscr_in%spsmear 1105 Hscr_cp%zcut = Hscr_in%zcut 1106 1107 ! Copy the abinit Header 1108 call hdr_copy(Hscr_in%Hdr,Hscr_cp%Hdr) 1109 1110 Hscr_cp%titles(:) = Hscr_in%titles(:) 1111 1112 ! Copy allocatable arrays. 1113 call alloc_copy(Hscr_in%gvec , Hscr_cp%gvec) 1114 call alloc_copy(Hscr_in%qibz , Hscr_cp%qibz) 1115 call alloc_copy(Hscr_in%qlwl , Hscr_cp%qlwl) 1116 call alloc_copy(Hscr_in%omega, Hscr_cp%omega) 1117 1118 ! HSCR_NEW 1119 hscr_cp%awtr = hscr_in%awtr 1120 hscr_cp%icutcoul = hscr_in%icutcoul 1121 hscr_cp%vcutgeo = hscr_in%vcutgeo 1122 hscr_cp%gwcomp = hscr_in%gwcomp 1123 hscr_cp%gwgamma = hscr_in%gwgamma 1124 hscr_cp%gwencomp = hscr_in%gwencomp 1125 hscr_cp%kind_cdata = hscr_in%kind_cdata 1126 ! HSCR_NEW 1127 1128 end subroutine hscr_copy
m_io_screening/hscr_free [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_free
FUNCTION
Deallocate the components of the header structured datatype
INPUTS
hdr <type(hdr_type)>=the header
OUTPUT
(only deallocate)
SOURCE
1043 subroutine hscr_free(hscr) 1044 1045 !Arguments ------------------------------------ 1046 class(hscr_t),intent(inout) :: hscr 1047 1048 ! ************************************************************************* 1049 1050 ABI_SFREE(hscr%gvec) 1051 ABI_SFREE(hscr%qibz) 1052 ABI_SFREE(hscr%qlwl) 1053 ABI_SFREE(hscr%omega) 1054 1055 call hscr%Hdr%free() 1056 1057 end subroutine hscr_free
m_io_screening/hscr_from_file [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_from_file
FUNCTION
Read the header of the (SCR/SUS) file
INPUTS
path=File name comm = MPI communicator.
OUTPUT
hscr<hscr_t>=The header. fform=Kind of the array in the file (0 signals an error)
SOURCE
309 subroutine hscr_from_file(hscr, path, fform, comm) 310 311 !Arguments ------------------------------------ 312 !scalars 313 class(hscr_t),intent(out) :: hscr 314 character(len=*),intent(in) :: path 315 integer,intent(in) :: comm 316 integer,intent(out) :: fform 317 318 !Local variables------------------------------- 319 !scalars 320 integer,parameter :: rdwr5 = 5, master = 0 321 integer :: unt,my_rank,ierr 322 character(len=500) :: msg 323 324 ! ************************************************************************* 325 326 my_rank = xmpi_comm_rank(comm) 327 328 ! Master reads and broadcasts. 329 if (my_rank == master) then 330 if (.not. endswith(path, ".nc")) then 331 ! Fortran-IO 332 if (open_file(path,msg,newunit=unt,form="unformatted", status="old",action="read") /= 0) then 333 ABI_ERROR(msg) 334 end if 335 call hscr_io(hscr,fform,rdwr5,unt,xmpi_comm_self,master,IO_MODE_FORTRAN) 336 close(unt) 337 else 338 ! Netcdf format 339 NCF_CHECK(nctk_open_read(unt, path, xmpi_comm_self)) 340 call hscr_io(hscr,fform,rdwr5,unt,xmpi_comm_self,master,IO_MODE_ETSF) 341 NCF_CHECK(nf90_close(unt)) 342 end if 343 344 ABI_CHECK(fform /= 0, sjoin("hscr_io returned fform == 0 while reading:", path)) 345 end if 346 347 ! Broadcast data. 348 if (xmpi_comm_size(comm) > 1) then 349 call hscr%bcast(master, my_rank, comm) 350 call xmpi_bcast(fform,master,comm,ierr) 351 end if 352 353 end subroutine hscr_from_file
m_io_screening/hscr_io [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_io
FUNCTION
This subroutine deals with the I/O of the hscr_t structured variables (read/write/echo). According to the value of rdwr, it reads the header of a file, writes it, or echo the value of the structured variable to a file. Note that, when reading, different records of hscr_t are allocated here, according to the values of the read variables. Records of hscr_t should be deallocated correctly by a call to hdr_free when hscr_t is not used anymore.
INPUTS
iomode=Option defining the file format of the external file. comm=MPI communicator. master=rank of the master node in comm, usually 0 rdwr= if 1, read the hscr_t structured variable from the header of the file, if 2, write the header to unformatted file if 3, echo part of the header to formatted file (records 1 and 2) if 4, echo the header to formatted file if 5, read the hscr_t without rewinding (unformatted) if 6, write the hscr_t without rewinding (unformatted) unt=unit number of the file (unformatted if rdwr=1, 2, 5 or 6 formatted if rdwr=3,4)
OUTPUT
(see side effects)
SIDE EFFECTS
The following variables are both input or output : fform=kind of the array in the file if rdwr=1,5 : will be output ; if the reading fail, return fform=0 if rdwr=2,3,4,6 : should be input, will be written or echo to file hscr_t <type(hscr_t)>=the header structured variable if rdwr=1,5 : will be output if rdwr=2,3,4,6 : should be input, will be written or echo to file
NOTES
In all cases, the file is supposed to be open already When reading (rdwr=1) or writing (rdwr=2), rewind the file When echoing (rdwr=3) does not rewind the file. When reading (rdwr=5) or writing (rdwr=6), DOES NOT rewind the file In writing mode, the routine is supposed to called by the master node. no check is done, it is up to the developer.
SOURCE
402 subroutine hscr_io(hscr, fform, rdwr, unt, comm, master, iomode) 403 404 !Arguments ------------------------------------ 405 !scalars 406 class(hscr_t),target,intent(inout) :: hscr 407 integer,intent(inout) :: fform 408 integer,intent(in) :: rdwr,unt,iomode,comm,master 409 410 !Local variables------------------------------- 411 !scalars 412 integer :: my_rank,nprocs,ncerr,ncid,varid,ierr !ii 413 character(len=500) :: errmsg 414 character(len=nctk_slen) :: varname !,head_shape,wing_shape 415 !arrays 416 real(dp),allocatable :: real_omega(:,:) 417 real(dp), ABI_CONTIGUOUS pointer :: r2vals(:,:) !,rvals3(:,:,:) 418 ! ************************************************************************* 419 420 DBG_ENTER("COLL") 421 !@hscr_t 422 ABI_UNUSED(master) ! FIXME 423 424 ! Initialize MPI info for comm === 425 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 426 427 if (rdwr==1 .or. rdwr==5) then 428 429 if (.True.) then 430 ! TODO: only master should read but then I have to skip the header. 431 !if (my_rank == master) then 432 ! Read the abinit header, rewinding of the file (if any) is done here. 433 if (iomode==IO_MODE_FORTRAN) then 434 call hdr_fort_read(hscr%hdr, unt, fform, rewind=(rdwr==1)) 435 else if (iomode==IO_MODE_ETSF) then 436 call hdr_ncread(hscr%hdr, unt, fform) 437 end if 438 439 ! Reset the variables absent in old versions. 440 Hscr%fform=fform 441 442 if (iomode==IO_MODE_FORTRAN .or. iomode==IO_MODE_MPI) then 443 select case (fform) 444 case (1003, 1004) 445 ! File format for epsilon^-1, espilon, chi0 446 read(unt, err=10, iomsg=errmsg)hscr%titles 447 read(unt, err=10, iomsg=errmsg)& 448 hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp,& 449 hscr%nI, hscr%nJ, hscr%nqibz, hscr%nqlwl, hscr%nomega, hscr%nbnds_used,& 450 hscr%npwe, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering,& 451 hscr%awtr, hscr%icutcoul, hscr%gwgamma, hscr%vcutgeo(1:3) 452 453 ! Read real scalars 454 read(unt, err=10, iomsg=errmsg)& 455 hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp,hscr%kind_cdata 456 457 ! Allocate arrays and read them 458 call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl) 459 460 read(unt, err=10, iomsg=errmsg)hscr%gvec(:,:) 461 read(unt, err=10, iomsg=errmsg)hscr%qibz(:,:) 462 read(unt, err=10, iomsg=errmsg)hscr%omega(:) 463 464 ! Read data for q-->0 limit. 465 if (hscr%nqlwl>0) then 466 read(unt, err=10, iomsg=errmsg)hscr%qlwl(:,:) 467 end if 468 469 case default 470 ABI_BUG(sjoin('Wrong fform read:', itoa(fform))) 471 end select 472 473 else if (iomode == IO_MODE_ETSF) then 474 ncid = unt 475 476 select case (fform) 477 case (1003, 1004) 478 ! Get dimensions and allocate arrays. 479 NCF_CHECK(nctk_get_dim(ncid, "number_of_coefficients_dielectric_function", hscr%npwe)) 480 NCF_CHECK(nctk_get_dim(ncid, "number_of_qpoints_dielectric_function", hscr%nqibz)) 481 NCF_CHECK(nctk_get_dim(ncid, "number_of_frequencies_dielectric_function", hscr%nomega)) 482 NCF_CHECK(nctk_get_dim(ncid, "number_of_qpoints_gamma_limit", hscr%nqlwl)) 483 NCF_CHECK(nctk_get_dim(ncid, "nI", hscr%ni)) 484 NCF_CHECK(nctk_get_dim(ncid, "nJ", hscr%nj)) 485 call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl) 486 487 varid = nctk_idname(ncid, 'reduced_coordinates_plane_waves_dielectric_function') 488 NCF_CHECK(nf90_get_var(ncid, varid, hscr%gvec, start=[1,1,1])) 489 NCF_CHECK(nf90_get_var(ncid, vid('qpoints_dielectric_function'), hscr%qibz)) 490 491 ABI_MALLOC(real_omega, (2, hscr%nomega)) 492 NCF_CHECK(nf90_get_var(ncid, vid('frequencies_dielectric_function'), real_omega)) 493 hscr%omega = dcmplx(real_omega(1,:), real_omega(2,:)) 494 ABI_FREE(real_omega) 495 496 ! Read extra data added in new header. 497 NCF_CHECK(nf90_get_var(ncid, vid("vcutgeo"), hscr%vcutgeo)) 498 NCF_CHECK(nf90_get_var(ncid, vid("id"), hscr%id)) 499 NCF_CHECK(nf90_get_var(ncid, vid("ikxc"), hscr%ikxc)) 500 NCF_CHECK(nf90_get_var(ncid, vid("inclvkb"), hscr%inclvkb)) 501 NCF_CHECK(nf90_get_var(ncid, vid("headform"), hscr%headform)) 502 NCF_CHECK(nf90_get_var(ncid, vid("fform"), hscr%fform)) 503 NCF_CHECK(nf90_get_var(ncid, vid("gwcalctyp"), hscr%gwcalctyp)) 504 NCF_CHECK(nf90_get_var(ncid, vid("nbands_used"), hscr%nbnds_used)) 505 NCF_CHECK(nf90_get_var(ncid, vid("npwwfn_used"), hscr%npwwfn_used)) 506 NCF_CHECK(nf90_get_var(ncid, vid("spmeth"), hscr%spmeth)) 507 NCF_CHECK(nf90_get_var(ncid, vid("test_type"), hscr%test_type)) 508 NCF_CHECK(nf90_get_var(ncid, vid("tordering"), hscr%tordering)) 509 NCF_CHECK(nf90_get_var(ncid, vid("awtr"), hscr%awtr)) 510 NCF_CHECK(nf90_get_var(ncid, vid("gw_icutcoul"), hscr%icutcoul)) 511 NCF_CHECK(nf90_get_var(ncid, vid("gwcomp"), hscr%gwcomp)) 512 NCF_CHECK(nf90_get_var(ncid, vid("gwgamma"), hscr%gwgamma)) 513 NCF_CHECK(nf90_get_var(ncid, vid("mbpt_sciss"), hscr%mbpt_sciss)) 514 NCF_CHECK(nf90_get_var(ncid, vid("spsmear"), hscr%spsmear)) 515 NCF_CHECK(nf90_get_var(ncid, vid("zcut"), hscr%zcut)) 516 NCF_CHECK(nf90_get_var(ncid, vid("gwencomp"), hscr%gwencomp)) 517 NCF_CHECK(nf90_get_var(ncid, vid("kind_cdata"), hscr%kind_cdata)) 518 call replace_ch0(hscr%kind_cdata) 519 520 NCF_CHECK(nf90_get_var(ncid, vid("titles"), hscr%titles)) 521 call replace_ch0(hscr%titles(:)) 522 523 ! TODO Read it 524 if (hscr%nqlwl /= 0) then 525 NCF_CHECK(nf90_get_var(ncid, vid("qpoints_gamma_limit"), hscr%qlwl)) 526 end if 527 528 case default 529 ABI_BUG(sjoin('Unsupported fform read:',itoa(fform))) 530 end select 531 else 532 ABI_ERROR(sjoin("Unsupported value of iomode:", iomode2str(iomode))) 533 end if 534 535 end if ! master 536 537 !call hscr%bcast(master, my_rank, comm) 538 !call hscr_mpio_skip(mpio_fh,fform,offset) 539 540 else if (rdwr == 2 .or. rdwr == 6) then 541 ! Writing the header of an unformatted file. 542 ! Always use the latest version. 543 544 if (iomode==IO_MODE_FORTRAN .or. iomode==IO_MODE_MPI) then 545 ! Write the abinit header. 546 call hscr%hdr%fort_write(unt, fform, ierr) 547 ABI_CHECK(ierr == 0, "hdr_fort_write retured ierr != 0") 548 549 write(unt, err=10, iomsg=errmsg)hscr%titles 550 551 ! Write integers 552 write(unt, err=10, iomsg=errmsg)& 553 hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp,& 554 hscr%nI, hscr%nJ, hscr%nqibz, hscr%nqlwl, hscr%nomega, hscr%nbnds_used,& 555 hscr%npwe, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering,& 556 hscr%awtr, hscr%icutcoul, hscr%gwgamma, hscr%vcutgeo(1:3) 557 558 ! Write real scalars 559 write(unt, err=10, iomsg=errmsg)& 560 hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp,hscr%kind_cdata 561 562 ! Write arrays 563 write(unt, err=10, iomsg=errmsg)hscr%gvec(:,:) 564 write(unt, err=10, iomsg=errmsg)hscr%qibz(:,:) 565 write(unt, err=10, iomsg=errmsg)hscr%omega(:) 566 567 ! Add q-points for heads and wings for q-->0. 568 if (hscr%nqlwl > 0) then 569 write(unt, err=10, iomsg=errmsg)hscr%qlwl(:,:) 570 end if 571 572 else if (iomode == IO_MODE_ETSF) then 573 ncid = unt 574 ! Write the abinit header, rewinding of the file (if any) is done here. 575 NCF_CHECK(hscr%hdr%ncwrite(ncid, fform, nc_define=.True.)) 576 577 ! Define dimensions 578 ! Part 2) of etsf-io specifications 579 ! FIXME: Spin is only used in particular cases, We usuall get the trace of W in spin space 580 ! and I'm not gonna allocate extra memory just to have up up, down down 581 ! Besides number_of_spins should be replaced by `number_of_spins_dielectric_function` 582 ! Should add spin_dependent attribute. 583 ncerr = nctk_def_dims(ncid, [ & 584 nctkdim_t("complex", 2), nctkdim_t("number_of_reduced_dimensions", 3), & 585 nctkdim_t("number_of_frequencies_dielectric_function", hscr%nomega), & 586 nctkdim_t("number_of_qpoints_dielectric_function", hscr%nqibz), & 587 nctkdim_t("number_of_qpoints_gamma_limit", hscr%nqlwl), & 588 nctkdim_t("number_of_spins", hscr%hdr%nsppol), & 589 nctkdim_t("nI", hscr%nI), nctkdim_t("nJ", hscr%nJ), & 590 nctkdim_t("number_of_coefficients_dielectric_function", hscr%npwe)], defmode=.True.) 591 NCF_CHECK(ncerr) 592 593 ! Part 3) of the specs 594 ! (note that, in the specs, the Gs depend on the q-point but npwe is a scalar 595 ! basis_set is added by the abinit header. 596 ! FIXME: g-vectors are not written properly. 597 ncerr = nctk_def_arrays(ncid, [& 598 ! Standard 599 nctkarr_t('frequencies_dielectric_function', "dp", 'complex, number_of_frequencies_dielectric_function'), & 600 nctkarr_t('qpoints_dielectric_function', "dp", 'number_of_reduced_dimensions, number_of_qpoints_dielectric_function'),& 601 nctkarr_t('qpoints_gamma_limit', "dp", 'number_of_reduced_dimensions, number_of_qpoints_gamma_limit'), & 602 nctkarr_t('reduced_coordinates_plane_waves_dielectric_function', "i", & 603 'number_of_reduced_dimensions, number_of_coefficients_dielectric_function, number_of_qpoints_dielectric_function'), & 604 ! Abinit 605 nctkarr_t('vcutgeo', "dp", 'number_of_reduced_dimensions'), & 606 nctkarr_t('kind_cdata', "char", 'character_string_length'), & 607 nctkarr_t('titles', "char", 'character_string_length, two') & 608 ]) 609 NCF_CHECK(ncerr) 610 611 ! FIXME problem with q-points, heads and wings? 612 ! The order in P in the specs is wrong, q should be the last dimension here I use the "correct" version 613 ! TODO: Remove. Use abifile_t 614 varname = ncname_from_id(hscr%id) 615 ncerr = nctk_def_arrays(ncid, & 616 nctkarr_t(varname, "dp", & 617 &"complex, number_of_coefficients_dielectric_function, number_of_coefficients_dielectric_function,& 618 &number_of_spins, number_of_spins, number_of_frequencies_dielectric_function, number_of_qpoints_dielectric_function")) 619 NCF_CHECK(ncerr) 620 621 !write(std_out,*)"nqlwl",hscr%nqlwl 622 NCF_CHECK(nctk_set_datamode(ncid)) 623 call c_f_pointer(c_loc(hscr%omega(1)), r2vals, shape=[2, size(hscr%omega)]) 624 NCF_CHECK(nf90_put_var(ncid, vid('frequencies_dielectric_function'), r2vals)) 625 NCF_CHECK(nf90_put_var(ncid, vid('qpoints_dielectric_function'), hscr%qibz)) 626 NCF_CHECK(nf90_put_var(ncid, vid('reduced_coordinates_plane_waves_dielectric_function'), hscr%gvec)) 627 628 NCF_CHECK(nf90_put_var(ncid, vid("titles"), hscr%titles)) 629 NCF_CHECK(nf90_put_var(ncid, vid("kind_cdata"), hscr%kind_cdata)) 630 NCF_CHECK(nf90_put_var(ncid, vid("vcutgeo"), hscr%vcutgeo)) 631 632 ncerr = nctk_defnwrite_ivars(ncid, [character(len=nctk_slen) :: & 633 "id", "ikxc", "inclvkb", "headform", "fform", "gwcalctyp", & 634 "nbands_used", "npwwfn_used", "spmeth", "test_type", "tordering", "awtr", "gw_icutcoul", & 635 "gwcomp", "gwgamma" & 636 ],& 637 [ hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp, & 638 hscr%nbnds_used, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering, hscr%awtr, hscr%icutcoul, & 639 hscr%gwcomp, hscr%gwgamma & 640 ]) 641 NCF_CHECK(ncerr) 642 643 ncerr = nctk_defnwrite_dpvars(ncid, [character(len=nctk_slen) :: & 644 "mbpt_sciss", "spsmear", "zcut", "gwencomp"], & 645 [hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp & 646 ]) 647 NCF_CHECK(ncerr) 648 649 ! Add q-points for heads and wings for q-->0. 650 if (hscr%nqlwl > 0) then 651 ! MG: This part has been commented out as it's not used 652 ! head_shape = "complex, number_of_spins, number_of_spins, number_of_frequencies_dielectric_function" 653 ! head_shape = trim(head_shape)//", number_of_qpoints_gamma_limit" 654 655 ! wing_shape = "complex, number_of_coefficients_dielectric_function, number_of_spins, number_of_spins" 656 ! wing_shape = trim(wing_shape)//", number_of_frequencies_dielectric_function, number_of_qpoints_gamma_limit" 657 658 ! ncerr = nctk_def_arrays(ncid, [& 659 ! nctkarr_t("dielectric_function_head", "dp", head_shape),& 660 ! nctkarr_t("dielectric_function_upper_wing", "dp", wing_shape),& 661 ! nctkarr_t("dielectric_function_lower_wing", "dp", wing_shape)], defmode=.True.) 662 ! NCF_CHECK(ncerr) 663 664 NCF_CHECK(nctk_set_datamode(ncid)) 665 NCF_CHECK(nf90_put_var(ncid, vid('qpoints_gamma_limit'), hscr%qlwl)) 666 end if 667 else 668 ABI_ERROR(sjoin('Unsupported iomode:',iomode2str(iomode))) 669 end if 670 671 else 672 ABI_BUG(sjoin("Wrong value for rdwr:", itoa(rdwr))) 673 end if ! read/write/echo 674 675 DBG_EXIT("COLL") 676 677 return 678 679 ! Handle Fortran IO error 680 10 continue 681 ABI_ERROR(errmsg) 682 683 contains 684 integer function vid(vname) 685 character(len=*),intent(in) :: vname 686 vid = nctk_idname(ncid, vname) 687 end function vid 688 689 end subroutine hscr_io
m_io_screening/hscr_malloc [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_malloc
FUNCTION
Allocate the components of the header structured datatype except for hscr%hdr
SOURCE
1008 subroutine hscr_malloc(hscr, npwe, nqibz, nomega, nqlwl) 1009 1010 !Arguments ------------------------------------ 1011 !scalars 1012 class(hscr_t),intent(inout) :: Hscr 1013 integer,intent(in) :: npwe, nqibz, nomega, nqlwl 1014 1015 ! ************************************************************************* 1016 1017 !@hscr_t 1018 ABI_MALLOC(hscr%gvec, (3, npwe)) 1019 ABI_MALLOC(hscr%qibz, (3, nqibz)) 1020 ABI_MALLOC(hscr%qlwl, (3, nqlwl)) 1021 ABI_MALLOC(hscr%omega, (nomega)) 1022 1023 end subroutine hscr_malloc
m_io_screening/hscr_merge [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_merge
FUNCTION
This subroutine merges diffrent header structured variable (hscr_t)
INPUTS
Hscr_in(:) <hscr_t)>=List of headers to be merged.
OUTPUT
Hscr_out<hscr_t>=The output merged header.
SOURCE
1148 subroutine hscr_merge(Hscr_in, Hscr_out) 1149 1150 !Arguments ------------------------------------ 1151 !scalars 1152 type(hscr_t),intent(in) :: Hscr_in(:) 1153 type(hscr_t),intent(out) :: Hscr_out 1154 1155 !Local variables------------------------------- 1156 !scalars 1157 integer :: nhds,restart,restartpaw,ihd,ii,nqtot,nqneq 1158 logical :: isok 1159 character(len=500) :: msg 1160 !arrays 1161 real(dp),allocatable :: qset(:,:) 1162 ! ************************************************************************* 1163 1164 !@hscr_t 1165 nhds=SIZE(Hscr_in) 1166 1167 ! Initial copy of the header === 1168 ! If multiple headers, select the header containing q-->0 so that we copy also heads and wings 1169 ii = imax_loc(Hscr_in(:)%nqlwl) 1170 call hscr_copy(Hscr_in(ii),Hscr_out) 1171 if (nhds==1) return 1172 1173 ! Check consistency of the abinit Headers. 1174 ! FFT grid might be q-point dependent so we stop only when restart==0 1175 isok=.TRUE. 1176 do ihd=2,nhds 1177 call hdr_check(Hscr_in(1)%fform,Hscr_in(ihd)%fform,Hscr_in(1)%Hdr,Hscr_in(ihd)%Hdr,'COLL',restart,restartpaw) 1178 if (restart==0) then 1179 isok=.FALSE. 1180 write(msg,'(a,i0,a)')' Abinit header no.',ihd,' is not consistent with the first header ' 1181 ABI_WARNING(msg) 1182 end if 1183 end do 1184 if (.not.isok) then 1185 ABI_ERROR('Cannot continue, Check headers') 1186 end if 1187 1188 ! Now check variables related to polarizability|epsilon^{-1}. 1189 ! 1) Tests quantities that must be equal 1190 ii = assert_eq(Hscr_in(:)%ID, 'Headers have different Identifiers') 1191 ii = assert_eq(Hscr_in(:)%ikxc, 'Headers have different ikxc' ) 1192 ii = assert_eq(Hscr_in(:)%headform, 'Headers have different headform' ) 1193 ii = assert_eq(Hscr_in(:)%fform, 'Headers have different fform' ) 1194 ii = assert_eq(Hscr_in(:)%gwcalctyp,'Headers have different gwcalctyp' ) 1195 ii = assert_eq(Hscr_in(:)%nI, 'Headers have different nI' ) 1196 ii = assert_eq(Hscr_in(:)%nJ, 'Headers have different nJ' ) 1197 ii = assert_eq(Hscr_in(:)%nomega, 'Headers have different nomega' ) 1198 ii = assert_eq(Hscr_in(:)%test_type,'Headers have different test_type' ) 1199 ii = assert_eq(Hscr_in(:)%tordering,'Headers have different tordering' ) 1200 1201 ! This is not mandatory but makes life easier! 1202 ii = assert_eq(Hscr_in(:)%npwe,'Headers have different number of G-vectors' ) 1203 1204 do ihd=2,nhds 1205 if (ANY(ABS(Hscr_in(ihd)%omega-Hscr_in(1)%omega)>tol6)) then 1206 write(msg,'(a,i0,a)')' Frequencies in the first and the ',ihd,'-th header differ' 1207 ABI_ERROR(msg) 1208 end if 1209 if (ANY(Hscr_in(ihd)%gvec(:,:)-Hscr_in(1)%gvec(:,:)/=0)) then 1210 write(msg,'(a,i0,a)')' Incompatible G-vector list found in the ',ihd,'-th header' 1211 ABI_ERROR(msg) 1212 end if 1213 if (hscr_in(ihd)%kind_cdata /= hscr_in(1)%kind_cdata) then 1214 write(msg,'(3a,i0,2a)')' Files contain data with different precisions.',ch10,& 1215 "In particular the ",ihd,'-th header has precision:',trim(hscr_in(ihd)%kind_cdata) 1216 ABI_ERROR(msg) 1217 end if 1218 end do !ihd 1219 1220 ! If error is not fatal, just warn === 1221 if (ANY(Hscr_in(:)%npwwfn_used/=Hscr_in(1)%npwwfn_used)) then 1222 ABI_COMMENT('Files have been produced with a different number of planewaves for the wavefunctions.') 1223 end if 1224 if (ANY(Hscr_in(:)%nbnds_used/=Hscr_in(1)%nbnds_used)) then 1225 ABI_COMMENT('Files have been produced with a different number of bands.') 1226 end if 1227 if (ANY(Hscr_in(:)%spmeth/=Hscr_in(1)%spmeth)) then 1228 ABI_COMMENT('Files have been produced with different algorithms.') 1229 end if 1230 if (ANY(ABS(Hscr_in(:)%mbpt_sciss-Hscr_in(1)%mbpt_sciss)>tol6)) then 1231 ABI_COMMENT('Files have benn produced with different values of mbpt_sciss.') 1232 end if 1233 if (ANY(ABS(Hscr_in(:)%spsmear-Hscr_in(1)%spsmear)>tol6)) then 1234 ABI_COMMENT('Files have been produced with different values of spsmear.') 1235 end if 1236 if (ANY(ABS(Hscr_in(:)%zcut-Hscr_in(1)%zcut)>tol6)) then 1237 ABI_COMMENT('Files have been produced with different values of zcut.') 1238 end if 1239 1240 ! Now merge the list of q-points. 1241 ! Take the union of the q-points, remove possible duplicated 1242 ! are change the parameters in hscr_out that depends on q-points. 1243 nqtot=SUM(Hscr_in(:)%nqibz) 1244 ABI_MALLOC(qset,(3,nqtot)) 1245 1246 ii=0 1247 do ihd=1,nhds 1248 qset(:,ii+1:ii+Hscr_in(ihd)%nqibz)=Hscr_in(ihd)%qibz(:,:) 1249 ii=ii+Hscr_in(ihd)%nqibz 1250 end do 1251 1252 call remove_copies(nqtot,qset,nqneq,isequalk) 1253 1254 if (nqneq /= nqtot) then 1255 write(msg,'(3a,2(i0,a))')& 1256 'COMMENT: Headers contain duplicated q-points ',ch10,& 1257 'Found ',nqneq,' distinct q-points among the total ',nqtot,' points reported in the headers. ' 1258 call wrtout(std_out, msg) 1259 end if 1260 1261 Hscr_out%nqibz = nqneq 1262 ABI_FREE(Hscr_out%qibz) 1263 ABI_MALLOC(Hscr_out%qibz,(3,nqneq)) 1264 Hscr_out%qibz(:,:)=qset(:,1:nqneq) 1265 ABI_FREE(qset) 1266 1267 end subroutine hscr_merge
m_io_screening/hscr_mpio_skip [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_mpio_skip
FUNCTION
Skip the header of the (SCR|SUSC) file in MPI-IO mode. This routine uses local MPI-IO calls hence it can be safely called by master node only. Note however that in this case the offset has to be communicated to the other nodes.
INPUTS
mpio_fh=MPI-IO file handler fmarker_bsize = Byte length of Fortran record marker. fmarker_mpi_type= MPI type of the Fortran record marker
OUTPUT
fform=kind of the array in the file offset=The offset of the Fortran record located immediately below the Abinit header.
SOURCE
1710 subroutine hscr_mpio_skip(mpio_fh, fform, offset) 1711 1712 !Arguments ------------------------------------ 1713 integer,intent(in) :: mpio_fh 1714 integer,intent(out) :: fform 1715 integer(kind=XMPI_OFFSET_KIND),intent(out) :: offset 1716 1717 !Local variables------------------------------- 1718 !scalars 1719 integer :: bsize_frm,mpi_type_frm 1720 #ifdef HAVE_MPI_IO 1721 integer :: ierr,isk,nqlwl 1722 !character(len=500) :: msg 1723 !arrays 1724 integer(kind=MPI_OFFSET_KIND) :: fmarker,positloc 1725 integer :: statux(MPI_STATUS_SIZE) 1726 #endif 1727 1728 ! ************************************************************************* 1729 1730 offset = 0 1731 bsize_frm = xmpio_bsize_frm ! Byte size of the Fortran record marker. 1732 mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker. 1733 1734 call hdr_mpio_skip(mpio_fh,fform,offset) 1735 1736 call wrtout(std_out, sjoin("in hdr_mpio_skip with fform = ",itoa(fform))) 1737 1738 #ifdef HAVE_MPI_IO 1739 select case (fform) 1740 case (1003, 1004) 1741 ! Ship the titles 1742 call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr) 1743 1744 ! read nqlwl from the 2d record. 1745 positloc = offset + bsize_frm + 9*xmpi_bsize_int 1746 call MPI_FILE_READ_AT(mpio_fh,positloc,nqlwl,1,MPI_INTEGER,statux,ierr) 1747 call wrtout(std_out, sjoin("nqlwl = ",itoa(nqlwl))) 1748 1749 do isk=1,5 1750 call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr) 1751 end do 1752 1753 if (nqlwl>0) then ! skip qlwl 1754 call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr) 1755 end if 1756 1757 case default 1758 ABI_BUG(sjoin('Wrong fform read:', itoa(fform))) 1759 end select 1760 1761 #else 1762 ABI_ERROR("hscr_mpio_skip cannot be used when MPI-IO is not enabled") 1763 #endif 1764 1765 end subroutine hscr_mpio_skip
m_io_screening/hscr_new [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_new
FUNCTION
Initialize the Hscr datatype and most of its content from the em1params_t data type Ep.
INPUTS
varname=Name of the netcdf variable (used to get fform and ID). ikxc=Integer flag definining the type of XC kernel (0 if None i.e RPA) test_type=Integer flag defining the type of probing charge (0 for None) tordering=The time-ordering of the Response function. gvec(3,Ep%npwe)=The G-vectors used. Ep<em1params_t>=Parameters defining the calculation of the screening. hdr_abinit<hdr_type>=The abinit header.
OUTPUT
Hscr<type(hscr_t)>=the header, initialized.
SOURCE
823 type(hscr_t) function hscr_new(varname,dtset,ep,hdr_abinit,ikxc,test_type,tordering,titles,ngvec,gvec) result(hscr) 824 825 !Arguments ------------------------------------ 826 !scalars 827 integer,intent(in) :: ikxc,test_type,tordering,ngvec 828 character(len=*),intent(in) :: varname 829 type(dataset_type),intent(in) :: dtset 830 type(em1params_t),intent(in) :: Ep 831 type(hdr_type),intent(in) :: hdr_abinit 832 !arrays 833 integer,intent(in) :: gvec(3,ngvec) 834 character(len=80),intent(in) :: titles(2) 835 836 !Local variables------------------------------- 837 integer :: id 838 type(abifile_t) :: abifile 839 ! ************************************************************************* 840 841 !@hscr_t 842 ABI_CHECK(ngvec == Ep%npwe, 'ngvec/=Ep%npwe') 843 844 ! Identifier used to define the type of Response function (e^-1, chi0) 845 id = 0 846 if (varname == "polarizability") id = 1 847 if (varname == "inverse_dielectric_function") id = 4 848 ABI_CHECK(id /= 0, sjoin("Invalid varname: ",varname)) 849 850 ! Get fform from abifile. 851 abifile = abifile_from_varname(varname) 852 if (abifile%fform == 0) then 853 ABI_ERROR(sjoin("Cannot find any abifile object associated to varname:", varname)) 854 end if 855 856 ! Copy the abinit header. 857 call hdr_copy(hdr_abinit, Hscr%Hdr) 858 859 ! Initialize quantities related to the screening file 860 hscr%id =id 861 hscr%ikxc =ikxc 862 hscr%inclvkb =Ep%inclvkb 863 hscr%headform =HSCR_LATEST_HEADFORM 864 hscr%fform =abifile%fform 865 hscr%gwcalctyp =Ep%gwcalctyp 866 hscr%nI =Ep%nI 867 hscr%nJ =Ep%nJ 868 hscr%nqibz =Ep%nqcalc ! nqcalc == nqibz except if we split the calculation with nqptdm 869 hscr%nqlwl =Ep%nqlwl 870 hscr%nomega =Ep%nomega 871 hscr%nbnds_used =Ep%nbnds 872 hscr%npwe =Ep%npwe 873 hscr%npwwfn_used=Ep%npwwfn 874 hscr%spmeth =Ep%spmeth 875 hscr%test_type =test_type 876 hscr%tordering =tordering 877 hscr%mbpt_sciss =Ep%mbpt_sciss 878 hscr%spsmear =Ep%spsmear 879 hscr%zcut =Ep%zcut 880 881 hscr%titles(:)=titles(:) 882 883 call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl) 884 hscr%gvec(:,:) = gvec(1:3,1:Ep%npwe) 885 hscr%qibz(:,:) = Ep%qcalc 886 hscr%qlwl(:,:) = Ep%qlwl 887 hscr%omega(:) = Ep%omega 888 889 ! HSCR_NEW 890 hscr%awtr = dtset%awtr 891 hscr%icutcoul = dtset%gw_icutcoul 892 hscr%vcutgeo = dtset%vcutgeo 893 hscr%gwcomp = dtset%gwcomp 894 hscr%gwgamma = dtset%gwgamma 895 hscr%gwencomp = dtset%gwencomp 896 hscr%kind_cdata = "dpc" ! For the time being, data is always written in double-precision 897 ! HSCR_NEW 898 899 end function hscr_new
m_io_screening/hscr_print [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
hscr_print
FUNCTION
Prints info on the header of the SCR|SUSC file.
INPUTS
OUTPUT
SOURCE
707 subroutine hscr_print(Hscr, header, unit, prtvol, mode_paral) 708 709 !Arguments ------------------------------------ 710 !scalars 711 class(hscr_t),intent(in) :: hscr 712 integer,intent(in),optional :: prtvol,unit 713 character(len=4),intent(in),optional :: mode_paral 714 character(len=*),intent(in),optional :: header 715 716 !Local variables------------------------------- 717 !scalars 718 integer :: iomega,iqibz,unt,verbose 719 character(len=4) :: mode 720 character(len=500) :: msg 721 722 ! ************************************************************************* 723 724 unt=std_out; if (PRESENT(unit )) unt =unit 725 verbose=0 ; if (PRESENT(prtvol )) verbose=prtvol 726 mode='COLL'; if (PRESENT(mode_paral)) mode =mode_paral 727 728 if (PRESENT(header)) then 729 msg=' ==== '//TRIM(ADJUSTL(header))//' ==== ' 730 call wrtout(unt,msg,mode) 731 end if 732 733 write(msg,'(1x,a)')TRIM(hscr%titles(1)) 734 call wrtout(unt,msg,mode) 735 write(msg,'(1x,a)')TRIM(hscr%titles(2)) 736 call wrtout(unt,msg,mode) 737 write(msg,'(a,i8)') ' Identifier ',hscr%ID 738 call wrtout(unt,msg,mode) 739 write(msg,'(a,i8)') ' Kxc kernel ',hscr%ikxc 740 call wrtout(unt,msg,mode) 741 write(msg,'(a,i8)') ' Treatment of q-->0 limit ',hscr%inclvkb 742 call wrtout(unt,msg,mode) 743 write(msg,'(a,i8)') ' headform ',hscr%headform 744 call wrtout(unt,msg,mode) 745 write(msg,'(a,i8)') ' fform ',hscr%fform 746 call wrtout(unt,msg,mode) 747 write(msg,'(a,i8)') ' gwcalctyp ',hscr%gwcalctyp 748 call wrtout(unt,msg,mode) 749 write(msg,'(a,2i8)')' Number of components ',hscr%nI,hscr%nJ 750 call wrtout(unt,msg,mode) 751 write(msg,'(a,i8)') ' Number of q-points ',hscr%nqibz 752 call wrtout(unt,msg,mode) 753 write(msg,'(a,i8)') ' Number of q-directions ',hscr%nqlwl 754 call wrtout(unt,msg,mode) 755 write(msg,'(a,i8)') ' Number of frequencies ',hscr%nomega 756 call wrtout(unt,msg,mode) 757 write(msg,'(a,i8)') ' Number of bands used ',hscr%nbnds_used 758 call wrtout(unt,msg,mode) 759 write(msg,'(a,i8)') ' Dimension of matrix ',hscr%npwe 760 call wrtout(unt,msg,mode) 761 write(msg,'(a,i8)') ' Number of planewaves used ',hscr%npwwfn_used 762 call wrtout(unt,msg,mode) 763 write(msg,'(a,i8)') ' Spectral method ',hscr%spmeth 764 call wrtout(unt,msg,mode) 765 write(msg,'(a,i8)') ' Test_type ',hscr%test_type 766 call wrtout(unt,msg,mode) 767 write(msg,'(a,i8)') ' Time-ordering ',hscr%tordering 768 call wrtout(unt,msg,mode) 769 write(msg,'(a,es16.6)')' Scissor Energy ',hscr%mbpt_sciss 770 call wrtout(unt,msg,mode) 771 write(msg,'(a,es16.6)')' Spectral smearing ',hscr%spsmear 772 call wrtout(unt,msg,mode) 773 write(msg,'(a,es16.6)')' Complex Imaginary Shift ',hscr%zcut 774 call wrtout(unt,msg,mode) 775 776 if (verbose==0) then 777 call wrtout(unt,' The header contains additional records.',mode) 778 else 779 write(msg,'(2a)')ch10,' q-points [r.l.u.]:' 780 call wrtout(unt,msg,mode) 781 do iqibz=1,hscr%nqibz 782 write(msg,'(i5,3f12.6)')iqibz,hscr%qibz(:,iqibz) 783 call wrtout(unt,msg,mode) 784 end do 785 786 write(msg,'(2a)')ch10,' Frequencies used [eV]:' 787 call wrtout(unt,msg,mode) 788 do iomega=1,hscr%nomega 789 write(msg,'(i3,2f7.2)')iomega,REAL(hscr%omega(iomega))*Ha_eV,AIMAG(hscr%omega(iomega))*Ha_eV 790 call wrtout(unt,msg,mode) 791 end do 792 end if 793 794 ! Echo the abinit header. 795 !if (prtvol>0) call hdr_echo(hscr%hdr,fform,rdwr,unit=unt) 796 797 end subroutine hscr_print
m_io_screening/hscr_t [ Types ]
[ Top ] [ m_io_screening ] [ Types ]
NAME
hscr_t
FUNCTION
The structure defining the header of the SCR/SUSC file. hscr_t contains the most important dimensions associated to the SCR/SUSC matrix, important GW metadata and the Abint header. The SCR/SUS matrices are saved with the same format. There are nqibz blocks, each block contains (npwe,npwe,nomega) matrices. SCR and SUS files mainly differ for what concerns the treatment of the q-->0 limit. The treatment of the non-analytic behaviour is not yet implemented but the main ideas are sketched in the NOTES below.
NOTES
On the treatment of the q-->0 limit 1) q=Gamma should be the first q-point 2) This point contains an initial section with data used to treat the q-->0 limit, followed by the SCR/SUS matrix evaluated for the small q-point (qlwl). The header should contains enough info so that we can skip this section and use the pre-existing routines to read the matrices. 3) The data stored in the q-->0 section depends on the type of matrix stored in the file: SUS file: we store the tensor and the wings needed to reconstruct the q-dependence around Gamma. This section if followed by the X(G,G') matrix evaluated at qlwl. Note that the content of the SUS file does not depend on a possible cutoff in vcoul. SCR file: two cases must be considered: No cutoff in vcoul: The q-->0 section stores the tensor, the wings as well as the inverse of the body B^{-1}. This info is used to integrate W(q) around q==Gamma. cutoff in vcoul: The content of the q-->0 section depends on dimensionality of the system. In 2D we need info on chi0 as well as A and a. See http://arxiv.org/pdf/1511.00129v1.pdf I think that this kind of calculations are easy to implement if we start from the SUS file. Ok, we have to recompute e-1 at each run but the logic is easier to implement.
SOURCE
105 type,public :: hscr_t 106 107 integer :: id = -1 108 ! Matrix identifier: 1 for chi0, 2 for chi, 3 for epsilon, 4 for espilon^{-1} 109 110 integer :: ikxc = 0 111 ! Kxc kernel used, 112 ! 0 for None (RPA), >0 for static TDDFT (=ixc), <0 for frequency-dependent TDDFT 113 114 integer :: inclvkb = 2 115 ! q-->0 treatment, 0 for None, 1-2 for transversal gauge, 3 for longitudinal 116 117 integer :: headform 118 ! format of the SCR header 119 120 integer :: fform 121 ! File format 122 123 integer :: gwcalctyp = 0 124 ! Calculation type (G0W0, G0W, GW ...) 125 126 integer :: nI = 1, nJ = 1 127 ! Number of spin components (rows,columns) in chi|eps^-1. (1,1) if collinear. 128 ! The internal representation of the matrix is eps(nI*npwe,nJ*npwe) 129 130 integer :: nqibz = -1 131 ! Number of q-points in the IBZ. 132 133 integer :: nqlwl = 1 134 ! Number of points for the treatment of the long wavelength limit. 135 136 integer :: nomega = -1 137 ! Total number of frequencies. 138 139 integer :: nbnds_used = -1 140 ! Number of bands used during the screening calculation (only for info) 141 142 integer :: npwe = -1 143 ! Number of G vectors reported on the file. 144 145 integer :: npwwfn_used = -1 146 ! Number of G vectors for wavefunctions used during the screening calculation (only for info) 147 148 integer :: spmeth = 0 149 ! Method used to approximate the delta function in the expression for Im Chi_0 150 151 integer :: test_type 152 ! 1 for TEST-PARTICLE, 2 for TEST-ELECTRON. 153 154 integer :: tordering = 1 155 ! 1 for Time-Ordered, 2 for Advanced, 3 for Retarded. 156 157 ! HSCR_NEW 158 integer :: awtr = 1 159 ! Input variable (time-reversal symmetry in RPA expression) 160 161 integer :: icutcoul = 0 162 ! Input variable (Coulomb singularity treatment) 163 164 integer :: gwcomp = 0 165 ! Input variable (GW compensation energy technique) 166 167 integer :: gwgamma = 0 168 ! Input variable Vertex correction 169 ! HSCR_NEW 170 171 real(dp) :: mbpt_sciss = zero 172 ! Scissor Energy, zero if not used 173 174 real(dp) :: spsmear = zero 175 ! Smearing of the delta in case of spmeth==2 176 177 real(dp) :: zcut = -one 178 ! Imaginary shift to avoid the poles along the real axis. 179 180 ! HSCR_NEW 181 real(dp) :: gwencomp = -one 182 ! Input variable (GW compensation energy technique) 183 184 character(len=3) :: kind_cdata 185 ! Flag to signal whether the data is in single or double precision ("spc" or "dpc") 186 ! For the time being, we always write/read in double precision. 187 ! This flag could be use to reduce the memory requirements if spc: 188 ! we run calculations in single precision dump the results with the same precision without 189 ! having to allocate extra memory. 190 ! HSCR_NEW 191 192 !arrays 193 194 ! HSCR_NEW 195 real(dp) :: vcutgeo(3) = zero 196 ! Input variable (defines coulomb cutoff) 197 ! HSCR_NEW 198 199 character(len=80) :: titles(2) 200 ! Titles describing the content of the file. 201 202 integer,allocatable :: gvec(:,:) 203 ! gvec(3,npwe) 204 ! G vectors in reduced coordinates. 205 206 real(dp),allocatable :: qibz(:,:) 207 ! qibz(3,nqibz) 208 ! q-points in the IBZ in reduced coordinates. 209 210 real(dp),allocatable :: qlwl(:,:) 211 ! qlwl(3,nqlwl) 212 ! q-points for the long wave-length limit treatment (r.l.u) 213 214 complex(dpc),allocatable :: omega(:) 215 ! omega(nomega) 216 ! All frequencies calculated both along the real and the imaginary axis. 217 ! Real frequencies are packed in the first section. 218 ! TODO: Add frequency mesh type? 219 220 type(hdr_type) :: hdr 221 ! The abinit header. 222 223 contains 224 225 procedure :: from_file => hscr_from_file ! Read the header from file. 226 procedure :: print => hscr_print ! Print the SCR-related part of the header. 227 procedure :: bcast => hscr_bcast ! Broadcast the header. 228 procedure :: free => hscr_free ! Free the header. 229 procedure :: copy => hscr_copy ! Copy the SCR|SUSC header. 230 231 end type hscr_t
m_io_screening/ioscr_qmerge [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
ioscr_qmerge
FUNCTION
Produce new file by merging the q-points stored in other files. This routine should be called by a single MPI process.
INPUTS
nfiles=Number of files to be merged. filenames(nfiles)=Paths of files to be merged. hscr_files(nfiles)<hscr_t>=Heades of the files to be merged. fname_out=Name of the file to be produced.
OUTPUT
ohscr<hscr_t>=The header of the output file.
SOURCE
1789 subroutine ioscr_qmerge(nfiles, filenames, hscr_files, fname_out, ohscr) 1790 1791 !Arguments ------------------------------------ 1792 !scalars 1793 integer,intent(in) :: nfiles 1794 character(len=*),intent(in) :: fname_out 1795 type(hscr_t),intent(out) :: ohscr 1796 !arrays 1797 character(len=*),intent(in) :: filenames(nfiles) 1798 type(hscr_t),intent(in) :: hscr_files(nfiles) 1799 1800 !Local variables------------------------------- 1801 !scalars 1802 integer,parameter :: rdwr2=2,master=0 1803 integer :: iqibz,ifound,ifile,iqf,ount,iomode,fform_merge,comm,nomega4m,npwe4m,iqiA,ierr 1804 character(len=500) :: msg 1805 character(len=nctk_slen) :: varname 1806 type(abifile_t) :: abifile 1807 !arrays 1808 integer,allocatable :: merge_table(:,:) 1809 real(dp) :: qdiff(3) 1810 complex(gwpc),allocatable :: epsm1(:,:,:,:) 1811 1812 ! ************************************************************************* 1813 1814 comm = xmpi_comm_self 1815 1816 if (file_exists(fname_out)) then 1817 ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out)) 1818 end if 1819 1820 ! Merge the headers creating the full list of q-points. 1821 call hscr_merge(Hscr_files(1:nfiles), ohscr) 1822 call hscr_print(ohscr, header='Header of the final file', unit=std_out, prtvol=1) 1823 1824 ! For each q to be merged, save the index of the file where q is stored as well as its sequential index. 1825 ! Useful to do the merge point-by-point thus avoiding the allocation of the entire epsm1 array. 1826 ABI_MALLOC(merge_table,(ohscr%nqibz,2)) 1827 do iqibz=1,ohscr%nqibz 1828 ifound=0 1829 fl: do ifile=1,nfiles 1830 do iqf=1,Hscr_files(ifile)%nqibz 1831 qdiff(:)=ohscr%qibz(:,iqibz)-Hscr_files(ifile)%qibz(:,iqf) 1832 if (all(abs(qdiff) < GW_TOLQ)) then 1833 merge_table(iqibz,1)=ifile 1834 merge_table(iqibz,2)=iqf 1835 ifound=ifound+1 1836 write(msg,'(a,3f12.6,2a)')'. q-point:',ohscr%qibz(:,iqibz),' will be taken from ',TRIM(filenames(ifile)) 1837 call wrtout(std_out, msg) 1838 EXIT fl 1839 end if 1840 end do 1841 end do fl 1842 ! Check if q-point has been found, multiple q-points not allowed. 1843 ABI_CHECK(ifound == 1, 'ifound/=1') 1844 end do 1845 1846 iomode = IO_MODE_FORTRAN; if (endswith(fname_out, ".nc")) iomode = IO_MODE_ETSF 1847 if (iomode == IO_MODE_FORTRAN) then 1848 if (open_file(fname_out,msg,newunit=ount,status='new',form='unformatted') /= 0) then 1849 ABI_ERROR(msg) 1850 end if 1851 else 1852 NCF_CHECK(nctk_open_create(ount, fname_out, comm)) 1853 end if 1854 1855 ! Write the header. 1856 fform_merge = hscr_files(1)%fform 1857 abifile = abifile_from_fform(fform_merge) 1858 if (abifile%fform == 0) then 1859 ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge))) 1860 end if 1861 varname = abifile%varname 1862 1863 if (any(hscr_files(:)%fform /= hscr_files(1)%fform)) then 1864 write(std_out,*)"fforms: ",hscr_files(:)%fform 1865 ABI_ERROR("Files to be merged have different fform. Cannot merge data") 1866 end if 1867 1868 call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode) 1869 1870 npwe4m = ohscr%npwe 1871 nomega4m = ohscr%nomega 1872 1873 ABI_MALLOC_OR_DIE(epsm1,(npwe4m,npwe4m,nomega4m,1), ierr) 1874 1875 do iqibz=1,ohscr%nqibz 1876 ifile = merge_table(iqibz,1) 1877 iqiA = merge_table(iqibz,2) 1878 call read_screening(varname,filenames(ifile),npwe4m,1,nomega4m,epsm1,iomode,comm,iqiA=iqiA) 1879 call write_screening(varname,ount,iomode,npwe4m,nomega4m,iqibz,epsm1) 1880 end do 1881 1882 ABI_FREE(epsm1) 1883 ABI_FREE(merge_table) 1884 1885 if (iomode == IO_MODE_FORTRAN) then 1886 close(ount) 1887 else 1888 NCF_CHECK(nf90_close(ount)) 1889 end if 1890 1891 write(msg,'(3a)')ch10,' ==== Files have been merged successfully === ',ch10 1892 call wrtout(std_out, msg) 1893 1894 end subroutine ioscr_qmerge
m_io_screening/ioscr_qrecover [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
ioscr_qrecover
FUNCTION
Recover q-points from a corrupted file produced e.g. from an interrupted run This routine should be called by a single MPI process.
INPUTS
path=Corrupted file. nqrec=Number of q-points to recover. fname_out=Name of the file to be produced.
OUTPUT
Output is written to file.
SOURCE
1917 subroutine ioscr_qrecover(ipath, nqrec, fname_out) 1918 1919 !Arguments ------------------------------------ 1920 !scalars 1921 integer,intent(in) :: nqrec 1922 character(len=*),intent(in) :: ipath,fname_out 1923 1924 !Local variables------------------------------- 1925 !scalars 1926 integer,parameter :: rdwr2=2,master=0 1927 integer :: iqiA,nqibzA,nomega_asked,unt,npwe_asked,iomode,comm,fform1,ifform,ierr 1928 character(len=500) :: msg 1929 character(len=nctk_slen) :: varname 1930 type(hscr_t) :: hscr_recov,hscr 1931 type(abifile_t) :: abifile 1932 !arrays 1933 complex(gwpc),allocatable :: epsm1(:,:,:,:) 1934 1935 ! ************************************************************************* 1936 1937 comm = xmpi_comm_self 1938 1939 call wrtout(std_out, sjoin(". Recovering q-points in file:", ipath)) 1940 call wrtout(std_out, sjoin(". Data written to file:", fname_out)) 1941 1942 if (file_exists(fname_out)) then 1943 ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out)) 1944 end if 1945 1946 ! Find iomode from file extension and open output file. 1947 if (endswith(fname_out, ".nc")) then 1948 iomode = IO_MODE_ETSF 1949 NCF_CHECK(nctk_open_create(unt, fname_out, comm)) 1950 else 1951 iomode = IO_MODE_FORTRAN 1952 if (open_file(fname_out, msg, newunit=unt, status='new', form='unformatted') /= 0) then 1953 ABI_ERROR(msg) 1954 end if 1955 end if 1956 1957 ! Read header. 1958 call hscr_from_file(hscr, ipath, ifform, comm) 1959 ABI_CHECK(ifform /= 0, sjoin("fform = 0 while reading:", ipath)) 1960 1961 if (nqrec < 1 .or. nqrec > hscr%nqibz) then 1962 ABI_ERROR(sjoin("Wrong input. nqibz on file:", itoa(hscr%nqibz))) 1963 end if 1964 1965 ! Copy header 1966 call hscr_copy(hscr, hscr_recov) 1967 1968 ! Change dimensions and arrays associated to nqibz. 1969 hscr_recov%nqibz = nqrec 1970 ABI_FREE(hscr_recov%qibz) 1971 ABI_MALLOC(hscr_recov%qibz, (3,nqrec)) 1972 hscr_recov%qibz = hscr%qibz(:,1:nqrec) 1973 1974 call hscr_print(hscr_recov,header="Header of the new SCR file",unit=std_out,prtvol=1) 1975 1976 ! Write the header of the recovered file. 1977 fform1 = hscr%fform 1978 1979 abifile = abifile_from_fform(fform1) 1980 if (abifile%fform == 0) then 1981 ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform1))) 1982 end if 1983 varname = abifile%varname 1984 1985 call hscr_io(hscr_recov,fform1,rdwr2,unt,comm,master,iomode) 1986 1987 nqibzA=1; nomega_asked=hscr%nomega; npwe_asked=hscr%npwe 1988 1989 ABI_MALLOC_OR_DIE(epsm1,(npwe_asked,npwe_asked,nomega_asked,1), ierr) 1990 1991 do iqiA=1,hscr_recov%nqibz 1992 call read_screening(varname,ipath,npwe_asked,nqibzA,nomega_asked,epsm1,iomode,comm,iqiA=iqiA) 1993 call write_screening(varname,unt,iomode,npwe_asked,nomega_asked,iqiA,epsm1) 1994 end do 1995 1996 if (iomode == IO_MODE_FORTRAN) close(unt) 1997 if (iomode == IO_MODE_ETSF) then 1998 NCF_CHECK(nf90_close(unt)) 1999 end if 2000 2001 ABI_FREE(epsm1) 2002 call hscr%free() 2003 call hscr_recov%free() 2004 2005 call wrtout(std_out, "Recovery completed") 2006 2007 end subroutine ioscr_qrecover
m_io_screening/ioscr_wmerge [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
ioscr_wmerge
FUNCTION
Produce new file by merging the frequencies stored in other files. This routine should be called by a single MPI process.
INPUTS
nfiles=Number of files to be merged. filenames(nfiles)=Paths of files to be merged. hscr_files(nfiles)<hscr_t>=Heades of the files to be merged. fname_out=Name of the file to be produced.
OUTPUT
ohscr<hscr_t>=The header of the output file.
SOURCE
2031 subroutine ioscr_wmerge(nfiles, filenames, hscr_file, freqremax, fname_out, ohscr) 2032 2033 !Arguments ------------------------------------ 2034 !scalars 2035 integer,intent(in) :: nfiles 2036 real(dp),intent(in) :: freqremax 2037 character(len=*),intent(in) :: fname_out 2038 type(hscr_t),intent(out) :: ohscr 2039 !arrays 2040 character(len=*),intent(in) :: filenames(nfiles) 2041 type(hscr_t),intent(in) :: hscr_file(nfiles) 2042 2043 !Local variables------------------------------- 2044 !scalars 2045 integer,parameter :: rdwr2=2,master=0 2046 integer :: ii,iqibz,ifile,ount,iomode,fform_merge,comm,nomega4m 2047 integer :: nfreq_tot,nfreqre,nfreqim,ifrq,npwe4mI,npwe4mJ,ierr 2048 character(len=500) :: msg 2049 logical :: skip 2050 character(len=nctk_slen) :: varname 2051 type(abifile_t) :: abifile 2052 !arrays 2053 integer,allocatable :: freq_indx(:,:),ifile_indx(:),pos_indx(:),i_temp(:),i2_temp(:,:) 2054 real(dp),allocatable :: real_omega(:),imag_omega(:) 2055 complex(gwpc),allocatable :: epsm1(:,:,:,:),epsm1_temp(:,:,:,:) 2056 complex(dpc),allocatable :: omega_storage(:) 2057 2058 ! ************************************************************************* 2059 2060 comm = xmpi_comm_self 2061 2062 ! Check that q-point sets are the same 2063 do ifile=1,nfiles 2064 do ii=1,nfiles 2065 if (ii == ifile) CYCLE 2066 if (Hscr_file(ifile)%nqibz /= Hscr_file(ii)%nqibz) then 2067 ABI_ERROR(' One or more files do not have the same number of q-points!') 2068 end if 2069 do iqibz=1,Hscr_file(1)%nqibz 2070 if (ABS(SUM(Hscr_file(ifile)%qibz(:,iqibz)-Hscr_file(ii)%qibz(:,iqibz))) > tol6) then 2071 ABI_ERROR('Q-point set differs between one or more files!') 2072 end if 2073 end do 2074 end do 2075 end do 2076 2077 ! nfreq_tot here is the total *possible* number of freq. 2078 nfreq_tot=0 2079 do ifile=1,nfiles 2080 nfreq_tot = nfreq_tot + Hscr_file(ifile)%nomega 2081 end do 2082 2083 ABI_MALLOC(omega_storage, (nfreq_tot)) 2084 ABI_MALLOC(freq_indx, (nfreq_tot,nfiles)) 2085 ABI_MALLOC(ifile_indx, (nfreq_tot)) 2086 omega_storage = CMPLX(-one,-one); freq_indx = 0; ifile_indx = 0 2087 2088 ! Calculate the total number of real freq and store 2089 nfreqre = 0 2090 do ifile=1,nfiles 2091 do ifrq=1,Hscr_file(ifile)%nomega 2092 skip = .FALSE. 2093 ! Check whether to skip this point 2094 if (AIMAG(Hscr_file(ifile)%omega(ifrq)) > tol16) skip = .TRUE. 2095 if (REAL(Hscr_file(ifile)%omega(ifrq)) > freqremax) skip = .TRUE. 2096 ! Check for repetition or non-monotonic points 2097 if (nfreqre>1) then 2098 do ii=1,nfreqre 2099 if (ABS(REAL(Hscr_file(ifile)%omega(ifrq)) - REAL(omega_storage(ii))) < tol6) skip = .TRUE. 2100 end do 2101 end if 2102 if (skip) CYCLE 2103 nfreqre = nfreqre + 1 2104 ! Store (complex) frequency and index 2105 omega_storage(nfreqre) = Hscr_file(ifile)%omega(ifrq) 2106 ifile_indx(nfreqre) = ifile 2107 freq_indx(nfreqre,ifile) = ifrq 2108 write(std_out,'(a,es16.6,a,i0,2(a,i0))')& 2109 ' Found real frequency: ',REAL(omega_storage(nfreqre))*Ha_eV,' [eV], number: ',nfreqre,& 2110 ', in file: ',ifile,' local index: ',ifrq 2111 end do 2112 end do 2113 2114 if (nfreqre > 0) then 2115 ! Sort real frequencies 2116 ABI_MALLOC(real_omega,(nfreqre)) 2117 ABI_MALLOC(pos_indx,(nfreqre)) 2118 ABI_MALLOC(i_temp,(nfreqre)) 2119 ABI_MALLOC(i2_temp,(nfreqre,nfiles)) 2120 real_omega(1:nfreqre) = REAL(omega_storage(1:nfreqre)) ! Copy real frequencies to temp. sorting array 2121 do ii=1,nfreqre ! Set up indexing array 2122 pos_indx(ii) = ii 2123 end do 2124 2125 ! Sort frequencies while keeping track of index 2126 call sort_dp(nfreqre,real_omega,pos_indx,tol16) 2127 i_temp(1:nfreqre) = ifile_indx(1:nfreqre) 2128 i2_temp(1:nfreqre,1:nfiles) = freq_indx(1:nfreqre,1:nfiles) 2129 2130 ! Copy sorted frequencies plus file and frequency index 2131 do ii=1,nfreqre 2132 omega_storage(ii) = CMPLX(real_omega(ii),zero) 2133 ifile_indx(ii) = i_temp(pos_indx(ii)) 2134 freq_indx(ii,1:nfiles) = i2_temp(pos_indx(ii),1:nfiles) 2135 end do 2136 2137 ABI_FREE(real_omega) 2138 ABI_FREE(pos_indx) 2139 ABI_FREE(i_temp) 2140 ABI_FREE(i2_temp) 2141 end if 2142 2143 ! Check imaginary frequencies and store them 2144 nfreqim = 0 2145 do ifile=1,nfiles 2146 do ifrq=1,Hscr_file(ifile)%nomega 2147 if (REAL(Hscr_file(ifile)%omega(ifrq)) > tol8) CYCLE 2148 if (AIMAG(Hscr_file(ifile)%omega(ifrq)) < tol8) CYCLE 2149 nfreqim = nfreqim + 1 2150 omega_storage(nfreqre+nfreqim) = Hscr_file(ifile)%omega(ifrq) 2151 ifile_indx(nfreqre+nfreqim) = ifile 2152 freq_indx(nfreqre+nfreqim,ifile) = ifrq 2153 write(std_out,'(a,es16.6,a,i0,2(a,i0))')& 2154 ' Found imaginary frequency: ',AIMAG(omega_storage(nfreqre+nfreqim))*Ha_eV,' [eV], number: ',nfreqim,& 2155 ', in file: ',ifile,' local index: ',ifrq 2156 end do 2157 end do 2158 2159 ! Sort imaginary frequencies 2160 ABI_MALLOC(imag_omega,(nfreqim)) 2161 ABI_MALLOC(pos_indx,(nfreqim)) 2162 ABI_MALLOC(i_temp,(nfreqim)) 2163 ABI_MALLOC(i2_temp,(nfreqim,nfiles)) 2164 2165 ! Copy imaginary frequencies to temp. sorting array 2166 imag_omega(1:nfreqim) = AIMAG(omega_storage(nfreqre+1:nfreqre+nfreqim)) 2167 do ii=1,nfreqim ! Set up indexing array 2168 pos_indx(ii) = ii 2169 end do 2170 2171 ! Sort frequencies while keeping track of index 2172 call sort_dp(nfreqim,imag_omega,pos_indx,tol16) 2173 i_temp(1:nfreqim) = ifile_indx(nfreqre+1:nfreqre+nfreqim) 2174 i2_temp(1:nfreqim,1:nfiles) = freq_indx(nfreqre+1:nfreqre+nfreqim,1:nfiles) 2175 2176 ! Copy sorted frequencies plus file and frequency index 2177 do ii=1,nfreqim 2178 omega_storage(nfreqre+ii) = CMPLX(zero,imag_omega(ii)) 2179 ifile_indx(nfreqre+ii) = i_temp(pos_indx(ii)) 2180 freq_indx(nfreqre+ii,1:nfiles) = i2_temp(pos_indx(ii),1:nfiles) 2181 end do 2182 2183 ABI_FREE(imag_omega) 2184 ABI_FREE(pos_indx) 2185 ABI_FREE(i_temp) 2186 ABI_FREE(i2_temp) 2187 2188 nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq 2189 write(std_out,'(2a,i0,a)') ch10,' Merging ',nfreq_tot,' frequencies.' 2190 write(std_out,'(2(a,i0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10 2191 2192 ! Copy old header 2193 call hscr_copy(Hscr_file(1),ohscr) 2194 2195 ! TODO: hscr_wmerge 2196 ! Then modify entries for new frequency grid 2197 ohscr%nomega = nfreq_tot 2198 ABI_FREE(ohscr%omega) 2199 ABI_MALLOC(ohscr%omega,(nfreq_tot)) 2200 ohscr%omega(1:nfreq_tot) = omega_storage(1:nfreq_tot) 2201 2202 npwe4mI = ohscr%npwe*ohscr%nI 2203 npwe4mJ = ohscr%npwe*ohscr%nJ 2204 2205 ! Print new header for info 2206 call hscr_print(ohscr,header='Header of the final file',unit=std_out,prtvol=1) 2207 2208 if (file_exists(fname_out)) then 2209 ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out)) 2210 end if 2211 2212 if (endswith(fname_out, ".nc")) then 2213 iomode = IO_MODE_ETSF 2214 NCF_CHECK(nctk_open_create(ount, fname_out, comm)) 2215 else 2216 iomode = IO_MODE_FORTRAN 2217 if (open_file(fname_out, msg, newunit=ount, status='new',form='unformatted') /= 0) then 2218 ABI_ERROR(msg) 2219 end if 2220 end if 2221 2222 ! * Write the header. 2223 fform_merge = ohscr%fform 2224 2225 abifile = abifile_from_fform(fform_merge) 2226 if (abifile%fform == 0) then 2227 ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge))) 2228 end if 2229 varname = abifile%varname 2230 2231 call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode) 2232 2233 npwe4mI = ohscr%npwe*ohscr%nI 2234 npwe4mJ = ohscr%npwe*ohscr%nJ 2235 nomega4m = ohscr%nomega 2236 ABI_MALLOC_OR_DIE(epsm1,(npwe4mI,npwe4mJ,nomega4m,1), ierr) 2237 2238 do iqibz=1,ohscr%nqibz 2239 2240 do ifile=1,nfiles 2241 ! allocate temporary array 2242 npwe4mI = Hscr_file(ifile)%npwe*Hscr_file(ifile)%nI 2243 npwe4mJ = Hscr_file(ifile)%npwe*Hscr_file(ifile)%nJ 2244 nomega4m = Hscr_file(ifile)%nomega 2245 ABI_MALLOC_OR_DIE(epsm1_temp,(npwe4mI,npwe4mJ,nomega4m,1), ierr) 2246 2247 ! read screening 2248 call read_screening(varname,filenames(ifile),npwe4mI,1,nomega4m,epsm1_temp,iomode,comm,iqiA=iqibz) 2249 2250 ! Copy matrices for relevant frequencies 2251 do ifrq=1,nfreq_tot 2252 if (ifile_indx(ifrq)==ifile) then 2253 epsm1(:,:,ifrq,1)=epsm1_temp(:,:,freq_indx(ifrq,ifile),1) 2254 end if 2255 end do 2256 2257 !Add imaginary part if needed 2258 !if (indx_imfreq_file==ifile) then 2259 !do ifrq=nfreqre+1,ohscr%nomega 2260 !epsm1(:,:,ifrq,1)=epsm1_temp(:,:,freq_indx(ifrq,ifile),1) 2261 !end do 2262 !end if 2263 2264 ABI_FREE(epsm1_temp) 2265 end do !ifile 2266 2267 ! Write data. 2268 npwe4mI = ohscr%npwe*ohscr%nI 2269 nomega4m = ohscr%nomega 2270 call write_screening(varname,ount,iomode,npwe4mI,nomega4m,iqibz,epsm1) 2271 end do ! iqibz 2272 2273 ABI_FREE(epsm1) 2274 ABI_FREE(omega_storage) 2275 ABI_FREE(freq_indx) 2276 ABI_FREE(ifile_indx) 2277 2278 if (iomode == IO_MODE_FORTRAN) then 2279 close(ount) 2280 else 2281 NCF_CHECK(nf90_close(ount)) 2282 end if 2283 2284 write(msg,'(3a)')ch10,' ==== Files have been merged successfully === ',ch10 2285 call wrtout(std_out, msg) 2286 2287 end subroutine ioscr_wmerge
m_io_screening/ioscr_wremove [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
ioscr_wremove
FUNCTION
Produce new file by removing selected frequencies in the initial file `inpath`. This routine should be called by a single MPI process.
INPUTS
inpath=Input file ihscr<hscr_t>=Headerf of the input file. fname_out=Output file. nfreq_tot=Number of frequencies in new file. freq_indx(nfreq_tot)=Index of frequency to be kept in input file.
OUTPUT
ohscr<hscr_t>=The header of the output file.
SOURCE
2312 subroutine ioscr_wremove(inpath, ihscr, fname_out, nfreq_tot, freq_indx, ohscr) 2313 2314 !Arguments ------------------------------------ 2315 !scalars 2316 integer,intent(in) :: nfreq_tot 2317 character(len=*),intent(in) :: inpath,fname_out 2318 type(hscr_t),intent(in) :: ihscr 2319 type(hscr_t),intent(out) :: ohscr 2320 !arrays 2321 integer,intent(in) :: freq_indx(nfreq_tot) 2322 2323 !Local variables------------------------------- 2324 !scalars 2325 integer,parameter :: rdwr2=2,master=0 2326 integer :: iqibz,fform_merge,comm,nomega4m,ierr 2327 integer :: ifrq,npwe4mI,npwe4mJ,iomode,ount 2328 character(len=500) :: msg 2329 character(len=nctk_slen) :: varname 2330 type(abifile_t) :: abifile 2331 !arrays 2332 complex(gwpc),allocatable :: epsm1(:,:,:),epsm1_temp(:,:,:) 2333 2334 ! ************************************************************************* 2335 2336 comm = xmpi_comm_self 2337 2338 ! check ifreq_idx 2339 ABI_CHECK(nfreq_tot > 0, "nfreq_tot <= 0!") 2340 if (all(freq_indx == 0)) ABI_ERROR("all(freq_indx == 0)") 2341 2342 ! Copy the old header 2343 call hscr_copy(ihscr, ohscr) 2344 2345 ! Then modify entries for new frequency grid. 2346 ohscr%nomega = nfreq_tot 2347 ABI_FREE(ohscr%omega) 2348 ABI_MALLOC(ohscr%omega,(nfreq_tot)) 2349 do ifrq=1,nfreq_tot 2350 ohscr%omega(ifrq) = ihscr%omega(freq_indx(ifrq)) 2351 end do 2352 2353 npwe4mI = ohscr%npwe*ohscr%nI 2354 npwe4mJ = ohscr%npwe*ohscr%nJ 2355 2356 ! Print new header for info 2357 call hscr_print(ohscr,header='Header of the final file',unit=std_out,prtvol=1) 2358 2359 ! Open output file. 2360 if (endswith(fname_out, ".nc")) then 2361 iomode = IO_MODE_ETSF 2362 NCF_CHECK(nctk_open_create(ount, fname_out, comm)) 2363 else 2364 iomode = IO_MODE_FORTRAN 2365 if (open_file(fname_out, msg, newunit=ount, status='new', form='unformatted') /= 0) then 2366 ABI_ERROR(msg) 2367 end if 2368 end if 2369 2370 ! Write the header. 2371 fform_merge = ohscr%fform 2372 abifile = abifile_from_fform(fform_merge) 2373 if (abifile%fform == 0) then 2374 ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge))) 2375 end if 2376 varname = abifile%varname 2377 2378 call hscr_io(ohscr,fform_merge,rdwr2,ount,comm,master,iomode) 2379 2380 npwe4mI = ohscr%npwe*ohscr%nI; npwe4mJ = ohscr%npwe*ohscr%nJ 2381 nomega4m = ohscr%nomega 2382 2383 ABI_MALLOC_OR_DIE(epsm1, (npwe4mI,npwe4mJ,nomega4m), ierr) 2384 2385 do iqibz=1,ohscr%nqibz 2386 ! allocate temporary array 2387 npwe4mI = ihscr%npwe * ihscr%nI 2388 npwe4mJ = ihscr%npwe * ihscr%nJ 2389 nomega4m = ihscr%nomega 2390 ABI_MALLOC_OR_DIE(epsm1_temp,(npwe4mI,npwe4mJ,nomega4m), ierr) 2391 2392 ! read full screening matrix for this q-point 2393 call read_screening(varname,inpath,npwe4mI,1,nomega4m,epsm1_temp,iomode,comm,iqiA=iqibz) 2394 2395 ! Copy relevant frequencies 2396 do ifrq=1,nfreq_tot 2397 epsm1(:,:,ifrq) = epsm1_temp(:,:,freq_indx(ifrq)) 2398 end do 2399 2400 ABI_FREE(epsm1_temp) 2401 2402 npwe4mI = ohscr%npwe*ohscr%nI; nomega4m = ohscr%nomega 2403 call write_screening(varname,ount,iomode,npwe4mI,nomega4m,iqibz,epsm1) 2404 end do ! iqibz 2405 2406 ABI_FREE(epsm1) 2407 2408 if (iomode == IO_MODE_FORTRAN) then 2409 close(ount) 2410 else 2411 NCF_CHECK(nf90_close(ount)) 2412 end if 2413 2414 write(msg,'(3a)')ch10,' ==== Frequencies have been removed successfully === ',ch10 2415 call wrtout(std_out, msg) 2416 2417 end subroutine ioscr_wremove
m_io_screening/ncname_from_id [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
ncname_from_id
FUNCTION
Return the name of the netcdf variable (chi0, espilon, em1...) from the id.
SOURCE
279 character(len=nctk_slen) function ncname_from_id(id) result(varname) 280 281 integer,intent(in) :: id 282 283 varname = "None" 284 if (id == 1) varname = chi0_ncname 285 if (id == 3) varname = e_ncname 286 if (id == 4) varname = em1_ncname 287 ABI_CHECK(varname /= "None", "Wrong id") 288 289 end function ncname_from_id
m_io_screening/read_screening [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
read_screening
FUNCTION
Read either a screening (\tilde epsilon^{-1}) file in the SCR format or the irreducible polarizability (chi0) in the SUSC format.
INPUTS
varname=Name of the array to read. Used for ETSF-IO files. iomode=Integer flag defining the format of the output file. Available options: IO_MODE_FORTRAN--> Plain Fortran file IO_MODE_ETSF--> ETSF format iqiA[optional]=Used if only a particular q-point is required. In this case iqiA define the index of the required q-point in the array qibz(3,Hscr%nqibz) nqibzA=number of asked q-points (used to dimension the output arrays). Equal to Hscr%nqibz if the full matrix is required comm=MPI communicator. npweA=number of asked planewaves nomegaA=number of asked frequencies
OUTPUT
epsm1(npweA,npweA,nomegaA,nqibzA) = \tilde\epsilon^{-1}(Ng,Ng,Nw,Nq)
NOTES
* If the epsilon matrix read is bigger than npweA x npweA, it will be truncated; if it is smaller, an error will occur * If the number of frequencies asked for is smaller than that reported in the file, the matrix will be truncated. If nomegaA > Hscr%nomega an error will occur
SOURCE
1400 subroutine read_screening(varname,fname,npweA,nqibzA,nomegaA,epsm1,iomode,comm, & 1401 & iqiA) ! Optional 1402 1403 !Arguments ------------------------------------ 1404 !scalars 1405 integer,intent(in) :: iomode,nomegaA,npweA,nqibzA,comm 1406 integer,optional,intent(in) :: iqiA 1407 character(len=*),intent(in) :: varname,fname 1408 !arrays 1409 complex(gwpc),target,intent(inout) :: epsm1(npweA,npweA,nomegaA,nqibzA) 1410 1411 !Local variables------------------------------- 1412 !scalars 1413 integer,parameter :: master = 0 1414 integer :: ipwe,fform,iomega,iqibz,unt,rdwr,my_rank,nprocs,my_iomode 1415 integer :: varid,ncerr 1416 #ifdef HAVE_MPI_IO 1417 integer :: test_fform,mpi_err,ierr,sc_mode 1418 integer :: bsize_frm,mpi_type_frm 1419 integer :: mpi_fh,buf_dim !,mat_ggw,mat_ggwq 1420 integer(XMPI_OFFSET_KIND) :: offset,displ_wq !,my_offpad 1421 !complex(dpc) :: ctmp 1422 #endif 1423 character(len=500) :: msg,errmsg 1424 logical :: read_qslice 1425 type(hscr_t) :: Hscr 1426 !arrays 1427 #ifdef HAVE_MPI_IO 1428 integer(MPI_OFFSET_KIND),allocatable :: offset_wq(:,:) 1429 #endif 1430 complex(dpc),allocatable :: bufdc2d(:,:),bufdc3d(:,:,:) 1431 ! pointers passed to netcdf4 routines (complex datatypes are not supported). 1432 #ifdef HAVE_GW_DPC 1433 real(dp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:) 1434 #else 1435 real(sp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:) 1436 #endif 1437 integer :: spins(2),s1,s2 1438 1439 ! ************************************************************************* 1440 1441 DBG_ENTER("COLL") 1442 1443 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1444 my_iomode = iomode 1445 if (endswith(fname, ".nc")) my_iomode = IO_MODE_ETSF 1446 !my_iomode = IO_MODE_MPI 1447 !if (my_iomode == IO_MODE_MPI) my_iomode = IO_MODE_FORTRAN 1448 1449 rdwr=1 1450 select case (my_iomode) 1451 case (IO_MODE_MPI) 1452 #ifdef HAVE_MPI_IO 1453 bsize_frm = xmpio_bsize_frm ! bsize_frm= Byte length of the Fortran record marker. 1454 mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker. 1455 sc_mode = xmpio_collective 1456 1457 ! Master reads the header via Fortran IO then bcast the data. 1458 call hscr_from_file(hscr, fname, fform, comm) 1459 1460 ! Open the file with MPI-IO 1461 call MPI_FILE_OPEN(comm, fname, MPI_MODE_RDONLY, xmpio_info ,mpi_fh, mpi_err) 1462 ABI_CHECK_MPI(mpi_err, sjoin("MPI_FILE_OPEN:", fname)) 1463 1464 ! Retrieve the offset of the section immediately below the header. 1465 call hscr_mpio_skip(mpi_fh,test_fform,offset) 1466 ABI_CHECK(test_fform == fform, "mismatch in fform!") 1467 1468 ! Offsets of the Fortran markers corresponding to the (w,q) slices. 1469 ABI_MALLOC(offset_wq,(HScr%nomega,HScr%nqibz)) 1470 displ_wq = offset 1471 do iqibz=1,Hscr%nqibz 1472 do iomega=1,Hscr%nomega 1473 ABI_CHECK(displ_wq > 0, "displ_wq < 0, your SCR|SUSC file is too big for MPI-IO!") 1474 offset_wq(iomega,iqibz) = displ_wq 1475 displ_wq = displ_wq + Hscr%npwe**2 * xmpi_bsize_dpc + Hscr%npwe * 2 * bsize_frm 1476 end do 1477 end do 1478 #else 1479 ABI_ERROR("MPI-IO support not enabled at configure-time") 1480 #endif 1481 1482 case (IO_MODE_FORTRAN) 1483 ! Plain Fortran IO, all nodes read. 1484 if (open_file(fname,msg,newunit=unt,form="unformatted",status="old",action="read") /= 0) then 1485 ABI_ERROR(msg) 1486 end if 1487 call hscr_io(hscr,fform,rdwr,unt,comm,master,my_iomode) 1488 1489 case (IO_MODE_ETSF) 1490 NCF_CHECK(nctk_open_read(unt, fname, xmpi_comm_self)) 1491 call hscr_io(hscr,fform,rdwr,unt,comm,master,my_iomode) 1492 1493 case default 1494 ABI_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode))) 1495 end select 1496 1497 ! Slice or full array? 1498 read_qslice = .False. 1499 if (PRESENT(iqiA)) then 1500 read_qslice = .True. 1501 !call wrtout(std_out, sjoin('. Reading q-slice for iq = ',itoa(iqiA),' from: ', fname)) 1502 if (iqiA <= 0 .or. iqiA > Hscr%nqibz) then 1503 ABI_BUG('iqiA out of range') 1504 end if 1505 end if 1506 1507 ! Do some check 1508 if (Hscr%npwe>npweA) then 1509 write(msg,'(a,i0,2a,i0)')& 1510 'Total number of G-vectors reported on file = ',Hscr%npwe,ch10,& 1511 'Reading a smaller matrix of dimension = ',npweA 1512 ABI_COMMENT(msg) 1513 end if 1514 1515 if (npweA>Hscr%npwe) then 1516 write(msg,'(2(a,i0))')' Dimension of matrix = ',Hscr%npwe," requiring a too big matrix = ",npweA 1517 ABI_ERROR(msg) 1518 end if 1519 1520 ABI_CHECK(nqibzA <= Hscr%nqibz, 'Requiring too many q-points') 1521 ABI_CHECK(nomegaA <= Hscr%nomega,'Requiring too many frequencies') 1522 1523 select case (my_iomode) 1524 case (IO_MODE_MPI) 1525 #ifdef HAVE_MPI_IO 1526 if (read_qslice) then 1527 call wrtout(std_out, "calling mpiotk to read_qslice") 1528 buf_dim = (npweA)**2 * nomegaA 1529 offset = offset_wq(1,iqiA) 1530 sc_mode = xmpio_collective 1531 1532 #ifdef HAVE_GW_DPC 1533 ! Read in-place. 1534 call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],& 1535 buf_dim,epsm1,xmpio_chunk_bsize,sc_mode,comm,ierr) 1536 ABI_CHECK(ierr==0,"Fortran matrix too big") 1537 #else 1538 ! Have to allocate workspace for dpc data. 1539 ! FIXME: Change the file format of the SCR and SUC file so that 1540 ! they are written in single precision if not HAVE_GW_DPC 1541 ABI_MALLOC_OR_DIE(bufdc3d,(npweA,npweA,nomegaA), ierr) 1542 1543 call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],& 1544 buf_dim,bufdc3d,xmpio_chunk_bsize,sc_mode,comm,ierr) 1545 ABI_CHECK(ierr==0,"Fortran matrix too big") 1546 1547 epsm1(:,:,:,1) = bufdc3d 1548 ABI_FREE(bufdc3d) 1549 #endif 1550 1551 else 1552 ! Full matrix (G,G',w,q) is needed. 1553 call wrtout(std_out, "calling mpiotk: Full matrix (G,G',w,q) is needed.") 1554 1555 #ifdef HAVE_GW_DPC 1556 ! Can read all data at once. 1557 buf_dim = (npweA)**2 * nomegaA * HScr%nqibz 1558 offset = offset_wq(1,1) 1559 sc_mode = xmpio_collective 1560 1561 call mpiotk_read_fsuba_dpc4D(mpi_fh,offset,& 1562 [HScr%npwe,HScr%npwe,HScr%nomega,HScr%nqibz], [npweA,npweA,nomegaA,HScr%nqibz], [1,1,1,1],& 1563 buf_dim,epsm1,xmpio_chunk_bsize,sc_mode,comm,ierr) 1564 ABI_CHECK(ierr==0,"Fortran record too big") 1565 #else 1566 ! Have to allocate workspace for dpc data. 1567 ABI_MALLOC_OR_DIE(bufdc3d,(npweA,npweA,nomegaA), ierr) 1568 sc_mode = xmpio_collective 1569 1570 do iqibz=1,Hscr%nqibz 1571 offset = offset_wq(1,iqibz) 1572 buf_dim = (2*npweA)**2 * nomegaA 1573 1574 call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, & 1575 [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],& 1576 buf_dim,bufdc3d,xmpio_chunk_bsize,sc_mode,comm,ierr) 1577 ABI_CHECK(ierr==0,"Fortran matrix too big") 1578 1579 epsm1(:,:,:,iqibz) = bufdc3d 1580 end do 1581 1582 ABI_FREE(bufdc3d) 1583 #endif 1584 end if 1585 1586 call MPI_FILE_CLOSE(mpi_fh,mpi_err) 1587 ABI_FREE(offset_wq) 1588 #endif 1589 1590 case (IO_MODE_FORTRAN) 1591 ! Read epsilon^-1 with Fortran IO 1592 ! Allocate a single column to save memory. 1593 ! TODO re-merge the two cases. 1594 ABI_MALLOC(bufdc2d,(Hscr%npwe,1)) 1595 1596 ! Two coding for different case just to keep it readable. 1597 select case (read_qslice) 1598 case (.True.) 1599 ! Read only a slice of the full array (useful if the entire array is huge). 1600 !if (dim_wings==1) STOP 'not implemented' 1601 !TODO this has to be done in a cleaner way. 1602 qread_loop: & 1603 & do iqibz=1,Hscr%nqibz 1604 if (iqibz==iqiA) then 1605 do iomega=1,nomegaA 1606 do ipwe=1,Hscr%npwe 1607 read(unt, err=10, iomsg=errmsg) bufdc2d(1:Hscr%npwe,1) 1608 if (ipwe<=npweA) epsm1(1:npweA,ipwe,iomega,1)=bufdc2d(1:npweA,1) 1609 end do 1610 end do 1611 EXIT qread_loop ! Got data. Do not need to read file till the end. 1612 else 1613 ! Skip other q-points i.e bufdc2d(1:Hscr%npwe,1:Hscr%npwe) 1614 do iomega=1,Hscr%nomega 1615 do ipwe=1,Hscr%npwe 1616 read(unt, err=10, iomsg=errmsg) 1617 end do 1618 end do 1619 end if ! iqibz==iqiA 1620 end do qread_loop ! iqibz 1621 1622 case (.False.) 1623 ! Read the entire array. 1624 do iqibz=1,Hscr%nqibz 1625 do iomega=1,nomegaA 1626 do ipwe=1,Hscr%npwe 1627 read(unt, err=10, iomsg=errmsg) bufdc2d(1:Hscr%npwe,1) 1628 if (ipwe<=npweA) epsm1(1:npweA,ipwe,iomega,iqibz)=bufdc2d(1:npweA,1) 1629 end do 1630 end do 1631 ! Skip other frequencies 1632 do iomega=nomegaA+1,Hscr%nomega 1633 do ipwe=1,Hscr%npwe 1634 read(unt, err=10, iomsg=errmsg) 1635 end do 1636 end do 1637 end do !iqibz 1638 end select 1639 1640 close(unt) 1641 1642 case (IO_MODE_ETSF) 1643 ! netcdf does not support complex datatypes. Here I use some C-magic to associate the memory 1644 ! to a Fortran real pointer with the correct type and shape. Note that the data on file is always in double precision. 1645 ! nf90_get_var will automatically convert from double to single if the GW code is in single precision mode. 1646 ! This is the reason why I'm using CPP option in the declaration of real_epsm1. 1647 1648 ! FIXME: Need to know the type to read 1649 !write(std_out,*)"in read_screening" 1650 varid = nctk_idname(unt, varname) 1651 1652 ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt] 1653 call c_f_pointer(c_loc(epsm1(1,1,1,1)), real_epsm1, [2,npweA,npweA,1,1,nomegaA,nqibzA]) 1654 spins = 1; s1 = spins(1); s2 = spins(2) 1655 if (read_qslice) then 1656 ncerr = nf90_get_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,iqia], count=[2,npweA,npweA,1,1,nomegaA,1]) 1657 else 1658 ncerr = nf90_get_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,1], count=[2,npweA,npweA,1,1,nomegaA,nqibzA]) 1659 !do iqibz=1,nqibzA 1660 ! !write(*,*)"epsm1: in read ",iqibz,epsm1(1:3,1,1,iqibz) 1661 !end do 1662 end if 1663 NCF_CHECK_MSG(ncerr, sjoin("getting var:", varname)) 1664 NCF_CHECK(nf90_close(unt)) 1665 !write(std_out,*)"read_screening done" 1666 1667 case default 1668 ABI_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode))) 1669 end select 1670 1671 ! Free memory 1672 ABI_SFREE(bufdc2d) 1673 ABI_SFREE(bufdc3d) 1674 1675 call Hscr%free() 1676 1677 DBG_EXIT("COLL") 1678 1679 return 1680 1681 ! Handle Fortran IO error. 1682 10 continue 1683 ABI_ERROR(errmsg) 1684 1685 end subroutine read_screening
m_io_screening/write_screening [ Functions ]
[ Top ] [ m_io_screening ] [ Functions ]
NAME
write_screening
FUNCTION
For a single q-point, write either \tilde epsilon^{-1} on the _SCR file or chi0 on the _SUSC file. The file is supposed to have been open in the calling routine.
INPUTS
varname=The name of the array to write (used if etsf-io format). unt=The unit number of the file to be written (supposed to be already open) iomode=Integer flag defining the format of the output file. Available options: IO_MODE_FORTRAN--> Plain Fortran file IO_MODE_ETSF--> ETSF format npwe=Number of plane waves in epsm1. nomega=Number of frequencies iqibz=Index of the q-points in the IBZ. epsm1(npwe,npwe,nomega)=The matrix to be written, for different frequencies, and a single q-point.
NOTES
On some architecture, the code crashes when trying to write or read a record containing the entire (G1,G2) matrix thus we use smaller records containing the columns of the two-point function.
OUTPUT
(only writing on file)
SOURCE
1300 subroutine write_screening(varname, unt, iomode, npwe, nomega, iqibz, epsm1) 1301 1302 !Arguments ------------------------------------ 1303 !scalars 1304 character(len=*),intent(in) :: varname 1305 integer,intent(in) :: nomega,npwe,iqibz,unt,iomode 1306 !arrays 1307 complex(gwpc),target,intent(in) :: epsm1(npwe,npwe,nomega) 1308 1309 !Local variables------------------------------- 1310 !scalars 1311 integer :: ipwe,iomega,spins(2),s1,s2 1312 character(len=500) :: errmsg 1313 !arrays 1314 complex(dpc),allocatable :: epsm1d(:,:) 1315 integer :: varid,ncerr 1316 #ifdef HAVE_GW_DPC 1317 real(dp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:) 1318 #else 1319 real(sp), ABI_CONTIGUOUS pointer :: real_epsm1(:,:,:,:,:,:,:) 1320 #endif 1321 1322 ! ************************************************************************* 1323 1324 DBG_ENTER("COLL") 1325 1326 select case (iomode) 1327 case (IO_MODE_FORTRAN, IO_MODE_MPI) 1328 ! Write a record for each omega, Always use double precision. 1329 ABI_MALLOC(epsm1d,(npwe,1)) 1330 1331 do iomega=1,nomega 1332 do ipwe=1,npwe 1333 epsm1d(:,1) = epsm1(:,ipwe,iomega) !spc ==> dpc 1334 write(unt, err=10, iomsg=errmsg)epsm1d(1:npwe,1) 1335 end do 1336 end do 1337 ABI_FREE(epsm1d) 1338 1339 case (IO_MODE_ETSF) 1340 ! netcdf does not support complex datatypes. Here I use some C-magic to associate the memory 1341 ! to a Fortran real pointer with the correct type and shape. Note that the data on file is always in double precision. 1342 ! but this is ok since: if the type of data differs from the netCDF variable type, type conversion will occur 1343 ! inside nf90_put_var 1344 varid = nctk_idname(unt, varname) 1345 call c_f_pointer(c_loc(epsm1(1,1,1)), real_epsm1, [2, npwe, npwe, 1, 1, nomega, 1]) 1346 ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt] 1347 spins = 1; s1 = spins(1); s2 = spins(2) 1348 ncerr = nf90_put_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,iqibz], count=[2,npwe,npwe,1,1,nomega,1]) 1349 NCF_CHECK_MSG(ncerr, sjoin("putting var:", varname)) 1350 1351 case default 1352 ABI_ERROR(sjoin("Wrong iomode:", iomode2str(iomode))) 1353 end select 1354 1355 DBG_EXIT("COLL") 1356 1357 return 1358 1359 ! Handle IO error 1360 10 continue 1361 ABI_ERROR(errmsg) 1362 1363 end subroutine write_screening