Line data Source code
1 : !!****m* ABINIT/m_io_screening
2 : !! NAME
3 : !! m_io_screening
4 : !!
5 : !! FUNCTION
6 : !! This module contains the definition of the header of the
7 : !! _SCR and _SUSC file as well as methods used to read/write/echo.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
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, intrinsic :: iso_c_binding
26 : use defs_basis
27 : use m_abicore
28 : USE_MPI
29 : use m_xmpi
30 : use m_mpiotk
31 : use m_nctk
32 : use m_errors
33 : use m_dtset
34 : use netcdf
35 : use m_hdr
36 : use m_sort
37 : use m_crystal
38 : use m_gwdefs
39 :
40 : use m_time, only : cwtime, cwtime_report !, sec2str, timab
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, kmesh_t
46 : use m_gsphere, only : gsphere_t
47 :
48 : implicit none
49 :
50 : private
51 : !!***
52 :
53 : #if defined HAVE_MPI1
54 : include 'mpif.h'
55 : #endif
56 :
57 : character(len=nctk_slen),public,parameter :: e_ncname = "dielectric_function"
58 : character(len=nctk_slen),public,parameter :: em1_ncname = "inverse_dielectric_function"
59 : character(len=nctk_slen),public,parameter :: chi0_ncname = "polarizability"
60 :
61 : public :: ncname_from_id ! return the name of the netcdf variable from the id
62 :
63 : !!****t* m_io_screening/hscr_t
64 : !! NAME
65 : !! hscr_t
66 : !!
67 : !! FUNCTION
68 : !! The structure defining the header of the SCR/SUSC file.
69 : !! hscr_t contains the most important dimensions associated to the SCR/SUSC matrix,
70 : !! important GW metadata and the Abint header. The SCR/SUS matrices are saved
71 : !! with the same format. There are nqibz blocks, each block contains (npwe,npwe,nomega) matrices.
72 : !! SCR and SUS files mainly differ for what concerns the treatment of the q-->0 limit.
73 : !! The treatment of the non-analytic behaviour is not yet implemented but the main ideas are
74 : !! sketched in the NOTES below.
75 : !!
76 : !! NOTES
77 : !! On the treatment of the q-->0 limit
78 : !!
79 : !! 1) q=Gamma should be the first q-point
80 : !!
81 : !! 2) This point contains an initial section with data used to treat the q-->0 limit, followed
82 : !! by the SCR/SUS matrix evaluated for the small q-point (qlwl). The header should contains
83 : !! enough info so that we can skip this section and use the pre-existing routines to read
84 : !! the matrices.
85 : !!
86 : !! 3) The data stored in the q-->0 section depends on the type of matrix stored in the file:
87 : !!
88 : !! SUS file: we store the tensor and the wings needed to reconstruct the q-dependence around Gamma.
89 : !! This section if followed by the X(G,G') matrix evaluated at qlwl.
90 : !! Note that the content of the SUS file does not depend on a possible cutoff in vcoul.
91 : !!
92 : !! SCR file: two cases must be considered:
93 : !!
94 : !! No cutoff in vcoul:
95 : !! The q-->0 section stores the tensor, the wings as well as the inverse of the body B^{-1}.
96 : !! This info is used to integrate W(q) around q==Gamma.
97 : !!
98 : !! cutoff in vcoul:
99 : !! The content of the q-->0 section depends on dimensionality of the system.
100 : !! In 2D we need info on chi0 as well as A and a. See http://arxiv.org/pdf/1511.00129v1.pdf
101 : !! I think that this kind of calculations are easy to implement if we start from the SUS file.
102 : !! Ok, we have to recompute e-1 at each run but the logic is easier to implement.
103 : !!
104 : !! SOURCE
105 :
106 : type,public :: hscr_t
107 :
108 : integer :: id = -1
109 : ! Matrix identifier: 1 for chi0, 2 for chi, 3 for epsilon, 4 for espilon^{-1}
110 :
111 : integer :: ikxc = 0
112 : ! Kxc kernel used,
113 : ! 0 for None (RPA), >0 for static TDDFT (=ixc), <0 for frequency-dependent TDDFT
114 :
115 : integer :: inclvkb = 2
116 : ! q-->0 treatment, 0 for None, 1-2 for transversal gauge, 3 for longitudinal
117 :
118 : integer :: headform
119 : ! format of the SCR header
120 :
121 : integer :: fform
122 : ! File format
123 :
124 : integer :: gwcalctyp = 0
125 : ! Calculation type (G0W0, G0W, GW ...)
126 :
127 : integer :: nI = 1, nJ = 1
128 : ! Number of spin components (rows,columns) in chi|eps^-1. (1,1) if collinear.
129 : ! The internal representation of the matrix is eps(nI*npwe,nJ*npwe)
130 :
131 : integer :: nqibz = -1
132 : ! Number of q-points in the IBZ.
133 :
134 : integer :: nqlwl = 1
135 : ! Number of points for the treatment of the long wavelength limit.
136 :
137 : integer :: nomega = -1
138 : ! Total number of frequencies.
139 :
140 : integer :: nbnds_used = -1
141 : ! Number of bands used during the screening calculation (only for info)
142 :
143 : integer :: npwe = -1
144 : ! Number of G vectors reported on the file.
145 :
146 : integer :: npwwfn_used = -1
147 : ! Number of G vectors for wavefunctions used during the screening calculation (only for info)
148 :
149 : integer :: spmeth = 0
150 : ! Method used to approximate the delta function in the expression for Im Chi_0
151 :
152 : integer :: test_type
153 : ! 1 for TEST-PARTICLE, 2 for TEST-ELECTRON.
154 :
155 : integer :: tordering = 1
156 : ! 1 for Time-Ordered, 2 for Advanced, 3 for Retarded.
157 :
158 : ! HSCR_NEW
159 : integer :: awtr = 1
160 : ! Input variable (time-reversal symmetry in RPA expression)
161 :
162 : integer :: icutcoul = 0
163 : ! Input variable (Coulomb singularity treatment)
164 :
165 : integer :: gwcomp = 0
166 : ! Input variable (GW compensation energy technique)
167 :
168 : integer :: gwgamma = 0
169 : ! Input variable Vertex correction
170 : ! HSCR_NEW
171 :
172 : real(dp) :: mbpt_sciss = zero
173 : ! Scissor Energy, zero if not used
174 :
175 : real(dp) :: spsmear = zero
176 : ! Smearing of the delta in case of spmeth==2
177 :
178 : real(dp) :: zcut = -one
179 : ! Imaginary shift to avoid the poles along the real axis.
180 :
181 : ! HSCR_NEW
182 : real(dp) :: gwencomp = -one
183 : ! Input variable (GW compensation energy technique)
184 :
185 : character(len=3) :: kind_cdata
186 : ! Flag to signal whether the data is in single or double precision ("spc" or "dp")
187 : ! For the time being, we always write/read in double precision.
188 : ! This flag could be use to reduce the memory requirements if spc:
189 : ! we run calculations in single precision dump the results with the same precision without
190 : ! having to allocate extra memory.
191 : ! HSCR_NEW
192 :
193 : !arrays
194 :
195 : ! HSCR_NEW
196 : real(dp) :: vcutgeo(3) = zero
197 : ! Input variable (defines coulomb cutoff)
198 : ! HSCR_NEW
199 :
200 : character(len=80) :: titles(2)
201 : ! Titles describing the content of the file.
202 :
203 : integer,allocatable :: gvec(:,:)
204 : ! gvec(3,npwe)
205 : ! G vectors in reduced coordinates.
206 :
207 : real(dp),allocatable :: qibz(:,:)
208 : ! qibz(3,nqibz)
209 : ! q-points in the IBZ in reduced coordinates.
210 :
211 : real(dp),allocatable :: qlwl(:,:)
212 : ! qlwl(3,nqlwl)
213 : ! q-points for the long wave-length limit treatment (r.l.u)
214 :
215 : complex(dp),allocatable :: omega(:)
216 : ! omega(nomega)
217 : ! All frequencies calculated both along the real and the imaginary axis.
218 : ! Real frequencies are packed in the first section.
219 :
220 : real(dp),allocatable :: omega_wgs(:)
221 : ! (nomega)
222 : ! Weights for numerical integration. Used for instance for minimax meshes.
223 :
224 : character(len=etsfio_charlen) :: iw_mesh_type="None", rw_mesh_type="None", cw_mesh_type="None"
225 : ! String defining the kind of sampling for imaginary (iw), real (rw) and complex (cw) frequencies
226 :
227 : type(hdr_type) :: hdr
228 : ! The abinit header.
229 :
230 : contains
231 : procedure :: from_file => hscr_from_file ! Read the header from file.
232 : procedure :: print => hscr_print ! Print the SCR-related part of the header.
233 : procedure :: bcast => hscr_bcast ! Broadcast the header.
234 : procedure :: free => hscr_free ! Free the header.
235 : procedure :: copy => hscr_copy ! Copy the SCR|SUSC header.
236 : procedure :: io => hscr_io ! I/O of the header (read/write/echo).
237 : end type hscr_t
238 : !!***
239 :
240 : !integer,private,parameter :: HSCR_KNOWN_HEADFORMS(1) = [80]
241 : integer,private,parameter :: HSCR_KNOWN_HEADFORMS(2) = [80, 81]
242 : ! The list of headforms used for SCR/SUSC so far.
243 :
244 : integer,private,parameter :: size_hscr_known_headforms = size(HSCR_KNOWN_HEADFORMS) ! Need this for Flang
245 : integer,public,parameter :: HSCR_LATEST_HEADFORM = HSCR_KNOWN_HEADFORMS(size_hscr_known_headforms)
246 : ! The latest headform used when writing.
247 :
248 : public :: hscr_new ! Create header.
249 : public :: hscr_merge ! Merge two or more headers.
250 : public :: write_screening ! Write a q-slice of the matrix in G-space.
251 : public :: read_screening ! Read the content of the (SCR|SUSC) file placed after the header.
252 : public :: get_hscr_qmesh_gsph
253 :
254 : ! =====================
255 : ! Tools used in mrgscr
256 : ! =====================
257 : public :: ioscr_qmerge ! Produce new file by merging the q-points stored in other files.
258 : public :: ioscr_qrecover ! Recover q-points from a corrupted file produced e.g. from an interrupted run
259 : public :: ioscr_wmerge ! Produce new file by merging the frequencies stored in other files.
260 : public :: ioscr_wremove ! Produce new file by removing selected frequencies in the initial file.
261 :
262 : CONTAINS !================================================================================================
263 : !!***
264 :
265 : !!****f* m_io_screening/ncname_from_id
266 : !! NAME
267 : !! ncname_from_id
268 : !!
269 : !! FUNCTION
270 : !! Return the name of the netcdf variable (chi0, espilon, em1...) from the id.
271 : !!
272 : !! SOURCE
273 :
274 498 : character(len=nctk_slen) function ncname_from_id(id) result(varname)
275 :
276 : integer,intent(in) :: id
277 :
278 498 : varname = "None"
279 498 : if (id == 1) varname = chi0_ncname
280 498 : if (id == 3) varname = e_ncname
281 498 : if (id == 4) varname = em1_ncname
282 498 : ABI_CHECK(varname /= "None", "Wrong id")
283 :
284 498 : end function ncname_from_id
285 : !!***
286 :
287 : !!****f* m_io_screening/hscr_from_file
288 : !! NAME
289 : !! hscr_from_file
290 : !!
291 : !! FUNCTION
292 : !! Read the header of the (SCR/SUS) file
293 : !!
294 : !! INPUTS
295 : !! path=File name
296 : !! comm = MPI communicator.
297 : !!
298 : !! OUTPUT
299 : !! hscr<hscr_t>=The header.
300 : !! fform=Kind of the array in the file (0 signals an error)
301 : !!
302 : !! SOURCE
303 :
304 212 : subroutine hscr_from_file(hscr, path, fform, comm)
305 :
306 : !Arguments ------------------------------------
307 : !scalars
308 : class(hscr_t),intent(out) :: hscr
309 : character(len=*),intent(in) :: path
310 : integer,intent(out) :: fform
311 : integer,intent(in) :: comm
312 :
313 : !Local variables-------------------------------
314 : !scalars
315 : integer,parameter :: rdwr5 = 5, master = 0
316 : integer :: unt,my_rank,ierr
317 : character(len=500) :: msg
318 : ! *************************************************************************
319 :
320 212 : my_rank = xmpi_comm_rank(comm)
321 :
322 : ! Master reads and broadcasts.
323 212 : if (my_rank == master) then
324 184 : if (.not. endswith(path, ".nc")) then
325 : ! Fortran-IO
326 0 : if (open_file(path,msg,newunit=unt,form="unformatted", status="old",action="read") /= 0) then
327 0 : ABI_ERROR(msg)
328 : end if
329 0 : call hscr%io(fform, rdwr5, unt, xmpi_comm_self, master, IO_MODE_FORTRAN)
330 0 : close(unt)
331 : else
332 : ! Netcdf format
333 184 : NCF_CHECK(nctk_open_read(unt, path, xmpi_comm_self))
334 184 : call hscr%io(fform, rdwr5, unt, xmpi_comm_self, master, IO_MODE_ETSF)
335 184 : NCF_CHECK(nf90_close(unt))
336 : end if
337 :
338 184 : ABI_CHECK(fform /= 0, sjoin("hscr_io returned fform == 0 while reading:", path))
339 : end if
340 :
341 : ! Broadcast data.
342 212 : if (xmpi_comm_size(comm) > 1) then
343 42 : call hscr%bcast(master, my_rank, comm)
344 42 : call xmpi_bcast(fform,master,comm,ierr)
345 : end if
346 :
347 212 : end subroutine hscr_from_file
348 : !!***
349 :
350 : !!****f* m_io_screening/hscr_io
351 : !! NAME
352 : !! hscr_io
353 : !!
354 : !! FUNCTION
355 : !! This subroutine deals with the I/O of the hscr_t structured variables (read/write/echo).
356 : !! According to the value of rdwr, it reads the header of a file, writes it, or echo the value
357 : !! of the structured variable to a file. Note that, when reading, different records of hscr_t
358 : !! are allocated here, according to the values of the read variables. Records of hscr_t should be
359 : !! deallocated correctly by a call to hdr_free when hscr_t is not used anymore.
360 : !!
361 : !! INPUTS
362 : !! iomode=Option defining the file format of the external file.
363 : !! comm=MPI communicator.
364 : !! master=rank of the master node in comm, usually 0
365 : !! rdwr= if 1, read the hscr_t structured variable from the header of the file,
366 : !! if 2, write the header to unformatted file
367 : !! if 3, echo part of the header to formatted file (records 1 and 2)
368 : !! if 4, echo the header to formatted file
369 : !! if 5, read the hscr_t without rewinding (unformatted)
370 : !! if 6, write the hscr_t without rewinding (unformatted)
371 : !! unt=unit number of the file (unformatted if rdwr=1, 2, 5 or 6 formatted if rdwr=3,4)
372 : !!
373 : !! OUTPUT
374 : !! (see side effects)
375 : !!
376 : !! SIDE EFFECTS
377 : !! The following variables are both input or output :
378 : !! fform=kind of the array in the file
379 : !! if rdwr=1,5 : will be output ; if the reading fail, return fform=0
380 : !! if rdwr=2,3,4,6 : should be input, will be written or echo to file
381 : !! hscr_t <type(hscr_t)>=the header structured variable
382 : !! if rdwr=1,5 : will be output
383 : !! if rdwr=2,3,4,6 : should be input, will be written or echo to file
384 : !!
385 : !! NOTES
386 : !! In all cases, the file is supposed to be open already
387 : !! When reading (rdwr=1) or writing (rdwr=2), rewind the file
388 : !! When echoing (rdwr=3) does not rewind the file.
389 : !! When reading (rdwr=5) or writing (rdwr=6), DOES NOT rewind the file
390 : !!
391 : !! In writing mode, the routine is supposed to called by the master node.
392 : !! no check is done, it is up to the developer.
393 : !!
394 : !! SOURCE
395 :
396 901 : subroutine hscr_io(hscr, fform, rdwr, unt, comm, master, iomode)
397 :
398 : !Arguments ------------------------------------
399 : !scalars
400 : class(hscr_t),target,intent(inout) :: hscr
401 : integer,intent(inout) :: fform
402 : integer,intent(in) :: rdwr,unt,iomode,comm,master
403 :
404 : !Local variables-------------------------------
405 : !scalars
406 : integer :: my_rank,nprocs,ncerr,ncid,varid,ierr !ii
407 : character(len=500) :: errmsg
408 : character(len=nctk_slen) :: varname !,head_shape,wing_shape
409 : !arrays
410 901 : real(dp),allocatable :: real_omega(:,:)
411 901 : real(dp), contiguous, pointer :: r2vals(:,:) !,rvals3(:,:,:)
412 : ! *************************************************************************
413 :
414 : DBG_ENTER("COLL")
415 : !@hscr_t
416 : ABI_UNUSED(master) ! FIXME
417 :
418 : ! Initialize MPI info for comm ===
419 901 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
420 :
421 901 : if (rdwr==1 .or. rdwr==5) then
422 :
423 : if (.True.) then
424 : ! TODO: only master should read but then I have to skip the header.
425 : !if (my_rank == master) then
426 : ! Read the abinit header, rewinding of the file (if any) is done here.
427 764 : if (iomode==IO_MODE_FORTRAN) then
428 0 : call hscr%hdr%fort_read(unt, fform, rewind=(rdwr==1))
429 764 : else if (iomode==IO_MODE_ETSF) then
430 764 : call hscr%hdr%ncread(unt, fform)
431 : end if
432 :
433 : ! Reset the variables absent in old versions.
434 764 : Hscr%fform=fform
435 :
436 764 : if (iomode==IO_MODE_FORTRAN .or. iomode==IO_MODE_MPI) then
437 0 : select case (fform)
438 : case (1003, 1004)
439 : ! File format for epsilon^-1, espilon, chi0
440 0 : read(unt, err=10, iomsg=errmsg)hscr%titles
441 : read(unt, err=10, iomsg=errmsg)&
442 0 : hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp,&
443 0 : hscr%nI, hscr%nJ, hscr%nqibz, hscr%nqlwl, hscr%nomega, hscr%nbnds_used,&
444 0 : hscr%npwe, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering,&
445 0 : hscr%awtr, hscr%icutcoul, hscr%gwgamma, hscr%vcutgeo(1:3)
446 :
447 : ! Read real scalars
448 : read(unt, err=10, iomsg=errmsg)&
449 0 : hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp, hscr%kind_cdata
450 :
451 : ! Allocate arrays and read them
452 0 : call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
453 :
454 0 : read(unt, err=10, iomsg=errmsg)hscr%gvec(:,:)
455 0 : read(unt, err=10, iomsg=errmsg)hscr%qibz(:,:)
456 0 : read(unt, err=10, iomsg=errmsg)hscr%omega(:)
457 :
458 : ! Read data for q-->0 limit.
459 0 : if (hscr%nqlwl>0) then
460 0 : read(unt, err=10, iomsg=errmsg)hscr%qlwl(:,:)
461 : end if
462 :
463 0 : call set_defaults_from_gwcalctyp()
464 0 : if (hscr%headform > 80) then
465 0 : read(unt, err=10, iomsg=errmsg)hscr%rw_mesh_type
466 0 : read(unt, err=10, iomsg=errmsg)hscr%iw_mesh_type
467 0 : read(unt, err=10, iomsg=errmsg)hscr%cw_mesh_type
468 0 : read(unt, err=10, iomsg=errmsg)hscr%omega_wgs
469 : end if
470 :
471 : case default
472 0 : ABI_BUG(sjoin('Wrong fform read:', itoa(fform)))
473 : end select
474 :
475 764 : else if (iomode == IO_MODE_ETSF) then
476 764 : ncid = unt
477 :
478 764 : select case (fform)
479 : case (1003, 1004)
480 : ! Get dimensions and allocate arrays.
481 764 : NCF_CHECK(nctk_get_dim(ncid, "number_of_coefficients_dielectric_function", hscr%npwe))
482 764 : NCF_CHECK(nctk_get_dim(ncid, "number_of_qpoints_dielectric_function", hscr%nqibz))
483 764 : NCF_CHECK(nctk_get_dim(ncid, "number_of_frequencies_dielectric_function", hscr%nomega))
484 764 : NCF_CHECK(nctk_get_dim(ncid, "number_of_qpoints_gamma_limit", hscr%nqlwl))
485 764 : NCF_CHECK(nctk_get_dim(ncid, "nI", hscr%ni))
486 764 : NCF_CHECK(nctk_get_dim(ncid, "nJ", hscr%nj))
487 764 : call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
488 :
489 764 : varid = nctk_idname(ncid, 'reduced_coordinates_plane_waves_dielectric_function')
490 764 : NCF_CHECK(nf90_get_var(ncid, varid, hscr%gvec, start=[1,1,1]))
491 764 : NCF_CHECK(nf90_get_var(ncid, vid('qpoints_dielectric_function'), hscr%qibz))
492 :
493 2292 : ABI_MALLOC(real_omega, (2, hscr%nomega))
494 764 : NCF_CHECK(nf90_get_var(ncid, vid('frequencies_dielectric_function'), real_omega))
495 5283 : hscr%omega = dcmplx(real_omega(1,:), real_omega(2,:))
496 764 : ABI_FREE(real_omega)
497 :
498 : ! Read extra data added in new header.
499 764 : NCF_CHECK(nf90_get_var(ncid, vid("vcutgeo"), hscr%vcutgeo))
500 764 : NCF_CHECK(nf90_get_var(ncid, vid("id"), hscr%id))
501 764 : NCF_CHECK(nf90_get_var(ncid, vid("ikxc"), hscr%ikxc))
502 764 : NCF_CHECK(nf90_get_var(ncid, vid("inclvkb"), hscr%inclvkb))
503 764 : NCF_CHECK(nf90_get_var(ncid, vid("headform"), hscr%headform))
504 764 : NCF_CHECK(nf90_get_var(ncid, vid("fform"), hscr%fform))
505 764 : NCF_CHECK(nf90_get_var(ncid, vid("gwcalctyp"), hscr%gwcalctyp))
506 764 : NCF_CHECK(nf90_get_var(ncid, vid("nbands_used"), hscr%nbnds_used))
507 764 : NCF_CHECK(nf90_get_var(ncid, vid("npwwfn_used"), hscr%npwwfn_used))
508 764 : NCF_CHECK(nf90_get_var(ncid, vid("spmeth"), hscr%spmeth))
509 764 : NCF_CHECK(nf90_get_var(ncid, vid("test_type"), hscr%test_type))
510 764 : NCF_CHECK(nf90_get_var(ncid, vid("tordering"), hscr%tordering))
511 764 : NCF_CHECK(nf90_get_var(ncid, vid("awtr"), hscr%awtr))
512 764 : NCF_CHECK(nf90_get_var(ncid, vid("gw_icutcoul"), hscr%icutcoul))
513 764 : NCF_CHECK(nf90_get_var(ncid, vid("gwcomp"), hscr%gwcomp))
514 764 : NCF_CHECK(nf90_get_var(ncid, vid("gwgamma"), hscr%gwgamma))
515 764 : NCF_CHECK(nf90_get_var(ncid, vid("mbpt_sciss"), hscr%mbpt_sciss))
516 764 : NCF_CHECK(nf90_get_var(ncid, vid("spsmear"), hscr%spsmear))
517 764 : NCF_CHECK(nf90_get_var(ncid, vid("zcut"), hscr%zcut))
518 764 : NCF_CHECK(nf90_get_var(ncid, vid("gwencomp"), hscr%gwencomp))
519 764 : NCF_CHECK(nf90_get_var(ncid, vid("kind_cdata"), hscr%kind_cdata))
520 764 : call replace_ch0(hscr%kind_cdata)
521 :
522 764 : NCF_CHECK(nf90_get_var(ncid, vid("titles"), hscr%titles))
523 2292 : call replace_ch0(hscr%titles(:))
524 :
525 : ! TODO Read it
526 764 : if (hscr%nqlwl /= 0) then
527 764 : NCF_CHECK(nf90_get_var(ncid, vid("qpoints_gamma_limit"), hscr%qlwl))
528 : end if
529 :
530 764 : call set_defaults_from_gwcalctyp()
531 764 : if (hscr%headform > 80) then
532 764 : NCF_CHECK(nf90_get_var(ncid, vid("rw_mesh_type"), hscr%rw_mesh_type))
533 764 : call replace_ch0(hscr%rw_mesh_type)
534 764 : NCF_CHECK(nf90_get_var(ncid, vid("iw_mesh_type"), hscr%iw_mesh_type))
535 764 : call replace_ch0(hscr%iw_mesh_type)
536 764 : NCF_CHECK(nf90_get_var(ncid, vid("cw_mesh_type"), hscr%cw_mesh_type))
537 764 : call replace_ch0(hscr%cw_mesh_type)
538 764 : NCF_CHECK(nf90_get_var(ncid, vid("omega_wgs"), hscr%omega_wgs))
539 : end if
540 :
541 : case default
542 764 : ABI_BUG(sjoin('Unsupported fform read:',itoa(fform)))
543 : end select
544 : else
545 0 : ABI_ERROR(sjoin("Unsupported value of iomode:", iomode2str(iomode)))
546 : end if
547 :
548 : end if ! master
549 :
550 : !call hscr%bcast(master, my_rank, comm)
551 : !call hscr_mpio_skip(mpio_fh,fform,offset)
552 :
553 137 : else if (rdwr == 2 .or. rdwr == 6) then
554 : ! Writing the header of an unformatted file.
555 : ! Always use the latest version.
556 :
557 137 : if (iomode==IO_MODE_FORTRAN .or. iomode==IO_MODE_MPI) then
558 : ! Write the abinit header.
559 0 : call hscr%hdr%fort_write(unt, fform, ierr)
560 0 : ABI_CHECK(ierr == 0, "hdr_fort_write returned ierr != 0")
561 :
562 0 : write(unt, err=10, iomsg=errmsg)hscr%titles
563 :
564 : ! Write integers
565 : write(unt, err=10, iomsg=errmsg)&
566 0 : hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp,&
567 0 : hscr%nI, hscr%nJ, hscr%nqibz, hscr%nqlwl, hscr%nomega, hscr%nbnds_used,&
568 0 : hscr%npwe, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering,&
569 0 : hscr%awtr, hscr%icutcoul, hscr%gwgamma, hscr%vcutgeo(1:3)
570 :
571 : ! Write real scalars
572 : write(unt, err=10, iomsg=errmsg)&
573 0 : hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp,hscr%kind_cdata
574 :
575 : ! Write arrays
576 0 : write(unt, err=10, iomsg=errmsg)hscr%gvec(:,:)
577 0 : write(unt, err=10, iomsg=errmsg)hscr%qibz(:,:)
578 0 : write(unt, err=10, iomsg=errmsg)hscr%omega(:)
579 :
580 : ! Add q-points for heads and wings for q-->0.
581 0 : if (hscr%nqlwl > 0) then
582 0 : write(unt, err=10, iomsg=errmsg)hscr%qlwl(:,:)
583 : end if
584 :
585 0 : if (hscr%headform > 80) then
586 0 : write(unt, err=10, iomsg=errmsg)hscr%rw_mesh_type
587 0 : write(unt, err=10, iomsg=errmsg)hscr%iw_mesh_type
588 0 : write(unt, err=10, iomsg=errmsg)hscr%cw_mesh_type
589 0 : write(unt, err=10, iomsg=errmsg)hscr%omega_wgs
590 : end if
591 :
592 137 : else if (iomode == IO_MODE_ETSF) then
593 137 : ncid = unt
594 : ! Write the abinit header, rewinding of the file (if any) is done here.
595 137 : NCF_CHECK(hscr%hdr%ncwrite(ncid, fform, nc_define=.True.))
596 :
597 : ! Define dimensions
598 : ! Part 2) of etsf-io specifications
599 : ! FIXME: Spin is only used in particular cases, We usually get the trace of W in spin space
600 : ! and I'm not gonna allocate extra memory just to have up up, down down
601 : ! Besides number_of_spins should be replaced by `number_of_spins_dielectric_function`
602 : ! Should add spin_dependent attribute.
603 : ncerr = nctk_def_dims(ncid, [ &
604 : nctkdim_t("complex", 2), nctkdim_t("number_of_reduced_dimensions", 3), &
605 : nctkdim_t("number_of_frequencies_dielectric_function", hscr%nomega), &
606 : nctkdim_t("number_of_qpoints_dielectric_function", hscr%nqibz), &
607 : nctkdim_t("number_of_qpoints_gamma_limit", hscr%nqlwl), &
608 : nctkdim_t("number_of_spins", hscr%hdr%nsppol), &
609 : nctkdim_t("nI", hscr%nI), nctkdim_t("nJ", hscr%nJ), &
610 1370 : nctkdim_t("number_of_coefficients_dielectric_function", hscr%npwe)], defmode=.True.)
611 137 : NCF_CHECK(ncerr)
612 :
613 : ! Part 3) of the specs
614 : ! (note that, in the specs, the Gs depend on the q-point but npwe is a scalar
615 : ! basis_set is added by the abinit header.
616 : ! FIXME: g-vectors are not written properly.
617 : ncerr = nctk_def_arrays(ncid, [&
618 : ! Standard
619 : nctkarr_t('frequencies_dielectric_function', "dp", 'complex, number_of_frequencies_dielectric_function'), &
620 : nctkarr_t('qpoints_dielectric_function', "dp", 'number_of_reduced_dimensions, number_of_qpoints_dielectric_function'),&
621 : nctkarr_t('qpoints_gamma_limit', "dp", 'number_of_reduced_dimensions, number_of_qpoints_gamma_limit'), &
622 : nctkarr_t('reduced_coordinates_plane_waves_dielectric_function', "i", &
623 : 'number_of_reduced_dimensions, number_of_coefficients_dielectric_function, number_of_qpoints_dielectric_function'), &
624 : ! Abinit
625 : nctkarr_t('vcutgeo', "dp", 'number_of_reduced_dimensions'), &
626 : nctkarr_t('kind_cdata', "char", 'character_string_length'), &
627 : nctkarr_t('titles', "char", 'character_string_length, two'), &
628 : nctkarr_t('rw_mesh_type', "char", 'character_string_length'), &
629 : nctkarr_t('iw_mesh_type', "char", 'character_string_length'), &
630 : nctkarr_t('cw_mesh_type', "char", 'character_string_length'), &
631 : nctkarr_t('omega_wgs', "dp", 'number_of_frequencies_dielectric_function') &
632 1644 : ])
633 137 : NCF_CHECK(ncerr)
634 :
635 : ! FIXME problem with q-points, heads and wings?
636 : ! The order in P in the specs is wrong, q should be the last dimension here I use the "correct" version
637 : ! TODO: Remove. Use abifile_t
638 137 : varname = ncname_from_id(hscr%id)
639 : ncerr = nctk_def_arrays(ncid, &
640 : nctkarr_t(varname, "dp", &
641 : &"complex, number_of_coefficients_dielectric_function, number_of_coefficients_dielectric_function,&
642 137 : &number_of_spins, number_of_spins, number_of_frequencies_dielectric_function, number_of_qpoints_dielectric_function"))
643 137 : NCF_CHECK(ncerr)
644 :
645 : !write(std_out,*)"nqlwl",hscr%nqlwl
646 137 : NCF_CHECK(nctk_set_datamode(ncid))
647 411 : call c_f_pointer(c_loc(hscr%omega(1)), r2vals, shape=[2, size(hscr%omega)])
648 137 : NCF_CHECK(nf90_put_var(ncid, vid('frequencies_dielectric_function'), r2vals))
649 137 : NCF_CHECK(nf90_put_var(ncid, vid('qpoints_dielectric_function'), hscr%qibz))
650 137 : NCF_CHECK(nf90_put_var(ncid, vid('reduced_coordinates_plane_waves_dielectric_function'), hscr%gvec))
651 :
652 137 : NCF_CHECK(nf90_put_var(ncid, vid("titles"), hscr%titles))
653 137 : NCF_CHECK(nf90_put_var(ncid, vid("kind_cdata"), hscr%kind_cdata))
654 137 : NCF_CHECK(nf90_put_var(ncid, vid("vcutgeo"), hscr%vcutgeo))
655 :
656 : ncerr = nctk_defnwrite_ivars(ncid, [character(len=nctk_slen) :: &
657 : "id", "ikxc", "inclvkb", "headform", "fform", "gwcalctyp", &
658 : "nbands_used", "npwwfn_used", "spmeth", "test_type", "tordering", "awtr", "gw_icutcoul", &
659 : "gwcomp", "gwgamma" &
660 : ],&
661 : [ hscr%id, hscr%ikxc, hscr%inclvkb, hscr%headform, hscr%fform, hscr%gwcalctyp, &
662 : hscr%nbnds_used, hscr%npwwfn_used, hscr%spmeth, hscr%test_type, hscr%tordering, hscr%awtr, hscr%icutcoul, &
663 : hscr%gwcomp, hscr%gwgamma &
664 4247 : ])
665 137 : NCF_CHECK(ncerr)
666 :
667 : ncerr = nctk_defnwrite_dpvars(ncid, [character(len=nctk_slen) :: &
668 : "mbpt_sciss", "spsmear", "zcut", "gwencomp"], &
669 : [hscr%mbpt_sciss, hscr%spsmear, hscr%zcut, hscr%gwencomp &
670 1233 : ])
671 137 : NCF_CHECK(ncerr)
672 :
673 : ! Add q-points for heads and wings for q-->0.
674 137 : if (hscr%nqlwl > 0) then
675 : ! MG: This part has been commented out as it's not used
676 : ! head_shape = "complex, number_of_spins, number_of_spins, number_of_frequencies_dielectric_function"
677 : ! head_shape = trim(head_shape)//", number_of_qpoints_gamma_limit"
678 :
679 : ! wing_shape = "complex, number_of_coefficients_dielectric_function, number_of_spins, number_of_spins"
680 : ! wing_shape = trim(wing_shape)//", number_of_frequencies_dielectric_function, number_of_qpoints_gamma_limit"
681 :
682 : ! ncerr = nctk_def_arrays(ncid, [&
683 : ! nctkarr_t("dielectric_function_head", "dp", head_shape),&
684 : ! nctkarr_t("dielectric_function_upper_wing", "dp", wing_shape),&
685 : ! nctkarr_t("dielectric_function_lower_wing", "dp", wing_shape)], defmode=.True.)
686 : ! NCF_CHECK(ncerr)
687 :
688 137 : NCF_CHECK(nctk_set_datamode(ncid))
689 137 : NCF_CHECK(nf90_put_var(ncid, vid('qpoints_gamma_limit'), hscr%qlwl))
690 :
691 : !if (hscr%headform > 80) then
692 137 : NCF_CHECK(nf90_put_var(ncid, vid("rw_mesh_type"), hscr%rw_mesh_type))
693 137 : NCF_CHECK(nf90_put_var(ncid, vid("iw_mesh_type"), hscr%iw_mesh_type))
694 137 : NCF_CHECK(nf90_put_var(ncid, vid("cw_mesh_type"), hscr%cw_mesh_type))
695 137 : NCF_CHECK(nf90_put_var(ncid, vid("omega_wgs"), hscr%omega_wgs))
696 : !end if
697 :
698 : end if
699 : else
700 0 : ABI_ERROR(sjoin('Unsupported iomode:',iomode2str(iomode)))
701 : end if
702 :
703 : else
704 0 : ABI_BUG(sjoin("Wrong value for rdwr:", itoa(rdwr)))
705 : end if ! read/write/echo
706 :
707 : DBG_EXIT("COLL")
708 :
709 901 : return
710 :
711 : ! Handle Fortran IO error
712 : 10 continue
713 1802 : ABI_ERROR(errmsg)
714 :
715 : contains
716 23663 : integer function vid(vname)
717 : character(len=*),intent(in) :: vname
718 23663 : vid = nctk_idname(ncid, vname)
719 : end function vid
720 :
721 764 : subroutine set_defaults_from_gwcalctyp()
722 : integer :: mod10
723 : ! Ppovide default values of mesh type on the basis of gw_calctyp if old format is read.
724 764 : mod10 = mod(hscr%gwcalctyp, 10)
725 764 : hscr%iw_mesh_type = "None"; hscr%rw_mesh_type = "None"; hscr%cw_mesh_type = "None"
726 :
727 764 : if (mod10 == SIG_GW_AC) then
728 32 : hscr%iw_mesh_type = "gauss_legendre"
729 732 : else if (any(mod10 == [SIG_QPGW_CD, SIG_QPGW_CD])) then
730 8 : hscr%rw_mesh_type = "linear"
731 8 : hscr%iw_mesh_type = "logarithmic"
732 : end if
733 764 : end subroutine set_defaults_from_gwcalctyp
734 :
735 : end subroutine hscr_io
736 : !!***
737 :
738 : !----------------------------------------------------------------------
739 :
740 : !!****f* m_io_screening/hscr_print
741 : !! NAME
742 : !! hscr_print
743 : !!
744 : !! FUNCTION
745 : !! Prints info on the header of the SCR|SUSC file.
746 : !!
747 : !! SOURCE
748 :
749 332 : subroutine hscr_print(Hscr, units, prtvol, header)
750 :
751 : !Arguments ------------------------------------
752 : !scalars
753 : class(hscr_t),intent(in) :: hscr
754 : integer,intent(in) :: units(:)
755 : integer,intent(in) :: prtvol
756 : character(len=*),intent(in),optional :: header
757 :
758 : !Local variables-------------------------------
759 : !scalars
760 : integer :: iomega,iq_ibz
761 : character(len=500) :: msg
762 : ! *************************************************************************
763 :
764 166 : if (PRESENT(header)) then
765 11 : msg = ' ==== '//TRIM(ADJUSTL(header))//' ==== '
766 11 : call wrtout(units, msg)
767 : end if
768 :
769 166 : write(msg,'(1x,a)')TRIM(hscr%titles(1))
770 166 : call wrtout(units,msg)
771 166 : write(msg,'(1x,a)')TRIM(hscr%titles(2))
772 166 : call wrtout(units,msg)
773 166 : write(msg,'(a,i8)') ' Identifier ',hscr%ID
774 166 : call wrtout(units,msg)
775 166 : write(msg,'(a,i8)') ' Kxc kernel ',hscr%ikxc
776 166 : call wrtout(units,msg)
777 166 : write(msg,'(a,i8)') ' Treatment of q-->0 limit ',hscr%inclvkb
778 166 : call wrtout(units,msg)
779 166 : write(msg,'(a,i8)') '- headform ',hscr%headform
780 166 : call wrtout(units,msg)
781 166 : write(msg,'(a,i8)') ' fform ',hscr%fform
782 166 : call wrtout(units,msg)
783 166 : write(msg,'(a,i8)') ' gwcalctyp ',hscr%gwcalctyp
784 166 : call wrtout(units,msg)
785 166 : write(msg,'(a,2i8)')' Number of components ',hscr%nI,hscr%nJ
786 166 : call wrtout(units,msg)
787 166 : write(msg,'(a,i8)') ' Number of q-points ',hscr%nqibz
788 166 : call wrtout(units,msg)
789 166 : write(msg,'(a,i8)') ' Number of q-directions ',hscr%nqlwl
790 166 : call wrtout(units,msg)
791 166 : write(msg,'(a,i8)') ' Number of frequencies ',hscr%nomega
792 166 : call wrtout(units,msg)
793 166 : write(msg,'(a,i8)') ' Number of bands used ',hscr%nbnds_used
794 166 : call wrtout(units,msg)
795 166 : write(msg,'(a,i8)') ' Dimension of matrix ',hscr%npwe
796 166 : call wrtout(units,msg)
797 166 : write(msg,'(a,i8)') ' Number of planewaves used ',hscr%npwwfn_used
798 166 : call wrtout(units,msg)
799 166 : write(msg,'(a,i8)') ' Spectral method ',hscr%spmeth
800 166 : call wrtout(units,msg)
801 166 : write(msg,'(a,i8)') ' Test_type ',hscr%test_type
802 166 : call wrtout(units,msg)
803 166 : write(msg,'(a,i8)') ' Time-ordering ',hscr%tordering
804 166 : call wrtout(units,msg)
805 166 : write(msg,'(a,es16.6)')' Scissor Energy ',hscr%mbpt_sciss
806 166 : call wrtout(units,msg)
807 166 : write(msg,'(a,es16.6)')' Spectral smearing ',hscr%spsmear
808 166 : call wrtout(units,msg)
809 166 : write(msg,'(a,es16.6)')' Complex Imaginary Shift ',hscr%zcut
810 166 : call wrtout(units,msg)
811 166 : write(msg,'(2a)') ' rw_mesh_type ',trim(hscr%rw_mesh_type)
812 166 : call wrtout(units,msg)
813 166 : write(msg,'(2a)') ' iw_mesh_type ',trim(hscr%iw_mesh_type)
814 166 : call wrtout(units,msg)
815 166 : write(msg,'(2a)') ' cw_mesh_type ',trim(hscr%cw_mesh_type)
816 166 : call wrtout(units,msg)
817 :
818 166 : if (prtvol == 0) then
819 142 : call wrtout(units,' The header contains additional records.')
820 : else
821 24 : write(msg,'(2a)')ch10,' q-points [r.l.u.]:'
822 24 : call wrtout(units,msg)
823 148 : do iq_ibz=1,hscr%nqibz
824 124 : write(msg,'(i5,3f12.6)')iq_ibz,hscr%qibz(:,iq_ibz)
825 148 : call wrtout(units,msg)
826 : end do
827 :
828 24 : write(msg,'(2a)')ch10,' Frequencies used [eV]:'
829 24 : call wrtout(units,msg)
830 248 : do iomega=1,hscr%nomega
831 224 : write(msg,'(i3,2f7.2)')iomega,REAL(hscr%omega(iomega))*Ha_eV,AIMAG(hscr%omega(iomega))*Ha_eV
832 248 : call wrtout(units,msg)
833 : end do
834 : end if
835 :
836 : ! Echo the abinit header.
837 : !if (prtvol>0) call hdr_echo(hscr%hdr,fform,rdwr,unit=unt)
838 :
839 166 : end subroutine hscr_print
840 : !!***
841 :
842 : !----------------------------------------------------------------------
843 :
844 : !!****f* m_io_screening/hscr_new
845 : !! NAME
846 : !! hscr_new
847 : !!
848 : !! FUNCTION
849 : !! Initialize the Hscr datatype and most of its content from the em1params_t data type Ep.
850 : !!
851 : !! INPUTS
852 : !! varname=Name of the netcdf variable (used to get fform and ID).
853 : !! ikxc=Integer flag defining the type of XC kernel (0 if None i.e RPA)
854 : !! test_type=Integer flag defining the type of probing charge (0 for None)
855 : !! tordering=The time-ordering of the Response function.
856 : !! gvec(3,Ep%npwe)=The G-vectors used.
857 : !! Ep<em1params_t>=Parameters defining the calculation of the screening.
858 : !! hdr_abinit<hdr_type>=The abinit header.
859 : !!
860 : !! OUTPUT
861 : !! Hscr<type(hscr_t)>=the header, initialized.
862 : !!
863 : !! SOURCE
864 :
865 256 : type(hscr_t) function hscr_new(varname, dtset, ep, hdr_abinit, ikxc, test_type, tordering, titles, ngvec, gvec) result(hscr)
866 :
867 : !Arguments ------------------------------------
868 : !scalars
869 : integer,intent(in) :: ikxc,test_type,tordering,ngvec
870 : character(len=*),intent(in) :: varname
871 : type(dataset_type),intent(in) :: dtset
872 : type(em1params_t),intent(in) :: Ep
873 : type(hdr_type),intent(in) :: hdr_abinit
874 : !arrays
875 : integer,intent(in) :: gvec(3,ngvec)
876 : character(len=80),intent(in) :: titles(2)
877 :
878 : !Local variables-------------------------------
879 : integer :: id
880 : type(abifile_t) :: abifile
881 : ! *************************************************************************
882 :
883 : !@hscr_t
884 0 : ABI_CHECK(ngvec == Ep%npwe, 'ngvec/=Ep%npwe')
885 :
886 : ! Identifier used to define the type of response function (e^-1, chi0)
887 128 : id = 0
888 128 : if (varname == "polarizability") id = 1
889 128 : if (varname == "inverse_dielectric_function") id = 4
890 8 : ABI_CHECK(id /= 0, sjoin("Invalid varname: ",varname))
891 :
892 : ! Get fform from abifile.
893 128 : abifile = abifile_from_varname(varname)
894 128 : if (abifile%fform == 0) then
895 0 : ABI_ERROR(sjoin("Cannot find any abifile object associated to varname:", varname))
896 : end if
897 :
898 : ! Copy the abinit header.
899 128 : call hdr_abinit%copy(Hscr%Hdr)
900 :
901 : ! Initialize quantities related to the screening file
902 128 : hscr%id =id
903 128 : hscr%ikxc =ikxc
904 128 : hscr%inclvkb =Ep%inclvkb
905 128 : hscr%headform =HSCR_LATEST_HEADFORM
906 128 : hscr%fform =abifile%fform
907 128 : hscr%gwcalctyp =Ep%gwcalctyp
908 128 : hscr%nI =Ep%nI
909 128 : hscr%nJ =Ep%nJ
910 128 : hscr%nqibz =Ep%nqcalc ! nqcalc == nqibz except if we split the calculation with nqptdm
911 128 : hscr%nqlwl =Ep%nqlwl
912 128 : hscr%nomega =Ep%nomega
913 128 : hscr%nbnds_used =Ep%nbnds
914 128 : hscr%npwe =Ep%npwe
915 128 : hscr%npwwfn_used=Ep%npwwfn
916 128 : hscr%spmeth =Ep%spmeth
917 128 : hscr%test_type =test_type
918 128 : hscr%tordering =tordering
919 128 : hscr%mbpt_sciss =Ep%mbpt_sciss
920 128 : hscr%spsmear =Ep%spsmear
921 128 : hscr%zcut =Ep%zcut
922 :
923 384 : hscr%titles(:)=titles(:)
924 :
925 128 : call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
926 29808 : hscr%gvec(:,:) = gvec(:,1:Ep%npwe)
927 3156 : hscr%qibz(:,:) = Ep%qcalc
928 640 : hscr%qlwl(:,:) = Ep%qlwl
929 995 : hscr%omega(:) = Ep%omega
930 :
931 995 : hscr%omega_wgs(:) = Ep%omega_wgs
932 128 : hscr%iw_mesh_type = Ep%iw_mesh_type
933 128 : hscr%rw_mesh_type = Ep%rw_mesh_type
934 128 : hscr%cw_mesh_type = Ep%cw_mesh_type
935 :
936 : ! HSCR_NEW
937 128 : hscr%awtr = dtset%awtr
938 128 : hscr%icutcoul = dtset%gw_icutcoul
939 512 : hscr%vcutgeo = dtset%vcutgeo
940 128 : hscr%gwcomp = dtset%gwcomp
941 128 : hscr%gwgamma = dtset%gwgamma
942 128 : hscr%gwencomp = dtset%gwencomp
943 128 : hscr%kind_cdata = "dpc" ! For the time being, data is always written in double-precision
944 : ! HSCR_NEW
945 :
946 512 : end function hscr_new
947 : !!***
948 :
949 : !----------------------------------------------------------------------
950 :
951 : !!****f* m_io_screening/hscr_bcast
952 : !! NAME
953 : !! hscr_bcast
954 : !!
955 : !! FUNCTION
956 : !! This subroutine transmit the header structured datatype initialized
957 : !! on one processor (or a group of processor), to the other processors.
958 : !! It also allocates the needed part of the header.
959 : !!
960 : !! INPUTS
961 : !! master=ID of the master node.
962 : !! my_rank=ID of the node that receives the data.
963 : !! comm=MPI communicator.
964 : !!
965 : !! OUTPUT
966 : !! (no output)
967 : !!
968 : !! SIDE EFFECTS
969 : !! Hscr<type(hscr_t)>=the SCR header. For the master, it is already
970 : !! initialized entirely, while for the other procs, everything has
971 : !! to be transmitted.
972 : !!
973 : !! SOURCE
974 :
975 1487 : subroutine hscr_bcast(hscr, master, my_rank, comm)
976 :
977 : !Arguments ------------------------------------
978 : class(hscr_t),intent(inout) :: hscr
979 : integer, intent(in) :: master, my_rank, comm
980 :
981 : !Local variables-------------------------------
982 : integer :: ierr
983 : ! *************************************************************************
984 :
985 : DBG_ENTER("COLL")
986 59 : if (xmpi_comm_size(comm) == 1) return ! Nothing to do
987 :
988 : ! integer
989 42 : call xmpi_bcast(hscr%id, master,comm,ierr)
990 42 : call xmpi_bcast(hscr%ikxc, master,comm,ierr)
991 42 : call xmpi_bcast(hscr%inclvkb, master,comm,ierr)
992 42 : call xmpi_bcast(hscr%headform, master,comm,ierr)
993 42 : call xmpi_bcast(hscr%fform, master,comm,ierr)
994 42 : call xmpi_bcast(hscr%gwcalctyp, master,comm,ierr)
995 42 : call xmpi_bcast(hscr%nI, master,comm,ierr)
996 42 : call xmpi_bcast(hscr%nJ, master,comm,ierr)
997 42 : call xmpi_bcast(hscr%nqibz, master,comm,ierr)
998 42 : call xmpi_bcast(hscr%nqlwl, master,comm,ierr)
999 42 : call xmpi_bcast(hscr%nomega, master,comm,ierr)
1000 42 : call xmpi_bcast(hscr%nbnds_used, master,comm,ierr)
1001 42 : call xmpi_bcast(hscr%npwe, master,comm,ierr)
1002 42 : call xmpi_bcast(hscr%npwwfn_used,master,comm,ierr)
1003 42 : call xmpi_bcast(hscr%spmeth, master,comm,ierr)
1004 42 : call xmpi_bcast(hscr%test_type, master,comm,ierr)
1005 42 : call xmpi_bcast(hscr%tordering, master,comm,ierr)
1006 :
1007 : ! Real
1008 42 : call xmpi_bcast(hscr%mbpt_sciss, master,comm,ierr)
1009 42 : call xmpi_bcast(hscr%spsmear, master,comm,ierr)
1010 42 : call xmpi_bcast(hscr%zcut, master,comm,ierr)
1011 :
1012 : ! arrays
1013 42 : call xmpi_bcast(hscr%titles, master,comm,ierr)
1014 :
1015 42 : if (my_rank /= master) call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
1016 :
1017 42 : call xmpi_bcast(hscr%gvec, master,comm,ierr)
1018 42 : call xmpi_bcast(hscr%qibz, master,comm,ierr)
1019 42 : call xmpi_bcast(hscr%qlwl, master,comm,ierr)
1020 42 : call xmpi_bcast(hscr%omega,master,comm,ierr)
1021 :
1022 : ! Communicate the Abinit header.
1023 42 : call hscr%Hdr%bcast(master, my_rank, comm)
1024 :
1025 : ! HSCR_NEW
1026 42 : call xmpi_bcast(hscr%awtr, master, comm, ierr)
1027 42 : call xmpi_bcast(hscr%icutcoul, master, comm, ierr)
1028 42 : call xmpi_bcast(hscr%vcutgeo, master, comm, ierr)
1029 42 : call xmpi_bcast(hscr%gwcomp, master, comm, ierr)
1030 42 : call xmpi_bcast(hscr%gwgamma, master, comm, ierr)
1031 42 : call xmpi_bcast(hscr%gwencomp, master, comm, ierr)
1032 42 : call xmpi_bcast(hscr%kind_cdata, master, comm, ierr)
1033 : ! HSCR_NEW
1034 :
1035 42 : call xmpi_bcast(hscr%omega_wgs, master, comm, ierr)
1036 42 : call xmpi_bcast(hscr%iw_mesh_type, master, comm, ierr)
1037 42 : call xmpi_bcast(hscr%rw_mesh_type, master, comm, ierr)
1038 42 : call xmpi_bcast(hscr%cw_mesh_type, master, comm, ierr)
1039 :
1040 : DBG_EXIT("COLL")
1041 :
1042 : end subroutine hscr_bcast
1043 : !!***
1044 :
1045 : !----------------------------------------------------------------------
1046 :
1047 : !!****f* m_io_screening/hscr_malloc
1048 : !! NAME
1049 : !! hscr_malloc
1050 : !!
1051 : !! FUNCTION
1052 : !! Allocate the components of the header structured datatype except for hscr%hdr
1053 : !!
1054 : !! SOURCE
1055 :
1056 920 : subroutine hscr_malloc(hscr, npwe, nqibz, nomega, nqlwl)
1057 :
1058 : !Arguments ------------------------------------
1059 : !scalars
1060 : class(hscr_t),intent(inout) :: Hscr
1061 : integer,intent(in) :: npwe, nqibz, nomega, nqlwl
1062 : ! *************************************************************************
1063 :
1064 : !@hscr_t
1065 2760 : ABI_MALLOC(hscr%gvec, (3, npwe))
1066 2760 : ABI_MALLOC(hscr%qibz, (3, nqibz))
1067 2760 : ABI_MALLOC(hscr%qlwl, (3, nqlwl))
1068 2760 : ABI_MALLOC(hscr%omega, (nomega))
1069 7478 : ABI_CALLOC(hscr%omega_wgs, (nomega))
1070 :
1071 920 : end subroutine hscr_malloc
1072 : !!***
1073 :
1074 : !----------------------------------------------------------------------
1075 :
1076 : !!****f* m_io_screening/hscr_free
1077 : !! NAME
1078 : !! hscr_free
1079 : !!
1080 : !! FUNCTION
1081 : !! Deallocate the components of the header structured datatype
1082 : !!
1083 : !! SOURCE
1084 :
1085 965 : subroutine hscr_free(hscr)
1086 :
1087 : !Arguments ------------------------------------
1088 : class(hscr_t),intent(inout) :: hscr
1089 : ! *************************************************************************
1090 :
1091 965 : ABI_SFREE(hscr%gvec)
1092 965 : ABI_SFREE(hscr%qibz)
1093 965 : ABI_SFREE(hscr%qlwl)
1094 965 : ABI_SFREE(hscr%omega)
1095 965 : ABI_SFREE(hscr%omega_wgs)
1096 :
1097 965 : call hscr%Hdr%free()
1098 :
1099 965 : end subroutine hscr_free
1100 : !!***
1101 :
1102 : !----------------------------------------------------------------------
1103 :
1104 : !!****f* m_io_screening/hscr_copy
1105 : !! NAME
1106 : !! hscr_copy
1107 : !!
1108 : !! FUNCTION
1109 : !! Deep copy of the header of the _SCR or _SUSC file.
1110 : !!
1111 : !! INPUTS
1112 : !!
1113 : !! SOURCE
1114 :
1115 9 : subroutine hscr_copy(Hscr_in, Hscr_cp)
1116 :
1117 : !Arguments ------------------------------------
1118 : !scalars
1119 : class(hscr_t),intent(in) :: Hscr_in
1120 : class(hscr_t),intent(inout) :: Hscr_cp
1121 : ! *************************************************************************
1122 :
1123 : !@hscr_t
1124 : ! Integer values.
1125 9 : Hscr_cp%id = Hscr_in%id
1126 9 : Hscr_cp%ikxc = Hscr_in%ikxc
1127 9 : Hscr_cp%inclvkb = Hscr_in%inclvkb
1128 9 : Hscr_cp%headform = Hscr_in%headform
1129 9 : Hscr_cp%fform = Hscr_in%fform
1130 9 : Hscr_cp%gwcalctyp = Hscr_in%gwcalctyp
1131 9 : Hscr_cp%nI = Hscr_in%nI
1132 9 : Hscr_cp%nJ = Hscr_in%nJ
1133 9 : Hscr_cp%nqibz = Hscr_in%nqibz
1134 9 : Hscr_cp%nqlwl = Hscr_in%nqlwl
1135 9 : Hscr_cp%nomega = Hscr_in%nomega
1136 9 : Hscr_cp%nbnds_used = Hscr_in%nbnds_used
1137 9 : Hscr_cp%npwe = Hscr_in%npwe
1138 9 : Hscr_cp%npwwfn_used = Hscr_in%npwwfn_used
1139 9 : Hscr_cp%spmeth = Hscr_in%spmeth
1140 9 : Hscr_cp%test_type = Hscr_in%test_type
1141 9 : Hscr_cp%tordering = Hscr_in%tordering
1142 :
1143 : ! Real variables
1144 9 : Hscr_cp%mbpt_sciss = Hscr_in%mbpt_sciss
1145 9 : Hscr_cp%spsmear = Hscr_in%spsmear
1146 9 : Hscr_cp%zcut = Hscr_in%zcut
1147 :
1148 : ! Copy the abinit Header
1149 9 : call Hscr_in%Hdr%copy(Hscr_cp%Hdr)
1150 :
1151 27 : Hscr_cp%titles(:) = Hscr_in%titles(:)
1152 :
1153 : ! Copy allocatable arrays.
1154 9 : call alloc_copy(Hscr_in%gvec , Hscr_cp%gvec)
1155 9 : call alloc_copy(Hscr_in%qibz , Hscr_cp%qibz)
1156 9 : call alloc_copy(Hscr_in%qlwl , Hscr_cp%qlwl)
1157 9 : call alloc_copy(Hscr_in%omega, Hscr_cp%omega)
1158 9 : call alloc_copy(Hscr_in%omega_wgs, Hscr_cp%omega_wgs)
1159 :
1160 : ! HSCR_NEW
1161 9 : hscr_cp%awtr = hscr_in%awtr
1162 9 : hscr_cp%icutcoul = hscr_in%icutcoul
1163 36 : hscr_cp%vcutgeo = hscr_in%vcutgeo
1164 9 : hscr_cp%gwcomp = hscr_in%gwcomp
1165 9 : hscr_cp%gwgamma = hscr_in%gwgamma
1166 9 : hscr_cp%gwencomp = hscr_in%gwencomp
1167 9 : hscr_cp%kind_cdata = hscr_in%kind_cdata
1168 :
1169 9 : hscr_cp%iw_mesh_type = hscr_in%iw_mesh_type
1170 9 : hscr_cp%rw_mesh_type = hscr_in%rw_mesh_type
1171 9 : hscr_cp%cw_mesh_type = hscr_in%cw_mesh_type
1172 : ! HSCR_NEW
1173 :
1174 9 : end subroutine hscr_copy
1175 : !!***
1176 :
1177 : !----------------------------------------------------------------------
1178 :
1179 : !!****f* m_io_screening/hscr_merge
1180 : !! NAME
1181 : !! hscr_merge
1182 : !!
1183 : !! FUNCTION
1184 : !! This subroutine merges different header structured variable (hscr_t)
1185 : !!
1186 : !! INPUTS
1187 : !! Hscr_in(:) <hscr_t)>=List of headers to be merged.
1188 : !!
1189 : !! OUTPUT
1190 : !! Hscr_out<hscr_t>=The output merged header.
1191 : !!
1192 : !! SOURCE
1193 :
1194 18 : subroutine hscr_merge(Hscr_in, Hscr_out)
1195 :
1196 : !Arguments ------------------------------------
1197 : !scalars
1198 : type(hscr_t),intent(in) :: Hscr_in(:)
1199 : type(hscr_t),intent(out) :: Hscr_out
1200 :
1201 : !Local variables-------------------------------
1202 : !scalars
1203 : integer :: nhds,restart,restartpaw,ihd,ii,nqtot,nqneq
1204 : logical :: isok
1205 : character(len=500) :: msg
1206 : !arrays
1207 3 : real(dp),allocatable :: qset(:,:)
1208 : ! *************************************************************************
1209 :
1210 : ! TODO: This should be called hscr_merge_qpoints
1211 :
1212 : !@hscr_t
1213 3 : nhds=SIZE(Hscr_in)
1214 :
1215 : ! Initial copy of the header ===
1216 : ! If multiple headers, select the header containing q-->0 so that we copy also heads and wings
1217 11 : ii = imax_loc(Hscr_in(:)%nqlwl)
1218 3 : call Hscr_in(ii)%copy(Hscr_out)
1219 3 : if (nhds==1) return
1220 :
1221 : ! Check consistency of the abinit Headers.
1222 : ! FFT grid might be q-point dependent so we stop only when restart==0
1223 3 : isok=.TRUE.
1224 8 : do ihd=2,nhds
1225 5 : call hdr_check(Hscr_in(1)%fform,Hscr_in(ihd)%fform,Hscr_in(1)%Hdr,Hscr_in(ihd)%Hdr,'COLL',restart,restartpaw)
1226 8 : if (restart==0) then
1227 0 : isok=.FALSE.
1228 0 : write(msg,'(a,i0,a)')' Abinit header no.',ihd,' is not consistent with the first header '
1229 0 : ABI_WARNING(msg)
1230 : end if
1231 : end do
1232 3 : if (.not.isok) then
1233 0 : ABI_ERROR('Cannot continue, Check headers')
1234 : end if
1235 :
1236 : ! Now check variables related to polarizability|epsilon^{-1}.
1237 : ! 1) Tests quantities that must be equal
1238 11 : ii = assert_eq(Hscr_in(:)%ID, 'Headers have different Identifiers')
1239 11 : ii = assert_eq(Hscr_in(:)%ikxc, 'Headers have different ikxc' )
1240 11 : ii = assert_eq(Hscr_in(:)%headform, 'Headers have different headform' )
1241 11 : ii = assert_eq(Hscr_in(:)%fform, 'Headers have different fform' )
1242 11 : ii = assert_eq(Hscr_in(:)%gwcalctyp,'Headers have different gwcalctyp' )
1243 11 : ii = assert_eq(Hscr_in(:)%nI, 'Headers have different nI' )
1244 11 : ii = assert_eq(Hscr_in(:)%nJ, 'Headers have different nJ' )
1245 11 : ii = assert_eq(Hscr_in(:)%nomega, 'Headers have different nomega' )
1246 11 : ii = assert_eq(Hscr_in(:)%test_type,'Headers have different test_type' )
1247 11 : ii = assert_eq(Hscr_in(:)%tordering,'Headers have different tordering' )
1248 :
1249 : ! This is not mandatory but makes life easier!
1250 11 : ii = assert_eq(Hscr_in(:)%npwe,'Headers have different number of G-vectors' )
1251 :
1252 8 : do ihd=2,nhds
1253 15 : if (ANY(ABS(Hscr_in(ihd)%omega-Hscr_in(1)%omega)>tol6)) then
1254 0 : write(msg,'(a,i0,a)')' Frequencies in the first and the ',ihd,'-th header differ'
1255 0 : ABI_ERROR(msg)
1256 : end if
1257 545 : if (ANY(Hscr_in(ihd)%gvec(:,:)-Hscr_in(1)%gvec(:,:)/=0)) then
1258 0 : write(msg,'(a,i0,a)')' Incompatible G-vector list found in the ',ihd,'-th header'
1259 0 : ABI_ERROR(msg)
1260 : end if
1261 8 : if (hscr_in(ihd)%kind_cdata /= hscr_in(1)%kind_cdata) then
1262 0 : write(msg,'(3a,i0,2a)')' Files contain data with different precisions.',ch10,&
1263 0 : "In particular the ",ihd,'-th header has precision:',trim(hscr_in(ihd)%kind_cdata)
1264 0 : ABI_ERROR(msg)
1265 : end if
1266 : end do !ihd
1267 :
1268 : ! If error is not fatal, just warn ===
1269 11 : if (ANY(Hscr_in(:)%npwwfn_used/=Hscr_in(1)%npwwfn_used)) then
1270 0 : ABI_COMMENT('Files have been produced with a different number of planewaves for the wavefunctions.')
1271 : end if
1272 11 : if (ANY(Hscr_in(:)%nbnds_used/=Hscr_in(1)%nbnds_used)) then
1273 0 : ABI_COMMENT('Files have been produced with a different number of bands.')
1274 : end if
1275 11 : if (ANY(Hscr_in(:)%spmeth/=Hscr_in(1)%spmeth)) then
1276 0 : ABI_COMMENT('Files have been produced with different algorithms.')
1277 : end if
1278 11 : if (ANY(ABS(Hscr_in(:)%mbpt_sciss-Hscr_in(1)%mbpt_sciss)>tol6)) then
1279 0 : ABI_COMMENT('Files have benn produced with different values of mbpt_sciss.')
1280 : end if
1281 11 : if (ANY(ABS(Hscr_in(:)%spsmear-Hscr_in(1)%spsmear)>tol6)) then
1282 0 : ABI_COMMENT('Files have been produced with different values of spsmear.')
1283 : end if
1284 11 : if (ANY(ABS(Hscr_in(:)%zcut-Hscr_in(1)%zcut)>tol6)) then
1285 0 : ABI_COMMENT('Files have been produced with different values of zcut.')
1286 : end if
1287 :
1288 : ! Now merge the list of q-points.
1289 : ! Take the union of the q-points, remove possible duplicated
1290 : ! are change the parameters in hscr_out that depends on q-points.
1291 11 : nqtot=SUM(Hscr_in(:)%nqibz)
1292 9 : ABI_MALLOC(qset,(3,nqtot))
1293 :
1294 3 : ii=0
1295 11 : do ihd=1,nhds
1296 92 : qset(:,ii+1:ii+Hscr_in(ihd)%nqibz)=Hscr_in(ihd)%qibz(:,:)
1297 11 : ii=ii+Hscr_in(ihd)%nqibz
1298 : end do
1299 :
1300 3 : call remove_copies(nqtot,qset,nqneq,isequalk)
1301 :
1302 3 : if (nqneq /= nqtot) then
1303 : write(msg,'(3a,2(i0,a))')&
1304 1 : 'COMMENT: Headers contain duplicated q-points ',ch10,&
1305 2 : 'Found ',nqneq,' distinct q-points among the total ',nqtot,' points reported in the headers. '
1306 1 : call wrtout(std_out, msg)
1307 : end if
1308 :
1309 3 : Hscr_out%nqibz = nqneq
1310 3 : ABI_FREE(Hscr_out%qibz)
1311 9 : ABI_MALLOC(Hscr_out%qibz,(3,nqneq))
1312 75 : Hscr_out%qibz(:,:)=qset(:,1:nqneq)
1313 3 : ABI_FREE(qset)
1314 :
1315 3 : end subroutine hscr_merge
1316 : !!***
1317 :
1318 : !----------------------------------------------------------------------
1319 :
1320 : !!****f* m_io_screening/write_screening
1321 : !! NAME
1322 : !! write_screening
1323 : !!
1324 : !! FUNCTION
1325 : !! For a single q-point, write either \tilde epsilon^{-1} on the _SCR file
1326 : !! or chi0 on the _SUSC file. The file is supposed to have been open in the calling routine.
1327 : !!
1328 : !! INPUTS
1329 : !! varname=The name of the array to write (used if etsf-io format).
1330 : !! unt=The unit number of the file to be written (supposed to be already open)
1331 : !! iomode=Integer flag defining the format of the output file. Available options:
1332 : !! IO_MODE_FORTRAN--> Plain Fortran file
1333 : !! IO_MODE_ETSF--> ETSF format
1334 : !! npwe=Number of plane waves in epsm1.
1335 : !! nomega=Number of frequencies
1336 : !! iq_ibz=Index of the q-points in the IBZ.
1337 : !! epsm1(npwe,npwe,nomega)=The matrix to be written, for different frequencies, and a single q-point.
1338 : !!
1339 : !! NOTES
1340 : !! On some architecture, the code crashes when trying to write or read a record containing the
1341 : !! entire (G1,G2) matrix thus we use smaller records containing the columns of the two-point function.
1342 : !!
1343 : !! OUTPUT
1344 : !! (only writing on file)
1345 : !!
1346 : !! SOURCE
1347 :
1348 812 : subroutine write_screening(varname, unt, iomode, npwe, nomega, iq_ibz, epsm1)
1349 :
1350 : !Arguments ------------------------------------
1351 : !scalars
1352 : character(len=*),intent(in) :: varname
1353 : integer,intent(in) :: nomega,npwe,iq_ibz,unt,iomode
1354 : !arrays
1355 : complex(gwp),target,intent(in) :: epsm1(npwe,npwe,nomega)
1356 :
1357 : !Local variables-------------------------------
1358 : !scalars
1359 : integer :: ipwe,iomega,spins(2),s1,s2
1360 : character(len=500) :: errmsg
1361 : !arrays
1362 812 : complex(dp),allocatable :: epsm1d(:,:)
1363 : integer :: varid,ncerr
1364 : #ifdef HAVE_GW_DPC
1365 812 : real(dp), contiguous, pointer :: real_epsm1(:,:,:,:,:,:,:)
1366 : #else
1367 : real(sp), contiguous, pointer :: real_epsm1(:,:,:,:,:,:,:)
1368 : #endif
1369 : ! *************************************************************************
1370 :
1371 : DBG_ENTER("COLL")
1372 :
1373 812 : select case (iomode)
1374 : case (IO_MODE_FORTRAN, IO_MODE_MPI)
1375 : ! Write a record for each omega, Always use double precision.
1376 0 : ABI_MALLOC(epsm1d,(npwe,1))
1377 :
1378 0 : do iomega=1,nomega
1379 0 : do ipwe=1,npwe
1380 0 : epsm1d(:,1) = epsm1(:,ipwe,iomega) !spc ==> dpc
1381 0 : write(unt, err=10, iomsg=errmsg)epsm1d(1:npwe,1)
1382 : end do
1383 : end do
1384 0 : ABI_FREE(epsm1d)
1385 :
1386 : case (IO_MODE_ETSF)
1387 : ! netcdf does not support complex datatypes. Here I use some C-magic to associate the memory
1388 : ! to a Fortran real pointer with the correct type and shape. Note that the data on file is always in double precision.
1389 : ! but this is ok since: if the type of data differs from the netCDF variable type, type conversion will occur
1390 : ! inside nf90_put_var
1391 812 : varid = nctk_idname(unt, varname)
1392 6496 : call c_f_pointer(c_loc(epsm1(1,1,1)), real_epsm1, [2, npwe, npwe, 1, 1, nomega, 1])
1393 : ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt]
1394 2436 : spins = 1; s1 = spins(1); s2 = spins(2)
1395 12180 : ncerr = nf90_put_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,iq_ibz], count=[2,npwe,npwe,1,1,nomega,1])
1396 812 : NCF_CHECK_MSG(ncerr, sjoin("putting var:", varname))
1397 :
1398 : case default
1399 0 : ABI_ERROR(sjoin("Wrong iomode:", iomode2str(iomode)))
1400 : end select
1401 :
1402 : DBG_EXIT("COLL")
1403 :
1404 : return
1405 :
1406 : ! Handle IO error
1407 : 10 continue
1408 0 : ABI_ERROR(errmsg)
1409 :
1410 812 : end subroutine write_screening
1411 : !!***
1412 :
1413 : !----------------------------------------------------------------------
1414 :
1415 : !!****f* m_io_screening/read_screening
1416 : !! NAME
1417 : !! read_screening
1418 : !!
1419 : !! FUNCTION
1420 : !! Read either a screening (\tilde epsilon^{-1}) file in the SCR format or
1421 : !! the irreducible polarizability (chi0) in the SUSC format.
1422 : !!
1423 : !! INPUTS
1424 : !! varname=Name of the array to read. Used for ETSF-IO files.
1425 : !! iomode=Integer flag defining the format of the output file. Available options:
1426 : !! IO_MODE_FORTRAN--> Plain Fortran file
1427 : !! IO_MODE_ETSF--> ETSF format
1428 : !! iqiA[optional]=Used if only a particular q-point is required. In this case iqiA define the index
1429 : !! of the required q-point in the array qibz(3,Hscr%nqibz)
1430 : !! nqibzA=number of asked q-points (used to dimension the output arrays).
1431 : !! Equal to Hscr%nqibz if the full matrix is required
1432 : !! comm=MPI communicator.
1433 : !! npweA=number of asked planewaves
1434 : !! nomegaA=number of asked frequencies
1435 : !!
1436 : !! OUTPUT
1437 : !! epsm1(npweA,npweA,nomegaA,nqibzA) = \tilde\epsilon^{-1}(Ng,Ng,Nw,Nq)
1438 : !!
1439 : !! NOTES
1440 : !! * If the epsilon matrix read is bigger than npweA x npweA, it will be truncated;
1441 : !! if it is smaller, an error will occur
1442 : !! * If the number of frequencies asked for is smaller than that reported in the file, the matrix
1443 : !! will be truncated. If nomegaA > Hscr%nomega an error will occur
1444 : !!
1445 : !! SOURCE
1446 :
1447 580 : subroutine read_screening(varname, fname, npweA, nqibzA, nomegaA, epsm1, iomode, comm, &
1448 : iqiA) ! Optional
1449 :
1450 : !Arguments ------------------------------------
1451 : !scalars
1452 : integer,intent(in) :: iomode,nomegaA,npweA,nqibzA,comm
1453 : integer,optional,intent(in) :: iqiA
1454 : character(len=*),intent(in) :: varname,fname
1455 : !arrays
1456 : complex(gwp),target,intent(inout) :: epsm1(npweA,npweA,nomegaA,nqibzA)
1457 :
1458 : !Local variables-------------------------------
1459 : !scalars
1460 : integer,parameter :: master = 0
1461 : integer :: ipwe,fform,iomega,iq_ibz,unt,rdwr,my_rank,nprocs,my_iomode, varid, ncerr
1462 : #ifdef HAVE_MPI_IO
1463 : integer :: test_fform,mpi_err,ierr,sc_mode, bsize_frm,mpi_type_frm
1464 : integer :: mpi_fh,buf_dim !,mat_ggw,mat_ggwq
1465 : integer(XMPI_OFFSET_KIND) :: offset,displ_wq !,my_offpad
1466 : !complex(dp) :: ctmp
1467 : #endif
1468 : real(dp) :: cpu, wall, gflops
1469 : character(len=500) :: msg,errmsg
1470 : logical :: read_qslice
1471 2320 : type(hscr_t) :: Hscr
1472 : !arrays
1473 : #ifdef HAVE_MPI_IO
1474 580 : integer(MPI_OFFSET_KIND),allocatable :: offset_wq(:,:)
1475 : #endif
1476 580 : complex(dp),allocatable :: bufdc2d(:,:),bufdc3d(:,:,:)
1477 : ! pointers passed to netcdf4 routines (complex datatypes are not supported).
1478 : #ifdef HAVE_GW_DPC
1479 580 : real(dp), contiguous, pointer :: real_epsm1(:,:,:,:,:,:,:)
1480 : #else
1481 : real(sp), contiguous, pointer :: real_epsm1(:,:,:,:,:,:,:)
1482 : #endif
1483 : integer :: spins(2),s1,s2
1484 : ! *************************************************************************
1485 :
1486 : DBG_ENTER("COLL")
1487 580 : call cwtime(cpu, wall, gflops, "start")
1488 :
1489 580 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1490 580 : my_iomode = iomode
1491 580 : if (endswith(fname, ".nc")) my_iomode = IO_MODE_ETSF
1492 : !my_iomode = IO_MODE_MPI
1493 : !if (my_iomode == IO_MODE_MPI) my_iomode = IO_MODE_FORTRAN
1494 :
1495 580 : rdwr=1
1496 0 : select case (my_iomode)
1497 : case (IO_MODE_MPI)
1498 : #ifdef HAVE_MPI_IO
1499 0 : bsize_frm = xmpio_bsize_frm ! bsize_frm= Byte length of the Fortran record marker.
1500 0 : mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker.
1501 0 : sc_mode = xmpio_collective
1502 :
1503 : ! Master reads the header via Fortran IO then bcast the data.
1504 0 : call hscr%from_file(fname, fform, comm)
1505 :
1506 : ! Open the file with MPI-IO
1507 0 : call MPI_FILE_OPEN(comm, fname, MPI_MODE_RDONLY, xmpio_info ,mpi_fh, mpi_err)
1508 0 : ABI_CHECK_MPI(mpi_err, sjoin("MPI_FILE_OPEN:", fname))
1509 :
1510 : ! Retrieve the offset of the section immediately below the header.
1511 0 : call hscr_mpio_skip(mpi_fh, test_fform, offset)
1512 0 : ABI_CHECK(test_fform == fform, "mismatch in fform!")
1513 :
1514 : ! Offsets of the Fortran markers corresponding to the (w,q) slices.
1515 0 : ABI_MALLOC(offset_wq,(HScr%nomega,HScr%nqibz))
1516 0 : displ_wq = offset
1517 0 : do iq_ibz=1,Hscr%nqibz
1518 0 : do iomega=1,Hscr%nomega
1519 0 : ABI_CHECK(displ_wq > 0, "displ_wq < 0, your SCR|SUSC file is too big for MPI-IO!")
1520 0 : offset_wq(iomega,iq_ibz) = displ_wq
1521 0 : displ_wq = displ_wq + Hscr%npwe**2 * xmpi_bsize_dpc + Hscr%npwe * 2 * bsize_frm
1522 : end do
1523 : end do
1524 : #else
1525 : ABI_ERROR("MPI-IO support not enabled at configure-time")
1526 : #endif
1527 :
1528 : case (IO_MODE_FORTRAN)
1529 : ! Plain Fortran IO, all nodes read.
1530 0 : if (open_file(fname,msg,newunit=unt,form="unformatted",status="old",action="read") /= 0) then
1531 0 : ABI_ERROR(msg)
1532 : end if
1533 0 : call hscr%io(fform, rdwr, unt, comm, master, my_iomode)
1534 :
1535 : case (IO_MODE_ETSF)
1536 580 : NCF_CHECK(nctk_open_read(unt, fname, xmpi_comm_self))
1537 580 : call hscr%io(fform, rdwr, unt, comm, master, my_iomode)
1538 :
1539 : case default
1540 580 : ABI_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode)))
1541 : end select
1542 :
1543 : ! Slice or full array?
1544 580 : read_qslice = .False.
1545 580 : if (PRESENT(iqiA)) then
1546 451 : read_qslice = .True.
1547 : !call wrtout(std_out, sjoin('. Reading q-slice for iq = ',itoa(iqiA),' from: ', fname))
1548 451 : if (iqiA <= 0 .or. iqiA > Hscr%nqibz) then
1549 0 : ABI_BUG('iqiA out of range')
1550 : end if
1551 : end if
1552 :
1553 : ! Do some check
1554 580 : if (Hscr%npwe>npweA) then
1555 : write(msg,'(a,i0,2a,i0)')&
1556 65 : 'Total number of G-vectors reported on file = ',Hscr%npwe,ch10,&
1557 130 : 'Reading a smaller matrix of dimension = ',npweA
1558 65 : ABI_COMMENT(msg)
1559 : end if
1560 :
1561 580 : if (npweA > Hscr%npwe) then
1562 0 : write(msg,'(2(a,i0))')' Dimension of matrix = ',Hscr%npwe," requiring a too big matrix = ",npweA
1563 0 : ABI_ERROR(msg)
1564 : end if
1565 :
1566 580 : ABI_CHECK(nqibzA <= Hscr%nqibz, 'Requiring too many q-points')
1567 580 : ABI_CHECK(nomegaA <= Hscr%nomega,'Requiring too many frequencies')
1568 :
1569 0 : select case (my_iomode)
1570 : case (IO_MODE_MPI)
1571 : #ifdef HAVE_MPI_IO
1572 0 : if (read_qslice) then
1573 : !call wrtout(std_out, "calling mpiotk to read_qslice")
1574 0 : buf_dim = (npweA)**2 * nomegaA
1575 0 : offset = offset_wq(1,iqiA)
1576 0 : sc_mode = xmpio_collective
1577 :
1578 : #ifdef HAVE_GW_DPC
1579 : ! Read in-place.
1580 : call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1581 0 : buf_dim,epsm1,xmpio_chunk_bsize,sc_mode,comm,ierr)
1582 0 : ABI_CHECK(ierr==0,"Fortran matrix too big")
1583 : #else
1584 : ! Have to allocate workspace for dp data.
1585 : ! FIXME: Change the file format of the SCR and SUC file so that
1586 : ! they are written in single precision if not HAVE_GW_DPC
1587 : ABI_MALLOC_OR_DIE(bufdc3d, (npweA,npweA,nomegaA), ierr)
1588 :
1589 : call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1590 : buf_dim,bufdc3d,xmpio_chunk_bsize,sc_mode,comm,ierr)
1591 : ABI_CHECK(ierr == 0,"Fortran matrix too big")
1592 :
1593 : epsm1(:,:,:,1) = bufdc3d
1594 : ABI_FREE(bufdc3d)
1595 : #endif
1596 :
1597 : else
1598 : ! Full matrix (G,G',w,q) is needed.
1599 : !call wrtout(std_out, "calling mpiotk: Full matrix (G,G',w,q) is needed.")
1600 :
1601 : #ifdef HAVE_GW_DPC
1602 : ! Can read all data at once.
1603 0 : buf_dim = (npweA)**2 * nomegaA * HScr%nqibz
1604 0 : offset = offset_wq(1,1)
1605 0 : sc_mode = xmpio_collective
1606 :
1607 : call mpiotk_read_fsuba_dpc4D(mpi_fh,offset,&
1608 : [HScr%npwe,HScr%npwe,HScr%nomega,HScr%nqibz], [npweA,npweA,nomegaA,HScr%nqibz], [1,1,1,1],&
1609 0 : buf_dim,epsm1,xmpio_chunk_bsize,sc_mode,comm,ierr)
1610 0 : ABI_CHECK(ierr==0,"Fortran record too big")
1611 : #else
1612 : ! Have to allocate workspace for dp data.
1613 : ABI_MALLOC_OR_DIE(bufdc3d,(npweA,npweA,nomegaA), ierr)
1614 : sc_mode = xmpio_collective
1615 :
1616 : do iq_ibz=1,Hscr%nqibz
1617 : offset = offset_wq(1,iq_ibz)
1618 : buf_dim = (2*npweA)**2 * nomegaA
1619 :
1620 : call mpiotk_read_fsuba_dpc3D(mpi_fh,offset, &
1621 : [HScr%npwe,HScr%npwe,HScr%nomega], [npweA,npweA,nomegaA], [1,1,1],&
1622 : buf_dim,bufdc3d,xmpio_chunk_bsize,sc_mode,comm,ierr)
1623 : ABI_CHECK(ierr==0,"Fortran matrix too big")
1624 :
1625 : epsm1(:,:,:,iq_ibz) = bufdc3d
1626 : end do
1627 :
1628 : ABI_FREE(bufdc3d)
1629 : #endif
1630 : end if
1631 :
1632 0 : call MPI_FILE_CLOSE(mpi_fh,mpi_err)
1633 0 : ABI_FREE(offset_wq)
1634 : #endif
1635 :
1636 : case (IO_MODE_FORTRAN)
1637 : ! Read epsilon^-1 with Fortran IO
1638 : ! Allocate a single column to save memory.
1639 : ! TODO re-merge the two cases.
1640 0 : ABI_MALLOC(bufdc2d,(Hscr%npwe,1))
1641 :
1642 : ! Two coding for different case just to keep it readable.
1643 : select case (read_qslice)
1644 : case (.True.)
1645 : ! Read only a slice of the full array (useful if the entire array is huge).
1646 : !if (dim_wings==1) STOP 'not implemented'
1647 : !TODO this has to be done in a cleaner way.
1648 : qread_loop: &
1649 0 : & do iq_ibz=1,Hscr%nqibz
1650 0 : if (iq_ibz==iqiA) then
1651 0 : do iomega=1,nomegaA
1652 0 : do ipwe=1,Hscr%npwe
1653 0 : read(unt, err=10, iomsg=errmsg) bufdc2d(1:Hscr%npwe,1)
1654 0 : if (ipwe<=npweA) epsm1(1:npweA,ipwe,iomega,1)=bufdc2d(1:npweA,1)
1655 : end do
1656 : end do
1657 : EXIT qread_loop ! Got data. Do not need to read file till the end.
1658 : else
1659 : ! Skip other q-points i.e bufdc2d(1:Hscr%npwe,1:Hscr%npwe)
1660 0 : do iomega=1,Hscr%nomega
1661 0 : do ipwe=1,Hscr%npwe
1662 0 : read(unt, err=10, iomsg=errmsg)
1663 : end do
1664 : end do
1665 : end if ! iq_ibz==iqiA
1666 : end do qread_loop ! iq_ibz
1667 :
1668 : case (.False.)
1669 : ! Read the entire array.
1670 0 : do iq_ibz=1,Hscr%nqibz
1671 0 : do iomega=1,nomegaA
1672 0 : do ipwe=1,Hscr%npwe
1673 0 : read(unt, err=10, iomsg=errmsg) bufdc2d(1:Hscr%npwe,1)
1674 0 : if (ipwe<=npweA) epsm1(1:npweA,ipwe,iomega,iq_ibz)=bufdc2d(1:npweA,1)
1675 : end do
1676 : end do
1677 : ! Skip other frequencies
1678 0 : do iomega=nomegaA+1,Hscr%nomega
1679 0 : do ipwe=1,Hscr%npwe
1680 0 : read(unt, err=10, iomsg=errmsg)
1681 : end do
1682 : end do
1683 : end do !iq_ibz
1684 : end select
1685 :
1686 0 : close(unt)
1687 :
1688 : case (IO_MODE_ETSF)
1689 : ! netcdf does not support complex datatypes. Here I use some C-magic to associate the memory
1690 : ! to a Fortran real pointer with the correct type and shape. Note that the data on file is always in double precision.
1691 : ! nf90_get_var will automatically convert from double to single if the GW code is in single precision mode.
1692 : ! This is the reason why I'm using CPP option in the declaration of real_epsm1.
1693 :
1694 : ! FIXME: Need to know the type to read
1695 : !write(std_out,*)"in read_screening"
1696 580 : varid = nctk_idname(unt, varname)
1697 :
1698 : ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt]
1699 4640 : call c_f_pointer(c_loc(epsm1(1,1,1,1)), real_epsm1, [2,npweA,npweA,1,1,nomegaA,nqibzA])
1700 1740 : spins = 1; s1 = spins(1); s2 = spins(2)
1701 580 : if (read_qslice) then
1702 6765 : ncerr = nf90_get_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,iqia], count=[2,npweA,npweA,1,1,nomegaA,1])
1703 : else
1704 1935 : ncerr = nf90_get_var(unt, varid, real_epsm1, start=[1,1,1,s1,s2,1,1], count=[2,npweA,npweA,1,1,nomegaA,nqibzA])
1705 : !do iq_ibz=1,nqibzA
1706 : ! !write(*,*)"epsm1: in read ",iq_ibz,epsm1(1:3,1,1,iq_ibz)
1707 : !end do
1708 : end if
1709 580 : NCF_CHECK_MSG(ncerr, sjoin("getting var:", varname))
1710 580 : NCF_CHECK(nf90_close(unt))
1711 : !write(std_out,*)"read_screening done"
1712 :
1713 : case default
1714 580 : ABI_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode)))
1715 : end select
1716 :
1717 : ! Free memory
1718 580 : ABI_SFREE(bufdc2d)
1719 : ABI_SFREE(bufdc3d)
1720 :
1721 580 : call Hscr%free()
1722 :
1723 : !call cwtime_report("- read_screening:", cpu, wall, gflops)
1724 : DBG_EXIT("COLL")
1725 :
1726 580 : return
1727 :
1728 : ! Handle Fortran IO error.
1729 : 10 continue
1730 0 : ABI_ERROR(errmsg)
1731 :
1732 1160 : end subroutine read_screening
1733 : !!***
1734 :
1735 : !----------------------------------------------------------------------
1736 :
1737 : !!****f* m_io_screening/hscr_mpio_skip
1738 : !! NAME
1739 : !! hscr_mpio_skip
1740 : !!
1741 : !! FUNCTION
1742 : !! Skip the header of the (SCR|SUSC) file in MPI-IO mode. This routine uses local MPI-IO calls hence
1743 : !! it can be safely called by master node only. Note however that in this case the
1744 : !! offset has to be communicated to the other nodes.
1745 : !!
1746 : !! INPUTS
1747 : !! mpio_fh=MPI-IO file handler
1748 : !! fmarker_bsize = Byte length of Fortran record marker.
1749 : !! fmarker_mpi_type= MPI type of the Fortran record marker
1750 : !!
1751 : !! OUTPUT
1752 : !! fform=kind of the array in the file
1753 : !! offset=The offset of the Fortran record located immediately below the Abinit header.
1754 : !!
1755 : !! SOURCE
1756 :
1757 0 : subroutine hscr_mpio_skip(mpio_fh, fform, offset)
1758 :
1759 : !Arguments ------------------------------------
1760 : integer,intent(in) :: mpio_fh
1761 : integer,intent(out) :: fform
1762 : integer(kind=XMPI_OFFSET_KIND),intent(out) :: offset
1763 :
1764 : !Local variables-------------------------------
1765 : !scalars
1766 : integer :: bsize_frm,mpi_type_frm, headform
1767 : #ifdef HAVE_MPI_IO
1768 : integer :: ierr,isk
1769 : !character(len=500) :: msg
1770 : !arrays
1771 : integer(kind=MPI_OFFSET_KIND) :: fmarker,positloc
1772 : integer :: nqlwl(1),statux(MPI_STATUS_SIZE)
1773 : #endif
1774 : ! *************************************************************************
1775 :
1776 : offset = 0
1777 0 : bsize_frm = xmpio_bsize_frm ! Byte size of the Fortran record marker.
1778 0 : mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker.
1779 :
1780 0 : call hdr_mpio_skip(mpio_fh,fform,offset)
1781 :
1782 : !call wrtout(std_out, sjoin("in hdr_mpio_skip with fform = ",itoa(fform)))
1783 :
1784 : #ifdef HAVE_MPI_IO
1785 0 : select case (fform)
1786 : case (1003, 1004)
1787 : ! Skip the titles
1788 0 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1789 :
1790 : ! read headform from the 2d record.
1791 0 : positloc = offset + bsize_frm + 3*xmpi_bsize_int
1792 0 : call MPI_FILE_READ_AT(mpio_fh,positloc,headform,1,MPI_INTEGER,statux,ierr)
1793 : !call wrtout(std_out, sjoin("headform = ",itoa(headform)))
1794 :
1795 : ! read nqlwl from the 2d record.
1796 0 : positloc = offset + bsize_frm + 9*xmpi_bsize_int
1797 0 : call MPI_FILE_READ_AT(mpio_fh,positloc,nqlwl,1,MPI_INTEGER,statux,ierr)
1798 : !call wrtout(std_out, sjoin("nqlwl = ",itoa(nqlwl(1))))
1799 :
1800 0 : do isk=1,5
1801 0 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1802 : end do
1803 :
1804 0 : if (nqlwl(1)>0) then ! skip qlwl
1805 0 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1806 : end if
1807 :
1808 0 : if (headform > 80) then ! skip new records
1809 0 : do isk=1,4
1810 0 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1811 : end do
1812 : end if
1813 :
1814 : case default
1815 0 : ABI_BUG(sjoin('Wrong fform read:', itoa(fform)))
1816 : end select
1817 :
1818 : #else
1819 : ABI_ERROR("hscr_mpio_skip cannot be used when MPI-IO is not enabled")
1820 : #endif
1821 :
1822 0 : end subroutine hscr_mpio_skip
1823 : !!***
1824 :
1825 : !----------------------------------------------------------------------
1826 :
1827 : !!****f* m_io_screening/ioscr_qmerge
1828 : !! NAME
1829 : !! ioscr_qmerge
1830 : !!
1831 : !! FUNCTION
1832 : !! Produce new file by merging the q-points stored in other files.
1833 : !! This routine should be called by a single MPI process.
1834 : !!
1835 : !! INPUTS
1836 : !! nfiles=Number of files to be merged.
1837 : !! filenames(nfiles)=Paths of files to be merged.
1838 : !! hscr_files(nfiles)<hscr_t>=Headers of the files to be merged.
1839 : !! fname_out=Name of the file to be produced.
1840 : !!
1841 : !! OUTPUT
1842 : !! ohscr<hscr_t>=The header of the output file.
1843 : !!
1844 : !! SOURCE
1845 :
1846 15 : subroutine ioscr_qmerge(nfiles, filenames, hscr_files, fname_out, ohscr)
1847 :
1848 : !Arguments ------------------------------------
1849 : !scalars
1850 : integer,intent(in) :: nfiles
1851 : character(len=*),intent(in) :: fname_out
1852 : type(hscr_t),intent(out) :: ohscr
1853 : !arrays
1854 : character(len=*),intent(in) :: filenames(nfiles)
1855 : type(hscr_t),intent(in) :: hscr_files(nfiles)
1856 :
1857 : !Local variables-------------------------------
1858 : !scalars
1859 : integer,parameter :: rdwr2=2,master=0
1860 : integer :: iq_ibz,ifound,ifile,iqf,ount,iomode,fform_merge,comm,nomega4m,npwe4m,iqiA,ierr
1861 : character(len=500) :: msg
1862 : character(len=nctk_slen) :: varname
1863 : type(abifile_t) :: abifile
1864 : !arrays
1865 3 : integer,allocatable :: merge_table(:,:)
1866 : real(dp) :: qdiff(3)
1867 3 : complex(gwp),allocatable :: epsm1(:,:,:,:)
1868 : ! *************************************************************************
1869 :
1870 3 : comm = xmpi_comm_self
1871 :
1872 3 : if (file_exists(fname_out)) then
1873 0 : ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
1874 : end if
1875 :
1876 : ! Merge the headers creating the full list of q-points.
1877 3 : call hscr_merge(Hscr_files(1:nfiles), ohscr)
1878 6 : call ohscr%print([std_out], 1, header='Header of the final file')
1879 :
1880 : ! For each q to be merged, save the index of the file where q is stored as well as its sequential index.
1881 : ! Useful to do the merge point-by-point thus avoiding the allocation of the entire epsm1 array.
1882 9 : ABI_MALLOC(merge_table,(ohscr%nqibz,2))
1883 21 : do iq_ibz=1,ohscr%nqibz
1884 34 : ifound=0
1885 34 : fl: do ifile=1,nfiles
1886 79 : do iqf=1,Hscr_files(ifile)%nqibz
1887 252 : qdiff(:)=ohscr%qibz(:,iq_ibz)-Hscr_files(ifile)%qibz(:,iqf)
1888 148 : if (all(abs(qdiff) < GW_TOLQ)) then
1889 18 : merge_table(iq_ibz,1)=ifile
1890 18 : merge_table(iq_ibz,2)=iqf
1891 18 : ifound=ifound+1
1892 18 : write(msg,'(a,3f12.6,2a)')'. q-point:',ohscr%qibz(:,iq_ibz),' will be taken from ',TRIM(filenames(ifile))
1893 18 : call wrtout(std_out, msg)
1894 : EXIT fl
1895 : end if
1896 : end do
1897 : end do fl
1898 : ! Check if q-point has been found, multiple q-points not allowed.
1899 3 : ABI_CHECK(ifound == 1, 'ifound/=1')
1900 : end do
1901 :
1902 3 : iomode = IO_MODE_FORTRAN; if (endswith(fname_out, ".nc")) iomode = IO_MODE_ETSF
1903 3 : if (iomode == IO_MODE_FORTRAN) then
1904 0 : if (open_file(fname_out,msg,newunit=ount,status='new',form='unformatted') /= 0) then
1905 0 : ABI_ERROR(msg)
1906 : end if
1907 : else
1908 3 : NCF_CHECK(nctk_open_create(ount, fname_out, comm))
1909 : end if
1910 :
1911 : ! Write the header.
1912 3 : fform_merge = hscr_files(1)%fform
1913 3 : abifile = abifile_from_fform(fform_merge)
1914 3 : if (abifile%fform == 0) then
1915 0 : ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
1916 : end if
1917 3 : varname = abifile%varname
1918 :
1919 11 : if (any(hscr_files(:)%fform /= hscr_files(1)%fform)) then
1920 0 : write(std_out,*)"fforms: ",hscr_files(:)%fform
1921 0 : ABI_ERROR("Files to be merged have different fform. Cannot merge data")
1922 : end if
1923 :
1924 3 : call ohscr%io(fform_merge, rdwr2, ount, comm, master, iomode)
1925 :
1926 3 : npwe4m = ohscr%npwe
1927 3 : nomega4m = ohscr%nomega
1928 :
1929 15 : ABI_MALLOC_OR_DIE(epsm1,(npwe4m,npwe4m,nomega4m,1), ierr)
1930 :
1931 21 : do iq_ibz=1,ohscr%nqibz
1932 18 : ifile = merge_table(iq_ibz,1)
1933 18 : iqiA = merge_table(iq_ibz,2)
1934 18 : call read_screening(varname, filenames(ifile), npwe4m, 1, nomega4m,epsm1, iomode, comm, iqiA=iqiA)
1935 21 : call write_screening(varname, ount, iomode, npwe4m, nomega4m, iq_ibz, epsm1)
1936 : end do
1937 :
1938 3 : ABI_FREE(epsm1)
1939 3 : ABI_FREE(merge_table)
1940 :
1941 3 : if (iomode == IO_MODE_FORTRAN) then
1942 0 : close(ount)
1943 : else
1944 3 : NCF_CHECK(nf90_close(ount))
1945 : end if
1946 :
1947 3 : write(msg,'(3a)')ch10,' ==== Files have been merged successfully === ',ch10
1948 3 : call wrtout(std_out, msg)
1949 :
1950 3 : end subroutine ioscr_qmerge
1951 : !!***
1952 :
1953 : !----------------------------------------------------------------------
1954 :
1955 : !!****f* m_io_screening/ioscr_qrecover
1956 : !! NAME
1957 : !! ioscr_qrecover
1958 : !!
1959 : !! FUNCTION
1960 : !! Recover q-points from a corrupted file produced e.g. from an interrupted run
1961 : !! This routine should be called by a single MPI process.
1962 : !!
1963 : !! INPUTS
1964 : !! path=Corrupted file.
1965 : !! nqrec=Number of q-points to recover.
1966 : !! fname_out=Name of the file to be produced.
1967 : !!
1968 : !! OUTPUT
1969 : !! Output is written to file.
1970 : !!
1971 : !! SOURCE
1972 :
1973 1 : subroutine ioscr_qrecover(ipath, nqrec, fname_out)
1974 :
1975 : !Arguments ------------------------------------
1976 : !scalars
1977 : integer,intent(in) :: nqrec
1978 : character(len=*),intent(in) :: ipath,fname_out
1979 :
1980 : !Local variables-------------------------------
1981 : !scalars
1982 : integer,parameter :: rdwr2=2,master=0
1983 : integer :: iqiA,nqibzA,nomega_asked,unt,npwe_asked,iomode,comm,fform1,ifform,ierr
1984 : character(len=500) :: msg
1985 : character(len=nctk_slen) :: varname
1986 7 : type(hscr_t) :: hscr_recov,hscr
1987 : type(abifile_t) :: abifile
1988 : !arrays
1989 1 : complex(gwp),allocatable :: epsm1(:,:,:,:)
1990 : ! *************************************************************************
1991 :
1992 1 : comm = xmpi_comm_self
1993 :
1994 1 : call wrtout(std_out, sjoin(". Recovering q-points in file:", ipath))
1995 1 : call wrtout(std_out, sjoin(". Data written to file:", fname_out))
1996 :
1997 1 : if (file_exists(fname_out)) then
1998 0 : ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
1999 : end if
2000 :
2001 : ! Find iomode from file extension and open output file.
2002 1 : if (endswith(fname_out, ".nc")) then
2003 1 : iomode = IO_MODE_ETSF
2004 1 : NCF_CHECK(nctk_open_create(unt, fname_out, comm))
2005 : else
2006 0 : iomode = IO_MODE_FORTRAN
2007 0 : if (open_file(fname_out, msg, newunit=unt, status='new', form='unformatted') /= 0) then
2008 0 : ABI_ERROR(msg)
2009 : end if
2010 : end if
2011 :
2012 : ! Read header.
2013 1 : call hscr%from_file(ipath, ifform, comm)
2014 1 : ABI_CHECK(ifform /= 0, sjoin("fform = 0 while reading:", ipath))
2015 :
2016 1 : if (nqrec < 1 .or. nqrec > hscr%nqibz) then
2017 0 : ABI_ERROR(sjoin("Wrong input. nqibz on file:", itoa(hscr%nqibz)))
2018 : end if
2019 :
2020 : ! Copy header
2021 1 : call hscr%copy(hscr_recov)
2022 :
2023 : ! Change dimensions and arrays associated to nqibz.
2024 1 : hscr_recov%nqibz = nqrec
2025 1 : ABI_FREE(hscr_recov%qibz)
2026 3 : ABI_MALLOC(hscr_recov%qibz, (3,nqrec))
2027 14 : hscr_recov%qibz = hscr%qibz(:,1:nqrec)
2028 :
2029 2 : call hscr_recov%print([std_out], 1, header="Header of the new SCR file")
2030 :
2031 : ! Write the header of the recovered file.
2032 1 : fform1 = hscr%fform
2033 :
2034 1 : abifile = abifile_from_fform(fform1)
2035 1 : if (abifile%fform == 0) then
2036 0 : ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform1)))
2037 : end if
2038 1 : varname = abifile%varname
2039 :
2040 1 : call hscr_recov%io(fform1,rdwr2,unt,comm,master,iomode)
2041 :
2042 1 : nqibzA=1; nomega_asked=hscr%nomega; npwe_asked=hscr%npwe
2043 :
2044 5 : ABI_MALLOC_OR_DIE(epsm1,(npwe_asked,npwe_asked,nomega_asked,1), ierr)
2045 :
2046 4 : do iqiA=1,hscr_recov%nqibz
2047 3 : call read_screening(varname,ipath,npwe_asked,nqibzA,nomega_asked,epsm1,iomode,comm,iqiA=iqiA)
2048 4 : call write_screening(varname,unt,iomode,npwe_asked,nomega_asked,iqiA,epsm1)
2049 : end do
2050 :
2051 1 : if (iomode == IO_MODE_FORTRAN) close(unt)
2052 1 : if (iomode == IO_MODE_ETSF) then
2053 1 : NCF_CHECK(nf90_close(unt))
2054 : end if
2055 :
2056 1 : ABI_FREE(epsm1)
2057 1 : call hscr%free()
2058 1 : call hscr_recov%free()
2059 :
2060 1 : call wrtout(std_out, "Recovery completed")
2061 :
2062 2 : end subroutine ioscr_qrecover
2063 : !!***
2064 :
2065 : !----------------------------------------------------------------------
2066 :
2067 : !!****f* m_io_screening/ioscr_wmerge
2068 : !! NAME
2069 : !! ioscr_wmerge
2070 : !!
2071 : !! FUNCTION
2072 : !! Produce new file by merging the frequencies stored in other files.
2073 : !! This routine should be called by a single MPI process.
2074 : !!
2075 : !! INPUTS
2076 : !! nfiles=Number of files to be merged.
2077 : !! filenames(nfiles)=Paths of files to be merged.
2078 : !! hscr_files(nfiles)<hscr_t>=Headers of the files to be merged.
2079 : !! fname_out=Name of the file to be produced.
2080 : !!
2081 : !! OUTPUT
2082 : !! ohscr<hscr_t>=The header of the output file.
2083 : !!
2084 : !! SOURCE
2085 :
2086 15 : subroutine ioscr_wmerge(nfiles, filenames, hscr_file, freqremax, fname_out, ohscr)
2087 :
2088 : !Arguments ------------------------------------
2089 : !scalars
2090 : integer,intent(in) :: nfiles
2091 : real(dp),intent(in) :: freqremax
2092 : character(len=*),intent(in) :: fname_out
2093 : type(hscr_t),intent(out) :: ohscr
2094 : !arrays
2095 : character(len=*),intent(in) :: filenames(nfiles)
2096 : type(hscr_t),intent(in) :: hscr_file(nfiles)
2097 :
2098 : !Local variables-------------------------------
2099 : !scalars
2100 : integer,parameter :: rdwr2=2,master=0
2101 : integer :: ii,iq_ibz,ifile,ount,iomode,fform_merge,comm,nomega4m
2102 : integer :: nfreq_tot,nfreqre,nfreqim,ifrq,npwe4mI,npwe4mJ,ierr, mod10, order
2103 : character(len=500) :: msg
2104 : logical :: skip
2105 : character(len=nctk_slen) :: varname
2106 : character(len=etsfio_charlen) :: rw_mesh_type, iw_mesh_type, cw_mesh_type
2107 : type(abifile_t) :: abifile
2108 : !arrays
2109 3 : integer,allocatable :: freq_indx(:,:),ifile_indx(:),pos_indx(:),i_temp(:),i2_temp(:,:)
2110 3 : real(dp),allocatable :: real_omega(:), real_omega_wgs(:), imag_omega(:), imag_omega_wgs(:) ,omega_wgs_storage(:)
2111 3 : complex(gwp),allocatable :: epsm1(:,:,:,:),epsm1_temp(:,:,:,:)
2112 3 : complex(dp),allocatable :: omega_storage(:)
2113 : ! *************************************************************************
2114 :
2115 3 : comm = xmpi_comm_self
2116 :
2117 3 : if (file_exists(fname_out)) then
2118 0 : ABI_ERROR(sjoin("Cannot overwrite existing file:", fname_out))
2119 : end if
2120 :
2121 3 : mod10 = mod(hscr_file(1)%gwcalctyp, 10)
2122 :
2123 : ! These are the meshes that will be reported in the final SCR file.
2124 : ! None means that we have a single point or no point at all along that axis
2125 : ! Clearly we cannot merge files with different meshes.
2126 3 : rw_mesh_type = "None"; iw_mesh_type = "None"; cw_mesh_type = "None"
2127 9 : do ifile=1,nfiles
2128 6 : if (hscr_file(ifile)%rw_mesh_type /= "None") then
2129 4 : if (rw_mesh_type /= "None" .and. rw_mesh_type /= hscr_file(ifile)%rw_mesh_type) then
2130 0 : ABI_ERROR(sjoin("Cannot merge rw_mesh_type:", rw_mesh_type, " with: ", hscr_file(ifile)%rw_mesh_type))
2131 : end if
2132 4 : rw_mesh_type = hscr_file(ifile)%rw_mesh_type
2133 : end if
2134 :
2135 6 : if (hscr_file(ifile)%iw_mesh_type /= "None") then
2136 5 : if (iw_mesh_type /= "None" .and. iw_mesh_type /= hscr_file(ifile)%iw_mesh_type) then
2137 0 : ABI_ERROR(sjoin("Cannot merge iw_mesh_type:", iw_mesh_type, " with: ", hscr_file(ifile)%iw_mesh_type))
2138 : end if
2139 5 : iw_mesh_type = hscr_file(ifile)%iw_mesh_type
2140 : end if
2141 :
2142 9 : if (hscr_file(ifile)%cw_mesh_type /= "None") then
2143 0 : if (cw_mesh_type /= "None" .and. cw_mesh_type /= hscr_file(ifile)%cw_mesh_type) then
2144 0 : ABI_ERROR(sjoin("Cannot merge cw_mesh_type:", cw_mesh_type, " with: ", hscr_file(ifile)%cw_mesh_type))
2145 : end if
2146 0 : cw_mesh_type = hscr_file(ifile)%cw_mesh_type
2147 : end if
2148 : end do ! ifile
2149 :
2150 6 : do ifile=2,nfiles
2151 :
2152 : ! Check that q-points are the same
2153 3 : ABI_CHECK_IEQ(Hscr_file(ifile)%nqibz, Hscr_file(1)%nqibz, 'files do not have the same number of q-points!')
2154 23 : do iq_ibz=1,Hscr_file(1)%nqibz
2155 83 : if (ABS(SUM(Hscr_file(ifile)%qibz(:,iq_ibz) - Hscr_file(1)%qibz(:,iq_ibz))) > tol6) then
2156 0 : ABI_ERROR('Q-point set differs between one or more files!')
2157 : end if
2158 : end do
2159 :
2160 : ! Check gwcalctyp and other basic parameters.
2161 3 : ABI_CHECK_IEQ(hscr_file(ifile)%gwcalctyp, hscr_file(1)%gwcalctyp, "Different gwcalctyp")
2162 3 : ABI_CHECK_IEQ(hscr_file(ifile)%npwe, hscr_file(1)%npwe, "Different npwe")
2163 3 : ABI_CHECK_IEQ(hscr_file(ifile)%id, hscr_file(1)%id, "Different IDs")
2164 3 : ABI_CHECK_IEQ(hscr_file(ifile)%ikxc, hscr_file(1)%ikxc, "Different IXCx")
2165 3 : ABI_CHECK_IEQ(hscr_file(ifile)%test_type, hscr_file(1)%test_type, "Different test_type")
2166 3 : ABI_CHECK_IEQ(hscr_file(ifile)%tordering, hscr_file(1)%tordering, "Different tordering")
2167 6 : ABI_CHECK_IEQ(hscr_file(ifile)%gwgamma, hscr_file(1)%gwgamma, "Different gwgamma")
2168 : end do ! ifile
2169 :
2170 : ! nfreq_tot here is the total *possible* number of freqs.
2171 3 : nfreq_tot = 0
2172 9 : do ifile=1,nfiles
2173 9 : nfreq_tot = nfreq_tot + Hscr_file(ifile)%nomega
2174 : end do
2175 :
2176 : ! freq_indx: mapping (global frequency index, file index) to local file index.
2177 : ! ifile_indx: mapping global frequency index to file index.
2178 : !
2179 9 : ABI_MALLOC(omega_storage, (nfreq_tot))
2180 9 : ABI_MALLOC(omega_wgs_storage, (nfreq_tot))
2181 12 : ABI_MALLOC(freq_indx, (nfreq_tot, nfiles))
2182 9 : ABI_MALLOC(ifile_indx, (nfreq_tot))
2183 321 : omega_storage = CMPLX(-one,-one); freq_indx = 0; ifile_indx = 0
2184 :
2185 : ! TODO: One should handle the case of complex frequencies as well.
2186 : ! Calculate the total number of real freqs and store
2187 3 : nfreqre = 0
2188 9 : do ifile=1,nfiles
2189 87 : do ifrq=1,Hscr_file(ifile)%nomega
2190 78 : skip = .FALSE.
2191 : ! Check whether to skip this point
2192 78 : if (AIMAG(Hscr_file(ifile)%omega(ifrq)) > tol16) skip = .TRUE.
2193 78 : if (REAL(Hscr_file(ifile)%omega(ifrq)) > freqremax) skip = .TRUE.
2194 : ! Check for repetition or non-monotonic points
2195 78 : if (nfreqre > 1) then
2196 832 : do ii=1,nfreqre
2197 832 : if (ABS(REAL(Hscr_file(ifile)%omega(ifrq)) - REAL(omega_storage(ii))) < tol6) skip = .TRUE.
2198 : end do
2199 : end if
2200 78 : if (skip) CYCLE
2201 :
2202 48 : nfreqre = nfreqre + 1
2203 :
2204 : ! Store (complex) frequency and index
2205 48 : omega_storage(nfreqre) = Hscr_file(ifile)%omega(ifrq)
2206 48 : omega_wgs_storage(nfreqre) = Hscr_file(ifile)%omega_wgs(ifrq)
2207 :
2208 48 : ifile_indx(nfreqre) = ifile
2209 48 : freq_indx(nfreqre, ifile) = ifrq
2210 : write(std_out,'(a,es16.6,a,i0,2(a,i0))')&
2211 48 : ' Found real frequency: ',REAL(omega_storage(nfreqre))*Ha_eV,' [eV], number: ',nfreqre,&
2212 132 : ', in file: ',ifile,' local index: ',ifrq
2213 : end do
2214 : end do
2215 :
2216 3 : if (nfreqre > 0) then
2217 : ! Sort real frequencies and rearrange weights as well.
2218 9 : ABI_MALLOC(real_omega, (nfreqre))
2219 6 : ABI_MALLOC(real_omega_wgs, (nfreqre))
2220 9 : ABI_MALLOC(pos_indx, (nfreqre))
2221 6 : ABI_MALLOC(i_temp, (nfreqre))
2222 12 : ABI_MALLOC(i2_temp, (nfreqre,nfiles))
2223 :
2224 : ! Copy real frequencies to temp. sorting array.
2225 51 : real_omega(1:nfreqre) = REAL(omega_storage(1:nfreqre))
2226 :
2227 : ! Sort frequencies while keeping track of index
2228 : ! Set up indexing array
2229 150 : pos_indx = [(ii, ii=1,nfreqre)]
2230 3 : call sort_dp(nfreqre, real_omega, pos_indx, tol16)
2231 :
2232 51 : i_temp(1:nfreqre) = ifile_indx(1:nfreqre)
2233 105 : i2_temp(1:nfreqre,1:nfiles) = freq_indx(1:nfreqre,1:nfiles)
2234 :
2235 : ! Copy sorted frequencies plus file and frequency index
2236 51 : do ii=1,nfreqre
2237 48 : omega_storage(ii) = CMPLX(real_omega(ii),zero)
2238 48 : real_omega_wgs(ii) = omega_wgs_storage(pos_indx(ii))
2239 48 : ifile_indx(ii) = i_temp(pos_indx(ii))
2240 147 : freq_indx(ii,1:nfiles) = i2_temp(pos_indx(ii),1:nfiles)
2241 : end do
2242 51 : omega_wgs_storage(1:nfreqre) = real_omega_wgs(1:nfreqre)
2243 :
2244 3 : ABI_FREE(real_omega)
2245 3 : ABI_FREE(real_omega_wgs)
2246 3 : ABI_FREE(pos_indx)
2247 3 : ABI_FREE(i_temp)
2248 3 : ABI_FREE(i2_temp)
2249 : end if
2250 :
2251 : ! Check imaginary frequencies and store them
2252 3 : nfreqim = 0
2253 9 : do ifile=1,nfiles
2254 87 : do ifrq=1,Hscr_file(ifile)%nomega
2255 78 : if (REAL(Hscr_file(ifile)%omega(ifrq)) > tol8) CYCLE
2256 33 : if (AIMAG(Hscr_file(ifile)%omega(ifrq)) < tol8) CYCLE
2257 30 : nfreqim = nfreqim + 1
2258 30 : omega_storage(nfreqre+nfreqim) = Hscr_file(ifile)%omega(ifrq)
2259 30 : omega_wgs_storage(nfreqre+nfreqim) = Hscr_file(ifile)%omega_wgs(ifrq)
2260 30 : ifile_indx(nfreqre+nfreqim) = ifile
2261 30 : freq_indx(nfreqre+nfreqim,ifile) = ifrq
2262 : write(std_out,'(a,es16.6,a,i0,2(a,i0))')&
2263 30 : ' Found imag frequency: ',AIMAG(omega_storage(nfreqre+nfreqim))*Ha_eV,' [eV], number: ',nfreqim,&
2264 114 : ', in file: ',ifile,' local index: ',ifrq
2265 : end do
2266 : end do
2267 :
2268 : ! Sort imaginary frequencies
2269 9 : ABI_MALLOC(imag_omega, (nfreqim))
2270 6 : ABI_MALLOC(imag_omega_wgs, (nfreqim))
2271 9 : ABI_MALLOC(pos_indx, (nfreqim))
2272 6 : ABI_MALLOC(i_temp, (nfreqim))
2273 12 : ABI_MALLOC(i2_temp, (nfreqim, nfiles))
2274 :
2275 : ! Copy imaginary frequencies to temp. sorting array
2276 33 : imag_omega(1:nfreqim) = AIMAG(omega_storage(nfreqre+1:nfreqre+nfreqim))
2277 :
2278 : ! Sort frequencies while keeping track of index (+1 for ascending order, -1 for descending order)
2279 3 : order = +1
2280 3 : if (mod10 == SIG_GW_AC .and. Hscr_file(1)%iw_mesh_type == "gauss_legendre") order = -1
2281 :
2282 : ! Set up indexing array
2283 96 : pos_indx = [(ii, ii=1,nfreqim)]
2284 3 : call sort_dp(nfreqim, imag_omega, pos_indx, tol16, order=order)
2285 :
2286 33 : i_temp(1:nfreqim) = ifile_indx(nfreqre+1:nfreqre+nfreqim)
2287 69 : i2_temp(1:nfreqim,1:nfiles) = freq_indx(nfreqre+1:nfreqre+nfreqim,1:nfiles)
2288 :
2289 : ! Copy sorted frequencies plus file and frequency index
2290 33 : do ii=1,nfreqim
2291 30 : omega_storage(nfreqre+ii) = CMPLX(zero,imag_omega(ii))
2292 30 : ifile = i_temp(pos_indx(ii))
2293 :
2294 30 : ifile_indx(nfreqre+ii) = ifile
2295 90 : freq_indx(nfreqre+ii,1:nfiles) = i2_temp(pos_indx(ii),1:nfiles)
2296 30 : ifrq = i2_temp(pos_indx(ii), ifile)
2297 33 : imag_omega_wgs(ii) = Hscr_file(ifile)%omega_wgs(ifrq) ! TODO
2298 : end do
2299 :
2300 33 : omega_wgs_storage(nfreqre+1:nfreqre+nfreqim) = imag_omega_wgs
2301 :
2302 3 : ABI_FREE(imag_omega)
2303 3 : ABI_FREE(imag_omega_wgs)
2304 3 : ABI_FREE(pos_indx)
2305 3 : ABI_FREE(i_temp)
2306 3 : ABI_FREE(i2_temp)
2307 :
2308 : ! Here nfreq_tot becomes the *true* number of freq
2309 3 : nfreq_tot = nfreqre + nfreqim
2310 3 : write(std_out,'(2a,i0,a)') ch10,' Merging ',nfreq_tot,' frequencies.'
2311 3 : write(std_out,'(2(a,i0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10
2312 :
2313 : ! Copy old header
2314 3 : call Hscr_file(1)%copy(ohscr)
2315 :
2316 : ! TODO: hscr_wmerge
2317 : ! Then modify entries for new frequency grid.
2318 3 : ohscr%nomega = nfreq_tot
2319 :
2320 9 : ABI_REMALLOC(ohscr%omega, (nfreq_tot))
2321 9 : ABI_REMALLOC(ohscr%omega_wgs, (nfreq_tot))
2322 84 : ohscr%omega = omega_storage(1:nfreq_tot)
2323 84 : ohscr%omega_wgs = omega_wgs_storage
2324 :
2325 3 : ohscr%rw_mesh_type = rw_mesh_type
2326 3 : ohscr%iw_mesh_type = iw_mesh_type
2327 3 : ohscr%cw_mesh_type = cw_mesh_type
2328 :
2329 3 : npwe4mI = ohscr%npwe*ohscr%nI
2330 3 : npwe4mJ = ohscr%npwe*ohscr%nJ
2331 :
2332 : ! Print new header for info
2333 6 : call ohscr%print([std_out], 1, header='Header of the final file')
2334 :
2335 3 : if (endswith(fname_out, ".nc")) then
2336 3 : iomode = IO_MODE_ETSF
2337 3 : NCF_CHECK(nctk_open_create(ount, fname_out, comm))
2338 : else
2339 0 : iomode = IO_MODE_FORTRAN
2340 0 : if (open_file(fname_out, msg, newunit=ount, status='new',form='unformatted') /= 0) then
2341 0 : ABI_ERROR(msg)
2342 : end if
2343 : end if
2344 :
2345 : ! Write the header.
2346 3 : fform_merge = ohscr%fform
2347 3 : abifile = abifile_from_fform(fform_merge)
2348 3 : if (abifile%fform == 0) then
2349 0 : ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
2350 : end if
2351 3 : varname = abifile%varname
2352 :
2353 3 : call ohscr%io(fform_merge, rdwr2, ount, comm, master, iomode)
2354 :
2355 3 : npwe4mI = ohscr%npwe*ohscr%nI
2356 3 : npwe4mJ = ohscr%npwe*ohscr%nJ
2357 3 : nomega4m = ohscr%nomega
2358 :
2359 : ! FIXME: Lot of memory allocated here. One can read and write inside the ifile loop
2360 : ! provided that partial frequencies are always contiguous and files are in the proper order
2361 15 : ABI_MALLOC_OR_DIE(epsm1, (npwe4mI, npwe4mJ, nomega4m, 1), ierr)
2362 :
2363 23 : do iq_ibz=1,ohscr%nqibz
2364 60 : do ifile=1,nfiles
2365 : ! allocate temporary array
2366 40 : npwe4mI = Hscr_file(ifile)%npwe*Hscr_file(ifile)%nI
2367 40 : npwe4mJ = Hscr_file(ifile)%npwe*Hscr_file(ifile)%nJ
2368 40 : nomega4m = Hscr_file(ifile)%nomega
2369 200 : ABI_MALLOC_OR_DIE(epsm1_temp, (npwe4mI, npwe4mJ, nomega4m, 1), ierr)
2370 :
2371 : ! read screening
2372 40 : call read_screening(varname,filenames(ifile),npwe4mI,1,nomega4m,epsm1_temp,iomode,comm,iqiA=iq_ibz)
2373 :
2374 : ! Copy matrices for relevant frequencies
2375 1124 : do ifrq=1,nfreq_tot
2376 1124 : if (ifile_indx(ifrq)==ifile) then
2377 1095158 : epsm1(:,:,ifrq,1)=epsm1_temp(:,:,freq_indx(ifrq,ifile),1)
2378 : end if
2379 : end do
2380 :
2381 60 : ABI_FREE(epsm1_temp)
2382 : end do !ifile
2383 :
2384 : ! Write data.
2385 20 : npwe4mI = ohscr%npwe*ohscr%nI
2386 20 : nomega4m = ohscr%nomega
2387 23 : call write_screening(varname,ount,iomode,npwe4mI,nomega4m,iq_ibz,epsm1)
2388 : end do ! iq_ibz
2389 :
2390 3 : ABI_FREE(epsm1)
2391 3 : ABI_FREE(omega_wgs_storage)
2392 3 : ABI_FREE(omega_storage)
2393 3 : ABI_FREE(freq_indx)
2394 3 : ABI_FREE(ifile_indx)
2395 :
2396 3 : if (iomode == IO_MODE_FORTRAN) then
2397 0 : close(ount)
2398 : else
2399 3 : NCF_CHECK(nf90_close(ount))
2400 : end if
2401 :
2402 3 : write(msg,'(3a)')ch10,' ==== Files have been merged successfully === ',ch10
2403 3 : call wrtout(std_out, msg)
2404 :
2405 3 : end subroutine ioscr_wmerge
2406 : !!***
2407 :
2408 : !----------------------------------------------------------------------
2409 :
2410 : !!****f* m_io_screening/ioscr_wremove
2411 : !! NAME
2412 : !! ioscr_wremove
2413 : !!
2414 : !! FUNCTION
2415 : !! Produce new file by removing selected frequencies in the initial file `inpath`.
2416 : !! This routine should be called by a single MPI process.
2417 : !!
2418 : !! INPUTS
2419 : !! inpath=Input file
2420 : !! ihscr<hscr_t>=Headerf of the input file.
2421 : !! fname_out=Output file.
2422 : !! nfreq_tot=Number of frequencies in new file.
2423 : !! freq_indx(nfreq_tot)=Index of frequency to be kept in input file.
2424 : !!
2425 : !! OUTPUT
2426 : !! ohscr<hscr_t>=The header of the output file.
2427 : !!
2428 : !! SOURCE
2429 :
2430 5 : subroutine ioscr_wremove(inpath, ihscr, fname_out, nfreq_tot, freq_indx, ohscr)
2431 :
2432 : !Arguments ------------------------------------
2433 : !scalars
2434 : integer,intent(in) :: nfreq_tot
2435 : character(len=*),intent(in) :: inpath,fname_out
2436 : type(hscr_t),intent(in) :: ihscr
2437 : type(hscr_t),intent(out) :: ohscr
2438 : !arrays
2439 : integer,intent(in) :: freq_indx(nfreq_tot)
2440 :
2441 : !Local variables-------------------------------
2442 : !scalars
2443 : integer,parameter :: rdwr2=2,master=0
2444 : integer :: iq_ibz,fform_merge,comm,nomega4m,ierr
2445 : integer :: ifrq,npwe4mI,npwe4mJ,iomode,ount
2446 : character(len=500) :: msg
2447 : character(len=nctk_slen) :: varname
2448 : type(abifile_t) :: abifile
2449 : !arrays
2450 1 : complex(gwp),allocatable :: epsm1(:,:,:),epsm1_temp(:,:,:)
2451 : ! *************************************************************************
2452 :
2453 1 : comm = xmpi_comm_self
2454 :
2455 : ! check ifreq_idx
2456 1 : ABI_CHECK(nfreq_tot > 0, "nfreq_tot <= 0!")
2457 1 : if (all(freq_indx == 0)) ABI_ERROR("all(freq_indx == 0)")
2458 :
2459 : ! Copy the old header
2460 1 : call ihscr%copy(ohscr)
2461 :
2462 : ! Then modify entries for new frequency grid.
2463 1 : ohscr%nomega = nfreq_tot
2464 3 : ABI_REMALLOC(ohscr%omega, (nfreq_tot))
2465 3 : ABI_REMALLOC(ohscr%omega_wgs, (nfreq_tot))
2466 17 : do ifrq=1,nfreq_tot
2467 16 : ohscr%omega(ifrq) = ihscr%omega(freq_indx(ifrq))
2468 17 : ohscr%omega_wgs(ifrq) = ihscr%omega_wgs(freq_indx(ifrq))
2469 : end do
2470 :
2471 1 : npwe4mI = ohscr%npwe*ohscr%nI
2472 1 : npwe4mJ = ohscr%npwe*ohscr%nJ
2473 :
2474 : ! Print new header for info
2475 2 : call ohscr%print([std_out], 1, header='Header of the final file')
2476 :
2477 : ! Open output file.
2478 1 : if (endswith(fname_out, ".nc")) then
2479 1 : iomode = IO_MODE_ETSF
2480 1 : NCF_CHECK(nctk_open_create(ount, fname_out, comm))
2481 : else
2482 0 : iomode = IO_MODE_FORTRAN
2483 0 : if (open_file(fname_out, msg, newunit=ount, status='new', form='unformatted') /= 0) then
2484 0 : ABI_ERROR(msg)
2485 : end if
2486 : end if
2487 :
2488 : ! Write the header.
2489 1 : fform_merge = ohscr%fform
2490 1 : abifile = abifile_from_fform(fform_merge)
2491 1 : if (abifile%fform == 0) then
2492 0 : ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform_merge)))
2493 : end if
2494 1 : varname = abifile%varname
2495 :
2496 1 : call ohscr%io(fform_merge,rdwr2,ount,comm,master,iomode)
2497 :
2498 1 : npwe4mI = ohscr%npwe*ohscr%nI; npwe4mJ = ohscr%npwe*ohscr%nJ
2499 1 : nomega4m = ohscr%nomega
2500 :
2501 5 : ABI_MALLOC_OR_DIE(epsm1, (npwe4mI,npwe4mJ,nomega4m), ierr)
2502 :
2503 9 : do iq_ibz=1,ohscr%nqibz
2504 : ! allocate temporary array
2505 8 : npwe4mI = ihscr%npwe * ihscr%nI
2506 8 : npwe4mJ = ihscr%npwe * ihscr%nJ
2507 8 : nomega4m = ihscr%nomega
2508 40 : ABI_MALLOC_OR_DIE(epsm1_temp,(npwe4mI,npwe4mJ,nomega4m), ierr)
2509 :
2510 : ! read full screening matrix for this q-point
2511 8 : call read_screening(varname,inpath,npwe4mI,1,nomega4m,epsm1_temp,iomode,comm,iqiA=iq_ibz)
2512 :
2513 : ! Copy relevant frequencies
2514 136 : do ifrq=1,nfreq_tot
2515 96904 : epsm1(:,:,ifrq) = epsm1_temp(:,:,freq_indx(ifrq))
2516 : end do
2517 :
2518 8 : ABI_FREE(epsm1_temp)
2519 :
2520 8 : npwe4mI = ohscr%npwe*ohscr%nI; nomega4m = ohscr%nomega
2521 9 : call write_screening(varname,ount,iomode,npwe4mI,nomega4m,iq_ibz,epsm1)
2522 : end do ! iq_ibz
2523 :
2524 1 : ABI_FREE(epsm1)
2525 :
2526 1 : if (iomode == IO_MODE_FORTRAN) then
2527 0 : close(ount)
2528 : else
2529 1 : NCF_CHECK(nf90_close(ount))
2530 : end if
2531 :
2532 1 : write(msg,'(3a)')ch10,' ==== Frequencies have been removed successfully === ',ch10
2533 1 : call wrtout(std_out, msg)
2534 :
2535 1 : end subroutine ioscr_wremove
2536 : !!***
2537 :
2538 : !!****f* m_io_screening/get_hscr_qmesh_gsph
2539 : !! NAME
2540 : !! get_hscr_qmesh_gsph
2541 : !!
2542 : !! FUNCTION
2543 : !!
2544 : !! INPUTS
2545 : !!
2546 : !! OUTPUT
2547 : !!
2548 : !! SOURCE
2549 :
2550 306 : subroutine get_hscr_qmesh_gsph(w_fname, dtset, cryst, hscr, qmesh, gsph_c, qlwl, comm)
2551 :
2552 : !Arguments ------------------------------------
2553 : !scalars
2554 : character(len=*),intent(inout) :: w_fname
2555 : type(dataset_type),intent(inout) :: dtset
2556 : type(crystal_t),intent(in) :: cryst
2557 : type(hscr_t),intent(out) :: hscr
2558 : type(kmesh_t),intent(out) :: qmesh
2559 : type(gsphere_t),intent(out) :: gsph_c
2560 : real(dp),allocatable,intent(out) :: qlwl(:,:)
2561 : integer,intent(in) :: comm
2562 :
2563 : !Local variables-------------------------------
2564 : integer,parameter :: master = 0
2565 : integer :: my_rank, fform, npwe_file, nqlwl, ierr
2566 : character(len=500) :: msg
2567 : ! *************************************************************************
2568 :
2569 17 : my_rank = xmpi_comm_rank(comm)
2570 :
2571 17 : if (my_rank == master) then
2572 : ! Read dimensions from the external file.
2573 17 : if (.not. file_exists(w_fname)) then
2574 13 : w_fname = nctk_ncify(w_fname)
2575 13 : ABI_COMMENT(sjoin("File not found. Will try netcdf file: ", w_fname))
2576 : end if
2577 : ! Master reads npw and nqlwl from the SCR file.
2578 17 : call wrtout(std_out, sjoin('Testing SCR file: ', w_fname))
2579 17 : call hscr%from_file(w_fname, fform, xmpi_comm_self)
2580 :
2581 : ! Have to change %npweps if it was larger than dim on disk.
2582 17 : npwe_file = hscr%npwe
2583 17 : nqlwl = hscr%nqlwl
2584 :
2585 17 : if (dtset%npweps > npwe_file) then
2586 : write(msg,'(2(a,i0),2a,i0)')&
2587 0 : "The number of G-vectors stored on file (",npwe_file,") is smaller than input dtset%npweps: ",dtset%npweps,ch10,&
2588 0 : "Calculation will proceed with the maximum available set, npwe_file: ",npwe_file
2589 0 : ABI_WARNING(msg)
2590 0 : dtset%npweps = npwe_file
2591 17 : else if (dtset%npweps < npwe_file .and. dtset%npweps /= 0) then
2592 : write(msg,'(2(a,i0),2a,i0)')&
2593 8 : "The number of G-vectors stored on file (",npwe_file,") is larger than dtset%npweps: ",dtset%npweps,ch10,&
2594 16 : "Calculation will proceed with dtset%npweps: ",dtset%npweps
2595 8 : ABI_COMMENT(msg)
2596 : else
2597 9 : call Gsph_c%init(cryst, 0, ecut=dtset%ecuteps)
2598 9 : if (Gsph_c%ng > npwe_file) then
2599 0 : dtset%npweps = npwe_file
2600 : write(msg,'(2a,f4.1,a,i0,2a,i0)')&
2601 0 : "npweps was not set in input.",&
2602 0 : ch10//"The number of G-vectors generated according to ecuteps (",dtset%ecuteps,") is larger than that stored on file (",npwe_file,")",&
2603 0 : ch10//"Calculation will proceed with the maximum available set: ",npwe_file
2604 0 : ABI_COMMENT(msg)
2605 : else
2606 9 : dtset%npweps = Gsph_c%ng
2607 : write(msg,'(2a,f4.1,a,i0,2a,f4.1)')&
2608 9 : "npweps was not set in input.",&
2609 9 : ch10//"The number of G-vectors generated according to ecuteps (",dtset%ecuteps,") is smaller than that stored on file (",npwe_file,")",&
2610 18 : ch10//"Calculation will proceed with ecuteps: ",dtset%ecuteps
2611 9 : ABI_COMMENT(msg)
2612 : end if
2613 9 : call Gsph_c%free()
2614 : end if
2615 : end if
2616 :
2617 17 : call xmpi_bcast(w_fname, master, comm, ierr)
2618 17 : call hscr%bcast(master, my_rank, comm)
2619 17 : call xmpi_bcast(dtset%npweps, master, comm, ierr)
2620 17 : call xmpi_bcast(nqlwl, master, comm, ierr)
2621 :
2622 17 : if (nqlwl > 0) then
2623 51 : ABI_MALLOC(qlwl, (3, nqlwl))
2624 102 : qlwl = hscr%qlwl
2625 : end if
2626 :
2627 : ! Init qmesh from the SCR file.
2628 17 : call qmesh%init(cryst, hscr%nqibz, hscr%qibz, dtset%kptopt)
2629 :
2630 17 : call gsph_c%init(cryst, dtset%npweps, gvec=hscr%gvec)
2631 :
2632 17 : end subroutine get_hscr_qmesh_gsph
2633 : !!***
2634 :
2635 636 : end module m_io_screening
2636 : !!***
|
