Line data Source code
1 : !!****m* ABINIT/m_hdr
2 : !! NAME
3 : !! m_hdr
4 : !!
5 : !! FUNCTION
6 : !! This module contains the definition of the abinit header and its methods
7 : !! If you have to change the hdr, pay attention to the following subroutines:
8 : !!
9 : !! hdr_malloc, hdr_init_lowlvl, hdr_free, hdr_bcast and the IO routines
10 : !! hdr_mpio_skip, hdr_fort_read, hdr_fort_write, hdr_ncread, hdr_ncwrite
11 : !!
12 : !! COPYRIGHT
13 : !! Copyright (C) 2008-2026 ABINIT group (XG, MB, MT, DC, MG)
14 : !! This file is distributed under the terms of the
15 : !! GNU General Public License, see ~abinit/COPYING
16 : !! or http://www.gnu.org/copyleft/gpl.txt .
17 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
18 : !!
19 : !! SOURCE
20 :
21 : #if defined HAVE_CONFIG_H
22 : #include "config.h"
23 : #endif
24 :
25 : #include "abi_common.h"
26 :
27 : !#define DEBUG_MODE
28 :
29 : ! This option enable the output of the new hdr entries in hdr%echo
30 : ! Reference files should be updated
31 : !#define DEV_NEW_HDR
32 :
33 : module m_hdr
34 :
35 : use defs_basis
36 : USE_MPI
37 : use m_xmpi
38 : use m_abicore
39 : use m_errors
40 : use m_crystal
41 : use m_wffile
42 : use m_sort
43 : use netcdf
44 : use m_nctk
45 : use m_dtset
46 :
47 : use m_build_info, only : abinit_version
48 : use m_copy, only : alloc_copy
49 : use m_io_tools, only : flush_unit, isncfile, file_exists, open_file
50 : use m_fstrings, only : sjoin, itoa, ftoa, ltoa, replace_ch0, startswith, endswith, ljust, strcat, atoi
51 : use m_symtk, only : print_symmetries
52 : !use m_kpts, only : kpts_timrev_from_kptopt
53 : use defs_wvltypes, only : wvl_internal_type
54 : use defs_datatypes, only : ebands_base_t, pseudopotential_type
55 : use m_pawtab, only : pawtab_type
56 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, &
57 : pawrhoij_io, pawrhoij_inquire_dim
58 :
59 : implicit none
60 :
61 : private
62 : !!***
63 :
64 : #if defined HAVE_MPI1
65 : include 'mpif.h'
66 : #endif
67 :
68 : !----------------------------------------------------------------------
69 :
70 : !!****t* m_hdr/hdr_type
71 : !! NAME
72 : !! hdr_type
73 : !!
74 : !! FUNCTION
75 : !! It contains all the information needed to write a header for a wf, den or pot file.
76 : !! The structure of the header is explained in the abinit_help.html and other associated html files.
77 : !! The datatype is considered as an object, to which are attached a whole
78 : !! set of "methods", actually, different subroutines.
79 : !! A few of these subroutines are: hdr_init, hdr_update, hdr_free, hdr_check, hdr_io, hdr_skip.
80 : !!
81 : !! SOURCE
82 :
83 : type, public :: hdr_type
84 :
85 : ! WARNING : if you modify this datatype, please check whether there might be creation/destruction/copy routines,
86 : ! declared in another part of ABINIT, that might need to take into account your modification.
87 :
88 : integer :: bantot ! total number of bands (sum of nband on all kpts and spins)
89 : integer :: date ! starting date
90 : integer :: headform ! format of the header
91 : integer :: intxc ! input variable
92 : integer :: ivalence=1 ! occopt=9 variable
93 : integer :: ixc ! input variable
94 : integer :: mband ! maxval(hdr%nband)
95 : integer :: natom ! input variable
96 : integer :: nkpt ! input variable
97 : integer :: npsp ! input variable
98 : integer :: nspden ! input variable
99 : integer :: nspinor ! input variable
100 : integer :: nsppol ! input variable
101 : integer :: nsym ! input variable
102 : integer :: ntypat ! input variable
103 : integer :: occopt ! input variable
104 : integer :: pertcase ! the index of the perturbation, 0 if GS calculation
105 : integer :: usepaw ! input variable (0=norm-conserving psps, 1=paw)
106 : integer :: usewvl ! input variable (0=plane-waves, 1=wavelets)
107 :
108 : integer :: kptopt ! input variable (defines symmetries used for k-point sampling)
109 : integer :: pawcpxocc ! input variable
110 : integer :: nshiftk_orig=1 ! original number of shifts given in input (changed in inkpts, the actual value is nshiftk)
111 : integer :: nshiftk=1 ! number of shifts after inkpts.
112 : integer :: icoulomb ! input variable.
113 :
114 : real(dp) :: ecut ! input variable
115 : real(dp) :: ecutdg ! input variable (ecut for NC psps, pawecutdg for paw)
116 : real(dp) :: ecutsm ! input variable
117 : real(dp) :: ecut_eff ! ecut*dilatmx**2 (dilatmx is an input variable)
118 : real(dp) :: etot ! EVOLVING variable
119 : real(dp) :: fermie ! EVOLVING variable
120 : real(dp) :: fermih=zero ! EVOLVING variable
121 : real(dp) :: residm ! EVOLVING variable
122 : real(dp) :: stmbias ! input variable
123 : real(dp) :: tphysel ! input variable
124 : real(dp) :: tsmear ! input variable
125 : real(dp) :: nelect ! number of electrons (computed from pseudos and cellcharge)
126 : real(dp) :: ne_qFD=zero ! number of excited electrons (input variable)
127 : real(dp) :: nh_qFD=zero ! number of excited holes (input variable)
128 : real(dp) :: cellcharge ! input variable (for the first image if more than one)
129 : real(dp) :: extfpmd_eshift=zero ! Energy shift of the extended for high temperature
130 :
131 : ! This record is not a part of the hdr_type, although it is present in the
132 : ! header of the files. This is because it depends on the kind of file
133 : ! that is written, while all other information does not depend on it.
134 : ! It was preferred to let it be initialized or defined outside of hdr_type.
135 : ! integer :: fform ! file format
136 :
137 : real(dp) :: qptn(3)
138 : ! the wavevector, in case of a perturbation
139 :
140 : real(dp) :: rprimd(3,3)
141 : ! EVOLVING variables
142 :
143 : integer :: ngfft(3)
144 : ! input variable
145 :
146 : integer :: nwvlarr(2)
147 : ! nwvlarr(2) array holding the number of wavelets for each resolution.
148 :
149 : integer :: kptrlatt_orig(3,3)
150 : ! Original kptrlatt
151 :
152 : integer :: kptrlatt(3,3)
153 : ! kptrlatt after inkpts.
154 :
155 : integer, allocatable :: istwfk(:)
156 : ! input variable istwfk(nkpt)
157 :
158 : integer, allocatable :: lmn_size(:)
159 : ! lmn_size(npsp) from psps
160 :
161 : integer, allocatable :: nband(:)
162 : ! input variable nband(nkpt*nsppol)
163 :
164 : integer, allocatable :: npwarr(:)
165 : ! npwarr(nkpt) array holding npw for each k point
166 :
167 : integer, allocatable :: pspcod(:)
168 : ! pscod(npsp) from psps
169 :
170 : integer, allocatable :: pspdat(:)
171 : ! psdat(npsp) from psps
172 :
173 : integer, allocatable :: pspso(:)
174 : ! pspso(npsp) from psps
175 :
176 : integer, allocatable :: pspxc(:)
177 : ! pspxc(npsp) from psps
178 :
179 : integer, allocatable :: so_psp(:)
180 : ! input variable so_psp(npsp)
181 :
182 : integer, allocatable :: symafm(:)
183 : ! input variable symafm(nsym)
184 :
185 : integer, allocatable :: symrel(:,:,:)
186 : ! input variable symrel(3,3,nsym)
187 :
188 : integer, allocatable :: typat(:)
189 : ! input variable typat(natom)
190 :
191 : real(dp), allocatable :: kptns(:,:)
192 : ! input variable kptns(3,nkpt)
193 :
194 : real(dp), allocatable :: occ(:)
195 : ! EVOLVING variable occ(bantot)
196 :
197 : real(dp), allocatable :: tnons(:,:)
198 : ! input variable tnons(3,nsym)
199 :
200 : real(dp), allocatable :: wtk(:)
201 : ! weight of kpoints wtk(nkpt)
202 :
203 : real(dp),allocatable :: shiftk_orig(:,:)
204 : ! original shifts given in input (changed in inkpts).
205 :
206 : real(dp),allocatable :: shiftk(:,:)
207 : ! shiftk(3,nshiftk), shiftks after inkpts
208 :
209 : real(dp),allocatable :: amu(:)
210 : ! amu(ntypat) ! EVOLVING variable
211 :
212 : real(dp), allocatable :: xred(:,:)
213 : ! EVOLVING variable xred(3,natom)
214 :
215 : real(dp), allocatable :: zionpsp(:)
216 : ! zionpsp(npsp) from psps
217 :
218 : real(dp), allocatable :: znuclpsp(:)
219 : ! znuclpsp(npsp) from psps
220 : ! Note the difference between (znucl|znucltypat) and znuclpsp !
221 :
222 : real(dp), allocatable :: znucltypat(:)
223 : ! znucltypat(ntypat) from alchemy
224 :
225 : character(len=8) :: codvsn
226 : ! version of the code
227 :
228 : character(len=132), allocatable :: title(:)
229 : ! title(npsp) from psps
230 :
231 : character(len=md5_slen),allocatable :: md5_pseudos(:)
232 : ! md5pseudos(npsp)
233 : ! md5 checksums associated to pseudos (read from file)
234 :
235 : ! EVOLVING variable, only for paw
236 : type(pawrhoij_type), allocatable :: pawrhoij(:)
237 :
238 : contains
239 :
240 : procedure :: init => hdr_init
241 : ! Initialize the header and most of its content from dtset and psps.
242 :
243 : procedure :: init_lowlvl => hdr_init_lowlvl
244 : ! Low level initialization method for Hdr (no dtset).
245 :
246 : procedure :: free => hdr_free
247 : ! Deallocates the components of the header.
248 :
249 : procedure :: get_nelect_from_occ => hdr_get_nelect_from_occ
250 : ! Returns the number of electrons calculated from the occupation factors Hdr%occ
251 :
252 : procedure :: from_fname => hdr_read_from_fname
253 : ! Read the header from file. Use Fortran IO or Netcdf depending on the extension of the file
254 :
255 : procedure :: fort_read => hdr_fort_read
256 : ! Reads the header from a logical unit associated to a Fortran unformatted file.
257 :
258 : procedure :: ncread => hdr_ncread
259 : ! Reads the header from a Netcdf file.
260 :
261 : procedure :: ncwrite => hdr_ncwrite
262 : ! Writes the header and fform to a Netcdf file.
263 :
264 : procedure :: vs_dtset => hdr_vs_dtset
265 : ! Check the compatibility of header with dtset.
266 :
267 : procedure :: get_crystal => hdr_get_crystal
268 : ! Return the crystal structure stored in the header.
269 :
270 : procedure :: bcast => hdr_bcast
271 : ! Broadcast the header.
272 :
273 : procedure :: compare => hdr_compare
274 : ! Compare two headers
275 :
276 : procedure :: update => hdr_update
277 : ! Update the header.
278 :
279 : procedure :: write_to_fname => hdr_write_to_fname
280 : ! Write the header (requires a string with the file name).
281 :
282 : procedure :: fort_write => hdr_fort_write
283 : ! Writes the header and fform to unformatted file
284 :
285 : procedure :: backspace => hdr_backspace
286 : ! Backspace the header (Fortran IO).
287 :
288 : procedure :: echo => hdr_echo
289 : ! Echo the header.
290 :
291 : procedure :: copy => hdr_copy
292 : ! Copy of the Header.
293 :
294 : end type hdr_type
295 : !!***
296 :
297 : public :: hdr_mpio_skip ! Skip the abinit header using MPI-IO routines.
298 : ! Return the offset of the first Fortran record after the header.
299 : public :: hdr_bsize_frecords ! Compute the size of the Fortran records from the header and formeig.
300 : public :: hdr_skip ! Skip the header.
301 : public :: hdr_io ! IO of the header.
302 : public :: hdr_check ! Compare two headers.
303 :
304 : public :: abifile_from_varname
305 : public :: abifile_from_fform
306 : public :: fform_from_ext ! Return the value of fform to be used from the file extension.
307 : public :: fform_contains ! Test whether fform provides varname.
308 : public :: varname_from_fname ! Return the name of the netcdf variable stored in a file from the file extension.
309 :
310 : ! Generic interface of the routines hdr_skip
311 : interface hdr_skip
312 : module procedure hdr_skip_int
313 : module procedure hdr_skip_wfftype
314 : end interface hdr_skip
315 :
316 : ! Generic interface of the routines hdr_io
317 : interface hdr_io
318 : module procedure hdr_io_int
319 : module procedure hdr_io_wfftype
320 : end interface hdr_io
321 :
322 : integer,private,parameter :: HDR_KNOWN_HEADFORMS(1) = [80]
323 : ! The list of headforms used so far.
324 :
325 : integer,private,parameter :: size_hdr_known_headforms = size(HDR_KNOWN_HEADFORMS) ! Need this for Flang
326 :
327 : ! The latest headform to be used for writing.
328 : integer,public,parameter :: HDR_LATEST_HEADFORM = HDR_KNOWN_HEADFORMS(size_hdr_known_headforms)
329 : !!***
330 :
331 : !!****t* m_hdr/abifile_t
332 : !! NAME
333 : !! abifile_t
334 : !!
335 : !! FUNCTION
336 : !! Gather information about a binary file with header.
337 : !! Every file with header must be registered in all_abifiles, see below.
338 : !!
339 : !! SOURCE
340 :
341 : type,public :: abifile_t
342 :
343 : character(len=nctk_slen) :: varname
344 : ! Name of the netcdf variable associated to the file.
345 : ! This string is used in fftdatar_write to find the value of fform to be written to file
346 :
347 : integer :: fform
348 : ! The value of fform associated to this file
349 :
350 : character(len=24) :: ext
351 : ! Abinit File extension (`.nc` is not included)
352 :
353 : character(len=24) :: class
354 : ! Each file belongs to a class e.g. wf_planewave, den, pot, data...
355 :
356 : logical :: has_pawrhoij=.True.
357 : ! True if this file contains pawrhoij when hdr%usepaw == 1.
358 :
359 : end type abifile_t
360 :
361 : ! Notes about abifiles:
362 : !
363 : ! *) fform are positive integers >0 and they must be unique inside the list.
364 : ! One might have used strings xxx.yyy.zzz instead of integers but there's a lot of code
365 : ! around that relies on integral fforms so we have to live with it.
366 : !
367 : ! *) class is used in postprocessing tools e.g. cut3d when we need to know if we are dealing
368 : ! with wavefunctions or density-like or potential-like data.
369 : ! Possible values are: "wf_planewave" for wavefunction files, "density" for density-like,
370 : ! "potential" for potential files, "data" for generic data a.k.a. internal files e.g. GKK matrix elements.
371 : !
372 : ! *) varname can appear multiple times, in this case the entries should be ordered chronologically
373 : ! i.e. the most recent format should follow the older ones. This could be useful if we decide to
374 : ! remove pawrhoij from a particular file. Let's assume, for example, that we've decided to remove
375 : ! pawrhoij from the POT file. In this case, abifiles should contain:
376 : !
377 : ! abifile_t(varname="potential", fform=102), & ! old file with pawrhoij
378 : ! abifile_t(varname="potential", fform=202, has_pawrhoij=.False.) & ! new file wo pawrhoij
379 : !
380 : ! *) The file extensions is used in fform_from_ext and varname_from_fname.
381 : ! fform_from_ext returns the most recent fform associated to a file extension.
382 : ! varname_from_fname is used in post-processing tools e.g. cut3d
383 : ! to read data from netcdf file without having to prompt the user for the variable name.
384 : ! In principle, the extension should be unique but there are exceptions e.g. the WFK produced in bigdft mode.
385 : ! Moreover the files produced by the DFPT code do not have a well-defined extension and, as a consequence,
386 : ! they require a special treatment. In python I would use regexp but Fortran is not python!
387 :
388 : type(abifile_t),private,parameter :: all_abifiles(52) = [ &
389 :
390 : ! Files with wavefunctions:
391 : abifile_t(varname="coefficients_of_wavefunctions", fform=2, ext="WFK", class="wf_planewave"), &
392 : abifile_t(varname="real_space_wavefunctions", fform=200, ext="WFK", class="wf_wavelet"), & ! Used by wavelets.
393 : abifile_t(varname="ur_ae", fform=602, ext="PAWAVES", class="wf_rspace"), & ! Used in pawmkaewf.
394 : abifile_t(varname="coefficients_of_wavefunctions", fform=502, ext="KSS", class="wf_planewave"), &
395 :
396 : ! Files with density-like data.
397 : abifile_t(varname="density", fform=52, ext="DEN", class="density"), & ! Official
398 : abifile_t(varname="positron_density", fform=53, ext="POSITRON", class="density"), &
399 : abifile_t(varname="first_order_density", fform=54, ext="DEN(\d+)", class="density"), &
400 : abifile_t(varname="pawrhor", fform=55, ext="PAWDEN", class="density"), &
401 : abifile_t(varname="pawrhor_core", fform=56, ext="ATMDEN_CORE", class="density"), &
402 : abifile_t(varname="pawrhor_val", fform=57, ext="ATMDEN_VAL", class="density"), &
403 : abifile_t(varname="pawrhor_full", fform=58, ext="ATMDEN_FULL", class="density"), &
404 : abifile_t(varname="pawrhor_ntilde_minus_nhat", fform=59, ext="N_TILDE", class="density"), &
405 : abifile_t(varname="pawrhor_n_one", fform=60, ext="N_ONE", class="density"), &
406 : abifile_t(varname="pawrhor_nt_one", fform=61, ext="NT_ONE", class="density"), &
407 : abifile_t(varname="qp_rhor", fform=62, ext="QP_DEN", class="density"), &
408 : abifile_t(varname="qp_pawrhor", fform=63, ext="QP_PAWDEN", class="density"), &
409 : abifile_t(varname="grhor_1", fform=67, ext="GDEN1", class="density"), &
410 : abifile_t(varname="grhor_2", fform=68, ext="GDEN2", class="density"), &
411 : abifile_t(varname="grhor_3", fform=69, ext="GDEN3", class="density"), &
412 :
413 : !???
414 : abifile_t(varname="stm", fform=110, ext="STM", class="density"), &
415 : abifile_t(varname="kinedr", fform=70, ext="KDEN", class="density"), &
416 : abifile_t(varname="elfr", fform=64, ext="ELF", class="density"), &
417 : abifile_t(varname="elfr_up", fform=65, ext="ELF_UP", class="density"), &
418 : abifile_t(varname="elfr_down", fform=66, ext="ELF_DOWN", class="density"), &
419 : abifile_t(varname="laprhor", fform=71, ext="LDEN", class="density"), &
420 :
421 : ! Files with potentials
422 : ! Official
423 : abifile_t(varname="potential", fform=102, ext="POT", class="potential"), & ! CHECK THESE TWO FILES
424 : abifile_t(varname="vtrial", fform=103, ext="POT", class="potential"), &
425 : abifile_t(varname="vhartree", fform=104, ext="VHA", class="potential"), &
426 : abifile_t(varname="vpsp", fform=105, ext="VPSP", class="potential"), &
427 : abifile_t(varname="vhartree_vloc", fform=106, ext="VCLMB", class="potential"), &
428 : abifile_t(varname="vhxc", fform=107, ext="VHXC", class="potential"), &
429 : abifile_t(varname="exchange_correlation_potential", fform=108, ext="VXC", class="potential"), &
430 :
431 : abifile_t(varname="first_order_potential", fform=109, ext="POT(\d+)", class="potential"), &
432 : ! fform 111 contains an extra record with rhog1_q(G=0) after the DFPT potential(r).
433 : abifile_t(varname="first_order_potential", fform=111, ext="POT(\d+)", class="potential"), &
434 :
435 : abifile_t(varname="first_order_vhartree", fform=112, ext="VHA(\d+)", class="potential"), &
436 : abifile_t(varname="first_order_vpsp", fform=113, ext="VPSP(\d+)", class="potential"), &
437 : abifile_t(varname="first_order_vxc", fform=114, ext="VXC(\d+)", class="potential"), &
438 :
439 : ! Data used in conducti
440 : abifile_t(varname="pawnabla", fform=610, ext="OPT1", class="data"), &
441 : abifile_t(varname="pawnabla_core", fform=611, ext="OPT2", class="data"), &
442 : abifile_t(varname="pawnabla_loc", fform=612, ext="OPT", class="data"), &
443 :
444 : ! Data used in E-PH code
445 : abifile_t(varname="gkk_elements", fform=42, ext="GKK", class="data"), &
446 :
447 : ! DKK matrix elements in netcdf format (optic, eph)
448 : abifile_t(varname="h1_matrix_elements", fform=43, ext="DKK", class="data"), &
449 :
450 : ! output files that are not supposed to be read by abinit.
451 : abifile_t(varname="this_file_is_not_read_by_abinit", fform=666, ext="666", class="data"), &
452 :
453 : ! GW files: old 1002, 1102
454 : !character(len=nctk_slen),public,parameter :: e_ncname="dielectric_function"
455 : ! FIXME This one should be rewritten
456 : abifile_t(varname="polarizability", fform=1003, ext="SUS", class="polariz"), &
457 : abifile_t(varname="inverse_dielectric_function", fform=1004, ext="SCR", class="epsm1"), &
458 : !abifile_t(varname="dielectric_function", fform=1002, ext="EPS", class="eps"), &
459 : !
460 : ! BSE: TODO. see m_bse_io
461 : !abifile_t(varname="bse_uresonant_q0", fform=1002, ext="BSR", class="bsreso"), &
462 : !abifile_t(varname="bse_ucoupling_q0", fform=1002, ext="BSC", class="bscoup"), &
463 :
464 : ! Miscellaneous
465 : abifile_t(varname="dos_fractions", fform=3000, ext="FATBANDS", class="data"), &
466 : abifile_t(varname="orbmag_mesh", fform=3100, ext="ORBMAG", class="data"), &
467 : abifile_t(varname="spectral_weights", fform=5000, ext="FOLD2BLOCH", class="data"), &
468 : abifile_t(varname="no_fftdatar_write", fform=6000, ext="ABIWAN", class="data"), &
469 : abifile_t(varname="None", fform=6001, ext="KERANGE", class="data"), &
470 : abifile_t(varname="None", fform=6002, ext="SIGEPH", class="data"), &
471 : abifile_t(varname="None", fform=6003, ext="GSTORE", class="data") &
472 : ]
473 :
474 : type(abifile_t),public,parameter :: abifile_none = abifile_t(varname="None", fform=0, ext="None", class="None")
475 : ! This object is returned when we cannot find the file in abifiles.
476 :
477 : CONTAINS !===========================================================
478 : !!***
479 :
480 : !----------------------------------------------------------------------
481 :
482 : !!****f* m_hdr/fform_from_ext
483 : !! NAME
484 : !! fform_from_ext
485 : !!
486 : !! FUNCTION
487 : !! Return the value of fform to be used from the file extension. If a file has multiple fforms,
488 : !! the most recent one is returned. Returns 0 if the extension is not registered.
489 : !!
490 : !! SOURCE
491 :
492 135 : integer function fform_from_ext(abiext) result(fform)
493 :
494 : !Arguments ---------------------------------------------
495 : character(len=*),intent(in) :: abiext
496 :
497 : !Local variables-------------------------------
498 : !scalars
499 : integer :: ii,ind,ierr,pertcase
500 135 : character(len=len(abiext)) :: ext
501 :
502 : ! *********************************************************************
503 : ! Remove .nc (if any) and work with ext
504 135 : ext = abiext
505 135 : if (endswith(abiext, ".nc")) then
506 120 : ind = index(abiext, ".nc", back=.True.); ext = abiext(:ind-1)
507 : end if
508 :
509 135 : fform = 0
510 7155 : do ii=1,size(all_abifiles)
511 7155 : if (ext == all_abifiles(ii)%ext) fform = all_abifiles(ii)%fform
512 : end do
513 135 : if (fform /= 0) return
514 : ! Here we handle special cases.
515 :
516 : ! Handle DEN[pertcase]
517 0 : if (startswith(ext, "DEN")) then
518 0 : read(ext(4:), *, iostat=ierr) pertcase
519 0 : if (ierr == 0) then
520 0 : do ii=1,size(all_abifiles)
521 0 : if (all_abifiles(ii)%ext == "DEN(\d+)") fform = all_abifiles(ii)%fform
522 : end do
523 : return
524 : end if
525 : end if
526 :
527 : ! Handle POT[pertcase]
528 0 : if (startswith(ext, "POT")) then
529 0 : read(ext(4:), *, iostat=ierr) pertcase
530 0 : if (ierr == 0) then
531 0 : do ii=1,size(all_abifiles)
532 0 : if (all_abifiles(ii)%ext == "POT(\d+)") fform = all_abifiles(ii)%fform
533 : end do
534 : return
535 : end if
536 : end if
537 :
538 0 : ABI_ERROR(sjoin("Cannot find fform associated to extension:", abiext))
539 :
540 135 : end function fform_from_ext
541 : !!***
542 :
543 : !!****f* m_hdr/fform_contains
544 : !! NAME
545 : !! fform_contains
546 : !!
547 : !! FUNCTION
548 : !! Test whether fform provides varname.
549 : !!
550 : !! SOURCE
551 :
552 2 : logical function fform_contains(fform, varname, msg) result(ok)
553 :
554 : !Arguments ---------------------------------------------
555 : integer,intent(in) :: fform
556 : character(len=*),intent(in) :: varname
557 : character(len=*),intent(out) :: msg
558 :
559 : !Local variables-------------------------------
560 : integer :: ii
561 : ! *********************************************************************
562 :
563 2 : ok = .False.; msg = ""
564 :
565 2 : if (fform == 0) then
566 0 : msg = "fform == 0 --> IO errror when reading the header from file"
567 0 : return
568 : end if
569 :
570 54 : do ii=1,size(all_abifiles)
571 54 : if (fform == all_abifiles(ii)%fform .and. varname == all_abifiles(ii)%varname) then
572 2 : ok = .True.; return
573 : end if
574 : end do
575 :
576 0 : do ii=1,size(all_abifiles)
577 0 : if (fform == all_abifiles(ii)%fform) then
578 0 : msg = sjoin("fform:", itoa(fform), " contains", all_abifiles(ii)%varname, "instead of: ", varname)
579 0 : return
580 : end if
581 : end do
582 :
583 0 : msg = sjoin("Cannot find fform:", itoa(fform), " in all_abifiles")
584 :
585 2 : end function fform_contains
586 : !!***
587 :
588 : !----------------------------------------------------------------------
589 :
590 : !!****f* m_hdr/varname_from_fname
591 : !! NAME
592 : !! varname_from_fname
593 : !!
594 : !! FUNCTION
595 : !! Return the name of the netcdf variable stored in a file from the file extension.
596 : !!
597 : !! NOTES
598 : !! The variable names should be consistent with the ones used in outscf.F90
599 : !!
600 : !! SOURCE
601 :
602 2002 : character(len=nctk_slen) function varname_from_fname(filename) result(varname)
603 :
604 : !Arguments ---------------------------------------------
605 : character(len=*),intent(in) :: filename
606 :
607 : !Local variables-------------------------------
608 : !scalars
609 : integer :: ind,pertcase,ierr
610 : logical :: found
611 2002 : character(len=len(filename)) :: ext
612 : ! *********************************************************************
613 :
614 : ! TODO: This should be a recursive function because we have
615 : ! to scan the string from left to right (extensions could have the same termination)
616 :
617 : ! Find the Abinit file extension. Examples: t43_VHXC.nc
618 2002 : if (endswith(filename, ".nc")) then
619 2002 : ind = index(filename, ".nc", back=.True.)
620 : else
621 : !ABI_ERROR(sjoin("Don't know how to handle: ", filename))
622 0 : ind = len_trim(filename) + 1
623 : end if
624 :
625 2002 : ext = filename(:ind-1)
626 2002 : ind = index(ext, "_", back=.True.)
627 2002 : ABI_CHECK(ind /= 0, "Cannot find `_` in file name!")
628 2002 : ABI_CHECK(ind /= len_trim(ext), sjoin("Wrong string: ", ext))
629 2002 : ext = ext(ind+1:)
630 :
631 2002 : found = .True.
632 597 : select case (ext)
633 : case ("DEN")
634 597 : varname = "density"
635 : !case ("DEN1")
636 : ! varname = "first_order_density"
637 : case ("POSITRON")
638 1 : varname = "positron_density"
639 : case ("PAWDEN")
640 1 : varname = "pawrhor"
641 : ! TODO: Other paw densities
642 : case ("ELF")
643 0 : varname = "elfr"
644 : case ("ELF_UP")
645 0 : varname = "elfr_up"
646 : case ("ELF_DOWN")
647 0 : varname = "elfr_down"
648 : case ("GDEN1")
649 0 : varname = "grhor_1"
650 : case ("GDEN2")
651 0 : varname = "grhor_2"
652 : case ("GDEN3")
653 0 : varname = "grhor_3"
654 : case ("KDEN")
655 3 : varname = "kinedr"
656 : case ("LDEN")
657 0 : varname = "laprhor"
658 : case ("POT")
659 30 : varname = "vtrial"
660 : case ("STM")
661 0 : varname = "stm"
662 : case ("VHA")
663 1 : varname = "vhartree"
664 : case ("VPSP")
665 0 : varname = "vpsp"
666 : case ("VHXC")
667 3 : varname = "vhxc"
668 : case ("VXC")
669 0 : varname = "exchange_correlation_potential"
670 : case ("VCLMB")
671 0 : varname = "vhartree_vloc"
672 : case default
673 2002 : found = .False.
674 : end select
675 :
676 2002 : if (found) return
677 :
678 : ! Handle DEN[pertcase]
679 1366 : if (startswith(ext, "DEN")) then
680 1366 : read(ext(4:), *, iostat=ierr) pertcase
681 1366 : if (ierr == 0) then
682 1366 : varname = "first_order_density"; return
683 : end if
684 : end if
685 :
686 : ! Handle POT[pertcase]
687 0 : if (startswith(ext, "POT")) then
688 0 : read(ext(4:), *, iostat=ierr) pertcase
689 0 : if (ierr == 0) then
690 0 : varname = "first_order_potential"; return
691 : end if
692 : end if
693 :
694 : ! Handle VXC[pertcase]
695 0 : if (startswith(ext, "VXC")) then
696 0 : read(ext(4:), *, iostat=ierr) pertcase
697 0 : if (ierr == 0) then
698 0 : varname = "first_order_vxc"; return
699 : end if
700 : end if
701 :
702 : ! Handle VHA[pertcase]
703 0 : if (startswith(ext, "VHA")) then
704 0 : read(ext(4:), *, iostat=ierr) pertcase
705 0 : if (ierr == 0) then
706 0 : varname = "first_order_vhartree"; return
707 : end if
708 : end if
709 :
710 : ! Handle VPSP[pertcase]
711 0 : if (startswith(ext, "VPSP")) then
712 0 : read(ext(4:), *, iostat=ierr) pertcase
713 0 : if (ierr == 0) then
714 0 : varname = "first_order_vpsp"; return
715 : end if
716 : end if
717 :
718 0 : ABI_ERROR(sjoin("Unknown abinit extension:", ext))
719 :
720 : end function varname_from_fname
721 : !!***
722 :
723 : !----------------------------------------------------------------------
724 :
725 : !!****f* m_hdr/abifile_from_varname
726 : !! NAME
727 : !! abifile_from_varname
728 : !!
729 : !! FUNCTION
730 : !! Return the abifile_t object corresponding to the given variable varname
731 : !! Return abifile_none if not found. This function is used to find the last
732 : !! value of fform when we write data to file.
733 : !!
734 : !! SOURCE
735 :
736 7962 : type(abifile_t) function abifile_from_varname(varname) result(afile)
737 :
738 : !Arguments ---------------------------------------------
739 : character(len=*),intent(in) :: varname
740 :
741 : !Local variables-------------------------------
742 : !scalars
743 : integer :: ii
744 : ! *************************************************************************
745 :
746 7962 : afile = abifile_none
747 421986 : do ii=1,size(all_abifiles)
748 421986 : if (all_abifiles(ii)%varname == varname) afile = all_abifiles(ii)
749 : end do
750 :
751 7962 : end function abifile_from_varname
752 : !!***
753 :
754 : !----------------------------------------------------------------------
755 :
756 : !!****f* m_hdr/abifile_from_fform
757 : !! NAME
758 : !! abifile_from_fform
759 : !!
760 : !! FUNCTION
761 : !! Return the abifile_t object corresponding to the given fform
762 : !! Return abifile_none if not found. This function is used to
763 : !! find the name of the netcdf variable from the fform and
764 : !! detect whether the file contains pawrhoij.
765 : !!
766 : !! SOURCE
767 :
768 8512 : type(abifile_t) function abifile_from_fform(fform) result(afile)
769 :
770 : !Arguments ---------------------------------------------
771 : integer,intent(in) :: fform
772 :
773 : !Local variables-------------------------------
774 : integer :: ii
775 : ! *************************************************************************
776 :
777 8512 : afile = abifile_none
778 451136 : do ii=1,size(all_abifiles)
779 451136 : if (all_abifiles(ii)%fform == fform) afile = all_abifiles(ii)
780 : end do
781 :
782 8512 : end function abifile_from_fform
783 : !!***
784 :
785 : !!****f* m_hdr/check_fform
786 : !! NAME
787 : !! check_fform
788 : !!
789 : !! FUNCTION
790 : !! This function is used ifdef DEBUG_MODE. It tests whether the value of fform is registered in all_abifiles.
791 : !!
792 : !! SOURCE
793 :
794 : subroutine check_fform(fform)
795 :
796 : !Local variables-------------------------------
797 : integer,intent(in) :: fform
798 : #ifdef DEBUG_MODE
799 : type(abifile_t) :: abifile
800 : character(len=500) :: msg
801 :
802 : ! *********************************************************************
803 : if (fform == 666) return
804 : abifile = abifile_from_fform(fform)
805 :
806 : if (abifile%fform == 0) then
807 : ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform)))
808 : end if
809 : if (abifile%fform /= fform) then
810 : write(msg,"(2a,2(a,i0))") &
811 : "Input fform does not agree with the one registered in abifile.",ch10,&
812 : "hdr%fform= ",fform,", abifile%fform= ",abifile%fform
813 : ABI_ERROR(msg)
814 : end if
815 :
816 : #else
817 : ABI_UNUSED(fform)
818 : #endif
819 :
820 : end subroutine check_fform
821 : !!***
822 :
823 : !----------------------------------------------------------------------
824 :
825 : !!****f* m_hdr/test_abifiles
826 : !! NAME
827 : !! test_abifiles
828 : !!
829 : !! FUNCTION
830 : !! Check the consistency of the internal abifiles table.
831 : !!
832 : !! SOURCE
833 :
834 : subroutine test_abifiles()
835 :
836 : !Arguments ---------------------------------------------
837 :
838 : !Local variables-------------------------------
839 : !scalars
840 : integer :: ii,nn,ierr
841 : integer :: all_fforms(size(all_abifiles)),iperm(size(all_abifiles))
842 : ! *************************************************************************
843 : nn = size(all_abifiles)
844 :
845 : do ii=1,nn
846 : all_fforms(ii) = all_abifiles(ii)%fform
847 : end do
848 : iperm = [(ii, ii=1,nn)]
849 : call sort_int(nn, all_fforms, iperm)
850 :
851 : ierr = 0
852 : do ii=1,nn-1
853 : if (all_fforms(ii) == all_fforms(ii+1)) then
854 : ABI_WARNING(sjoin("fform: ", itoa(all_fforms(ii+1)), "is already in the abifiles list"))
855 : ierr = ierr + 1
856 : end if
857 : end do
858 :
859 : if (ierr /= 0) then
860 : ABI_ERROR("test_abifiles gave ierr != 0. Aborting now")
861 : end if
862 :
863 : end subroutine test_abifiles
864 : !!***
865 :
866 : !----------------------------------------------------------------------
867 :
868 : !!****f* m_hdr/hdr_malloc
869 : !! NAME
870 : !! hdr_malloc
871 : !!
872 : !! FUNCTION
873 : !! Allocate memory from dimensions with the exception of pawrhoij.
874 : !! This is a private routine. Client code should use hdr_init or hdr_fort_read.
875 : !!
876 : !! SOURCE
877 :
878 38801 : subroutine hdr_malloc(hdr, bantot, nkpt, nsppol, npsp, natom, ntypat, nsym, nshiftk_orig, nshiftk)
879 :
880 : !Arguments ---------------------------------------------
881 : !scalars
882 : class(hdr_type),intent(inout) :: hdr
883 : integer,intent(in) :: bantot,nkpt,nsppol,npsp,natom,ntypat,nsym,nshiftk,nshiftk_orig
884 : ! *************************************************************************
885 :
886 : !@hdt_type
887 38801 : call hdr%free()
888 :
889 116403 : ABI_MALLOC(hdr%istwfk, (nkpt))
890 116403 : ABI_MALLOC(hdr%nband, (nkpt*nsppol))
891 77602 : ABI_MALLOC(hdr%npwarr, (nkpt))
892 116403 : ABI_MALLOC(hdr%pspcod, (npsp))
893 77602 : ABI_MALLOC(hdr%pspdat, (npsp))
894 77602 : ABI_MALLOC(hdr%pspso, (npsp))
895 77602 : ABI_MALLOC(hdr%pspxc, (npsp))
896 77602 : ABI_MALLOC(hdr%lmn_size, (npsp))
897 77602 : ABI_MALLOC(hdr%so_psp, (npsp))
898 116403 : ABI_MALLOC(hdr%symafm, (nsym))
899 116403 : ABI_MALLOC(hdr%symrel, (3,3,nsym))
900 116403 : ABI_MALLOC(hdr%typat, (natom))
901 116403 : ABI_MALLOC(hdr%kptns, (3,nkpt))
902 116403 : ABI_MALLOC(hdr%occ, (bantot))
903 116403 : ABI_MALLOC(hdr%tnons, (3,nsym))
904 116403 : ABI_MALLOC(hdr%wtk, (nkpt))
905 116403 : ABI_MALLOC(hdr%xred, (3,natom))
906 116403 : ABI_MALLOC(hdr%zionpsp, (npsp))
907 77602 : ABI_MALLOC(hdr%znuclpsp, (npsp))
908 116403 : ABI_MALLOC(hdr%znucltypat, (ntypat))
909 116403 : ABI_MALLOC(hdr%title, (npsp))
910 116403 : ABI_MALLOC(hdr%shiftk, (3,nshiftk))
911 116403 : ABI_MALLOC(hdr%shiftk_orig, (3,nshiftk_orig))
912 116403 : ABI_MALLOC(hdr%md5_pseudos, (npsp))
913 77602 : ABI_MALLOC(hdr%amu, (ntypat))
914 :
915 38801 : end subroutine hdr_malloc
916 : !!***
917 :
918 : !----------------------------------------------------------------------
919 :
920 : !!****f* m_hdr/hdr_init
921 : !! NAME
922 : !! hdr_init
923 : !!
924 : !! FUNCTION
925 : !! Initialize the header structured datatype
926 : !! and most of its content from dtset and psps, and put default values for evolving variables.
927 : !!
928 : !! INPUTS
929 : !! ebands <type(ebands_base_t)>=band structure information including Brillouin zone description
930 : !! codvsn=code version
931 : !! dtset <type(dataset_type)>=all input variables for this dataset
932 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
933 : !! pertcase=index of the perturbation, or 0 if GS calculation
934 : !! psps <type(pseudopotential_type)>=all the information about psps
935 : !! [my_atomtab(:)]=Index of the atoms (in global numbering ) treated by current proc.
936 : !! [mpi_atmtab(:)]= indexes of the atoms treated by current proc
937 : !! [comm_atom]= MPI communicator over atoms
938 : !!
939 : !! OUTPUT
940 : !! hdr <type(hdr_type)>=the header, initialized, and for most part of
941 : !! it, contain its definite values, except for evolving variables
942 : !!
943 : !! SOURCE
944 :
945 23066 : subroutine hdr_init(hdr, ebands, codvsn, dtset, pawtab, pertcase, psps, wvl, &
946 11533 : mpi_atmtab, comm_atom) ! optional arguments (parallelism)
947 :
948 : !Arguments ------------------------------------
949 : !scalars
950 : class(hdr_type),intent(inout) :: hdr !vz_i
951 : class(ebands_base_t),intent(in) :: ebands
952 : integer,intent(in) :: pertcase
953 : integer,intent(in),optional :: comm_atom
954 : character(len=8),intent(in) :: codvsn
955 : type(dataset_type),intent(in) :: dtset
956 : type(pseudopotential_type),intent(in) :: psps
957 : type(wvl_internal_type),intent(in) :: wvl
958 : !arrays
959 : integer,optional,target,intent(in) :: mpi_atmtab(:)
960 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
961 :
962 : !Local variables-------------------------------
963 : !scalars
964 : integer,parameter :: image=1
965 : character(len=500) :: msg
966 : ! *************************************************************************
967 :
968 : #ifdef DEBUG_MODE
969 : call test_abifiles()
970 : #endif
971 :
972 : !@hdr_type
973 : ! More checking would be needed ...
974 11533 : if (dtset%ntypat/=psps%ntypat) then
975 0 : write(msg,'(a,2(i0,2x))')' dtset%ntypat and psps%ntypat differs. They are: ',dtset%ntypat,psps%ntypat
976 0 : ABI_ERROR(msg)
977 : end if
978 :
979 11533 : if (dtset%npsp/=psps%npsp) then
980 0 : write(msg,'(a,2(i0,2x))')' dtset%npsp and psps%npsp differs. They are: ',dtset%npsp,psps%npsp
981 0 : ABI_ERROR(msg)
982 : end if
983 :
984 : ! Note: The structure parameters are taken from the first image, also cellcharge !
985 11533 : if (present(comm_atom)) then
986 11155 : if (present(mpi_atmtab)) then
987 : call hdr%init_lowlvl(ebands,psps,pawtab,wvl,codvsn,pertcase,&
988 : dtset%natom,dtset%nsym,dtset%nspden,dtset%ecut,dtset%pawecutdg,dtset%ecutsm,dtset%dilatmx,&
989 : dtset%intxc,dtset%ixc,dtset%stmbias,dtset%usewvl,dtset%pawcpxocc,dtset%pawspnorb,dtset%ngfft,dtset%ngfftdg,&
990 : dtset%so_psp,dtset%qptn, dtset%rprimd_orig(:,:,image),dtset%xred_orig(:,:,image),&
991 : dtset%symrel,dtset%tnons,dtset%symafm,dtset%typat,dtset%amu_orig(:,image),dtset%icoulomb,&
992 : dtset%kptopt,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,dtset%cellcharge(image),&
993 : dtset%kptrlatt_orig,dtset%kptrlatt,dtset%nshiftk_orig,dtset%nshiftk,dtset%shiftk_orig,dtset%shiftk,&
994 738 : comm_atom=comm_atom,mpi_atmtab=mpi_atmtab)
995 : else
996 : call hdr%init_lowlvl(ebands,psps,pawtab,wvl,codvsn,pertcase,&
997 : dtset%natom,dtset%nsym,dtset%nspden,dtset%ecut,dtset%pawecutdg,dtset%ecutsm,dtset%dilatmx,&
998 : dtset%intxc,dtset%ixc,dtset%stmbias,dtset%usewvl,dtset%pawcpxocc,dtset%pawspnorb,dtset%ngfft,dtset%ngfftdg,&
999 : dtset%so_psp,dtset%qptn, dtset%rprimd_orig(:,:,image),dtset%xred_orig(:,:,image),&
1000 : dtset%symrel,dtset%tnons,dtset%symafm,dtset%typat,dtset%amu_orig(:,image),dtset%icoulomb,&
1001 : dtset%kptopt,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,dtset%cellcharge(image),&
1002 : dtset%kptrlatt_orig,dtset%kptrlatt,dtset%nshiftk_orig,dtset%nshiftk,dtset%shiftk_orig,dtset%shiftk,&
1003 10417 : comm_atom=comm_atom)
1004 : end if
1005 : else
1006 : call hdr%init_lowlvl(ebands,psps,pawtab,wvl,codvsn,pertcase,&
1007 : dtset%natom,dtset%nsym,dtset%nspden,dtset%ecut,dtset%pawecutdg,dtset%ecutsm,dtset%dilatmx,&
1008 : dtset%intxc,dtset%ixc,dtset%stmbias,dtset%usewvl,dtset%pawcpxocc,dtset%pawspnorb,dtset%ngfft,dtset%ngfftdg,&
1009 : dtset%so_psp,dtset%qptn, dtset%rprimd_orig(:,:,image),dtset%xred_orig(:,:,image),dtset%symrel,&
1010 : dtset%tnons,dtset%symafm,dtset%typat,dtset%amu_orig(:,image),dtset%icoulomb,&
1011 : dtset%kptopt,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,dtset%cellcharge(image),&
1012 378 : dtset%kptrlatt_orig,dtset%kptrlatt,dtset%nshiftk_orig,dtset%nshiftk,dtset%shiftk_orig,dtset%shiftk)
1013 : end if
1014 :
1015 11533 : end subroutine hdr_init
1016 : !!***
1017 :
1018 : !----------------------------------------------------------------------
1019 :
1020 : !!****f* m_hdr/hdr_free
1021 : !! NAME
1022 : !! hdr_free
1023 : !!
1024 : !! FUNCTION
1025 : !! This subroutine deallocates the components of the header structured datatype
1026 : !!
1027 : !! SOURCE
1028 :
1029 86946 : subroutine hdr_free(hdr)
1030 :
1031 : !Arguments ------------------------------------
1032 : class(hdr_type),intent(inout) :: hdr
1033 : ! *************************************************************************
1034 :
1035 : !integer
1036 86946 : ABI_SFREE(hdr%istwfk)
1037 86946 : ABI_SFREE(hdr%lmn_size)
1038 86946 : ABI_SFREE(hdr%nband)
1039 86946 : ABI_SFREE(hdr%npwarr)
1040 86946 : ABI_SFREE(hdr%pspcod)
1041 86946 : ABI_SFREE(hdr%pspdat)
1042 86946 : ABI_SFREE(hdr%pspso)
1043 86946 : ABI_SFREE(hdr%pspxc)
1044 86946 : ABI_SFREE(hdr%so_psp)
1045 86946 : ABI_SFREE(hdr%symafm)
1046 86946 : ABI_SFREE(hdr%symrel)
1047 86946 : ABI_SFREE(hdr%typat)
1048 :
1049 : !real
1050 86946 : ABI_SFREE(hdr%amu)
1051 86946 : ABI_SFREE(hdr%kptns)
1052 86946 : ABI_SFREE(hdr%occ)
1053 86946 : ABI_SFREE(hdr%tnons)
1054 86946 : ABI_SFREE(hdr%wtk)
1055 86946 : ABI_SFREE(hdr%shiftk)
1056 86946 : ABI_SFREE(hdr%shiftk_orig)
1057 86946 : ABI_SFREE(hdr%xred)
1058 86946 : ABI_SFREE(hdr%zionpsp)
1059 86946 : ABI_SFREE(hdr%znuclpsp)
1060 86946 : ABI_SFREE(hdr%znucltypat)
1061 :
1062 : !string arrays
1063 86946 : ABI_SFREE(hdr%md5_pseudos)
1064 86946 : ABI_SFREE(hdr%title)
1065 :
1066 86946 : if (hdr%usepaw==1 .and. allocated(hdr%pawrhoij) ) then
1067 5489 : call pawrhoij_free(hdr%pawrhoij)
1068 19483 : ABI_FREE(hdr%pawrhoij)
1069 : end if
1070 :
1071 86946 : end subroutine hdr_free
1072 : !!***
1073 :
1074 : !----------------------------------------------------------------------
1075 :
1076 : !!****f* m_hdr/hdr_copy
1077 : !! NAME
1078 : !! hdr_copy
1079 : !!
1080 : !! FUNCTION
1081 : !! Deep copy of the abinit header.
1082 : !!
1083 : !! NOTES
1084 : !! The present version deals with versions of the header up to 56.
1085 : !!
1086 : !! SOURCE
1087 :
1088 3337 : subroutine hdr_copy(Hdr_in, Hdr_cp)
1089 :
1090 : !Arguments ------------------------------------
1091 : !scalars
1092 : class(hdr_type),intent(in) :: Hdr_in
1093 : class(hdr_type),intent(out) :: Hdr_cp
1094 :
1095 : !Local variables-------------------------------
1096 : !scalars
1097 : integer :: cplex_rhoij,nspden_rhoij,qphase_rhoij
1098 : ! *************************************************************************
1099 :
1100 : !@hdr_type
1101 : ! Integer values
1102 3337 : Hdr_cp%bantot = Hdr_in%bantot
1103 3337 : Hdr_cp%date = Hdr_in%date
1104 3337 : Hdr_cp%headform = Hdr_in%headform
1105 3337 : hdr_cp%icoulomb = hdr_in%icoulomb
1106 3337 : Hdr_cp%intxc = Hdr_in%intxc
1107 3337 : Hdr_cp%ivalence = Hdr_in%ivalence
1108 3337 : Hdr_cp%ixc = Hdr_in%ixc
1109 3337 : Hdr_cp%natom = Hdr_in%natom
1110 3337 : Hdr_cp%nkpt = Hdr_in%nkpt
1111 3337 : Hdr_cp%npsp = Hdr_in%npsp
1112 3337 : Hdr_cp%nspden = Hdr_in%nspden
1113 3337 : Hdr_cp%nspinor = Hdr_in%nspinor
1114 3337 : Hdr_cp%nsppol = Hdr_in%nsppol
1115 3337 : Hdr_cp%nsym = Hdr_in%nsym
1116 3337 : Hdr_cp%ntypat = Hdr_in%ntypat
1117 3337 : Hdr_cp%occopt = Hdr_in%occopt
1118 3337 : Hdr_cp%pertcase = Hdr_in%pertcase
1119 3337 : Hdr_cp%usepaw = Hdr_in%usepaw
1120 3337 : Hdr_cp%usewvl = Hdr_in%usewvl
1121 3337 : Hdr_cp%mband = Hdr_in%mband
1122 179021 : ABI_CHECK(hdr_in%mband == maxval(hdr_in%nband), "mband != maxval(hdr_in%nband)")
1123 3337 : hdr_cp%kptopt = hdr_in%kptopt
1124 3337 : hdr_cp%pawcpxocc = hdr_in%pawcpxocc
1125 3337 : hdr_cp%nshiftk_orig = hdr_in%nshiftk_orig
1126 3337 : hdr_cp%nshiftk = hdr_in%nshiftk
1127 :
1128 : ! Integer arrays
1129 13348 : Hdr_cp%ngfft = Hdr_in%ngfft
1130 10011 : Hdr_cp%nwvlarr = Hdr_in%nwvlarr
1131 43381 : hdr_cp%kptrlatt = hdr_in%kptrlatt
1132 43381 : hdr_cp%kptrlatt_orig = hdr_in%kptrlatt_orig
1133 :
1134 : ! Integer allocatable arrays
1135 3337 : call alloc_copy( Hdr_in%istwfk, Hdr_cp%istwfk )
1136 3337 : call alloc_copy( Hdr_in%lmn_size,Hdr_cp%lmn_size )
1137 3337 : call alloc_copy( Hdr_in%nband, Hdr_cp%nband )
1138 3337 : call alloc_copy( Hdr_in%npwarr, Hdr_cp%npwarr )
1139 3337 : call alloc_copy( Hdr_in%pspcod, Hdr_cp%pspcod )
1140 3337 : call alloc_copy( Hdr_in%pspdat, Hdr_cp%pspdat )
1141 3337 : call alloc_copy( Hdr_in%pspso , Hdr_cp%pspso )
1142 3337 : call alloc_copy( Hdr_in%pspxc , Hdr_cp%pspxc )
1143 3337 : call alloc_copy( Hdr_in%so_psp, Hdr_cp%so_psp )
1144 3337 : call alloc_copy( Hdr_in%symafm, Hdr_cp%symafm )
1145 3337 : call alloc_copy( Hdr_in%symrel, Hdr_cp%symrel )
1146 3337 : call alloc_copy( Hdr_in%typat , Hdr_cp%typat )
1147 :
1148 : ! Real variables
1149 3337 : Hdr_cp%ecut = Hdr_in%ecut
1150 3337 : Hdr_cp%ecutdg = Hdr_in%ecutdg
1151 3337 : Hdr_cp%ecutsm = Hdr_in%ecutsm
1152 3337 : Hdr_cp%ecut_eff = Hdr_in%ecut_eff
1153 3337 : Hdr_cp%etot = Hdr_in%etot
1154 3337 : Hdr_cp%fermie = Hdr_in%fermie
1155 3337 : Hdr_cp%fermih = Hdr_in%fermih
1156 3337 : Hdr_cp%residm = Hdr_in%residm
1157 3337 : Hdr_cp%stmbias = Hdr_in%stmbias
1158 3337 : Hdr_cp%tphysel = Hdr_in%tphysel
1159 3337 : Hdr_cp%tsmear = Hdr_in%tsmear
1160 3337 : hdr_cp%nelect = hdr_in%nelect
1161 3337 : hdr_cp%ne_qFD = hdr_in%ne_qFD
1162 3337 : hdr_cp%nh_qFD = hdr_in%nh_qFD
1163 3337 : hdr_cp%cellcharge = hdr_in%cellcharge
1164 3337 : Hdr_cp%extfpmd_eshift = Hdr_in%extfpmd_eshift
1165 :
1166 13348 : Hdr_cp%qptn(:) = Hdr_in%qptn(:)
1167 43381 : Hdr_cp%rprimd(:,:) = Hdr_in%rprimd(:,:)
1168 :
1169 : ! Real allocatable arrays
1170 3337 : call alloc_copy(Hdr_in%amu, Hdr_cp%amu)
1171 3337 : call alloc_copy( Hdr_in%kptns ,Hdr_cp%kptns )
1172 3337 : call alloc_copy( Hdr_in%occ ,Hdr_cp%occ )
1173 3337 : call alloc_copy( Hdr_in%tnons ,Hdr_cp%tnons )
1174 3337 : call alloc_copy( Hdr_in%wtk ,Hdr_cp%wtk )
1175 3337 : call alloc_copy( Hdr_in%xred ,Hdr_cp%xred )
1176 3337 : call alloc_copy( Hdr_in%zionpsp ,Hdr_cp%zionpsp )
1177 3337 : call alloc_copy( Hdr_in%znuclpsp ,Hdr_cp%znuclpsp )
1178 3337 : call alloc_copy( Hdr_in%znucltypat,Hdr_cp%znucltypat)
1179 3337 : call alloc_copy(Hdr_in%shiftk, Hdr_cp%shiftk)
1180 3337 : call alloc_copy(Hdr_in%shiftk_orig, Hdr_cp%shiftk_orig)
1181 :
1182 : ! Character arrays
1183 3337 : Hdr_cp%codvsn = Hdr_in%codvsn
1184 : ! THIS DOES NOT WORK ON XLF: Hdr_cp%title string length becomes huge and segfaults
1185 : ! call alloc_copy( Hdr_in%title,Hdr_cp%title )
1186 10011 : ABI_MALLOC(Hdr_cp%title,(Hdr_cp%npsp))
1187 11289 : Hdr_cp%title = Hdr_in%title
1188 :
1189 10011 : ABI_MALLOC(hdr_cp%md5_pseudos, (hdr_cp%npsp))
1190 11289 : hdr_cp%md5_pseudos = hdr_in%md5_pseudos
1191 :
1192 : ! For PAW have to copy Pawrhoij ====
1193 : ! NOTE alchemy requires a different treatment but for the moment it is not available within PAW.
1194 3337 : if (Hdr_in%usepaw==1) then
1195 252 : cplex_rhoij = Hdr_in%Pawrhoij(1)%cplex_rhoij
1196 252 : qphase_rhoij = Hdr_in%Pawrhoij(1)%qphase
1197 252 : nspden_rhoij = Hdr_in%Pawrhoij(1)%nspden
1198 1479 : ABI_MALLOC(Hdr_cp%Pawrhoij,(Hdr_in%natom))
1199 : call pawrhoij_alloc(Hdr_cp%Pawrhoij,cplex_rhoij,nspden_rhoij,Hdr_in%nspinor,Hdr_in%nsppol,Hdr_in%typat,&
1200 252 : lmnsize=Hdr_in%lmn_size(1:Hdr_in%ntypat),qphase=qphase_rhoij)
1201 252 : call pawrhoij_copy(Hdr_in%Pawrhoij,Hdr_cp%Pawrhoij)
1202 : end if
1203 :
1204 3337 : end subroutine hdr_copy
1205 : !!***
1206 :
1207 : !----------------------------------------------------------------------
1208 :
1209 : !!****f* m_hdr/hdr_get_nelect_from_occ
1210 : !! NAME
1211 : !! hdr_get_nelect_from_occ
1212 : !!
1213 : !! FUNCTION
1214 : !! Return the number of electrons from the occupation numbers
1215 : !! This function is mainly used for debugging purposes, use hdr%nelect and hdr%cellcharge
1216 : !!
1217 : !! INPUTS
1218 : !! Hdr<hdr_type>
1219 : !!
1220 : !! OUTPUT
1221 : !! nelect=Number of electrons in the unit cell.
1222 : !!
1223 : !! SOURCE
1224 :
1225 0 : real(dp) pure function hdr_get_nelect_from_occ(Hdr) result(nelect)
1226 :
1227 : !Arguments ---------------------------------------------
1228 : !scalars
1229 : class(hdr_type),intent(in) :: Hdr
1230 :
1231 : !Local variables ---------------------------------------
1232 : !scalars
1233 : integer :: idx,isppol,ikibz,nband_k
1234 : ! *************************************************************************
1235 :
1236 : ! Cannot use znucl because we might have additional cellcharge or alchemy.
1237 0 : nelect=zero ; idx=0
1238 0 : do isppol=1,Hdr%nsppol
1239 0 : do ikibz=1,Hdr%nkpt
1240 0 : nband_k=Hdr%nband(ikibz+(isppol-1)*Hdr%nkpt)
1241 0 : nelect = nelect + Hdr%wtk(ikibz)*SUM(Hdr%occ(idx+1:idx+nband_k))
1242 0 : idx=idx+nband_k
1243 : end do
1244 : end do
1245 :
1246 0 : end function hdr_get_nelect_from_occ
1247 : !!***
1248 :
1249 : !----------------------------------------------------------------------
1250 :
1251 : !!****f* m_hdr/hdr_init_lowlvl
1252 : !! NAME
1253 : !! hdr_init_lowlvl
1254 : !!
1255 : !! FUNCTION
1256 : !! This subroutine initializes the header structured datatype
1257 : !! and most of its content from psps and other input variables that
1258 : !! are passed explicitly. It also use default values for evolving variables.
1259 : !! Note that Dtset is not required thus rendering the initialization of the header much easier.
1260 : !!
1261 : !! INPUTS
1262 : !! ebands <type(ebands_t)>=band structure information including Brillouin zone description
1263 : !! codvsn=code version
1264 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
1265 : !! comm_atom=--optional-- MPI communicator over atoms
1266 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
1267 : !! pertcase=index of the perturbation, or 0 if GS calculation
1268 : !! psps <type(pseudopotential_type)>=all the information about psps
1269 : !! For the meaning of the other varialble see the definition of dataset_type.
1270 : !!
1271 : !! OUTPUT
1272 : !! hdr <type(hdr_type)>=the header, initialized, and for most part of
1273 : !! it, contain its definite values, except for evolving variables
1274 : !!
1275 : !! SOURCE
1276 :
1277 23086 : subroutine hdr_init_lowlvl(hdr,ebands,psps,pawtab,wvl,&
1278 : codvsn,pertcase,natom,nsym,nspden,ecut,pawecutdg,ecutsm,dilatmx,&
1279 11543 : intxc,ixc,stmbias,usewvl,pawcpxocc,pawspnorb,ngfft,ngfftdg,so_psp,qptn,&
1280 11543 : rprimd,xred,symrel,tnons,symafm,typat,amu,icoulomb,&
1281 : kptopt,nelect,ne_qFD,nh_qFD,ivalence,cellcharge,kptrlatt_orig,kptrlatt,&
1282 11543 : nshiftk_orig,nshiftk,shiftk_orig,shiftk,&
1283 11543 : mpi_atmtab,comm_atom) ! optional arguments (parallelism)
1284 :
1285 : !Arguments ------------------------------------
1286 : !scalars
1287 : class(hdr_type),intent(inout) :: hdr
1288 : class(ebands_base_t),intent(in) :: ebands
1289 : integer,intent(in) :: natom,nsym,nspden,intxc,ixc,usewvl,pawcpxocc,pawspnorb,pertcase
1290 : integer,intent(in) :: ivalence,kptopt,nshiftk_orig,nshiftk,icoulomb
1291 : integer, intent(in),optional :: comm_atom
1292 : real(dp),intent(in) :: ecut,ecutsm,dilatmx,stmbias,pawecutdg,nelect,ne_qFD,nh_qFD,cellcharge
1293 : character(len=8),intent(in) :: codvsn
1294 : type(pseudopotential_type),intent(in) :: psps
1295 : type(wvl_internal_type),intent(in) :: wvl
1296 : !arrays
1297 : integer,intent(in) :: typat(natom)
1298 : integer,intent(in) :: so_psp(psps%npsp)
1299 : integer,intent(in) :: symrel(3,3,nsym),symafm(nsym)
1300 : integer,intent(in) :: ngfft(18),ngfftdg(18),kptrlatt_orig(3,3),kptrlatt(3,3)
1301 : integer,optional,target,intent(in) :: mpi_atmtab(:)
1302 : real(dp),intent(in) :: tnons(3,nsym),amu(psps%ntypat)
1303 : real(dp),intent(in) :: qptn(3) ! the wavevector, in case of a perturbation
1304 : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
1305 : real(dp),intent(in) :: shiftk_orig(3,nshiftk_orig),shiftk(3,nshiftk)
1306 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
1307 :
1308 : !Local variables-------------------------------
1309 : !scalars
1310 : integer :: bantot,date,nkpt,npsp,ntypat,nsppol,nspinor
1311 : integer :: cplex_rhoij,nspden_rhoij,qphase_rhoij
1312 : integer :: idx,isppol,ikpt,iband,ipsp
1313 : character(len=8) :: date_time
1314 : ! *************************************************************************
1315 :
1316 : !@hdr_type
1317 11543 : call date_and_time(date_time)
1318 11543 : read(date_time,'(i8)')date
1319 :
1320 11543 : npsp = psps%npsp
1321 11543 : ntypat = psps%ntypat
1322 11543 : nkpt = ebands%nkpt
1323 11543 : nsppol = ebands%nsppol
1324 11543 : nspinor= ebands%nspinor
1325 11543 : bantot = ebands%bantot
1326 :
1327 : !Transfer dimensions and other scalars to hdr.
1328 11543 : hdr%intxc =intxc
1329 11543 : hdr%ixc =ixc
1330 11543 : hdr%natom =natom
1331 11543 : hdr%npsp =npsp
1332 11543 : hdr%nspden =nspden
1333 11543 : hdr%nspinor =nspinor
1334 11543 : hdr%nsym =nsym
1335 11543 : hdr%ntypat =ntypat
1336 11543 : hdr%bantot =bantot
1337 11543 : hdr%nkpt =nkpt
1338 11543 : hdr%nshiftk_orig = nshiftk_orig
1339 11543 : hdr%nshiftk = nshiftk
1340 11543 : hdr%nsppol =nsppol
1341 11543 : hdr%usepaw =psps%usepaw
1342 11543 : hdr%usewvl =usewvl !hdr%nwvlarr will be set later since the number !of wavelets have not yet been computed.
1343 11543 : hdr%occopt =ebands%occopt
1344 11543 : hdr%codvsn =codvsn
1345 11543 : hdr%date =date
1346 11543 : hdr%headform =HDR_LATEST_HEADFORM ! Initialize with the latest headform
1347 11543 : hdr%pertcase =pertcase
1348 11543 : hdr%ecut =ecut
1349 11543 : hdr%ecutsm =ecutsm
1350 11543 : hdr%ecut_eff =ecut * (dilatmx)**2
1351 11543 : hdr%stmbias =stmbias
1352 11543 : hdr%tphysel =ebands%tphysel
1353 11543 : hdr%tsmear =ebands%tsmear
1354 46172 : hdr%qptn =qptn
1355 150059 : hdr%rprimd =rprimd ! Evolving data
1356 :
1357 : !Default for other data (all evolving data)
1358 11543 : hdr%etot =1.0d20
1359 11543 : hdr%fermie =1.0d20
1360 11543 : hdr%fermih =1.0d20
1361 11543 : hdr%residm =1.0d20
1362 :
1363 : ! Allocate all components of hdr
1364 : call hdr_malloc(hdr, hdr%bantot, hdr%nkpt, hdr%nsppol, hdr%npsp, hdr%natom, hdr%ntypat,&
1365 11543 : hdr%nsym, hdr%nshiftk_orig, hdr%nshiftk)
1366 :
1367 : !Transfer data from ebands
1368 439380 : hdr%istwfk(1:nkpt) = ebands%istwfk(1:nkpt)
1369 1722891 : hdr%kptns(:,:) = ebands%kptns(:,:)
1370 880189 : hdr%nband(1:nkpt*nsppol)=ebands%nband(1:nkpt*nsppol); hdr%mband = maxval(hdr%nband)
1371 439380 : hdr%npwarr(:) = ebands%npwarr(:)
1372 439380 : hdr%wtk(:) = ebands%wtk(:)
1373 :
1374 : !Transfer data from psps
1375 39919 : hdr%pspcod =psps%pspcod
1376 39919 : hdr%pspdat =psps%pspdat
1377 39919 : hdr%pspso =psps%pspso
1378 39919 : hdr%pspxc =psps%pspxc
1379 39919 : hdr%znuclpsp =psps%znuclpsp
1380 39891 : hdr%znucltypat=psps%znucltypat
1381 39919 : hdr%zionpsp =psps%zionpsp
1382 28376 : do ipsp=1,psps%npsp
1383 28376 : write(hdr%title(ipsp), "(A)") psps%title(ipsp)(1:132)
1384 : end do
1385 39919 : hdr%md5_pseudos = psps%md5_pseudos
1386 :
1387 39919 : hdr%so_psp=so_psp
1388 245178 : hdr%symafm(1:min(size(symafm),size(hdr%symafm)))=symafm(1:min(size(symafm),size(hdr%symafm)))
1389 3048798 : hdr%symrel(:,:,1:min(size(symrel,3),size(hdr%symrel,3))) =symrel(:,:,1:min(size(symrel,3),size(hdr%symrel,3)))
1390 946083 : hdr%tnons(:,1:min(size(tnons,2),size(hdr%tnons,2)))=tnons(:,1:min(size(tnons,2),size(hdr%tnons,2)))
1391 :
1392 38354 : hdr%typat(1:natom) =typat(1:natom) ! PMA : in tests/v2/t11 size(dtset%typat) is bigger dtset%natom
1393 118787 : hdr%xred(:,1:natom)=xred(:,1:natom) ! Evolving data
1394 :
1395 11543 : hdr%kptopt = kptopt
1396 11543 : hdr%pawcpxocc = pawcpxocc
1397 11543 : hdr%nelect = nelect
1398 11543 : hdr%ne_qFD = ne_qFD
1399 11543 : hdr%nh_qFD = nh_qFD
1400 11543 : hdr%ivalence = ivalence
1401 11543 : hdr%cellcharge = cellcharge
1402 150059 : hdr%kptrlatt_orig = kptrlatt_orig
1403 150059 : hdr%kptrlatt = kptrlatt
1404 108078 : hdr%shiftk_orig = shiftk_orig(:, 1:hdr%nshiftk_orig)
1405 69258 : hdr%shiftk = shiftk
1406 11543 : hdr%icoulomb = icoulomb
1407 39891 : hdr%amu = amu
1408 :
1409 11543 : if (psps%usepaw==1)then
1410 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
1411 2239 : nspden=nspden,spnorb=pawspnorb,cpxocc=pawcpxocc,qpt=qptn)
1412 12780 : ABI_MALLOC(hdr%pawrhoij,(natom))
1413 : ! Values of nspden/nspinor/nsppol are dummy ones; they are overwritten later (by hdr_update)
1414 4478 : if (present(comm_atom)) then
1415 2225 : if (present(mpi_atmtab)) then
1416 : call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,nsppol,typat,qphase=qphase_rhoij,&
1417 738 : pawtab=pawtab,comm_atom=comm_atom,mpi_atmtab=mpi_atmtab)
1418 : else
1419 : call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,nsppol,typat,qphase=qphase_rhoij,&
1420 1487 : pawtab=pawtab,comm_atom=comm_atom)
1421 : end if
1422 : else
1423 : call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,nsppol,typat,qphase=qphase_rhoij,&
1424 14 : pawtab=pawtab)
1425 : end if
1426 : end if
1427 :
1428 11543 : if (psps%usepaw==1 .and. usewvl ==0 ) then
1429 8956 : hdr%ngfft(:) =ngfftdg(1:3)
1430 9304 : else if (usewvl==1) then
1431 : #if defined HAVE_BIGDFT
1432 : hdr%ngfft(:) = (/ wvl%Glr%d%n1i, wvl%Glr%d%n2i, wvl%Glr%d%n3i /)
1433 : #else
1434 0 : BIGDFT_NOTENABLED_ERROR()
1435 : #endif
1436 : else
1437 37216 : hdr%ngfft(:) =ngfft(1:3)
1438 : end if
1439 :
1440 : !Transfer paw data
1441 11543 : if(psps%usepaw==1) then
1442 2239 : hdr%ecutdg =pawecutdg
1443 5766 : hdr%lmn_size(1:npsp)=pawtab(1:npsp)%lmn_size
1444 : else
1445 9304 : hdr%ecutdg=hdr%ecut
1446 22610 : hdr%lmn_size(:)=psps%lmnmax
1447 : end if
1448 :
1449 2576792 : hdr%occ(:)=zero; idx=0
1450 24036 : do isppol=1,nsppol
1451 458359 : do ikpt=1,nkpt
1452 3012065 : do iband=1,hdr%nband(ikpt+(isppol-1)*nkpt)
1453 2565249 : idx=idx+1
1454 2999572 : hdr%occ(idx)=ebands%occ(iband,ikpt,isppol)
1455 : end do
1456 : end do
1457 : end do
1458 :
1459 : ABI_UNUSED(wvl%h(1))
1460 :
1461 11543 : end subroutine hdr_init_lowlvl
1462 : !!***
1463 :
1464 : !----------------------------------------------------------------------
1465 :
1466 : !!****f* m_hdr/hdr_read_from_fname
1467 : !! NAME
1468 : !! hdr_read_from_fname
1469 : !!
1470 : !! FUNCTION
1471 : !! Read the header from file fname. Use Fortran IO or Netcdf depending on the extension of the file
1472 : !! Only rank0 process reads the header and then broadcast data to the other processes inside comm.
1473 : !!
1474 : !! INPUTS
1475 : !! fname=String with the name of the file.
1476 : !! comm = MPI communicator.
1477 : !!
1478 : !! OUTPUT
1479 : !! fform=Kind of the array in the file (0 signals an error)
1480 : !!
1481 : !! SOURCE
1482 :
1483 9100 : subroutine hdr_read_from_fname(Hdr, fname, fform, comm)
1484 :
1485 : !Arguments ------------------------------------
1486 : class(hdr_type),intent(out) :: Hdr
1487 : character(len=*),intent(in) :: fname
1488 : integer,intent(out) :: fform
1489 : integer,intent(in) :: comm
1490 :
1491 : !Local variables-------------------------------
1492 : integer,parameter :: rdwr1=1,master=0
1493 : integer :: fh,my_rank,mpierr
1494 : character(len=500) :: msg
1495 9100 : character(len=len(fname)) :: my_fname
1496 : ! *************************************************************************
1497 :
1498 9100 : my_rank = xmpi_comm_rank(comm)
1499 :
1500 9100 : if (my_rank == master) then
1501 9068 : my_fname = fname
1502 9068 : if (nctk_try_fort_or_ncfile(my_fname, msg) /= 0) then
1503 0 : ABI_ERROR(msg)
1504 : end if
1505 :
1506 9068 : if (.not. isncfile(my_fname)) then
1507 : ! Use Fortran IO to open the file and read the header.
1508 0 : if (open_file(my_fname, msg, newunit=fh, form="unformatted", status="old") /= 0) then
1509 0 : ABI_ERROR(msg)
1510 : end if
1511 :
1512 0 : call hdr%fort_read(fh,fform,rewind=(rdwr1==1))
1513 0 : ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname))
1514 0 : close(fh)
1515 :
1516 : else
1517 : ! Use Netcdf to open the file and read the header.
1518 9068 : NCF_CHECK(nctk_open_read(fh, my_fname, xmpi_comm_self))
1519 9068 : call hdr%ncread(fh, fform)
1520 9068 : ABI_CHECK(fform /= 0, sjoin("Error while reading:", my_fname))
1521 9068 : NCF_CHECK(nf90_close(fh))
1522 : end if
1523 : end if ! my_rank
1524 :
1525 : ! Broadcast fform and the header.
1526 9100 : if (xmpi_comm_size(comm) > 1) then
1527 48 : call hdr%bcast(master, my_rank, comm)
1528 48 : call xmpi_bcast(fform, master, comm, mpierr)
1529 : end if
1530 :
1531 9100 : end subroutine hdr_read_from_fname
1532 : !!***
1533 :
1534 : !----------------------------------------------------------------------
1535 :
1536 : !!****f* m_hdr/hdr_write_to_fname
1537 : !! NAME
1538 : !! hdr_write_to_fname
1539 : !!
1540 : !! FUNCTION
1541 : !! Write the header and fform to file fname.
1542 : !! Use Fortran IO or Netcdf depending on the extension of the file
1543 : !!
1544 : !! INPUTS
1545 : !! fname=String with the name of the file.
1546 : !! fform=Kind of the array in the file
1547 : !! Hdr<hdr_type>=The abinit header.
1548 : !!
1549 : !! OUTPUT
1550 : !! Only writing.
1551 : !!
1552 : !! SOURCE
1553 :
1554 72 : subroutine hdr_write_to_fname(Hdr, fname, fform)
1555 :
1556 : !Arguments ------------------------------------
1557 : class(hdr_type),intent(inout) :: Hdr
1558 : character(len=*),intent(in) :: fname
1559 : integer,intent(in) :: fform
1560 :
1561 : !Local variables-------------------------------
1562 : integer :: fh,ierr
1563 : character(len=500) :: msg
1564 : ! *************************************************************************
1565 :
1566 72 : if (.not.isncfile(fname)) then
1567 : ! Use Fortran IO to write the header.
1568 72 : if (open_file(fname,msg,newunit=fh,form="unformatted", status="unknown") /= 0) then
1569 0 : ABI_ERROR(msg)
1570 : end if
1571 72 : call hdr%fort_write(fh, fform, ierr)
1572 72 : ABI_CHECK(ierr==0, sjoin("Error while writing Abinit header to file:", fname))
1573 72 : close(fh)
1574 :
1575 : else
1576 : ! Use Netcdf to open the file and write the header.
1577 0 : if (file_exists(fname)) then
1578 0 : NCF_CHECK(nctk_open_modify(fh,fname, xmpi_comm_self))
1579 : else
1580 0 : NCF_CHECK_MSG(nctk_open_create(fh, fname, xmpi_comm_self), sjoin("Creating file:", fname))
1581 : end if
1582 :
1583 0 : NCF_CHECK(hdr%ncwrite(fh, fform, nc_define=.True.))
1584 0 : NCF_CHECK(nf90_close(fh))
1585 : end if
1586 :
1587 72 : end subroutine hdr_write_to_fname
1588 : !!***
1589 :
1590 : !----------------------------------------------------------------------
1591 :
1592 : !!****f* m_hdr/hdr_mpio_skip
1593 : !! NAME
1594 : !! hdr_mio_skip
1595 : !!
1596 : !! FUNCTION
1597 : !! Skip the abinit header in MPI-IO mode. This routine uses local MPI-IO calls hence
1598 : !! it can be safely called by master node only. Note however that in this case the
1599 : !! offset has to be communicated to the other nodes.
1600 : !!
1601 : !! INPUTS
1602 : !! mpio_fh=MPI-IO file handler
1603 : !!
1604 : !! TODO
1605 : !! We don't need to read record to skip. We just need to compute the offset from the dimensions.
1606 : !! The algorithm is as follows:
1607 : !!
1608 : !! 1) master reads and broadcast the header.
1609 : !! 2) The offset is computed from the header
1610 : !! 3) Open the file with MPI and use the offset to point the data to be read.
1611 : !!
1612 : !! OUTPUT
1613 : !! fform=kind of the array in the file
1614 : !! offset=The offset of the Fortran record located immediately below the Abinit header.
1615 : !!
1616 : !! SOURCE
1617 :
1618 139 : subroutine hdr_mpio_skip(mpio_fh, fform, offset)
1619 :
1620 : !Arguments ------------------------------------
1621 : integer,intent(in) :: mpio_fh
1622 : integer,intent(out) :: fform
1623 : integer(kind=XMPI_OFFSET_KIND),intent(out) :: offset
1624 :
1625 : !Local variables-------------------------------
1626 : !scalars
1627 : integer :: bsize_frm,mpi_type_frm
1628 : #ifdef HAVE_MPI_IO
1629 : integer :: headform,ierr,mu,usepaw,npsp
1630 : !arrays
1631 : integer(kind=MPI_OFFSET_KIND) :: fmarker,positloc
1632 : integer :: iread(1),statux(MPI_STATUS_SIZE)
1633 : #endif
1634 : character(len=500) :: msg
1635 : ! *************************************************************************
1636 :
1637 : !@hdr_type
1638 139 : offset = 0; fform = 0
1639 :
1640 139 : bsize_frm = xmpio_bsize_frm ! bsize_frm= Byte length of the Fortran record marker.
1641 139 : mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker.
1642 :
1643 : #ifdef HAVE_MPI_IO
1644 : !Reading the first record of the file -------------------------------------
1645 : !read (unitfi) codvsn,headform,..............
1646 139 : positloc = bsize_frm + 8*xmpi_bsize_ch
1647 139 : call MPI_FILE_READ_AT(mpio_fh,positloc,iread,1,MPI_INTEGER,statux,ierr)
1648 139 : fform=iread(1)
1649 :
1650 973 : if (ANY(fform == [1,2,51,52,101,102] )) then
1651 : ! This is the old format !read (unitfi) codvsn,fform
1652 0 : headform=22
1653 : write(msg,'(3a,i0,4a)') &
1654 0 : "ABINIT version: ",trim(abinit_version)," cannot read old files with headform: ",headform,ch10,&
1655 0 : "produced by previous versions. Use an old ABINIT version to read this file or ",ch10,&
1656 0 : "regenerate your files with version >= 8.0."
1657 0 : ABI_ERROR(msg)
1658 :
1659 : else
1660 : !read (unitfi)codvsn,headform,fform
1661 139 : call MPI_FILE_READ_AT(mpio_fh,positloc,iread,1,MPI_INTEGER,statux,ierr)
1662 139 : headform=iread(1)
1663 139 : positloc = positloc + xmpi_bsize_int
1664 139 : call MPI_FILE_READ_AT(mpio_fh,positloc,iread,1,MPI_INTEGER,statux,ierr)
1665 139 : fform=iread(1)
1666 : end if
1667 :
1668 139 : if (headform < 80) then
1669 : write(msg,'(3a,i0,4a)') &
1670 0 : "ABINIT version: ",trim(abinit_version)," cannot read old files with headform: ",headform,ch10,&
1671 0 : "produced by previous versions. Use an old ABINIT version to read this file or ",ch10,&
1672 0 : "regenerate your files with version >= 8.0."
1673 0 : ABI_ERROR(msg)
1674 : end if
1675 :
1676 : ! Skip first record.
1677 139 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1678 :
1679 : !Read npsp and usepaw from the second record and skip it
1680 139 : positloc = offset + bsize_frm + xmpi_bsize_int*13
1681 139 : call MPI_FILE_READ_AT(mpio_fh,positloc,iread,1,MPI_INTEGER,statux,ierr)
1682 139 : npsp=iread(1)
1683 139 : positloc = positloc + xmpi_bsize_int*4
1684 139 : call MPI_FILE_READ_AT(mpio_fh,positloc,iread,1,MPI_INTEGER,statux,ierr)
1685 139 : usepaw=iread(1)
1686 139 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1687 :
1688 : ! Skip the rest of the file ---------------------------------------------
1689 723 : do mu=1,3+npsp
1690 723 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1691 : end do
1692 :
1693 139 : if (usepaw == 1) then ! skip rhoij records.
1694 90 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1695 90 : call xmpio_read_frm(mpio_fh,offset,xmpio_single,fmarker,ierr)
1696 : end if
1697 :
1698 : #else
1699 : ABI_ERROR("hdr_mpio_skip cannot be used when MPI-IO is not enabled")
1700 : ABI_UNUSED(mpio_fh)
1701 : #endif
1702 :
1703 139 : end subroutine hdr_mpio_skip
1704 : !!***
1705 :
1706 : !----------------------------------------------------------------------
1707 :
1708 : !!****f* m_hdr/hdr_bsize_frecords
1709 : !! NAME
1710 : !! hdr_bsize_frecords
1711 : !!
1712 : !! FUNCTION
1713 : !! Compute the size of the Fortran records of the WFK file from the header and formeig.
1714 : !!
1715 : !! INPUTS
1716 : !! Hdr<hdr_type>=The abinit header.
1717 : !! formeig = 0 for GS WFK, 1 for response function WFK.
1718 : !!
1719 : !! OUTPUTS
1720 : !! nfrec = Number fof Fortran records
1721 : !! bsize_frecords(nfrec) = Byte size of each records. Allocated inside this routine.
1722 : !!
1723 : !! SOURCE
1724 :
1725 5 : subroutine hdr_bsize_frecords(Hdr, formeig, nfrec, bsize_frecords)
1726 :
1727 : !Arguments ------------------------------------
1728 : !scalars
1729 : class(hdr_type),intent(in) :: Hdr
1730 : integer,intent(in) :: formeig
1731 : integer,intent(out) :: nfrec
1732 : !arrays
1733 : integer(XMPI_OFFSET_KIND),allocatable,intent(out) :: bsize_frecords(:)
1734 :
1735 : !Local variables-------------------------------
1736 : !scalars
1737 : integer :: max_nfrec,ik_ibz,spin,mband,nband_k,npw_k,band
1738 : !arrays
1739 5 : integer(XMPI_OFFSET_KIND),allocatable :: bsz_frec(:)
1740 : !************************************************************************
1741 :
1742 : !@hdr_type
1743 35 : mband = MAXVAL(Hdr%nband)
1744 5 : max_nfrec = Hdr%nkpt*Hdr%nsppol * (3 + mband)
1745 :
1746 5 : if (formeig==1) max_nfrec = max_nfrec + Hdr%nkpt*Hdr%nsppol*mband
1747 15 : ABI_MALLOC(bsz_frec, (max_nfrec))
1748 :
1749 5 : nfrec = 0
1750 10 : do spin=1,Hdr%nsppol
1751 40 : do ik_ibz=1,Hdr%nkpt
1752 30 : nband_k = Hdr%nband(ik_ibz + (spin-1)*Hdr%nkpt)
1753 30 : npw_k = Hdr%npwarr(ik_ibz)
1754 :
1755 : ! First record: npw, nspinor, nband_disk
1756 30 : nfrec = nfrec + 1
1757 30 : bsz_frec(nfrec) = 3*xmpi_bsize_int
1758 :
1759 : ! Record with kg_k(3,npw_k) vectors
1760 30 : nfrec = nfrec + 1
1761 30 : bsz_frec(nfrec) = 3*npw_k*xmpi_bsize_int
1762 :
1763 35 : if (formeig==0) then
1764 : ! Record with the eigenvalues
1765 : ! eig_k(nband_k), occ_k(nband_k)
1766 30 : nfrec = nfrec + 1
1767 30 : bsz_frec(nfrec) = 2*nband_k*xmpi_bsize_dp
1768 :
1769 : ! cg_k record
1770 270 : do band=1,nband_k
1771 240 : nfrec = nfrec + 1
1772 270 : bsz_frec(nfrec) = 2*npw_k*Hdr%nspinor*xmpi_bsize_dp
1773 : end do
1774 :
1775 0 : else if (formeig==1) then
1776 0 : do band=1,nband_k
1777 : ! Record with the eigenvalues
1778 0 : nfrec = nfrec + 1
1779 0 : bsz_frec(nfrec) = 2*nband_k*xmpi_bsize_dp
1780 :
1781 : ! cg_k record
1782 0 : nfrec = nfrec + 1
1783 0 : bsz_frec(nfrec) = 2*npw_k*Hdr%nspinor*xmpi_bsize_dp
1784 : end do
1785 : else
1786 0 : ABI_ERROR("Wrong formeig")
1787 : end if
1788 :
1789 : end do
1790 : end do
1791 :
1792 15 : ABI_MALLOC(bsize_frecords, (nfrec))
1793 340 : bsize_frecords = bsz_frec(1:nfrec)
1794 :
1795 5 : ABI_FREE(bsz_frec)
1796 :
1797 5 : end subroutine hdr_bsize_frecords
1798 : !!***
1799 :
1800 : !----------------------------------------------------------------------
1801 :
1802 : !!****f* m_hdr/hdr_io_wfftype
1803 : !! NAME
1804 : !! hdr_io_wfftype
1805 : !!
1806 : !! FUNCTION
1807 : !! This subroutine deals with the I/O of the hdr_type
1808 : !! structured variables (read/write/echo).
1809 : !! According to the value of rdwr, it reads the header
1810 : !! of a file, writes it, or echo the value of the structured
1811 : !! variable to a file.
1812 : !! Note that, when reading, different records of hdr
1813 : !! are allocated here, according to the values of the
1814 : !! read variables. Records of hdr should be deallocated
1815 : !! correctly by a call to hdr_free when hdr is not used anymore.
1816 : !! Two instances of the hdr_io routines are defined :
1817 : !! hdr_io_int to which only the unit number is given
1818 : !! hdr_io_wfftype to which a wffil datatype is given
1819 : !!
1820 : !! INPUTS
1821 : !! rdwr= if 1, read the hdr structured variable from the header of the file,
1822 : !! if 2, write the header to unformatted file
1823 : !! if 3, echo part of the header to formatted file (records 1 and 2)
1824 : !! if 4, echo the header to formatted file
1825 : !! if 5, read the hdr without rewinding (unformatted)
1826 : !! if 6, write the hdr without rewinding (unformatted)
1827 : !! unitfi=unit number of the file (unformatted if rdwr=1, 2, 5 or 6 formatted if rdwr=3,4)
1828 : !!
1829 : !! OUTPUT
1830 : !! (see side effects)
1831 : !!
1832 : !! SIDE EFFECTS
1833 : !! The following variables are both input or output :
1834 : !! fform=kind of the array in the file
1835 : !! if rdwr=1,5 : will be output ; if the reading fail, return fform=0
1836 : !! if rdwr=2,3,4,6 : should be input, will be written or echo to file
1837 : !! hdr <type(hdr_type)>=the header structured variable
1838 : !! if rdwr=1,5 : will be output
1839 : !! if rdwr=2,3,4,6 : should be input, will be written or echo to file
1840 : !!
1841 : !! NOTES
1842 : !! In all cases, the file is supposed to be open already
1843 : !! When reading (rdwr=1) or writing (rdwr=2), rewind the file
1844 : !! When echoing (rdwr=3) does not rewind the file.
1845 : !! When reading (rdwr=5) or writing (rdwr=6), DOES NOT rewind the file
1846 : !!
1847 : !! SOURCE
1848 :
1849 164 : subroutine hdr_io_wfftype(fform,hdr,rdwr,wff)
1850 :
1851 : !Arguments ------------------------------------
1852 : integer,intent(inout) :: fform
1853 : integer,intent(in) :: rdwr
1854 : type(hdr_type),intent(inout) :: hdr
1855 : type(wffile_type),intent(inout) :: wff
1856 :
1857 : !Local variables-------------------------------
1858 : #if defined HAVE_MPI
1859 : integer :: ierr
1860 : #endif
1861 : ! *************************************************************************
1862 :
1863 : DBG_ENTER("COLL")
1864 :
1865 : if ( wff%iomode==IO_MODE_FORTRAN .or. &
1866 164 : (wff%iomode==IO_MODE_FORTRAN_MASTER .and.wff%master==wff%me).or. &
1867 : (wff%iomode==IO_MODE_MPI .and.wff%master==wff%me) ) then
1868 139 : call hdr_io_int(fform,hdr,rdwr,wff%unwff)
1869 : ! Master node **MUST** flush the output buffer so that the
1870 : ! other nodes can read headform and therefore the Fortran marker length when MPI-IO is used
1871 139 : if (rdwr == 2) call flush_unit(wff%unwff)
1872 : end if
1873 :
1874 : #if defined HAVE_MPI
1875 : !In the parallel case, if the files were not local, need to bcast the data
1876 164 : if(rdwr==1)then
1877 4 : if (wff%iomode==IO_MODE_FORTRAN_MASTER .or. wff%iomode==IO_MODE_MPI) then
1878 4 : if (wff%spaceComm/=MPI_COMM_SELF) then
1879 4 : call xmpi_bcast(fform,wff%master,wff%spaceComm,ierr)
1880 4 : call hdr%bcast(wff%master, wff%me, wff%spaceComm)
1881 : end if
1882 4 : wff%headform=hdr%headform
1883 4 : if(wff%iomode==IO_MODE_MPI)then
1884 4 : call hdr_skip_wfftype(wff,ierr)
1885 : end if
1886 : end if
1887 : end if
1888 : #if defined HAVE_MPI_IO
1889 164 : if (rdwr == 2 .and. wff%iomode==IO_MODE_MPI) then
1890 28 : if (wff%spaceComm/=MPI_COMM_SELF) then
1891 28 : call xmpi_barrier(wff%spaceComm)
1892 : end if
1893 28 : wff%headform=hdr%headform
1894 28 : call hdr_skip_wfftype(wff,ierr)
1895 : end if
1896 : #endif
1897 164 : if (rdwr==5) wff%headform=hdr%headform
1898 : #else
1899 : if (rdwr==1.or.rdwr==5) wff%headform=hdr%headform
1900 : #endif
1901 :
1902 : DBG_EXIT("COLL")
1903 :
1904 164 : end subroutine hdr_io_wfftype
1905 : !!***
1906 :
1907 : !----------------------------------------------------------------------
1908 :
1909 : !!****f* m_hdr/hdr_io_int
1910 : !! NAME
1911 : !! hdr_io_int
1912 : !!
1913 : !! FUNCTION
1914 : !! This subroutine deals with the I/O of the hdr_type structured variables (read/write/echo).
1915 : !! According to the value of rdwr, it reads the header of a file, writes it, or echo the value of the structured
1916 : !! variable to a file. Note that, when reading, different records of hdr are allocated here, according to the values of the
1917 : !! read variables. Records of hdr should be deallocated correctly by a call to hdr_free when hdr is not used anymore.
1918 : !! Two instances of the hdr_io routines are defined :
1919 : !! hdr_io_int to which only the unit number is given
1920 : !! hdr_io_wfftype to which a wffil datatype is given
1921 : !!
1922 : !! INPUTS
1923 : !! rdwr= if 1, read the hdr structured variable from the header of the file,
1924 : !! if 2, write the header to unformatted file
1925 : !! if 3, echo part of the header to formatted file (records 1 and 2)
1926 : !! if 4, echo the header to formatted file
1927 : !! if 5, read the hdr without rewinding (unformatted)
1928 : !! if 6, write the hdr without rewinding (unformatted)
1929 : !! unitfi=unit number of the file (unformatted if rdwr=1, 2, 5 or 6 formatted if rdwr=3,4)
1930 : !!
1931 : !! OUTPUT
1932 : !! (see side effects)
1933 : !!
1934 : !! SIDE EFFECTS
1935 : !! The following variables are both input or output :
1936 : !! fform=kind of the array in the file
1937 : !! if rdwr=1,5 : will be output ; if the reading fail, return fform=0
1938 : !! if rdwr=2,3,4,6 : should be input, will be written or echo to file
1939 : !! hdr <type(hdr_type)>=the header structured variable
1940 : !! if rdwr=1,5 : will be output
1941 : !! if rdwr=2,3,4,6 : should be input, will be written or echo to file
1942 : !!
1943 : !! NOTES
1944 : !! In all cases, the file is supposed to be open already
1945 : !! When reading (rdwr=1) or writing (rdwr=2), rewind the file
1946 : !! When echoing (rdwr=3) does not rewind the file.
1947 : !! When reading (rdwr=5) or writing (rdwr=6), DOES NOT rewind the file
1948 : !!
1949 : !! SOURCE
1950 :
1951 139 : subroutine hdr_io_int(fform,hdr,rdwr,unitfi)
1952 :
1953 : !Arguments ------------------------------------
1954 : integer,intent(inout) :: fform
1955 : integer,intent(in) :: rdwr,unitfi
1956 : type(hdr_type),intent(inout) :: hdr
1957 :
1958 : !Local variables-------------------------------
1959 : integer :: ierr
1960 : !*************************************************************************
1961 :
1962 : DBG_ENTER("COLL")
1963 :
1964 140 : select case(rdwr)
1965 : case (1, 5)
1966 : ! Reading the header of an unformatted file
1967 1 : call hdr%fort_read(unitfi,fform,rewind=(rdwr==1))
1968 :
1969 : case (2, 6)
1970 : ! Writing the header of an unformatted file
1971 138 : call hdr%fort_write(unitfi, fform, ierr, rewind=(rdwr==2))
1972 :
1973 : case (3, 4)
1974 : ! Writing the header of a formatted file
1975 0 : call hdr%echo(fform,rdwr,unit=unitfi)
1976 : case default
1977 139 : ABI_ERROR(sjoin("Wrong value for rdwr: ",itoa(rdwr)))
1978 : end select
1979 :
1980 : DBG_EXIT("COLL")
1981 :
1982 139 : end subroutine hdr_io_int
1983 : !!***
1984 : !----------------------------------------------------------------------
1985 :
1986 : !!****f* m_hdr/hdr_echo
1987 : !! NAME
1988 : !! hdr_echo
1989 : !!
1990 : !! FUNCTION
1991 : !! Echo the header
1992 : !!
1993 : !! INPUTS
1994 : !! hdr <type(hdr_type)>=the header structured variable
1995 : !! rdwr= if 3, echo part of the header to formatted file (records 1 and 2)
1996 : !! if 4, echo the header to formatted file
1997 : !! fform=kind of the array in the file
1998 : !! [unit]=unit number of the formatted file [DEFAULT: std_out]
1999 : !! [header]=Optional title.
2000 : !!
2001 : !! OUTPUT
2002 : !! Only writing
2003 : !!
2004 : !! TODO
2005 : !! Activate new header, avoid printing tons of lines with occupations.
2006 : !!
2007 : !! SOURCE
2008 :
2009 119 : subroutine hdr_echo(hdr, fform, rdwr, unit, header)
2010 :
2011 : !Arguments ------------------------------------
2012 : integer,intent(inout) :: fform
2013 : integer,intent(in) :: rdwr
2014 : integer,optional,intent(in) :: unit
2015 : class(hdr_type),intent(inout) :: hdr
2016 : character(len=*),optional,intent(in) :: header
2017 :
2018 : !Local variables-------------------------------
2019 : integer,parameter :: max_ns=6
2020 : integer :: iatom,ii,ikpt,ipsp,isym,ount !,ns
2021 : !character(len=500) :: msg
2022 : !*************************************************************************
2023 :
2024 119 : ount = std_out; if (present(unit)) ount = unit; if (ount == dev_null) return
2025 :
2026 119 : write(ount,'(a)')' ==============================================================================='
2027 119 : if (present(header)) write(ount, "(a)")ch10//' === '//trim(adjustl(header))//' === '
2028 119 : if (rdwr==3) write(ount, '(a)' ) ' ECHO of part of the ABINIT file header '
2029 119 : if (rdwr==4) write(ount, '(a)' ) ' ECHO of the ABINIT file header '
2030 119 : write(ount, '(a)' ) ' '
2031 119 : write(ount, '(a)' ) ' First record :'
2032 119 : write(ount, '(a,a8,2i5)' ) '.codvsn,headform,fform = ',hdr%codvsn, hdr%headform, fform
2033 119 : write(ount, '(a)' ) ' '
2034 119 : write(ount, '(a)' ) ' Second record :'
2035 119 : write(ount, '(a,4i6)') ' bantot,intxc,ixc,natom =',hdr%bantot, hdr%intxc, hdr%ixc, hdr%natom
2036 119 : write(ount, '(a,4i6)') ' ngfft(1:3),nkpt =',hdr%ngfft(1:3), hdr%nkpt
2037 119 : write(ount, '(a,2i6)') ' nspden,nspinor =',hdr%nspden, hdr%nspinor
2038 119 : write(ount, '(a,4i6)' ) ' nsppol,nsym,npsp,ntypat =',hdr%nsppol,hdr%nsym,hdr%npsp,hdr%ntypat
2039 119 : write(ount, '(a,3i6)' ) ' occopt,pertcase,usepaw =',hdr%occopt,hdr%pertcase,hdr%usepaw
2040 119 : write(ount, '(a,3es18.10)') ' ecut,ecutdg,ecutsm =',hdr%ecut, hdr%ecutdg, hdr%ecutsm
2041 119 : write(ount, '(a, es18.10)' ) ' ecut_eff =',hdr%ecut_eff
2042 119 : write(ount, '(a,3es18.10)') ' qptn(1:3) =',hdr%qptn(1:3)
2043 119 : write(ount, '(a,3es18.10)' ) ' rprimd(1:3,1) =',hdr%rprimd(1:3,1)
2044 119 : write(ount, '(a,3es18.10)' ) ' rprimd(1:3,2) =',hdr%rprimd(1:3,2)
2045 119 : write(ount, '(a,3es18.10)' ) ' rprimd(1:3,3) =',hdr%rprimd(1:3,3)
2046 119 : write(ount, '(a,3es18.10)') ' stmbias,tphysel,tsmear =',hdr%stmbias,hdr%tphysel, hdr%tsmear
2047 :
2048 : #ifdef DEV_NEW_HDR
2049 : write(ount, "(a,2es18.10,i0)") ' nelect,cellcharge,icoulomb =',hdr%nelect, hdr%cellcharge, hdr%icoulomb
2050 : write(ount, "(a,2i6)") ' kptopt,pawcpxocc =',hdr%kptopt, hdr%pawcpxocc
2051 : write(ount, '(a,9(i0,1x))') ' kptrlatt_orig = ',hdr%kptrlatt_orig
2052 : write(ount, '(a,9(i0,1x))' ) ' kptrlatt = ',hdr%kptrlatt
2053 :
2054 : ns = min(size(hdr%shiftk_orig, dim=2), max_ns)
2055 : write(msg, sjoin("(a,",itoa(3*ns),"(f4.2,1x))")) ' shiftk_orig = ',hdr%shiftk_orig(:,1:ns)
2056 : if (size(hdr%shiftk_orig, dim=2) > max_ns) msg = sjoin(msg, "...")
2057 : write(ount,"(a)")trim(msg)
2058 :
2059 : ns = min(size(hdr%shiftk, dim=2), max_ns)
2060 : write(msg, sjoin("(a,",itoa(3*ns),"(f4.2,1x))")) ' shiftk = ',hdr%shiftk(:,1:ns)
2061 : if (size(hdr%shiftk, dim=2) > max_ns) msg = sjoin(msg, "...")
2062 : write(ount,"(a)")trim(msg)
2063 : #endif
2064 :
2065 119 : write(ount, '(a)' )
2066 119 : if (rdwr==3)then
2067 0 : write(ount, '(a,i3,a)' ) ' The header contain ',hdr%npsp+2,' additional records.'
2068 : else
2069 119 : write(ount, '(a)' ) ' Third record :'
2070 119 : write(ount, '(a,(12i5,8x))') ' istwfk=',hdr%istwfk
2071 119 : write(ount, '(a,(12i5,8x))') ' nband =',hdr%nband
2072 119 : write(ount, '(a,(10i5,8x))') ' npwarr=',hdr%npwarr
2073 :
2074 119 : write(ount, '(a,(12i4,8x))') ' so_psp=',hdr%so_psp(:)
2075 : !write(ount,'(a,(12f6.2,1x))' )' amu =',hdr%amu
2076 :
2077 119 : write(ount, '(a)') ' symafm='
2078 119 : write(ount, '(8x,24i3,8x)') hdr%symafm
2079 :
2080 119 : write(ount, '(a)' ) ' symrel='
2081 2622 : do isym=1,hdr%nsym/2
2082 2622 : write(ount, '(a,9i4,a,9i4)' ) ' ',hdr%symrel(:,:,2*isym-1),' ',hdr%symrel(:,:,2*isym)
2083 : end do
2084 119 : if(2*(hdr%nsym/2)/=hdr%nsym)write(ount, '(a,9i4)' ) ' ',hdr%symrel(:,:,hdr%nsym)
2085 :
2086 119 : write(ount, '(a,(12i4,8x))') ' type =',hdr%typat(:)
2087 119 : write(ount, '(a)' ) ' kptns = (max 50 k-points will be written)'
2088 1552 : do ikpt=1,min(hdr%nkpt,50)
2089 1552 : write(ount, '(a,3es16.6)' ) ' ',hdr%kptns(:,ikpt)
2090 : end do
2091 119 : write(ount, '(a)' ) ' wtk ='
2092 386 : do ikpt=1,hdr%nkpt,10
2093 386 : write(ount, '(a,10f6.2)' ) ' ',hdr%wtk(ikpt:min(hdr%nkpt,ikpt + 10 - 1))
2094 : end do
2095 119 : write(ount, '(a)' ) ' occ ='
2096 3084 : do ii=1,hdr%bantot,10
2097 3084 : write(ount, '(a,10f6.2)') ' ',hdr%occ(ii:min(hdr%bantot,ii+10-1))
2098 : end do
2099 119 : write(ount, '(a)' ) ' tnons ='
2100 2622 : do isym=1,hdr%nsym/2
2101 2622 : write(ount, '(a,3f10.6,a,3f10.6)' ) ' ',hdr%tnons(:,2*isym-1),' ',hdr%tnons(:,2*isym)
2102 : end do
2103 119 : if(2*(hdr%nsym/2)/=hdr%nsym)write(ount, '(a,3f10.6)' ) ' ',hdr%tnons(:,hdr%nsym)
2104 119 : write(ount, '(a,(10f6.2,8x))') ' znucl=',hdr%znucltypat(:)
2105 119 : write(ount,'(a)')
2106 :
2107 119 : write(ount, '(a)' ) ' Pseudopotential info :'
2108 262 : do ipsp=1,hdr%npsp
2109 143 : write(ount,'(a,a)' ) ' title=',trim(hdr%title(ipsp))
2110 : ! TODO: This part should always be printed.
2111 : !write(ount,'(a,a)' ) ' md5=',trim(hdr%md5_pseudos(ipsp))
2112 : write(ount,'(a,f6.2,a,f6.2,a,i3,a,i6,a,i3,a,i3)' ) &
2113 143 : ' znuclpsp=',hdr%znuclpsp(ipsp), ', zionpsp=', hdr%zionpsp(ipsp),&
2114 143 : ', pspso=' , hdr%pspso(ipsp), ', pspdat=',hdr%pspdat(ipsp), &
2115 286 : ', pspcod=', hdr%pspcod(ipsp), ', pspxc=', hdr%pspxc(ipsp)
2116 :
2117 262 : if(hdr%usepaw==1)then
2118 2 : write(ount,'(a,i3)' ) ' lmn_size=', hdr%lmn_size(ipsp)
2119 : else
2120 141 : write(ount,'(a,i3)' ) ' lmnmax =', hdr%lmn_size(ipsp)
2121 : end if
2122 : end do
2123 :
2124 119 : write(ount, '(a)' ) ' '
2125 119 : write(ount, '(a)' ) ' Last record :'
2126 119 : write(ount, '(a,es16.6,es22.12,es16.6)' )' residm,etot,fermie=',hdr%residm, hdr%etot, hdr%fermie
2127 119 : write(ount, '(a)' ) ' xred ='
2128 380 : do iatom=1,hdr%natom
2129 380 : write(ount, '(a,3es16.6)' ) ' ',hdr%xred(:,iatom)
2130 : end do
2131 :
2132 119 : if (hdr%usepaw==1)then
2133 2 : call pawrhoij_io(hdr%pawrhoij,ount,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Echo")
2134 : end if
2135 :
2136 : if (rdwr==3)write(ount, '(a)' ) ' End the ECHO of part of the ABINIT file header '
2137 119 : if (rdwr==4)write(ount, '(a)' ) ' End the ECHO of the ABINIT file header '
2138 119 : write(ount,'(a)')' ==============================================================================='
2139 : end if ! rdwr is 3 or 4
2140 :
2141 119 : call flush_unit(ount)
2142 :
2143 119 : end subroutine hdr_echo
2144 : !!***
2145 :
2146 : !----------------------------------------------------------------------
2147 :
2148 : !!****f* m_hdr/hdr_skip_int
2149 : !! NAME
2150 : !! hdr_skip_int
2151 : !!
2152 : !! FUNCTION
2153 : !! Skip wavefunction or density file header, after having rewound the file.
2154 : !! Two instances of the hdr_skip routines are defined:
2155 : !! hdr_skip_int to which only the unit number is given
2156 : !! hdr_skip_wfftype to which a wffil datatype is given
2157 : !!
2158 : !! INPUTS
2159 : !! unit = number of unit to be read
2160 : !!
2161 : !! OUTPUT
2162 : !! ierr = error code returned by the MPI calls
2163 : !!
2164 : !! SIDE EFFECTS
2165 : !!
2166 : !! NOTES
2167 : !! No checking performed, since hdr_skip is assumed to be used only
2168 : !! on temporary wavefunction files.
2169 : !! This initialize further reading and checking by rwwf
2170 : !!
2171 : !! SOURCE
2172 :
2173 22 : subroutine hdr_skip_int(unitfi,ierr)
2174 :
2175 : !Arguments ------------------------------------
2176 : integer,intent(in) :: unitfi
2177 : integer,intent(out) :: ierr
2178 :
2179 : !Local variables-------------------------------
2180 : type(wffile_type) :: wff
2181 : ! *************************************************************************
2182 :
2183 : !Use default values for wff
2184 22 : wff%unwff=unitfi; wff%iomode=IO_MODE_FORTRAN
2185 22 : wff%me=0; wff%master=0
2186 : !Then, transmit to hdr_skip_wfftype
2187 22 : call hdr_skip_wfftype(wff,ierr)
2188 :
2189 22 : end subroutine hdr_skip_int
2190 : !!***
2191 :
2192 : !----------------------------------------------------------------------
2193 :
2194 : !!****f* m_hdr/hdr_skip_wfftype
2195 : !! NAME
2196 : !! hdr_skip_wfftype
2197 : !!
2198 : !! FUNCTION
2199 : !! Skip wavefunction or density file header, after having rewound the file.
2200 : !! Two instances of the hdr_skip routines are defined :
2201 : !! hdr_skip_int to which only the unit number is given
2202 : !! hdr_skip_wfftype to which a wffil datatype is given
2203 : !!
2204 : !! INPUTS
2205 : !! unit = number of unit to be read
2206 : !!
2207 : !! OUTPUT
2208 : !! ierr = error code returned by the MPI calls
2209 : !!
2210 : !! NOTES
2211 : !! No checking performed, since hdr_skip is assumed to be used only
2212 : !! on temporary wavefunction files.
2213 : !! This initialize further reading and checking by rwwf
2214 : !!
2215 : !! SOURCE
2216 :
2217 54 : subroutine hdr_skip_wfftype(wff,ierr)
2218 :
2219 : !Arguments ------------------------------------
2220 : type(wffile_type),intent(inout) :: wff
2221 : integer, intent(out) :: ierr
2222 :
2223 : !Local variables-------------------------------
2224 : integer :: headform,mu,npsp,unit,usepaw,fform
2225 : integer :: integers(17)
2226 : character(len=8) :: codvsn
2227 : character(len=500) :: msg,errmsg
2228 : #if defined HAVE_MPI_IO
2229 : integer(kind=MPI_OFFSET_KIND) :: delim_record,posit,positloc,off(1)
2230 : integer :: iread(1),statux(MPI_STATUS_SIZE)
2231 : #endif
2232 : !*************************************************************************
2233 :
2234 : !@hdr_type
2235 54 : unit=wff%unwff; ierr=0
2236 :
2237 54 : if( wff%iomode==IO_MODE_FORTRAN .or. (wff%iomode==IO_MODE_FORTRAN_MASTER.and.wff%master==wff%me) ) then
2238 :
2239 22 : rewind(unit, err=10, iomsg=errmsg)
2240 :
2241 : ! Pick off headform from WF file. Support for pre-v9 (length of codvsn was changed from 6 to 8) is implemented.
2242 22 : ABI_CHECK(read_first_record(unit, codvsn, headform, fform, errmsg) == 0, errmsg)
2243 :
2244 : if (headform==1 .or. headform==2 .or. &
2245 : headform==51 .or. headform==52 .or. &
2246 22 : headform==101 .or. headform==102 ) headform=22
2247 :
2248 22 : if (headform < 80) then
2249 : write(msg,'(3a,i0,4a)')&
2250 0 : "ABINIT version: ",trim(abinit_version)," cannot read old files with headform: ",headform,ch10,&
2251 0 : "produced by previous versions. Use an old ABINIT version to read this file or ",ch10,&
2252 0 : "regenerate your files with version >= 8.0."
2253 0 : ABI_ERROR(msg)
2254 : end if
2255 :
2256 22 : read (unit, err=10, iomsg=errmsg) integers(1:13),npsp,integers(15:17),usepaw
2257 :
2258 : ! Skip rest of header records
2259 111 : do mu=1,3+npsp
2260 111 : read (unit, err=10, iomsg=errmsg)
2261 : end do
2262 :
2263 22 : if (usepaw==1) then
2264 2 : read (unit, err=10, iomsg=errmsg)
2265 2 : read (unit, err=10, iomsg=errmsg)
2266 : end if
2267 :
2268 : #if defined HAVE_MPI_IO
2269 32 : else if(wff%iomode==IO_MODE_MPI)then
2270 :
2271 32 : headform=wff%headform
2272 : if (headform==1 .or. headform==2 .or. &
2273 : headform==51 .or. headform==52 .or. &
2274 32 : headform==101 .or. headform==102) headform=22
2275 :
2276 32 : if (headform < 80) then
2277 : write(msg,'(3a,i0,4a)')&
2278 0 : "ABINIT version: ",trim(abinit_version)," cannot read old files with headform: ",headform,ch10,&
2279 0 : "produced by previous versions. Use an old ABINIT version to read this file or ",ch10,&
2280 0 : "regenerate your files with version >= 8.0."
2281 0 : ABI_ERROR(msg)
2282 : end if
2283 :
2284 : ! Causes all previous writes to be transferred to the storage device
2285 32 : call flush_unit(wff%unwff)
2286 32 : call MPI_FILE_SYNC(wff%fhwff,ierr)
2287 :
2288 : ! Check FORTRAN record marker length (only at first call)
2289 32 : if (wff%nbOct_recMarker<=0) then
2290 0 : call getRecordMarkerLength_wffile(wff)
2291 : end if
2292 :
2293 32 : if (wff%master==wff%me) then
2294 :
2295 : ! Reading the first record of the file -------------------------------------
2296 : ! read (unitfi) codvsn,headform,..............
2297 27 : posit = 0
2298 27 : call rwRecordMarker(1,posit,delim_record,wff,ierr)
2299 :
2300 : ! Reading the second record of the file ------------------------------------
2301 : ! read(unitfi) bantot, hdr%date, hdr%intxc.................
2302 : ! Pick off npsp and usepaw from WF file
2303 27 : positloc = posit + wff%nbOct_recMarker + wff%nbOct_int*13
2304 27 : call MPI_FILE_READ_AT(wff%fhwff,positloc,iread,1,MPI_INTEGER,statux,ierr)
2305 27 : npsp=iread(1)
2306 :
2307 : ! Read usepaw and skip the fortran record
2308 27 : positloc = positloc + wff%nbOct_int*4
2309 27 : call MPI_FILE_READ_AT(wff%fhwff,positloc,iread,1,MPI_INTEGER,statux,ierr)
2310 27 : usepaw=iread(1)
2311 27 : call rwRecordMarker(1,posit,delim_record,wff,ierr)
2312 :
2313 : ! Skip the rest of the file ---------------------------------------------
2314 146 : do mu=1,3+npsp
2315 146 : call rwRecordMarker(1,posit,delim_record,wff,ierr)
2316 : end do
2317 :
2318 27 : if (usepaw==1) then
2319 7 : call rwRecordMarker(1,posit,delim_record,wff,ierr)
2320 7 : call rwRecordMarker(1,posit,delim_record,wff,ierr)
2321 : end if
2322 :
2323 27 : wff%offwff=posit
2324 : end if
2325 :
2326 32 : if (wff%spaceComm/=MPI_COMM_SELF) then
2327 32 : off(1)=wff%offwff
2328 32 : call MPI_BCAST(off,1,wff%offset_mpi_type,wff%master,wff%spaceComm,ierr)
2329 32 : wff%offwff=off(1)
2330 : end if
2331 : #endif
2332 : end if
2333 :
2334 : ! Handle IO-error: write warning and let the caller handle the exception.
2335 54 : return
2336 0 : 10 ierr=1
2337 0 : ABI_WARNING(errmsg)
2338 :
2339 : end subroutine hdr_skip_wfftype
2340 : !!***
2341 :
2342 : !----------------------------------------------------------------------
2343 :
2344 : !!****f* m_hdr/hdr_update
2345 : !! NAME
2346 : !! hdr_update
2347 : !!
2348 : !! FUNCTION
2349 : !! This subroutine update the header structured datatype.
2350 : !! Most of its records had been initialized correctly, but some corresponds
2351 : !! to evolving variables, or change with the context (like fform),
2352 : !! This routine is to be called before writing the header
2353 : !! to a file, in order to have up-to-date information.
2354 : !!
2355 : !! INPUTS
2356 : !! bantot=total number of bands
2357 : !! etot=total energy (Hartree)
2358 : !! extfpmd_eshift=--optional-- extfpmd energy shift (Hartree)
2359 : !! fermie=Fermi energy (Hartree)
2360 : !! fermih=Fermi energy for holes (Hartree), useful when occopt = 9
2361 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
2362 : !! comm_atom=--optional-- MPI communicator over atoms
2363 : !! residm=maximal residual
2364 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
2365 : !! occ(bantot)=occupancies for each band and k point
2366 : !! pawrhoij(natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
2367 : !! xred(3,natom)= relative coords of atoms in unit cell (dimensionless)
2368 : !! amu(ntypat)=masses in atomic mass units for each kind of atom in cell.
2369 : !!
2370 : !! OUTPUT
2371 : !! hdr <type(hdr_type)>=the header, initialized, and for most part of
2372 : !! it, contain its definite values, except for evolving variables
2373 : !!
2374 : !! SOURCE
2375 25943 : subroutine hdr_update(hdr,bantot,etot,fermie,fermih,residm,rprimd,occ,pawrhoij,xred,amu, &
2376 25943 : comm_atom,extfpmd_eshift,mpi_atmtab) ! optional arguments (parallelism and extfpmd)
2377 :
2378 : !Arguments ------------------------------------
2379 : !scalars
2380 : integer,intent(in) :: bantot
2381 : integer,optional,intent(in) :: comm_atom
2382 : real(dp),intent(in) :: etot,fermie,fermih,residm
2383 : real(dp),optional,intent(in) :: extfpmd_eshift
2384 : class(hdr_type),intent(inout) :: hdr
2385 : !arrays
2386 : integer,optional,target,intent(in) :: mpi_atmtab(:)
2387 : real(dp),intent(in) :: occ(bantot),rprimd(3,3),xred(3,hdr%natom),amu(hdr%ntypat)
2388 : type(pawrhoij_type),intent(inout) :: pawrhoij(:)
2389 : ! *************************************************************************
2390 :
2391 : !@hdr_type
2392 : !Update of the "evolving" data
2393 25943 : hdr%etot =etot
2394 25943 : hdr%fermie =fermie
2395 25943 : hdr%fermih =fermih
2396 25943 : hdr%residm =residm
2397 337259 : hdr%rprimd(:,:)=rprimd(:,:)
2398 4070879 : hdr%occ(:) =occ(:)
2399 283075 : hdr%xred(:,:)=xred(:,:)
2400 62137 : hdr%amu(:) = amu
2401 25943 : if(present(extfpmd_eshift)) hdr%extfpmd_eshift = extfpmd_eshift
2402 :
2403 25943 : if (hdr%usepaw==1) then
2404 6598 : if (present(comm_atom)) then
2405 6522 : if (present(mpi_atmtab)) then
2406 1750 : call pawrhoij_copy(pawrhoij,hdr%pawrhoij,comm_atom=comm_atom,mpi_atmtab=mpi_atmtab)
2407 : else
2408 4772 : call pawrhoij_copy(pawrhoij,hdr%pawrhoij,comm_atom=comm_atom)
2409 : end if
2410 : else
2411 76 : call pawrhoij_copy(pawrhoij,hdr%pawrhoij)
2412 : end if
2413 : end if
2414 :
2415 25943 : end subroutine hdr_update
2416 : !!***
2417 :
2418 : !----------------------------------------------------------------------
2419 :
2420 : !!****f* m_hdr/hdr_bcast
2421 : !! NAME
2422 : !! hdr_bcast
2423 : !!
2424 : !! FUNCTION
2425 : !! This subroutine transmits the header structured datatype
2426 : !! initialized on one processor (or a group of processor),
2427 : !! to the other processors. It also allocate the needed
2428 : !! part of the header.
2429 : !!
2430 : !! INPUTS
2431 : !! master = id of the master process
2432 : !! me = id of the current process
2433 : !! comm = id of the space communicator handler
2434 : !!
2435 : !! OUTPUT
2436 : !! (no output)
2437 : !!
2438 : !! SIDE EFFECTS
2439 : !! hdr <type(hdr_type)>=the header. For the master, it is already
2440 : !! initialized entirely, while for the other procs, everything has
2441 : !! to be transmitted.
2442 : !!
2443 : !! NOTES
2444 : !! This routine is called only in the case of MPI version of the code.
2445 : !!
2446 : !! SOURCE
2447 :
2448 1439 : subroutine hdr_bcast(hdr, master, me, comm)
2449 :
2450 : !Arguments ------------------------------------
2451 : class(hdr_type),intent(inout) :: hdr
2452 : integer, intent(in) :: master,me,comm
2453 :
2454 : !Local variables-------------------------------
2455 : !scalars
2456 : integer :: bantot,cplex_rhoij,iatom,ierr,index,index2,ipsp,iq,iq0,ispden,list_size,list_size2
2457 : integer :: lmn2_size,natom,nkpt,npsp,nsel,nspden,nsppol,nsym,nrhoij,ntypat,qphase
2458 : character(len=fnlen) :: list_tmp
2459 : !arrays
2460 1439 : integer,allocatable :: list_int(:)
2461 1439 : real(dp),allocatable :: list_dpr(:)
2462 1439 : character(len=fnlen),allocatable :: list_char(:)
2463 : ! *************************************************************************
2464 :
2465 : !@hdr_type
2466 4070 : if (xmpi_comm_size(comm) == 1) return ! Nothing to do
2467 :
2468 : DBG_ENTER("COLL")
2469 :
2470 : !Transmit the integer scalars
2471 562 : list_size = 44
2472 562 : ABI_MALLOC(list_int,(list_size))
2473 562 : if (master==me)then
2474 183 : list_int(1)=hdr%bantot
2475 183 : list_int(2)=hdr%date
2476 183 : list_int(3)=hdr%headform
2477 183 : list_int(4)=hdr%intxc
2478 183 : list_int(5)=hdr%ixc
2479 183 : list_int(6)=hdr%natom
2480 183 : list_int(7)=hdr%nkpt
2481 183 : list_int(8)=hdr%npsp
2482 183 : list_int(9)=hdr%nspden
2483 183 : list_int(10)=hdr%nspinor
2484 183 : list_int(11)=hdr%nsppol
2485 183 : list_int(12)=hdr%nsym
2486 183 : list_int(13)=hdr%ntypat
2487 183 : list_int(14)=hdr%occopt
2488 183 : list_int(15)=hdr%pertcase
2489 183 : list_int(16)=hdr%usepaw
2490 732 : list_int(17:19)=hdr%ngfft(1:3)
2491 183 : list_int(20)=hdr%usewvl
2492 183 : list_int(21)=hdr%kptopt
2493 183 : list_int(22)=hdr%pawcpxocc
2494 183 : list_int(23)=hdr%nshiftk_orig
2495 183 : list_int(24)=hdr%nshiftk
2496 183 : list_int(25:33)=reshape(hdr%kptrlatt_orig, [3*3])
2497 183 : list_int(34:42)=reshape(hdr%kptrlatt, [3*3])
2498 183 : list_int(43)=hdr%icoulomb
2499 183 : list_int(44)=hdr%ivalence
2500 : end if
2501 :
2502 562 : call xmpi_bcast(list_int,master,comm,ierr)
2503 :
2504 562 : if(master/=me)then
2505 379 : hdr%bantot =list_int(1)
2506 379 : hdr%date =list_int(2)
2507 379 : hdr%headform=list_int(3)
2508 379 : hdr%intxc =list_int(4)
2509 379 : hdr%ixc =list_int(5)
2510 379 : hdr%natom =list_int(6)
2511 379 : hdr%nkpt =list_int(7)
2512 379 : hdr%npsp =list_int(8)
2513 379 : hdr%nspden =list_int(9)
2514 379 : hdr%nspinor =list_int(10)
2515 379 : hdr%nsppol =list_int(11)
2516 379 : hdr%nsym =list_int(12)
2517 379 : hdr%ntypat =list_int(13)
2518 379 : hdr%occopt =list_int(14)
2519 379 : hdr%pertcase=list_int(15)
2520 379 : hdr%usepaw =list_int(16)
2521 1516 : hdr%ngfft(1:3)=list_int(17:19)
2522 379 : hdr%usewvl =list_int(20)
2523 379 : hdr%kptopt = list_int(21)
2524 379 : hdr%pawcpxocc = list_int(22)
2525 379 : hdr%nshiftk_orig = list_int(23)
2526 379 : hdr%nshiftk = list_int(24)
2527 4927 : hdr%kptrlatt_orig = reshape(list_int(25:33), [3,3])
2528 4927 : hdr%kptrlatt = reshape(list_int(34:42), [3,3])
2529 379 : hdr%icoulomb = list_int(43)
2530 379 : hdr%ivalence = list_int(44)
2531 : end if
2532 562 : ABI_FREE(list_int)
2533 :
2534 562 : bantot=hdr%bantot
2535 562 : natom =hdr%natom
2536 562 : nkpt =hdr%nkpt
2537 562 : npsp =hdr%npsp
2538 562 : nspden=hdr%nspden
2539 562 : nsppol=hdr%nsppol
2540 562 : nsym =hdr%nsym
2541 562 : ntypat=hdr%ntypat
2542 :
2543 562 : if (master/=me) then
2544 : ! Allocate all components of hdr
2545 : call hdr_malloc(hdr, bantot, nkpt, nsppol, npsp, natom, ntypat,&
2546 379 : nsym, hdr%nshiftk_orig, hdr%nshiftk)
2547 : end if
2548 :
2549 : !Transmit the integer arrays
2550 562 : list_size=nkpt*(2+nsppol)+6*npsp+10*nsym+natom
2551 1686 : ABI_MALLOC(list_int,(list_size))
2552 562 : if (master==me)then
2553 2405 : list_int(1 :nkpt )=hdr%istwfk ; index=nkpt
2554 2468 : list_int(1+index:nkpt*nsppol+index)=hdr%nband ; index=index+nkpt*nsppol
2555 2405 : list_int(1+index:nkpt +index)=hdr%npwarr ; index=index+nkpt
2556 445 : list_int(1+index:npsp +index)=hdr%pspcod ; index=index+npsp
2557 445 : list_int(1+index:npsp +index)=hdr%pspdat ; index=index+npsp
2558 445 : list_int(1+index:npsp +index)=hdr%pspso ; index=index+npsp
2559 445 : list_int(1+index:npsp +index)=hdr%pspxc ; index=index+npsp
2560 445 : list_int(1+index:npsp +index)=hdr%lmn_size ; index=index+npsp
2561 445 : list_int(1+index:npsp +index)=hdr%so_psp ; index=index+npsp
2562 5909 : list_int(1+index:nsym +index)=hdr%symafm ; index=index+nsym
2563 366 : list_int(1+index:nsym*3*3 +index)=reshape(hdr%symrel,(/3*3*nsym/))
2564 183 : index=index+nsym*3*3
2565 555 : list_int(1+index:natom +index)=hdr%typat ; index=index+natom
2566 : end if
2567 :
2568 562 : call xmpi_bcast(list_int,master,comm,ierr)
2569 :
2570 562 : if(master/=me)then
2571 5446 : hdr%istwfk=list_int(1 :nkpt ) ; index=nkpt
2572 5633 : hdr%nband =list_int(1+index:nkpt*nsppol+index) ; index=index+nkpt*nsppol
2573 5446 : hdr%npwarr=list_int(1+index:nkpt +index) ; index=index+nkpt
2574 1298 : hdr%pspcod=list_int(1+index:npsp +index) ; index=index+npsp
2575 1298 : hdr%pspdat=list_int(1+index:npsp +index) ; index=index+npsp
2576 1298 : hdr%pspso =list_int(1+index:npsp +index) ; index=index+npsp
2577 1298 : hdr%pspxc =list_int(1+index:npsp +index) ; index=index+npsp
2578 1298 : hdr%lmn_size=list_int(1+index:npsp +index) ; index=index+npsp
2579 1298 : hdr%so_psp =list_int(1+index:npsp +index) ; index=index+npsp
2580 12714 : hdr%symafm=list_int(1+index:nsym +index) ; index=index+nsym
2581 157323 : hdr%symrel=reshape(list_int(1+index:nsym*3*3 +index),(/3,3,nsym/))
2582 379 : index=index+nsym*3*3
2583 1530 : hdr%typat =list_int(1+index:natom +index) ; index=index+natom
2584 : end if
2585 562 : ABI_FREE(list_int)
2586 :
2587 : !Transmit the double precision scalars and arrays
2588 : list_size = 22+ 3*nkpt+nkpt+bantot + 3*nsym + 3*natom + 2*npsp+ntypat + &
2589 562 : 5 + 3*hdr%nshiftk_orig + 3*hdr%nshiftk + hdr%ntypat
2590 1686 : ABI_MALLOC(list_dpr,(list_size))
2591 :
2592 562 : if (master==me)then
2593 183 : list_dpr(1)=hdr%ecut_eff
2594 183 : list_dpr(2)=hdr%etot
2595 183 : list_dpr(3)=hdr%fermie
2596 183 : list_dpr(4)=hdr%fermih
2597 183 : list_dpr(5)=hdr%residm
2598 183 : list_dpr(6:14)=reshape(hdr%rprimd(1:3,1:3),(/9/))
2599 183 : list_dpr(15)=hdr%ecut
2600 183 : list_dpr(16)=hdr%ecutdg
2601 183 : list_dpr(17)=hdr%ecutsm
2602 183 : list_dpr(18)=hdr%tphysel
2603 183 : list_dpr(19)=hdr%tsmear
2604 732 : list_dpr(20:22)=hdr%qptn(1:3) ; index=22
2605 366 : list_dpr(1+index:3*nkpt +index)=reshape(hdr%kptns,(/3*nkpt/)) ; index=index+3*nkpt
2606 2405 : list_dpr(1+index:nkpt +index)=hdr%wtk ; index=index+nkpt
2607 28187 : list_dpr(1+index:bantot +index)=hdr%occ ; index=index+bantot
2608 366 : list_dpr(1+index:3*nsym +index)=reshape(hdr%tnons,(/3*nsym/)) ; index=index+3*nsym
2609 366 : list_dpr(1+index:3*natom+index)=reshape(hdr%xred,(/3*natom/)) ; index=index+3*natom
2610 445 : list_dpr(1+index:npsp +index)=hdr%zionpsp ; index=index+npsp
2611 445 : list_dpr(1+index:npsp +index)=hdr%znuclpsp ; index=index+npsp
2612 445 : list_dpr(1+index:ntypat +index)=hdr%znucltypat ; index=index+ntypat
2613 183 : list_dpr(1+index)=hdr%nelect; index=index+1
2614 183 : list_dpr(1+index)=hdr%ne_qFD; index=index+1
2615 183 : list_dpr(1+index)=hdr%nh_qFD; index=index+1
2616 183 : list_dpr(1+index)=hdr%cellcharge; index=index+1
2617 183 : list_dpr(1+index)=hdr%extfpmd_eshift; index=index+1
2618 366 : list_dpr(1+index:index+3*hdr%nshiftk_orig) = reshape(hdr%shiftk_orig, [3*hdr%nshiftk_orig])
2619 183 : index=index+3*hdr%nshiftk_orig
2620 366 : list_dpr(1+index:index+3*hdr%nshiftk) = reshape(hdr%shiftk, [3*hdr%nshiftk])
2621 183 : index=index+3*hdr%nshiftk
2622 445 : list_dpr(1+index:index+hdr%ntypat) = hdr%amu(1:hdr%ntypat)
2623 : end if
2624 :
2625 562 : call xmpi_bcast(list_dpr,master,comm,ierr)
2626 :
2627 562 : if(master/=me)then
2628 379 : hdr%ecut_eff=list_dpr(1)
2629 379 : hdr%etot =list_dpr(2)
2630 379 : hdr%fermie =list_dpr(3)
2631 379 : hdr%fermih =list_dpr(4)
2632 379 : hdr%residm =list_dpr(5)
2633 4927 : hdr%rprimd =reshape(list_dpr(6:14),(/3,3/))
2634 379 : hdr%ecut =list_dpr(15)
2635 379 : hdr%ecutdg =list_dpr(16)
2636 379 : hdr%ecutsm =list_dpr(17)
2637 379 : hdr%tphysel =list_dpr(18)
2638 379 : hdr%tsmear =list_dpr(19)
2639 1516 : hdr%qptn(1:3)=list_dpr(20:22) ; index=22
2640 20268 : hdr%kptns =reshape(list_dpr(1+index:3*nkpt +index),(/3,nkpt/)) ; index=index+3*nkpt
2641 5446 : hdr%wtk =list_dpr(1+index:nkpt +index) ; index=index+nkpt
2642 61070 : hdr%occ =list_dpr(1+index:bantot +index) ; index=index+bantot
2643 49340 : hdr%tnons =reshape(list_dpr(1+index:3*nsym +index),(/3,nsym/)) ; index=index+3*nsym
2644 4604 : hdr%xred =reshape(list_dpr(1+index:3*natom+index),(/3,natom/)); index=index+3*natom
2645 1298 : hdr%zionpsp =list_dpr(1+index:npsp +index) ; index=index+npsp
2646 1298 : hdr%znuclpsp=list_dpr(1+index:npsp +index) ; index=index+npsp
2647 1298 : hdr%znucltypat=list_dpr(1+index:ntypat +index) ; index=index+ntypat
2648 379 : hdr%nelect = list_dpr(1+index); index=index+1
2649 379 : hdr%ne_qFD = list_dpr(1+index); index=index+1
2650 379 : hdr%nh_qFD = list_dpr(1+index); index=index+1
2651 379 : hdr%cellcharge = list_dpr(1+index); index=index+1
2652 379 : hdr%extfpmd_eshift = list_dpr(1+index); index=index+1
2653 5672 : hdr%shiftk_orig = reshape(list_dpr(1+index:index+3*hdr%nshiftk_orig), [3, hdr%nshiftk_orig])
2654 379 : index=index+3*hdr%nshiftk_orig
2655 3032 : hdr%shiftk = reshape(list_dpr(1+index:index+3*hdr%nshiftk), [3, hdr%nshiftk])
2656 379 : index=index+3*hdr%nshiftk
2657 1298 : hdr%amu = list_dpr(1+index:index+hdr%ntypat)
2658 : end if
2659 562 : ABI_FREE(list_dpr)
2660 :
2661 : !Transmit the characters
2662 562 : list_size=npsp+1 + npsp
2663 1686 : ABI_MALLOC(list_char,(list_size))
2664 562 : if (master==me)then
2665 183 : list_char(1) =hdr%codvsn ! Only 8 characters are stored in list_char(1)
2666 445 : list_char(2:npsp+1)=hdr%title
2667 445 : list_char(npsp+2:) =hdr%md5_pseudos
2668 : end if
2669 :
2670 562 : call xmpi_bcast(list_char,master,comm,ierr)
2671 :
2672 562 : if(master/=me)then
2673 379 : list_tmp=list_char(1)
2674 379 : hdr%codvsn=list_tmp(1:8)
2675 919 : do ipsp=2,npsp+1
2676 540 : list_tmp =list_char(ipsp)
2677 919 : hdr%title(ipsp-1) =list_tmp(1:min(fnlen,132))
2678 : end do
2679 919 : do ipsp=npsp+2,2*npsp+1
2680 919 : hdr%md5_pseudos(ipsp-npsp-1) = list_char(ipsp)(1:md5_slen)
2681 : end do
2682 : end if
2683 562 : ABI_FREE(list_char)
2684 :
2685 : !Transmit the structured variables in case of PAW
2686 562 : if (hdr%usepaw==1) then
2687 :
2688 48 : nrhoij=0
2689 48 : if (master==me)then
2690 13 : cplex_rhoij=hdr%pawrhoij(1)%cplex_rhoij
2691 13 : qphase=hdr%pawrhoij(1)%qphase
2692 13 : nspden=hdr%pawrhoij(1)%nspden
2693 43 : do iatom=1,natom
2694 43 : nrhoij=nrhoij+hdr%pawrhoij(iatom)%nrhoijsel
2695 : end do
2696 : end if
2697 :
2698 48 : call xmpi_bcast(nrhoij,master,comm,ierr)
2699 48 : call xmpi_bcast(cplex_rhoij,master,comm,ierr)
2700 48 : call xmpi_bcast(qphase,master,comm,ierr)
2701 48 : call xmpi_bcast(nspden,master,comm,ierr)
2702 :
2703 48 : list_size=natom+nrhoij;list_size2=nspden*nrhoij*cplex_rhoij*qphase
2704 144 : ABI_MALLOC(list_int,(list_size))
2705 144 : ABI_MALLOC(list_dpr,(list_size2))
2706 48 : if (master==me)then
2707 : index=0;index2=0
2708 43 : do iatom=1,natom
2709 30 : nsel=hdr%pawrhoij(iatom)%nrhoijsel
2710 30 : lmn2_size=hdr%pawrhoij(iatom)%lmn2_size
2711 30 : list_int(1+index)=nsel
2712 582 : list_int(2+index:1+nsel+index)=hdr%pawrhoij(iatom)%rhoijselect(1:nsel)
2713 30 : index=index+1+nsel
2714 84 : do ispden=1,nspden
2715 112 : do iq=1,qphase
2716 41 : iq0=merge(0,lmn2_size*cplex_rhoij,iq==1)
2717 3026 : list_dpr(1+index2:nsel*cplex_rhoij+index2)=hdr%pawrhoij(iatom)%rhoijp(iq0+1:iq0+nsel*cplex_rhoij,ispden)
2718 82 : index2=index2+nsel*cplex_rhoij
2719 : end do
2720 : end do
2721 : end do
2722 : end if
2723 :
2724 48 : call xmpi_bcast(list_int,master,comm,ierr)
2725 48 : call xmpi_bcast(list_dpr,master,comm,ierr)
2726 :
2727 48 : if(master/=me)then
2728 35 : index=0;index2=0
2729 179 : ABI_MALLOC(hdr%pawrhoij,(natom))
2730 : call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden,hdr%nspinor,hdr%nsppol,hdr%typat,&
2731 35 : lmnsize=hdr%lmn_size,qphase=qphase)
2732 109 : do iatom=1,natom
2733 74 : nsel=list_int(1+index)
2734 74 : lmn2_size=hdr%pawrhoij(iatom)%lmn2_size
2735 74 : hdr%pawrhoij(iatom)%nrhoijsel=nsel
2736 1730 : hdr%pawrhoij(iatom)%rhoijselect(1:nsel)=list_int(2+index:1+nsel+index)
2737 74 : index=index+1+nsel
2738 216 : do ispden=1,nspden
2739 288 : do iq=1,qphase
2740 107 : iq0=merge(0,lmn2_size*cplex_rhoij,iq==1)
2741 9062 : hdr%pawrhoij(iatom)%rhoijp(iq0+1:iq0+nsel*cplex_rhoij,ispden)=list_dpr(1+index2:nsel*cplex_rhoij+index2)
2742 214 : index2=index2+nsel*cplex_rhoij
2743 : end do
2744 : end do
2745 : end do
2746 : end if
2747 48 : ABI_FREE(list_int)
2748 144 : ABI_FREE(list_dpr)
2749 : end if
2750 :
2751 7722 : hdr%mband = maxval(hdr%nband)
2752 :
2753 : DBG_EXIT("COLL")
2754 :
2755 5441 : end subroutine hdr_bcast
2756 : !!***
2757 :
2758 : !----------------------------------------------------------------------
2759 :
2760 : !!****f* m_hdr/read_first_record
2761 : !! NAME
2762 : !! read_first_record
2763 : !!
2764 : !! FUNCTION
2765 : !! Read the first record of the header.
2766 : !! This function is needed to support pre-Abinitv9 headers
2767 : !! since length of codvsn was changed from 6 to 8 in v9
2768 : !!
2769 : !! SOURCE
2770 :
2771 12504 : integer function read_first_record(unit, codvsn8, headform, fform, errmsg) result(ierr)
2772 :
2773 : !Arguments ------------------------------------
2774 : integer,intent(in) :: unit
2775 : integer,intent(out) :: headform, fform
2776 : character(len=8),intent(out) :: codvsn8
2777 : character(len=*),intent(out) :: errmsg
2778 :
2779 : !Local variables-------------------------------
2780 : integer :: major, ii
2781 : character(len=6) :: codvsn6
2782 : !*************************************************************************
2783 :
2784 : ! Try pre-v9 first. This read should not fail as we have enough space in the record
2785 : ! Obviously headform and fform are wrong in > Abinit9.
2786 12504 : read(unit, iostat=ierr, iomsg=errmsg) codvsn6, headform, fform
2787 12504 : if (ierr /= 0) then
2788 0 : call wrtout(std_out, "Fatal error while reading the first record of the Abinit header!")
2789 0 : return
2790 : end if
2791 :
2792 12504 : ii = index(codvsn6, ".")
2793 12504 : if (ii == 0 .or. ii == 1) then
2794 0 : errmsg = sjoin("Cannot find major.minor pattern in codvsn:", codvsn6)
2795 0 : ierr = 1; return
2796 : end if
2797 :
2798 12504 : major = atoi(codvsn6(:ii-1))
2799 : !call wrtout(std_out, sjoin("Reading HDR file generated by major version:", itoa(major)))
2800 12504 : if (major > 8) then
2801 12504 : backspace(unit)
2802 12504 : read(unit, iostat=ierr, iomsg=errmsg) codvsn8, headform, fform
2803 12504 : if (ierr /= 0) then
2804 0 : call wrtout(std_out, "Fatal error while reading the first record of the Abinit header version > 8!")
2805 0 : return
2806 : end if
2807 : else
2808 0 : codvsn8 = ""
2809 0 : codvsn8(1:6) = codvsn6
2810 : end if
2811 :
2812 25008 : end function read_first_record
2813 : !!***
2814 :
2815 : !----------------------------------------------------------------------
2816 :
2817 : !!****f* m_hdr/hdr_fort_read
2818 : !! NAME
2819 : !! hdr_fort_read
2820 : !!
2821 : !! FUNCTION
2822 : !! Reads the header from a logical unit associated to a unformatted file.
2823 : !! Note that, when reading, different records of hdr are allocated here, according to the values of the
2824 : !! read variables. Records of hdr should be deallocated correctly by a call to hdr_free when hdr is not used anymore.
2825 : !!
2826 : !! INPUTS
2827 : !! unit=unit number of the unformatted file
2828 : !! [rewind]=True to rewind the file. Default: False
2829 : !!
2830 : !! OUTPUT
2831 : !! Hdr<hdr_type>=The header of the file fully initialized (if fform /=0)
2832 : !! fform=kind of the array in the file. if the reading fail, return fform=0
2833 : !!
2834 : !! NOTES
2835 : !! The file is supposed to be open already
2836 : !!
2837 : !! SOURCE
2838 :
2839 12482 : subroutine hdr_fort_read(Hdr, unit, fform, rewind)
2840 :
2841 : !Arguments ------------------------------------
2842 : class(hdr_type),intent(out) :: hdr
2843 : integer,intent(out) :: fform
2844 : integer,intent(in) :: unit
2845 : logical,optional,intent(in) :: rewind
2846 :
2847 : !Local variables-------------------------------
2848 : !integer :: ierr
2849 : integer :: ipsp
2850 : character(len=500) :: msg,errmsg
2851 12482 : real(dp),allocatable :: occ3d(:,:,:)
2852 : !*************************************************************************
2853 :
2854 : !@hdr_type
2855 : DBG_ENTER("COLL")
2856 :
2857 12482 : if (present(rewind)) then
2858 52 : if (rewind) rewind(unit, err=10, iomsg=errmsg)
2859 : end if
2860 :
2861 : ! Reading the first record of the file ------------------------------------
2862 : ! fform is not a record of hdr_type
2863 12482 : ABI_CHECK(read_first_record(unit, hdr%codvsn, hdr%headform, fform, errmsg) == 0, errmsg)
2864 :
2865 : !write(std_out,*) 'In hdr_fort_read, l. 3032, headform, fform = ', hdr%headform, ', ', fform
2866 :
2867 12482 : if (hdr%headform < 80) then
2868 : write(msg,'(3a,i0,4a)') &
2869 0 : "ABINIT version: ",trim(abinit_version)," cannot read old files with headform: ",hdr%headform,ch10,&
2870 0 : "produced by previous versions. Use an old ABINIT version to read this file or ",ch10,&
2871 0 : "regenerate your files with version >= 8.0."
2872 0 : ABI_ERROR(msg)
2873 : end if
2874 :
2875 : call check_fform(fform)
2876 :
2877 : !Reading the second record of the file ------------------------------------
2878 : read(unit, err=10, iomsg=errmsg) &
2879 12482 : hdr%bantot, hdr%date, hdr%intxc, hdr%ixc, hdr%natom, hdr%ngfft(1:3),&
2880 12482 : hdr%nkpt, hdr%nspden, hdr%nspinor, hdr%nsppol, hdr%nsym, hdr%npsp, hdr%ntypat, hdr%occopt, hdr%pertcase,&
2881 12482 : hdr%usepaw, hdr%ecut, hdr%ecutdg, hdr%ecutsm, hdr%ecut_eff, hdr%qptn(1:3), hdr%rprimd,&
2882 24964 : hdr%stmbias, hdr%tphysel, hdr%tsmear, hdr%usewvl, hdr%nshiftk_orig, hdr%nshiftk, hdr%mband
2883 :
2884 : !Allocate all parts of hdr that need to be --------------------------------
2885 : call hdr_malloc(hdr, hdr%bantot, hdr%nkpt, hdr%nsppol, hdr%npsp, hdr%natom, hdr%ntypat,&
2886 12482 : hdr%nsym, hdr%nshiftk_orig, hdr%nshiftk)
2887 :
2888 12482 : if (hdr%usepaw==1) then
2889 29 : ABI_MALLOC(hdr%pawrhoij,(hdr%natom))
2890 : end if
2891 :
2892 : ! Reading the third record of the file ------------------------------------
2893 :
2894 : ! Take into account future migration to occ(:,:,:) in the Format
2895 : ! read 3d matrix with stride and transfer to (stupid) 1d hdr%occ in packed form.
2896 62410 : ABI_MALLOC(occ3d, (hdr%mband,hdr%nkpt,hdr%nsppol))
2897 :
2898 : read(unit, err=10, iomsg=errmsg) &
2899 12482 : hdr%istwfk(:), hdr%nband(:), hdr%npwarr(:), &
2900 12482 : hdr%so_psp(:), hdr%symafm(:), hdr%symrel(:,:,:), &
2901 12482 : hdr%typat(:), hdr%kptns(:,:), occ3d, &
2902 24964 : hdr%tnons(:,:), hdr%znucltypat(:), hdr%wtk(:)
2903 1059335 : ABI_CHECK(hdr%mband == maxval(hdr%nband), "mband != max(hdr%nband). Are you reading an Abinit8 file with Abinit9?")
2904 :
2905 12482 : call hdr_set_occ(hdr, occ3d)
2906 12482 : ABI_FREE(occ3d)
2907 :
2908 : ! Reading the final record of the header ---------------------------------
2909 12482 : read(unit, err=10, iomsg=errmsg) hdr%residm, hdr%xred(:,:), hdr%etot, hdr%fermie, hdr%amu(:)
2910 :
2911 : read(unit, err=10, iomsg=errmsg)&
2912 12482 : hdr%kptopt,hdr%pawcpxocc,hdr%nelect,hdr%cellcharge,hdr%icoulomb,&
2913 24964 : hdr%kptrlatt,hdr%kptrlatt_orig, hdr%shiftk_orig,hdr%shiftk
2914 :
2915 12482 : hdr%ivalence = hdr%nelect / 2 ! Read in case occopt = 9
2916 12482 : hdr%ne_qFD = zero
2917 12482 : hdr%nh_qFD = zero
2918 12482 : hdr%fermih = zero
2919 12482 : if (hdr%occopt == 9) then
2920 : ! This was erroneous
2921 : ! write(unit,err=10, iomsg=errmsg) hdr%ivalence, hdr%ne_qFD, hdr%nh_qFD, hdr%fermie, hdr%fermih
2922 : ! But this induced problems on bob_gnu_7.5_openmp
2923 0 : read(unit,err=10, iomsg=errmsg) hdr%ivalence, hdr%ne_qFD, hdr%nh_qFD, hdr%fermie, hdr%fermih
2924 : !
2925 : end if
2926 :
2927 : ! Reading the records with psp information ---------------------------------
2928 34353 : do ipsp=1,hdr%npsp
2929 : read(unit, err=10, iomsg=errmsg) &
2930 21871 : hdr%title(ipsp), hdr%znuclpsp(ipsp), hdr%zionpsp(ipsp), hdr%pspso(ipsp), hdr%pspdat(ipsp), &
2931 56224 : hdr%pspcod(ipsp), hdr%pspxc(ipsp), hdr%lmn_size(ipsp), hdr%md5_pseudos(ipsp)
2932 : end do
2933 :
2934 12482 : if (hdr%usepaw==1) then ! Reading the Rhoij tab if the PAW method was used.
2935 6 : call pawrhoij_io(hdr%pawrhoij,unit,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Read")
2936 : end if
2937 :
2938 : DBG_EXIT("COLL")
2939 : return
2940 :
2941 : ! Handle IO-error: write warning and let the caller handle the exception.
2942 0 : 10 fform=0
2943 0 : ABI_WARNING(errmsg)
2944 :
2945 0 : end subroutine hdr_fort_read
2946 : !!***
2947 :
2948 : !----------------------------------------------------------------------
2949 :
2950 : !!****f* m_hdr/hdr_ncread
2951 : !! NAME
2952 : !! hdr_ncread
2953 : !!
2954 : !! FUNCTION
2955 : !! This subroutine deals with the reading of the hdr_type structured variables
2956 : !! It handles variables according to the ETSF format, whenever
2957 : !! possible and uses new variables when not available in the ETSF format.
2958 : !! Note that, when reading, different records of hdr are allocated here,
2959 : !! Records of hdr should be deallocated
2960 : !! correctly by a call to hdr_free when hdr is not used anymore.
2961 : !!
2962 : !! INPUTS
2963 : !! ncid=the unit of the open NetCDF file.
2964 : !!
2965 : !! OUTPUT
2966 : !! fform=kind of the array in the file. if the reading fails, return fform=0
2967 : !!
2968 : !! SOURCE
2969 :
2970 14397 : subroutine hdr_ncread(Hdr, ncid, fform)
2971 :
2972 : !Arguments ------------------------------------
2973 : !scalars
2974 : class(hdr_type),target,intent(out) :: hdr
2975 : integer,intent(in) :: ncid
2976 : integer,intent(out) :: fform
2977 :
2978 : !Local variables-------------------------------
2979 : !scalars
2980 : integer :: nresolution, itypat, ii, varid, ncerr
2981 : character(len=500) :: msg
2982 : !arrays
2983 14397 : integer,allocatable :: nband2d(:,:)
2984 14397 : real(dp),allocatable :: occ3d(:,:,:)
2985 : ! *************************************************************************
2986 :
2987 : !@hdr_type
2988 14397 : NCF_CHECK(nctk_set_datamode(ncid))
2989 14397 : NCF_CHECK(nf90_get_var(ncid, vid("fform"), fform))
2990 14397 : NCF_CHECK(nf90_get_var(ncid, vid("headform"), hdr%headform))
2991 :
2992 14397 : if (hdr%headform < 80) then
2993 : write(msg,'(3a,i0,4a)')&
2994 0 : "ABINIT version: ",trim(abinit_version)," cannot read old files with headform: ",hdr%headform,ch10,&
2995 0 : "produced by previous versions. Use an old ABINIT version to read this file or ",ch10,&
2996 0 : "regenerate your files with version >= 8.0."
2997 0 : ABI_ERROR(msg)
2998 : end if
2999 :
3000 : call check_fform(fform)
3001 :
3002 : ! First, we read the declaration of code, fform ...
3003 : ! pad the returned string with " " instead of "\0"
3004 : !
3005 : ! Support for pre-v9 (length of codvsn was changed from 6 to 8)
3006 14397 : NCF_CHECK(nctk_get_dim(ncid, "codvsnlen", ii))
3007 14397 : NCF_CHECK(nf90_get_var(ncid, vid("codvsn"), hdr%codvsn(1:ii)))
3008 14397 : call replace_ch0(hdr%codvsn)
3009 :
3010 : ! Get ETSF dimensions
3011 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_atoms", hdr%natom))
3012 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_kpoints", hdr%nkpt))
3013 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_components", hdr%nspden))
3014 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_spinor_components", hdr%nspinor))
3015 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_spins", hdr%nsppol))
3016 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_symmetry_operations", hdr%nsym))
3017 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_atom_species", hdr%ntypat))
3018 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_grid_points_vector1", hdr%ngfft(1)))
3019 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_grid_points_vector2", hdr%ngfft(2)))
3020 14397 : NCF_CHECK(nctk_get_dim(ncid, "number_of_grid_points_vector3", hdr%ngfft(3)))
3021 14397 : NCF_CHECK(nctk_get_dim(ncid, "max_number_of_states", hdr%mband))
3022 : ! bantot is used to dimension %occ in hdr_malloc. Note that hdr%bantot != sum(nband) because states
3023 : ! are packed in hdr%occ and therefore bantot <= hdr%mband * hdr%nkpt * hdr%nsppol
3024 14397 : NCF_CHECK(nctk_get_dim(ncid, "bantot", hdr%bantot))
3025 :
3026 : ! Read other dimensions, not handled by ETSF format.
3027 14397 : NCF_CHECK(nctk_get_dim(ncid, "npsp", hdr%npsp))
3028 14397 : NCF_CHECK(nctk_get_dim(ncid, "nshiftk_orig", hdr%nshiftk_orig))
3029 14397 : NCF_CHECK(nctk_get_dim(ncid, "nshiftk", hdr%nshiftk))
3030 :
3031 : ! Read other important scalar variables.
3032 14397 : NCF_CHECK(nf90_get_var(ncid, vid("usepaw"), hdr%usepaw))
3033 14397 : NCF_CHECK(nf90_get_var(ncid, vid("usewvl"), hdr%usewvl))
3034 :
3035 14397 : nresolution=0
3036 14397 : if (hdr%usewvl == 1) then
3037 : ! This value must be 2...
3038 0 : NCF_CHECK(nctk_get_dim(ncid, "number_of_wavelet_resolutions", nresolution))
3039 : ! We set the right ngfft, adding the padding space for wavelets.
3040 0 : hdr%ngfft = hdr%ngfft + 31
3041 : end if
3042 :
3043 : ! Allocate all parts of hdr that need to be
3044 : call hdr_malloc(hdr, hdr%bantot, hdr%nkpt, hdr%nsppol, hdr%npsp, hdr%natom, hdr%ntypat,&
3045 14397 : hdr%nsym, hdr%nshiftk_orig, hdr%nshiftk)
3046 :
3047 57588 : ABI_MALLOC(nband2d, (hdr%nkpt, hdr%nsppol))
3048 14397 : NCF_CHECK(nf90_get_var(ncid, vid("number_of_states"), nband2d))
3049 671199 : hdr%nband(:) = reshape(nband2d, [hdr%nkpt*hdr%nsppol])
3050 14397 : ABI_FREE(nband2d)
3051 656802 : ABI_CHECK(hdr%mband == maxval(hdr%nband), "mband != maxval(hdr%nband)")
3052 :
3053 14397 : if (hdr%usepaw==1) then
3054 15999 : ABI_MALLOC(hdr%pawrhoij,(hdr%natom))
3055 : end if
3056 :
3057 : !We get then all variables included in ETSF
3058 14397 : if (hdr%usewvl==0) then
3059 14397 : NCF_CHECK(nf90_get_var(ncid, vid("kinetic_energy_cutoff"), hdr%ecut))
3060 14397 : NCF_CHECK(nf90_get_var(ncid, vid("number_of_coefficients"), hdr%npwarr))
3061 : else
3062 0 : NCF_CHECK(nf90_get_var(ncid, vid("number_of_wavelets"), hdr%nwvlarr))
3063 : end if
3064 :
3065 : ! read 3d matrix with stride and transfer to (stupid) 1d hdr%occ in packed form.
3066 4724651 : ABI_CALLOC(occ3d, (hdr%mband, hdr%nkpt, hdr%nsppol))
3067 14397 : NCF_CHECK(nf90_get_var(ncid, vid("occupations"), occ3d))
3068 14397 : call hdr_set_occ(hdr, occ3d)
3069 14397 : ABI_FREE(occ3d)
3070 :
3071 14397 : NCF_CHECK(nf90_get_var(ncid, vid("fermi_energy"), hdr%fermie))
3072 14397 : NCF_CHECK(nf90_get_var(ncid, vid("primitive_vectors"), hdr%rprimd))
3073 14397 : NCF_CHECK(nf90_get_var(ncid, vid("reduced_symmetry_matrices"), hdr%symrel))
3074 14397 : NCF_CHECK(nf90_get_var(ncid, vid("atom_species"), hdr%typat))
3075 14397 : NCF_CHECK(nf90_get_var(ncid, vid("reduced_symmetry_translations"), hdr%tnons))
3076 14397 : NCF_CHECK(nf90_get_var(ncid, vid("reduced_atom_positions"), hdr%xred))
3077 14397 : NCF_CHECK(nf90_get_var(ncid, vid("atomic_numbers"), hdr%znucltypat))
3078 14397 : NCF_CHECK(nf90_get_var(ncid, vid("reduced_coordinates_of_kpoints"), hdr%kptns))
3079 14397 : NCF_CHECK(nf90_get_var(ncid, vid("kpoint_weights"), hdr%wtk))
3080 14397 : NCF_CHECK(nf90_get_var(ncid, vid("date"), hdr%date))
3081 14397 : NCF_CHECK(nf90_get_var(ncid, vid("ecut_eff"), hdr%ecut_eff))
3082 14397 : NCF_CHECK(nf90_get_var(ncid, vid("ecutsm"), hdr%ecutsm))
3083 14397 : NCF_CHECK(nf90_get_var(ncid, vid("etot"), hdr%etot))
3084 14397 : NCF_CHECK(nf90_get_var(ncid, vid("intxc"), hdr%intxc))
3085 14397 : NCF_CHECK(nf90_get_var(ncid, vid("ixc"), hdr%ixc))
3086 14397 : NCF_CHECK(nf90_get_var(ncid, vid("occopt"), hdr%occopt))
3087 14397 : NCF_CHECK(nf90_get_var(ncid, vid("pertcase"), hdr%pertcase))
3088 14397 : NCF_CHECK(nf90_get_var(ncid, vid("qptn"), hdr%qptn))
3089 14397 : NCF_CHECK(nf90_get_var(ncid, vid("residm"), hdr%residm))
3090 14397 : NCF_CHECK(nf90_get_var(ncid, vid("stmbias"), hdr%stmbias))
3091 14397 : NCF_CHECK(nf90_get_var(ncid, vid("tphysel"), hdr%tphysel))
3092 14397 : NCF_CHECK(nf90_get_var(ncid, vid("tsmear"), hdr%tsmear))
3093 14397 : NCF_CHECK(nf90_get_var(ncid, vid("ecutdg"), hdr%ecutdg))
3094 :
3095 : ! Multidimensional variables. Be careful with zionpsp if alchemical mixing!
3096 14397 : NCF_CHECK(nf90_get_var(ncid, vid("istwfk"), hdr%istwfk))
3097 14397 : NCF_CHECK(nf90_get_var(ncid, vid("pspcod"), hdr%pspcod))
3098 14397 : NCF_CHECK(nf90_get_var(ncid, vid("pspdat"), hdr%pspdat))
3099 14397 : NCF_CHECK(nf90_get_var(ncid, vid("pspso"), hdr%pspso))
3100 14397 : NCF_CHECK(nf90_get_var(ncid, vid("pspxc"), hdr%pspxc))
3101 14397 : NCF_CHECK(nf90_get_var(ncid, vid("so_psp"), hdr%so_psp))
3102 14397 : NCF_CHECK(nf90_get_var(ncid, vid("symafm"), hdr%symafm))
3103 14397 : NCF_CHECK(nf90_get_var(ncid, vid("zionpsp"), hdr%zionpsp))
3104 14397 : NCF_CHECK(nf90_get_var(ncid, vid("znuclpsp"), hdr%znuclpsp))
3105 14397 : NCF_CHECK(nf90_get_var(ncid, vid("kptopt"), hdr%kptopt))
3106 14397 : NCF_CHECK(nf90_get_var(ncid, vid("pawcpxocc"), hdr%pawcpxocc))
3107 14397 : NCF_CHECK(nf90_get_var(ncid, vid("nelect"), hdr%nelect))
3108 14397 : NCF_CHECK(nf90_get_var(ncid, vid("charge"), hdr%cellcharge))
3109 14397 : NCF_CHECK(nf90_get_var(ncid, vid("kptrlatt_orig"), hdr%kptrlatt_orig))
3110 14397 : NCF_CHECK(nf90_get_var(ncid, vid("kptrlatt"), hdr%kptrlatt))
3111 14397 : NCF_CHECK(nf90_get_var(ncid, vid("shiftk_orig"), hdr%shiftk_orig))
3112 14397 : NCF_CHECK(nf90_get_var(ncid, vid("shiftk"), hdr%shiftk))
3113 14397 : NCF_CHECK(nf90_get_var(ncid, vid("md5_pseudos"), hdr%md5_pseudos))
3114 14397 : NCF_CHECK(nf90_get_var(ncid, vid("amu"), hdr%amu))
3115 14397 : NCF_CHECK(nf90_get_var(ncid, vid("icoulomb"), hdr%icoulomb))
3116 14397 : NCF_CHECK(nf90_get_var(ncid, vid("title"), hdr%title))
3117 :
3118 : ! Pad the returned string with " " instead of "\0"
3119 36030 : do itypat=1,size(hdr%title)
3120 36030 : call replace_ch0(hdr%title(itypat))
3121 : end do
3122 :
3123 14397 : NCF_CHECK(nf90_get_var(ncid, vid("lmn_size"), hdr%lmn_size))
3124 14397 : if (hdr%usepaw==1) then
3125 : call pawrhoij_io(hdr%pawrhoij,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,&
3126 2958 : hdr%headform,"Read",form="netcdf")
3127 : end if
3128 :
3129 : ! Reading the values of fermih, ne_qFD, nh_qFD and ivalence if occopt = 9
3130 14397 : hdr%fermih = zero
3131 14397 : hdr%ne_qFD = zero
3132 14397 : hdr%nh_qFD = zero
3133 14397 : hdr%ivalence = hdr%nelect / 2
3134 14397 : if (hdr%occopt == 9) then
3135 0 : ncerr = nf90_inq_varid(ncid, "hole_fermi_energy", varid)
3136 0 : if (ncerr /= nf90_noerr) then
3137 0 : NCF_CHECK(nf90_get_var(ncid, vid("hole_fermi_energy"), hdr%fermih))
3138 : end if
3139 0 : ncerr = nf90_inq_varid(ncid, "ne_qFD", varid)
3140 0 : if (ncerr /= nf90_noerr) then
3141 0 : NCF_CHECK(nf90_get_var(ncid, vid("ne_qFD"), hdr%ne_qFD))
3142 : end if
3143 0 : ncerr = nf90_inq_varid(ncid, "nh_qFD", varid)
3144 0 : if (ncerr /= nf90_noerr) then
3145 0 : NCF_CHECK(nf90_get_var(ncid, vid("nh_qFD"), hdr%nh_qFD))
3146 : end if
3147 0 : ncerr = nf90_inq_varid(ncid, "ivalence", varid)
3148 0 : if (ncerr /= nf90_noerr) then
3149 0 : NCF_CHECK(nf90_get_var(ncid, vid("ivalence"), hdr%ivalence))
3150 : end if
3151 :
3152 : endif
3153 :
3154 : contains
3155 777438 : integer function vid(vname)
3156 :
3157 : character(len=*),intent(in) :: vname
3158 777438 : vid = nctk_idname(ncid, vname)
3159 : end function vid
3160 :
3161 : end subroutine hdr_ncread
3162 : !!***
3163 :
3164 : !----------------------------------------------------------------------
3165 :
3166 : !!****f* m_hdr/hdr_fort_write
3167 : !! NAME
3168 : !! hdr_fort_write
3169 : !!
3170 : !! FUNCTION
3171 : !! Writes the header and fform to unformatted file
3172 : !!
3173 : !! INPUTS
3174 : !! Hdr<hdr_type>=The header of the file.
3175 : !! fform=kind of the array in the file
3176 : !! unit=unit number of the unformatted file
3177 : !! [rewind]=True to rewind the file. Default: False
3178 : !!
3179 : !! OUTPUT
3180 : !! ierr=Exit status
3181 : !!
3182 : !! NOTES
3183 : !! The file is supposed to be open already
3184 : !!
3185 : !! SOURCE
3186 :
3187 1110 : subroutine hdr_fort_write(Hdr,unit,fform,ierr,rewind)
3188 :
3189 : !Arguments ------------------------------------
3190 : integer,intent(out) :: ierr
3191 : integer,intent(in) :: unit,fform
3192 : logical,optional,intent(in) :: rewind
3193 : class(hdr_type),intent(inout) :: hdr
3194 :
3195 : !Local variables-------------------------------
3196 : integer :: headform,ipsp,major,ii
3197 : character(len=500) :: errmsg
3198 1110 : real(dp),allocatable :: occ3d(:,:,:)
3199 : !*************************************************************************
3200 :
3201 : ! TODO: Change intent to in. Change pawrhoij_io first!
3202 1110 : ierr = 0
3203 1110 : if (present(rewind)) then
3204 138 : if (rewind) rewind(unit, err=10, iomsg=errmsg)
3205 : end if
3206 :
3207 : call check_fform(fform)
3208 :
3209 1110 : ii = index(hdr%codvsn, ".")
3210 1110 : if (ii == 0 .or. ii == 1) then
3211 0 : ABI_WARNING(sjoin("Cannot find major.minor pattern in codvsn:", hdr%codvsn))
3212 1110 : ierr = 1; return
3213 : end if
3214 :
3215 1110 : major = atoi(hdr%codvsn(:ii-1))
3216 :
3217 : ! Writing always use last format version
3218 1110 : headform = HDR_LATEST_HEADFORM
3219 :
3220 1110 : if (major > 8) then
3221 1110 : write(unit, err=10, iomsg=errmsg) hdr%codvsn, headform, fform
3222 : else
3223 0 : write(unit, err=10, iomsg=errmsg) hdr%codvsn(1:6), headform, fform
3224 : end if
3225 :
3226 : write(unit, err=10, iomsg=errmsg) &
3227 1110 : hdr%bantot, hdr%date, hdr%intxc, hdr%ixc, hdr%natom, hdr%ngfft(1:3), &
3228 1110 : hdr%nkpt, hdr%nspden, hdr%nspinor, hdr%nsppol, hdr%nsym, hdr%npsp, hdr%ntypat, hdr%occopt, hdr%pertcase,&
3229 1110 : hdr%usepaw, hdr%ecut, hdr%ecutdg, hdr%ecutsm, hdr%ecut_eff, hdr%qptn, hdr%rprimd, &
3230 2220 : hdr%stmbias, hdr%tphysel, hdr%tsmear, hdr%usewvl, hdr%nshiftk_orig, hdr%nshiftk, hdr%mband
3231 109072 : ABI_CHECK(hdr%mband == maxval(hdr%nband), "mband != maxval(hdr%nband)")
3232 :
3233 5550 : ABI_MALLOC(occ3d, (hdr%mband,hdr%nkpt,hdr%nsppol))
3234 1110 : call hdr_get_occ3d(hdr, occ3d)
3235 1110 : write(unit,err=10, iomsg=errmsg) hdr%istwfk(:), hdr%nband(:), hdr%npwarr(:),&
3236 1110 : hdr%so_psp(:), hdr%symafm(:), hdr%symrel(:,:,:), hdr%typat(:), hdr%kptns(:,:), occ3d, &
3237 2220 : hdr%tnons(:,:), hdr%znucltypat(:), hdr%wtk(:)
3238 1110 : ABI_FREE(occ3d)
3239 :
3240 1110 : write(unit,err=10, iomsg=errmsg) hdr%residm, hdr%xred(:,:), hdr%etot, hdr%fermie, hdr%amu(:)
3241 : write(unit,err=10, iomsg=errmsg) &
3242 1110 : hdr%kptopt, hdr%pawcpxocc, hdr%nelect, hdr%cellcharge, hdr%icoulomb,&
3243 2220 : hdr%kptrlatt,hdr%kptrlatt_orig, hdr%shiftk_orig(:,1:hdr%nshiftk_orig),hdr%shiftk(:,1:hdr%nshiftk)
3244 :
3245 : ! Write record for occopt 9 option if needed
3246 1110 : if (hdr%occopt == 9) then
3247 0 : write(unit,err=10, iomsg=errmsg) hdr%ivalence, hdr%ne_qFD, hdr%nh_qFD, hdr%fermie, hdr%fermih
3248 : end if
3249 :
3250 : ! Write the records with psp information ---------------------------------
3251 2641 : do ipsp=1,hdr%npsp
3252 : write(unit, err=10, iomsg=errmsg) &
3253 1531 : hdr%title(ipsp), hdr%znuclpsp(ipsp), hdr%zionpsp(ipsp), hdr%pspso(ipsp), hdr%pspdat(ipsp), &
3254 4172 : hdr%pspcod(ipsp), hdr%pspxc(ipsp), hdr%lmn_size(ipsp), hdr%md5_pseudos(ipsp)
3255 : end do
3256 :
3257 1110 : if (hdr%usepaw==1) then
3258 150 : call pawrhoij_io(hdr%pawrhoij,unit,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,headform,"Write")
3259 : end if
3260 :
3261 : return
3262 :
3263 : ! Handle IO-error: write warning and let the caller handle the exception.
3264 0 : 10 ierr=1
3265 0 : ABI_WARNING(errmsg)
3266 :
3267 1110 : end subroutine hdr_fort_write
3268 : !!***
3269 :
3270 : !----------------------------------------------------------------------
3271 :
3272 : !!****f* m_hdr/hdr_backspace
3273 : !! NAME
3274 : !! hdr_backspace
3275 : !!
3276 : !! FUNCTION
3277 : !! Backspace the header. Return exit status and error message
3278 : !! The file is supposed to be open already
3279 : !!
3280 : !! INPUTS
3281 : !! Hdr<hdr_type>=The header of the file.
3282 : !! unit=unit number of the unformatted file
3283 : !!
3284 : !! SOURCE
3285 :
3286 2974 : integer function hdr_backspace(hdr, unit, msg) result(ierr)
3287 :
3288 : !Arguments ------------------------------------
3289 : class(hdr_type),intent(in) :: hdr
3290 : integer,intent(in) :: unit
3291 : character(len=*),intent(out) :: msg
3292 :
3293 : !Local variables-------------------------------
3294 : integer :: irec
3295 : !*************************************************************************
3296 :
3297 2974 : ierr = 0
3298 23378 : do irec=1,5 + hdr%npsp
3299 23378 : backspace(unit=unit, err=10, iomsg=msg)
3300 : end do
3301 :
3302 2974 : if (hdr%usepaw == 1) then
3303 0 : do irec=1,2
3304 0 : backspace(unit=unit, err=10, iomsg=msg)
3305 : end do
3306 : end if
3307 :
3308 : return
3309 :
3310 : ! Handle IO-error
3311 2974 : 10 ierr = 1
3312 :
3313 : end function hdr_backspace
3314 : !!***
3315 :
3316 : !!****f* m_hdr/hdr_ncwrite
3317 : !! NAME
3318 : !! hdr_ncwrite
3319 : !!
3320 : !! FUNCTION
3321 : !! This subroutine deals with the output of the hdr_type structured variables in ETSF+NETCDF fornat.
3322 : !! It handles variables according to the ETSF format, whenever possible and uses new variables
3323 : !! when not available in the ETSF format.
3324 : !!
3325 : !! INPUTS
3326 : !! fform=kind of the array in the file
3327 : !! ncid=the unit of the open NetCDF file.
3328 : !! [spinat(3, natom)]= Spinat input variable. In principle it should be stored in the hdr
3329 : !! as it actually affects the real space group but this requires changing hdr_type and fform and dealing
3330 : !! with backward compatibility issues especially in Fortran-IO.
3331 : !! [nc_define]=Optional flag. If True, the basic dimensions required by the ETSF specification
3332 : !! are written. Default: False.
3333 : !!
3334 : !! OUTPUT
3335 : !! Only writing
3336 : !!
3337 : !! SOURCE
3338 :
3339 4155 : integer function hdr_ncwrite(hdr, ncid, fform, spinat, nc_define) result(ncerr)
3340 :
3341 : !Arguments ------------------------------------
3342 : !scalars
3343 : integer,intent(in) :: ncid,fform
3344 : logical,optional,intent(in) :: nc_define
3345 : class(hdr_type),target,intent(in) :: hdr
3346 : real(dp),optional,intent(in) :: spinat(3, hdr%natom)
3347 :
3348 : !Local variables-------------------------------
3349 : !scalars
3350 : logical :: my_define
3351 : character(len=etsfio_charlen) :: basis_set,k_dependent,symmorphic
3352 : !character(len=500) :: msg
3353 : !arrays
3354 16533 : integer,allocatable :: arr2d(:,:)
3355 16533 : real(dp),allocatable :: arr3d(:,:,:)
3356 16533 : type(pawrhoij_type),pointer :: rhoij_ptr(:)
3357 : ! *************************************************************************
3358 :
3359 : call check_fform(fform)
3360 :
3361 : !@hdr_type
3362 16533 : my_define = .False.; if (present(nc_define)) my_define = nc_define
3363 16533 : ncerr = nf90_noerr
3364 :
3365 474869 : k_dependent = "no"; if (any(hdr%nband(1) /= hdr%nband)) k_dependent = "yes"
3366 878734 : symmorphic = "no"; if (all(abs(hdr%tnons) < tol6)) symmorphic = "yes"
3367 :
3368 16533 : if (my_define) then
3369 : !call wrtout(std_out, "hdr_ncwrite: defining variables")
3370 16533 : NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
3371 :
3372 : ! Write ETSF-dims
3373 : ncerr = nctk_def_dims(ncid, [ &
3374 : nctkdim_t("max_number_of_states", hdr%mband), &
3375 : nctkdim_t("number_of_atoms", hdr%natom), &
3376 : nctkdim_t("number_of_atom_species", hdr%ntypat), &
3377 : nctkdim_t("number_of_components", hdr%nspden), &
3378 : nctkdim_t("number_of_kpoints", hdr%nkpt), &
3379 : nctkdim_t("number_of_spinor_components", hdr%nspinor), &
3380 : nctkdim_t("number_of_spins", hdr%nsppol), &
3381 : nctkdim_t("number_of_symmetry_operations", hdr%nsym) &
3382 148797 : ])
3383 : !nctkdim_t("nshiftk_orig", ebands%nshiftk_orig), &
3384 : !nctkdim_t("nshiftk", ebands%nshiftk)], &
3385 16533 : NCF_CHECK(ncerr)
3386 :
3387 : ! Define part of geometry section contained in the header.
3388 : ncerr = nctk_def_arrays(ncid, [ &
3389 : ! Atomic structure and symmetry operations
3390 : nctkarr_t("primitive_vectors", "dp", "number_of_cartesian_directions, number_of_vectors"), &
3391 : nctkarr_t("reduced_symmetry_matrices", "int", &
3392 : "number_of_reduced_dimensions, number_of_reduced_dimensions, number_of_symmetry_operations"), &
3393 : nctkarr_t("reduced_symmetry_translations", "dp", "number_of_reduced_dimensions, number_of_symmetry_operations"), &
3394 : nctkarr_t("atom_species", "int", "number_of_atoms"), &
3395 : nctkarr_t("reduced_atom_positions", "dp", "number_of_reduced_dimensions, number_of_atoms"), &
3396 : nctkarr_t("atomic_numbers", "dp", "number_of_atom_species") &
3397 : !nctkarr_t("atom_species_names", "char", "character_string_length, number_of_atom_species"), &
3398 : !nctkarr_t("chemical_symbols", "char", "symbol_length, number_of_atom_species"), &
3399 : ! Atomic information.
3400 : !nctkarr_t("valence_charges", "dp", "number_of_atom_species"), & ! NB: This variable is not written if alchemical
3401 : !nctkarr_t("pseudopotential_types", "char", "character_string_length, number_of_atom_species") &
3402 115731 : ])
3403 16533 : NCF_CHECK(ncerr)
3404 :
3405 : ! Some variables require the "symmorphic" attribute.
3406 16533 : NCF_CHECK(nf90_put_att(ncid, vid("reduced_symmetry_matrices"), "symmorphic", symmorphic))
3407 16533 : NCF_CHECK(nf90_put_att(ncid, vid("reduced_symmetry_translations"), "symmorphic", symmorphic))
3408 :
3409 : ! At this point we have an ETSF-compliant file. Add additional data for internal use in abinit.
3410 16533 : ncerr = nctk_def_arrays(ncid, nctkarr_t('symafm', "int", "number_of_symmetry_operations"))
3411 16533 : NCF_CHECK(ncerr)
3412 :
3413 : ! TODO spinat should be added to the header
3414 16533 : if (present(spinat)) then
3415 4155 : ncerr = nctk_def_arrays(ncid, nctkarr_t("spinat", "dp", "three, number_of_atoms"))
3416 4155 : NCF_CHECK(ncerr)
3417 : end if
3418 :
3419 : ! Define k-points. Note: monkhorst_pack_folding is replaced by kptrlatt and shiftk
3420 : ncerr = nctk_def_arrays(ncid, [ &
3421 : nctkarr_t("reduced_coordinates_of_kpoints", "dp", "number_of_reduced_dimensions, number_of_kpoints"), &
3422 : nctkarr_t("kpoint_weights", "dp", "number_of_kpoints") &
3423 : !nctkarr_t("monkhorst_pack_folding", "int", "number_of_vectors") &
3424 49599 : ])
3425 16533 : NCF_CHECK(ncerr)
3426 :
3427 : ! Define states section.
3428 : ncerr = nctk_def_arrays(ncid, [ &
3429 : nctkarr_t("number_of_states", "int", "number_of_kpoints, number_of_spins"), &
3430 : nctkarr_t("eigenvalues", "dp", "max_number_of_states, number_of_kpoints, number_of_spins"), &
3431 : nctkarr_t("occupations", "dp", "max_number_of_states, number_of_kpoints, number_of_spins"), &
3432 : nctkarr_t("smearing_scheme", "char", "character_string_length") &
3433 82665 : ])
3434 16533 : NCF_CHECK(ncerr)
3435 :
3436 33066 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "number_of_electrons"])
3437 16533 : NCF_CHECK(ncerr)
3438 49599 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "fermi_energy", "smearing_width"])
3439 16533 : NCF_CHECK(ncerr)
3440 16533 : NCF_CHECK(nctk_set_atomic_units(ncid, "smearing_width"))
3441 :
3442 : ! Some variables require the specifications of units.
3443 16533 : NCF_CHECK(nctk_set_atomic_units(ncid, "eigenvalues"))
3444 16533 : NCF_CHECK(nctk_set_atomic_units(ncid, "fermi_energy"))
3445 16533 : NCF_CHECK(nf90_put_att(ncid, vid("number_of_states"), "k_dependent", k_dependent))
3446 :
3447 : ! Define dimensions.
3448 : ncerr = nctk_def_dims(ncid, [&
3449 : nctkdim_t("npsp", hdr%npsp), nctkdim_t("codvsnlen", 8), nctkdim_t("psptitlen", 132)&
3450 66132 : ])
3451 16533 : NCF_CHECK(ncerr)
3452 :
3453 16533 : if (hdr%usewvl==1) then ! Add the BigDFT private dimensions.
3454 0 : ncerr = nctk_def_dims(ncid, nctkdim_t("number_of_wavelet_resolutions", 2))
3455 0 : NCF_CHECK(ncerr)
3456 : end if
3457 :
3458 : ! Define scalars.
3459 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
3460 165330 : "date", "ixc", "intxc", "occopt", "pertcase", "headform", "fform", "usepaw", "usewvl"])
3461 16533 : NCF_CHECK(ncerr)
3462 :
3463 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
3464 165330 : "ecut_eff", "ecutdg", "ecutsm", "etot", "extfpmd_eshift", "residm", "stmbias", "tphysel", "tsmear"])
3465 16533 : NCF_CHECK(ncerr)
3466 :
3467 : ! Multi-dimensional variables.
3468 : ncerr = nctk_def_arrays(ncid, [&
3469 : nctkarr_t("istwfk", "i", "number_of_kpoints"),&
3470 : nctkarr_t("codvsn", "c", "codvsnlen"),&
3471 : nctkarr_t("pspcod", "i", "npsp"),&
3472 : nctkarr_t("pspdat", "i", "npsp"),&
3473 : nctkarr_t("pspso", "i", "npsp"),&
3474 : nctkarr_t("pspxc", "i", "npsp"),&
3475 : nctkarr_t("qptn", "dp", "number_of_reduced_dimensions"),&
3476 : nctkarr_t("so_psp", "i", "npsp"),&
3477 : !nctkarr_t("symafm", "i", "number_of_symmetry_operations"),&
3478 : nctkarr_t("title", "c", "psptitlen, npsp"),&
3479 : nctkarr_t("zionpsp", "dp", "npsp"),&
3480 : nctkarr_t("znuclpsp", "dp", "npsp"),&
3481 214929 : nctkarr_t("lmn_size", "i", "npsp")])
3482 16533 : NCF_CHECK(ncerr)
3483 :
3484 : ! Add the BigDFT private variables.
3485 16533 : if (hdr%usewvl == 1) then
3486 0 : ncerr = nctk_def_arrays(ncid, nctkarr_t("number_of_wavelets", "i", "number_of_wavelet_resolutions"))
3487 0 : NCF_CHECK(ncerr)
3488 : end if
3489 :
3490 16533 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("basis_set", "char", "character_string_length")))
3491 16533 : if (hdr%usewvl == 0) then
3492 33066 : NCF_CHECK(nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "kinetic_energy_cutoff"]))
3493 16533 : NCF_CHECK(nctk_set_atomic_units(ncid, "kinetic_energy_cutoff"))
3494 16533 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("number_of_coefficients", "int", "number_of_kpoints")))
3495 : end if
3496 :
3497 16533 : NCF_CHECK(nf90_put_att(ncid, vid("number_of_states"), "k_dependent", k_dependent))
3498 :
3499 16533 : if (hdr%usewvl == 0) then
3500 : ! Note that here we always use the coarse FFT mesh even if usepaw == 1
3501 : ncerr = nctk_def_dims(ncid, [&
3502 : nctkdim_t("number_of_grid_points_vector1", hdr%ngfft(1)),&
3503 : nctkdim_t("number_of_grid_points_vector2", hdr%ngfft(2)),&
3504 66132 : nctkdim_t("number_of_grid_points_vector3", hdr%ngfft(3))], defmode=.True.)
3505 16533 : NCF_CHECK(ncerr)
3506 : else
3507 0 : ABI_WARNING("Don't know how to define grid_points for wavelets!")
3508 : end if
3509 :
3510 : !write(std_out,*)"hdr%nshiftk_orig,hdr%nshiftk",hdr%nshiftk_orig,hdr%nshiftk
3511 : ncerr = nctk_def_dims(ncid, [&
3512 : nctkdim_t("nshiftk_orig", hdr%nshiftk_orig),&
3513 : nctkdim_t("nshiftk", hdr%nshiftk), &
3514 : nctkdim_t("bantot", hdr%bantot), &
3515 82665 : nctkdim_t("md5_slen", md5_slen)], defmode=.True.)
3516 16533 : NCF_CHECK(ncerr)
3517 :
3518 66132 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "kptopt", "pawcpxocc", "icoulomb"])
3519 16533 : NCF_CHECK(ncerr)
3520 49599 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "nelect", "charge"])
3521 16533 : NCF_CHECK(ncerr)
3522 :
3523 : ncerr = nctk_def_arrays(ncid, [&
3524 : nctkarr_t("kptrlatt_orig", "i", "number_of_reduced_dimensions, number_of_reduced_dimensions"),&
3525 : nctkarr_t("kptrlatt", "i", "number_of_reduced_dimensions, number_of_reduced_dimensions"),&
3526 : nctkarr_t("shiftk_orig", "dp", "number_of_reduced_dimensions, nshiftk_orig"),&
3527 : nctkarr_t("shiftk", "dp", "number_of_reduced_dimensions, nshiftk"), &
3528 : nctkarr_t("amu", "dp", "number_of_atom_species"), &
3529 115731 : nctkarr_t("md5_pseudos", "ch", "md5_slen, npsp") ])
3530 16533 : NCF_CHECK(ncerr)
3531 :
3532 : !call wrtout(std_out, "hdr_ncwrite completed define mode")
3533 : end if ! my_define
3534 :
3535 : ! Switch to write mode.
3536 16533 : NCF_CHECK(nctk_set_datamode(ncid))
3537 :
3538 : ! write ETSF variables.
3539 16533 : if (hdr%usewvl == 0) then
3540 : ! Plane wave case.
3541 16533 : basis_set = "plane_waves"
3542 16533 : NCF_CHECK(nf90_put_var(ncid, vid("basis_set"), basis_set))
3543 16533 : NCF_CHECK(nf90_put_var(ncid, vid("kinetic_energy_cutoff"), hdr%ecut))
3544 16533 : NCF_CHECK(nf90_put_var(ncid, vid("number_of_coefficients"), hdr%npwarr))
3545 : else
3546 : ! Wavelet case.
3547 0 : basis_set = "daubechies_wavelets"
3548 0 : NCF_CHECK(nf90_put_var(ncid, vid("basis_set"), basis_set))
3549 : ! Required variable than should enter the standard.
3550 0 : NCF_CHECK(nf90_put_var(ncid, vid("number_of_wavelets"), hdr%nwvlarr))
3551 : end if
3552 :
3553 : ! Write electrons
3554 16533 : NCF_CHECK(nf90_put_var(ncid, vid("fermi_energy"), hdr%fermie))
3555 16533 : NCF_CHECK(nf90_put_var(ncid, vid("smearing_width"), hdr%tsmear))
3556 16533 : NCF_CHECK(nf90_put_var(ncid, vid("smearing_scheme"), nctk_string_from_occopt(hdr%occopt)))
3557 :
3558 : ! transfer data from (stupid) 1d hdr%nband and hdr%occ in packed form to 2d - 3d matrix with stride
3559 : ! native support for array and array syntax is one of the reasons why we still use Fortran
3560 : ! and we program like in C but without the power of C!o
3561 :
3562 : ! Also, strange problem with Petrus + Nag5: had to explicitly specify nf90_put_var,
3563 : ! with explicit definition of start, count and stride .
3564 : ! Direct calls to NCF_CHECK, see below, were working for selected tests, but not all tests
3565 66132 : ABI_MALLOC(arr2d, (hdr%nkpt, hdr%nsppol))
3566 49599 : arr2d(:,:) = reshape(hdr%nband, [hdr%nkpt, hdr%nsppol])
3567 49599 : ncerr = nf90_put_var(ncid, vid("number_of_states"), arr2d, start=[1,1], count=[hdr%nkpt,hdr%nsppol], stride=[1,1])
3568 16533 : NCF_CHECK(ncerr)
3569 16533 : ABI_FREE(arr2d)
3570 :
3571 82665 : ABI_MALLOC(arr3d, (hdr%mband, hdr%nkpt, hdr%nsppol))
3572 16533 : call hdr_get_occ3d(hdr, arr3d)
3573 16533 : NCF_CHECK(nf90_put_var(ncid, vid("occupations"), arr3d))
3574 16533 : ABI_FREE(arr3d)
3575 :
3576 : ! Write geometry
3577 16533 : NCF_CHECK(nf90_put_var(ncid, vid("primitive_vectors"), hdr%rprimd))
3578 16533 : NCF_CHECK(nf90_put_var(ncid, vid("reduced_symmetry_matrices"), hdr%symrel))
3579 16533 : NCF_CHECK(nf90_put_var(ncid, vid("atom_species"), hdr%typat))
3580 16533 : NCF_CHECK(nf90_put_var(ncid, vid("reduced_symmetry_translations"), hdr%tnons))
3581 16533 : NCF_CHECK(nf90_put_var(ncid, vid("reduced_atom_positions"), hdr%xred))
3582 16533 : NCF_CHECK(nf90_put_var(ncid, vid("atomic_numbers"), hdr%znucltypat))
3583 :
3584 : ! Write k-points.
3585 16533 : NCF_CHECK(nf90_put_var(ncid, vid("reduced_coordinates_of_kpoints"), hdr%kptns))
3586 16533 : NCF_CHECK(nf90_put_var(ncid, vid("kpoint_weights"), hdr%wtk))
3587 :
3588 : ! Write non-ETSF variables.
3589 16533 : NCF_CHECK(nf90_put_var(ncid, vid("codvsn"), hdr%codvsn))
3590 16533 : NCF_CHECK(nf90_put_var(ncid, vid("title"), hdr%title))
3591 :
3592 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
3593 : & "date", "ixc", "intxc", "occopt", "pertcase", "headform", "fform", "usepaw", "icoulomb"],&
3594 314127 : & [hdr%date, hdr%ixc ,hdr%intxc ,hdr%occopt, hdr%pertcase, HDR_LATEST_HEADFORM, fform, hdr%usepaw, hdr%icoulomb])
3595 16533 : NCF_CHECK(ncerr)
3596 :
3597 : ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
3598 : & "ecut_eff", "ecutdg", "ecutsm", "etot", "extfpmd_eshift", "residm", "stmbias", "tphysel", "tsmear"],&
3599 314127 : & [hdr%ecut_eff, hdr%ecutdg, hdr%ecutsm, hdr%etot, hdr%extfpmd_eshift, hdr%residm, hdr%stmbias, hdr%tphysel, hdr%tsmear])
3600 16533 : NCF_CHECK(ncerr)
3601 :
3602 : ! Write Abinit array variables.
3603 :
3604 : ! FIXME Be careful with zionpsp if alchemical mixing!
3605 16533 : NCF_CHECK(nf90_put_var(ncid, vid("istwfk"), hdr%istwfk))
3606 16533 : NCF_CHECK(nf90_put_var(ncid, vid("pspcod"), hdr%pspcod))
3607 16533 : NCF_CHECK(nf90_put_var(ncid, vid("pspdat"), hdr%pspdat))
3608 16533 : NCF_CHECK(nf90_put_var(ncid, vid("pspso"), hdr%pspso))
3609 16533 : NCF_CHECK(nf90_put_var(ncid, vid("pspxc"), hdr%pspxc))
3610 16533 : NCF_CHECK(nf90_put_var(ncid, vid("qptn"), hdr%qptn))
3611 16533 : NCF_CHECK(nf90_put_var(ncid, vid("so_psp"), hdr%so_psp))
3612 16533 : NCF_CHECK(nf90_put_var(ncid, vid("symafm"), hdr%symafm))
3613 16533 : NCF_CHECK(nf90_put_var(ncid, vid("znuclpsp"), hdr%znuclpsp))
3614 16533 : NCF_CHECK(nf90_put_var(ncid, vid("zionpsp"), hdr%zionpsp))
3615 16533 : NCF_CHECK(nf90_put_var(ncid, vid("lmn_size"), hdr%lmn_size))
3616 16533 : NCF_CHECK(nf90_put_var(ncid, vid("usewvl"), hdr%usewvl))
3617 :
3618 16533 : if (present(spinat)) then
3619 4155 : NCF_CHECK(nf90_put_var(ncid, vid("spinat"), spinat))
3620 : end if
3621 :
3622 : ! Write hdr%pawrhoij.
3623 16533 : if (hdr%usepaw == 1) then
3624 : ! Dirty trick to bypass check on the intent, but the problem is in the intent(inout) of pawrhoij_io
3625 2305 : rhoij_ptr => hdr%pawrhoij
3626 : call pawrhoij_io(rhoij_ptr,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,&
3627 2305 : HDR_LATEST_HEADFORM,"Write",form="netcdf")
3628 : end if
3629 :
3630 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
3631 82665 : "kptopt", "pawcpxocc"], [hdr%kptopt, hdr%pawcpxocc])
3632 16533 : NCF_CHECK(ncerr)
3633 :
3634 : ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
3635 82665 : "nelect", "charge"], [hdr%nelect, hdr%cellcharge])
3636 16533 : NCF_CHECK(ncerr)
3637 :
3638 : ! NB: In etsf_io the number of electrons is declared as integer.
3639 : ! We use abinit nelect to store the value as real(dp).
3640 16533 : NCF_CHECK(nf90_put_var(ncid, vid("number_of_electrons"), nint(hdr%nelect)))
3641 16533 : NCF_CHECK(nf90_put_var(ncid, vid("kptrlatt_orig"), hdr%kptrlatt_orig))
3642 16533 : NCF_CHECK(nf90_put_var(ncid, vid("kptrlatt"), hdr%kptrlatt))
3643 16533 : NCF_CHECK(nf90_put_var(ncid, vid("shiftk_orig"), hdr%shiftk_orig))
3644 16533 : NCF_CHECK(nf90_put_var(ncid, vid("shiftk"), hdr%shiftk))
3645 16533 : NCF_CHECK(nf90_put_var(ncid, vid("md5_pseudos"), hdr%md5_pseudos))
3646 16533 : NCF_CHECK(nf90_put_var(ncid, vid("amu"), hdr%amu))
3647 :
3648 37221 : if (hdr%occopt == 9) then
3649 4 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "hole_fermi_energy", "ne_qFD", "nh_qFD"])
3650 1 : NCF_CHECK(ncerr)
3651 1 : NCF_CHECK(nctk_set_atomic_units(ncid, "hole_fermi_energy"))
3652 1 : NCF_CHECK(nf90_put_var(ncid, vid("hole_fermi_energy"), hdr%fermih))
3653 2 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "ivalence"])
3654 1 : NCF_CHECK(ncerr)
3655 3 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: "ivalence"], [hdr%ivalence])
3656 1 : NCF_CHECK(ncerr)
3657 5 : ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: "ne_qFD", "nh_qFD"],[hdr%ne_qFD, hdr%nh_qFD])
3658 : end if
3659 :
3660 : contains
3661 682009 : integer function vid(vname)
3662 : character(len=*),intent(in) :: vname
3663 682009 : vid = nctk_idname(ncid, vname)
3664 : end function vid
3665 :
3666 : end function hdr_ncwrite
3667 : !!***
3668 :
3669 : !!****f* m_hdr/hdr_set_occ
3670 : !! NAME
3671 : !! hdr_set_occ
3672 : !!
3673 : !! FUNCTION
3674 : !! Set the occuations hdr%occ(:) from a 3d array with stride.
3675 : !!
3676 : !! SOURCE
3677 :
3678 26879 : subroutine hdr_set_occ(hdr, occ3d)
3679 :
3680 : !Arguments ------------------------------------
3681 : class(hdr_type),intent(inout) :: hdr
3682 : real(dp),intent(in) :: occ3d(hdr%mband,hdr%nkpt,hdr%nsppol)
3683 :
3684 : !Local variables-------------------------------
3685 : !scalars
3686 : integer :: ii,band,ikpt,spin
3687 : !*************************************************************************
3688 :
3689 26879 : ii = 0
3690 56167 : do spin=1,hdr%nsppol
3691 1745425 : do ikpt=1,hdr%nkpt
3692 14084952 : do band=1,hdr%nband(ikpt + (spin-1) * hdr%nkpt)
3693 12366406 : ii = ii +1
3694 14055664 : hdr%occ(ii) = occ3d(band,ikpt,spin)
3695 : end do
3696 : end do
3697 : end do
3698 :
3699 26879 : end subroutine hdr_set_occ
3700 : !!***
3701 :
3702 : !!****f* m_hdr/hdr_get_occ3d
3703 : !! NAME
3704 : !! hdr_get_occ3d
3705 : !!
3706 : !! FUNCTION
3707 : !! Return occupations in a 3d array with stride.
3708 : !!
3709 : !! SOURCE
3710 :
3711 17643 : subroutine hdr_get_occ3d(hdr, occ3d)
3712 :
3713 : !Arguments ------------------------------------
3714 : class(hdr_type),intent(in) :: hdr
3715 : real(dp),intent(out) :: occ3d(hdr%mband,hdr%nkpt,hdr%nsppol)
3716 :
3717 : !Local variables-------------------------------
3718 : !scalars
3719 : integer :: ii,band,ikpt,spin
3720 : !*************************************************************************
3721 :
3722 4468029 : ii = 0; occ3d = huge(one)
3723 36986 : do spin=1,hdr%nsppol
3724 603929 : do ikpt=1,hdr%nkpt
3725 4449510 : do band=1,hdr%nband(ikpt + (spin-1) * hdr%nkpt)
3726 3863224 : ii = ii +1
3727 4430167 : occ3d(band,ikpt,spin) = hdr%occ(ii)
3728 : end do
3729 : end do
3730 : end do
3731 :
3732 17643 : end subroutine hdr_get_occ3d
3733 : !!***
3734 :
3735 : !----------------------------------------------------------------------
3736 :
3737 : !!****f* m_hdr/hdr_check
3738 : !! NAME
3739 : !! hdr_check
3740 : !!
3741 : !! FUNCTION
3742 : !! This subroutine compare the header structured variable (hdr)
3743 : !! from input data (mostly dtset and psps) with the one (hdr0) of
3744 : !! an input data file (e.g. wf, density, potential).
3745 : !! Various values are checked for agreement or near agreement in the
3746 : !! case of floating point numbers. The program will exit or produce
3747 : !! warning messages when unexpected values are found.
3748 : !! A record of the comparison of the headers is written to stdout.
3749 : !!
3750 : !! Decisions have been taken about whether a restart is allowed.
3751 : !! In the self-consistent case, a restart will always be allowed, but
3752 : !! one has to distinguish between a direct restart and a restart with
3753 : !! translation of wavefunction.
3754 : !! In the non-self-consistent case, the conditions below
3755 : !! must be fulfilled to allow a restart.
3756 : !!
3757 : !! INPUTS
3758 : !! fform=integer specification of data type (expected)
3759 : !! fform0=integer specification of data type (from disk file)
3760 : !! mode_paral: COLL or PERS, for all wrtout calls
3761 : !! hdr <type(hdr_type)>=the header structured variable from dtset and psps
3762 : !! hdr0<type(hdr_type)>=the header structured variable from the disk file
3763 : !!
3764 : !! OUTPUT
3765 : !! restart=1 if direct restart, =2 if translation is needed, =0 if no
3766 : !! restart is possible.
3767 : !! restartpaw= deals with the additional information in the PAW method
3768 : !! =1 if direct restart, =0 if no restart from spherical data is possible.
3769 : !! also 0 if no restart is possible
3770 : !!
3771 : !! NOTES
3772 : !! In the current version of the user interface restarts are allowed from
3773 : !! wavefunction files for self-consistent runs and from densities for
3774 : !! non-self-consistent runs. The precise conditions under which we will
3775 : !! allow a restart in this release are as follows.
3776 : !!
3777 : !! self-consistent case : direct restarts
3778 : !! ======================================
3779 : !!
3780 : !! A direct restart will be allowed provided the following quantities in
3781 : !! old and new calculations are the same:
3782 : !!
3783 : !! (A) the primitive vectors (tprim)
3784 : !! (B) the plane-wave cutoff (tecut)
3785 : !! (C) nkpt, kpt(3,nkpt), wtk(nkpt) (tkpt)
3786 : !! (D) istwfk(nkpt), the format of wavefunctions (twfk)
3787 : !! (E) nspinor, the scalar or spinor wf characteristics (tspinor)
3788 : !! For PAW calculations:
3789 : !! (F) the use of PAW method (tpaw)
3790 : !! (G) the number of lmn elements for the paw basis (tlmn)
3791 : !! (H) the energy cutoff for the double (fine) grid (tdg)
3792 : !! For WVL calculations:
3793 : !! (I) the number of wavelets differs (twvl)
3794 : !! (J) the space-grid size differs (tgrid)
3795 : !!
3796 : !! non-self-consistent restarts
3797 : !! ============================
3798 : !!
3799 : !! A restart will be allowed provided the following quantities in
3800 : !! old and new calculation are the same
3801 : !!
3802 : !! (A) the primitive vectors (tprim)
3803 : !! (B) the number of atoms of each type (tatty)
3804 : !! (C) xred(3,natom) (txred)
3805 : !! (D) pseudopotentials (not just pseudocharges) (tpseu)
3806 : !! (E) the plane-wave cutoff (tecut)
3807 : !! (F) ngfft(1:3) (tng)
3808 : !! For PAW calculations:
3809 : !! (G) the use of PAW method (tpaw)
3810 : !! (H) the number of lmn elements for the paw basis (tlmn)
3811 : !! (I) the energy cutoff for the double (fine) grid (tdg)
3812 : !!
3813 : !! SOURCE
3814 :
3815 4235 : subroutine hdr_check(fform, fform0, hdr, hdr0, mode_paral, restart, restartpaw)
3816 :
3817 : !Arguments ------------------------------------
3818 : !scalars
3819 : class(hdr_type),intent(in) :: hdr,hdr0
3820 : integer,intent(in) :: fform,fform0
3821 : integer,intent(out) :: restart,restartpaw
3822 : character(len=4),intent(in) :: mode_paral
3823 :
3824 : !Local variables-------------------------------
3825 : character(len=500) :: bndfmt, occfmt, wtkfmt, zatfmt, typfmt
3826 : !scalars
3827 : integer,parameter :: mwarning=3,nkpt_max=5
3828 : integer :: bantot,bantot_eff,ii,ipsp,isppol,istart,istop,isym,itest,iwarning
3829 : integer :: jj,mu,natom,nelm,nkpt,npsp,nsppol,nsym,ntypat,tatty,tband,tdg
3830 : integer :: tecut,tgrid,tkpt,tlmn,tng,tpaw,tprim,tpsch,tpseu,tspinor,tsym,twfk
3831 : integer :: twvl,txred,enough
3832 : real(dp) :: rms
3833 : logical :: tfform2,tfform52
3834 : character(len=500) :: msg
3835 : type(abifile_t) :: abifile,abifile0
3836 : ! *************************************************************************
3837 :
3838 : !@hdr_type
3839 : DBG_ENTER("COLL")
3840 :
3841 : ! We will adopt convention that if things agree between restart
3842 : ! and current calculation then the tflag is 0. Begin by assuming
3843 : ! that there is complete agreement between the files
3844 :
3845 4235 : tatty = 0; tband = 0; tdg = 0 ; tecut = 0; tkpt = 0;
3846 4235 : tlmn = 0; tng = 0; tpaw = 0; tprim = 0; tpsch = 0; tpseu = 0;
3847 4235 : tspinor=0; tsym = 0; twfk = 0 ; txred = 0 ; twvl = 0 ; tgrid = 0
3848 :
3849 : ! Write out a header
3850 : write(msg,'(a1,80a,2a1,10x,a,3a1,10x,a,27x,a,a1,10x,19a,27x,12a,a1)' )&
3851 347270 : ch10,('=',ii=1,80),ch10,ch10,&
3852 4235 : '- hdr_check: checking restart file header for consistency -',&
3853 491260 : (ch10,ii=1,3),'current calculation','restart file',ch10,('-',ii=1,19),('-',ii=1,12),ch10
3854 4235 : call wrtout(std_out,msg,mode_paral)
3855 :
3856 : ! Check validity of fform, and find filetype
3857 4235 : abifile = abifile_from_fform(fform)
3858 4235 : if (abifile%fform == 0) then
3859 0 : ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform)))
3860 : end if
3861 :
3862 : ! Check validity of fform0, and find filetype
3863 4235 : abifile0 = abifile_from_fform(fform0)
3864 4235 : if (abifile0%fform == 0) then
3865 0 : ABI_ERROR(sjoin("Cannot find any abifile object associated to fform:", itoa(fform0)))
3866 : end if
3867 :
3868 : write(msg,'(a,a17,3x,2a,a17)') &
3869 4235 : ' calculation expects a ',ljust(abifile%class, 17),'|',' input file contains a ',ljust(abifile0%class, 17)
3870 4235 : call wrtout(std_out,msg,mode_paral)
3871 :
3872 : write(msg,'(a,a,13x,a,a,a)')&
3873 4235 : '. ABINIT code version ',hdr%codvsn,'|',' ABINIT code version ',hdr0%codvsn
3874 4235 : call wrtout(std_out,msg,mode_paral)
3875 :
3876 : ! Check fform from input, not from header file
3877 4235 : if ( fform /= fform0) then
3878 0 : write(msg,'(a,i0,a,i0,a)')'input fform=',fform,' differs from disk file fform=',fform0,'.'
3879 0 : ABI_ERROR(msg)
3880 : end if
3881 :
3882 : write(msg, '(a,i8,a,i8,a,i4,2x,a,a,i8,a,i8,a,i4)' ) &
3883 4235 : '. date ',hdr %date,' bantot ',hdr %bantot,' natom ',hdr %natom,'|',&
3884 8470 : ' date ',hdr0%date,' bantot ',hdr0%bantot,' natom ',hdr0%natom
3885 4235 : call wrtout(std_out,msg,mode_paral)
3886 :
3887 : write(msg, '(a,i8,a,i3,3(a,i4),2x,a,a,i8,a,i3,3(a,i4))' )&
3888 4235 : ' nkpt',hdr %nkpt,' nsym',hdr %nsym,' ngfft',hdr %ngfft(1),',',hdr %ngfft(2),',',hdr %ngfft(3),'|',&
3889 8470 : ' nkpt',hdr0%nkpt,' nsym',hdr0%nsym,' ngfft',hdr0%ngfft(1),',',hdr0%ngfft(2),',',hdr0%ngfft(3)
3890 4235 : call wrtout(std_out,msg,mode_paral)
3891 :
3892 4235 : if (hdr%usewvl == 0) then
3893 : ! Note that the header actually contains ecut_eff=ecut*dilatmx**2
3894 : write(msg,'(a,i3,a,f12.7,12x,a,a,i3,a,f12.7)')&
3895 4235 : ' ntypat',hdr %ntypat,' ecut_eff',hdr %ecut_eff,'|',&
3896 8470 : ' ntypat',hdr0%ntypat,' ecut_eff',hdr0%ecut_eff
3897 4235 : call wrtout(std_out,msg,mode_paral)
3898 : else
3899 : write(msg,'(a,i3,a,f12.7,12x,a,a,i3,a,f12.7)')&
3900 0 : ' ntypat',hdr %ntypat,' hgrid ', 2. * hdr %rprimd(1,1) / (hdr %ngfft(1) - 31),'|',&
3901 0 : ' ntypat',hdr0%ntypat,' hgrid ', 2. * hdr0%rprimd(1,1) / (hdr0%ngfft(1) - 31)
3902 0 : call wrtout(std_out,msg,mode_paral)
3903 : ! Check hgrid and rprimd values.
3904 : if (hdr0%rprimd(1,2) /= zero .or. hdr0%rprimd(1,3) /= zero .or. &
3905 : hdr0%rprimd(2,1) /= zero .or. hdr0%rprimd(2,3) /= zero .or. &
3906 0 : hdr0%rprimd(3,1) /= zero .or. hdr0%rprimd(3,2) /= zero) then
3907 0 : ABI_ERROR('disk file rprimd is not parallelepipedic.')
3908 : end if
3909 0 : if (abs(hdr0%rprimd(1,1) / hdr0%ngfft(1) - hdr %rprimd(1,1) / hdr %ngfft(1)) > tol8) then
3910 : write(msg,'(a,F7.4,a,F7.4)')&
3911 0 : 'input wvl_hgrid=', 2. * hdr%rprimd(1,1) / hdr%ngfft(1), &
3912 0 : 'not equal disk file wvl_hgrid=', 2. * hdr0%rprimd(1,1) / hdr0%ngfft(1)
3913 0 : ABI_COMMENT(msg)
3914 0 : tgrid = 1
3915 : end if
3916 : end if
3917 :
3918 4235 : write(msg, '(a,i3,33x,a,a,i3)' )' usepaw',hdr %usepaw,'|',' usepaw',hdr0%usepaw
3919 4235 : call wrtout(std_out,msg,mode_paral)
3920 :
3921 4235 : write(msg, '(a,i3,33x,a,a,i3)' )' usewvl',hdr %usewvl,'|',' usewvl',hdr0%usewvl
3922 4235 : call wrtout(std_out,msg,mode_paral)
3923 :
3924 : write(msg,'(a,35x,a,a,3(a1,2x,3f12.7,6x,a,2x,3f12.7))')&
3925 4235 : ' rprimd:','|',' rprimd:',ch10,&
3926 4235 : hdr%rprimd(:,1),'|',hdr0%rprimd(:,1),ch10,&
3927 4235 : hdr%rprimd(:,2),'|',hdr0%rprimd(:,2),ch10,&
3928 8470 : hdr%rprimd(:,3),'|',hdr0%rprimd(:,3)
3929 4235 : call wrtout(std_out,msg,mode_paral)
3930 :
3931 4235 : if (hdr%bantot/=hdr0%bantot) tband=1
3932 :
3933 4235 : if (hdr%intxc/=hdr0%intxc) then
3934 2 : write(msg,'(a,i0,a,i0)')'input intxc=',hdr%intxc,' not equal disk file intxc=',hdr0%intxc
3935 2 : ABI_WARNING(msg)
3936 : end if
3937 :
3938 4235 : if (hdr%ixc/=hdr0%ixc) then
3939 134 : write(msg,'(a,i0,a,i0)')'input ixc=',hdr%ixc,' not equal disk file ixc=',hdr0%ixc
3940 134 : ABI_WARNING(msg)
3941 : end if
3942 :
3943 4235 : if (hdr%natom/=hdr0%natom) then
3944 0 : write(msg,'(a,i0,a,i0)')'input natom=',hdr%natom,' not equal disk file natom=',hdr0%natom
3945 0 : ABI_WARNING(msg)
3946 0 : tatty=1
3947 : end if
3948 :
3949 16818 : if ( ANY(hdr%ngfft/=hdr0%ngfft) ) then
3950 : ! For sensible rho(r) or V(r) data, fft grid must be identical
3951 : ! Note, however, that we allow for different FFT meshes and we interpolate the density in the
3952 : ! caller when we are restarting a SCF calculation.
3953 41 : if (abifile%class == "density" .or. abifile%class == "potential") then
3954 : write(msg, '(10a)' )&
3955 0 : 'FFT grids must be the same to restart from a ',trim(abifile%class),' file.',ch10,&
3956 0 : "ngfft from file: ", trim(ltoa(hdr0%ngfft(1:3))), ", from input: ", trim(ltoa(hdr%ngfft(1:3))), ch10, &
3957 0 : 'Action: change the FFT grid in the input via ngfft or change the restart file.'
3958 0 : ABI_ERROR(msg)
3959 : end if
3960 : tng=1
3961 : end if
3962 :
3963 4235 : if (hdr%nkpt/=hdr0%nkpt) then
3964 650 : if (abifile%class == "wf_planewave") then
3965 197 : write(msg,'(a,i0,a,i0)' )'input nkpt=',hdr%nkpt,' not equal disk file nkpt=',hdr0%nkpt
3966 197 : ABI_COMMENT(msg)
3967 : end if
3968 : tkpt=1; twfk=1
3969 : end if
3970 :
3971 4235 : if (hdr%nspinor/=hdr0%nspinor) then
3972 16 : if (abifile%class == "wf_planewave") then
3973 13 : write(msg,'(a,i0,a,i0)')'input nspinor=',hdr%nspinor,' not equal disk file nspinor=',hdr0%nspinor
3974 13 : ABI_WARNING(msg)
3975 : end if
3976 : tspinor=1
3977 : end if
3978 :
3979 : ! No check is present for nspden
3980 4235 : if (hdr%nsppol/=hdr0%nsppol) then
3981 6 : write(msg,'(a,i0,a,i0)')'input nsppol=',hdr%nsppol,' not equal disk file nsppol=',hdr0%nsppol
3982 6 : ABI_WARNING(msg)
3983 : end if
3984 :
3985 4235 : if (hdr%nsym/=hdr0%nsym) then
3986 110 : write(msg, '(a,i0,a,i0)' )'input nsym=',hdr%nsym,' not equal disk file nsym=',hdr0%nsym
3987 110 : ABI_WARNING(msg)
3988 110 : tsym=1
3989 : end if
3990 :
3991 4235 : if (hdr%ntypat/=hdr0%ntypat) then
3992 0 : write(msg,'(a,i0,a,i0)')'input ntypat=',hdr%ntypat,' not equal disk file ntypat=',hdr0%ntypat
3993 0 : call wrtout(std_out,msg,mode_paral)
3994 0 : ABI_WARNING(msg)
3995 0 : tatty=1
3996 : end if
3997 :
3998 4235 : if (hdr%usepaw/=hdr0%usepaw) then
3999 0 : write(msg,'(a,i0,a,i0)')'input usepaw=',hdr%usepaw,' not equal disk file usepaw=',hdr0%usepaw
4000 0 : ABI_WARNING(msg)
4001 0 : tpaw=1
4002 : end if
4003 :
4004 4235 : if (hdr%usewvl/=hdr0%usewvl) then
4005 : write(msg, '(a,i6,a,i6,a,a)' )&
4006 0 : 'input usewvl=',hdr%usewvl,' not equal disk file usewvl=',hdr0%usewvl, ch10, &
4007 0 : 'Action: change usewvl input variable or your restart file.'
4008 0 : ABI_ERROR(msg)
4009 : end if
4010 :
4011 : ! Also examine agreement of floating point data
4012 4235 : if (hdr%usewvl == 0 .and. abs(hdr%ecut_eff-hdr0%ecut_eff)>tol8) then
4013 28 : write(msg,'(a,f12.6,a,f12.6,a)')'input ecut_eff=',hdr%ecut_eff,' /= disk file ecut_eff=',hdr0%ecut_eff,'.'
4014 28 : ABI_WARNING(msg)
4015 28 : tecut=1
4016 : end if
4017 :
4018 16940 : do ii=1,3
4019 55055 : do jj=1,3
4020 50820 : if (abs(hdr%rprimd(ii,jj)-hdr0%rprimd(ii,jj))>tol6) then
4021 : write(msg, '(a,i1,a,i1,a,1p,e17.9,a,i1,a,i1,a,e17.9)' )&
4022 369 : 'input rprimd(',ii,',',jj,')=',hdr%rprimd(ii,jj),' /= disk file rprimd(',ii,',',jj,')=',hdr0%rprimd(ii,jj)
4023 369 : ABI_WARNING(msg)
4024 369 : tprim=1
4025 : end if
4026 : end do
4027 : end do
4028 :
4029 : ! Below this point many comparisons only make sense if
4030 : ! certain things agree, e.g. nkpt, natom. Also have to
4031 : ! accomodate different amounts of data in general.
4032 :
4033 4235 : if (hdr%usepaw==1 .and. hdr0%usepaw==1) then
4034 :
4035 : ! Compare ecutdg (PAW)
4036 972 : write(msg, '(a,f12.6,19x,a,a,f12.6)' )' PAW: ecutdg',hdr %ecutdg,'|',' PAW: ecutdg',hdr0%ecutdg
4037 972 : call wrtout(std_out,msg,mode_paral)
4038 972 : if (hdr%ecutdg/=hdr0%ecutdg) then
4039 14 : write(msg, '(a,f12.6,a,f12.6)' )'input ecutdg=',hdr%ecutdg,'not equal disk file ecutdg=',hdr0%ecutdg
4040 14 : ABI_WARNING(msg)
4041 14 : tdg=1
4042 : end if
4043 : end if
4044 :
4045 : ! Compare nband(nkpt*nsppol) (cannot compare if nkpt and nsppol not same)
4046 4235 : if (hdr%nkpt==hdr0%nkpt .and. hdr%nsppol==hdr0%nsppol) then
4047 3579 : nkpt=hdr%nkpt ; nsppol=hdr%nsppol
4048 3579 : write(msg,'(a,36x,a,a)') ' nband:','|',' nband:'
4049 3579 : call wrtout(std_out,msg,mode_paral)
4050 20559 : do istart = 1,nsppol*nkpt,9
4051 16980 : istop = min(istart + 8,nsppol*nkpt)
4052 16980 : mu = istop - istart + 1
4053 : ! generate a format specifier
4054 16980 : bndfmt = strcat('(2x,',itoa(mu),'i4,t41,a,2x,',itoa(mu),'i4)')
4055 20559 : if (istart<=100) then
4056 11817 : write(msg,fmt=bndfmt) hdr%nband(istart:istop),' |',hdr0%nband(istart:istop)
4057 11817 : call wrtout(std_out,msg,mode_paral)
4058 11817 : if (istop>100) call wrtout(std_out, '=> stop printing nband after 100 values', mode_paral)
4059 : end if
4060 : end do
4061 :
4062 : enough = 0
4063 8214 : do isppol=1,nsppol
4064 144663 : do ii=1,nkpt
4065 141084 : if (hdr%nband(ii)/=hdr0%nband(ii)) then
4066 361 : tband=1
4067 361 : enough = enough + 1
4068 361 : if (abifile%class == "wf_planewave") then
4069 119 : if (enough > 5) then
4070 8 : write(std_out, "(a)")"Stop writing warnings after 5 values"
4071 8 : exit
4072 : else
4073 : write(msg,'(a,i0,a,i0,a,i0)' )&
4074 111 : 'kpt num ',ii,' input nband= ',hdr%nband(ii),' not equal disk file nband=',hdr0%nband(ii)
4075 111 : ABI_WARNING(msg)
4076 : end if
4077 : end if
4078 : end if
4079 : end do
4080 : end do
4081 : end if
4082 :
4083 : ! Compare the number of wavelets in each resolution.
4084 : if (hdr%usewvl == 1) then
4085 : if (size(hdr%nwvlarr) /= size(hdr0%nwvlarr) .or. size(hdr%nwvlarr) /= 2) then
4086 : write(msg, '(a,i0,a,i0,a,a)' )&
4087 : 'input nwvlres= ',size(hdr%nwvlarr),' not equal disk file nwvlres= ',size(hdr0%nwvlarr),' or 2',&
4088 : ' ABINIT is not implemented for wavelet resolutions different from 2.'
4089 : ABI_ERROR(msg)
4090 : end if
4091 : end if
4092 :
4093 : ! Compare symmetry arrays (integers) symafm(nsym)
4094 : ! only for same number of symmetries nsym
4095 4235 : itest=0
4096 4235 : if (hdr%nsym==hdr0%nsym) then
4097 4125 : nsym=hdr%nsym
4098 4125 : write(msg,'(a,35x,a,a)') ' symafm:','|',' symafm:'
4099 4125 : call wrtout(std_out,msg,mode_paral)
4100 14131 : do istart = 1,nsym,12
4101 10006 : istop=min(istart+11,nsym)
4102 10006 : nelm = istop - istart + 1
4103 10006 : typfmt = strcat('(2x,',itoa(nelm),'i3,t41,a,2x,',itoa(nelm),'i3)')
4104 10006 : write(msg,fmt=typfmt) hdr%symafm(istart:istop),' |',hdr0%symafm(istart:istop)
4105 14131 : call wrtout(std_out,msg,mode_paral)
4106 : end do
4107 : end if
4108 :
4109 4235 : if (itest/=0) then
4110 0 : write(msg,'(a,i0,a)' )'For symmetry number',itest,' input symafm not equal disk file symafm'
4111 0 : ABI_WARNING(msg)
4112 0 : tsym=1
4113 : end if
4114 :
4115 : ! Compare symmetry arrays (integers) symrel(3,3,nsym)
4116 : ! only for same number of symmetries nsym
4117 4235 : itest=0
4118 4235 : if (hdr%nsym==hdr0%nsym) then
4119 4125 : nsym=hdr%nsym
4120 4125 : write(msg,'(a,35x,a,a)') ' symrel:','|',' symrel:'
4121 4125 : call wrtout(std_out,msg,mode_paral)
4122 106574 : do isym=1,nsym
4123 102450 : write(msg,'(2x,9i3,15x,a,2x,9i3)')hdr%symrel(:,:,isym),'|',hdr0%symrel(:,:,isym)
4124 102450 : call wrtout(std_out,msg,mode_paral)
4125 1335974 : if(sum(abs(hdr%symrel(:,:,isym)-hdr0%symrel(:,:,isym)))/=0)then
4126 1 : itest=isym
4127 1 : exit
4128 : end if
4129 : end do
4130 : end if
4131 :
4132 4235 : if (itest/=0) then
4133 1 : write(msg,'(a,i0,a)')'For symmetry number',itest,' input symrel not equal disk file symrel'
4134 1 : ABI_WARNING(msg)
4135 1 : tsym=1
4136 : end if
4137 :
4138 : ! Compare typat(natom)
4139 4235 : if (hdr%natom==hdr0%natom) then
4140 4235 : natom=hdr%natom
4141 4235 : write(msg,'(a,36x,a,a)') ' typat:','|',' typat:'
4142 4235 : call wrtout(std_out,msg,mode_paral)
4143 8471 : do istart = 1,natom,12
4144 4236 : istop=min(istart+11,natom)
4145 4236 : nelm = istop - istart + 1
4146 4236 : typfmt = strcat('(2x,',itoa(nelm),'i3,t41,a,2x,',itoa(nelm),'i3)')
4147 4236 : write(msg,fmt=typfmt) hdr%typat(istart:istop),' |',hdr0%typat(istart:istop)
4148 8471 : call wrtout(std_out,msg,mode_paral)
4149 : end do
4150 13080 : do ii=1,natom
4151 13080 : if (hdr%typat(ii)/=hdr0%typat(ii)) then
4152 : write(msg, '(a,i0,a,i0,a,i0)' )&
4153 2 : 'For atom number ',ii,' input typat=',hdr%typat(ii),' not equal disk file typat=',hdr0%typat(ii)
4154 2 : ABI_WARNING(msg)
4155 2 : tatty=1
4156 : end if
4157 : end do
4158 : end if
4159 :
4160 : ! Compare so_psp(npsp)
4161 4235 : if (hdr%npsp==hdr0%npsp) then
4162 4235 : npsp=hdr%npsp
4163 4235 : write(msg,'(a,33x,a,a)') ' so_psp :','|',' so_psp :'
4164 4235 : call wrtout(std_out,msg,mode_paral)
4165 8470 : do istart = 1,npsp ,12
4166 4235 : istop=min(istart+11,npsp )
4167 4235 : nelm = istop - istart + 1
4168 4235 : typfmt = strcat('(2x,',itoa(nelm),'i3,t41,a,2x,',itoa(nelm),'i3)')
4169 4235 : write(msg,fmt=typfmt) hdr%so_psp (istart:istop),' |',hdr0%so_psp (istart:istop)
4170 8470 : call wrtout(std_out,msg,mode_paral)
4171 : end do
4172 10285 : do ii=1,npsp
4173 10285 : if (hdr%so_psp (ii)/=hdr0%so_psp (ii)) then
4174 : write(msg,'(a,i0,a,i0,a,i0)')&
4175 5 : 'For pseudopotential number ',ii,' input so_psp =',hdr%so_psp(ii),' not equal disk file so_psp=',hdr0%so_psp(ii)
4176 5 : ABI_WARNING(msg)
4177 : end if
4178 : end do
4179 : end if
4180 :
4181 : ! Compare istwfk(nkpt)
4182 4235 : if (hdr%nkpt==hdr0%nkpt) then
4183 3585 : nkpt=hdr%nkpt
4184 3585 : write(msg,'(a,35x,a,a)') ' istwfk:','|',' istwfk:'
4185 3585 : call wrtout(std_out,msg,mode_paral)
4186 16265 : do istart = 1,nkpt,12
4187 12680 : istop=min(istart+11,nkpt)
4188 12680 : nelm = istop - istart + 1
4189 12680 : typfmt = strcat('(2x,',itoa(nelm),'i3,t41,a,2x,',itoa(nelm),'i3)')
4190 16265 : if (istart<=100) then
4191 8915 : write(msg,fmt=typfmt) hdr%istwfk(istart:istop),' |',hdr0%istwfk(istart:istop)
4192 8915 : call wrtout(std_out,msg,mode_paral)
4193 8915 : if (istop>100) then
4194 551 : call wrtout(std_out, '=> stop printing istwfk after 100 values' ,mode_paral)
4195 : end if
4196 : end if
4197 : end do
4198 131904 : do ii=1,nkpt
4199 131904 : if (hdr%istwfk(ii)/=hdr0%istwfk(ii)) then
4200 : write(msg, '(a,i0,a,i0,a,i0)' )&
4201 440 : 'For k point number ',ii,' input istwfk=',hdr%istwfk(ii),' not equal disk file istwfk=',hdr0%istwfk(ii)
4202 440 : ABI_COMMENT(msg)
4203 440 : twfk=1
4204 : end if
4205 : end do
4206 : end if
4207 :
4208 : !NEW_HDR
4209 51208 : if (any(hdr%kptrlatt /= hdr0%kptrlatt)) then
4210 323 : write(msg,"(2(a,9(i0,1x)))")"input kptrlatt = ",hdr%kptrlatt," /= disk file kptrlatt = ",hdr0%kptrlatt
4211 323 : ABI_COMMENT(msg)
4212 : end if
4213 4235 : if (hdr%kptopt /= hdr0%kptopt) then
4214 488 : ABI_COMMENT(sjoin("input kptopt = ", itoa(hdr%kptopt)," /= disk file kptopt = ", itoa(hdr0%kptopt)))
4215 : end if
4216 4235 : if (hdr%pawcpxocc /= hdr0%pawcpxocc) then
4217 3 : ABI_WARNING(sjoin("input pawcpxocc = ", itoa(hdr%pawcpxocc)," /= disk file pawcpxocc = ", itoa(hdr0%pawcpxocc)))
4218 : end if
4219 4235 : if (hdr%icoulomb /= hdr0%icoulomb) then
4220 0 : ABI_WARNING(sjoin("input icoulomb = ", itoa(hdr%icoulomb)," /= disk file icoulomb = ", itoa(hdr0%icoulomb)))
4221 : end if
4222 :
4223 4235 : if (abs(hdr%nelect - hdr0%nelect) > tol6) then
4224 6 : ABI_WARNING(sjoin("input nelect = ", ftoa(hdr%nelect)," /= disk file nelect = ",ftoa(hdr0%nelect)))
4225 : end if
4226 :
4227 4235 : if (abs(hdr%ne_qFD - hdr0%ne_qFD) > tol6) then
4228 0 : ABI_WARNING(sjoin("input ne_qFD = ", ftoa(hdr%ne_qFD)," /= disk file nelect = ",ftoa(hdr0%ne_qFD)))
4229 : end if
4230 4235 : if (abs(hdr%nh_qFD - hdr0%nh_qFD) > tol6) then
4231 0 : ABI_WARNING(sjoin("input nh_qFD = ", ftoa(hdr%nh_qFD)," /= disk file nelect = ",ftoa(hdr0%nh_qFD)))
4232 : end if
4233 4235 : if (hdr%ivalence/=hdr0%ivalence) then
4234 4151 : write(msg,'(a,i0,a,i0)')'input ival=',hdr%ivalence,' not equal disk file ival=',hdr0%ivalence
4235 4151 : ABI_WARNING(msg)
4236 : end if
4237 :
4238 4235 : if (abs(hdr%cellcharge - hdr0%cellcharge) > tol6) then
4239 0 : ABI_WARNING(sjoin("input cellcharge = ", ftoa(hdr%cellcharge)," /= disk file cellcharge = ", ftoa(hdr0%cellcharge)))
4240 : end if
4241 :
4242 4235 : if (hdr%ntypat==hdr0%ntypat) then
4243 10281 : if (any(abs(hdr%amu - hdr0%amu) > tol6)) then
4244 2 : ABI_WARNING(sjoin("input amu = ",ltoa(hdr%amu)," /= disk file amu = ",ltoa(hdr0%amu)))
4245 : end if
4246 : end if
4247 : !end NEW_HDR
4248 :
4249 : ! Compare kpt(3,nkpt)
4250 4235 : if (hdr%nkpt==hdr0%nkpt) then
4251 3585 : nkpt=hdr%nkpt
4252 3585 : write(msg,'(a,38x,a,a)') ' kpt:','|',' kpt:'
4253 3585 : call wrtout(std_out,msg,mode_paral)
4254 18038 : do ii = 1,min(nkpt,nkpt_max)
4255 14453 : write(msg,'(2x,3f12.7,2x,a,2x,3f12.7)')hdr%kptns(:,ii),' |',hdr0%kptns(:,ii)
4256 14453 : call wrtout(std_out,msg,mode_paral)
4257 3585 : if(ii>nkpt_max)then
4258 : call wrtout(std_out,'The number of printed k points is sufficient... stop writing them.',mode_paral)
4259 : exit
4260 : end if
4261 : end do
4262 3585 : iwarning=0
4263 130129 : do ii=1,nkpt
4264 126636 : itest=0
4265 506544 : do mu=1,3
4266 506544 : if(abs( hdr%kptns(mu,ii)-hdr0%kptns(mu,ii) )>tol6)itest=1
4267 : end do
4268 130129 : if (itest==1) then
4269 : write(msg, '(a,i5,a,3es17.7,a,a,3es17.7)' )&
4270 299 : 'kpt num',ii,', input kpt=',hdr%kptns(:,ii),ch10,&
4271 598 : 'not equal disk file kpt=',hdr0%kptns(:,ii)
4272 299 : ABI_WARNING(msg)
4273 299 : tkpt=1 ; iwarning=iwarning+1
4274 299 : if(iwarning>=mwarning)then
4275 92 : call wrtout(std_out,'The number of comment messages is sufficient ... stop writing them.',mode_paral)
4276 92 : exit
4277 : end if
4278 : end if
4279 : end do
4280 : end if
4281 :
4282 : ! Compare wtk(nkpt)
4283 4235 : if (hdr%nkpt==hdr0%nkpt) then
4284 3585 : nkpt=hdr%nkpt
4285 :
4286 3585 : write(msg,'(a,38x,a,a)') ' wtk:','|',' wtk:'
4287 3585 : call wrtout(std_out,msg,mode_paral)
4288 3585 : istop = min(nkpt,nkpt_max)
4289 7170 : do ii = 1, istop, 5
4290 3585 : mu = min(5, istop - ii + 1)
4291 3585 : wtkfmt = strcat('(2x,',itoa(mu),'f7.3,t41,a,2x,',itoa(mu),'f7.3)')
4292 3585 : write(msg, wtkfmt)hdr%wtk(ii:min(istop, ii + 5 - 1)),' |',hdr0%wtk(ii:min(istop, ii + 5 - 1))
4293 7170 : call wrtout(std_out,msg,mode_paral)
4294 : end do
4295 3585 : iwarning=0
4296 131882 : do ii=1,nkpt
4297 128308 : itest=0
4298 131882 : if (abs( hdr%wtk(ii)-hdr0%wtk(ii) )>tol6) then
4299 : write(msg,'(a,i5,a,es17.7,a,a,es17.7)')&
4300 59 : 'kpt num',ii,', input weight=',hdr%wtk(ii),ch10,&
4301 118 : 'not equal to disk file weight=',hdr0%wtk(ii)
4302 59 : ABI_WARNING(msg)
4303 :
4304 59 : tkpt=1 ; iwarning=iwarning+1
4305 59 : if(iwarning>=mwarning)then
4306 11 : call wrtout(std_out,'The number of warning messages is sufficient ... stop writing them.',mode_paral)
4307 11 : exit
4308 : end if
4309 : end if
4310 : end do
4311 : end if
4312 :
4313 : ! Compare occ(bantot)
4314 4235 : if (hdr%nkpt==hdr0%nkpt.and. hdr%bantot==hdr0%bantot) then
4315 3472 : nkpt=hdr%nkpt
4316 3472 : bantot=hdr%bantot
4317 :
4318 3472 : write(msg,'(a,38x,a,a)') ' occ:','|',' occ:'
4319 3472 : call wrtout(std_out,msg,mode_paral)
4320 3472 : bantot_eff=min(bantot,9*nkpt_max)
4321 16974 : do istart = 1,bantot_eff,9
4322 13502 : istop = min(istart+8,bantot_eff)
4323 13502 : mu = istop - istart + 1
4324 13502 : occfmt = strcat('(2x,',itoa(mu),'f4.1,t41,a,2x,',itoa(mu),'f4.1)')
4325 13502 : write(msg,fmt=occfmt)hdr%occ(istart:istop),' |', hdr0%occ(istart:istop)
4326 13502 : call wrtout(std_out,msg,mode_paral)
4327 16974 : if(istart>9*nkpt_max)then
4328 0 : call wrtout(std_out,'The number of printed occupation numbers is sufficient ... stop writing them.',mode_paral)
4329 0 : exit
4330 : end if
4331 : end do
4332 3472 : iwarning=0
4333 697904 : do ii=1,bantot
4334 697904 : if (abs( hdr%occ(ii)-hdr0%occ(ii) )>tol6) then
4335 1037 : write(msg,'(a,i0,a,1p,e15.7,a,e15.7)')'band,k: ',ii,', input occ=',hdr%occ(ii),' disk occ=',hdr0%occ(ii)
4336 1037 : ABI_COMMENT(msg)
4337 1037 : tband=1 ; iwarning=iwarning+1
4338 1037 : if(iwarning>=mwarning)then
4339 339 : call wrtout(std_out,'The number of warning msgs is sufficient ... stop writing them.',mode_paral)
4340 339 : exit
4341 : end if
4342 : end if
4343 : end do
4344 : end if
4345 :
4346 : ! Compare tnons(3,nsym)
4347 4235 : if (hdr%nsym==hdr0%nsym) then
4348 4125 : nsym=hdr%nsym
4349 4125 : itest=0
4350 4125 : write(msg,'(a,36x,a,a)') ' tnons:','|',' tnons:'
4351 4125 : call wrtout(std_out,msg,mode_paral)
4352 106582 : do isym=1,nsym
4353 102457 : write(msg,'(2x,3f12.7,2x,a,2x,3f12.7)') hdr%tnons(:,isym),' |',hdr0%tnons(:,isym)
4354 106582 : call wrtout(std_out,msg,mode_paral)
4355 : end do
4356 :
4357 96908 : do isym=1,nsym
4358 375899 : if( sum(abs( hdr%tnons(:,isym)-hdr0%tnons(:,isym) )) > tol6) then
4359 214 : itest=isym
4360 214 : exit
4361 : end if
4362 : end do
4363 4125 : if (itest/=0) then
4364 214 : write(msg, '(a,i0,a)' )'For symmetry number ',itest,' input tnons not equal disk file tnons'
4365 214 : ABI_WARNING(msg)
4366 : end if
4367 : end if
4368 :
4369 : ! Compare znucltypat(ntypat)
4370 4235 : if (hdr%ntypat==hdr0%ntypat) then
4371 4235 : ntypat=hdr%ntypat
4372 :
4373 4235 : write(msg,'(a,35x,a,a)') ' znucl:','|',' znucl:'
4374 4235 : call wrtout(std_out,msg,mode_paral)
4375 8470 : do istart = 1,ntypat,6
4376 4235 : istop = min(istart+5,ntypat)
4377 4235 : mu = istop-istart+1
4378 4235 : zatfmt = strcat('(2x,',itoa(mu),'f6.2,t41,a,6x,',itoa(mu),'f6.2)')
4379 4235 : write(msg,fmt=zatfmt) hdr%znucltypat(istart:istop),' |',hdr0%znucltypat(istart:istop)
4380 8470 : call wrtout(std_out,msg,mode_paral)
4381 : end do
4382 :
4383 10283 : do ii=1,ntypat
4384 10283 : if (abs(hdr%znucltypat(ii)-hdr0%znucltypat(ii))>tol6) then
4385 : write(msg, '(a,i5,a,f12.6,a,f12.6)' )&
4386 2 : ' For atom number ',ii,' input znucl=',hdr%znucltypat(ii),' not equal disk file znucl=',hdr0%znucltypat(ii)
4387 2 : ABI_WARNING(msg)
4388 : end if
4389 : end do
4390 : end if
4391 :
4392 : ! Should perform some checks related to pertcase and qptn,
4393 : ! that have been introduced in the header in v4.1
4394 : ! Warning: a GS file might be read, while the hdr corresponds
4395 : ! to a RF file (to initialize k+q), and vice-versa (in nonlinear).
4396 :
4397 : ! Now check agreement of psp headers too
4398 4235 : if (hdr%npsp==hdr0%npsp) then
4399 4235 : npsp=hdr%npsp
4400 4235 : itest=0
4401 :
4402 10285 : do ipsp=1,npsp
4403 : write(msg,'(a,i3,a,13x,a,a,i3,a)')&
4404 6050 : ' pseudopotential atom type',ipsp,':','|',' pseudopotential atom type',ipsp,':'
4405 6050 : call wrtout(std_out,msg,mode_paral)
4406 :
4407 6050 : if (hdr%usepaw==1 .and. hdr0%usepaw==1) then
4408 : write(msg,'(a,i3,a,i7,a,i3,5x,a,a,i3,a,i7,a,i3)')&
4409 1769 : ' pspso ',hdr %pspso(ipsp),' pspxc ',hdr %pspxc(ipsp),&
4410 1769 : ' lmn_size ',hdr%lmn_size(ipsp),'|',&
4411 1769 : ' pspso ',hdr0%pspso(ipsp),' pspxc ',hdr0%pspxc(ipsp),&
4412 3538 : ' lmn_size ',hdr0%lmn_size(ipsp)
4413 1769 : call wrtout(std_out,msg,mode_paral)
4414 1769 : if (hdr%lmn_size(ipsp)/=hdr0%lmn_size(ipsp)) then
4415 : write(msg, '(a,i3,a,i3,a,i3)' )&
4416 0 : 'For atom type ',ipsp,' input lmn_size=',hdr%lmn_size(ipsp),&
4417 0 : 'not equal disk file lmn_size=',hdr0%lmn_size(ipsp)
4418 0 : ABI_WARNING(msg)
4419 0 : tlmn=1
4420 : end if
4421 : else
4422 : write(msg,'(a,i3,a,i3,23x,a,a,i3,a,i3)')&
4423 4281 : ' pspso ',hdr %pspso(ipsp),' pspxc ',hdr %pspxc(ipsp),'|',&
4424 8562 : ' pspso ',hdr0%pspso(ipsp),' pspxc ',hdr0%pspxc(ipsp)
4425 4281 : call wrtout(std_out,msg,mode_paral)
4426 : end if
4427 : write(msg,'(a,i8,a,i4,a,f5.1,4x,a,a,i8,a,i4,a,f5.1)')&
4428 6050 : ' pspdat ',hdr %pspdat(ipsp),' pspcod ',hdr %pspcod(ipsp),&
4429 6050 : ' zion ',hdr %zionpsp(ipsp),'|',&
4430 6050 : ' pspdat ',hdr0%pspdat(ipsp),' pspcod ',hdr0%pspcod(ipsp),&
4431 12100 : ' zion ',hdr0%zionpsp(ipsp)
4432 6050 : call wrtout(std_out,msg,mode_paral)
4433 :
4434 : ! Check on md5 values.
4435 6050 : if (hdr%md5_pseudos(ipsp) /= hdr0%md5_pseudos(ipsp)) then
4436 : write(msg, '(a,i0,6a)' )&
4437 3 : ' Different md5 checksum for pseudo ',ipsp,ch10,&
4438 3 : ' input md5= ',hdr%md5_pseudos(ipsp),ch10,&
4439 6 : ' disk md5= ',hdr0%md5_pseudos(ipsp)
4440 3 : ABI_WARNING(msg)
4441 3 : itest=1; tpsch=1
4442 : end if
4443 :
4444 : ! Second, test
4445 : ! NOTE, XG 000719: should do something about pspso
4446 : ! NOTE, XG 020716: znucl and zion are not written
4447 6050 : if (abs(hdr%znuclpsp(ipsp)-hdr0%znuclpsp(ipsp))>tol6) itest=1
4448 6050 : if (abs(hdr%zionpsp(ipsp)-hdr0%zionpsp(ipsp))>tol6) then
4449 0 : itest=1; tpsch=1
4450 : end if
4451 6050 : if (hdr%pspdat(ipsp)/= hdr0%pspdat(ipsp)) itest=1
4452 6050 : if (hdr%pspcod(ipsp)/= hdr0%pspcod(ipsp)) itest=1
4453 10285 : if (hdr%pspxc(ipsp) /= hdr0%pspxc(ipsp) ) itest=1
4454 : end do
4455 :
4456 4235 : if (itest==1) then
4457 3 : ABI_WARNING('input psp header does not agree perfectly with disk file psp header.')
4458 3 : tpseu=1
4459 : end if
4460 : end if
4461 :
4462 : ! Finally, read residm and etotal ("current value" not known), and check xred.
4463 4235 : if (hdr%natom==hdr0%natom) then
4464 4235 : natom=hdr%natom
4465 4235 : write(msg,'(a,37x,a,a)') ' xred:','|',' xred:'
4466 4235 : call wrtout(std_out,msg,mode_paral)
4467 13080 : do ii=1,natom
4468 8845 : write(msg,'(2x,3f12.7,6x,a,2x,3f12.7)') hdr%xred(:,ii),'|',hdr0%xred(:,ii)
4469 13080 : call wrtout(std_out,msg,mode_paral)
4470 : end do
4471 :
4472 : ! check atom positions one atom at a time and allow possibility
4473 : ! that there is a harmless translation of atoms by a cell vector.
4474 13080 : do ii=1,natom
4475 8845 : rms=0.0_dp
4476 35380 : do jj=1,3
4477 35380 : rms=rms+(hdr%xred(jj,ii)-hdr0%xred(jj,ii) - dble(nint((hdr%xred(jj,ii)-hdr0%xred(jj,ii)))) )**2
4478 : end do
4479 8845 : rms=sqrt(rms/3.0_dp)
4480 13080 : if (rms>tol6) txred=1
4481 : end do
4482 : end if
4483 :
4484 : ! Run tests here to establish whether this is a valid restart
4485 :
4486 : ! tfform2 will be true if there is a problem for the wavefunctions
4487 : tfform2 = (hdr%usewvl == 0 .and. &
4488 : (tprim /= 0 .or. tecut /= 0 .or. tkpt /= 0 .or. &
4489 : twfk /=0 .or. tspinor /= 0)) .or. &
4490 : (hdr%usepaw == 1 .and. &
4491 : (tpaw /= 0 .or. tlmn /= 0 .or. tdg /= 0)) .or. &
4492 : (hdr%usewvl == 1 .and. &
4493 4235 : (tatty /= 0 .or. tband /= 0))
4494 :
4495 : ! tfform52 will be true if there is a problem for the format 52
4496 : tfform52 = tprim /= 0 .or. tatty /= 0 .or. txred /= 0 .or.&
4497 : tpseu /= 0 .or. tecut /= 0 .or. tng /= 0 .or. &
4498 4235 : (hdr%usepaw == 1 .and. (tpaw /= 0 .or. tlmn /= 0 .or. tdg /= 0))
4499 :
4500 4235 : restart=1; restartpaw=hdr%usepaw
4501 :
4502 : ! If there is a problem somewhere
4503 : if ( (abifile%class == "wf_planewave" .and. tfform2 ) .or. &
4504 4235 : (abifile%class == "density" .and. tfform52 ) .or. &
4505 : (abifile%class == "wf_wavelet" .and. tfform2 ) ) then
4506 :
4507 518 : if (abifile%class == "wf_planewave") then
4508 517 : restart=2
4509 517 : ABI_COMMENT('Restart of self-consistent calculation need translated wavefunctions.')
4510 1 : else if (abifile%class == "density") then
4511 1 : restart=0
4512 1 : ABI_WARNING('Illegal restart of non-self-consistent calculation')
4513 : end if
4514 :
4515 : write(msg,'(a,a1,a)') &
4516 518 : ' Indeed, critical differences between current calculation and',ch10,&
4517 1036 : ' restart file have been detected in:'
4518 518 : call wrtout(std_out,msg,mode_paral)
4519 :
4520 518 : if ( (abifile%class == "density" .or. abifile%class == "wf_wavelet") .and. tatty /= 0 ) then
4521 0 : write(msg, '(8x,a)' ) '* the number of atoms of each type'
4522 0 : call wrtout(std_out,msg,mode_paral)
4523 : end if
4524 518 : if ( abifile%class /= "wf_wavelet" .and. tecut /= 0 ) then
4525 28 : write(msg, '(8x,a)' ) '* the plane-wave cutoff'
4526 28 : call wrtout(std_out,msg,mode_paral)
4527 : end if
4528 518 : if ( abifile%class == "wf_wavelent" .and. tband /= 0 ) then
4529 0 : write(msg, '(8x,a)' ) '* the band and their occupation'
4530 0 : call wrtout(std_out,msg,mode_paral)
4531 : end if
4532 518 : if ( abifile%class == "wf_planewave" .and. tkpt /= 0 ) then
4533 248 : write(msg, '(8x,a)' ) '* the number, position, or weight of k-points'
4534 248 : call wrtout(std_out,msg,mode_paral)
4535 : end if
4536 518 : if ( abifile%class == "wf_planewave" .and. twfk /= 0 ) then
4537 369 : write(msg, '(8x,a)' ) '* the format of wavefunctions (istwfk)'
4538 369 : call wrtout(std_out,msg,mode_paral)
4539 : end if
4540 518 : if ( abifile%class == "wf_planewave" .and. tspinor /= 0 ) then
4541 13 : write(msg, '(8x,a)' ) '* the scalar/spinor character of the wf (nspinor)'
4542 13 : call wrtout(std_out,msg,mode_paral)
4543 : end if
4544 518 : if ( abifile%class == "density" .and. tng /= 0 ) then
4545 0 : write(msg, '(8x,a)' ) '* the Fourier transform box dimensions'
4546 0 : call wrtout(std_out,msg,mode_paral)
4547 : end if
4548 518 : if ( tprim /= 0 ) then
4549 74 : write(msg, '(8x,a)' )'* the vectors defining the unit cell (obtained from rprim and acell)'
4550 74 : call wrtout(std_out,msg,mode_paral)
4551 : end if
4552 518 : if ( abifile%class == "density" .and. tpseu /= 0 ) then
4553 1 : write(msg, '(8x,a)' )'* the pseudopotential files'
4554 1 : call wrtout(std_out,msg,mode_paral)
4555 : end if
4556 518 : if ( abifile%class == "density" .and. txred /= 0 ) then
4557 0 : write(msg, '(8x,a)' ) '* the positions of the ions in the basis'
4558 0 : call wrtout(std_out,msg,mode_paral)
4559 : end if
4560 :
4561 : ! Tests for a restart in the framework of the PAW method
4562 518 : if (hdr%usepaw/=0 .or. hdr0%usepaw/=0) then
4563 56 : if (tpaw /= 0 .or. tlmn /= 0) restartpaw=0
4564 56 : if (restartpaw == 0) then
4565 0 : write(msg,'(8x,a)') 'Critical differences for a restart within PAW method:'
4566 0 : call wrtout(std_out,msg,mode_paral)
4567 0 : if ( tpaw /= 0 ) then
4568 0 : write(msg, '(8x,a)' ) '* the use of the PAW method'
4569 0 : call wrtout(std_out,msg,mode_paral)
4570 : else
4571 0 : if(tlmn/=0)then
4572 0 : write(msg, '(8x,a)' ) '* the number of lmn elements for the paw basis'
4573 0 : call wrtout(std_out,msg,mode_paral)
4574 : end if
4575 : end if
4576 56 : else if (tdg/=0) then
4577 14 : write(msg,'(a,a,a,a,a,a)') ch10,&
4578 14 : ' hdr_check: WARNING -',ch10,&
4579 28 : ' Restart of calculation within PAW may be inconsistent because of:"'
4580 14 : call wrtout(std_out,msg,mode_paral)
4581 : if(tdg/=0)then
4582 14 : write(msg, '(8x,a)' )'* the cutoff energy of the paw double (fine) grid'
4583 14 : call wrtout(std_out,msg,mode_paral)
4584 : end if
4585 : end if
4586 : end if
4587 :
4588 : else
4589 :
4590 3717 : if (abifile%class == "wf_planewave" .or. abifile%class == "wf_wavelet") then
4591 1755 : write(msg,'(a,a)') ' hdr_check: ',' Wavefunction file is OK for direct restart of calculation'
4592 1755 : call wrtout(std_out,msg,mode_paral)
4593 1962 : else if (abifile%class == "density") then
4594 1957 : write(msg,'(a,a)') ' hdr_check: ',' Density/Potential file is OK for restart of calculation'
4595 1957 : call wrtout(std_out,msg,mode_paral)
4596 : end if
4597 : end if
4598 :
4599 343035 : write(msg,'(80a)') ('=',ii=1,80)
4600 4235 : call wrtout(std_out,msg,mode_paral)
4601 :
4602 4235 : end subroutine hdr_check
4603 : !!***
4604 :
4605 : !----------------------------------------------------------------------
4606 :
4607 : !!****f* m_wfk/hdr_compare
4608 : !! NAME
4609 : !! hdr_compare
4610 : !!
4611 : !! FUNCTION
4612 : !! Test two hdr_t objects for consistency. Return non-zero value if test fails.
4613 : !!
4614 : !! INPUTS
4615 : !! hdr1, hdr2 <class(hdr_t)> = hdr handlers to be compared
4616 : !!
4617 : !! OUTPUT
4618 : !! ierr
4619 : !!
4620 : !! SOURCE
4621 :
4622 29 : integer function hdr_compare(hdr1, hdr2) result(ierr)
4623 :
4624 : !Arguments ------------------------------------
4625 : !scalars
4626 : class(hdr_type),intent(in) :: hdr1, hdr2
4627 :
4628 : !Local variables-------------------------------
4629 : character(len=500) :: msg
4630 : !************************************************************************
4631 :
4632 29 : ierr = 0
4633 :
4634 : ! Test basic dimensions
4635 29 : if (hdr1%nsppol /= hdr2%nsppol) then
4636 0 : write(msg,'(a,i0,a,i0)')'Different nsppol : ',hdr1%nsppol,' and ',hdr2%nsppol
4637 0 : ierr = ierr + 1; ABI_WARNING(msg)
4638 : end if
4639 29 : if (hdr1%nspinor /= hdr2%nspinor) then
4640 0 : write(msg,'(a,i0,a,i0)')'Different nspinor : ',hdr1%nspinor,' and ',hdr2%nspinor
4641 0 : ierr = ierr + 1; ABI_WARNING(msg)
4642 : end if
4643 29 : if (hdr1%nspden /= hdr2%nspden) then
4644 0 : write(msg,'(a,i0,a,i0)')'Different nspden : ',hdr1%nspden,' and ',hdr2%nspden
4645 0 : ierr = ierr + 1; ABI_WARNING(msg)
4646 : end if
4647 29 : if (hdr1%nkpt /= hdr2%nkpt) then
4648 0 : write(msg,'(a,i0,a,i0)')'Different nkpt : ',hdr1%nkpt,' and ',hdr2%nkpt
4649 0 : ierr = ierr + 1; ABI_WARNING(msg)
4650 : end if
4651 29 : if (hdr1%usepaw /= hdr2%usepaw) then
4652 0 : write(msg,'(a,i0,a,i0)')'Different usepaw : ',hdr1%usepaw,' and ',hdr2%usepaw
4653 0 : ierr = ierr + 1; ABI_WARNING(msg)
4654 : end if
4655 29 : if (hdr1%ntypat /= hdr2%ntypat) then
4656 0 : write(msg,'(a,i0,a,i0)')'Different ntypat : ',hdr1%ntypat,' and ',hdr2%ntypat
4657 0 : ierr = ierr + 1; ABI_WARNING(msg)
4658 : end if
4659 29 : if (hdr1%natom /= hdr2%natom) then
4660 0 : write(msg,'(a,i0,a,i0)')'Different natom : ',hdr1%natom,' and ',hdr2%natom
4661 0 : ierr = ierr + 1; ABI_WARNING(msg)
4662 : end if
4663 :
4664 : ! Return immediately if important dimensions are not equal.
4665 29 : if (ierr /= 0) return
4666 :
4667 : ! Test important arrays (rprimd is not tested)
4668 79 : if (any(hdr1%typat /= hdr2%typat)) then
4669 0 : write(msg,'(a,i0,a,i0)')'Different ntypat array : ',hdr1%typat(1),' ... and ',hdr2%typat(1)
4670 0 : ierr = ierr + 1; ABI_WARNING(msg)
4671 : end if
4672 : !Should test npwarr, however taking into account differences due to istwfk !
4673 : !if (any(hdr1%npwarr /= hdr2%npwarr)) then
4674 : ! write(msg,'(a,i0,a,i0)')'Different npwarr array : ',hdr1%npwarr(1),' ... and ',hdr2%npwarr(1)
4675 : ! ierr = ierr + 1; ABI_WARNING(msg)
4676 : !end if
4677 9157 : if (any(abs(hdr1%kptns - hdr2%kptns) > tol6)) then
4678 0 : write(msg,'(a,i0,a,i0)')'Different kptns array '
4679 0 : ierr = ierr + 1; ABI_WARNING(msg)
4680 : end if
4681 :
4682 : end function hdr_compare
4683 : !!***
4684 :
4685 : !----------------------------------------------------------------------
4686 :
4687 : !!****f* m_hdr/hdr_vs_dtset
4688 : !! NAME
4689 : !! hdr_vs_dtset
4690 : !!
4691 : !! FUNCTION
4692 : !! Check the compatibility of the Abinit header with respect to the
4693 : !! input variables defined in the input file.
4694 : !!
4695 : !! INPUTS
4696 : !! Dtset<type(dataset_type)>=all input variables for this dataset
4697 : !! Hdr <type(hdr_type)>=the header structured variable
4698 : !!
4699 : !! OUTPUT
4700 : !! Only check
4701 : !!
4702 : !! SOURCE
4703 :
4704 527 : subroutine hdr_vs_dtset(hdr, dtset)
4705 :
4706 : !Arguments ------------------------------------
4707 : class(Hdr_type),intent(in) :: hdr
4708 : type(Dataset_type),intent(in) :: dtset
4709 :
4710 : !Local variables-------------------------------
4711 : integer :: ik, jj, ierr
4712 : logical :: test, tsymrel,ttnons, tsymafm
4713 : character(len=5000) :: msg
4714 : ! *************************************************************************
4715 :
4716 : ! Check basic dimensions
4717 527 : ierr = 0
4718 527 : call compare_int('natom', hdr%natom, dtset%natom, ierr)
4719 527 : call compare_int('nkpt', hdr%nkpt, dtset%nkpt, ierr)
4720 527 : call compare_int('npsp', hdr%npsp, dtset%npsp, ierr)
4721 527 : call compare_int('nspden', hdr%nspden, dtset%nspden, ierr)
4722 527 : call compare_int('nspinor',hdr%nspinor,dtset%nspinor,ierr)
4723 527 : call compare_int('nsppol', hdr%nsppol, dtset%nsppol, ierr)
4724 527 : call compare_int('nsym', hdr%nsym, dtset%nsym, ierr)
4725 527 : call compare_int('ntypat', hdr%ntypat, dtset%ntypat, ierr)
4726 527 : call compare_int('usepaw', hdr%usepaw, dtset%usepaw, ierr)
4727 527 : call compare_int('usewvl', hdr%usewvl, dtset%usewvl, ierr)
4728 527 : call compare_int('kptopt', hdr%kptopt, dtset%kptopt, ierr)
4729 527 : call compare_int('pawcpxocc', hdr%pawcpxocc, dtset%pawcpxocc, ierr)
4730 527 : call compare_int('nshiftk_orig', hdr%nshiftk_orig, dtset%nshiftk_orig, ierr)
4731 527 : call compare_int('nshiftk', hdr%nshiftk, dtset%nshiftk, ierr)
4732 : !call compare_int("ixc", hdr%ixc, dtset%ixc, ierr)
4733 :
4734 : ! The number of fatal errors must be zero.
4735 527 : if (ierr /= 0) then
4736 : write(msg,'(3a)')&
4737 0 : 'Cannot continue, basic dimensions/parameters reported in the header do not agree with input file. ',ch10,&
4738 0 : 'Check consistency between the content of the WFK external file and the input file.'
4739 0 : ABI_ERROR(msg)
4740 : end if
4741 :
4742 5099 : test = ALL(ABS(Hdr%xred - Dtset%xred_orig(:,1:Dtset%natom,1)) < tol3)
4743 527 : ABI_CHECK(test, 'Mismatch in xred')
4744 :
4745 1670 : test=ALL(Hdr%typat == Dtset%typat(1:Dtset%natom))
4746 527 : ABI_CHECK(test,'Mismatch in typat')
4747 :
4748 : ! Check if the lattice from the input file agrees with the one read from the WFK file
4749 6851 : if ( (ANY(ABS(Hdr%rprimd - Dtset%rprimd_orig(1:3,1:3,1)) > tol6)) ) then
4750 0 : write(msg,'(5a,3(3es16.6),3a,3(3es16.6),3a)')ch10,&
4751 0 : ' real lattice vectors read from Header differ from the values specified in the input file', ch10, &
4752 0 : ' rprimd from Hdr file = ',ch10,(Hdr%rprimd(:,jj),jj=1,3),ch10,&
4753 0 : ' rprimd from input file = ',ch10,(Dtset%rprimd_orig(:,jj,1),jj=1,3),ch10,ch10,&
4754 0 : ' Modify the lattice vectors in the input file '
4755 0 : ABI_ERROR(msg)
4756 : end if
4757 :
4758 : ! Check symmetry operations.
4759 283277 : tsymrel=(ALL(Hdr%symrel==Dtset%symrel(:,:,1:Dtset%nsym)))
4760 527 : if (.not.tsymrel) then
4761 : write(msg,'(3a)')&
4762 0 : ' real space symmetries read from Header ',ch10,&
4763 0 : ' differ from the values inferred from the input file'
4764 0 : ABI_WARNING(msg)
4765 0 : tsymrel=.FALSE.
4766 : end if
4767 :
4768 87527 : ttnons=ALL(ABS(Hdr%tnons-Dtset%tnons(:,1:Dtset%nsym))<tol6)
4769 527 : if (.not.ttnons) then
4770 : write(msg,'(3a)')&
4771 0 : ' fractional translations read from Header ',ch10,&
4772 0 : ' differ from the values inferred from the input file'
4773 0 : ABI_WARNING(msg)
4774 0 : ttnons=.FALSE.
4775 : end if
4776 :
4777 22277 : tsymafm=ALL(Hdr%symafm==Dtset%symafm(1:Dtset%nsym))
4778 527 : if (.not.tsymafm) then
4779 : write(msg,'(3a)')&
4780 0 : ' AFM symmetries read from Header ',ch10,&
4781 0 : ' differ from the values inferred from the input file'
4782 0 : ABI_WARNING(msg)
4783 0 : tsymafm=.FALSE.
4784 : end if
4785 :
4786 527 : if (.not. (tsymrel.and.ttnons.and.tsymafm)) then
4787 0 : write(msg,'(a)')' Header '
4788 0 : call wrtout(std_out,msg)
4789 0 : call print_symmetries([std_out], Hdr%nsym,Hdr%symrel,Hdr%tnons,Hdr%symafm)
4790 0 : write(msg,'(a)')' Dtset '
4791 0 : call wrtout(std_out,msg)
4792 0 : call print_symmetries([std_out], Dtset%nsym,Dtset%symrel,Dtset%tnons,Dtset%symafm)
4793 0 : ABI_ERROR('Check symmetry operations')
4794 : end if
4795 :
4796 527 : if (abs(Dtset%nelect-hdr%nelect)>tol6) then
4797 0 : write(msg,'(2(a,f8.2))')"File contains ", hdr%nelect," electrons but nelect initialized from input is ",Dtset%nelect
4798 0 : ABI_ERROR(msg)
4799 : end if
4800 :
4801 527 : if (abs(Dtset%ne_qFD-hdr%ne_qFD)>tol6) then
4802 0 : write(msg,'(2(a,f8.2))')"File contains ", hdr%ne_qFD, &
4803 0 : " electrons in the conduction bands but nelect initialized from input is ",Dtset%ne_qFD
4804 0 : ABI_ERROR(msg)
4805 : end if
4806 :
4807 527 : if (abs(Dtset%nh_qFD-hdr%nh_qFD)>tol6) then
4808 0 : write(msg,'(2(a,f8.2))')"File contains ", hdr%nh_qFD,&
4809 0 : " electrons in the valence bands but nelect initialized from input is ",Dtset%nh_qFD
4810 0 : ABI_ERROR(msg)
4811 : end if
4812 :
4813 527 : if (abs(Dtset%cellcharge(1)-hdr%cellcharge)>tol6) then
4814 0 : write(msg,'(2(a,f8.2))')"File contains cellcharge ", hdr%cellcharge," but cellcharge from input is ",Dtset%cellcharge
4815 0 : ABI_ERROR(msg)
4816 : end if
4817 :
4818 6851 : if (any(hdr%kptrlatt_orig /= dtset%kptrlatt_orig)) then
4819 : write(msg,"(5a)")&
4820 0 : "hdr%kptrlatt_orig: ",trim(ltoa(reshape(hdr%kptrlatt_orig,[9]))),ch10,&
4821 0 : "dtset%kptrlatt_orig: ",trim(ltoa(reshape(dtset%kptrlatt_orig, [9])))
4822 0 : ABI_ERROR(msg)
4823 : end if
4824 :
4825 6851 : if (any(hdr%kptrlatt /= dtset%kptrlatt)) then
4826 : write(msg,"(5a)") &
4827 0 : "hdr%kptrlatt: ",trim(ltoa(reshape(hdr%kptrlatt, [9]))),ch10,&
4828 0 : "dtset%kptrlatt: ",trim(ltoa(reshape(dtset%kptrlatt, [9])))
4829 0 : ABI_ERROR(msg)
4830 : end if
4831 :
4832 4363 : if (any(abs(hdr%shiftk_orig - dtset%shiftk_orig(:,1:dtset%nshiftk_orig)) > tol6)) then
4833 : write(msg,"(5a)")&
4834 0 : "hdr%shiftk_orig: ",trim(ltoa(reshape(hdr%shiftk_orig, [3*hdr%nshiftk_orig]))),ch10,&
4835 0 : "dtset%shiftk_orig: ",trim(ltoa(reshape(dtset%shiftk_orig, [3*dtset%nshiftk_orig])))
4836 0 : ABI_ERROR(msg)
4837 : end if
4838 :
4839 2635 : if (any(abs(hdr%shiftk - dtset%shiftk(:,1:dtset%nshiftk)) > tol6)) then
4840 : write(msg,"(5a)")&
4841 0 : "hdr%shiftk: ",trim(ltoa(reshape(hdr%shiftk, [3*hdr%nshiftk]))),ch10,&
4842 0 : "dtset%shiftk: ",trim(ltoa(reshape(dtset%shiftk, [3*dtset%nshiftk])))
4843 0 : ABI_ERROR(msg)
4844 : end if
4845 :
4846 : ! Check if the k-points from the input file agrees with that read from the WFK file
4847 : !
4848 : ! Note that the EPH code frees several dtset arrays that depend on nkpt (see dtset_free_nkpt_arrays)
4849 : ! in order to reduce memory when large meshes are used (e.g. 200x200x200)
4850 : ! so we have to test whether the dtset array is allocated before testing.
4851 : !
4852 527 : if (allocated (dtset%kpt)) then
4853 36395 : if ((any(abs(Hdr%kptns(:,:) - Dtset%kpt(:,1:Dtset%nkpt)) > tol6))) then
4854 0 : write(msg,'(9a)')ch10,&
4855 0 : ' hdr_vs_dtset: ERROR - ',ch10,&
4856 0 : ' k-points read from Header ',ch10,&
4857 0 : ' differ from the values specified in the input file',ch10,&
4858 0 : ' k-points from Hdr file | k-points from input file ',ch10
4859 0 : call wrtout(std_out,msg)
4860 0 : do ik=1,Dtset%nkpt
4861 0 : if (any(abs(Hdr%kptns(:,ik) - Dtset%kpt(:,ik)) > tol6)) then
4862 0 : write(msg,'(3(3es16.6,3x))')Hdr%kptns(:,ik),Dtset%kpt(:,ik)
4863 0 : call wrtout(std_out,msg)
4864 : end if
4865 : end do
4866 0 : ABI_ERROR('Modify the k-mesh in the input file')
4867 : end if
4868 : end if
4869 :
4870 527 : if (allocated(dtset%wtk)) then
4871 9482 : if (ANY(ABS(Hdr%wtk(:) - Dtset%wtk(1:Dtset%nkpt)) > tol6)) then
4872 0 : write(msg,'(9a)')ch10,&
4873 0 : ' hdr_vs_dtset : ERROR - ',ch10,&
4874 0 : ' k-point weights read from Header ',ch10,&
4875 0 : ' differ from the values specified in the input file',ch10,&
4876 0 : ' Hdr file | File ',ch10
4877 0 : call wrtout(std_out,msg)
4878 0 : do ik=1,Dtset%nkpt
4879 0 : if (abs(Hdr%wtk(ik) - Dtset%wtk(ik)) > tol6) then
4880 0 : write(msg,'(2(f11.5,1x))')Hdr%wtk(ik),Dtset%wtk(ik)
4881 0 : call wrtout(std_out,msg)
4882 : end if
4883 : end do
4884 0 : ABI_ERROR('Check the k-mesh and the symmetries of the system. ')
4885 : end if
4886 : end if
4887 :
4888 : ! Check istwfk storage
4889 527 : if (allocated(dtset%istwfk)) then
4890 9308 : if ((any(Hdr%istwfk(:) /= Dtset%istwfk(1:Dtset%nkpt))) ) then
4891 2 : ABI_COMMENT('istwfk from header differs from the values specified in the input file (not critical)')
4892 : !call wrtout(std_out, " Hdr | input ")
4893 : !do ik=1,Dtset%nkpt
4894 : ! write(msg,'(i5,3x,i5)')Hdr%istwfk(ik),Dtset%istwfk(ik)
4895 : ! call wrtout(std_out,msg)
4896 : !end do
4897 : !ABI_ERROR('Modify istwfk in the input file.')
4898 : end if
4899 : end if
4900 :
4901 : contains
4902 : !!***
4903 :
4904 : !!****f* hdr_vs_dtset/compare_int
4905 : !! NAME
4906 : !! compare_int
4907 : !!
4908 : !! FUNCTION
4909 : !! Compare two int value and may raise an exception on error.
4910 : !!
4911 : !! INPUTS
4912 : !! vname=Name of the variable
4913 : !! iexp= expected value.
4914 : !! ifound=the actuval value
4915 : !!
4916 : !! SIDE EFFECTS
4917 : !! ierr=increased by one if values differ
4918 : !!
4919 : !! SOURCE
4920 :
4921 7378 : subroutine compare_int(vname, iexp, ifound, ierr)
4922 :
4923 : !Arguments ------------------------------------
4924 : integer,intent(in) :: iexp,ifound
4925 : integer,intent(inout) :: ierr
4926 : character(len=*),intent(in) :: vname
4927 :
4928 : !Local variables-------------------------------
4929 : character(len=500) :: msg
4930 : ! *************************************************************************
4931 :
4932 7378 : if (.not. iexp == ifound) then
4933 0 : write(msg,'(2a,i0,a,i0)')' Mismatch in '//trim(vname),': Expected = ', iexp, ' Found = ', ifound
4934 0 : call wrtout(std_out, msg)
4935 : ! Increase ierr to signal we should stop in the caller.
4936 0 : ierr = ierr + 1
4937 : end if
4938 :
4939 7378 : end subroutine compare_int
4940 : !!***
4941 :
4942 : end subroutine hdr_vs_dtset
4943 : !!***
4944 :
4945 : !!****f* m_hdr/hdr_get_crystal
4946 : !! NAME
4947 : !! hdr_get_crystal
4948 : !!
4949 : !! FUNCTION
4950 : !! Initializes a crystal_t data type starting from the abinit header.
4951 : !!
4952 : !! INPUTS
4953 : !! hdr<hdr_type>=the abinit header
4954 : !! [gw_timrev] ==2 => take advantage of time-reversal symmetry
4955 : !! ==1 ==> do not use time-reversal symmetry
4956 : !! Default: 2
4957 : !! NOTE THAT HERE WE USE THE GW CONVENTIONS I.E ABINIT_TIMREV + !
4958 : !! [remove_inv] = if .TRUE. the inversion symmetry is removed from the set of operations
4959 : !! even if it is present in the header
4960 : !!
4961 : !! OUTPUT
4962 : !! cryst<crystal_t>= the data type filled with data reported in the abinit header
4963 : !!
4964 : !! TODO
4965 : !! Add information on the use of time-reversal in the Abinit header.
4966 : !!
4967 : !! SOURCE
4968 :
4969 11498 : type(crystal_t) function hdr_get_crystal(hdr, gw_timrev, remove_inv) result(cryst)
4970 :
4971 : !Arguments ------------------------------------
4972 : class(hdr_type),intent(in) :: hdr
4973 : integer,optional,intent(in) :: gw_timrev
4974 : logical,optional,intent(in) :: remove_inv
4975 :
4976 : !Local variables-------------------------------
4977 : integer :: my_timrev, space_group
4978 : logical :: rinv, use_antiferro
4979 : ! *********************************************************************
4980 :
4981 386 : rinv=.FALSE.; if (PRESENT(remove_inv)) rinv=remove_inv
4982 11498 : use_antiferro = hdr%nspden == 2 .and. hdr%nsppol ==1
4983 :
4984 11498 : if (.not. present(gw_timrev)) then
4985 : ! Get it from kptopt
4986 : !my_timrev = kpts_timrev_from_kptopt(hdr%kptopt) + 1
4987 7228 : my_timrev = 1; if (any(hdr%kptopt == [3, 4])) my_timrev = 0
4988 7228 : my_timrev = my_timrev + 1
4989 : !print *, "my_timrev", my_timrev
4990 : else
4991 4270 : my_timrev = gw_timrev
4992 : end if
4993 :
4994 : ! Consistency check
4995 11498 : ABI_CHECK(any(my_timrev == [1, 2]), "timrev should be in (1|2)")
4996 11498 : if (use_antiferro) then
4997 78 : ABI_CHECK(ANY(hdr%symafm == -1), "Wrong nspden, nsppol, symafm.")
4998 : end if
4999 :
5000 11498 : space_group = 0 ! FIXME not known at this level.
5001 :
5002 : call cryst%init(hdr%amu,space_group,hdr%natom,hdr%npsp,hdr%ntypat,hdr%nsym,hdr%rprimd,hdr%typat,hdr%xred,&
5003 : hdr%zionpsp,hdr%znuclpsp,my_timrev,use_antiferro,rinv,hdr%title,&
5004 11498 : symrel=hdr%symrel,tnons=hdr%tnons,symafm=hdr%symafm) ! Optional
5005 :
5006 643888 : end function hdr_get_crystal
5007 : !!***
5008 :
5009 117948 : end module m_hdr
5010 : !!***
|
