Branch data 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 : 512 : character(len=nctk_slen) function ncname_from_id(id) result(varname)
275 : :
276 : : integer,intent(in) :: id
277 : :
278 : 512 : varname = "None"
279 [ + + ]: 512 : if (id == 1) varname = chi0_ncname
280 [ - + ]: 512 : if (id == 3) varname = e_ncname
281 [ + + ]: 512 : if (id == 4) varname = em1_ncname
282 [ - + ]: 512 : ABI_CHECK(varname /= "None", "Wrong id")
283 : :
284 : 512 : 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 [ + - ]: 218 : 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 : 218 : my_rank = xmpi_comm_rank(comm)
321 : :
322 : : ! Master reads and broadcasts.
323 [ + + ]: 218 : if (my_rank == master) then
324 [ - + ]: 190 : 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 [ - + ]: 190 : NCF_CHECK(nctk_open_read(unt, path, xmpi_comm_self))
334 : 190 : call hscr%io(fform, rdwr5, unt, xmpi_comm_self, master, IO_MODE_ETSF)
335 [ - + ]: 190 : NCF_CHECK(nf90_close(unt))
336 : : end if
337 : :
338 [ - + ]: 190 : ABI_CHECK(fform /= 0, sjoin("hscr_io returned fform == 0 while reading:", path))
339 : : end if
340 : :
341 : : ! Broadcast data.
342 [ + + ]: 218 : 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 : 218 : 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 : 915 : 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 : 915 : real(dp),allocatable :: real_omega(:,:)
411 : 915 : 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 : 915 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
420 : :
421 [ + + ]: 915 : 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 [ - + ]: 776 : if (iomode==IO_MODE_FORTRAN) then
428 : 0 : call hscr%hdr%fort_read(unt, fform, rewind=(rdwr==1))
429 [ + - ]: 776 : else if (iomode==IO_MODE_ETSF) then
430 : 776 : call hscr%hdr%ncread(unt, fform)
431 : : end if
432 : :
433 : : ! Reset the variables absent in old versions.
434 : 776 : Hscr%fform=fform
435 : :
436 [ - + ]: 776 : 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 [ + - ]: 776 : else if (iomode == IO_MODE_ETSF) then
476 : 776 : ncid = unt
477 : :
478 : 776 : select case (fform)
479 : : case (1003, 1004)
480 : : ! Get dimensions and allocate arrays.
481 [ - + ]: 776 : NCF_CHECK(nctk_get_dim(ncid, "number_of_coefficients_dielectric_function", hscr%npwe))
482 [ - + ]: 776 : NCF_CHECK(nctk_get_dim(ncid, "number_of_qpoints_dielectric_function", hscr%nqibz))
483 [ - + ]: 776 : NCF_CHECK(nctk_get_dim(ncid, "number_of_frequencies_dielectric_function", hscr%nomega))
484 [ - + ]: 776 : NCF_CHECK(nctk_get_dim(ncid, "number_of_qpoints_gamma_limit", hscr%nqlwl))
485 [ - + ]: 776 : NCF_CHECK(nctk_get_dim(ncid, "nI", hscr%ni))
486 [ - + ]: 776 : NCF_CHECK(nctk_get_dim(ncid, "nJ", hscr%nj))
487 : 776 : call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
488 : :
489 : 776 : varid = nctk_idname(ncid, 'reduced_coordinates_plane_waves_dielectric_function')
490 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, varid, hscr%gvec, start=[1,1,1]))
491 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid('qpoints_dielectric_function'), hscr%qibz))
492 : :
493 [ - + + - : 2328 : ABI_MALLOC(real_omega, (2, hscr%nomega))
+ - + - -
+ - + ]
494 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid('frequencies_dielectric_function'), real_omega))
495 [ + - - + : 5341 : hscr%omega = dcmplx(real_omega(1,:), real_omega(2,:))
- - + + ]
496 [ - + ]: 776 : ABI_FREE(real_omega)
497 : :
498 : : ! Read extra data added in new header.
499 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("vcutgeo"), hscr%vcutgeo))
500 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("id"), hscr%id))
501 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("ikxc"), hscr%ikxc))
502 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("inclvkb"), hscr%inclvkb))
503 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("headform"), hscr%headform))
504 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("fform"), hscr%fform))
505 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("gwcalctyp"), hscr%gwcalctyp))
506 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("nbands_used"), hscr%nbnds_used))
507 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("npwwfn_used"), hscr%npwwfn_used))
508 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("spmeth"), hscr%spmeth))
509 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("test_type"), hscr%test_type))
510 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("tordering"), hscr%tordering))
511 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("awtr"), hscr%awtr))
512 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("gw_icutcoul"), hscr%icutcoul))
513 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("gwcomp"), hscr%gwcomp))
514 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("gwgamma"), hscr%gwgamma))
515 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("mbpt_sciss"), hscr%mbpt_sciss))
516 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("spsmear"), hscr%spsmear))
517 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("zcut"), hscr%zcut))
518 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("gwencomp"), hscr%gwencomp))
519 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("kind_cdata"), hscr%kind_cdata))
520 : 776 : call replace_ch0(hscr%kind_cdata)
521 : :
522 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("titles"), hscr%titles))
523 [ + + ]: 2328 : call replace_ch0(hscr%titles(:))
524 : :
525 : : ! TODO Read it
526 [ + - ]: 776 : if (hscr%nqlwl /= 0) then
527 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("qpoints_gamma_limit"), hscr%qlwl))
528 : : end if
529 : :
530 : 776 : call set_defaults_from_gwcalctyp()
531 [ + - ]: 776 : if (hscr%headform > 80) then
532 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("rw_mesh_type"), hscr%rw_mesh_type))
533 : 776 : call replace_ch0(hscr%rw_mesh_type)
534 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("iw_mesh_type"), hscr%iw_mesh_type))
535 : 776 : call replace_ch0(hscr%iw_mesh_type)
536 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("cw_mesh_type"), hscr%cw_mesh_type))
537 : 776 : call replace_ch0(hscr%cw_mesh_type)
538 [ - + ]: 776 : NCF_CHECK(nf90_get_var(ncid, vid("omega_wgs"), hscr%omega_wgs))
539 : : end if
540 : :
541 : : case default
542 [ + - ]: 776 : 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 [ + - ]: 139 : else if (rdwr == 2 .or. rdwr == 6) then
554 : : ! Writing the header of an unformatted file.
555 : : ! Always use the latest version.
556 : :
557 [ - + ]: 139 : 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 [ + - ]: 139 : else if (iomode == IO_MODE_ETSF) then
593 : 139 : ncid = unt
594 : : ! Write the abinit header, rewinding of the file (if any) is done here.
595 [ - + ]: 139 : 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 [ + + ]: 1390 : nctkdim_t("number_of_coefficients_dielectric_function", hscr%npwe)], defmode=.True.)
611 [ - + ]: 139 : 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 [ + + ]: 1668 : ])
633 [ - + ]: 139 : 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 : 139 : 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 : 139 : &number_of_spins, number_of_spins, number_of_frequencies_dielectric_function, number_of_qpoints_dielectric_function"))
643 [ - + ]: 139 : NCF_CHECK(ncerr)
644 : :
645 : : !write(std_out,*)"nqlwl",hscr%nqlwl
646 [ - + ]: 139 : NCF_CHECK(nctk_set_datamode(ncid))
647 [ + + ]: 417 : call c_f_pointer(c_loc(hscr%omega(1)), r2vals, shape=[2, size(hscr%omega)])
648 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid('frequencies_dielectric_function'), r2vals))
649 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid('qpoints_dielectric_function'), hscr%qibz))
650 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid('reduced_coordinates_plane_waves_dielectric_function'), hscr%gvec))
651 : :
652 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid("titles"), hscr%titles))
653 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid("kind_cdata"), hscr%kind_cdata))
654 [ - + ]: 139 : 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 [ + + + + ]: 4309 : ])
665 [ - + ]: 139 : 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 [ + + + + ]: 1251 : ])
671 [ - + ]: 139 : NCF_CHECK(ncerr)
672 : :
673 : : ! Add q-points for heads and wings for q-->0.
674 [ + - ]: 139 : 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 [ - + ]: 139 : NCF_CHECK(nctk_set_datamode(ncid))
689 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid('qpoints_gamma_limit'), hscr%qlwl))
690 : :
691 : : !if (hscr%headform > 80) then
692 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid("rw_mesh_type"), hscr%rw_mesh_type))
693 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid("iw_mesh_type"), hscr%iw_mesh_type))
694 [ - + ]: 139 : NCF_CHECK(nf90_put_var(ncid, vid("cw_mesh_type"), hscr%cw_mesh_type))
695 [ - + ]: 139 : 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 : 915 : return
710 : :
711 : : ! Handle Fortran IO error
712 : : 10 continue
713 [ + - ]: 1830 : ABI_ERROR(errmsg)
714 : :
715 : : contains
716 : 24033 : integer function vid(vname)
717 : : character(len=*),intent(in) :: vname
718 : 24033 : vid = nctk_idname(ncid, vname)
719 : : end function vid
720 : :
721 : 776 : 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 : 776 : mod10 = mod(hscr%gwcalctyp, 10)
725 : 776 : hscr%iw_mesh_type = "None"; hscr%rw_mesh_type = "None"; hscr%cw_mesh_type = "None"
726 : :
727 [ + + ]: 776 : if (mod10 == SIG_GW_AC) then
728 : 34 : hscr%iw_mesh_type = "gauss_legendre"
729 [ + + ]: 742 : 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 : 776 : 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 : 344 : 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 [ + + ]: 172 : if (PRESENT(header)) then
765 [ + - + - ]: 11 : msg = ' ==== '//TRIM(ADJUSTL(header))//' ==== '
766 : 11 : call wrtout(units, msg)
767 : : end if
768 : :
769 [ + - ]: 172 : write(msg,'(1x,a)')TRIM(hscr%titles(1))
770 : 172 : call wrtout(units,msg)
771 [ + + ]: 172 : write(msg,'(1x,a)')TRIM(hscr%titles(2))
772 : 172 : call wrtout(units,msg)
773 : 172 : write(msg,'(a,i8)') ' Identifier ',hscr%ID
774 : 172 : call wrtout(units,msg)
775 : 172 : write(msg,'(a,i8)') ' Kxc kernel ',hscr%ikxc
776 : 172 : call wrtout(units,msg)
777 : 172 : write(msg,'(a,i8)') ' Treatment of q-->0 limit ',hscr%inclvkb
778 : 172 : call wrtout(units,msg)
779 : 172 : write(msg,'(a,i8)') '- headform ',hscr%headform
780 : 172 : call wrtout(units,msg)
781 : 172 : write(msg,'(a,i8)') ' fform ',hscr%fform
782 : 172 : call wrtout(units,msg)
783 : 172 : write(msg,'(a,i8)') ' gwcalctyp ',hscr%gwcalctyp
784 : 172 : call wrtout(units,msg)
785 : 172 : write(msg,'(a,2i8)')' Number of components ',hscr%nI,hscr%nJ
786 : 172 : call wrtout(units,msg)
787 : 172 : write(msg,'(a,i8)') ' Number of q-points ',hscr%nqibz
788 : 172 : call wrtout(units,msg)
789 : 172 : write(msg,'(a,i8)') ' Number of q-directions ',hscr%nqlwl
790 : 172 : call wrtout(units,msg)
791 : 172 : write(msg,'(a,i8)') ' Number of frequencies ',hscr%nomega
792 : 172 : call wrtout(units,msg)
793 : 172 : write(msg,'(a,i8)') ' Number of bands used ',hscr%nbnds_used
794 : 172 : call wrtout(units,msg)
795 : 172 : write(msg,'(a,i8)') ' Dimension of matrix ',hscr%npwe
796 : 172 : call wrtout(units,msg)
797 : 172 : write(msg,'(a,i8)') ' Number of planewaves used ',hscr%npwwfn_used
798 : 172 : call wrtout(units,msg)
799 : 172 : write(msg,'(a,i8)') ' Spectral method ',hscr%spmeth
800 : 172 : call wrtout(units,msg)
801 : 172 : write(msg,'(a,i8)') ' Test_type ',hscr%test_type
802 : 172 : call wrtout(units,msg)
803 : 172 : write(msg,'(a,i8)') ' Time-ordering ',hscr%tordering
804 : 172 : call wrtout(units,msg)
805 : 172 : write(msg,'(a,es16.6)')' Scissor Energy ',hscr%mbpt_sciss
806 : 172 : call wrtout(units,msg)
807 : 172 : write(msg,'(a,es16.6)')' Spectral smearing ',hscr%spsmear
808 : 172 : call wrtout(units,msg)
809 : 172 : write(msg,'(a,es16.6)')' Complex Imaginary Shift ',hscr%zcut
810 : 172 : call wrtout(units,msg)
811 : 172 : write(msg,'(2a)') ' rw_mesh_type ',trim(hscr%rw_mesh_type)
812 : 172 : call wrtout(units,msg)
813 : 172 : write(msg,'(2a)') ' iw_mesh_type ',trim(hscr%iw_mesh_type)
814 : 172 : call wrtout(units,msg)
815 : 172 : write(msg,'(2a)') ' cw_mesh_type ',trim(hscr%cw_mesh_type)
816 : 172 : call wrtout(units,msg)
817 : :
818 [ + + ]: 172 : if (prtvol == 0) then
819 : 148 : 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 [ - + ]: 172 : 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 : 260 : 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 : 130 : id = 0
888 [ + + ]: 130 : if (varname == "polarizability") id = 1
889 [ + + ]: 130 : if (varname == "inverse_dielectric_function") id = 4
890 [ - + ]: 8 : ABI_CHECK(id /= 0, sjoin("Invalid varname: ",varname))
891 : :
892 : : ! Get fform from abifile.
893 : 130 : abifile = abifile_from_varname(varname)
894 [ - + ]: 130 : 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 : 130 : call hdr_abinit%copy(Hscr%Hdr)
900 : :
901 : : ! Initialize quantities related to the screening file
902 : 130 : hscr%id =id
903 : 130 : hscr%ikxc =ikxc
904 : 130 : hscr%inclvkb =Ep%inclvkb
905 : 130 : hscr%headform =HSCR_LATEST_HEADFORM
906 : 130 : hscr%fform =abifile%fform
907 : 130 : hscr%gwcalctyp =Ep%gwcalctyp
908 : 130 : hscr%nI =Ep%nI
909 : 130 : hscr%nJ =Ep%nJ
910 : 130 : hscr%nqibz =Ep%nqcalc ! nqcalc == nqibz except if we split the calculation with nqptdm
911 : 130 : hscr%nqlwl =Ep%nqlwl
912 : 130 : hscr%nomega =Ep%nomega
913 : 130 : hscr%nbnds_used =Ep%nbnds
914 : 130 : hscr%npwe =Ep%npwe
915 : 130 : hscr%npwwfn_used=Ep%npwwfn
916 : 130 : hscr%spmeth =Ep%spmeth
917 : 130 : hscr%test_type =test_type
918 : 130 : hscr%tordering =tordering
919 : 130 : hscr%mbpt_sciss =Ep%mbpt_sciss
920 : 130 : hscr%spsmear =Ep%spsmear
921 : 130 : hscr%zcut =Ep%zcut
922 : :
923 [ + + ]: 390 : hscr%titles(:)=titles(:)
924 : :
925 : 130 : call hscr_malloc(hscr, hscr%npwe, hscr%nqibz, hscr%nomega, hscr%nqlwl)
926 [ + + + + ]: 29978 : hscr%gvec(:,:) = gvec(:,1:Ep%npwe)
927 [ + + + + ]: 3194 : hscr%qibz(:,:) = Ep%qcalc
928 [ + + + + ]: 650 : hscr%qlwl(:,:) = Ep%qlwl
929 [ + + ]: 1001 : hscr%omega(:) = Ep%omega
930 : :
931 [ + + ]: 1001 : hscr%omega_wgs(:) = Ep%omega_wgs
932 : 130 : hscr%iw_mesh_type = Ep%iw_mesh_type
933 : 130 : hscr%rw_mesh_type = Ep%rw_mesh_type
934 : 130 : hscr%cw_mesh_type = Ep%cw_mesh_type
935 : :
936 : : ! HSCR_NEW
937 : 130 : hscr%awtr = dtset%awtr
938 : 130 : hscr%icutcoul = dtset%gw_icutcoul
939 [ + + ]: 520 : hscr%vcutgeo = dtset%vcutgeo
940 : 130 : hscr%gwcomp = dtset%gwcomp
941 : 130 : hscr%gwgamma = dtset%gwgamma
942 : 130 : hscr%gwencomp = dtset%gwencomp
943 : 130 : hscr%kind_cdata = "dpc" ! For the time being, data is always written in double-precision
944 : : ! HSCR_NEW
945 : :
946 [ + + - + ]: 520 : 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 : 934 : 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 [ - + + - : 2802 : ABI_MALLOC(hscr%gvec, (3, npwe))
+ - + - -
+ - + -
+ ]
1066 [ - + + - : 2802 : ABI_MALLOC(hscr%qibz, (3, nqibz))
+ - + - -
+ - + -
+ ]
1067 [ - + + - : 2802 : ABI_MALLOC(hscr%qlwl, (3, nqlwl))
+ - + - -
+ - + -
+ ]
1068 [ - + + - : 2802 : ABI_MALLOC(hscr%omega, (nomega))
+ - + - -
+ - + -
+ ]
1069 [ - + + - : 7558 : ABI_CALLOC(hscr%omega_wgs, (nomega))
+ - + - -
+ - + - +
+ + ]
1070 : :
1071 : 934 : 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 : 979 : subroutine hscr_free(hscr)
1086 : :
1087 : : !Arguments ------------------------------------
1088 : : class(hscr_t),intent(inout) :: hscr
1089 : : ! *************************************************************************
1090 : :
1091 [ + + ]: 979 : ABI_SFREE(hscr%gvec)
1092 [ + + ]: 979 : ABI_SFREE(hscr%qibz)
1093 [ + + ]: 979 : ABI_SFREE(hscr%qlwl)
1094 [ + + ]: 979 : ABI_SFREE(hscr%omega)
1095 [ + + ]: 979 : ABI_SFREE(hscr%omega_wgs)
1096 : :
1097 : 979 : call hscr%Hdr%free()
1098 : :
1099 : 979 : 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 : 821 : 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 : 821 : complex(dp),allocatable :: epsm1d(:,:)
1363 : : integer :: varid,ncerr
1364 : : #ifdef HAVE_GW_DPC
1365 : 821 : 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 : 821 : 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 : 821 : varid = nctk_idname(unt, varname)
1392 [ + + ]: 6568 : 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 [ + + ]: 2463 : spins = 1; s1 = spins(1); s2 = spins(2)
1395 [ + + + + ]: 12315 : 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 [ - + ]: 821 : 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 [ - + - - : 821 : 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 : 586 : 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 [ + + ]: 2344 : type(hscr_t) :: Hscr
1472 : : !arrays
1473 : : #ifdef HAVE_MPI_IO
1474 : 586 : integer(MPI_OFFSET_KIND),allocatable :: offset_wq(:,:)
1475 : : #endif
1476 : 586 : complex(dp),allocatable :: bufdc2d(:,:),bufdc3d(:,:,:)
1477 : : ! pointers passed to netcdf4 routines (complex datatypes are not supported).
1478 : : #ifdef HAVE_GW_DPC
1479 : 586 : 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 : 586 : call cwtime(cpu, wall, gflops, "start")
1488 : :
1489 : 586 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1490 : 586 : my_iomode = iomode
1491 [ + - ]: 586 : 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 : 586 : 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 [ - + ]: 586 : NCF_CHECK(nctk_open_read(unt, fname, xmpi_comm_self))
1537 : 586 : call hscr%io(fform, rdwr, unt, comm, master, my_iomode)
1538 : :
1539 : : case default
1540 [ - - + - ]: 586 : ABI_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode)))
1541 : : end select
1542 : :
1543 : : ! Slice or full array?
1544 : 586 : read_qslice = .False.
1545 [ + + ]: 586 : 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 [ + + ]: 586 : 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 [ - + ]: 586 : 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 [ - + ]: 586 : ABI_CHECK(nqibzA <= Hscr%nqibz, 'Requiring too many q-points')
1567 [ - + ]: 586 : 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 : 586 : varid = nctk_idname(unt, varname)
1697 : :
1698 : : ! [cplex, npwe, npwe, nspin, nspin, nomega, nqpt]
1699 [ + + ]: 4688 : call c_f_pointer(c_loc(epsm1(1,1,1,1)), real_epsm1, [2,npweA,npweA,1,1,nomegaA,nqibzA])
1700 [ + + ]: 1758 : spins = 1; s1 = spins(1); s2 = spins(2)
1701 [ + + ]: 586 : 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 [ + + + + ]: 2025 : 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 [ - + ]: 586 : NCF_CHECK_MSG(ncerr, sjoin("getting var:", varname))
1710 [ - + ]: 586 : NCF_CHECK(nf90_close(unt))
1711 : : !write(std_out,*)"read_screening done"
1712 : :
1713 : : case default
1714 [ - - + - ]: 586 : ABI_ERROR(sjoin("Wrong iomode:", iomode2str(my_iomode)))
1715 : : end select
1716 : :
1717 : : ! Free memory
1718 [ - + ]: 586 : ABI_SFREE(bufdc2d)
1719 : : ABI_SFREE(bufdc3d)
1720 : :
1721 : 586 : call Hscr%free()
1722 : :
1723 : : !call cwtime_report("- read_screening:", cpu, wall, gflops)
1724 : : DBG_EXIT("COLL")
1725 : :
1726 : 586 : return
1727 : :
1728 : : ! Handle Fortran IO error.
1729 : : 10 continue
1730 : 0 : ABI_ERROR(errmsg)
1731 : :
1732 [ - + - + : 1172 : 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 [ - + - - : 654 : end module m_io_screening
+ + - + +
- - + - +
+ - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
- - - - -
- - - - -
+ ]
2636 : : !!***
|
