Line data Source code
1 : !!****m* ABINIT/m_ddb_hdr
2 : !! NAME
3 : !! m_ddb_hdr
4 : !!
5 : !! FUNCTION
6 : !! This module contains the declaration of data types and methods
7 : !! to handle the header of the DDB files.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2011-2026 ABINIT group (GA, MJV, XG, MT, MM, MVeithen, MG, PB, JCC, MMignolet)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_ddb_hdr
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_xmpi
29 : use m_dtset
30 : use m_crystal
31 : use m_nctk
32 : use netcdf
33 :
34 : use defs_datatypes, only : pseudopotential_type
35 : use m_copy, only : alloc_copy
36 : use m_pawtab, only : pawtab_type, pawtab_nullify, pawtab_free, pawtab_bcast !, pawtab_copy
37 : use m_psps, only : psps_copy, psps_free, psps_ncwrite, psps_ncread
38 : use m_io_tools, only : open_file, get_unit, file_exists
39 : use m_copy, only : alloc_copy
40 : use m_fstrings, only : sjoin, endswith
41 : use m_geometry, only : mkrdim
42 :
43 : implicit none
44 :
45 : private
46 :
47 : public :: is_type_d0E ! Is this block type of a d0E kind?
48 : public :: is_type_d1E ! Is this block type of a d1E kind?
49 : public :: is_type_d2E ! Is this block type of a d2E kind?
50 : public :: is_type_d3E ! Is this block type of a d3E kind?
51 : public :: is_type_d2eig ! Is this block type of a d2eig kind?
52 :
53 : ! GA: These 3 functions should no longer be public
54 : public :: ddb_getdims ! Open a DDB file and read basic dimensions and variables.
55 : public :: ioddb8_in ! Temporary
56 : public :: psddb8 ! Temporary
57 :
58 :
59 : ! Description of perturbations and block types
60 : ! --------------------------------------------
61 : !
62 : ! GA: Variable names should not end with _d0E, _d1E, d2E, ...
63 : ! Such names may affect the variable type, with certain compilers.
64 : ! This is why I appended _xx, lack of a better suffix.
65 : !
66 : !
67 : ! These parameters should be consistent with input variable rfmeth
68 : ! and input parameter rftyp (e.g. ddb_get_block, ddb_read_block_txt).
69 : ! Changing these values requires to change the documentation as well.
70 : !
71 : integer,public,parameter :: BLKTYP_d0E_xx=0 ! Total energy
72 : integer,public,parameter :: BLKTYP_d1E_xx=4 ! First-order derivatives of total energy
73 : integer,public,parameter :: BLKTYP_d2E_ns=1 ! Second-order derivatives of total energy, non-stationary
74 : integer,public,parameter :: BLKTYP_d2E_st=2 ! Second-order derivatives of total energy, stationary
75 : integer,public,parameter :: BLKTYP_d2E_mbc=85 ! Molecular Berry curvature (MBC)
76 : integer,public,parameter :: BLKTYP_d3E_xx=3 ! Third-order derivatives of total energy
77 : integer,public,parameter :: BLKTYP_d3E_lw=33 ! Long-wave third-order derivatives of total energy
78 : integer,public,parameter :: BLKTYP_d2eig_re=5 ! Second-order derivatives of eigenvalues
79 : integer,public,parameter :: BLKTYP_d2eig_im=6 ! Static broadening of eigenvalues
80 :
81 : integer,public,parameter :: descrlen=500
82 : integer,public,parameter :: ntypes=9 ! The number of different block types
83 : character(len=descrlen),public,parameter :: DESCR_d0E_xx = ' Total energy - '
84 : character(len=descrlen),public,parameter :: DESCR_d1E_xx = ' 1st derivatives - '
85 : character(len=descrlen),public,parameter :: DESCR_d2E_ns = ' 2nd derivatives (non-stat.) - '
86 : character(len=descrlen),public,parameter :: DESCR_d2E_st = ' 2nd derivatives (stationary) - '
87 : character(len=descrlen),public,parameter :: DESCR_d2E_mbc = ' 2nd derivatives (MBC) - '
88 : character(len=descrlen),public,parameter :: DESCR_d3E_xx = ' 3rd derivatives - '
89 : character(len=descrlen),public,parameter :: DESCR_d3E_lw = ' 3rd derivatives (long wave) - '
90 : character(len=descrlen),public,parameter :: DESCR_d2eig_re = ' 2nd eigenvalue derivatives - '
91 : character(len=descrlen),public,parameter :: DESCR_d2eig_im = ' 2nd eigenvalue derivatives (imaginary part) - '
92 :
93 : ! These describe the perturbation indices above natom
94 : ! GA: It would be nice if different ddb instance could order perturbations
95 : ! differently. This could potentially reduce the ddb size.
96 : character(len=descrlen),public,parameter :: DESCR_ipert_0 = 'displacement of atom '
97 : character(len=descrlen),public,parameter :: DESCR_ipert_1 = 'derivative wrt k'
98 : character(len=descrlen),public,parameter :: DESCR_ipert_2 = 'electric field'
99 : character(len=descrlen),public,parameter :: DESCR_ipert_3 = 'strain' ! diagonal stress components
100 : character(len=descrlen),public,parameter :: DESCR_ipert_4 = 'strain' ! non-diagonal stress
101 : character(len=descrlen),public,parameter :: DESCR_ipert_5 = 'magnetic field'
102 : character(len=descrlen),public,parameter :: DESCR_ipert_10 = '2nd derivative wrt to k'
103 : character(len=descrlen),public,parameter :: DESCR_ipert_11 = '2nd derivative wrt to k and electric field'
104 :
105 : integer,public,parameter :: DDB_VERSION=20240201 !
106 : !integer,public,parameter :: DDB_VERSION=20230401 ! TODO: check if we should update this with new G matrix stuff
107 : ! DDB Version number for text format.
108 :
109 : integer,public,parameter :: DDB_VERSION_NC=20240201 ! TODO:
110 : !integer,public,parameter :: DDB_VERSION_NC=20230219 ! TODO: check if we should update this with new G matrix stuff
111 : ! DDB NetCDF version number.
112 :
113 : type,public :: ddb_hdr_type
114 :
115 : logical :: has_open_file_txt=.false. ! Has an open file in text format
116 : logical :: has_open_file_nc=.false. ! Has an open file in netcdf format
117 :
118 : integer :: iomode=IO_MODE_FORTRAN
119 :
120 : integer :: unddb ! Unit for open file in text format
121 : integer :: ncid ! Unit for open file in netcdf format
122 :
123 : ! The variables with_psps / with_dfpt_vars
124 : ! Must be identical in text format, but may differ in netcdf.
125 : integer :: with_psps=0 ! Whether info on pseudopotentials is present.
126 : integer :: with_dfpt_vars=0 ! Whether this ddb comes from a dfpt calculation
127 : ! 0 -> comes from a ground state calculation.
128 : ! 1 -> comes from a dfpt calculation
129 : ! and k-point may or may not be reduced
130 : ! by time reversal symmetry.
131 :
132 : integer :: ddb_version ! Version of the DDB file
133 :
134 : integer :: matom
135 : integer :: mband
136 : integer :: mkpt
137 : integer :: msym ! GA: Do we need both msym and nsym ?
138 : ! ddb_hdr reading is first called with msym=192
139 : ! for dimensions only, then it is called with
140 : ! actual number of sym.
141 : ! See also comment in dfpt_phfrq
142 : ! TODO: try to eliminate msym as much as possible.
143 : integer :: mtypat
144 : integer :: mpert
145 : integer :: msize
146 :
147 : integer :: intxc
148 : integer :: iscf
149 : integer :: ixc
150 : integer :: natom
151 : integer :: nkpt
152 : integer :: nspden
153 : integer :: nspinor
154 : integer :: nsppol
155 : integer :: nsym
156 : integer :: ntypat
157 : integer :: occopt
158 : integer :: usepaw
159 :
160 : integer :: nblok ! Number of blocks in the ddb
161 : logical :: has_d0E_xx=.false.
162 : logical :: has_d1E_xx=.false.
163 : logical :: has_d2E_xx=.false.
164 : logical :: has_d3E_xx=.false.
165 : logical :: has_d3E_lw=.false. ! Note that d3E_xx is implied by d3E_lw
166 : logical :: has_d2eig=.false.
167 :
168 : real(dp) :: dilatmx
169 : real(dp) :: ecut
170 : real(dp) :: ecutsm
171 : real(dp) :: kptnrm
172 : real(dp) :: pawecutdg
173 : real(dp) :: dfpt_sciss
174 : real(dp) :: tolwfr
175 : real(dp) :: tphysel
176 : real(dp) :: tsmear
177 :
178 : character(len=fnlen) :: dscrpt
179 :
180 : integer :: ngfft(18)
181 : real(dp) :: acell(3)
182 : real(dp) :: rprim(3,3)
183 :
184 : integer,allocatable :: nband(:)
185 : ! nband(mkpt*nsppol)
186 :
187 : integer,allocatable :: symafm(:)
188 : ! symafm(msym)
189 :
190 : integer,allocatable :: symrel(:,:,:)
191 : ! symrel(3,3,msym)
192 :
193 : integer,allocatable :: typat(:)
194 : ! typat(matom)
195 :
196 : integer,allocatable :: typ(:)
197 : ! typ(nblok), type of the blocks: d2E_st, d3E_lw,...
198 :
199 : real(dp),allocatable :: amu(:)
200 : ! amu(mtypat)
201 :
202 : real(dp),allocatable :: kpt(:,:)
203 : ! kpt(3,mkpt)
204 :
205 : real(dp),allocatable :: occ(:)
206 : ! occ(mband*mkpt*nsppol)
207 :
208 : real(dp),allocatable :: spinat(:,:)
209 : ! spinat(3,matom)
210 :
211 : real(dp),allocatable :: tnons(:,:)
212 : ! tnons(3,msym)
213 :
214 : real(dp),allocatable :: wtk(:)
215 : ! wtk(mkpt)
216 :
217 : real(dp),allocatable :: xred(:,:)
218 : ! xred(3,matom)
219 :
220 : real(dp),allocatable :: zion(:)
221 : ! zion(mtypat)
222 :
223 : real(dp),allocatable :: znucl(:)
224 : ! znucl(mtypat)
225 :
226 : type(pawtab_type),allocatable :: pawtab(:)
227 : ! pawtab(psps%ntypat*psps%usepaw)
228 :
229 : type(pseudopotential_type) :: psps
230 :
231 : type(crystal_t) :: crystal
232 :
233 : contains
234 :
235 : procedure :: init => ddb_hdr_init
236 : ! Initialize object.
237 :
238 : procedure :: init_from_crystal => ddb_hdr_init_from_crystal
239 : ! Initialize object from a crystal.
240 :
241 : procedure :: malloc => ddb_hdr_malloc
242 : ! Allocate dynamic memory.
243 :
244 : procedure :: free => ddb_hdr_free
245 : ! Free dynamic memory.
246 :
247 : procedure :: bcast_dim => ddb_hdr_bcast_dim
248 : ! Master broadcasts header dimensions.
249 :
250 : procedure :: bcast => ddb_hdr_bcast
251 : ! Master broadcasts header data.
252 :
253 : procedure :: print => ddb_hdr_print
254 : ! Print out the content of the header.
255 :
256 : procedure :: compare => ddb_hdr_compare
257 : ! Compare two DDB headers.
258 :
259 : procedure :: copy_missing_variables => ddb_hdr_copy_missing_variables
260 : ! Copy some missing variables from one header to another
261 :
262 : procedure :: copy_psps_from => ddb_hdr_copy_psps_from
263 : ! Copy the information on pseudopotentials from one header to another
264 :
265 : procedure :: set_typ => ddb_hdr_set_typ
266 : ! Set the type of each blok in the ddb
267 :
268 : procedure :: get_block_dims => ddb_hdr_get_block_dims
269 : ! Compute mpert and msize
270 :
271 : procedure :: open_write => ddb_hdr_open_write
272 : ! Open the DDB file and write the header.
273 :
274 : procedure :: open_write_txt => ddb_hdr_open_write_txt
275 : ! Open the DDB text file and write the header.
276 :
277 : procedure :: open_write_nc => ddb_hdr_open_write_nc
278 : ! Open the DDB NetCDF file and write the header.
279 :
280 : procedure :: open_read => ddb_hdr_open_read
281 : ! Open the DDB file and read the header.
282 :
283 : procedure :: open_read_txt => ddb_hdr_open_read_txt
284 : ! Open the DDB text file and read the header.
285 :
286 : procedure :: open_read_nc => ddb_hdr_open_read_nc
287 : ! Open the DDB NetCDF file and read the header.
288 :
289 : procedure :: close => ddb_hdr_close
290 : ! Close any open file
291 :
292 : procedure :: get_iomode => ddb_hdr_get_iomode
293 : ! Decide on iomode based on file name and current iomode parameter
294 :
295 : end type ddb_hdr_type
296 :
297 : CONTAINS !===========================================================
298 : !!***
299 :
300 : !----------------------------------------------------------------------
301 :
302 : !!****f* m_ddb_hdr/ddb_hdr_init_from_crystal
303 : !! NAME
304 : !! ddb_hdr_init_from_crystal
305 : !!
306 : !! FUNCTION
307 : !! Initialize a ddb_hdr object from a crystal.
308 : !! This creates a minimal ddb_hdr with no information whatsoever
309 : !! about the parameters that created the ddb.
310 : !!
311 : !! INPUTS
312 : !! ddb_hdr=the new ddb_hdr object
313 : !! crystal=crystal object
314 : !!
315 : !! SOURCE
316 :
317 44 : subroutine ddb_hdr_init_from_crystal(ddb_hdr, crystal)
318 :
319 : !Arguments ------------------------------------
320 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
321 : type(crystal_t),intent(inout) :: crystal
322 : ! ************************************************************************
323 :
324 44 : ddb_hdr%nblok = 1
325 44 : ddb_hdr%ddb_version = DDB_VERSION
326 44 : ddb_hdr%dscrpt = ''
327 :
328 44 : ddb_hdr%with_psps = 0
329 44 : ddb_hdr%with_dfpt_vars = 0
330 :
331 44 : ddb_hdr%matom = crystal%natom
332 44 : ddb_hdr%natom = crystal%natom
333 44 : ddb_hdr%msym = crystal%nsym
334 44 : ddb_hdr%nsym = crystal%nsym
335 44 : ddb_hdr%mtypat = crystal%ntypat
336 44 : ddb_hdr%ntypat = crystal%ntypat
337 :
338 44 : ddb_hdr%mband = 1
339 44 : ddb_hdr%mkpt = 1
340 44 : ddb_hdr%nkpt = 1
341 44 : ddb_hdr%mpert = crystal%natom
342 44 : ddb_hdr%msize = 3*ddb_hdr%mpert*3*ddb_hdr%mpert
343 :
344 44 : ddb_hdr%intxc = 0
345 44 : ddb_hdr%iscf = 0
346 44 : ddb_hdr%ixc = 0
347 44 : ddb_hdr%nspden = 1
348 44 : ddb_hdr%nspinor = 1
349 44 : ddb_hdr%nsppol = 1
350 44 : ddb_hdr%occopt = 0
351 44 : ddb_hdr%usepaw = 0
352 :
353 44 : ddb_hdr%dilatmx = 0
354 44 : ddb_hdr%ecut = 0
355 44 : ddb_hdr%ecutsm = 0
356 44 : ddb_hdr%kptnrm = 1
357 44 : ddb_hdr%pawecutdg = 0
358 44 : ddb_hdr%dfpt_sciss = 0
359 44 : ddb_hdr%tolwfr = 0
360 44 : ddb_hdr%tphysel = 0
361 44 : ddb_hdr%tsmear = 0
362 :
363 836 : ddb_hdr%ngfft = zero
364 176 : ddb_hdr%acell = one
365 :
366 44 : call ddb_hdr%malloc()
367 :
368 572 : ddb_hdr%rprim(:,:) = crystal%rprimd(:,:)
369 1666 : ddb_hdr%symafm(:) = crystal%symafm(:)
370 21130 : ddb_hdr%symrel(:,:,:) = crystal%symrel(:,:,:)
371 6532 : ddb_hdr%tnons(:,:) = crystal%tnons(:,:)
372 150 : ddb_hdr%typat(:) = crystal%typat(:)
373 99 : ddb_hdr%amu(:) = crystal%amu(:)
374 468 : ddb_hdr%xred(:,:) = crystal%xred(:,:)
375 99 : ddb_hdr%zion(:) = crystal%zion(:)
376 99 : ddb_hdr%znucl(:) = crystal%znucl(:)
377 :
378 44 : call crystal%copy(ddb_hdr%crystal)
379 :
380 : ! Pseudopotential info will not be written by this header,
381 : ! but we still need to initialize certain dimensions and arrays
382 : ! because they will be passed to psddb8
383 44 : ddb_hdr%psps%dimekb = zero
384 44 : ddb_hdr%psps%lmnmax = zero
385 44 : ddb_hdr%psps%useylm = zero
386 44 : ddb_hdr%psps%usepaw = zero
387 88 : ABI_MALLOC(ddb_hdr%psps%ekb,(ddb_hdr%psps%dimekb,ddb_hdr%ntypat))
388 176 : ABI_MALLOC(ddb_hdr%psps%indlmn,(6,ddb_hdr%psps%lmnmax,ddb_hdr%ntypat))
389 132 : ABI_MALLOC(ddb_hdr%psps%pspso,(ddb_hdr%ntypat))
390 132 : ABI_MALLOC(ddb_hdr%pawtab,(ddb_hdr%ntypat*ddb_hdr%psps%usepaw))
391 44 : call pawtab_nullify(ddb_hdr%pawtab)
392 99 : ddb_hdr%psps%ekb = zero
393 99 : ddb_hdr%psps%indlmn = zero
394 99 : ddb_hdr%psps%pspso = zero
395 :
396 44 : end subroutine ddb_hdr_init_from_crystal
397 : !!***
398 :
399 : !----------------------------------------------------------------------
400 :
401 : !!****f* m_ddb_hdr/ddb_hdr_init
402 : !! NAME
403 : !! ddb_hdr_init
404 : !!
405 : !! FUNCTION
406 : !! Initialize a ddb_hdr object from a dataset.
407 : !!
408 : !! INPUTS
409 : !! ddb_hdr=the new ddb_hdr object
410 : !! dtset=dtset object of the current calculation
411 : !! psps=the pseudopotential objects
412 : !! pawtab=pawtab object
413 : !! dscrpt=string description of the ddb
414 : !! nblok=number of blocks
415 : !! xred=reduced coordinate positions of the atoms
416 : !! occ=occupation number of the bands
417 : !! ngfft=fft grid
418 : !! mband=maximum number of bands
419 : !! mkpt=maximum number of kpoints
420 : !!
421 : !! SOURCE
422 :
423 6798 : subroutine ddb_hdr_init(ddb_hdr, dtset, psps, pawtab, dscrpt, &
424 : nblok, msize, mpert, xred, occ, ngfft,&
425 137 : mband,nkpt,kpt)
426 :
427 : !Arguments ------------------------------------
428 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
429 : type(dataset_type),intent(in) :: dtset
430 : type(pseudopotential_type),intent(inout) :: psps
431 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
432 : character(len=*),intent(in) :: dscrpt ! TODO: Make this one optional.
433 : integer,intent(in) :: nblok
434 : !integer,intent(in),optional :: mpert ! GA: TODO
435 : integer,intent(in),optional :: mband,nkpt,msize,mpert
436 : integer,intent(in),optional :: ngfft(18)
437 : real(dp),intent(in),optional :: xred(3,dtset%natom)
438 : real(dp),intent(in),optional :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
439 : real(dp),intent(in),optional :: kpt(:,:)
440 :
441 : !Local variables -------------------------
442 : integer :: ii, nn, ikpt
443 : ! ************************************************************************
444 :
445 3399 : ddb_hdr%nblok = nblok
446 3399 : ddb_hdr%ddb_version = DDB_VERSION
447 3399 : ddb_hdr%dscrpt = trim(dscrpt)
448 :
449 3399 : call psps_copy(psps, ddb_hdr%psps)
450 :
451 : ! Copy scalars from dtset
452 3399 : ddb_hdr%matom = dtset%natom
453 3399 : ddb_hdr%natom = dtset%natom
454 3399 : ddb_hdr%mband = dtset%mband
455 3399 : ddb_hdr%mkpt = dtset%nkpt
456 3399 : ddb_hdr%nkpt = dtset%nkpt
457 3399 : ddb_hdr%msym = dtset%nsym
458 3399 : ddb_hdr%nsym = dtset%nsym
459 3399 : ddb_hdr%mtypat = dtset%ntypat
460 3399 : ddb_hdr%ntypat = dtset%ntypat
461 3399 : ddb_hdr%nspden = dtset%nspden
462 3399 : ddb_hdr%nspinor = dtset%nspinor
463 3399 : ddb_hdr%nsppol = dtset%nsppol
464 :
465 3399 : if (present(mband)) then
466 137 : ddb_hdr%mband = mband
467 : end if
468 :
469 3399 : if (present(nkpt)) then
470 137 : ddb_hdr%mkpt = nkpt
471 137 : ddb_hdr%nkpt = nkpt
472 : end if
473 :
474 3399 : if (present(mpert)) then
475 137 : ddb_hdr%mpert = mpert
476 : end if
477 :
478 3399 : if (present(msize)) then
479 0 : ddb_hdr%msize = msize
480 : else
481 3399 : if (present(mpert)) then
482 137 : call ddb_hdr%get_block_dims()
483 : end if
484 : end if
485 :
486 3399 : ddb_hdr%occopt = dtset%occopt
487 3399 : ddb_hdr%usepaw = dtset%usepaw
488 :
489 3399 : ddb_hdr%intxc = dtset%intxc
490 3399 : ddb_hdr%ixc = dtset%ixc
491 3399 : ddb_hdr%iscf = dtset%iscf
492 :
493 3399 : ddb_hdr%dilatmx = dtset%dilatmx
494 3399 : ddb_hdr%ecut = dtset%ecut
495 3399 : ddb_hdr%ecutsm = dtset%ecutsm
496 3399 : ddb_hdr%pawecutdg = dtset%pawecutdg
497 3399 : ddb_hdr%kptnrm = dtset%kptnrm
498 3399 : ddb_hdr%dfpt_sciss = dtset%dfpt_sciss
499 3399 : ddb_hdr%tolwfr = one
500 3399 : ddb_hdr%tphysel = dtset%tphysel
501 3399 : ddb_hdr%tsmear = dtset%tsmear
502 :
503 3399 : ddb_hdr%with_psps = 1
504 3399 : ddb_hdr%with_dfpt_vars = 1
505 :
506 3399 : ddb_hdr%iomode = dtset%iomode
507 3399 : ddb_hdr%has_open_file_txt = .false.
508 3399 : ddb_hdr%has_open_file_nc = .false.
509 :
510 3399 : call ddb_hdr%malloc()
511 :
512 3399 : if (present(ngfft)) then
513 61199 : ddb_hdr%ngfft = ngfft
514 : else
515 3382 : ddb_hdr%ngfft = dtset%ngfft
516 : end if
517 :
518 : ! Copy arrays from dtset
519 13596 : ddb_hdr%acell(:) = dtset%acell_orig(1:3,1)
520 44187 : ddb_hdr%rprim(:,:) = dtset%rprim_orig(1:3,1:3,1)
521 7689 : ddb_hdr%amu(:) = dtset%amu_orig(:,1)
522 63428 : ddb_hdr%nband(:) = dtset%nband(1:ddb_hdr%mkpt*ddb_hdr%nsppol)
523 83509 : ddb_hdr%symafm(:) = dtset%symafm(1:ddb_hdr%msym)
524 1044829 : ddb_hdr%symrel(:,:,:) = dtset%symrel(1:3,1:3,1:ddb_hdr%msym)
525 11191 : ddb_hdr%typat(:) = dtset%typat(1:ddb_hdr%matom)
526 61554 : ddb_hdr%wtk(:) = dtset%wtk(1:ddb_hdr%mkpt)
527 34567 : ddb_hdr%spinat(:,:) = dtset%spinat(1:3,1:ddb_hdr%matom)
528 323839 : ddb_hdr%tnons(:,:) = dtset%tnons(1:3,1:ddb_hdr%msym)
529 7689 : ddb_hdr%zion(:) = dtset%ziontypat(1:ddb_hdr%mtypat)
530 7689 : ddb_hdr%znucl(:) = dtset%znucl(1:ddb_hdr%mtypat)
531 :
532 3399 : if (present(kpt)) then
533 2200 : do ikpt=1,ddb_hdr%nkpt
534 8389 : do ii = 1,3
535 8252 : ddb_hdr%kpt(ii,ikpt) = kpt(ii,ikpt)
536 : end do
537 : end do
538 : else
539 227630 : ddb_hdr%kpt(:,:) = dtset%kpt(1:3,1:ddb_hdr%mkpt)
540 : end if
541 :
542 3399 : if (present(mband)) then
543 2200 : do ii=1,ddb_hdr%nkpt
544 2200 : ddb_hdr%nband(ii) = mband
545 : end do
546 : end if
547 :
548 3399 : if (present(xred)) then
549 34567 : ddb_hdr%xred(:,:) = xred(1:3,1:ddb_hdr%matom)
550 : else
551 0 : ddb_hdr%xred(:,:) = dtset%xred_orig(1:3,1:ddb_hdr%matom,1)
552 : end if
553 3399 : if (present(occ)) then
554 420823 : ddb_hdr%occ(:) = occ(1:ddb_hdr%mband*ddb_hdr%mkpt*ddb_hdr%nsppol)
555 : else
556 0 : ddb_hdr%occ(:) = dtset%occ_orig(1:ddb_hdr%mband*ddb_hdr%mkpt*ddb_hdr%nsppol,1)
557 : end if
558 :
559 : ! GA: I had way too many problems implementing pawtab_copy.
560 : ! The script check-libpaw would report all sorts of errors.
561 : ! Therefore, I do a cheap copy here, copying only the relevant info.
562 : !call pawtab_copy(pawtab, ddb_hdr%pawtab)
563 3399 : nn=size(pawtab)
564 11192 : ABI_MALLOC(ddb_hdr%pawtab,(ddb_hdr%psps%ntypat*ddb_hdr%psps%usepaw))
565 3399 : if (nn.gt.0) then
566 :
567 653 : call pawtab_nullify(ddb_hdr%pawtab)
568 :
569 1648 : do ii=1,nn
570 995 : ddb_hdr%pawtab(ii)%basis_size = pawtab(ii)%basis_size
571 995 : ddb_hdr%pawtab(ii)%lmn_size = pawtab(ii)%lmn_size
572 995 : ddb_hdr%pawtab(ii)%lmn2_size = pawtab(ii)%lmn2_size
573 995 : ddb_hdr%pawtab(ii)%rpaw = pawtab(ii)%rpaw
574 995 : ddb_hdr%pawtab(ii)%rshp = pawtab(ii)%rshp
575 995 : ddb_hdr%pawtab(ii)%shape_type = pawtab(ii)%shape_type
576 1648 : if (allocated(pawtab(ii)%dij0)) then
577 995 : call alloc_copy(pawtab(ii)%dij0, ddb_hdr%pawtab(ii)%dij0)
578 : end if
579 : end do
580 : end if
581 :
582 : call ddb_hdr%crystal%init(dtset%amu_orig(:,1), &
583 : & dtset%spgroup, dtset%natom, dtset%npsp, psps%ntypat, &
584 : & dtset%nsym, dtset%rprimd_orig(:,:,1), dtset%typat, &
585 : & ddb_hdr%xred, dtset%ziontypat, dtset%znucl, 1, &
586 : & dtset%nspden==2.and.dtset%nsppol==1, .false., '', &
587 6769 : & dtset%symrel, dtset%tnons, dtset%symafm)
588 :
589 3399 : end subroutine ddb_hdr_init
590 : !!***
591 :
592 : !----------------------------------------------------------------------
593 :
594 : !!****f* m_ddb_hdr/ddb_hdr_malloc
595 : !! NAME
596 : !! ddb_hdr_malloc
597 : !!
598 : !! FUNCTION
599 : !! Allocate dynamic memory.
600 : !!
601 : !! SOURCE
602 :
603 4242 : subroutine ddb_hdr_malloc(ddb_hdr)
604 :
605 : !Arguments ------------------------------------
606 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
607 : ! ************************************************************************
608 :
609 : ! integer
610 12726 : ABI_MALLOC(ddb_hdr%nband,(ddb_hdr%mkpt*ddb_hdr%nsppol))
611 12726 : ABI_MALLOC(ddb_hdr%symafm,(ddb_hdr%msym))
612 12726 : ABI_MALLOC(ddb_hdr%symrel,(3,3,ddb_hdr%msym))
613 12726 : ABI_MALLOC(ddb_hdr%typat,(ddb_hdr%matom))
614 :
615 : ! real
616 12726 : ABI_MALLOC(ddb_hdr%amu,(ddb_hdr%mtypat))
617 12726 : ABI_MALLOC(ddb_hdr%kpt,(3, ddb_hdr%mkpt))
618 12726 : ABI_MALLOC(ddb_hdr%occ,(ddb_hdr%mband*ddb_hdr%mkpt*ddb_hdr%nsppol))
619 12726 : ABI_MALLOC(ddb_hdr%spinat,(3, ddb_hdr%matom))
620 12726 : ABI_MALLOC(ddb_hdr%tnons,(3, ddb_hdr%msym))
621 12726 : ABI_MALLOC(ddb_hdr%wtk,(ddb_hdr%mkpt))
622 8484 : ABI_MALLOC(ddb_hdr%xred,(3, ddb_hdr%matom))
623 8484 : ABI_MALLOC(ddb_hdr%zion,(ddb_hdr%mtypat))
624 8484 : ABI_MALLOC(ddb_hdr%znucl,(ddb_hdr%mtypat))
625 :
626 : ! types
627 : !ABI_MALLOC(ddb_hdr%pawtab,(ddb_hdr%psps%ntypat*ddb_hdr%psps%usepaw))
628 : !call pawtab_nullify(ddb_hdr%pawtab)
629 :
630 : ! Do not allocate the pseudo
631 :
632 80598 : ddb_hdr%ngfft = zero
633 264955 : ddb_hdr%nband = zero
634 210175 : ddb_hdr%symafm = one
635 2681371 : ddb_hdr%symrel = zero
636 14109 : ddb_hdr%typat = zero
637 9905 : ddb_hdr%amu = zero
638 809582 : ddb_hdr%kpt = zero
639 4953449 : ddb_hdr%occ = zero
640 43710 : ddb_hdr%spinat = zero
641 827974 : ddb_hdr%tnons = zero
642 205577 : ddb_hdr%wtk = zero
643 43710 : ddb_hdr%xred = zero
644 9905 : ddb_hdr%zion = zero
645 9905 : ddb_hdr%znucl = zero
646 :
647 4242 : end subroutine ddb_hdr_malloc
648 : !!***
649 :
650 : !----------------------------------------------------------------------
651 :
652 : !!****f* m_ddb_hdr/ddb_hdr_set_typ
653 : !! NAME
654 : !! ddb_hdr_set_typ
655 : !!
656 : !! FUNCTION
657 : !! Set the type of each block.
658 : !!
659 : !! INPUTS
660 : !! nblok=number of blocks
661 : !! typ=type of each block
662 : !!
663 : !! SOURCE
664 :
665 3346 : subroutine ddb_hdr_set_typ(ddb_hdr, nblok, typ)
666 :
667 : !Arguments ------------------------------------
668 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
669 : integer,intent(in) :: nblok
670 : integer,intent(in) :: typ(nblok)
671 : ! ************************************************************************
672 :
673 3346 : ddb_hdr%nblok = nblok
674 3346 : ABI_SFREE(ddb_hdr%typ)
675 10038 : ABI_MALLOC(ddb_hdr%typ,(nblok))
676 9575 : ddb_hdr%typ(:) = typ(:)
677 :
678 3346 : end subroutine ddb_hdr_set_typ
679 : !!***
680 :
681 : !----------------------------------------------------------------------
682 :
683 : !!****f* m_ddb_hdr/ddb_hdr_get_block_dims
684 : !! NAME
685 : !! ddb_hdr_get_block_dims
686 : !!
687 : !! FUNCTION
688 : !! Get the maximum number of perturbations from natom,
689 : !! and get the maximum matrix size of one block, knowing the block types.
690 : !!
691 : !! SOURCE
692 :
693 1310 : subroutine ddb_hdr_get_block_dims(ddb_hdr)
694 :
695 : !Arguments ------------------------------------
696 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
697 : ! ************************************************************************
698 :
699 : ! Compute mpert
700 : !ddb_hdr%mpert = 2*ddb_hdr%natom+MPERT_MAX
701 : ! GA: mpert is stored in netcdf format but not in text format.
702 :
703 : ! Compute msize
704 1310 : if (ddb_hdr%has_d3E_xx) then
705 69 : ddb_hdr%msize=3*ddb_hdr%mpert*3*ddb_hdr%mpert*3*ddb_hdr%mpert
706 : else
707 1241 : ddb_hdr%msize=3*ddb_hdr%mpert*3*ddb_hdr%mpert
708 : end if
709 :
710 1310 : end subroutine ddb_hdr_get_block_dims
711 : !!***
712 :
713 : !----------------------------------------------------------------------
714 :
715 : !!****f* m_ddb_hdr/ddb_hdr_free
716 : !! NAME
717 : !! ddb_hdr_free
718 : !!
719 : !! FUNCTION
720 : !! Free memory.
721 : !!
722 : !! SOURCE
723 :
724 4999 : subroutine ddb_hdr_free(ddb_hdr)
725 :
726 : !Arguments ------------------------------------
727 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
728 : ! ************************************************************************
729 :
730 : ! integer
731 4999 : ABI_SFREE(ddb_hdr%nband)
732 4999 : ABI_SFREE(ddb_hdr%symafm)
733 4999 : ABI_SFREE(ddb_hdr%symrel)
734 4999 : ABI_SFREE(ddb_hdr%typat)
735 :
736 : ! real
737 4999 : ABI_SFREE(ddb_hdr%amu)
738 4999 : ABI_SFREE(ddb_hdr%kpt)
739 4999 : ABI_SFREE(ddb_hdr%occ)
740 4999 : ABI_SFREE(ddb_hdr%spinat)
741 4999 : ABI_SFREE(ddb_hdr%tnons)
742 4999 : ABI_SFREE(ddb_hdr%wtk)
743 4999 : ABI_SFREE(ddb_hdr%xred)
744 4999 : ABI_SFREE(ddb_hdr%zion)
745 4999 : ABI_SFREE(ddb_hdr%znucl)
746 4999 : ABI_SFREE(ddb_hdr%typ)
747 :
748 : ! types
749 4999 : call ddb_hdr%crystal%free()
750 4999 : call psps_free(ddb_hdr%psps)
751 :
752 4999 : if (allocated(ddb_hdr%pawtab)) then
753 4408 : call pawtab_free(ddb_hdr%pawtab)
754 5444 : ABI_FREE(ddb_hdr%pawtab)
755 : end if
756 :
757 4999 : end subroutine ddb_hdr_free
758 : !!***
759 :
760 : !----------------------------------------------------------------------
761 :
762 : !!****f* m_ddb_hdr/ddb_hdr_open_write_txt
763 : !! NAME
764 : !! ddb_hdr_open_write_txt
765 : !!
766 : !! FUNCTION
767 : !! Open the DDB text file and write the header.
768 : !!
769 : !! INPUTS
770 : !! filename=name of the file being written (abo_DS*_DDB)
771 : !! with_psps
772 : !! 1-> include information on pseudopoentials
773 : !! 0-> do not include information on pseudopoentials
774 : !!
775 : !! OUTPUT
776 : !! unddb=unit of the open ddb file in text format.
777 : !!
778 : !! SOURCE
779 :
780 150 : subroutine ddb_hdr_open_write_txt(ddb_hdr, filename, with_psps)
781 :
782 : !Arguments ------------------------------------
783 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
784 : character(len=*),intent(in) :: filename
785 : integer,intent(in),optional :: with_psps
786 :
787 : !Local variables -------------------------
788 : character(len=500) :: message
789 : integer :: unddb
790 : integer :: ierr
791 : ! ************************************************************************
792 :
793 150 : if (present(with_psps)) then
794 48 : ddb_hdr%with_psps = with_psps
795 48 : ddb_hdr%with_dfpt_vars = with_psps
796 : end if
797 :
798 150 : unddb = get_unit()
799 150 : ddb_hdr%unddb = unddb
800 150 : ddb_hdr%has_open_file_txt = .true.
801 :
802 : ! Open the output derivative database.
803 150 : ierr = open_file(filename,message,unit=unddb,status='unknown',form='formatted')
804 150 : if (ierr /= 0) then
805 0 : ABI_ERROR(message)
806 : end if
807 :
808 150 : call ddb_hdr%print(unddb)
809 :
810 150 : end subroutine ddb_hdr_open_write_txt
811 : !!***
812 :
813 : !----------------------------------------------------------------------
814 :
815 : !!****f* m_ddb_hdr/ddb_hdr_open_write_nc
816 : !! NAME
817 : !! ddb_hdr_open_write_nc
818 : !!
819 : !! FUNCTION
820 : !! Open the DDB NetCDF file and write the header.
821 : !! Create the groups and define all variables.
822 : !!
823 : !! INPUTS
824 : !! filename=name of the file being written
825 : !! with_psps
826 : !! 1-> include information on pseudopoentials
827 : !! 0-> do not include information on pseudopoentials
828 : !!
829 : !! OUTPUT
830 : !! ncid=netcdf identifier.
831 : !!
832 : !! SOURCE
833 :
834 3279 : subroutine ddb_hdr_open_write_nc(ddb_hdr, filename, with_psps, with_dfpt_vars)
835 :
836 : !Arguments ------------------------------------
837 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
838 : character(len=*),intent(in) :: filename
839 : integer,intent(in),optional :: with_psps
840 : integer,intent(in),optional :: with_dfpt_vars
841 :
842 : !Local variables -------------------------
843 : !scalars
844 : integer :: ncid
845 : integer :: ncid_crystal, ncid_psps, ncid_pawtab, ncid_d0E, ncid_d1E, ncid_d2E, ncid_d3E, ncid_d2eig
846 : integer :: ncerr,nkpt,bandtot,occopt,blktyp
847 : integer :: ii,iat,itypat,isppol,ikpt,iband,iblok,giblok
848 : integer :: nblok_d0E,nblok_d1E,nblok_d2E,nblok_d3E,nblok_d2eig
849 : integer :: iblok_d0E,iblok_d1E,iblok_d2E,iblok_d3E,iblok_d2eig
850 : integer :: basis_size, max_basis_size, lmn_size, lmn2_size, max_lmn2_size, shape_type
851 : integer :: natom, mpert
852 : !arrays
853 3279 : integer,allocatable :: available_block_types(:)
854 3279 : integer,allocatable :: blktyp_d0E(:),blktyp_d1E(:),blktyp_d2E(:),blktyp_d3E(:),blktyp_d2eig(:)
855 3279 : integer,allocatable :: nband(:,:)
856 : character(len=descrlen) :: descr
857 3279 : real(dp),allocatable :: occ(:,:,:)
858 3279 : integer,allocatable :: pawtab_basis_size(:)
859 3279 : integer,allocatable :: pawtab_lmn_size(:)
860 3279 : integer,allocatable :: pawtab_lmn2_size(:)
861 3279 : integer,allocatable :: pawtab_shape_type(:)
862 3279 : real(dp),allocatable :: pawtab_rshp(:)
863 3279 : real(dp),allocatable :: pawtab_rpaw(:)
864 :
865 : ! ************************************************************************
866 :
867 3279 : if (present(with_psps)) then
868 2835 : ddb_hdr%with_psps = with_psps
869 : end if
870 :
871 3279 : if (present(with_dfpt_vars)) then
872 0 : ddb_hdr%with_dfpt_vars = with_dfpt_vars
873 : else
874 3279 : ddb_hdr%with_dfpt_vars = ddb_hdr%with_psps
875 : end if
876 :
877 : ! Initialize NetCDF file.
878 3279 : NCF_CHECK(nctk_open_create(ncid, filename, xmpi_comm_self))
879 3279 : ddb_hdr%ncid = ncid
880 3279 : ddb_hdr%has_open_file_nc = .True.
881 :
882 : ! Write version number
883 6558 : NCF_CHECK(nctk_defnwrite_ivars(ncid, ['ddb_version'], [DDB_VERSION_NC]))
884 :
885 : ! ----------
886 : ! Dimensions
887 : ! ----------
888 :
889 : ! Dimensions specified by ETSF
890 : ncerr = nctk_def_dims(ncid, [&
891 : nctkdim_t('number_of_cartesian_directions', 3), &
892 : nctkdim_t('number_of_atoms', ddb_hdr%natom), &
893 : nctkdim_t('number_of_kpoints', ddb_hdr%nkpt), &
894 : nctkdim_t('number_of_spins', ddb_hdr%nsppol), &
895 : nctkdim_t('number_of_spinor_components', ddb_hdr%nspinor), &
896 : nctkdim_t('number_of_spin_densities', ddb_hdr%nspden), &
897 : nctkdim_t('one_dim', 1), &
898 : nctkdim_t('two_dim', 2), &
899 : nctkdim_t('three_dim', 3), &
900 36069 : nctkdim_t('cplex',2)], defmode=.True.)
901 3279 : NCF_CHECK(ncerr)
902 :
903 : ! Dimensions relevant to the DDB
904 : ncerr = nctk_def_dims(ncid, [&
905 : nctkdim_t('number_of_symmetries', ddb_hdr%nsym), &
906 : nctkdim_t('number_of_atom_types', ddb_hdr%ntypat), &
907 : nctkdim_t('maximum_number_of_bands', ddb_hdr%mband) &
908 13116 : ], defmode=.True.)
909 3279 : NCF_CHECK(ncerr)
910 :
911 : ! Compute total number of bands
912 3279 : occopt = ddb_hdr%occopt
913 3279 : nkpt = ddb_hdr%nkpt
914 3279 : if(ddb_hdr%occopt==2)then
915 : bandtot=0
916 444 : do ikpt=1,nkpt
917 444 : bandtot=bandtot+ddb_hdr%nband(ikpt)*ddb_hdr%nsppol
918 : end do
919 : else
920 3076 : bandtot=nkpt*ddb_hdr%nband(1)*ddb_hdr%nsppol
921 : end if
922 : ncerr = nctk_def_dims(ncid, [&
923 : nctkdim_t('total_number_of_states', bandtot) &
924 6558 : ], defmode=.True.)
925 3279 : NCF_CHECK(ncerr)
926 :
927 : ! Construct local array of occupations and nband
928 16395 : ABI_MALLOC(occ, (ddb_hdr%mband,ddb_hdr%nkpt,ddb_hdr%nsppol))
929 479832 : occ = zero
930 : ii = 0
931 6989 : do isppol=1,ddb_hdr%nsppol
932 65819 : do ikpt=1,ddb_hdr%nkpt
933 476553 : do iband=1,ddb_hdr%mband
934 414013 : ii = ii + 1
935 472843 : occ(iband,ikpt,isppol) = ddb_hdr%occ(ii)
936 : end do
937 : end do
938 : end do
939 13116 : ABI_MALLOC(nband, (ddb_hdr%nkpt,ddb_hdr%nsppol))
940 65819 : nband = zero
941 : ii = 0
942 6989 : do isppol=1,ddb_hdr%nsppol
943 65819 : do ikpt=1,ddb_hdr%nkpt
944 58830 : ii = ii + 1
945 62540 : nband(ikpt,isppol) = ddb_hdr%nband(ii)
946 : end do
947 : end do
948 :
949 :
950 9837 : ncerr = nctk_defnwrite_ivars(ncid, ["with_psps"], [ddb_hdr%with_psps])
951 3279 : NCF_CHECK(ncerr)
952 9837 : ncerr = nctk_defnwrite_ivars(ncid, ["with_dfpt_vars"], [ddb_hdr%with_dfpt_vars])
953 : !ncerr = nctk_defnwrite_ivars(ncid, &
954 : !& [ character(len=nctk_slen) :: &
955 : !& "with_psps", "with_dfpt_vars"], &
956 : !& [ddb_hdr%with_psps, ddb_hdr%with_dfpt_vars])
957 3279 : NCF_CHECK(ncerr)
958 :
959 : ! --------------------------
960 : ! Write objects in subgroups
961 : ! --------------------------
962 :
963 : ! Crystal info
964 3279 : NCF_CHECK(nf90_def_grp(ncid, 'Crystal', ncid_crystal))
965 3279 : NCF_CHECK(ddb_hdr%crystal%ncwrite(ncid_crystal))
966 :
967 : ! Pseudopotentials info
968 3279 : NCF_CHECK(nf90_def_grp(ncid, 'Pseudopotentials', ncid_psps))
969 3279 : call psps_ncwrite(ddb_hdr%psps, ncid_psps) ! We still need to write usepaw, for example.
970 :
971 : ! PAW info
972 : ! Note that PAW shares the namespace with Pseudopotentials.
973 3279 : if (ddb_hdr%psps%usepaw > 0) then
974 :
975 : ! -------------------------------------
976 : ! GA: TODO This should be a function pawtab_ncwrite
977 603 : NCF_CHECK(nf90_def_grp(ncid_psps, 'PAW', ncid_pawtab))
978 :
979 : ! get dimensions
980 1809 : ABI_MALLOC(pawtab_basis_size, (ddb_hdr%ntypat))
981 1206 : ABI_MALLOC(pawtab_lmn_size, (ddb_hdr%ntypat))
982 1206 : ABI_MALLOC(pawtab_lmn2_size, (ddb_hdr%ntypat))
983 1206 : ABI_MALLOC(pawtab_shape_type, (ddb_hdr%ntypat))
984 :
985 1536 : pawtab_basis_size = zero
986 1536 : pawtab_lmn_size = zero
987 1536 : pawtab_lmn2_size = zero
988 1536 : pawtab_shape_type = zero
989 :
990 : max_basis_size = 0 ; max_lmn2_size = 0
991 1536 : do itypat=1,ddb_hdr%ntypat
992 :
993 933 : basis_size = ddb_hdr%pawtab(itypat)%basis_size
994 933 : lmn_size = ddb_hdr%pawtab(itypat)%lmn_size
995 933 : lmn2_size = ddb_hdr%pawtab(itypat)%lmn2_size
996 933 : shape_type = ddb_hdr%pawtab(itypat)%shape_type
997 :
998 933 : pawtab_basis_size(itypat) = basis_size
999 933 : pawtab_lmn_size(itypat) = lmn_size
1000 933 : pawtab_lmn2_size(itypat) = lmn2_size
1001 933 : pawtab_shape_type(itypat) = shape_type
1002 :
1003 933 : if (lmn2_size > max_lmn2_size) then
1004 : max_lmn2_size = lmn2_size
1005 : end if
1006 1536 : if (basis_size > max_basis_size) then
1007 : max_basis_size = basis_size
1008 : end if
1009 : end do
1010 :
1011 : ncerr = nctk_def_dims(ncid_pawtab, [&
1012 : nctkdim_t('maximum_basis_size', max_basis_size), &
1013 : nctkdim_t('maximum_lmn2_size', max_lmn2_size)], &
1014 1809 : defmode=.True.)
1015 603 : NCF_CHECK(ncerr)
1016 :
1017 : ncerr = nctk_def_arrays(ncid_pawtab, [&
1018 : nctkarr_t('basis_size', "i", 'number_of_atom_types'), &
1019 : nctkarr_t('lmn_size', "i", 'number_of_atom_types'), &
1020 : nctkarr_t('lmn2_size', "i", 'number_of_atom_types'), &
1021 : nctkarr_t('shape_type', "i", 'number_of_atom_types') &
1022 3015 : &])
1023 603 : NCF_CHECK(ncerr)
1024 :
1025 : ncerr = nctk_def_arrays(ncid_pawtab, [&
1026 : nctkarr_t('Dij0', "dp", 'number_of_atom_types, maximum_lmn2_size'), &
1027 : nctkarr_t('r_PAW', "dp", 'number_of_atom_types'), &
1028 : nctkarr_t('r_shape', "dp", 'number_of_atom_types') &
1029 2412 : &])
1030 603 : NCF_CHECK(ncerr)
1031 :
1032 603 : NCF_CHECK(nf90_put_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'basis_size'), pawtab_basis_size))
1033 603 : NCF_CHECK(nf90_put_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'lmn_size'), pawtab_lmn_size))
1034 603 : NCF_CHECK(nf90_put_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'lmn2_size'), pawtab_lmn2_size))
1035 603 : NCF_CHECK(nf90_put_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'shape_type'), pawtab_shape_type))
1036 :
1037 603 : ABI_FREE(pawtab_basis_size)
1038 603 : ABI_FREE(pawtab_lmn_size)
1039 603 : ABI_FREE(pawtab_lmn2_size)
1040 603 : ABI_FREE(pawtab_shape_type)
1041 :
1042 1809 : ABI_MALLOC(pawtab_rpaw, (ddb_hdr%ntypat))
1043 1206 : ABI_MALLOC(pawtab_rshp, (ddb_hdr%ntypat))
1044 :
1045 1536 : do itypat=1,ddb_hdr%ntypat
1046 :
1047 933 : pawtab_rpaw(itypat) = ddb_hdr%pawtab(itypat)%rpaw
1048 933 : pawtab_rshp(itypat) = ddb_hdr%pawtab(itypat)%rshp
1049 :
1050 : ncerr = nf90_put_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'Dij0'), &
1051 : ddb_hdr%pawtab(itypat)%dij0, start=[itypat,1], &
1052 4665 : count=[1,ddb_hdr%pawtab(itypat)%lmn2_size])
1053 1536 : NCF_CHECK(ncerr)
1054 :
1055 : end do
1056 :
1057 603 : ncerr = nf90_put_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'r_PAW'), pawtab_rpaw)
1058 603 : NCF_CHECK(ncerr)
1059 :
1060 603 : ncerr = nf90_put_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'r_shape'), pawtab_rshp)
1061 603 : NCF_CHECK(ncerr)
1062 :
1063 603 : ABI_FREE(pawtab_rpaw)
1064 603 : ABI_FREE(pawtab_rshp)
1065 :
1066 : ! -------------------------------------
1067 :
1068 : end if
1069 :
1070 :
1071 : ! ----------------
1072 : ! Scalar variables
1073 : ! ----------------
1074 :
1075 : ! Scalar integer variables
1076 : ncerr = nctk_defnwrite_ivars(ncid, &
1077 : [character(len=nctk_slen) :: 'intxc', 'iscf', 'ixc', 'occopt'], &
1078 29511 : [ddb_hdr%intxc, ddb_hdr%iscf, ddb_hdr%ixc, ddb_hdr%occopt])
1079 3279 : NCF_CHECK(ncerr)
1080 :
1081 : ! Scalar float variables
1082 : ncerr = nctk_defnwrite_dpvars(ncid, &
1083 : [character(len=nctk_slen) :: 'dilatmx', 'ecut', 'pawecutdg', 'ecutsm', 'kptnrm'], &
1084 36069 : [ddb_hdr%dilatmx, ddb_hdr%ecut, ddb_hdr%pawecutdg, ddb_hdr%ecutsm, ddb_hdr%kptnrm])
1085 3279 : NCF_CHECK(ncerr)
1086 :
1087 : ncerr = nctk_defnwrite_dpvars(ncid, &
1088 : [character(len=nctk_slen) :: 'dfpt_sciss', 'tolwfr', 'tphysel', 'tsmear'], &
1089 29511 : [ddb_hdr%dfpt_sciss, ddb_hdr%tolwfr, ddb_hdr%tphysel, ddb_hdr%tsmear])
1090 3279 : NCF_CHECK(ncerr)
1091 :
1092 : ! ------
1093 : ! Arrays
1094 : ! ------
1095 :
1096 : ! FFT grid
1097 6558 : ncerr = nctk_def_arrays(ncid, [nctkarr_t('ngfft', "i", 'number_of_cartesian_directions')])
1098 3279 : NCF_CHECK(ncerr)
1099 :
1100 : ! Geometry
1101 : ncerr = nctk_def_arrays(ncid, [&
1102 : nctkarr_t('acell', "dp", 'number_of_cartesian_directions') &
1103 6558 : &])
1104 3279 : NCF_CHECK(ncerr)
1105 :
1106 : ! Basis
1107 : ! note: In case occopt!=2, only the first entry of number_of_bands matters
1108 : ncerr = nctk_def_arrays(ncid, [&
1109 : nctkarr_t('reduced_coordinates_of_kpoints', "dp", 'number_of_cartesian_directions, number_of_kpoints'), &
1110 : nctkarr_t('number_of_bands', "i", 'number_of_kpoints, number_of_spins'), &
1111 : nctkarr_t('spinat', "dp", 'number_of_cartesian_directions, number_of_atoms'), &
1112 : nctkarr_t('occupations', "dp", 'maximum_number_of_bands, number_of_kpoints, number_of_spins'), &
1113 : nctkarr_t('kpoints_weights', "dp", 'number_of_kpoints') &
1114 19674 : &])
1115 3279 : NCF_CHECK(ncerr)
1116 :
1117 : ! ----------------
1118 : ! Set array values
1119 : ! ----------------
1120 :
1121 : ! FFT grid
1122 3279 : NCF_CHECK(nf90_put_var(ncid, vid('ngfft'), ddb_hdr%ngfft(1:3)))
1123 :
1124 : ! Geometry
1125 3279 : NCF_CHECK(nf90_put_var(ncid, vid('acell'), ddb_hdr%acell(:)))
1126 :
1127 : ! Basis
1128 3279 : NCF_CHECK(nf90_put_var(ncid, vid('spinat'), ddb_hdr%spinat(:,:)))
1129 3279 : NCF_CHECK(nf90_put_var(ncid, vid('occupations'), occ))
1130 3279 : NCF_CHECK(nf90_put_var(ncid, vid('reduced_coordinates_of_kpoints'), ddb_hdr%kpt(:,:)))
1131 3279 : NCF_CHECK(nf90_put_var(ncid, vid('kpoints_weights'), ddb_hdr%wtk(:)))
1132 3279 : NCF_CHECK(nf90_put_var(ncid, vid('number_of_bands'), nband))
1133 :
1134 3279 : ABI_FREE(occ)
1135 3279 : ABI_FREE(nband)
1136 :
1137 : ! =============================================================
1138 : ! Write info on blocks, create groups, and define all variables
1139 : ! =============================================================
1140 :
1141 : ! --------------------------------------------
1142 : ! Description of perturbations and block types
1143 : ! --------------------------------------------
1144 3279 : mpert = ddb_hdr%mpert
1145 3279 : natom = ddb_hdr%natom
1146 :
1147 : ncerr = nctk_def_dims(ncid, [&
1148 : nctkdim_t('number_of_perturbations', ddb_hdr%mpert), &
1149 : nctkdim_t("number_of_available_block_types", ntypes), &
1150 : nctkdim_t("description_length", descrlen) &
1151 13116 : &], defmode=.True.)
1152 3279 : NCF_CHECK(ncerr)
1153 :
1154 : ncerr = nctk_def_arrays(ncid, [&
1155 : nctkarr_t('available_block_types', "i", 'number_of_available_block_types'), &
1156 : nctkarr_t('description_of_block_types', "char", 'description_length, number_of_available_block_types'), &
1157 : nctkarr_t('description_of_perturbations', "char", 'description_length, number_of_perturbations') &
1158 13116 : &])
1159 3279 : NCF_CHECK(ncerr)
1160 :
1161 :
1162 10799 : do iat=1,natom
1163 :
1164 7520 : write(descr, '(2a,i7)') trim(DESCR_ipert_0), ' ', iat
1165 :
1166 : ncerr = nf90_put_var(ncid, vid("description_of_perturbations"), &
1167 22560 : trim(descr) , start=[1, iat])
1168 10799 : NCF_CHECK(ncerr)
1169 : end do
1170 :
1171 :
1172 3279 : if (mpert >= natom+1) then
1173 : ncerr = nf90_put_var(ncid, vid("description_of_perturbations"), &
1174 9486 : trim(DESCR_ipert_1), start=[1, natom+1])
1175 3162 : NCF_CHECK(ncerr)
1176 : end if
1177 :
1178 3279 : if (mpert >= natom+2) then
1179 : ncerr = nf90_put_var(ncid, vid("description_of_perturbations"), &
1180 9486 : trim(DESCR_ipert_2), start=[1, natom+2])
1181 3162 : NCF_CHECK(ncerr)
1182 : end if
1183 :
1184 3279 : if (mpert >= natom+3) then
1185 : ncerr = nf90_put_var(ncid, vid("description_of_perturbations"), &
1186 9486 : trim(DESCR_ipert_3), start=[1, natom+3])
1187 3162 : NCF_CHECK(ncerr)
1188 : end if
1189 :
1190 3279 : if (mpert >= natom+4) then
1191 : ncerr = nf90_put_var(ncid, vid("description_of_perturbations"), &
1192 9486 : trim(DESCR_ipert_4), start=[1, natom+4])
1193 3162 : NCF_CHECK(ncerr)
1194 : end if
1195 :
1196 3279 : if (mpert >= natom+5) then
1197 : ncerr = nf90_put_var(ncid, vid("description_of_perturbations"), &
1198 9486 : trim(DESCR_ipert_5), start=[1, natom+5])
1199 3162 : NCF_CHECK(ncerr)
1200 : end if
1201 :
1202 : ! GA: The numbering natom+5, natom+10 , natom+11 comes from dfpt_looppert
1203 : ! However, we dont want to increase mpert needlessly.
1204 : ! Hence, we use natom+5, natom+6, natom+7.
1205 3279 : if (mpert >= natom+6) then
1206 : ncerr = nf90_put_var(ncid, vid("description_of_perturbations"), &
1207 9486 : trim(DESCR_ipert_10), start=[1, natom+6])
1208 3162 : NCF_CHECK(ncerr)
1209 : end if
1210 :
1211 3279 : if (mpert >= natom+7) then
1212 : ncerr = nf90_put_var(ncid, vid("description_of_perturbations"), &
1213 1260 : trim(DESCR_ipert_11), start=[1, natom+7])
1214 420 : NCF_CHECK(ncerr)
1215 : end if
1216 :
1217 :
1218 3279 : ABI_MALLOC(available_block_types, (ntypes))
1219 : ! The order is unimportant, as long as it matches the description order
1220 :
1221 3279 : available_block_types(1) = BLKTYP_d0E_xx
1222 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1223 6558 : trim(DESCR_d0E_xx), start=[1, 1])
1224 3279 : NCF_CHECK(ncerr)
1225 :
1226 3279 : available_block_types(2) = BLKTYP_d1E_xx
1227 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1228 6558 : trim(DESCR_d1E_xx), start=[1, 2])
1229 3279 : NCF_CHECK(ncerr)
1230 :
1231 3279 : available_block_types(3) = BLKTYP_d2E_ns
1232 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1233 6558 : trim(DESCR_d2E_ns), start=[1, 3])
1234 3279 : NCF_CHECK(ncerr)
1235 :
1236 3279 : available_block_types(4) = BLKTYP_d2E_st
1237 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1238 6558 : trim(DESCR_d2E_st), start=[1, 4])
1239 3279 : NCF_CHECK(ncerr)
1240 :
1241 3279 : available_block_types(5) = BLKTYP_d2E_mbc
1242 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1243 6558 : trim(DESCR_d2E_mbc), start=[1, 5])
1244 3279 : NCF_CHECK(ncerr)
1245 :
1246 3279 : available_block_types(6) = BLKTYP_d3E_xx
1247 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1248 6558 : trim(DESCR_d3E_xx), start=[1, 6])
1249 3279 : NCF_CHECK(ncerr)
1250 :
1251 3279 : available_block_types(7) = BLKTYP_d3E_lw
1252 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1253 6558 : trim(DESCR_d3E_lw), start=[1, 7])
1254 3279 : NCF_CHECK(ncerr)
1255 :
1256 3279 : available_block_types(8) = BLKTYP_d2eig_re
1257 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1258 6558 : trim(DESCR_d2eig_re), start=[1, 8])
1259 3279 : NCF_CHECK(ncerr)
1260 :
1261 3279 : available_block_types(9) = BLKTYP_d2eig_im
1262 : ncerr = nf90_put_var(ncid, vid("description_of_block_types"),&
1263 6558 : trim(DESCR_d2eig_im), start=[1, 9])
1264 3279 : NCF_CHECK(ncerr)
1265 :
1266 6558 : ncerr = nf90_put_var(ncid, vid("available_block_types"), available_block_types)
1267 3279 : NCF_CHECK(ncerr)
1268 3279 : ABI_FREE(available_block_types)
1269 :
1270 6558 : ncerr = nctk_def_dims(ncid, [nctkdim_t('number_of_blocks', ddb_hdr%nblok)])
1271 :
1272 : ! Info on blocks
1273 : ncerr = nctk_def_arrays(ncid, [&
1274 : nctkarr_t('block_types', "i", 'number_of_blocks') &
1275 6558 : &])
1276 3279 : NCF_CHECK(ncerr)
1277 :
1278 3279 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'block_types'), ddb_hdr%typ))
1279 :
1280 : ! TODO Write a chart of block types
1281 :
1282 : ! Count each types of blocks
1283 3279 : nblok_d0E = 0
1284 3279 : nblok_d1E = 0
1285 3279 : nblok_d2E = 0
1286 3279 : nblok_d3E = 0
1287 3279 : nblok_d2eig = 0
1288 :
1289 9287 : do iblok=1,ddb_hdr%nblok
1290 6008 : blktyp = ddb_hdr%typ(iblok)
1291 2712 : if (is_type_d0E(blktyp)) nblok_d0E = nblok_d0E + 1
1292 2719 : if (is_type_d1E(blktyp)) nblok_d1E = nblok_d1E + 1
1293 436 : if (is_type_d2E(blktyp)) nblok_d2E = nblok_d2E + 1
1294 24 : if (is_type_d3E(blktyp)) nblok_d3E = nblok_d3E + 1
1295 9287 : if (is_type_d2eig(blktyp)) nblok_d2eig = nblok_d2eig + 1
1296 : end do
1297 :
1298 : ! ------------------------------------------
1299 : ! Build local arrays of global block indices
1300 : ! to keep track of the initial order.
1301 : ! ------------------------------------------
1302 :
1303 9837 : ABI_MALLOC(blktyp_d0E, (nblok_d0E))
1304 9837 : ABI_MALLOC(blktyp_d1E, (nblok_d1E))
1305 9837 : ABI_MALLOC(blktyp_d2E, (nblok_d2E))
1306 9837 : ABI_MALLOC(blktyp_d3E, (nblok_d3E))
1307 9837 : ABI_MALLOC(blktyp_d2eig, (nblok_d2eig))
1308 :
1309 3279 : iblok_d0E = 0
1310 3279 : iblok_d1E = 0
1311 3279 : iblok_d2E = 0
1312 3279 : iblok_d3E = 0
1313 3279 : iblok_d2eig = 0
1314 :
1315 9287 : do giblok=1,ddb_hdr%nblok
1316 :
1317 6008 : blktyp = ddb_hdr%typ(giblok)
1318 :
1319 3279 : if (is_type_d0E(blktyp)) then
1320 2712 : iblok_d0E = iblok_d0E + 1
1321 2712 : blktyp_d0E(iblok_d0E) = blktyp
1322 :
1323 3296 : else if (is_type_d1E(blktyp)) then
1324 2719 : iblok_d1E = iblok_d1E + 1
1325 2719 : blktyp_d1E(iblok_d1E) = blktyp
1326 :
1327 577 : else if (is_type_d2E(blktyp)) then
1328 436 : iblok_d2E = iblok_d2E + 1
1329 436 : blktyp_d2E(iblok_d2E) = blktyp
1330 :
1331 141 : else if (is_type_d3E(blktyp)) then
1332 24 : iblok_d3E = iblok_d3E + 1
1333 24 : blktyp_d3E(iblok_d3E) = blktyp
1334 :
1335 6008 : else if (is_type_d2eig(blktyp)) then
1336 117 : iblok_d2eig = iblok_d2eig + 1
1337 117 : blktyp_d2eig(iblok_d2eig) = blktyp
1338 : end if
1339 :
1340 : end do
1341 :
1342 : ! ------------------------
1343 : ! Zeroth-order derivatives
1344 : ! ------------------------
1345 :
1346 3279 : NCF_CHECK(nf90_def_grp(ncid, 'd0E', ncid_d0E))
1347 :
1348 6558 : ncerr = nctk_def_dims(ncid_d0E, [nctkdim_t('number_of_d0E_blocks', nblok_d0E)], defmode=.True.)
1349 3279 : NCF_CHECK(ncerr)
1350 :
1351 : ! Info on blocks and matrix values
1352 : ncerr = nctk_def_arrays(ncid_d0E, [&
1353 : nctkarr_t('d0E_block_types', "i", 'number_of_d0E_blocks'), &
1354 : nctkarr_t('matrix_values', "dp", 'number_of_d0E_blocks'), &
1355 : nctkarr_t('matrix_mask', "i", 'number_of_d0E_blocks') &
1356 13116 : &])
1357 3279 : NCF_CHECK(ncerr)
1358 :
1359 3279 : NCF_CHECK(nf90_put_var(ncid_d0E, nctk_idname(ncid_d0E, 'd0E_block_types'), blktyp_d0E))
1360 :
1361 3279 : ABI_FREE(blktyp_d0E)
1362 :
1363 : ! -----------------------
1364 : ! First-order derivatives
1365 : ! -----------------------
1366 :
1367 3279 : NCF_CHECK(nf90_def_grp(ncid, 'd1E', ncid_d1E))
1368 :
1369 6558 : ncerr = nctk_def_dims(ncid_d1E, [nctkdim_t('number_of_d1E_blocks', nblok_d1E)], defmode=.True.)
1370 3279 : NCF_CHECK(ncerr)
1371 :
1372 : ! Info on blocks and matrix values
1373 : ncerr = nctk_def_arrays(ncid_d1E, [&
1374 : nctkarr_t('d1E_block_types', "i", 'number_of_d1E_blocks'),&
1375 : nctkarr_t('matrix_values', "dp",&
1376 : sjoin('cplex,',&
1377 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1378 : 'number_of_d1E_blocks')),&
1379 : nctkarr_t('matrix_mask', "i",&
1380 : sjoin('number_of_cartesian_directions, number_of_perturbations, ',&
1381 : 'number_of_d1E_blocks'))&
1382 13116 : &])
1383 3279 : NCF_CHECK(ncerr)
1384 :
1385 3279 : NCF_CHECK(nf90_put_var(ncid_d1E, nctk_idname(ncid_d1E, 'd1E_block_types'), blktyp_d1E))
1386 :
1387 3279 : ABI_FREE(blktyp_d1E)
1388 :
1389 : ! ------------------------
1390 : ! Second-order derivatives
1391 : ! ------------------------
1392 :
1393 3279 : NCF_CHECK(nf90_def_grp(ncid, 'd2E', ncid_d2E))
1394 :
1395 : ! Qpoints
1396 : ncerr = nctk_def_dims(ncid_d2E, [&
1397 : nctkdim_t('number_of_d2E_blocks', nblok_d2E) &
1398 6558 : ], defmode=.True.)
1399 3279 : NCF_CHECK(ncerr)
1400 :
1401 : ncerr = nctk_def_arrays(ncid_d2E, [&
1402 : nctkarr_t('reduced_coordinates_of_qpoints', "dp",&
1403 : 'number_of_cartesian_directions, number_of_d2E_blocks') &
1404 6558 : &])
1405 3279 : NCF_CHECK(ncerr)
1406 :
1407 : ncerr = nctk_def_arrays(ncid_d2E, [&
1408 : nctkarr_t('qpoints_normalization', "dp", 'number_of_d2E_blocks') &
1409 6558 : &])
1410 3279 : NCF_CHECK(ncerr)
1411 :
1412 : ncerr = nctk_def_arrays(ncid_d2E, [&
1413 : nctkarr_t('frequency', "dp", 'number_of_d2E_blocks') &
1414 6558 : &])
1415 3279 : NCF_CHECK(ncerr)
1416 :
1417 : ! Info on blocks and matrix values
1418 : ncerr = nctk_def_arrays(ncid_d2E, [&
1419 : nctkarr_t('matrix_values', "dp",&
1420 : sjoin('cplex, ',&
1421 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1422 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1423 : 'number_of_d2E_blocks')), &
1424 : nctkarr_t('matrix_mask', "i",&
1425 : sjoin('number_of_cartesian_directions, number_of_perturbations, ',&
1426 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1427 : 'number_of_d2E_blocks')), &
1428 : nctkarr_t('d2E_block_types', "i", 'number_of_d2E_blocks') &
1429 13116 : &])
1430 3279 : NCF_CHECK(ncerr)
1431 :
1432 3279 : NCF_CHECK(nf90_put_var(ncid_d2E, nctk_idname(ncid_d2E, 'd2E_block_types'), blktyp_d2E))
1433 :
1434 3279 : ABI_FREE(blktyp_d2E)
1435 :
1436 : ! -----------------------
1437 : ! Third-order derivatives
1438 : ! -----------------------
1439 :
1440 3279 : NCF_CHECK(nf90_def_grp(ncid, 'd3E', ncid_d3E))
1441 :
1442 : ! Qpoints
1443 : ncerr = nctk_def_dims(ncid_d3E, [&
1444 : nctkdim_t('number_of_d3E_blocks', nblok_d3E) &
1445 6558 : ], defmode=.True.)
1446 3279 : NCF_CHECK(ncerr)
1447 :
1448 : ncerr = nctk_def_arrays(ncid_d3E, [&
1449 : nctkarr_t('reduced_coordinates_of_qpoints', "dp",&
1450 : 'three_dim, number_of_cartesian_directions, number_of_d3E_blocks') &
1451 6558 : &])
1452 3279 : NCF_CHECK(ncerr)
1453 :
1454 : ncerr = nctk_def_arrays(ncid_d3E, [&
1455 : nctkarr_t('qpoints_normalization', "dp",&
1456 : 'three_dim, number_of_d3E_blocks') &
1457 6558 : &])
1458 3279 : NCF_CHECK(ncerr)
1459 :
1460 : ncerr = nctk_def_arrays(ncid_d3E, [&
1461 : nctkarr_t('frequency', "dp",&
1462 : 'three_dim, number_of_d3E_blocks') &
1463 6558 : &])
1464 3279 : NCF_CHECK(ncerr)
1465 :
1466 : ! Info on blocks and matrix values
1467 : ncerr = nctk_def_arrays(ncid_d3E, [&
1468 : nctkarr_t('matrix_values', "dp",&
1469 : sjoin('cplex, ',&
1470 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1471 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1472 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1473 : 'number_of_d3E_blocks')), &
1474 : nctkarr_t('matrix_mask', "i",&
1475 : sjoin('number_of_cartesian_directions, number_of_perturbations, ',&
1476 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1477 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1478 : 'number_of_d3E_blocks')), &
1479 : nctkarr_t('d3E_block_types', "i", 'number_of_d3E_blocks') &
1480 13116 : &])
1481 3279 : NCF_CHECK(ncerr)
1482 :
1483 3279 : NCF_CHECK(nf90_put_var(ncid_d3E, nctk_idname(ncid_d3E, 'd3E_block_types'), blktyp_d3E))
1484 :
1485 3279 : ABI_FREE(blktyp_d3E)
1486 :
1487 : ! ---------------------------------------
1488 : ! Second-order derivatives of eigenvalues
1489 : ! ---------------------------------------
1490 :
1491 3279 : NCF_CHECK(nf90_def_grp(ncid, 'd2eig', ncid_d2eig))
1492 :
1493 : ! Qpoints
1494 6558 : ncerr = nctk_def_dims(ncid_d2eig, [nctkdim_t('number_of_d2eig_blocks', nblok_d2eig)], defmode=.True.)
1495 3279 : NCF_CHECK(ncerr)
1496 :
1497 : ncerr = nctk_def_arrays(ncid_d2eig, [&
1498 : nctkarr_t('reduced_coordinates_of_qpoints', "dp",&
1499 : 'number_of_cartesian_directions, number_of_d2eig_blocks') &
1500 6558 : &])
1501 3279 : NCF_CHECK(ncerr)
1502 :
1503 : ncerr = nctk_def_arrays(ncid_d2eig, [&
1504 : nctkarr_t('qpoints_normalization', "dp", 'number_of_d2eig_blocks') &
1505 6558 : &])
1506 3279 : NCF_CHECK(ncerr)
1507 :
1508 : ! Info on blocks and matrix values
1509 : ncerr = nctk_def_arrays(ncid_d2eig, [&
1510 : nctkarr_t('matrix_values', "dp",&
1511 : sjoin('cplex, ',&
1512 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1513 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1514 : 'maximum_number_of_bands, number_of_kpoints, number_of_spins, ',&
1515 : 'number_of_d2eig_blocks')), &
1516 : nctkarr_t('matrix_mask', "i",&
1517 : sjoin('number_of_cartesian_directions, number_of_perturbations, ',&
1518 : 'number_of_cartesian_directions, number_of_perturbations, ',&
1519 : 'number_of_d2eig_blocks')), &
1520 : nctkarr_t('d2eig_block_types', "i", 'number_of_d2eig_blocks') &
1521 13116 : &])
1522 3279 : NCF_CHECK(ncerr)
1523 :
1524 3279 : NCF_CHECK(nf90_put_var(ncid_d2eig, nctk_idname(ncid_d2eig, 'd2eig_block_types'), blktyp_d2eig))
1525 :
1526 3279 : ABI_FREE(blktyp_d2eig)
1527 :
1528 :
1529 : contains
1530 82655 : integer function vid(vname)
1531 : character(len=*),intent(in) :: vname
1532 82655 : vid = nctk_idname(ncid, vname)
1533 : end function vid
1534 :
1535 : end subroutine ddb_hdr_open_write_nc
1536 : !!***
1537 :
1538 : !----------------------------------------------------------------------
1539 :
1540 : !!****f* m_ddb_hdr/ddb_hdr_open_write
1541 : !! NAME
1542 : !! ddb_hdr_open_write
1543 : !!
1544 : !! FUNCTION
1545 : !! Open the DDB file and write the header.
1546 : !!
1547 : !! INPUTS
1548 : !! filename=name of the file being written
1549 : !! with_psps
1550 : !! 1-> include information on pseudopoentials
1551 : !! 0-> do not include information on pseudopoentials
1552 : !! comm=MPI communicator.
1553 : !!
1554 : !! OUTPUT
1555 : !!
1556 : !! SOURCE
1557 :
1558 143 : subroutine ddb_hdr_open_write(ddb_hdr, filename, with_psps, comm)
1559 :
1560 : !Arguments ------------------------------------
1561 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
1562 : character(len=*),intent(in) :: filename
1563 : integer,intent(in),optional :: with_psps
1564 : integer,intent(in),optional :: comm
1565 :
1566 : !Local variables-------------------------------
1567 : !scalars
1568 : integer,parameter :: master=0
1569 : integer :: comm_, iomode
1570 : character(len=fnlen) :: filename_
1571 : ! ************************************************************************
1572 :
1573 143 : if (present(comm)) then
1574 118 : comm_ = comm
1575 : else
1576 25 : comm_ = xmpi_comm_self
1577 : end if
1578 :
1579 143 : if (xmpi_comm_rank(comm_) /= master) return
1580 :
1581 123 : call ddb_hdr%get_iomode(filename, 2, iomode, filename_)
1582 :
1583 123 : if (iomode==IO_MODE_ETSF) then
1584 117 : call ddb_hdr%open_write_nc(nctk_ncify(filename), with_psps=with_psps)
1585 6 : else if (iomode==IO_MODE_FORTRAN) then
1586 6 : call ddb_hdr%open_write_txt(filename, with_psps)
1587 : end if
1588 :
1589 : end subroutine ddb_hdr_open_write
1590 : !!***
1591 :
1592 : !----------------------------------------------------------------------
1593 :
1594 : !!****f* m_ddb_hdr/ddb_hdr_open_read
1595 : !! NAME
1596 : !! ddb_hdr_open_read
1597 : !!
1598 : !! FUNCTION
1599 : !! Open the DDB file and read the header.
1600 : !!
1601 : !! INPUTS
1602 : !! filename=name of the file being written (abo_DS*_DDB)
1603 : !! comm=MPI communicator
1604 : !! matom,mtypat,mband,mkpt,msym,dimekb,lmnmax,usepaw=dimensions to be set.
1605 : !! dimonly
1606 : !! 1-> only read dimensions and close the file
1607 : !! 0-> leave the file open
1608 : !!
1609 : !! SOURCE
1610 :
1611 706 : subroutine ddb_hdr_open_read(ddb_hdr, filename, comm, &
1612 : matom,mtypat,mband,mkpt,msym,dimekb,lmnmax,usepaw,dimonly)
1613 :
1614 : !Arguments ------------------------------------
1615 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
1616 : character(len=*),intent(in) :: filename
1617 : integer,intent(in) :: comm
1618 : integer,intent(in),optional :: matom,mtypat,mband,mkpt,msym
1619 : integer,intent(in),optional :: dimekb,lmnmax,usepaw
1620 : integer,intent(in),optional :: dimonly
1621 :
1622 : !Local variables-------------------------------
1623 : !scalars
1624 : integer,parameter :: master=0
1625 : integer :: iomode
1626 : character(len=fnlen) :: filename_
1627 :
1628 706 : call ddb_hdr%get_iomode(filename, 1, iomode, filename_)
1629 :
1630 706 : if (xmpi_comm_rank(comm) == master) then
1631 : ! GA: FIXME
1632 : ! Here we keep the old output behavior just to make test pass.
1633 : ! However, test should allow differences in filenames.
1634 : !call wrtout(std_out, sjoin(" Opening DDB file:", filename_), 'COLL')
1635 690 : call wrtout(std_out, sjoin(" Opening DDB file:", filename), 'COLL')
1636 : end if
1637 :
1638 706 : if (iomode==IO_MODE_ETSF) then
1639 :
1640 : call ddb_hdr%open_read_nc(filename_, comm, &
1641 392 : matom,mtypat,mband,mkpt,msym,dimekb,lmnmax,usepaw,dimonly)
1642 :
1643 314 : else if (iomode==IO_MODE_FORTRAN) then
1644 :
1645 : call ddb_hdr%open_read_txt(filename, comm, &
1646 314 : matom,mtypat,mband,mkpt,msym,dimekb,lmnmax,usepaw,dimonly)
1647 :
1648 : end if
1649 :
1650 :
1651 706 : end subroutine ddb_hdr_open_read
1652 : !!***
1653 :
1654 : !----------------------------------------------------------------------
1655 :
1656 : !!****f* m_ddb_hdr/ddb_hdr_open_read_txt
1657 : !! NAME
1658 : !! ddb_hdr_open_read_txt
1659 : !!
1660 : !! FUNCTION
1661 : !! Open the DDB text file and read the header.
1662 : !! If dimonly, only read dimensions and close the file,
1663 : !! otherwise leave the file open.
1664 : !!
1665 : !! INPUTS
1666 : !! filename=name of the file being written (abo_DS*_DDB)
1667 : !! comm=MPI communicator
1668 : !! matom,mtypat,mband,mkpt,msym,dimekb,lmnmax,usepaw=dimensions to be set.
1669 : !! dimonly
1670 : !! 1-> only read dimensions and close the file
1671 : !! 0-> leave the file open
1672 : !!
1673 : !! SOURCE
1674 :
1675 1220 : subroutine ddb_hdr_open_read_txt(ddb_hdr, filename, comm, &
1676 : & matom,mtypat,mband,mkpt,msym,dimekb,lmnmax,usepaw,dimonly)
1677 :
1678 : !Arguments ------------------------------------
1679 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
1680 : character(len=*),intent(in) :: filename
1681 : integer,intent(in),optional :: comm
1682 : integer,intent(in),optional :: matom,mtypat,mband,mkpt,msym
1683 : integer,intent(in),optional :: dimekb,lmnmax,usepaw
1684 : integer,intent(in),optional :: dimonly
1685 :
1686 : !Local variables -------------------------
1687 : !scalars
1688 : integer,parameter :: master=0
1689 : integer :: unddb
1690 : integer :: spgroup, timrev, npsp
1691 : integer :: choice
1692 : integer :: ii, nblok
1693 : integer :: matom_l,mtypat_l,mband_l,mkpt_l,msym_l,nsppol_l
1694 : integer :: dimekb_l,lmnmax_l,usepaw_l
1695 : integer :: comm_l
1696 : logical :: has_d0E_xx,has_d1E_xx,has_d2E_xx,has_d3E_xx,has_d3E_lw,has_d2eig
1697 : !arrays
1698 : real(dp):: rprimd(3,3)
1699 : ! ************************************************************************
1700 :
1701 610 : if (present(comm)) then
1702 610 : comm_l = comm
1703 : else
1704 0 : comm_l = xmpi_comm_self
1705 : end if
1706 :
1707 : ! =====================================================
1708 : ! Open the DDB once and read the dimensions from header
1709 : ! =====================================================
1710 :
1711 : call ddb_getdims(filename,comm_l,dimekb_l,lmnmax_l,mband_l, &
1712 : & msym_l,matom_l,nblok,mkpt_l,nsppol_l,mtypat_l,usepaw_l,&
1713 610 : & has_d0E_xx, has_d1E_xx, has_d2E_xx, has_d3E_xx, has_d3E_lw, has_d2eig)
1714 :
1715 :
1716 : ! When merging ddbs, we want to open a header with fixed dimensions.
1717 : ! GA: Might not be necessary if we handle the merging better.
1718 : !
1719 610 : if (present(matom)) matom_l = matom
1720 610 : if (present(mtypat)) mtypat_l = mtypat
1721 610 : if (present(mband)) mband_l = mband
1722 610 : if (present(mkpt)) mkpt_l = mkpt
1723 610 : if (present(msym)) msym_l = msym
1724 610 : if (present(dimekb)) dimekb_l = dimekb
1725 610 : if (present(lmnmax)) lmnmax_l = lmnmax
1726 610 : if (present(usepaw)) usepaw_l = usepaw
1727 :
1728 610 : ddb_hdr%ddb_version = DDB_VERSION
1729 610 : ddb_hdr%usepaw = usepaw_l
1730 610 : ddb_hdr%mband = mband_l
1731 610 : ddb_hdr%matom = matom_l
1732 610 : ddb_hdr%msym = msym_l
1733 610 : ddb_hdr%mtypat = mtypat_l
1734 610 : ddb_hdr%mkpt = mkpt_l
1735 610 : ddb_hdr%nsppol = nsppol_l
1736 :
1737 610 : ddb_hdr%nblok = nblok
1738 610 : ddb_hdr%has_d0E_xx = has_d0E_xx
1739 610 : ddb_hdr%has_d1E_xx = has_d1E_xx
1740 610 : ddb_hdr%has_d2E_xx = has_d2E_xx
1741 610 : ddb_hdr%has_d3E_xx = has_d3E_xx
1742 610 : ddb_hdr%has_d3E_lw = has_d3E_lw
1743 610 : ddb_hdr%has_d2eig = has_d2eig
1744 :
1745 610 : ddb_hdr%psps%dimekb = dimekb_l
1746 610 : ddb_hdr%psps%ntypat = mtypat_l
1747 610 : ddb_hdr%psps%lmnmax = lmnmax_l
1748 610 : ddb_hdr%psps%usepaw = usepaw_l
1749 610 : ddb_hdr%psps%useylm = usepaw_l ! yep.
1750 :
1751 : ! GA: Eventually, we might want to initialize other entries of psps here.
1752 : ! such as alchemical psps info.
1753 :
1754 610 : ddb_hdr%natom = ddb_hdr%matom
1755 610 : ddb_hdr%nkpt = ddb_hdr%mkpt
1756 610 : ddb_hdr%nsym = ddb_hdr%msym
1757 610 : ddb_hdr%ntypat = ddb_hdr%mtypat
1758 :
1759 610 : npsp = ddb_hdr%mtypat
1760 :
1761 : ! Set maximal value for mpert
1762 610 : ddb_hdr%mpert = 2*ddb_hdr%natom+MPERT_MAX
1763 :
1764 : ! Compute the block dimensions
1765 610 : call ddb_hdr%get_block_dims()
1766 :
1767 610 : unddb = get_unit()
1768 610 : ddb_hdr%unddb = unddb
1769 :
1770 610 : if (present(dimonly)) then
1771 216 : if (dimonly==1) return
1772 : end if
1773 :
1774 :
1775 : ! =========================================
1776 : ! Open the DDB a second time to read header
1777 : ! =========================================
1778 :
1779 394 : ddb_hdr%has_open_file_txt = .true.
1780 :
1781 394 : if (xmpi_comm_rank(comm) == master) then
1782 :
1783 : ! Allocate the memory
1784 386 : call ddb_hdr%malloc()
1785 :
1786 1544 : ABI_MALLOC(ddb_hdr%psps%indlmn,(6,ddb_hdr%psps%lmnmax,ddb_hdr%mtypat))
1787 1158 : ABI_MALLOC(ddb_hdr%psps%pspso,(ddb_hdr%mtypat))
1788 1544 : ABI_MALLOC(ddb_hdr%psps%ekb,(ddb_hdr%psps%dimekb,ddb_hdr%mtypat))
1789 :
1790 1169 : ABI_MALLOC(ddb_hdr%pawtab,(ddb_hdr%psps%ntypat*ddb_hdr%psps%usepaw))
1791 386 : call pawtab_nullify(ddb_hdr%pawtab)
1792 :
1793 : ! This is needed to read the DDBs in the old format
1794 : ! GA : Not sure if we really need this.
1795 74498 : ddb_hdr%symafm(:)=1
1796 386 : if(ddb_hdr%mtypat>=1)then
1797 1125 : ddb_hdr%psps%pspso(:)=0
1798 4848 : ddb_hdr%psps%ekb(:,:)=zero
1799 1125 : ddb_hdr%znucl(:)=zero
1800 : end if
1801 386 : if(ddb_hdr%matom>=1)then
1802 4862 : ddb_hdr%spinat(:,:)=zero
1803 : end if
1804 :
1805 : ! Note: the maximum parameters (matom, mkpt, etc.) are inputs to ioddb8_in
1806 : ! wile the actual parameters (natom, nkpt, etc.) are outputs
1807 : call ioddb8_in(filename,ddb_hdr%matom,ddb_hdr%mband,&
1808 : & ddb_hdr%mkpt,ddb_hdr%msym,ddb_hdr%mtypat,unddb,&
1809 : & ddb_hdr%acell,ddb_hdr%amu,ddb_hdr%ddb_version,ddb_hdr%dilatmx,ddb_hdr%ecut,ddb_hdr%ecutsm,&
1810 : & ddb_hdr%intxc,ddb_hdr%iscf,ddb_hdr%ixc,ddb_hdr%kpt,ddb_hdr%kptnrm,&
1811 : & ddb_hdr%natom,ddb_hdr%nband,ddb_hdr%ngfft,ddb_hdr%nkpt,ddb_hdr%nspden,&
1812 : & ddb_hdr%nspinor,ddb_hdr%nsppol,ddb_hdr%nsym,ddb_hdr%ntypat,&
1813 : & ddb_hdr%occ,ddb_hdr%occopt,ddb_hdr%pawecutdg,ddb_hdr%rprim,&
1814 : & ddb_hdr%dfpt_sciss,ddb_hdr%spinat,ddb_hdr%symafm,ddb_hdr%symrel,&
1815 : & ddb_hdr%tnons,ddb_hdr%tolwfr,ddb_hdr%tphysel,ddb_hdr%tsmear,&
1816 : & ddb_hdr%typat,ddb_hdr%usepaw,ddb_hdr%wtk,ddb_hdr%xred,ddb_hdr%zion,&
1817 386 : & ddb_hdr%znucl)
1818 :
1819 :
1820 : ! Read the psp information of the input DDB
1821 386 : choice=1 ! Read
1822 : call psddb8(choice,ddb_hdr%psps%dimekb,ddb_hdr%psps%ekb,ddb_hdr%with_psps,&
1823 : & ddb_hdr%psps%indlmn,ddb_hdr%psps%lmnmax,&
1824 : & ddb_hdr%nblok,ddb_hdr%ntypat,unddb,ddb_hdr%pawtab,ddb_hdr%psps%pspso,&
1825 386 : & ddb_hdr%psps%usepaw,ddb_hdr%psps%useylm)
1826 :
1827 :
1828 : ! This is a limitation of the text file format.
1829 : ! We might be able to get rid of it.
1830 386 : ddb_hdr%with_dfpt_vars = ddb_hdr%with_psps
1831 :
1832 386 : ddb_hdr%dscrpt = ''
1833 :
1834 : ! Complete the initialization of the psps object
1835 : ! so that we may copy it.
1836 386 : ddb_hdr%psps%npsp = npsp
1837 386 : ddb_hdr%psps%ntypat = ddb_hdr%ntypat
1838 386 : ddb_hdr%psps%ntypalch = 0 ! Alchemical potential not supported
1839 386 : ddb_hdr%psps%mtypalch = 0 ! Alchemical potential not supported
1840 386 : ddb_hdr%psps%n1xccc = 0
1841 386 : ddb_hdr%psps%usewvl = 0
1842 386 : ddb_hdr%psps%mqgrid_vl = 0
1843 386 : ddb_hdr%psps%mqgrid_ff = 0
1844 386 : ddb_hdr%psps%mpspso = 1 ! FIXME
1845 1125 : do ii=1,npsp
1846 1125 : ddb_hdr%psps%mpspso = max(ddb_hdr%psps%mpspso, ddb_hdr%psps%pspso(ii))
1847 : end do
1848 386 : ddb_hdr%psps%mpssoang = 1
1849 386 : ddb_hdr%psps%mpsang = 1
1850 386 : ddb_hdr%psps%mproj = 0
1851 386 : ddb_hdr%psps%lnmax = 0
1852 386 : ddb_hdr%psps%ntyppure = 0
1853 386 : ddb_hdr%psps%npspalch = 0
1854 386 : ddb_hdr%psps%nc_xccc_gspace = 0
1855 386 : ddb_hdr%psps%optnlxccc = 0
1856 386 : ddb_hdr%psps%positron = 0
1857 386 : ddb_hdr%psps%vlspl_recipSpace = .false.
1858 :
1859 1158 : ABI_MALLOC(ddb_hdr%psps%ziontypat,(ddb_hdr%psps%ntypat))
1860 1158 : ABI_MALLOC(ddb_hdr%psps%znucltypat,(ddb_hdr%psps%ntypat))
1861 1158 : ABI_MALLOC(ddb_hdr%psps%filpsp,(ddb_hdr%psps%npsp))
1862 1158 : ABI_MALLOC(ddb_hdr%psps%md5_pseudos,(ddb_hdr%psps%npsp))
1863 1158 : ABI_MALLOC(ddb_hdr%psps%title,(ddb_hdr%psps%npsp))
1864 1158 : ABI_MALLOC(ddb_hdr%psps%pspxc,(ddb_hdr%psps%ntypat))
1865 :
1866 1158 : ABI_MALLOC(ddb_hdr%psps%pspcod,(ddb_hdr%psps%npsp))
1867 1158 : ABI_MALLOC(ddb_hdr%psps%algalch,(ddb_hdr%psps%ntypalch))
1868 1158 : ABI_MALLOC(ddb_hdr%psps%pspdat,(ddb_hdr%psps%ntypat))
1869 1158 : ABI_MALLOC(ddb_hdr%psps%zionpsp,(ddb_hdr%psps%npsp))
1870 1158 : ABI_MALLOC(ddb_hdr%psps%znuclpsp,(ddb_hdr%psps%npsp))
1871 :
1872 1897 : ABI_MALLOC(ddb_hdr%psps%nctab,(ddb_hdr%psps%ntypat))
1873 :
1874 1125 : ddb_hdr%psps%ziontypat(:) = ddb_hdr%zion(:)
1875 1125 : ddb_hdr%psps%znucltypat(:) = ddb_hdr%znucl(:)
1876 1125 : do ii=1,npsp
1877 739 : ddb_hdr%psps%filpsp(ii) = ''
1878 739 : ddb_hdr%psps%md5_pseudos(ii) = ''
1879 1125 : ddb_hdr%psps%title(ii) = ''
1880 : !ddb_hdr%psps%pspxc(ii) = ddb_hdr%ixc
1881 : end do
1882 :
1883 1125 : ddb_hdr%psps%pspxc = zero ! Unknown
1884 :
1885 : ! -------------------------------
1886 : ! Initialize the crystal
1887 : ! -------------------------------
1888 :
1889 386 : call mkrdim(ddb_hdr%acell,ddb_hdr%rprim,rprimd)
1890 :
1891 : ! GA: space group and time reversal are not written in the text file
1892 : ! but they are written in the netcdf version.
1893 : ! It doesnt seem to make any difference to anaddb.
1894 386 : spgroup = 1
1895 386 : timrev = 2
1896 :
1897 : call ddb_hdr%crystal%init(ddb_hdr%amu, &
1898 : & spgroup, ddb_hdr%natom, npsp, ddb_hdr%ntypat, &
1899 : & ddb_hdr%nsym, rprimd, ddb_hdr%typat, &
1900 : & ddb_hdr%xred, ddb_hdr%zion, ddb_hdr%znucl, timrev, &
1901 : & ddb_hdr%nspden==2.and.ddb_hdr%nsppol==1, .false., '', &
1902 772 : & ddb_hdr%symrel, ddb_hdr%tnons, ddb_hdr%symafm)
1903 :
1904 :
1905 : end if
1906 :
1907 : ! Master broadcasts data
1908 394 : call ddb_hdr%bcast(comm_l)
1909 :
1910 : end subroutine ddb_hdr_open_read_txt
1911 : !!***
1912 :
1913 : !----------------------------------------------------------------------
1914 :
1915 : !!****f* m_ddb_hdr/ddb_hdr_open_read_nc
1916 : !! NAME
1917 : !! ddb_hdr_open_read
1918 : !!
1919 : !! FUNCTION
1920 : !! Open the DDB.nc file and read the header.
1921 : !!
1922 : !! INPUTS
1923 : !! filename=name of the file being written (abo_DS*_DDB.nc)
1924 : !! comm=MPI communicator
1925 : !! matom,mtypat,mband,mkpt,msym,dimekb,lmnmax,usepaw=dimensions to be set.
1926 : !! dimonly
1927 : !! 1-> only read dimensions and close the file
1928 : !! 0-> leave the file open
1929 : !!
1930 : !! SOURCE
1931 :
1932 563 : subroutine ddb_hdr_open_read_nc(ddb_hdr, filename, comm, &
1933 : & matom,mtypat,mband,mkpt,msym,dimekb,lmnmax,usepaw,dimonly)
1934 :
1935 : !Arguments ------------------------------------
1936 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
1937 : character(len=*),intent(in) :: filename
1938 : integer,intent(in),optional :: comm
1939 : integer,intent(in),optional :: matom,mtypat,mband,mkpt,msym
1940 : integer,intent(in),optional :: dimekb,lmnmax,usepaw
1941 : integer,intent(in),optional :: dimonly
1942 :
1943 : !Local variables-------------------------------
1944 : !scalars
1945 : integer,parameter :: master=0
1946 : integer :: ncid
1947 : integer :: ncid_crystal, ncid_psps, ncid_pawtab
1948 : integer :: ii,isppol,isym,ikpt,iband,itypat,iblok
1949 : integer :: comm_
1950 : integer :: natom_file, ntypat_file, mband_file, nkpt_file, nsym_file
1951 : integer :: mband_, nkpt_
1952 : integer :: blktyp
1953 : integer :: ncerr
1954 : !arrays
1955 : integer :: ngfft(3)
1956 563 : integer,allocatable :: nband(:,:)
1957 563 : integer,allocatable :: pawtab_basis_size(:)
1958 563 : integer,allocatable :: pawtab_lmn_size(:)
1959 563 : integer,allocatable :: pawtab_lmn2_size(:)
1960 563 : integer,allocatable :: pawtab_shape_type(:)
1961 563 : real(dp),allocatable :: pawtab_rpaw(:)
1962 563 : real(dp),allocatable :: pawtab_rshp(:)
1963 563 : real(dp),allocatable :: occ(:,:,:)
1964 : ! ************************************************************************
1965 :
1966 563 : if (present(comm)) then
1967 563 : comm_ = comm
1968 : else
1969 0 : comm_ = xmpi_comm_self
1970 : end if
1971 :
1972 563 : if (xmpi_comm_rank(comm_) == master) then
1973 563 : NCF_CHECK(nctk_open_read(ncid, filename, xmpi_comm_self))
1974 563 : ddb_hdr%ncid = ncid
1975 : end if
1976 :
1977 563 : if (xmpi_comm_size(comm_) > 1) then
1978 0 : call xmpi_bcast(ddb_hdr%ncid, master, comm_, ncerr)
1979 : end if
1980 :
1981 563 : ddb_hdr%has_open_file_nc = .True.
1982 :
1983 : ! Flush allocated arrays
1984 563 : call ddb_hdr%free()
1985 :
1986 : ! ===============
1987 : ! Read dimensions
1988 : ! ===============
1989 563 : if (xmpi_comm_rank(comm_) == master) then
1990 :
1991 : ! No check is performed on the ddb version number, yet.
1992 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'ddb_version'), ddb_hdr%ddb_version))
1993 :
1994 : ! Scalar integers
1995 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_atoms", natom_file))
1996 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_kpoints", nkpt_file))
1997 563 : ddb_hdr%natom = natom_file
1998 563 : ddb_hdr%nkpt = nkpt_file
1999 :
2000 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_spins", ddb_hdr%nsppol))
2001 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_spinor_components", ddb_hdr%nspinor))
2002 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_spin_densities", ddb_hdr%nspden))
2003 :
2004 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_symmetries", nsym_file))
2005 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_atom_types", ntypat_file))
2006 563 : ddb_hdr%nsym = nsym_file
2007 563 : ddb_hdr%ntypat = ntypat_file
2008 :
2009 563 : NCF_CHECK(nctk_get_dim(ncid, "maximum_number_of_bands", mband_file))
2010 563 : ddb_hdr%mband = mband_file
2011 :
2012 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_blocks", ddb_hdr%nblok))
2013 563 : NCF_CHECK(nctk_get_dim(ncid, "number_of_perturbations", ddb_hdr%mpert))
2014 :
2015 563 : ddb_hdr%msym = ddb_hdr%nsym
2016 563 : ddb_hdr%mkpt = ddb_hdr%nkpt
2017 563 : ddb_hdr%matom = ddb_hdr%natom
2018 563 : ddb_hdr%mtypat = ddb_hdr%ntypat
2019 :
2020 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "with_psps"), ddb_hdr%with_psps))
2021 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "with_dfpt_vars"), ddb_hdr%with_dfpt_vars))
2022 :
2023 : ! Read pseudopotentials info
2024 563 : ncid_psps = nctk_idgroup(ncid, 'Pseudopotentials')
2025 :
2026 563 : call psps_ncread(ddb_hdr%psps, ncid_psps)
2027 563 : ddb_hdr%usepaw = ddb_hdr%psps%usepaw
2028 :
2029 1703 : ABI_MALLOC(ddb_hdr%pawtab, (ddb_hdr%ntypat*ddb_hdr%usepaw))
2030 563 : if (ddb_hdr%usepaw > 0) then
2031 :
2032 : ! Read PAW info
2033 : ! -------------------------------------
2034 : ! GA: TODO There should be a function pawtab_ncread
2035 : ! Here, we only read the info needed to print the ddb in text format.
2036 :
2037 7 : ncid_pawtab = nctk_idgroup(ncid_psps, 'PAW')
2038 7 : call pawtab_nullify(ddb_hdr%pawtab)
2039 :
2040 21 : ABI_MALLOC(pawtab_basis_size, (ddb_hdr%ntypat))
2041 14 : ABI_MALLOC(pawtab_lmn_size, (ddb_hdr%ntypat))
2042 14 : ABI_MALLOC(pawtab_lmn2_size, (ddb_hdr%ntypat))
2043 14 : ABI_MALLOC(pawtab_shape_type, (ddb_hdr%ntypat))
2044 :
2045 7 : NCF_CHECK(nf90_get_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'basis_size'), pawtab_basis_size))
2046 7 : NCF_CHECK(nf90_get_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'lmn_size'), pawtab_lmn_size))
2047 7 : NCF_CHECK(nf90_get_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'lmn2_size'), pawtab_lmn2_size))
2048 7 : NCF_CHECK(nf90_get_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'shape_type'), pawtab_shape_type))
2049 :
2050 21 : ABI_MALLOC(pawtab_rpaw, (ddb_hdr%ntypat))
2051 7 : ncerr = nf90_get_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'r_PAW'), pawtab_rpaw)
2052 7 : NCF_CHECK(ncerr)
2053 21 : ABI_MALLOC(pawtab_rshp, (ddb_hdr%ntypat))
2054 7 : ncerr = nf90_get_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'r_shape'), pawtab_rshp)
2055 7 : NCF_CHECK(ncerr)
2056 :
2057 21 : do itypat=1,ddb_hdr%ntypat
2058 :
2059 14 : ddb_hdr%pawtab(itypat)%basis_size = pawtab_basis_size(itypat)
2060 14 : ddb_hdr%pawtab(itypat)%lmn_size = pawtab_lmn_size(itypat)
2061 14 : ddb_hdr%pawtab(itypat)%lmn2_size = pawtab_lmn2_size(itypat)
2062 14 : ddb_hdr%pawtab(itypat)%shape_type = pawtab_shape_type(itypat)
2063 :
2064 14 : ddb_hdr%pawtab(itypat)%rpaw = pawtab_rpaw(itypat)
2065 14 : ddb_hdr%pawtab(itypat)%rshp = pawtab_rshp(itypat)
2066 :
2067 42 : ABI_MALLOC(ddb_hdr%pawtab(itypat)%dij0, (ddb_hdr%pawtab(itypat)%lmn2_size))
2068 : ncerr = nf90_get_var(ncid_pawtab, nctk_idname(ncid_pawtab, 'Dij0'), &
2069 : ddb_hdr%pawtab(itypat)%dij0, &
2070 : start=[itypat,1], &
2071 70 : count=[1,ddb_hdr%pawtab(itypat)%lmn2_size])
2072 21 : NCF_CHECK(ncerr)
2073 :
2074 : end do
2075 :
2076 7 : ABI_FREE(pawtab_rpaw)
2077 7 : ABI_FREE(pawtab_rshp)
2078 7 : ABI_FREE(pawtab_basis_size)
2079 7 : ABI_FREE(pawtab_lmn_size)
2080 7 : ABI_FREE(pawtab_lmn2_size)
2081 7 : ABI_FREE(pawtab_shape_type)
2082 :
2083 : ! -------------------------------------
2084 :
2085 : end if
2086 :
2087 : ! Read block types and compute maximum block type
2088 563 : ABI_SFREE(ddb_hdr%typ)
2089 1689 : ABI_MALLOC(ddb_hdr%typ,(ddb_hdr%nblok))
2090 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'block_types'), ddb_hdr%typ))
2091 :
2092 1190 : do iblok=1,ddb_hdr%nblok
2093 627 : blktyp = ddb_hdr%typ(iblok)
2094 563 : if (is_type_d0E(blktyp)) then
2095 24 : ddb_hdr%has_d0E_xx=.true.
2096 603 : else if (is_type_d1E(blktyp)) then
2097 36 : ddb_hdr%has_d1E_xx=.true.
2098 567 : else if (is_type_d2E(blktyp)) then
2099 468 : ddb_hdr%has_d2E_xx=.true.
2100 99 : else if (is_type_d3E(blktyp)) then
2101 9 : ddb_hdr%has_d3E_xx=.true.
2102 9 : if (blktyp == BLKTYP_d3E_lw) then
2103 0 : ddb_hdr%has_d3E_lw=.true.
2104 : end if
2105 627 : else if (is_type_d2eig(blktyp)) then
2106 90 : ddb_hdr%has_d2eig=.true.
2107 : end if
2108 : end do
2109 :
2110 : ! Overwrite dimensions if requested
2111 : ! GA: TODO Try to remove these optional arguments
2112 : ! Should at least give a warning when inconsistencies are found
2113 563 : if (present(matom)) ddb_hdr%matom = matom
2114 563 : if (present(mtypat)) ddb_hdr%mtypat = mtypat
2115 563 : if (present(mband)) ddb_hdr%mband = mband
2116 563 : if (present(mkpt)) ddb_hdr%mkpt = mkpt
2117 563 : if (present(msym)) ddb_hdr%msym = msym
2118 563 : if (present(dimekb)) ddb_hdr%psps%dimekb = dimekb
2119 563 : if (present(lmnmax)) ddb_hdr%psps%lmnmax = lmnmax
2120 563 : if (present(usepaw)) ddb_hdr%usepaw = usepaw
2121 :
2122 : ! Compute the block dimensions (mpert and msize)
2123 563 : call ddb_hdr%get_block_dims()
2124 :
2125 : end if
2126 :
2127 563 : ddb_hdr%dscrpt = ''
2128 :
2129 563 : if (present(dimonly)) then
2130 166 : if (dimonly==1) then
2131 166 : call ddb_hdr%bcast_dim(comm_)
2132 166 : call ddb_hdr%close()
2133 166 : return
2134 : end if
2135 : end if
2136 :
2137 : ! ===============
2138 : ! Allocate arrays
2139 : ! ===============
2140 397 : if (xmpi_comm_rank(comm_) == master) then
2141 :
2142 397 : call ddb_hdr%malloc()
2143 :
2144 : !! GA : Not sure we still need this.
2145 : !ddb_hdr%symafm(:)=1
2146 : !if(ddb_hdr%mtypat>=1)then
2147 : ! ddb_hdr%znucl(:)=zero
2148 : !end if
2149 : !if(ddb_hdr%matom>=1)then
2150 : ! ddb_hdr%spinat(:,:)=zero
2151 : !end if
2152 :
2153 : end if
2154 :
2155 : ! ==============
2156 : ! Read variables
2157 : ! ==============
2158 397 : if (xmpi_comm_rank(comm_) == master) then
2159 :
2160 : ! -------
2161 : ! Scalars
2162 : ! -------
2163 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'intxc'), ddb_hdr%intxc))
2164 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'iscf'), ddb_hdr%iscf))
2165 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'ixc'), ddb_hdr%ixc))
2166 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'occopt'), ddb_hdr%occopt))
2167 :
2168 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'dilatmx'), ddb_hdr%dilatmx))
2169 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'ecut'), ddb_hdr%ecut))
2170 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'pawecutdg'), ddb_hdr%pawecutdg))
2171 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'ecutsm'), ddb_hdr%ecutsm))
2172 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'kptnrm'), ddb_hdr%kptnrm))
2173 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'dfpt_sciss'), ddb_hdr%dfpt_sciss))
2174 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'tolwfr'), ddb_hdr%tolwfr))
2175 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'tphysel'), ddb_hdr%tphysel))
2176 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'tsmear'), ddb_hdr%tsmear))
2177 :
2178 : ! ------
2179 : ! Arrays
2180 : ! ------
2181 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'acell'), ddb_hdr%acell))
2182 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'ngfft'), ngfft))
2183 7543 : ddb_hdr%ngfft = zero
2184 1588 : do ii=1,3
2185 1588 : ddb_hdr%ngfft(ii) = ngfft(ii)
2186 : end do
2187 :
2188 : ! Basis
2189 : ! Note that internal value of mband might differ from file when merging ddbs.
2190 1985 : ABI_MALLOC(occ, (mband_file,nkpt_file,ddb_hdr%nsppol))
2191 1588 : ABI_MALLOC(nband, (nkpt_file,ddb_hdr%nsppol))
2192 :
2193 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'spinat'), ddb_hdr%spinat))
2194 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'occupations'), occ))
2195 1191 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, 'reduced_coordinates_of_kpoints'), ddb_hdr%kpt, count=[3,nkpt_file])
2196 397 : NCF_CHECK(ncerr)
2197 794 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, 'kpoints_weights'), ddb_hdr%wtk, count=[nkpt_file])
2198 397 : NCF_CHECK(ncerr)
2199 397 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, 'number_of_bands'), nband))
2200 :
2201 :
2202 : ! Flatten some arrays for internal storage
2203 227390 : ddb_hdr%occ = zero
2204 397 : mband_ = min(mband_file, ddb_hdr%mband)
2205 397 : nkpt_ = min(nkpt_file, ddb_hdr%nkpt)
2206 397 : ii = 0
2207 813 : do isppol=1,ddb_hdr%nsppol
2208 34967 : do ikpt=1,nkpt_
2209 257032 : do iband=1,mband_
2210 222462 : ii = ii + 1
2211 256616 : ddb_hdr%occ(ii) = occ(iband,ikpt,isppol)
2212 : end do
2213 : end do
2214 : end do
2215 35633 : ddb_hdr%nband = zero
2216 397 : ii = 0
2217 813 : do isppol=1,ddb_hdr%nsppol
2218 34967 : do ikpt=1,nkpt_
2219 34154 : ii = ii + 1
2220 34570 : ddb_hdr%nband(ii) = nband(ikpt,isppol)
2221 : end do
2222 : end do
2223 397 : ABI_FREE(occ)
2224 397 : ABI_FREE(nband)
2225 :
2226 : ! ----------------------
2227 : ! Initialize the crystal
2228 : ! ----------------------
2229 397 : ncid_crystal = nctk_idgroup(ncid, 'Crystal')
2230 397 : call ddb_hdr%crystal%ncread(ncid_crystal)
2231 :
2232 : ! Copy symmetries
2233 10064 : do isym=1,ddb_hdr%nsym
2234 9667 : ddb_hdr%symafm(isym) = ddb_hdr%crystal%symafm(isym)
2235 38668 : ddb_hdr%tnons(:,isym) = ddb_hdr%crystal%tnons(:,isym)
2236 126068 : ddb_hdr%symrel(:,:,isym) = ddb_hdr%crystal%symrel(:,:,isym)
2237 : end do
2238 :
2239 : ! Copy geometry
2240 1588 : do ii=1,3
2241 5161 : ddb_hdr%rprim(:,ii) = ddb_hdr%crystal%rprimd(:,ii) / ddb_hdr%acell(ii)
2242 : end do
2243 :
2244 3906 : ddb_hdr%xred = ddb_hdr%crystal%xred
2245 1572 : ddb_hdr%typat = ddb_hdr%crystal%typat
2246 1333 : ddb_hdr%znucl = ddb_hdr%crystal%znucl
2247 1333 : ddb_hdr%zion = ddb_hdr%crystal%zion
2248 1333 : ddb_hdr%amu = ddb_hdr%crystal%amu
2249 :
2250 : end if
2251 :
2252 : ! ======================
2253 : ! Master broadcasts data
2254 : ! ======================
2255 397 : call ddb_hdr%bcast(comm_)
2256 :
2257 563 : end subroutine ddb_hdr_open_read_nc
2258 : !!***
2259 :
2260 : !----------------------------------------------------------------------
2261 :
2262 : !!****f* m_ddb_hdr/ddb_hdr_close
2263 : !! NAME
2264 : !! ddb_hdr_close
2265 : !!
2266 : !! FUNCTION
2267 : !! Close the open DDB file in either txt or nc format
2268 : !! using the unit or ncid stored in the ddb_hdr object.
2269 : !!
2270 : !! INPUTS
2271 : !! comm=MPI communicator.
2272 : !!
2273 : !! SOURCE
2274 :
2275 1244 : subroutine ddb_hdr_close(ddb_hdr, comm)
2276 : !Arguments ------------------------------------
2277 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
2278 : integer,intent(in),optional :: comm
2279 :
2280 : !Local variables -------------------------
2281 : !scalars
2282 : integer,parameter :: master=0
2283 : integer :: ncerr, ierr
2284 : ! ************************************************************************
2285 :
2286 1244 : if (present(comm)) then
2287 0 : if (xmpi_comm_rank(comm) /= master) return
2288 : end if
2289 :
2290 1244 : if (ddb_hdr%has_open_file_txt) then
2291 544 : close(ddb_hdr%unddb, iostat=ierr)
2292 544 : ddb_hdr%has_open_file_txt = .false.
2293 : end if
2294 1244 : if (ddb_hdr%has_open_file_nc) then
2295 680 : ncerr = nf90_close(ddb_hdr%ncid)
2296 680 : ddb_hdr%has_open_file_nc = .false.
2297 : end if
2298 :
2299 : end subroutine ddb_hdr_close
2300 : !!***
2301 :
2302 : !----------------------------------------------------------------------
2303 :
2304 : !!****f* m_ddb_hdr/ddb_hdr_compare
2305 : !! NAME
2306 : !! ddb_hdr_compare
2307 : !!
2308 : !! FUNCTION
2309 : !! Compare two DDB headers and raise error if they differ.
2310 : !! Also, complete psps information if one has more info than the other.
2311 : !!
2312 : !! INPUTS
2313 : !! ddb_hdr2: ddb header to compare with.
2314 : !!
2315 : !! OUTPUT
2316 : !!
2317 : !! SOURCE
2318 :
2319 243 : subroutine ddb_hdr_compare(ddb_hdr1, ddb_hdr2)
2320 :
2321 : !Arguments ------------------------------------
2322 : class(ddb_hdr_type),intent(inout) :: ddb_hdr1
2323 : type(ddb_hdr_type),intent(inout) :: ddb_hdr2
2324 :
2325 : !Local variables -------------------------
2326 : integer :: bantot,ii,ij,isym
2327 : integer :: fullinit,fullinit2
2328 : !integer :: itypat
2329 : !real(dp) :: ekbcm8,ekbcmp
2330 : real(dp),parameter :: tol=1.0d-6 ! Limited by the precision of text DDB
2331 : character(len=500) :: msg
2332 : ! ************************************************************************
2333 :
2334 243 : fullinit = ddb_hdr1%with_psps * ddb_hdr1%with_dfpt_vars
2335 243 : fullinit2 = ddb_hdr2%with_psps * ddb_hdr2%with_dfpt_vars
2336 :
2337 : ! natom
2338 243 : call chki8(ddb_hdr1%natom,ddb_hdr2%natom,' natom')
2339 :
2340 : ! nkpt
2341 243 : if(fullinit/=0)then
2342 229 : if(ddb_hdr1%nkpt/=2*ddb_hdr2%nkpt .and. 2*ddb_hdr1%nkpt/=ddb_hdr2%nkpt)then
2343 : ! GA: I've disabled this check, but I might want to put it back...
2344 : !call chki8(nkpt,nkpt8,' nkpt')
2345 : else
2346 6 : write(std_out,*)' compar8 : assume that one of the DDB to be',' merged use Time-Reversal to'
2347 6 : write(std_out,*)' decrease the number of k-points'
2348 : end if
2349 : end if
2350 :
2351 : ! nband
2352 : ! There can be the case of perturbation at Gamma, that
2353 : ! only need half of the number of k points.
2354 229 : if(fullinit==0 .or. ddb_hdr2%nkpt==2*ddb_hdr1%nkpt)then
2355 19 : bantot=0
2356 1392 : do ii=1,ddb_hdr2%nkpt
2357 1392 : bantot=bantot+ddb_hdr1%nband(ii)
2358 : end do
2359 : else
2360 : bantot=0
2361 19445 : do ii=1,ddb_hdr1%nkpt
2362 : !if(ddb_hdr1%nkpt==ddb_hdr2%nkpt)then
2363 : ! call chki8(ddb_hdr1%nband(ii),ddb_hdr2%nband(ii),' nband')
2364 : !end if
2365 19445 : bantot=bantot+ddb_hdr1%nband(ii)
2366 : end do
2367 : end if
2368 :
2369 : ! nsppol
2370 243 : call chki8(ddb_hdr1%nsppol,ddb_hdr2%nsppol,'nsppol')
2371 :
2372 : ! nsym
2373 : !if(ddb_hdr1%nsym/=1 .and. ddb_hdr2%nsym/=1)then
2374 : ! call chki8(ddb_hdr1%nsym,ddb_hdr2%nsym,' nsym')
2375 : !end if
2376 :
2377 : ! ntypat
2378 243 : call chki8(ddb_hdr1%ntypat,ddb_hdr2%ntypat,'ntypat')
2379 :
2380 : ! acell
2381 972 : do ii=1,3
2382 972 : call chkr8(ddb_hdr1%acell(ii),ddb_hdr2%acell(ii),' acell',tol)
2383 : end do
2384 :
2385 : ! amu
2386 641 : do ii=1,ddb_hdr1%ntypat
2387 641 : call chkr8(ddb_hdr1%amu(ii),ddb_hdr2%amu(ii),' amu',tol)
2388 : end do
2389 :
2390 : ! ecut
2391 243 : call chkr8(ddb_hdr1%ecut,ddb_hdr2%ecut,' ecut',tol)
2392 :
2393 : ! pawecutdg (PAW only)
2394 243 : if (ddb_hdr1%usepaw==1) then
2395 2 : call chkr8(ddb_hdr1%pawecutdg,ddb_hdr2%pawecutdg,' ecut',tol)
2396 : end if
2397 :
2398 : ! iscf
2399 243 : if(fullinit/=0)then
2400 229 : call chki8(ddb_hdr1%iscf,ddb_hdr2%iscf,' iscf')
2401 : end if
2402 :
2403 : ! ixc
2404 243 : call chki8(ddb_hdr1%ixc,ddb_hdr2%ixc,' ixc')
2405 :
2406 : ! kpt and 14. kptnrm
2407 243 : if(ddb_hdr2%nkpt == ddb_hdr1%nkpt .and. fullinit/=0)then
2408 17264 : do ij=1,ddb_hdr1%nkpt
2409 68543 : do ii=1,3
2410 68372 : call chkr8(ddb_hdr1%kpt(ii,ij)/ddb_hdr1%kptnrm,ddb_hdr2%kpt(ii,ij)/ddb_hdr2%kptnrm,' kpt',tol)
2411 : end do
2412 : end do
2413 : end if
2414 :
2415 : ! ngfft
2416 : ! Allow for different ngfft
2417 : ! E.g. DFPT calculation might use a different ngfft from GS WFK.
2418 972 : do ii=1,3
2419 729 : if (ddb_hdr1%ngfft(ii) == ddb_hdr2%ngfft(ii)) cycle
2420 : write(msg,'(3a,i10,3a,i10,a)') &
2421 0 : 'Comparing integers for variable ngfft.',ch10,&
2422 0 : 'Value from input DDB is',ddb_hdr1%ngfft(ii),' and',ch10,&
2423 0 : 'from transfer DDB is',ddb_hdr2%ngfft(ii),'.'
2424 972 : ABI_WARNING(msg)
2425 : end do
2426 :
2427 : ! occ
2428 243 : if(fullinit/=0 .and. ddb_hdr1%nkpt==ddb_hdr2%nkpt)then
2429 151011 : do ii=1,bantot
2430 151011 : call chkr8(ddb_hdr1%occ(ii),ddb_hdr2%occ(ii),' occ',tol)
2431 : end do
2432 : end if
2433 :
2434 : ! rprim
2435 972 : do ii=1,3
2436 3159 : do ij=1,3
2437 2916 : call chkr8(ddb_hdr1%rprim(ii,ij),ddb_hdr2%rprim(ii,ij),' rprim',tol)
2438 : end do
2439 : end do
2440 :
2441 : ! dfpt_sciss
2442 243 : if(fullinit/=0)then
2443 229 : call chkr8(ddb_hdr1%dfpt_sciss,ddb_hdr2%dfpt_sciss,' dfpt_sciss',tol)
2444 : end if
2445 :
2446 : ! symrel and tnons
2447 243 : if(ddb_hdr1%nsym==ddb_hdr2%nsym)then
2448 7363 : do isym=1,ddb_hdr1%nsym
2449 28762 : do ii=1,3
2450 21399 : call chkr8(ddb_hdr1%tnons(ii,isym),ddb_hdr2%tnons(ii,isym),' tnons',tol)
2451 92729 : do ij=1,3
2452 85596 : call chki8(ddb_hdr1%symrel(ii,ij,isym),ddb_hdr2%symrel(ii,ij,isym),'symrel')
2453 : end do
2454 : end do
2455 : end do
2456 : end if
2457 :
2458 : ! typat
2459 641 : do ii=1,ddb_hdr1%ntypat
2460 641 : call chki8(ddb_hdr1%typat(ii),ddb_hdr2%typat(ii),' typat')
2461 : end do
2462 :
2463 : ! wtk
2464 243 : if(ddb_hdr1%nkpt==ddb_hdr2%nkpt .and. fullinit/=0)then
2465 17264 : do ii=1,ddb_hdr1%nkpt
2466 17264 : call chkr8(ddb_hdr1%wtk(ii),ddb_hdr2%wtk(ii),' wtk',tol)
2467 : end do
2468 : end if
2469 :
2470 : ! xred
2471 800 : do ij=1,ddb_hdr1%natom
2472 2471 : do ii=1,3
2473 2228 : call chkr8(ddb_hdr1%xred(ii,ij),ddb_hdr2%xred(ii,ij),' xred',tol)
2474 : end do
2475 : end do
2476 : ! zion
2477 641 : do ii=1,ddb_hdr1%ntypat
2478 641 : call chkr8(ddb_hdr1%zion(ii),ddb_hdr2%zion(ii),' zion',tol)
2479 : end do
2480 :
2481 : ! Now compare the NC pseudopotential information
2482 : ! GA: I disabled this comparison because the ekb
2483 : ! are not properly written to the DDB in text format.
2484 : ! See comment in psddb8
2485 : !if (ddb_hdr1%usepaw==0) then
2486 : ! if(ddb_hdr1%psps%dimekb/=0 .and. fullinit/=0 .and. fullinit2/=0 )then
2487 : ! do ii=1,dimekb
2488 : ! do itypat=1,ntypat
2489 : ! ekbcmp=ekb(ii,itypat)
2490 : ! ekbcm8=ekb8(ii,itypat)
2491 : ! !call chkr8(ekbcmp,ekbcm8,' ekb',tol)
2492 : ! end do
2493 : ! end do
2494 : !end if
2495 :
2496 : !Now compare several PAW dataset information
2497 : ! if (ddb_hdr1%usepaw==1) then
2498 : ! if (fullinit/=0 .and. fullinit2/=0) then
2499 : ! do itypat=1,ddb_hdr1%ntypat
2500 : ! call chki8(ddb_hdr1%pawtab(itypat)%basis_size,ddb_hdr2%pawtab(itypat)%basis_size,'bas_sz')
2501 : ! call chki8(ddb_hdr1%pawtab(itypat)%lmn_size,ddb_hdr2%pawtab(itypat)%lmn_size,'lmn_sz')
2502 : ! call chki8(ddb_hdr1%pawtab(itypat)%lmn2_size,ddb_hdr2%pawtab(itypat)%lmn2_size,'lmn2sz')
2503 : ! call chkr8(ddb_hdr1%pawtab(itypat)%rpaw,ddb_hdr2%pawtab(itypat)%rpaw,' rpaw',tol3)
2504 : ! call chkr8(ddb_hdr1%pawtab(itypat)%rshp,ddb_hdr2%pawtab(itypat)%rshp,'rshape',tol3)
2505 : ! call chki8(ddb_hdr1%pawtab(itypat)%shape_type,ddb_hdr2%pawtab(itypat)%shape_type,'shp_tp')
2506 : ! if (ddb_hdr1%pawtab(itypat)%lmn2_size>0) then
2507 : ! do ii=1,ddb_hdr1%pawtab(itypat)%lmn2_size
2508 : ! call chkr8(ddb_hdr1%pawtab(itypat)%dij0(ii),ddb_hdr2%pawtab(itypat)%dij0(ii),' dij0',tol)
2509 : ! end do
2510 : ! end if
2511 : ! end do
2512 : ! end if
2513 : ! end if
2514 :
2515 :
2516 :
2517 : ! Should also compare indlmn and pspso ... but suppose that
2518 : ! the checking of ekb is enough for the psps.
2519 : ! Should also compare many other variables ...
2520 :
2521 :
2522 243 : end subroutine ddb_hdr_compare
2523 : !!***
2524 :
2525 : !----------------------------------------------------------------------
2526 :
2527 : !!****f* m_ddb_hdr/ddb_hdr_copy_missing_variables
2528 : !! NAME
2529 : !! ddb_hdr_copy_missing_variables
2530 : !!
2531 : !! FUNCTION
2532 : !! Copy missing variables from ddb_hdr2 to ddb_hdr1
2533 : !!
2534 : !! INPUTS
2535 : !!
2536 : !! OUTPUT
2537 : !!
2538 : !! SOURCE
2539 :
2540 243 : subroutine ddb_hdr_copy_missing_variables(ddb_hdr1, ddb_hdr2)
2541 :
2542 : !Arguments ------------------------------------
2543 : class(ddb_hdr_type),intent(inout) :: ddb_hdr1
2544 : type(ddb_hdr_type),intent(inout) :: ddb_hdr2
2545 :
2546 : !Local variables -------------------------
2547 : integer :: fullinit, fullinit2
2548 : integer :: bantot,ii,ij,itypat
2549 : ! ************************************************************************
2550 :
2551 243 : fullinit = ddb_hdr1%with_psps * ddb_hdr1%with_dfpt_vars
2552 243 : fullinit2 = ddb_hdr2%with_psps * ddb_hdr2%with_dfpt_vars
2553 :
2554 : ! nkpt
2555 243 : if(fullinit==0)then
2556 14 : ddb_hdr1%nkpt = ddb_hdr2%nkpt
2557 : end if
2558 :
2559 : ! occopt
2560 : ! Take here the most favorable case.
2561 243 : if(ddb_hdr2%occopt==0) ddb_hdr1%occopt=0
2562 :
2563 : ! nband
2564 : ! There can be the case of perturbation at Gamma, that
2565 : ! only need half of the number of k points.
2566 243 : if(fullinit==0 .or. ddb_hdr2%nkpt==2*ddb_hdr1%nkpt)then
2567 19 : bantot=0
2568 1392 : do ii=1,ddb_hdr2%nkpt
2569 1373 : ddb_hdr1%nband(ii) = ddb_hdr2%nband(ii)
2570 1392 : bantot=bantot+ddb_hdr1%nband(ii)
2571 : end do
2572 : else
2573 : bantot=0
2574 19445 : do ii=1,ddb_hdr1%nkpt
2575 19445 : bantot=bantot+ddb_hdr1%nband(ii)
2576 : end do
2577 : end if
2578 :
2579 : ! iscf
2580 243 : if(fullinit==0)then
2581 14 : ddb_hdr1%iscf = ddb_hdr2%iscf
2582 : end if
2583 :
2584 : ! kpt and kptnrm
2585 243 : if(ddb_hdr2%nkpt == 2*ddb_hdr1%nkpt .or. fullinit==0) then
2586 :
2587 : ! Copy the largest number of k points in the right place
2588 1392 : do ij=1,ddb_hdr2%nkpt
2589 5511 : do ii=1,3
2590 5492 : ddb_hdr1%kpt(ii,ij)=ddb_hdr2%kpt(ii,ij)
2591 : end do
2592 : end do
2593 :
2594 19 : ddb_hdr1%kptnrm=ddb_hdr2%kptnrm
2595 :
2596 : end if
2597 :
2598 : ! occ
2599 243 : if (fullinit==0 .or. ddb_hdr2%nkpt==2*ddb_hdr1%nkpt) then
2600 9151 : do ii=1,bantot
2601 9375 : ddb_hdr1%occ(ii)=ddb_hdr2%occ(ii)
2602 : end do
2603 : end if
2604 :
2605 : ! dfpt_sciss
2606 243 : if(fullinit==0)then
2607 14 : ddb_hdr1%dfpt_sciss = ddb_hdr2%dfpt_sciss
2608 : end if
2609 :
2610 : ! symrel
2611 243 : if(ddb_hdr1%nsym/=ddb_hdr2%nsym .and. ddb_hdr2%nsym/=1)then
2612 : ! GA: It is suspicious to assume the array are allocated with the same size
2613 : ! but this is only used in merge_ddb
2614 0 : ddb_hdr1%nsym = ddb_hdr2%nsym
2615 0 : ddb_hdr1%symrel(:,:,1:ddb_hdr2%nsym)=ddb_hdr2%symrel(:,:,1:ddb_hdr2%nsym)
2616 0 : ddb_hdr1%tnons(:,1:ddb_hdr2%nsym)=ddb_hdr2%tnons(:,1:ddb_hdr2%nsym)
2617 : end if
2618 :
2619 : ! tolwfr
2620 : !Take the less converged value...
2621 : !ddb_hdr1%tolwfr=max(ddb_hdr1%tolwfr,ddb_hdr2%tolwfr)
2622 :
2623 : ! wtk
2624 243 : if(ddb_hdr2%nkpt==2*ddb_hdr1%nkpt .or. fullinit==0)then
2625 1392 : do ii=1,ddb_hdr2%nkpt
2626 1616 : ddb_hdr1%wtk(ii)=ddb_hdr2%wtk(ii)
2627 : end do
2628 : end if
2629 :
2630 : ! Copy nkpt in the case of the use of the time-reversal symmetry
2631 243 : if(2*ddb_hdr1%nkpt==ddb_hdr2%nkpt)then
2632 5 : ddb_hdr1%nkpt=ddb_hdr2%nkpt
2633 : end if
2634 :
2635 : ! PAW dataset information
2636 243 : if (ddb_hdr1%usepaw==1) then
2637 2 : if (fullinit==0 .and. fullinit2/=0) then
2638 3 : do itypat=1,ddb_hdr1%ntypat
2639 2 : ddb_hdr1%pawtab(itypat)%basis_size = ddb_hdr2%pawtab(itypat)%basis_size
2640 2 : ddb_hdr1%pawtab(itypat)%lmn_size = ddb_hdr2%pawtab(itypat)%lmn_size
2641 2 : ddb_hdr1%pawtab(itypat)%rpaw = ddb_hdr2%pawtab(itypat)%rpaw
2642 2 : ddb_hdr1%pawtab(itypat)%rshp = ddb_hdr2%pawtab(itypat)%rshp
2643 2 : ddb_hdr1%pawtab(itypat)%shape_type = ddb_hdr2%pawtab(itypat)%shape_type
2644 2 : if (ddb_hdr2%pawtab(itypat)%lmn2_size>0) then
2645 2 : if (ddb_hdr1%pawtab(itypat)%lmn2_size==0) then
2646 0 : ABI_MALLOC(ddb_hdr1%pawtab(itypat)%dij0,(ddb_hdr2%pawtab(itypat)%lmn2_size))
2647 : end if
2648 74 : do ii=1,ddb_hdr2%pawtab(itypat)%lmn2_size
2649 74 : ddb_hdr1%pawtab(itypat)%dij0(ii) = ddb_hdr2%pawtab(itypat)%dij0(ii)
2650 : end do
2651 : end if
2652 3 : ddb_hdr1%pawtab(itypat)%lmn2_size = ddb_hdr2%pawtab(itypat)%lmn2_size
2653 : end do
2654 : end if
2655 : end if
2656 243 : end subroutine ddb_hdr_copy_missing_variables
2657 : !!***
2658 :
2659 : !----------------------------------------------------------------------
2660 :
2661 : !!****f* m_ddb_hdr/ddb_hdr_copy_psps_info
2662 : !! NAME
2663 : !! ddb_hdr_compare
2664 : !!
2665 : !! FUNCTION
2666 : !! Copy the information on pseudopotentials from one header to another.
2667 : !!
2668 : !! INPUTS
2669 : !! ddb_hdr1: the ddb header into which info is copied
2670 : !! ddb_hdr2: the ddb header from which info is copied
2671 : !!
2672 : !! OUTPUT
2673 : !!
2674 : !! SOURCE
2675 :
2676 47 : subroutine ddb_hdr_copy_psps_from(ddb_hdr1, ddb_hdr2)
2677 :
2678 : !Arguments ------------------------------------
2679 : class(ddb_hdr_type),intent(inout) :: ddb_hdr1
2680 : type(ddb_hdr_type),intent(inout) :: ddb_hdr2
2681 :
2682 : !Local variables -------------------------
2683 : integer :: itypat
2684 :
2685 : ! ************************************************************************
2686 :
2687 47 : if (ddb_hdr2%with_psps > 0) then
2688 46 : ddb_hdr1%with_psps = ddb_hdr2%with_psps
2689 :
2690 46 : call psps_free(ddb_hdr1%psps)
2691 46 : call psps_copy(ddb_hdr2%psps, ddb_hdr1%psps)
2692 :
2693 : end if
2694 :
2695 : ! TODO: Need a function pawtab_copy
2696 47 : if (ddb_hdr1%psps%usepaw > 0) then
2697 3 : do itypat=1,ddb_hdr1%mtypat
2698 2 : ddb_hdr1%pawtab(itypat)%basis_size = ddb_hdr2%pawtab(itypat)%basis_size
2699 2 : ddb_hdr1%pawtab(itypat)%lmn_size = ddb_hdr2%pawtab(itypat)%lmn_size
2700 2 : ddb_hdr1%pawtab(itypat)%lmn2_size = ddb_hdr2%pawtab(itypat)%lmn2_size
2701 2 : ddb_hdr1%pawtab(itypat)%shape_type = ddb_hdr2%pawtab(itypat)%shape_type
2702 2 : ddb_hdr1%pawtab(itypat)%rpaw = ddb_hdr2%pawtab(itypat)%rpaw
2703 2 : ddb_hdr1%pawtab(itypat)%rshp = ddb_hdr2%pawtab(itypat)%rshp
2704 3 : if (allocated(ddb_hdr2%pawtab(itypat)%dij0)) then
2705 2 : ABI_SFREE(ddb_hdr1%pawtab(itypat)%dij0)
2706 6 : ABI_MALLOC(ddb_hdr1%pawtab(itypat)%dij0, (ddb_hdr1%pawtab(itypat)%lmn2_size))
2707 74 : ddb_hdr1%pawtab(itypat)%dij0(:) = ddb_hdr2%pawtab(itypat)%dij0(:)
2708 : end if
2709 : end do
2710 : end if
2711 :
2712 47 : end subroutine ddb_hdr_copy_psps_from
2713 : !!***
2714 :
2715 : !----------------------------------------------------------------------
2716 :
2717 : !!****f* m_ddb_hdr/ddb_hdr_bcast_dim
2718 : !! NAME
2719 : !! ddb_hdr_bcast_dim
2720 : !!
2721 : !! FUNCTION
2722 : !! Master broadcasts the header dimensions.
2723 : !!
2724 : !! INPUTS
2725 : !! comm=MPI communicator.
2726 : !!
2727 : !! SOURCE
2728 :
2729 190 : subroutine ddb_hdr_bcast_dim(ddb_hdr, comm)
2730 :
2731 : !Arguments ------------------------------------
2732 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
2733 : integer, intent(in) :: comm
2734 :
2735 : !Local variables -------------------------
2736 : integer, parameter :: master=0
2737 : integer :: ierr
2738 : ! ************************************************************************
2739 :
2740 190 : if (xmpi_comm_size(comm) == 1) return
2741 :
2742 : DBG_ENTER("COLL")
2743 :
2744 : ! Integers
2745 24 : call xmpi_bcast(ddb_hdr%ddb_version, master, comm, ierr)
2746 :
2747 24 : call xmpi_bcast(ddb_hdr%matom, master, comm, ierr)
2748 24 : call xmpi_bcast(ddb_hdr%mband, master, comm, ierr)
2749 24 : call xmpi_bcast(ddb_hdr%mkpt, master, comm, ierr)
2750 24 : call xmpi_bcast(ddb_hdr%msym, master, comm, ierr)
2751 24 : call xmpi_bcast(ddb_hdr%mtypat, master, comm, ierr)
2752 :
2753 24 : call xmpi_bcast(ddb_hdr%natom, master, comm, ierr)
2754 24 : call xmpi_bcast(ddb_hdr%nkpt, master, comm, ierr)
2755 24 : call xmpi_bcast(ddb_hdr%nsym, master, comm, ierr)
2756 24 : call xmpi_bcast(ddb_hdr%ntypat, master, comm, ierr)
2757 :
2758 24 : call xmpi_bcast(ddb_hdr%nsppol, master, comm, ierr)
2759 24 : call xmpi_bcast(ddb_hdr%nspinor, master, comm, ierr)
2760 24 : call xmpi_bcast(ddb_hdr%nspden, master, comm, ierr)
2761 :
2762 24 : call xmpi_bcast(ddb_hdr%nblok, master, comm, ierr)
2763 24 : call xmpi_bcast(ddb_hdr%has_d0E_xx, master, comm, ierr)
2764 24 : call xmpi_bcast(ddb_hdr%has_d1E_xx, master, comm, ierr)
2765 24 : call xmpi_bcast(ddb_hdr%has_d2E_xx, master, comm, ierr)
2766 24 : call xmpi_bcast(ddb_hdr%has_d3E_xx, master, comm, ierr)
2767 24 : call xmpi_bcast(ddb_hdr%has_d3E_lw, master, comm, ierr)
2768 24 : call xmpi_bcast(ddb_hdr%has_d2eig, master, comm, ierr)
2769 :
2770 24 : call xmpi_bcast(ddb_hdr%mpert, master, comm, ierr)
2771 24 : call xmpi_bcast(ddb_hdr%msize, master, comm, ierr)
2772 :
2773 24 : call xmpi_bcast(ddb_hdr%with_psps, master, comm, ierr)
2774 24 : call xmpi_bcast(ddb_hdr%with_dfpt_vars, master, comm, ierr)
2775 :
2776 24 : call xmpi_bcast(ddb_hdr%usepaw, master, comm, ierr)
2777 24 : call xmpi_bcast(ddb_hdr%psps%dimekb, master, comm, ierr)
2778 24 : call xmpi_bcast(ddb_hdr%psps%ntypat, master, comm, ierr)
2779 24 : call xmpi_bcast(ddb_hdr%psps%lmnmax, master, comm, ierr)
2780 24 : call xmpi_bcast(ddb_hdr%psps%usepaw, master, comm, ierr)
2781 24 : call xmpi_bcast(ddb_hdr%psps%useylm, master, comm, ierr)
2782 :
2783 : DBG_EXIT("COLL")
2784 :
2785 : end subroutine ddb_hdr_bcast_dim
2786 : !!***
2787 :
2788 : !----------------------------------------------------------------------
2789 :
2790 : !!****f* m_ddb_hdr/ddb_hdr_bcast
2791 : !! NAME
2792 : !! ddb_hdr_bcast
2793 : !!
2794 : !! FUNCTION
2795 : !! Master broadcasts the data and others allocate their arrays
2796 : !!
2797 : !! INPUTS
2798 : !! comm=MPI communicator.
2799 : !!
2800 : !! SOURCE
2801 :
2802 803 : subroutine ddb_hdr_bcast(ddb_hdr, comm)
2803 :
2804 : !Arguments ------------------------------------
2805 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
2806 : integer, intent(in) :: comm
2807 :
2808 : !Local variables -------------------------
2809 : integer, parameter :: master=0
2810 : integer :: ierr, ii,nn
2811 : ! ************************************************************************
2812 :
2813 803 : if (xmpi_comm_size(comm) == 1) return
2814 :
2815 : DBG_ENTER("COLL")
2816 :
2817 24 : call ddb_hdr%bcast_dim(comm)
2818 :
2819 : ! Allocate arrays on the other nodes.
2820 24 : if (xmpi_comm_rank(comm) /= master) then
2821 :
2822 16 : call ddb_hdr%free()
2823 16 : call ddb_hdr%malloc()
2824 :
2825 64 : ABI_MALLOC(ddb_hdr%psps%indlmn,(6,ddb_hdr%psps%lmnmax,ddb_hdr%mtypat))
2826 48 : ABI_MALLOC(ddb_hdr%psps%pspso,(ddb_hdr%mtypat))
2827 64 : ABI_MALLOC(ddb_hdr%psps%ekb,(ddb_hdr%psps%dimekb,ddb_hdr%mtypat * (1 - ddb_hdr%psps%usepaw)))
2828 :
2829 64 : ABI_MALLOC(ddb_hdr%pawtab,(ddb_hdr%psps%ntypat*ddb_hdr%psps%usepaw))
2830 16 : call pawtab_nullify(ddb_hdr%pawtab)
2831 :
2832 16 : ddb_hdr%dscrpt = ''
2833 :
2834 : end if
2835 :
2836 : ! Floats
2837 24 : call xmpi_bcast(ddb_hdr%dilatmx, master, comm, ierr)
2838 24 : call xmpi_bcast(ddb_hdr%ecut, master, comm, ierr)
2839 24 : call xmpi_bcast(ddb_hdr%kptnrm, master, comm, ierr)
2840 24 : call xmpi_bcast(ddb_hdr%pawecutdg, master, comm, ierr)
2841 24 : call xmpi_bcast(ddb_hdr%dfpt_sciss, master, comm, ierr)
2842 24 : call xmpi_bcast(ddb_hdr%tolwfr, master, comm, ierr)
2843 24 : call xmpi_bcast(ddb_hdr%tphysel, master, comm, ierr)
2844 24 : call xmpi_bcast(ddb_hdr%tsmear, master, comm, ierr)
2845 :
2846 : ! Arrays
2847 24 : call xmpi_bcast(ddb_hdr%ngfft, master, comm, ierr)
2848 24 : call xmpi_bcast(ddb_hdr%acell, master, comm, ierr)
2849 24 : call xmpi_bcast(ddb_hdr%rprim, master, comm, ierr)
2850 :
2851 24 : call xmpi_bcast(ddb_hdr%nband, master, comm, ierr)
2852 24 : call xmpi_bcast(ddb_hdr%symafm, master, comm, ierr)
2853 24 : call xmpi_bcast(ddb_hdr%symrel, master, comm, ierr)
2854 24 : call xmpi_bcast(ddb_hdr%typat, master, comm, ierr)
2855 :
2856 24 : call xmpi_bcast(ddb_hdr%amu, master, comm, ierr)
2857 24 : call xmpi_bcast(ddb_hdr%kpt, master, comm, ierr)
2858 24 : call xmpi_bcast(ddb_hdr%occ, master, comm, ierr)
2859 24 : call xmpi_bcast(ddb_hdr%spinat, master, comm, ierr)
2860 24 : call xmpi_bcast(ddb_hdr%tnons, master, comm, ierr)
2861 24 : call xmpi_bcast(ddb_hdr%wtk, master, comm, ierr)
2862 24 : call xmpi_bcast(ddb_hdr%xred, master, comm, ierr)
2863 24 : call xmpi_bcast(ddb_hdr%zion, master, comm, ierr)
2864 24 : call xmpi_bcast(ddb_hdr%znucl, master, comm, ierr)
2865 :
2866 24 : call xmpi_bcast(ddb_hdr%psps%indlmn, master, comm, ierr)
2867 24 : call xmpi_bcast(ddb_hdr%psps%pspso, master, comm, ierr)
2868 :
2869 24 : if (ddb_hdr%psps%dimekb > 0 .and. ddb_hdr%psps%usepaw==0) then
2870 12 : call xmpi_bcast(ddb_hdr%psps%ekb, master, comm, ierr)
2871 : end if
2872 :
2873 24 : nn = ddb_hdr%psps%ntypat * ddb_hdr%psps%usepaw
2874 24 : if (nn.gt.0) then
2875 36 : do ii=1,nn
2876 36 : call pawtab_bcast(ddb_hdr%pawtab(ii), comm)
2877 : end do
2878 : end if
2879 :
2880 : ! Strings
2881 : !call xmpi_bcast(ddb_hdr%dscrpt, master, comm, ierr)
2882 :
2883 24 : call ddb_hdr%crystal%bcast(comm)
2884 :
2885 : end subroutine ddb_hdr_bcast
2886 : !!***
2887 :
2888 : !----------------------------------------------------------------------
2889 :
2890 : !!****f* m_ddb_hdr/ddb_hdr_get_iomode
2891 : !! NAME
2892 : !! ddb_hdr_get_iomode
2893 : !!
2894 : !! FUNCTION
2895 : !! Decide whether to use text or netcdf IO
2896 : !!
2897 : !! INPUTS
2898 : !! filename=name of the file to examine.
2899 : !! io=specify whether we want to read (io=1) or write (io=2)
2900 : !!
2901 : !! OUTPUT
2902 : !! iomode=The io mode.
2903 : !! filenameout=The filename, possibly appended by '.nc'
2904 : !!
2905 : !! SOURCE
2906 :
2907 4583 : subroutine ddb_hdr_get_iomode(ddb_hdr, filenamein, io, iomode, filenameout)
2908 :
2909 : !Arguments ------------------------------------
2910 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
2911 : character(len=fnlen),intent(in) :: filenamein
2912 : integer,intent(in) :: io
2913 : integer,intent(out) :: iomode
2914 : character(len=fnlen),intent(out) :: filenameout
2915 :
2916 : !Local variables ------------------------------
2917 : character(len=fnlen) :: filenamenc
2918 : ! ************************************************************************
2919 :
2920 : ! Simple case: file name was given with netcdf extension
2921 : #ifdef HAVE_NETCDF
2922 4583 : if (endswith(filenamein, ".nc")) then
2923 5 : filenameout = filenamein
2924 5 : iomode = IO_MODE_ETSF
2925 4583 : return
2926 : else
2927 : ! In most case, this function is called without filename extension.
2928 4578 : filenamenc = nctk_ncify(filenamein)
2929 : end if
2930 : #else
2931 : if (endswith(filenamein, ".nc")) then
2932 : ABI_ERROR("NETCDF support required to read or write DDB.nc file.")
2933 : else
2934 : filenameout = filenamein
2935 : iomode = IO_MODE_FORTRAN
2936 : return
2937 : end if
2938 : #endif
2939 :
2940 : ! Reading
2941 4578 : if (io==1) then
2942 :
2943 : ! First check for the existence of netcdf file
2944 1169 : if (file_exists(filenamenc)) then
2945 559 : filenameout = filenamenc
2946 559 : iomode = IO_MODE_ETSF
2947 559 : return
2948 :
2949 : ! Then check for the existence of text file
2950 610 : else if (file_exists(filenamein)) then
2951 610 : filenameout = filenamein
2952 610 : iomode = IO_MODE_FORTRAN
2953 610 : return
2954 :
2955 : ! File is not available for reading.
2956 : else
2957 :
2958 : ! This doesnt always have to be an error
2959 : !ABI_ERROR(sjoin("Cannot find DDB file:", filenamein))
2960 0 : ABI_WARNING(sjoin("Cannot find DDB file:", filenamein))
2961 :
2962 : ! Rely on default value of iomode.
2963 0 : iomode = ddb_hdr%iomode
2964 0 : if (iomode==IO_MODE_ETSF) then
2965 0 : filenameout = filenamenc
2966 : else
2967 0 : filenameout = filenamein
2968 : end if
2969 0 : return
2970 : end if
2971 :
2972 : ! Writing
2973 3409 : else if (io==2) then
2974 :
2975 : ! Rely on ddb_hdr%iomode, which should have been initialized from dtset
2976 3409 : if (ddb_hdr%iomode==IO_MODE_ETSF) then
2977 3277 : filenameout = filenamenc
2978 3277 : iomode = IO_MODE_ETSF
2979 3277 : return
2980 132 : else if ((ddb_hdr%iomode==IO_MODE_FORTRAN).or.(ddb_hdr%iomode==IO_MODE_FORTRAN_MASTER).or.(ddb_hdr%iomode==IO_MODE_MPI)) then
2981 132 : filenameout = filenamein
2982 132 : iomode = IO_MODE_FORTRAN
2983 132 : return
2984 : else
2985 0 : ABI_ERROR("Unexpected value for iomode.")
2986 : end if
2987 : else
2988 0 : ABI_ERROR("Unexpected value for io.")
2989 : end if
2990 :
2991 : end subroutine ddb_hdr_get_iomode
2992 : !!***
2993 :
2994 : !----------------------------------------------------------------------
2995 :
2996 : !!****f* m_ddb_hdr/psddb8
2997 : !!
2998 : !! NAME
2999 : !! psddb8
3000 : !!
3001 : !! FUNCTION
3002 : !! Take care of the i/o of pseudopotentials for the
3003 : !! Derivative DataBase, and also the number of data blocks.
3004 : !!
3005 : !! INPUTS
3006 : !! choice=(1 => read), (2=> write)
3007 : !! dimekb=dimension of ekb (contains Kleimann-Bylander energies)
3008 : !! used only for norm-conserving pseudopotentials
3009 : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
3010 : !! =if useylm=0, max number of (l,n) comp. over all type of psps
3011 : !! nunit=unit number for the Derivative DataBase.
3012 : !! ntypat=number of atom types
3013 : !! pspso(ntypat)=For each type of psp, 1 if no spin-orbit component is taken
3014 : !! into account, 2 if a spin-orbit component is used
3015 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
3016 : !! useylm=governs the way the nonlocal operator is to be applied:
3017 : !! 1=using Ylm, 0=using Legendre polynomials
3018 : !!
3019 : !! OUTPUT
3020 : !! (see side effects)
3021 : !!
3022 : !! SIDE EFFECTS
3023 : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
3024 : !! or i=lmn (if useylm=1)
3025 : !! ekb(dimekb,ntypat)= (norm-conserving psps only) (Real) Kleinman-Bylander energies (hartree)
3026 : !! Presently the only characteristics of the psp
3027 : !! with_psps=0 if the ekb are not available, at input as well as at output
3028 : !! pawtab(ntypat*usepaw)= (PAW only) PAW datasets characteristics
3029 : !! Presently only pawtab%basis_size,pawtab%lmn_size,pawtab%shape_type
3030 : !! pawtab%rpaw,pawtab%rshp,pawtab%dij0 are used
3031 : !! nblok=number of blocks
3032 : !!
3033 : !! NOTES
3034 : !! Only executed by one processor
3035 : !!
3036 : !! SOURCE
3037 :
3038 536 : subroutine psddb8 (choice,dimekb,ekb,with_psps,indlmn,lmnmax,&
3039 536 : & nblok,ntypat,nunit,pawtab,pspso,usepaw,useylm)
3040 :
3041 : !Arguments -------------------------------
3042 : !scalars
3043 : integer,intent(in) :: choice,dimekb,lmnmax,ntypat,nunit,usepaw,useylm
3044 : integer,intent(inout) :: with_psps,nblok
3045 : !arrays
3046 : integer,intent(in) :: pspso(ntypat)
3047 : integer,intent(inout) :: indlmn(6,lmnmax,ntypat)
3048 : real(dp),intent(inout) :: ekb(dimekb,ntypat)
3049 : type(pawtab_type),intent(inout) :: pawtab(ntypat*usepaw)
3050 :
3051 : !Local variables -------------------------
3052 : !Set the version number
3053 : !scalars
3054 : integer,parameter :: vrsio8=100401,vrsio8_old=010929,vrsio8_old_old=990527
3055 : integer :: basis_size0,dimekb0,iekb,ii,ij,il,ilm,ilmn,iln,iln0,im,ios,iproj,iproj0,itypat,itypat0
3056 : integer :: jekb,jlmn,jln,lmnmax0,lmn_size0,lmn2_size0,lpsang,nekb,nproj,npsang,pspso0,shape_type0
3057 : integer :: usepaw0,vrspsp8
3058 : real(dp) :: rpaw0,rshape0
3059 : character(len=12) :: string
3060 : character(len=500) :: message
3061 : !arrays
3062 536 : integer,allocatable :: i1(:),i2(:),nprj(:),orbitals(:)
3063 536 : real(dp),allocatable :: dij0(:),ekb0(:,:)
3064 : ! *********************************************************************
3065 :
3066 : !Check the value of choice
3067 536 : if (choice<=0.or.choice>=3) then
3068 : write(message, '(a,a,a,i10,a)' )&
3069 0 : 'The permitted values for choice are 1 or 2.',ch10,&
3070 0 : 'The calling routine asks ',choice,'.'
3071 0 : ABI_BUG(message)
3072 : end if
3073 :
3074 : !==================================================================================
3075 : !First option: read psp characteristic from file =================================
3076 : !==================================================================================
3077 536 : if (choice==1) then
3078 :
3079 386 : read(nunit,*)
3080 386 : read(nunit, '(a12)' )string
3081 28808 : with_psps=1 ; if (lmnmax>0) indlmn(:,:,:)=0
3082 :
3083 : ! --------------------------------------------
3084 : ! ----- NEW FORMAT (NCPP+PAW) ---------------
3085 : ! --------------------------------------------
3086 386 : if (string==' Descriptio')then
3087 : ! ==============================
3088 : ! ======== Common data =========
3089 : ! ==============================
3090 386 : read(nunit, '(32x,i6)' )vrspsp8
3091 386 : if (vrspsp8==vrsio8_old.or.vrspsp8==vrsio8_old_old) then
3092 55 : usepaw0=0
3093 : ! this format statement is inconsistent with test v2[13]
3094 : ! JWZ 31 May 2024
3095 : !read(nunit, '(10x,i3,14x,i3,11x,i3)', iostat=ios)dimekb0,lmnmax0,usepaw0
3096 55 : read(nunit, '(10x,i3,14x,i3,14x,i3)', iostat=ios)dimekb0,lmnmax0,usepaw0
3097 55 : if(ios/=0)then
3098 0 : backspace(nunit)
3099 0 : read (nunit, '(10x,i3,14x,i3)' )dimekb0,lmnmax0
3100 0 : usepaw0=0
3101 : end if
3102 331 : else if (vrspsp8==vrsio8) then
3103 331 : read(nunit, '(10x,i3)') usepaw0
3104 331 : if (usepaw/=usepaw0) then
3105 0 : write(message, '(a,i1,a,i1,a)' )'usepaw is announced to be ',usepaw,' but read usepaw is ',usepaw0,' !'
3106 0 : ABI_ERROR(message)
3107 : end if
3108 331 : if (usepaw==0) then
3109 326 : read (nunit, '(10x,i3,14x,i3)' )dimekb0,lmnmax0
3110 : end if
3111 : end if
3112 :
3113 : ! ==============================
3114 : ! === Norm-conserving psps =====
3115 : ! ==============================
3116 386 : if (usepaw==0) then
3117 4832 : ekb(:,:)=zero
3118 1524 : ABI_MALLOC(ekb0,(dimekb,dimekb))
3119 1109 : do itypat=1,ntypat
3120 728 : read(nunit, '(13x,i4,9x,i3,8x,i4)' )itypat0,pspso0,nekb
3121 : ! Check the compatibility with the main code dimensioning
3122 728 : if(nekb>dimekb)then
3123 : write(message, '(a,i8,a,a,a,i3,a)' )&
3124 0 : ' ',nekb,' components of ekb are announced',ch10,'but dimekb=',dimekb,'.'
3125 0 : ABI_BUG(message)
3126 : end if
3127 728 : read(nunit,*)
3128 728 : ilmn=0;iproj0=0
3129 4456 : do iekb=1,nekb
3130 3347 : read(nunit, '(3i6,3x,8d15.7)' ) iln,lpsang,iproj,(ekb0(ii,iekb),ii=1,min(nekb,4))
3131 3347 : if(nekb>4)then
3132 5448 : do jekb=5,nekb,4
3133 5448 : read(nunit, '(21x,8d15.7)' )(ekb0(ii,iekb),ii=jekb,min(nekb,jekb+3))
3134 : end do
3135 : end if
3136 3347 : if (lpsang==0.and.iproj>iproj0) iproj0=iproj
3137 3347 : if (useylm==1) then
3138 0 : do im=-lpsang,lpsang
3139 0 : ilmn=ilmn+1
3140 0 : indlmn(1,ilmn,itypat)=lpsang
3141 0 : indlmn(2,ilmn,itypat)=im
3142 0 : indlmn(3,ilmn,itypat)=iproj
3143 0 : indlmn(4,ilmn,itypat)=lpsang**2+lpsang+1+im
3144 0 : indlmn(5,ilmn,itypat)=iln
3145 0 : indlmn(6,ilmn,itypat)=1
3146 0 : if (pspso0/=1.and.iln>(nekb-iproj0)/2) indlmn(6,ilmn,itypat)=2
3147 : end do
3148 : else
3149 3347 : ilmn=ilmn+1
3150 3347 : indlmn(1,ilmn,itypat)=lpsang
3151 3347 : indlmn(2,ilmn,itypat)=lpsang
3152 3347 : indlmn(3,ilmn,itypat)=iproj
3153 3347 : indlmn(4,ilmn,itypat)=lpsang**2+lpsang+1
3154 3347 : indlmn(5,ilmn,itypat)=iln
3155 3347 : indlmn(6,ilmn,itypat)=1
3156 3347 : if (pspso0/=1.and.iln>(nekb-iproj0)/2) indlmn(6,ilmn,itypat)=2
3157 : end if
3158 : ! For the time being, only diagonal ekb are treated in abinit v3
3159 4075 : ekb(iekb,itypat)=ekb0(iekb,iekb)
3160 : ! For non-diagonal ekb, one could use:
3161 : ! do jekb=iekb to nekb
3162 : ! ekb(jekb+iekb*(iekb-1)/2,itypat)=ekb0(jekb,iekb)
3163 : ! end do
3164 : end do
3165 : end do
3166 381 : ABI_FREE(ekb0)
3167 :
3168 : ! ==============================
3169 : ! ============ PAW =============
3170 : ! ==============================
3171 : else
3172 16 : do itypat=1,ntypat
3173 11 : read(nunit, '(12x,i4,12x,i3,12x,i5)' )itypat0,basis_size0,lmn_size0
3174 11 : lmn2_size0=lmn_size0*(lmn_size0+1)/2
3175 33 : ABI_MALLOC(orbitals,(basis_size0))
3176 11 : read(nunit, '(20x,50i2)' ) orbitals(1:basis_size0)
3177 11 : read(nunit, '(11x,f6.3,13x,i2,11x,f6.3)' ) rpaw0,shape_type0,rshape0
3178 11 : read(nunit,'(24x,i3)') nekb
3179 11 : read(nunit,*)
3180 33 : ABI_MALLOC(dij0,(nekb))
3181 33 : ABI_MALLOC(i1,(nekb))
3182 22 : ABI_MALLOC(i2,(nekb))
3183 52 : do ii=1,nekb,4
3184 214 : read(nunit,'(3x,4(1x,i4,1x,i4,1x,d12.5))') (i1(ij),i2(ij),dij0(ij),ij=ii,min(ii+3,nekb))
3185 : end do
3186 11 : if (lmn_size0>lmnmax) then
3187 : write(message, '(a,i5,3a,i5,a)' )&
3188 0 : 'max. value of ',lmnmax,' for lmn_size is announced',ch10,'but ',lmn_size0,' is read.'
3189 0 : ABI_BUG(message)
3190 : end if
3191 11 : if (allocated(pawtab(itypat)%dij0)) then
3192 0 : if (lmn_size0>pawtab(itypat)%lmn_size) then
3193 : write(message, '(a,i5,3a,i5,a)' )&
3194 0 : 'lmn_size=,',pawtab(itypat)%lmn_size,' is announced',ch10,'but ',lmn_size0,' is read.'
3195 0 : ABI_BUG(message)
3196 : end if
3197 : end if
3198 81 : ABI_MALLOC(nprj,(0:maxval(orbitals)))
3199 35 : ilmn=0;nprj=0
3200 59 : do iln=1,basis_size0
3201 48 : il=orbitals(iln)
3202 48 : nprj(il)=nprj(il)+1
3203 156 : do ilm=1,2*il+1
3204 108 : indlmn(1,ilmn+ilm,itypat)=il
3205 108 : indlmn(2,ilmn+ilm,itypat)=ilm-(il+1)
3206 108 : indlmn(3,ilmn+ilm,itypat)=nprj(il)
3207 108 : indlmn(4,ilmn+ilm,itypat)=il*il+ilm
3208 108 : indlmn(5,ilmn+ilm,itypat)=iln
3209 156 : indlmn(6,ilmn+ilm,itypat)=1
3210 : end do
3211 59 : ilmn=ilmn+2*il+1
3212 : end do
3213 11 : pawtab(itypat)%basis_size=basis_size0
3214 11 : pawtab(itypat)%lmn_size =lmn_size0
3215 11 : pawtab(itypat)%lmn2_size =lmn2_size0
3216 11 : pawtab(itypat)%shape_type=shape_type0
3217 11 : pawtab(itypat)%rpaw =rpaw0
3218 11 : pawtab(itypat)%rshp =rshape0
3219 11 : if (.not.allocated(pawtab(itypat)%dij0)) then
3220 33 : ABI_MALLOC(pawtab(itypat)%dij0,(lmn2_size0))
3221 : end if
3222 677 : pawtab(itypat)%dij0(1:lmn2_size0)=zero
3223 173 : do ii=1,nekb
3224 162 : ij=i1(ii)+i2(ii)*(i2(ii)-1)/2
3225 173 : pawtab(itypat)%dij0(ij)=dij0(ii)
3226 : end do
3227 11 : ABI_FREE(nprj)
3228 11 : ABI_FREE(orbitals)
3229 11 : ABI_FREE(dij0)
3230 11 : ABI_FREE(i1)
3231 16 : ABI_FREE(i2)
3232 : end do
3233 :
3234 : end if ! NCPP or PAW
3235 :
3236 : ! --------------------------------------------
3237 : ! ----- OLD FORMAT (NCPP only) --------------
3238 : ! --------------------------------------------
3239 0 : else if (string==' Description')then
3240 0 : if (usepaw==1) then
3241 0 : ABI_BUG("old DDB pspformat not compatible with PAW")
3242 : end if
3243 :
3244 0 : read (nunit, '(10x,i3,10x,i3)' )nproj,npsang
3245 0 : nekb=nproj*npsang
3246 : ! Check the compatibility with the main code dimensioning
3247 0 : if(nekb>dimekb)then
3248 : write(message, '(a,i8,a,a,a,i3,a)' )&
3249 0 : ' ',nekb,' components of ekb are announced',ch10,'but the maximum is dimekb=',dimekb,'.'
3250 0 : ABI_BUG(message)
3251 : end if
3252 0 : if(useylm/=0)then
3253 0 : ABI_BUG('useylm must be 0 !')
3254 : end if
3255 : ! Read the data
3256 0 : ABI_MALLOC(ekb0,(dimekb,dimekb))
3257 0 : ekb0(:,:)=zero
3258 0 : do itypat=1,ntypat
3259 0 : read (nunit, '(13x,i4)' )ij
3260 0 : do iproj=1,nproj
3261 0 : read (nunit, '(6x,3d22.14)' )(ekb0(iproj+nproj*(ii-1),iproj+nproj*(ii-1)),ii=1,min(npsang,3))
3262 0 : if(npsang>3)read (nunit, '(6x,3d22.14)' )(ekb0(iproj+nproj*(ii-1),iproj+nproj*(ii-1)),ii=4,npsang)
3263 0 : do ii=1,npsang
3264 0 : iekb=iproj+nproj*(ii-1)
3265 0 : indlmn(1,iekb,itypat)=ii-1
3266 0 : indlmn(2,iekb,itypat)=ii-1
3267 0 : indlmn(3,iekb,itypat)=iproj
3268 0 : indlmn(4,iekb,itypat)=ii**2-ii+1
3269 0 : indlmn(5,iekb,itypat)=iekb
3270 0 : indlmn(6,iekb,itypat)=1
3271 : ! For the time being, only diagonal ekb are treated in abinit v3
3272 0 : ekb(iekb,itypat)=ekb0(iekb,iekb)
3273 : end do
3274 : end do
3275 : end do
3276 0 : ABI_FREE(ekb0)
3277 :
3278 : ! --------------------------------------------
3279 : ! ----- OTHER CASES -------------------------
3280 : ! --------------------------------------------
3281 0 : else if(string==' No informat')then
3282 0 : with_psps=0
3283 : else
3284 0 : ABI_BUG('Error when reading the psp information')
3285 : end if
3286 :
3287 : ! Now, the number of blocks
3288 386 : read(nunit,*)
3289 386 : read(nunit,*)
3290 386 : read(nunit, '(24x,i4)' )nblok
3291 :
3292 : ! ==================================================================================
3293 : ! Second option: write psp characteristic to file ==================================
3294 : ! ==================================================================================
3295 150 : else if(choice==2)then
3296 :
3297 150 : write(nunit, '(a)' )' '
3298 150 : if (with_psps==0)then
3299 : ! This possibility is used when the DDB is initialized,
3300 : ! and the ekb s are not available from the GS input file...
3301 86 : write(nunit, '(a)' )' No information on the potentials yet '
3302 : else
3303 :
3304 : ! ==============================
3305 : ! === Norm-conserving psps =====
3306 : ! ==============================
3307 : ! GA: Something is wrong here.
3308 : ! I noticed this with test tutorespfn[tnlo_2]
3309 : ! where the ekb are reported both in the main output
3310 : ! and in the ddb, but some components written to ddb are wrong.
3311 64 : if (usepaw==0) then
3312 63 : write(nunit, '(a)' )' Description of the potentials (KB energies)'
3313 63 : write(nunit, '(a,i6)' )' vrsio8 (for pseudopotentials)=',vrsio8
3314 63 : write(nunit, '(a,i3)' ) ' usepaw =',usepaw
3315 63 : write(nunit, '(a,i3,a,i3,a,i3)' )' dimekb =',dimekb,' lmnmax=',lmnmax
3316 252 : ABI_MALLOC(ekb0,(dimekb,dimekb))
3317 161 : do itypat=1,ntypat
3318 : ! Compute nekb
3319 98 : nekb=0
3320 788 : do jlmn=1,lmnmax
3321 690 : jln=indlmn(5,jlmn,itypat)
3322 788 : if(jln>nekb)then
3323 428 : nekb=jln
3324 : end if
3325 : end do
3326 98 : write(nunit, '(a,i4,a,i3,a,i4)' )' Atom type= ',itypat,' pspso=',pspso(itypat),' nekb=',nekb
3327 98 : write(nunit, '(a)' ) ' iln lpsang iproj ekb(:)'
3328 98 : iln0=0
3329 4918 : ekb0(:,:)=zero
3330 851 : do ilmn=1,lmnmax
3331 690 : iln =indlmn(5,ilmn,itypat)
3332 788 : if (iln>iln0) then
3333 428 : iln0=iln
3334 428 : lpsang=indlmn(1,ilmn,itypat)
3335 428 : iproj=indlmn(3,ilmn,itypat)
3336 : ! For the time being, only diagonal ekb are treated in abinit v3
3337 428 : ekb0(iln,iln)=ekb(iln,itypat)
3338 : ! For non-diagonal ekb, one could use:
3339 : ! do ii=iln to nekb
3340 : ! ekb0(ii,iln)=ekb(ii+iln*(iln-1)/2,itypat)
3341 : ! end do
3342 428 : write(nunit, '(3i6,3x,4es15.7)' ) iln,lpsang,iproj,(ekb0(ii,iln),ii=1,min(nekb,4))
3343 428 : if(nekb>4)then
3344 996 : do iekb=5,nekb,4
3345 996 : write(nunit, '(21x,4es15.7)' )(ekb0(ii,iekb),ii=iekb,min(nekb,iekb+3))
3346 : end do
3347 : end if
3348 : end if
3349 : end do
3350 : end do
3351 63 : ABI_FREE(ekb0)
3352 :
3353 : ! ==============================
3354 : ! ============ PAW =============
3355 : ! ==============================
3356 : else
3357 1 : write(nunit, '(a)' )' Description of the PAW dataset(s)'
3358 1 : write(nunit, '(a,i6)' )' vrsio8 (for pseudopotentials)=',vrsio8
3359 1 : write(nunit, '(a,i3)' ) ' usepaw =',usepaw
3360 3 : do itypat=1,ntypat
3361 2 : iln0=0
3362 6 : ABI_MALLOC(orbitals,(pawtab(itypat)%basis_size))
3363 18 : do ilmn=1,pawtab(itypat)%lmn_size
3364 16 : iln =indlmn(5,ilmn,itypat)
3365 18 : if (iln>iln0) then
3366 8 : iln0=iln;orbitals(iln)=indlmn(1,ilmn,itypat)
3367 : end if
3368 : end do
3369 : write(nunit, '(a,i4,a,i3,a,i5)' ) &
3370 2 : ' Atom type=',itypat,' basis_size=',pawtab(itypat)%basis_size,&
3371 4 : ' lmn_size=',pawtab(itypat)%lmn_size
3372 : write(nunit, '(a,50i2)' ) &
3373 2 : ' Basis functions=',orbitals(1:pawtab(itypat)%basis_size)
3374 : write(nunit, '(a,f6.3,a,i2,a,f6.3)' ) &
3375 2 : ' r_PAW= ',pawtab(itypat)%rpaw, &
3376 2 : ' shape_type= ',pawtab(itypat)%shape_type,&
3377 4 : ' r_shape= ',pawtab(itypat)%rshp
3378 2 : nekb=0
3379 6 : ABI_MALLOC(dij0,(pawtab(itypat)%lmn2_size))
3380 6 : ABI_MALLOC(i1,(pawtab(itypat)%lmn2_size))
3381 4 : ABI_MALLOC(i2,(pawtab(itypat)%lmn2_size))
3382 18 : do jlmn=1,pawtab(itypat)%lmn_size
3383 16 : ij=jlmn*(jlmn-1)/2
3384 90 : do ilmn=1,jlmn
3385 88 : if (abs(pawtab(itypat)%dij0(ij+ilmn))>tol16) then
3386 24 : nekb=nekb+1;i1(nekb)=ilmn;i2(nekb)=jlmn
3387 24 : dij0(nekb)=pawtab(itypat)%dij0(ij+ilmn)
3388 : end if
3389 : end do
3390 : end do
3391 2 : write(nunit,'(a,i3,a)') ' Dij0= (only the ',nekb,' values different from zero)'
3392 2 : write(nunit,'(2a)') ' i j Dij0 i j Dij0 ',&
3393 4 : ' i j Dij0 i j Dij0'
3394 8 : do ii=1,nekb,4
3395 32 : write(nunit,'(3x,4(1x,i4,1x,i4,1x,es12.5))') (i1(ij),i2(ij),dij0(ij),ij=ii,min(ii+3,nekb))
3396 : end do
3397 2 : ABI_FREE(dij0)
3398 2 : ABI_FREE(i1)
3399 2 : ABI_FREE(i2)
3400 3 : ABI_FREE(orbitals)
3401 : end do
3402 :
3403 : end if ! NCPP or PAW
3404 : end if ! with_psps==0
3405 :
3406 : ! Now, write the number of blocks
3407 150 : write(nunit, '(a)' )' '
3408 150 : write(nunit, '(a)' )' **** Database of total energy derivatives ****'
3409 150 : write(nunit, '(a,i4)' ) ' Number of data blocks= ',nblok
3410 : end if
3411 :
3412 536 : end subroutine psddb8
3413 : !!***
3414 :
3415 : !----------------------------------------------------------------------
3416 :
3417 : !!****f* m_ddb_hdr/ioddb8_in
3418 : !!
3419 : !! NAME
3420 : !! ioddb8_in
3421 : !!
3422 : !! FUNCTION
3423 : !! Open Derivative DataBase, and read preliminary information.
3424 : !! Note: only one processor reads the DDB.
3425 : !!
3426 : !! INPUTS
3427 : !! character(len=*)= name of input file
3428 : !! matom=maximum number of atoms
3429 : !! mband=maximum number of bands
3430 : !! mkpt=maximum number of special points
3431 : !! msym=maximum number of symetries
3432 : !! mtypat=maximum number of atom types
3433 : !! unddb=unit number for input
3434 : !!
3435 : !! OUTPUT
3436 : !! acell(3)=length scales of primitive translations (bohr)
3437 : !! amu(mtypat)=mass of the atoms (atomic mass unit)
3438 : !! ddb_version=version of the ddb file
3439 : !! dilatmx=the maximal dilatation factor
3440 : !! ecut=kinetic energy planewave cutoff (hartree)
3441 : !! ecutsm=smearing energy for plane wave kinetic energy (Ha)
3442 : !! intxc=control xc quadrature
3443 : !! iscf=parameter controlling scf or non-scf choice
3444 : !! ixc=exchange-correlation choice parameter
3445 : !! kpt(3,mkpt)=k point set (reduced coordinates)
3446 : !! kptnrm=normalisation of k points
3447 : !! natom=number of atoms in the unit cell
3448 : !! nband(mkpt)=number of bands at each k point, for each polarization
3449 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
3450 : !! nkpt=number of k points
3451 : !! nspden=number of spin-density components
3452 : !! nspinor=number of spinorial components of the wavefunctions
3453 : !! nsppol=1 for unpolarized, 2 for spin-polarized
3454 : !! nsym=number of symmetry elements in space group
3455 : !! ntypat=number of atom types
3456 : !! occ(mband*mkpt)=occupation number for each band and k
3457 : !! occopt=option for occupancies
3458 : !! pawecutdg=cut-off for fine "double grid" used in PAW calculations (unused for NCPP)
3459 : !! rprim(3,3)=dimensionless primitive translations in real space
3460 : !! dfpt_sciss=scissor shift (Ha)
3461 : !! spinat(3,matom)=initial spin of each atom, in unit of hbar/2
3462 : !! symafm(msym)=(anti)ferromagnetic part of symmetry operations
3463 : !! symrel(3,3,msym)=symmetry operations in real space
3464 : !! tnons(3,msym)=nonsymmorphic translations for symmetry operations
3465 : !! tolwfr=tolerance on largest wf residual
3466 : !! tphysel="physical" electronic temperature with FD occupations
3467 : !! tsmear=smearing width (or temperature) in Hartree
3468 : !! typat(matom)=type of each atom
3469 : !! usepaw=flag for PAW
3470 : !! wtk(mkpt)=weight assigned to each k point
3471 : !! xred(3,matom)=reduced atomic coordinates
3472 : !! zion(mtypat)=valence charge of each type of atom
3473 : !! znucl(mtypat)=atomic number of atom type
3474 : !!
3475 : !! SOURCE
3476 :
3477 386 : subroutine ioddb8_in(filename,matom,mband,mkpt,msym,mtypat,unddb,&
3478 386 : & acell,amu,ddb_version,dilatmx,ecut,ecutsm,intxc,iscf,ixc,kpt,kptnrm,&
3479 386 : & natom,nband,ngfft,nkpt,nspden,nspinor,nsppol,nsym,ntypat,occ,occopt,&
3480 386 : & pawecutdg,rprim,dfpt_sciss,spinat,symafm,symrel,tnons,tolwfr,tphysel,tsmear,&
3481 386 : & typat,usepaw,wtk,xred,zion,znucl)
3482 :
3483 : !Arguments -------------------------------
3484 : !scalars
3485 : integer,intent(in) :: matom,mband,mkpt,msym,mtypat,unddb,usepaw
3486 : integer,intent(out) :: ddb_version,intxc,iscf,ixc,natom,nkpt,nspden,nspinor,nsppol,nsym,ntypat,occopt
3487 : real(dp),intent(out) :: dilatmx,ecut,ecutsm,pawecutdg,kptnrm,dfpt_sciss,tolwfr,tphysel,tsmear
3488 : character(len=*),intent(in) :: filename
3489 : !arrays
3490 : integer,intent(out) :: nband(mkpt),ngfft(18),symafm(msym),symrel(3,3,msym),typat(matom)
3491 : real(dp),intent(out) :: acell(3),amu(mtypat),kpt(3,mkpt),occ(mband*mkpt)
3492 : real(dp),intent(out) :: rprim(3,3),spinat(3,matom),tnons(3,msym),wtk(mkpt)
3493 : real(dp),intent(out) :: xred(3,matom),zion(mtypat),znucl(mtypat)
3494 :
3495 : !Local variables -------------------------
3496 : !Set routine version number here:
3497 : !scalars
3498 : integer,parameter :: vrsio8=100401,vrsio8_old=010929,vrsio8_old_old=990527 ! should I modify this ? I guess not
3499 : integer,parameter :: cvrsio9=20230401,cvrsio8=20100401,cvrsio8_old=20010929,cvrsio8_old_old=19990527 ! should I modify this ? I guess not
3500 : integer,parameter :: cvrsio9_new=20240201
3501 : integer :: bantot,ddbvrs,iband,ii,ij,ikpt,iline,im,ndig,usepaw0
3502 : logical :: ddbvrs_is_current_or_old,testn,testv
3503 : character(len=500) :: message
3504 : character(len=6) :: name_old, ddbvrs6
3505 : character(len=8) :: ddbvrs8
3506 : character(len=3) :: prefix
3507 : !arrays
3508 : character(len=12) :: name(9)
3509 : ! *********************************************************************
3510 :
3511 : !Open the input derivative database.
3512 : !write(message,'(a,a)')' About to open file ',TRIM(filename)
3513 : !call wrtout(std_out,message,'COLL')
3514 386 : if (open_file(filename,message,unit=unddb,form="formatted",status="old",action="read") /= 0) then
3515 0 : ABI_ERROR(message)
3516 : end if
3517 :
3518 : !Check the compatibility of the input DDB with the DDB code
3519 386 : read (unddb,*)
3520 386 : read (unddb,*)
3521 386 : read (unddb, '(20x,i10)' )ddbvrs
3522 :
3523 : !write(std_out,'(a,i10)')' ddbvrs=',ddbvrs
3524 : if(ddbvrs/=cvrsio9_new .and. ddbvrs/=cvrsio9 .and. ddbvrs/=vrsio8 &
3525 386 : & .and. ddbvrs/=vrsio8_old .and. ddbvrs/=vrsio8_old_old)then
3526 : write(message, '(a,i10,2a,4(a,i10),a)' )&
3527 0 : 'The input DDB version number=',ddbvrs,' does not agree',ch10,&
3528 0 : 'with the allowed code DDB version numbers,',cvrsio9,', ',vrsio8,', ',vrsio8_old,' and ',vrsio8_old_old,' .'
3529 0 : ABI_BUG(message)
3530 : end if
3531 :
3532 : !Convert older version to 8 digit format
3533 386 : if (ddbvrs /= cvrsio9 .and. ddbvrs /= cvrsio9_new) then
3534 210 : ndig= int(log10(real(ddbvrs))) + 1
3535 210 : write(ddbvrs6,'(i0)') ddbvrs
3536 210 : if (ddbvrs==vrsio8 .or.ddbvrs==vrsio8_old) then
3537 210 : if (ndig==6) then
3538 155 : write(prefix,'(i2)') 20
3539 55 : else if (ndig==5) then
3540 55 : write(prefix,'(i3)') 200
3541 : end if
3542 0 : else if (ddbvrs==vrsio8_old_old) then
3543 0 : if (ndig==6) then
3544 0 : write(prefix,'(i2)') 19
3545 0 : else if (ndig==5) then
3546 0 : write(prefix,'(i3)') 199
3547 : end if
3548 : end if
3549 210 : ddbvrs8= trim(prefix) // trim(ddbvrs6)
3550 210 : read(ddbvrs8,'(i8)') ddbvrs
3551 : end if
3552 386 : ddb_version=ddbvrs
3553 :
3554 : !Read the 4 n-integers, also testing the names of data, and checking that their value is acceptable.
3555 : !This is important to insure that any array has a sufficient dimension.
3556 386 : read (unddb,*)
3557 386 : read (unddb,*)
3558 386 : read (unddb,*)
3559 386 : testn=.true.; testv=.true.
3560 : ! ddbvrs_is_current_or_old=(ddbvrs==vrsio8.or.ddbvrs==vrsio8_old)
3561 386 : ddbvrs_is_current_or_old=(ddbvrs>=cvrsio8_old)
3562 :
3563 : !1. usepaw
3564 386 : if(ddbvrs>=cvrsio8)then
3565 331 : read (unddb, '(1x,a9,i10)' )name(1),usepaw0
3566 : else
3567 55 : usepaw0=0;name(1)=' usepaw'
3568 : end if
3569 386 : if(name(1)/=' usepaw')testn=.false.
3570 386 : if(usepaw0/=usepaw)testv=.false.
3571 : !2. natom
3572 386 : if(ddbvrs_is_current_or_old)then
3573 386 : read (unddb, '(1x,a9,i10)' )name(2),natom
3574 : else
3575 0 : read (unddb, '(1x,a6,i10)' )name_old,natom ; name(2)=' '//name_old
3576 : end if
3577 386 : if(name(2)/=' natom')testn=.false.
3578 386 : if(natom<=0.or.natom>matom)testv=.false.
3579 : !3. nkpt
3580 386 : if(ddbvrs_is_current_or_old)then
3581 386 : read (unddb, '(1x,a9,i10)' )name(3),nkpt
3582 : else
3583 0 : read (unddb, '(1x,a6,i10)' )name_old,nkpt ; name(3)=' '//name_old
3584 : end if
3585 386 : if(name(3)/=' nkpt')testn=.false.
3586 386 : if(nkpt <=0.or.nkpt >mkpt )testv=.false.
3587 : !4. nsppol
3588 386 : if(ddbvrs_is_current_or_old)then
3589 386 : read (unddb, '(1x,a9,i10)' )name(4),nsppol
3590 : else
3591 0 : read (unddb, '(1x,a6,i10)' )name_old,nsppol ; name(4)=' '//name_old
3592 : end if
3593 386 : if(name(4)/=' nsppol')testn=.false.
3594 386 : if(nsppol<=0.or.nsppol>2)testv=.false.
3595 : !5. nsym
3596 386 : if(ddbvrs_is_current_or_old)then
3597 386 : read (unddb, '(1x,a9,i10)' )name(5),nsym
3598 : else
3599 0 : read (unddb, '(1x,a6,i10)' )name_old,nsym ; name(5)=' '//name_old
3600 : end if
3601 386 : if(name(5)/=' nsym')testn=.false.
3602 386 : if(nsym <=0.or.nsym >msym )testv=.false.
3603 : !6. ntypat
3604 386 : if(ddbvrs_is_current_or_old)then
3605 386 : read (unddb, '(1x,a9,i10)' )name(6),ntypat
3606 : else
3607 0 : read (unddb, '(1x,a6,i10)' )name_old,ntypat ; name(6)=' '//name_old
3608 : end if
3609 386 : if(name(6)/=' ntypat' .and. name(6)/=' ntype')testn=.false.
3610 386 : if(ntypat<=0.or.ntypat>mtypat)testv=.false.
3611 : !7. occopt
3612 : !Before reading nband, the last parameters that define
3613 : !the dimension of some array, need to know what is their
3614 : !representation, given by occopt
3615 386 : if(ddbvrs_is_current_or_old)then
3616 386 : read (unddb, '(1x,a9,i10)' )name(7),occopt
3617 : else
3618 0 : read (unddb, '(1x,a6,i10)' )name_old,occopt ; name(7)=' '//name_old
3619 : end if
3620 386 : if(name(7)/=' occopt')testn=.false.
3621 386 : if(occopt<0.or.occopt>8)testv=.false.
3622 : !Message if the names or values are not right
3623 386 : if (.not.testn.or..not.testv) then
3624 0 : write(message, '(a,a,a)' )' ioddb8_in : An error has been found in one',ch10,&
3625 0 : ' of the positive n-integers contained in the DDB : '
3626 0 : call wrtout(std_out,message,'COLL')
3627 : write(message, '(a)' )&
3628 0 : ' Expected Found '
3629 0 : call wrtout(std_out,message,'COLL')
3630 : write(message, '(a,i10,a,a,a,i10)' )&
3631 0 : ' usepaw equal to ',usepaw,' ',trim(name(1)),' =',usepaw0
3632 0 : call wrtout(std_out,message,'COLL')
3633 : write(message, '(a,i10,a,a,a,i10)' )&
3634 0 : ' natom , lower than',matom+1,' ',trim(name(2)),' =',natom
3635 0 : call wrtout(std_out,message,'COLL')
3636 : write(message, '(a,i10,a,a,a,i10)' )&
3637 0 : ' nkpt , lower than',mkpt+1 ,' ',trim(name(3)),' =',nkpt
3638 0 : call wrtout(std_out,message,'COLL')
3639 : write(message, '(a,i10,a,a,a,i10)' )&
3640 0 : ' nsppol, lower than',3 ,' ',trim(name(4)),' =',nsppol
3641 0 : call wrtout(std_out,message,'COLL')
3642 : write(message, '(a,i10,a,a,a,i10)' )&
3643 0 : ' nsym , lower than',msym+1 ,' ',trim(name(5)),' =',nsym
3644 0 : call wrtout(std_out,message,'COLL')
3645 : write(message, '(a,i10,a,a,a,i10)' )&
3646 0 : ' ntypat, lower than',mtypat+1,' ',trim(name(6)),' =',ntypat
3647 0 : call wrtout(std_out,message,'COLL')
3648 : write(message, '(a,a,a,i10)' )&
3649 0 : ' occopt, between 0 and 7 ',trim(name(7)),' =',occopt
3650 0 : call wrtout(std_out,message,'COLL')
3651 :
3652 0 : ABI_ERROR('See the error message above.')
3653 : end if
3654 :
3655 : !One more set of parameters define the dimensions of the
3656 : !array : nband. Morever, it depends on occopt and nkpt, and has to be
3657 : !tested after the test on nkpt is performed.
3658 : !8. nband
3659 386 : if(occopt==2)then
3660 9 : im=12
3661 74 : do iline=1,(nkpt+11)/12
3662 650 : name(:) = "" ! reset all line strings
3663 65 : if(iline==(nkpt+11)/12)im=nkpt-12*(iline-1)
3664 65 : if(ddbvrs_is_current_or_old)then
3665 833 : read (unddb, '(1x,a9,5x,12i5)' )name(1),(nband((iline-1)*12+ii),ii=1,im)
3666 : else
3667 0 : read (unddb, '(1x,a6,5x,12i5)' )name_old,(nband((iline-1)*12+ii),ii=1,im) ; name(1)=' '//name_old
3668 : end if
3669 74 : if (iline==1) then
3670 9 : call ddb_chkname(name(1),' nband')
3671 : else
3672 56 : call ddb_chkname(name(1),' ')
3673 : end if
3674 : end do
3675 : else
3676 3770 : name(:) = "" ! reset all line strings
3677 377 : if(ddbvrs_is_current_or_old)then
3678 377 : read (unddb, '(1x,a9,i10)' )name(1),nband(1)
3679 : else
3680 0 : read (unddb, '(1x,a6,i10)' )name_old,nband(1) ; name(1)=' '//name_old
3681 : end if
3682 377 : call ddb_chkname(name(1),' nband')
3683 377 : if(nkpt>1)then
3684 73501 : do ikpt=2,nkpt
3685 73501 : nband(ikpt)=nband(1)
3686 : end do
3687 : end if
3688 : end if
3689 :
3690 : !check all nband values, and sum them
3691 386 : bantot=0
3692 74656 : do ikpt=1,nkpt
3693 74270 : if(nband(ikpt)<0)then
3694 : write(message, '(a,i4,a,i4,3a)' )&
3695 0 : & 'For ikpt = ',ikpt,' nband = ',nband(ikpt),' is negative.',ch10,&
3696 0 : & 'Action: correct your DDB.'
3697 0 : ABI_ERROR(message)
3698 74270 : else if(nband(ikpt)>mband)then
3699 : write(message, '(a,i4,a,i4,a,a,i4,3a)' )&
3700 0 : & 'For ikpt = ',ikpt,', nband = ',nband(ikpt),ch10,&
3701 0 : & 'is larger than mband = ',mband,'.',ch10,&
3702 0 : & 'Action: recompile the calling code with a larger mband.'
3703 0 : ABI_ERROR(message)
3704 : end if
3705 74656 : bantot=bantot+nband(ikpt)
3706 : end do
3707 :
3708 : !Read the rest of variables, with check of the names
3709 : !9. acell
3710 386 : if(ddbvrs_is_current_or_old)then
3711 386 : read (unddb, '(1x,a9,3d22.14)' )name(1),acell
3712 : else
3713 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,acell ; name(1)=' '//name_old
3714 : end if
3715 386 : call ddb_chkname(name(1),' acell')
3716 : !9. amu
3717 386 : im=3
3718 773 : do iline=1,(ntypat+2)/3
3719 387 : if(iline==(ntypat+2)/3)im=ntypat-3*(iline-1)
3720 387 : if(ddbvrs_is_current_or_old)then
3721 1126 : read (unddb, '(1x,a9,3d22.14)' )name(1),(amu((iline-1)*3+ii),ii=1,im)
3722 : else
3723 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,(amu((iline-1)*3+ii),ii=1,im) ; name(1)=' '//name_old
3724 : end if
3725 773 : if (iline==1) then
3726 386 : call ddb_chkname(name(1),' amu')
3727 : else
3728 1 : call ddb_chkname(name(1),' ')
3729 : end if
3730 : end do
3731 : !11. dilatmx
3732 386 : if(ddbvrs_is_current_or_old)then
3733 386 : read (unddb, '(1x,a9,d22.14)' )name(1),dilatmx
3734 386 : call ddb_chkname(name(1),' dilatmx')
3735 : else
3736 0 : dilatmx=one
3737 : end if
3738 : !12. ecut
3739 : if(ddbvrs_is_current_or_old)then
3740 386 : read (unddb, '(1x,a9,d22.14)' )name(1),ecut
3741 : else
3742 0 : read (unddb, '(1x,a6,d22.14)' )name_old,ecut ; name(1)=' '//name_old
3743 : end if
3744 386 : call ddb_chkname(name(1),' ecut')
3745 : !12b. pawecutdg (PAW only)
3746 386 : if(ddbvrs>=cvrsio8.and.usepaw==1) then
3747 5 : read (unddb, '(1x,a9,d22.14)' )name(1),pawecutdg
3748 : else
3749 381 : pawecutdg=ecut;name(1)='pawecutdg'
3750 : end if
3751 386 : call ddb_chkname(name(1),'pawecutdg')
3752 : !13. ecutsm
3753 386 : if(ddbvrs_is_current_or_old)then
3754 386 : read (unddb, '(1x,a9,d22.14)' )name(1),ecutsm
3755 386 : call ddb_chkname(name(1),' ecutsm')
3756 : else
3757 0 : ecutsm=zero
3758 : end if
3759 : !14. intxc
3760 : if(ddbvrs_is_current_or_old)then
3761 386 : read (unddb, '(1x,a9,i10)' )name(1),intxc
3762 386 : call ddb_chkname(name(1),' intxc')
3763 : else
3764 0 : intxc=1
3765 : end if
3766 : !15. iscf
3767 : if(ddbvrs_is_current_or_old)then
3768 386 : read (unddb, '(1x,a9,i10)' )name(1),iscf
3769 : else
3770 0 : read (unddb, '(1x,a6,i10)' )name_old,iscf ; name(1)=' '//name_old
3771 : end if
3772 386 : call ddb_chkname(name(1),' iscf')
3773 : !16. ixc
3774 386 : if(ddbvrs_is_current_or_old)then
3775 386 : read (unddb, '(1x,a9,i10)' )name(1),ixc
3776 : else
3777 0 : read (unddb, '(1x,a6,i10)' )name_old,ixc ; name(1)=' '//name_old
3778 : end if
3779 386 : call ddb_chkname(name(1),' ixc')
3780 : !17. kpt
3781 74656 : do iline=1,nkpt
3782 74270 : if(ddbvrs_is_current_or_old)then
3783 74270 : read (unddb, '(1x,a9,3d22.14)' )name(1),(kpt(ii,iline),ii=1,3)
3784 : else
3785 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,(kpt(ii,iline),ii=1,3) ; name(1)=' '//name_old
3786 : end if
3787 74656 : if (iline==1) then
3788 386 : call ddb_chkname(name(1),' kpt')
3789 : else
3790 73884 : call ddb_chkname(name(1),' ')
3791 : end if
3792 : end do
3793 : !18. kptnrm
3794 386 : if(ddbvrs_is_current_or_old)then
3795 386 : read (unddb, '(1x,a9,d22.14)' )name(1),kptnrm
3796 : else
3797 0 : read (unddb, '(1x,a6,d22.14)' )name_old,kptnrm ; name(1)=' '//name_old
3798 : end if
3799 386 : call ddb_chkname(name(1),' kptnrm')
3800 : !19. ngfft
3801 386 : if(ddbvrs_is_current_or_old)then
3802 386 : read (unddb, '(1x,a9,5x,3i5)' )name(1),ngfft(1:3)
3803 : else
3804 0 : read (unddb, '(1x,a6,5x,3i5)' )name_old,ngfft(1:3) ; name(1)=' '//name_old
3805 : end if
3806 : !For the time being, do not check the validity of the name,
3807 : !in order to accept both ng and ngfft
3808 : !20. nspden
3809 386 : if(ddbvrs_is_current_or_old)then
3810 386 : read (unddb, '(1x,a9,i10)' )name(1),nspden
3811 386 : call ddb_chkname(name(1),' nspden')
3812 : else
3813 0 : nspden=0
3814 : end if
3815 : !21. nspinor
3816 : if(ddbvrs_is_current_or_old)then
3817 386 : read (unddb, '(1x,a9,i10)' )name(1),nspinor
3818 386 : call ddb_chkname(name(1),' nspinor')
3819 : else
3820 0 : nspinor=0
3821 : end if
3822 : !22. occ
3823 386 : if(occopt==2)then
3824 9 : im=3
3825 5145 : do iline=1,(bantot+2)/3
3826 5136 : if(iline==(bantot+2)/3)im=bantot-3*(iline-1)
3827 5136 : if(ddbvrs_is_current_or_old)then
3828 20544 : read (unddb, '(1x,a9,3d22.14)' )name(1),(occ((iline-1)*3+ii),ii=1,im)
3829 : else
3830 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,(occ((iline-1)*3+ii),ii=1,im) ; name(1)=' '//name_old
3831 : end if
3832 5145 : if (iline==1) then
3833 9 : call ddb_chkname(name(1),' occ')
3834 : else
3835 5127 : call ddb_chkname(name(1),' ')
3836 : end if
3837 : end do
3838 : else
3839 377 : im=3
3840 2163 : do iline=1,(nband(1)+2)/3
3841 1786 : if(iline==(nband(1)+2)/3)im=nband(1)-3*(iline-1)
3842 1786 : if(ddbvrs_is_current_or_old)then
3843 6836 : read (unddb, '(1x,a9,3d22.14)' )name(1),(occ((iline-1)*3+ii),ii=1,im)
3844 : else
3845 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,(occ((iline-1)*3+ii),ii=1,im) ; name(1)=' '//name_old
3846 : end if
3847 2163 : if (iline==1) then
3848 377 : call ddb_chkname(name(1),' occ')
3849 : else
3850 1409 : call ddb_chkname(name(1),' ')
3851 : end if
3852 : end do
3853 377 : if(nkpt>1)then
3854 73501 : do ikpt=2,nkpt
3855 1457821 : do iband=1,nband(1)
3856 1457445 : occ(iband+nband(1)*(ikpt-1))=occ(iband)
3857 : end do
3858 : end do
3859 : end if
3860 : end if
3861 : !23. rprim
3862 1544 : do iline=1,3
3863 1158 : if(ddbvrs_is_current_or_old)then
3864 1158 : read (unddb, '(1x,a9,3d22.14)' )name(1),(rprim(ii,iline),ii=1,3)
3865 : else
3866 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,(rprim(ii,iline),ii=1,3) ; name(1)=' '//name_old
3867 : end if
3868 1544 : if (iline==1) then
3869 386 : call ddb_chkname(name(1),' rprim')
3870 : else
3871 772 : call ddb_chkname(name(1),' ')
3872 : end if
3873 : end do
3874 : !24. dfpt_sciss
3875 386 : if(ddbvrs_is_current_or_old)then
3876 386 : read (unddb, '(1x,a11,d22.14)' )name(1),dfpt_sciss
3877 386 : call ddb_chkname(name(1),'dfpt_sciss', 'sciss')
3878 : else
3879 0 : read (unddb, '(1x,a6,d22.14)' )name_old,dfpt_sciss ; name(1)=name_old
3880 0 : call ddb_chkname(name(1),'sciss')
3881 : end if
3882 : !25. spinat
3883 386 : if(ddbvrs_is_current_or_old)then
3884 1505 : do iline=1,natom
3885 1119 : read (unddb, '(1x,a9,3d22.14)' )name(1),(spinat(ii,iline),ii=1,3)
3886 1505 : if (iline==1) then
3887 386 : call ddb_chkname(name(1),' spinat')
3888 : else
3889 733 : call ddb_chkname(name(1),' ')
3890 : end if
3891 : end do
3892 : else
3893 : ! spinat is set to zero by default in mrgddb.f
3894 : ! spinat(:,1:natom)=zero
3895 : end if
3896 : !26. symafm
3897 : if(ddbvrs_is_current_or_old)then
3898 386 : im=12
3899 1611 : do iline=1,(nsym+11)/12
3900 1225 : if(iline==(nsym+11)/12)im=nsym-12*(iline-1)
3901 15147 : read (unddb, '(1x,a9,5x,12i5)' )name(1),(symafm((iline-1)*12+ii),ii=1,im)
3902 1611 : if (iline==1) then
3903 386 : call ddb_chkname(name(1),' symafm')
3904 : else
3905 839 : call ddb_chkname(name(1),' ')
3906 : end if
3907 : end do
3908 : else
3909 : ! symafm is set to 1 by default in mrgddb.f
3910 : ! symafm(1:nsym)=1
3911 : end if
3912 : !27. symrel
3913 14308 : do iline=1,nsym
3914 13922 : if(ddbvrs_is_current_or_old)then
3915 13922 : read (unddb, '(1x,a9,5x,9i5)' )name(1),((symrel(ii,ij,iline),ii=1,3),ij=1,3)
3916 : else
3917 0 : read (unddb, '(1x,a6,5x,9i5)' )name_old,&
3918 0 : & ((symrel(ii,ij,iline),ii=1,3),ij=1,3) ; name(1)=' '//name_old
3919 : end if
3920 14308 : if (iline==1) then
3921 386 : call ddb_chkname(name(1),' symrel')
3922 : else
3923 13536 : call ddb_chkname(name(1),' ')
3924 : end if
3925 : end do
3926 : !28old. xred
3927 386 : if(.not.ddbvrs_is_current_or_old)then
3928 0 : do iline=1,natom
3929 0 : read (unddb, '(1x,a6,3d22.14)' )name(1),(xred(ii,iline),ii=1,3)
3930 : end do
3931 : ! No check of name, to allow the old tn
3932 : end if
3933 : !28. tnons
3934 14308 : do iline=1,nsym
3935 13922 : if(ddbvrs_is_current_or_old)then
3936 13922 : read (unddb, '(1x,a9,3d22.14)' )name(1),(tnons(ii,iline),ii=1,3)
3937 : else
3938 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,(tnons(ii,iline),ii=1,3) ; name(1)=' '//name_old
3939 : end if
3940 14308 : if (iline==1) then
3941 386 : call ddb_chkname(name(1),' tnons')
3942 : else
3943 13536 : call ddb_chkname(name(1),' ')
3944 : end if
3945 : end do
3946 : !29. tolwfr
3947 386 : if(ddbvrs_is_current_or_old)then
3948 386 : read (unddb, '(1x,a9,d22.14)' )name(1),tolwfr
3949 : end if
3950 : !Do not check the name, in order to allow both tolwfr and wftol
3951 : !30. tphysel
3952 : if(ddbvrs_is_current_or_old)then
3953 386 : read (unddb, '(1x,a9,d22.14)' )name(1),tphysel
3954 386 : call ddb_chkname(name(1),' tphysel')
3955 : else
3956 0 : tphysel=zero
3957 : end if
3958 : !31. tsmear
3959 : if(ddbvrs_is_current_or_old)then
3960 386 : read (unddb, '(1x,a9,d22.14)' )name(1),tsmear
3961 386 : call ddb_chkname(name(1),' tsmear')
3962 : else
3963 0 : tsmear=zero
3964 : end if
3965 : !32. typat
3966 386 : im=12
3967 773 : do iline=1,(natom+11)/12
3968 387 : if(iline==(natom+11)/12)im=natom-12*(iline-1)
3969 387 : if(ddbvrs_is_current_or_old)then
3970 1506 : read (unddb, '(1x,a9,5x,12i5)' )name(1),(typat((iline-1)*12+ii),ii=1,im)
3971 : else
3972 0 : read (unddb, '(1x,a6,5x,12i5)' )name_old,(typat((iline-1)*12+ii),ii=1,im) ; name(1)=' '//name_old
3973 : end if
3974 773 : if (iline==1) then
3975 : ! Both type and typat are allowed => no check
3976 : ! call ddb_chkname(name(1),' typat')
3977 : else
3978 1 : call ddb_chkname(name(1),' ')
3979 : end if
3980 : end do
3981 : !33old. tolwfr
3982 386 : if(.not.ddbvrs_is_current_or_old)then
3983 0 : read (unddb, '(1x,a6,d22.14)' )name(1),tolwfr
3984 : end if
3985 : !Do not check the name, in order to allow both tolwfr and wftol
3986 : !33. wtk
3987 386 : im=3
3988 25304 : do iline=1,(nkpt+2)/3
3989 24918 : if(iline==(nkpt+2)/3)im=nkpt-3*(iline-1)
3990 24918 : if(ddbvrs_is_current_or_old)then
3991 99188 : read (unddb, '(1x,a9,3d22.14)' )name(1),(wtk((iline-1)*3+ii),ii=1,im)
3992 : else
3993 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,(wtk((iline-1)*3+ii),ii=1,im) ; name(1)=' '//name_old
3994 : end if
3995 25304 : if (iline==1) then
3996 386 : call ddb_chkname(name(1),' wtk')
3997 : else
3998 24532 : call ddb_chkname(name(1),' ')
3999 : end if
4000 : end do
4001 : !34. xred
4002 386 : if(ddbvrs_is_current_or_old)then
4003 1505 : do iline=1,natom
4004 1119 : read (unddb, '(1x,a9,3d22.14)' )name(1),(xred(ii,iline),ii=1,3)
4005 1505 : if (iline==1) then
4006 386 : call ddb_chkname(name(1),' xred')
4007 : else
4008 733 : call ddb_chkname(name(1),' ')
4009 : end if
4010 : end do
4011 : end if
4012 : !35. znucl
4013 : if(ddbvrs_is_current_or_old)then
4014 386 : im=3
4015 773 : do iline=1,(ntypat+2)/3
4016 387 : if(iline==(ntypat+2)/3)im=ntypat-3*(iline-1)
4017 1126 : read (unddb, '(1x,a9,3d22.14)' )name(1),(znucl((iline-1)*3+ii),ii=1,im)
4018 773 : if (iline==1) then
4019 386 : call ddb_chkname(name(1),' znucl')
4020 : else
4021 1 : call ddb_chkname(name(1),' ')
4022 : end if
4023 : end do
4024 : else
4025 : ! znucl is set to zero by default in mrgddb.f
4026 : ! znucl(:)=zero
4027 : end if
4028 : !36. zion
4029 386 : im=3
4030 773 : do iline=1,(ntypat+2)/3
4031 387 : if(iline==(ntypat+2)/3)im=ntypat-3*(iline-1)
4032 387 : if(ddbvrs_is_current_or_old)then
4033 1126 : read (unddb, '(1x,a9,3d22.14)' )name(1),(zion((iline-1)*3+ii),ii=1,im)
4034 : else
4035 0 : read (unddb, '(1x,a6,3d22.14)' )name_old,(zion((iline-1)*3+ii),ii=1,im) ; name(1)=' '//name_old
4036 : end if
4037 773 : if (iline==1) then
4038 : ! Do not check the names, to allow both zion and znucl - the latter for 990527 format
4039 : ! call ddb_chkname(name(1),' zion')
4040 : else
4041 1 : call ddb_chkname(name(1),' ')
4042 : end if
4043 : end do
4044 :
4045 386 : end subroutine ioddb8_in
4046 : !!***
4047 :
4048 : !----------------------------------------------------------------------
4049 :
4050 : !!****f* m_ddb_hdr/ddb_getdims
4051 : !! NAME
4052 : !! ddb_getdims
4053 : !!
4054 : !! FUNCTION
4055 : !! Open Derivative DataBase, then reads the variables that
4056 : !! must be known in order to dimension the arrays, and close the file.
4057 : !!
4058 : !! INPUTS
4059 : !! character(len=*) filename: name of input or output file
4060 : !! unddb=unit number for input or output
4061 : !!
4062 : !! OUTPUT
4063 : !! dimekb=dimension of ekb (only used for norm-conserving psps)
4064 : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
4065 : !! =if useylm=0, max number of (l,n) comp. over all type of psps
4066 : !! mband=maximum number of bands
4067 : !! msym=maximum number of symmetries
4068 : !! natom=number of atoms
4069 : !! nblok=number of bloks in the DDB
4070 : !! nkpt=number of k points
4071 : !! ntypat=number of atom types
4072 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
4073 : !! comm=MPI communicator.
4074 : !!
4075 : !! SOURCE
4076 :
4077 610 : subroutine ddb_getdims(filename,comm,dimekb,lmnmax,mband,msym,natom,nblok,nkpt,nsppol,ntypat,usepaw,&
4078 : & has_d0E_xx, has_d1E_xx, has_d2E_xx, has_d3E_xx, has_d3E_lw, has_d2eig)
4079 :
4080 : !Arguments -------------------------------
4081 : !scalars
4082 : character(len=*),intent(in) :: filename
4083 : integer,intent(in) :: comm
4084 : integer,intent(out) :: msym,dimekb,lmnmax,mband,natom,nblok,nkpt,ntypat,nsppol,usepaw
4085 : logical,intent(out) :: has_d0E_xx,has_d1E_xx,has_d2E_xx,has_d3E_xx,has_d3E_lw,has_d2eig
4086 :
4087 : !Local variables-------------------------------
4088 : !scalars
4089 : integer,parameter :: master=0
4090 : integer :: ierr,unddb
4091 : ! *********************************************************************
4092 :
4093 : ! Master node reads dims from file and then broadcast.
4094 610 : if (xmpi_comm_rank(comm) == master) then
4095 586 : unddb = get_unit()
4096 : call inprep8(filename,unddb,dimekb,lmnmax,mband,msym,natom,nblok,nkpt,nsppol,ntypat,usepaw,&
4097 586 : has_d0E_xx, has_d1E_xx, has_d2E_xx, has_d3E_xx, has_d3E_lw, has_d2eig)
4098 : end if
4099 :
4100 : ! couldnt we use ddb_hdr_bcast_dim here?
4101 610 : if (xmpi_comm_size(comm) > 1) then
4102 36 : call xmpi_bcast(dimekb, master, comm, ierr)
4103 36 : call xmpi_bcast(lmnmax, master, comm, ierr)
4104 36 : call xmpi_bcast(mband, master, comm, ierr)
4105 36 : call xmpi_bcast(msym, master, comm, ierr)
4106 36 : call xmpi_bcast(natom, master, comm, ierr)
4107 36 : call xmpi_bcast(nblok, master, comm, ierr)
4108 36 : call xmpi_bcast(nkpt, master, comm, ierr)
4109 36 : call xmpi_bcast(ntypat, master, comm, ierr)
4110 36 : call xmpi_bcast(nsppol, master, comm, ierr)
4111 36 : call xmpi_bcast(usepaw, master, comm, ierr)
4112 36 : call xmpi_bcast(has_d0E_xx, master, comm, ierr)
4113 36 : call xmpi_bcast(has_d1E_xx, master, comm, ierr)
4114 36 : call xmpi_bcast(has_d2E_xx, master, comm, ierr)
4115 36 : call xmpi_bcast(has_d3E_xx, master, comm, ierr)
4116 36 : call xmpi_bcast(has_d3E_lw, master, comm, ierr)
4117 36 : call xmpi_bcast(has_d2eig, master, comm, ierr)
4118 : end if
4119 :
4120 : ! Maximum number of perturbations and size of matrix.
4121 : !mpert=natom+6
4122 : !msize=3*mpert*3*mpert; if (has_d3E_xx) msize=msize*3*mpert
4123 :
4124 610 : end subroutine ddb_getdims
4125 : !!***
4126 :
4127 : !----------------------------------------------------------------------
4128 :
4129 : !!****f* m_ddb_hdr/inprep8
4130 : !! NAME
4131 : !! inprep8
4132 : !!
4133 : !! FUNCTION
4134 : !! Open Derivative DataBase, then reads the variables that
4135 : !! must be known in order to dimension the arrays before complete reading
4136 : !! Note: only one processor read or write the DDB.
4137 : !!
4138 : !! INPUTS
4139 : !! character(len=*) filename: name of input or output file
4140 : !! unddb=unit number for input or output
4141 : !!
4142 : !! OUTPUT
4143 : !! dimekb=dimension of ekb (only used for norm-conserving psps)
4144 : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
4145 : !! =if useylm=0, max number of (l,n) comp. over all type of psps
4146 : !! mband=maximum number of bands
4147 : !! msym=maximum number of symmetries
4148 : !! natom=number of atoms
4149 : !! nblok=number of bloks in the DDB
4150 : !! nkpt=number of k points
4151 : !! ntypat=number of atom types
4152 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
4153 : !!
4154 : !! SOURCE
4155 :
4156 :
4157 586 : subroutine inprep8 (filename,unddb,dimekb,lmnmax,mband,msym,natom,nblok,nkpt,&
4158 : & nsppol,ntypat,usepaw,&
4159 : & has_d0E_xx, has_d1E_xx, has_d2E_xx, has_d3E_xx, has_d3E_lw, has_d2eig)
4160 :
4161 : !Arguments -------------------------------
4162 : !scalars
4163 : character(len=*),intent(in) :: filename
4164 : integer,intent(in) :: unddb
4165 : integer,intent(out) :: dimekb,lmnmax,mband,msym,natom,nblok,nkpt,ntypat,nsppol,usepaw
4166 : logical, intent(out) :: has_d0E_xx,has_d1E_xx,has_d2E_xx,has_d3E_xx,has_d3E_lw,has_d2eig
4167 :
4168 : !Local variables -------------------------
4169 : !scalars
4170 : !Set routine version number here:
4171 : integer,parameter :: vrsio8=100401,vrsio8_old=010929,vrsio8_old_old=990527
4172 : integer,parameter :: cvrsio9=20230401,cvrsio8=20100401,cvrsio8_old=20010929,cvrsio8_old_old=19990527
4173 : integer,parameter :: cvrsio9_new=20240201
4174 : integer :: bantot,basis_size0,blktyp,ddbvrs,iband,iblok,iekb,ii,ikpt,iline,im,ios,iproj
4175 : integer :: itypat,itypat0,jekb,lmn_size0,mproj,mpsang,nekb,ndig,nelmts
4176 : integer :: occopt,pspso0,nsym
4177 : logical :: ddbvrs_is_current_or_old,testn,testv
4178 : character(len=12) :: string
4179 : character(len=32) :: blkname
4180 : character(len=500) :: message
4181 : character(len=6) :: name_old, ddbvrs6
4182 : character(len=80) :: rdstring
4183 : character(len=8) :: ddbvrs8
4184 : character(len=3) :: prefix
4185 : !arrays
4186 586 : integer,allocatable :: nband(:)
4187 : character(len=12) :: name(9)
4188 : ! *********************************************************************
4189 :
4190 : !Open the input derivative database.
4191 586 : if (open_file(filename,message,unit=unddb,form="formatted",status="old",action="read") /= 0) then
4192 0 : ABI_ERROR(message)
4193 : end if
4194 :
4195 : !Check the compatibility of the input DDB with the DDB code
4196 586 : read (unddb,*)
4197 586 : read (unddb,*)
4198 586 : read (unddb, '(20x,i10)' )ddbvrs
4199 :
4200 1421 : if (all(ddbvrs/= [cvrsio9_new, cvrsio9, vrsio8, vrsio8_old, vrsio8_old_old]) )then
4201 : write(message, '(a,i10,2a,4(a,i10))' )&
4202 0 : & 'The input DDB version number=',ddbvrs,' does not agree',ch10,&
4203 0 : & 'with the allowed code DDB version numbers,',cvrsio9,', ',vrsio8,', ',vrsio8_old,' and ',vrsio8_old_old
4204 0 : ABI_ERROR(message)
4205 : end if
4206 :
4207 : !Convert older version to 8 digit format
4208 586 : if (ddbvrs /= cvrsio9 .and. ddbvrs /= cvrsio9_new) then
4209 359 : ndig= int(log10(real(ddbvrs))) + 1
4210 359 : write(ddbvrs6,'(i0)') ddbvrs
4211 359 : if (ddbvrs==vrsio8 .or.ddbvrs==vrsio8_old) then
4212 359 : if (ndig==6) then
4213 268 : write(prefix,'(i2)') 20
4214 91 : else if (ndig==5) then
4215 91 : write(prefix,'(i3)') 200
4216 : end if
4217 0 : else if (ddbvrs==vrsio8_old_old) then
4218 0 : if (ndig==6) then
4219 0 : write(prefix,'(i2)') 19
4220 0 : else if (ndig==5) then
4221 0 : write(prefix,'(i3)') 199
4222 : end if
4223 : end if
4224 359 : ddbvrs8= trim(prefix) // trim(ddbvrs6)
4225 359 : read(ddbvrs8,'(i8)') ddbvrs
4226 : end if
4227 :
4228 : !Read the 4 n-integers, also testing the names of data,
4229 : !and checking that their value is acceptable.
4230 : !This is important to insure that any array has a sufficient dimension.
4231 586 : read (unddb,*)
4232 586 : read (unddb,*)
4233 586 : read (unddb,*)
4234 586 : testn=.true.
4235 586 : testv=.true.
4236 : ! ddbvrs_is_current_or_old=(ddbvrs==vrsio8.or.ddbvrs==vrsio8_old)
4237 586 : ddbvrs_is_current_or_old=(ddbvrs>=cvrsio8_old)
4238 :
4239 : !1. usepaw
4240 586 : if(ddbvrs>=cvrsio8)then
4241 495 : read (unddb, '(1x,a9,i10)' )name(1),usepaw
4242 : else
4243 91 : usepaw=0;name(1)=' usepaw'
4244 : end if
4245 586 : if(name(1)/=' usepaw')testn=.false.
4246 : !2. natom
4247 586 : if(ddbvrs_is_current_or_old)then
4248 586 : read (unddb, '(1x,a9,i10)' )name(2),natom
4249 : else
4250 0 : read (unddb, '(1x,a6,i10)' )name_old,natom ; name(2)=' '//name_old
4251 : end if
4252 586 : if(name(2)/=' natom')testn=.false.
4253 586 : if(natom<=0)testv=.false.
4254 : !3. nkpt
4255 586 : if(ddbvrs_is_current_or_old)then
4256 586 : read (unddb, '(1x,a9,i10)' )name(3),nkpt
4257 : else
4258 0 : read (unddb, '(1x,a6,i10)' )name_old,nkpt ; name(3)=' '//name_old
4259 : end if
4260 586 : if(name(3)/=' nkpt')testn=.false.
4261 586 : if(nkpt <=0)testv=.false.
4262 : !4. nsppol
4263 586 : if(ddbvrs_is_current_or_old)then
4264 586 : read (unddb, '(1x,a9,i10)' )name(4),nsppol
4265 : else
4266 0 : read (unddb, '(1x,a6,i10)' )name_old,nsppol ; name(4)=' '//name_old
4267 : end if
4268 586 : if(name(4)/=' nsppol')testn=.false.
4269 586 : if(nsppol<=0.or.nsppol>2)testv=.false.
4270 : !5. nsym
4271 586 : if(ddbvrs_is_current_or_old)then
4272 586 : read (unddb, '(1x,a9,i10)' )name(5),nsym
4273 : else
4274 0 : read (unddb, '(1x,a6,i10)' )name_old,nsym ; name(5)=' '//name_old
4275 : end if
4276 586 : if(name(5)/=' nsym')testn=.false.
4277 : !MG FIXME Why this and why do we need msym?
4278 586 : msym = 192
4279 586 : if (nsym > msym) msym=nsym
4280 : !if(nsym <=0.or.nsym >msym )testv=.false.
4281 : !6. ntypat
4282 586 : if(ddbvrs_is_current_or_old)then
4283 586 : read (unddb, '(1x,a9,i10)' )name(6),ntypat
4284 : else
4285 0 : read (unddb, '(1x,a6,i10)' )name_old,ntypat ; name(6)=' '//name_old
4286 : end if
4287 586 : if(name(6)/=' ntypat' .and. name(6)/=' ntype')testn=.false.
4288 586 : if(ntypat<=0)testv=.false.
4289 : !7. occopt
4290 : !Before reading nband, the last parameters that define
4291 : !the dimension of some array, need to know what is their
4292 : !representation, given by occopt
4293 586 : if(ddbvrs_is_current_or_old)then
4294 586 : read (unddb, '(1x,a9,i10)' )name(7),occopt
4295 : else
4296 0 : read (unddb, '(1x,a6,i10)' )name_old,occopt ; name(7)=' '//name_old
4297 : end if
4298 586 : if(name(7)/=' occopt')testn=.false.
4299 586 : if(occopt<0.or.occopt>8)testv=.false.
4300 :
4301 : !Message if the names or values are not right
4302 586 : if (.not.testn.or..not.testv) then
4303 0 : write(message, '(a,a)' )' inprep8 : An error has been found in the',' positive n-integers contained in the DDB : '
4304 0 : call wrtout(std_out,message,'COLL')
4305 0 : write(message, '(a)' ) ' Expected Found '
4306 0 : call wrtout(std_out,message,'COLL')
4307 0 : write(message, '(a,i9,a,a,a,i10)' )' natom , larger than',0,' ',trim(name(2)),' =',natom
4308 0 : call wrtout(std_out,message,'COLL')
4309 0 : write(message, '(a,i9,a,a,a,i10)' )' nkpt , larger than',0,' ',trim(name(3)),' =',nkpt
4310 0 : call wrtout(std_out,message,'COLL')
4311 0 : write(message, '(a,i1,a,a,a,i10)' )' nsppol, either 1 or ',2,' ',trim(name(4)),' =',nsppol
4312 0 : call wrtout(std_out,message,'COLL')
4313 : ! write(message, '(a,i10,a,a,a,i10)' )& ' nsym , lower than',msym,' ',trim(name(5)),' =',nsym
4314 0 : call wrtout(std_out,message,'COLL')
4315 0 : write(message, '(a,i9,a,a,a,i10)' )' ntypat , larger than',0,' ',trim(name(6)),' =',ntypat
4316 0 : call wrtout(std_out,message,'COLL')
4317 0 : write(message, '(a,a,a,i10)' )' occopt, equal to 0,1 or 2 ',trim(name(7)),' =',occopt
4318 0 : call wrtout(std_out,message,'COLL')
4319 :
4320 0 : ABI_ERROR('See the error message above.')
4321 : end if
4322 :
4323 : !One more set of parameters define the dimensions of the
4324 : !array : nband. Morever, it depends on occopt and !nkpt, and has to be
4325 : !tested after the test on nkpt is performed.
4326 :
4327 : !8. nband
4328 1758 : ABI_MALLOC(nband,(nkpt))
4329 586 : if(occopt==2)then
4330 13 : im=12
4331 106 : do iline=1,(nkpt+11)/12
4332 93 : if(iline==(nkpt+11)/12)im=nkpt-12*(iline-1)
4333 93 : if(ddbvrs_is_current_or_old)then
4334 1191 : read (unddb, '(1x,a9,5x,12i5)' )name(1),(nband((iline-1)*12+ii),ii=1,im)
4335 : else
4336 0 : read (unddb, '(1x,a6,5x,12i5)' )name_old,&
4337 0 : & (nband((iline-1)*12+ii),ii=1,im) ; name(1)=' '//name_old
4338 : end if
4339 106 : if (iline==1) then
4340 13 : call ddb_chkname(name(1),' nband')
4341 : else
4342 80 : call ddb_chkname(name(1),' ')
4343 : end if
4344 : end do
4345 : else
4346 573 : if(ddbvrs_is_current_or_old)then
4347 573 : read (unddb, '(1x,a9,i10)' )name(1),nband(1)
4348 : else
4349 0 : read (unddb, '(1x,a6,i10)' )name_old,nband(1) ; name(1)=' '//name_old
4350 : end if
4351 573 : call ddb_chkname(name(1),' nband')
4352 573 : if(nkpt>1)then
4353 152260 : do ikpt=2,nkpt
4354 152260 : nband(ikpt)=nband(1)
4355 : end do
4356 : end if
4357 : end if
4358 :
4359 : !Check all nband values, and sum them
4360 586 : bantot=0
4361 153946 : do ikpt=1,nkpt
4362 153360 : if(nband(ikpt)<0)then
4363 : write(message, '(a,i0,a,i0,3a)' )&
4364 0 : & 'For ikpt = ',ikpt,' nband = ',nband(ikpt),' is negative.',ch10,&
4365 0 : & 'Action: correct your DDB.'
4366 0 : ABI_ERROR(message)
4367 : end if
4368 153946 : bantot=bantot+nband(ikpt)
4369 : end do
4370 :
4371 153946 : mband=maxval(nband(:))
4372 :
4373 : !Skip the rest of variables
4374 : !9. acell
4375 586 : read (unddb,*)
4376 : !10. amu
4377 1174 : do iline=1,(ntypat+2)/3
4378 1174 : read (unddb,*)
4379 : end do
4380 : !11. dilatmx
4381 586 : if(ddbvrs_is_current_or_old) read (unddb,*)
4382 : !12. ecut
4383 586 : read (unddb,*)
4384 : !12b. pawecutdg (PAW only)
4385 586 : if(ddbvrs>=cvrsio8.and.usepaw==1) read (unddb,*)
4386 : !13. ecutsm
4387 586 : if(ddbvrs_is_current_or_old) read (unddb,*)
4388 : !14. intxc
4389 586 : if(ddbvrs_is_current_or_old) read (unddb,*)
4390 : !15. iscf
4391 586 : read (unddb,*)
4392 : !16. ixc
4393 586 : read (unddb,*)
4394 : !17. kpt
4395 153946 : do iline=1,nkpt
4396 153946 : read (unddb,*)
4397 : end do
4398 : !18. kptnrm
4399 586 : read (unddb,*)
4400 : !19. ngfft
4401 586 : read (unddb,*)
4402 : !20. nspden
4403 586 : if(ddbvrs_is_current_or_old) read (unddb,*)
4404 : !21. nspinor
4405 586 : if(ddbvrs_is_current_or_old) read (unddb,*)
4406 : !22. occ
4407 586 : if(occopt==2)then
4408 7357 : do iline=1,(bantot+2)/3
4409 7357 : read (unddb,*)
4410 : end do
4411 : else
4412 : !write(message,*)' inprep8 : nband(1)=',nband(1)
4413 : !call wrtout(std_out,message,'COLL')
4414 3523 : do iline=1,(nband(1)+2)/3
4415 3523 : read (unddb,'(a80)')rdstring
4416 : !write(message,*)trim(rdstring)
4417 : !call wrtout(std_out,message,'COLL') ! GA: why are we printing this?
4418 : end do
4419 : end if
4420 : !23. rprim
4421 2344 : do iline=1,3
4422 2344 : read (unddb,*)
4423 : end do
4424 : !24. dfpt_sciss
4425 586 : read (unddb,*)
4426 : !25. spinat
4427 586 : if(ddbvrs_is_current_or_old)then
4428 2431 : do iline=1,natom
4429 2431 : read (unddb,*)
4430 : end do
4431 : end if
4432 : !26. symafm
4433 : if(ddbvrs_is_current_or_old)then
4434 2408 : do iline=1,(nsym+11)/12
4435 2408 : read (unddb,*)
4436 : end do
4437 : end if
4438 : !27. symrel
4439 21277 : do iline=1,nsym
4440 21277 : read (unddb,*)
4441 : end do
4442 : !28old. xred
4443 586 : if(.not.ddbvrs_is_current_or_old)then
4444 0 : do iline=1,natom
4445 0 : read (unddb,*)
4446 : end do
4447 : end if
4448 : !28. tnons
4449 21277 : do iline=1,nsym
4450 21277 : read (unddb,*)
4451 : end do
4452 : !29. tolwfr
4453 586 : if(ddbvrs_is_current_or_old) read (unddb,*)
4454 : !30. tphysel
4455 586 : if(ddbvrs_is_current_or_old) read (unddb,*)
4456 : !31. tsmear
4457 586 : if(ddbvrs_is_current_or_old) read (unddb,*)
4458 : !32. type
4459 1175 : do iline=1,(natom+11)/12
4460 1175 : read (unddb,*)
4461 : end do
4462 : !33old. tolwfr
4463 586 : if(.not.ddbvrs_is_current_or_old) read (unddb,*)
4464 : !33. wtk
4465 51950 : do iline=1,(nkpt+2)/3
4466 51950 : read (unddb,*)
4467 : end do
4468 : !34. xred
4469 586 : if(ddbvrs_is_current_or_old)then
4470 2431 : do iline=1,natom
4471 2431 : read (unddb,*)
4472 : end do
4473 : end if
4474 : !35. znucl
4475 : if(ddbvrs_is_current_or_old)then
4476 1174 : do iline=1,(ntypat+2)/3
4477 1174 : read (unddb,*)
4478 : end do
4479 : end if
4480 : !36. zion
4481 1174 : do iline=1,(ntypat+2)/3
4482 1174 : read (unddb,*)
4483 : end do
4484 :
4485 586 : read (unddb,*)
4486 :
4487 : !Now, take care of the pseudopotentials
4488 586 : read(unddb, '(a12)' )string
4489 :
4490 586 : if(string==' Descriptio')then
4491 :
4492 586 : read (unddb,*)
4493 586 : if (ddbvrs==cvrsio8_old.or.ddbvrs==cvrsio8_old_old) then
4494 : ! this format statement is inconsistent with the ddb for test v2[13]
4495 : ! JWZ 31 May 2024
4496 : !read (unddb, '(10x,i3,14x,i3,11x,i3)', iostat=ios )dimekb,lmnmax,usepaw
4497 91 : read (unddb, '(10x,i3,14x,i3,14x,i3)', iostat=ios )dimekb,lmnmax,usepaw
4498 91 : if(ios/=0)then
4499 0 : backspace(unddb)
4500 0 : read (unddb, '(10x,i3,14x,i3)')dimekb,lmnmax
4501 0 : usepaw=0
4502 : end if
4503 495 : else if (ddbvrs>=cvrsio8) then
4504 495 : read (unddb, '(10x,i3)') usepaw
4505 495 : if (usepaw==0) then
4506 480 : read (unddb, '(10x,i3,14x,i3)' ) dimekb,lmnmax
4507 : else
4508 15 : dimekb=0;lmnmax=0
4509 : end if
4510 : end if
4511 586 : if (usepaw==0) then
4512 1705 : do itypat=1,ntypat
4513 1134 : read(unddb, '(13x,i4,9x,i3,8x,i4)' )itypat0,pspso0,nekb
4514 1134 : read(unddb,*)
4515 7035 : do iekb=1,nekb
4516 16720 : do jekb=1,nekb,4
4517 15586 : read(unddb,*)
4518 : end do
4519 : end do
4520 : end do
4521 : else
4522 48 : do itypat=1,ntypat
4523 33 : read(unddb, '(12x,i4,12x,i3,12x,i5)' )itypat0,basis_size0,lmn_size0
4524 33 : lmnmax=max(lmnmax,lmn_size0)
4525 33 : read(unddb,*)
4526 33 : read(unddb,*)
4527 33 : read(unddb,'(24x,i3)') nekb
4528 33 : read(unddb,*)
4529 171 : do iekb=1,nekb,4
4530 156 : read(unddb,*)
4531 : end do
4532 : end do
4533 : end if
4534 :
4535 0 : else if(string==' Description')then
4536 0 : if (usepaw==1) then
4537 0 : ABI_BUG('old DDB pspformat not compatible with PAW 1')
4538 : end if
4539 :
4540 0 : read (unddb, '(10x,i3,10x,i3)' )mproj,mpsang
4541 0 : dimekb=mproj*mpsang
4542 0 : usepaw=0
4543 0 : do itypat=1,ntypat
4544 0 : read (unddb,*)
4545 : ! For f-electrons, one more line has been written
4546 0 : do iproj=1,mproj*max(1,(mpsang+2)/3)
4547 0 : read (unddb,*)
4548 : end do
4549 : end do
4550 0 : else if(string==' No informat')then
4551 :
4552 0 : dimekb=0
4553 0 : lmnmax=0
4554 : !usepaw=0 ! GA: usepaw is also declared earlier in the header
4555 : ! and it is that earlier value that usepaw will
4556 : ! be compared to in ioddb8_in, so there is no reason
4557 : ! to override the value here.
4558 :
4559 : else
4560 : write(message, '(a,a,a,a)' )&
4561 0 : & 'Error when reading the psp information',ch10,&
4562 0 : & 'String=',trim(string)
4563 0 : ABI_BUG(message)
4564 : end if
4565 :
4566 : !Now, the number of blocks
4567 586 : read(unddb,*)
4568 586 : read(unddb,*)
4569 586 : read(unddb, '(24x,i4)' )nblok
4570 :
4571 : !Now, the type of each blok, in turn
4572 : ! GA: Certain types of block are not expected to mix (3 and 5)
4573 586 : has_d0E_xx=.false.
4574 586 : has_d1E_xx=.false.
4575 586 : has_d2E_xx=.false.
4576 586 : has_d3E_xx=.false.
4577 586 : has_d3E_lw=.false.
4578 586 : has_d2eig=.false.
4579 :
4580 586 : if(nblok>=1)then
4581 4528 : do iblok=1,nblok
4582 :
4583 3942 : read(unddb,*)
4584 3942 : read(unddb, '(a32,12x,i12)' )blkname,nelmts
4585 :
4586 : ! wouldn't it be better to use conditions like:
4587 : ! else if(blkname==DESCR_d2E_mbc) then
4588 : ! GA: yes.
4589 3942 : if(blkname==' 2nd derivatives (non-stat.) - ' .or. blkname==' 2rd derivatives (non-stat.) - ')then
4590 3583 : blktyp=BLKTYP_d2E_ns
4591 3583 : has_d2E_xx=.true.
4592 359 : else if(blkname==' 2nd derivatives (stationary) - ' .or. blkname==' 2rd derivatives (stationary) - ')then
4593 123 : blktyp=BLKTYP_d2E_st
4594 123 : has_d2E_xx=.true.
4595 236 : else if(blkname==' 2nd derivatives (MBC) - ') then
4596 0 : blktyp=BLKTYP_d2E_mbc
4597 0 : has_d2E_xx=.true.
4598 236 : else if(blkname==' 3rd derivatives - ')then
4599 8 : blktyp=BLKTYP_d3E_xx
4600 8 : has_d3E_xx=.true.
4601 228 : else if(blkname==' Total energy - ')then
4602 73 : blktyp=BLKTYP_d0E_xx
4603 73 : has_d0E_xx=.true.
4604 155 : else if(blkname==' 1st derivatives - ')then
4605 73 : blktyp=BLKTYP_d1E_xx
4606 73 : has_d1E_xx=.true.
4607 82 : else if(blkname==' 2nd eigenvalue derivatives - ' .or. blkname==' 2rd eigenvalue derivatives - ')then
4608 38 : blktyp=BLKTYP_d2eig_re
4609 38 : has_d2eig=.true.
4610 44 : else if(blkname==' 3rd derivatives (long wave) - ')then
4611 44 : blktyp=BLKTYP_d3E_lw
4612 44 : has_d3E_lw=.true.
4613 44 : has_d3E_xx=.true. ! Note that d3E_xx is implied by d3E_lw
4614 : else
4615 : write(message, '(a,a,a,a,a,a,a,a,a)' )&
4616 0 : & 'The following string appears in the DDB in place of',' the block type description :',ch10,blkname,ch10,&
4617 0 : & 'Action: check your DDB.',ch10,&
4618 0 : & 'Note: If you did use an abinit version prior to 6.12 to generate your DDB',&
4619 0 : & 'pay attention to the change:: 2rd derivatives ==> 2nd derivatives'
4620 0 : ABI_ERROR(message)
4621 : end if
4622 :
4623 3942 : if (is_type_d2E(blktyp)) then
4624 : ! Read the phonon wavevector
4625 3706 : read(unddb,*)
4626 : ! Read the perturbation frequency
4627 3706 : if (ddbvrs >= cvrsio9_new) then
4628 704 : read(unddb,*)
4629 : end if
4630 : ! else if(blktyp==3.or.blktyp==33)then
4631 236 : else if (is_type_d3E(blktyp)) then
4632 : ! Read the perturbation wavevectors
4633 52 : read(unddb,*)
4634 52 : read(unddb,*)
4635 52 : read(unddb,*)
4636 : ! Read the perturbation frequency
4637 52 : if (ddbvrs >= cvrsio9_new) then
4638 44 : read(unddb,*)
4639 44 : read(unddb,*)
4640 44 : read(unddb,*)
4641 : end if
4642 : ! else if(blktyp==5)then
4643 3980 : else if (is_type_d2eig(blktyp)) then
4644 38 : read(unddb,*)
4645 : end if
4646 :
4647 : ! Read every element
4648 : !if(blktyp==5)then
4649 586 : if (is_type_d2eig(blktyp)) then
4650 406 : do ikpt=1,nkpt
4651 368 : read(unddb,*)
4652 2454 : do iband=1,nband(ikpt)
4653 2048 : read(unddb,*)
4654 62320 : do ii=1,nelmts
4655 61952 : read(unddb,*)
4656 : end do
4657 : end do
4658 : end do
4659 : else
4660 741576 : do ii=1,nelmts
4661 741576 : read(unddb,*)
4662 : end do
4663 : end if
4664 :
4665 : end do
4666 : end if
4667 :
4668 586 : ABI_FREE(nband)
4669 :
4670 : !Close the DDB
4671 586 : close(unddb)
4672 :
4673 586 : end subroutine inprep8
4674 : !!***
4675 :
4676 : !----------------------------------------------------------------------
4677 :
4678 : !!****f* m_ddb_hdr/ddb_chkname
4679 : !! NAME
4680 : !! ddb_chkname
4681 : !!
4682 : !! FUNCTION
4683 : !! This small subroutine check the identity of its argument,
4684 : !! who are a6 names, and eventually send a message and stop
4685 : !! if they are found unequal
4686 : !!
4687 : !! INPUTS
4688 : !! nmfond= name which has to be checked
4689 : !! nmxpct= name expected for nmfond
4690 : !! nmxpct2= eventual second optional name (backward compatibility)
4691 : !!
4692 : !! OUTPUT
4693 : !!
4694 : !! TODO
4695 : !! Describe the inputs
4696 : !!
4697 : !! SOURCE
4698 :
4699 145863 : subroutine ddb_chkname(nmfond,nmxpct,nmxpct2)
4700 :
4701 : !Arguments -------------------------------
4702 : !scalars
4703 : character(len=*),intent(in) :: nmfond,nmxpct
4704 : character(len=*),intent(in),optional :: nmxpct2
4705 :
4706 : !Local variables-------------------------------
4707 : !scalars
4708 : logical :: found
4709 : character(len=500) :: nmfond_,nmxpct_,nmxpct2_, message
4710 : ! *********************************************************************
4711 :
4712 145863 : nmxpct_ = trim(adjustl(nmxpct))
4713 145863 : nmfond_ = trim(adjustl(nmfond))
4714 :
4715 145863 : found = (nmxpct_ == nmfond_)
4716 :
4717 145863 : if (present(nmxpct2) .and. .not. found) then
4718 11 : nmxpct2_ = trim(adjustl(nmxpct2))
4719 11 : found = (nmxpct2_==nmfond_)
4720 : end if
4721 :
4722 145863 : if (.not. found) then
4723 : write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )&
4724 0 : & 'Reading DDB, expected name was "',trim(nmxpct_),'"',ch10,&
4725 0 : & ' and name found is "',trim(nmfond_),'"',ch10,&
4726 0 : & 'Likely your DDB is incorrect.',ch10,&
4727 0 : & 'Action: correct your DDB, or contact the ABINIT group.'
4728 0 : ABI_ERROR(message)
4729 : end if
4730 :
4731 145863 : end subroutine ddb_chkname
4732 : !!***
4733 :
4734 : !----------------------------------------------------------------------
4735 :
4736 : !!****f* m_ddb_hdr/ddb_hdr_print
4737 : !! NAME
4738 : !! ddb_hdr_print
4739 : !!
4740 : !! FUNCTION
4741 : !! Print out the content of the header.
4742 : !!
4743 : !! INPUTS
4744 : !! unddb=unit to print out the content.
4745 : !!
4746 : !! SOURCE
4747 :
4748 150 : subroutine ddb_hdr_print(ddb_hdr, unddb)
4749 :
4750 : !Arguments ------------------------------------
4751 : class(ddb_hdr_type),intent(inout) :: ddb_hdr
4752 : integer,intent(in) :: unddb
4753 :
4754 : !Local variables -------------------------
4755 : integer,parameter :: choice=2
4756 : ! ************************************************************************
4757 :
4758 : call ddb_io_out(unddb,ddb_hdr%dscrpt,ddb_hdr%matom,ddb_hdr%mband,&
4759 : & ddb_hdr%mkpt,ddb_hdr%msym,ddb_hdr%mtypat,&
4760 : & ddb_hdr%acell,ddb_hdr%amu,ddb_hdr%dilatmx,ddb_hdr%ecut,ddb_hdr%ecutsm,&
4761 : & ddb_hdr%intxc,ddb_hdr%iscf,ddb_hdr%ixc,ddb_hdr%kpt,ddb_hdr%kptnrm,&
4762 : & ddb_hdr%natom,ddb_hdr%nband,ddb_hdr%ngfft,ddb_hdr%nkpt,ddb_hdr%nspden,&
4763 : & ddb_hdr%nspinor,ddb_hdr%nsppol,ddb_hdr%nsym,ddb_hdr%ntypat,ddb_hdr%occ,&
4764 : & ddb_hdr%occopt,ddb_hdr%pawecutdg,ddb_hdr%rprim,ddb_hdr%dfpt_sciss,&
4765 : & ddb_hdr%spinat,ddb_hdr%symafm,ddb_hdr%symrel,ddb_hdr%tnons,ddb_hdr%tolwfr,&
4766 : & ddb_hdr%tphysel,ddb_hdr%tsmear,ddb_hdr%typat,ddb_hdr%usepaw,ddb_hdr%wtk,&
4767 150 : & ddb_hdr%xred,ddb_hdr%zion,ddb_hdr%znucl)
4768 :
4769 : call psddb8(choice,ddb_hdr%psps%dimekb,ddb_hdr%psps%ekb,ddb_hdr%with_psps,&
4770 : & ddb_hdr%psps%indlmn,ddb_hdr%psps%lmnmax,ddb_hdr%nblok,ddb_hdr%ntypat,unddb,&
4771 150 : & ddb_hdr%pawtab,ddb_hdr%psps%pspso,ddb_hdr%psps%usepaw,ddb_hdr%psps%useylm)
4772 :
4773 150 : end subroutine ddb_hdr_print
4774 : !!***
4775 :
4776 : !----------------------------------------------------------------------
4777 :
4778 : !!****f* m_ddb_hdr/chkr8
4779 : !!
4780 : !! NAME
4781 : !! chkr8
4782 : !!
4783 : !! FUNCTION
4784 : !! This small subroutine check the identity of reali and realt,
4785 : !! who are integers, and eventually send a message and stop
4786 : !! if they are found unequal by more than tol
4787 : !!
4788 : !! INPUTS
4789 : !! reali=first real number
4790 : !! intt=second real number
4791 : !! character(len=6) name=name of the variable in the calling routine,
4792 : !! to be echoed
4793 : !! tol=tolerance
4794 : !!
4795 : !! OUTPUT
4796 : !! (only checking)
4797 : !!
4798 : !! SOURCE
4799 :
4800 246468 : subroutine chkr8(reali,realt,name,tol)
4801 :
4802 : !Arguments -------------------------------
4803 : !scalars
4804 : real(dp),intent(in) :: reali,realt,tol
4805 : character(len=6),intent(in) :: name
4806 :
4807 : !Local variables-------------------------------
4808 : !scalars
4809 : character(len=500) :: message
4810 : ! *********************************************************************
4811 :
4812 246468 : if(abs(reali-realt)>tol) then
4813 : write(message, '(5a,es16.6,3a,es16.6,3a)' )&
4814 0 : 'Comparing reals for variable',trim(name),'.',ch10,&
4815 0 : 'Value from input DDB is',reali,' and',ch10,&
4816 0 : 'from transfer DDB is',realt,'.',ch10,&
4817 0 : 'Action: check your DDBs.'
4818 0 : ABI_ERROR(message)
4819 : end if
4820 :
4821 246468 : end subroutine chkr8
4822 : !!***
4823 :
4824 : !----------------------------------------------------------------------
4825 :
4826 : !!****f* m_ddb_hdr/chki8
4827 : !!
4828 : !! NAME
4829 : !! chki8
4830 : !!
4831 : !! FUNCTION
4832 : !! This small subroutine check the identity of inti and intt,
4833 : !! who are integers, and eventually send a message and stop
4834 : !! if they are found unequal
4835 : !!
4836 : !! INPUTS
4837 : !! inti=first integer
4838 : !! intt=second integer
4839 : !! character(len=6) name=name of the variable in the calling routine,
4840 : !! to be echoed
4841 : !!
4842 : !! OUTPUT
4843 : !! (only checking)
4844 : !!
4845 : !! SOURCE
4846 :
4847 65796 : subroutine chki8(inti,intt,name)
4848 :
4849 : !Arguments -------------------------------
4850 : !scalars
4851 : integer,intent(in) :: inti,intt
4852 : character(len=6),intent(in) :: name
4853 :
4854 : !Local variables-------------------------------
4855 : !scalars
4856 : character(len=500) :: message
4857 : ! *********************************************************************
4858 :
4859 65796 : if(inti/=intt) then
4860 : write(message, '(5a,i0,3a,i0,3a)' )&
4861 0 : 'Comparing integers for variable',trim(name),'.',ch10,&
4862 0 : 'Value from input DDB is ',inti,' and',ch10,&
4863 0 : 'from transfer DDB is ',intt,'.',ch10,&
4864 0 : 'Action: check your DDBs.'
4865 0 : ABI_ERROR(message)
4866 : end if
4867 :
4868 65796 : end subroutine chki8
4869 : !!***
4870 :
4871 : !----------------------------------------------------------------------
4872 :
4873 : !!****f* m_ddb_hdr/ddb_io_out
4874 : !!
4875 : !! NAME
4876 : !! ddb_io_out
4877 : !!
4878 : !! FUNCTION
4879 : !! Open Derivative DataBase, then
4880 : !! write Derivative DataBase preliminary information.
4881 : !! Note: only one processor writes the DDB.
4882 : !!
4883 : !! INPUTS
4884 : !! unddb=unit number for output
4885 : !! acell(3)=length scales of primitive translations (bohr)
4886 : !! amu(mtypat)=mass of the atoms (atomic mass unit)
4887 : !! dilatmx=the maximal dilatation factor
4888 : !! character(len=fnlen) dscrpt:string that describe the output database
4889 : !! ecut=kinetic energy planewave cutoff (hartree)
4890 : !! ecutsm=smearing energy for plane wave kinetic energy (Ha)
4891 : !! character(len=fnlen) filename: name of output file
4892 : !! intxc=control xc quadrature
4893 : !! iscf=parameter controlling scf or non-scf choice
4894 : !! ixc=exchange-correlation choice parameter
4895 : !! kpt(3,mkpt)=k point set (reduced coordinates)
4896 : !! kptnrm=normalisation of k points
4897 : !! matom=maximum number of atoms
4898 : !! mband=maximum number of bands
4899 : !! mkpt=maximum number of special points
4900 : !! msym=maximum number of symetries
4901 : !! mtypat=maximum number of atom types
4902 : !! natom=number of atoms in the unit cell
4903 : !! nband(mkpt)=number of bands at each k point, for each polarization
4904 : !! ngfft(18)=contain all needed information about 3D FFT,
4905 : !! see ~abinit/doc/variables/vargs.htm#ngfft
4906 : !! nkpt=number of k points
4907 : !! nspden=number of spin-density components
4908 : !! nspinor=number of spinorial components of the wavefunctions
4909 : !! nsppol=1 for unpolarized, 2 for spin-polarized
4910 : !! nsym=number of symmetry elements in space group
4911 : !! ntypat=number of atom types
4912 : !! occ(mband*mkpt)=occupation number for each band and k
4913 : !! occopt=option for occupancies
4914 : !! pawecutdg=cut-off for fine "double grid" used in PAW calculations (unused for NCPP)
4915 : !! rprim(3,3)=dimensionless primitive translations in real space
4916 : !! dfpt_sciss=scissor shift (Ha)
4917 : !! spinat(3,matom)=initial spin of each atom, in unit of hbar/2
4918 : !! symafm(msym)=(anti)ferromagnetic part of symmetry operations
4919 : !! symrel(3,3,msym)=symmetry operations in real space
4920 : !! tnons(3,msym)=nonsymmorphic translations for symmetry operations
4921 : !! tolwfr=tolerance on largest wf residual
4922 : !! tphysel="physical" electronic temperature with FD occupations
4923 : !! tsmear=smearing width (or temperature) in Hartree
4924 : !! typat(matom)=type of each atom
4925 : !! usepaw=flag for PAW
4926 : !! wtk(mkpt)=weight assigned to each k point
4927 : !! xred(3,matom)=reduced atomic coordinates
4928 : !! zion(mtypat)=valence charge of each type of atom
4929 : !! znucl(mtypat)=atomic number of atom type
4930 : !!
4931 : !! OUTPUT
4932 : !! Only writing
4933 : !!
4934 : !! SOURCE
4935 :
4936 150 : subroutine ddb_io_out (unddb,dscrpt,matom,mband,&
4937 : & mkpt,msym,mtypat,&
4938 150 : & acell,amu,dilatmx,ecut,ecutsm,intxc,iscf,ixc,kpt,kptnrm,&
4939 150 : & natom,nband,ngfft,nkpt,nspden,nspinor,nsppol,nsym,ntypat,occ,occopt,&
4940 150 : & pawecutdg,rprim,dfpt_sciss,spinat,symafm,symrel,tnons,tolwfr,tphysel,tsmear,&
4941 150 : & typat,usepaw,wtk,xred,zion,znucl)
4942 :
4943 : !Arguments -------------------------------
4944 : !scalars
4945 : integer,intent(in) :: unddb,matom,mband,mkpt,msym,mtypat
4946 : integer,intent(in) :: intxc,iscf,ixc,natom,nkpt,nspden,nspinor,nsppol,nsym
4947 : integer,intent(in) :: ntypat,occopt,usepaw
4948 : real(dp),intent(in) :: dilatmx,ecut,ecutsm,kptnrm,pawecutdg,dfpt_sciss,tolwfr,tphysel
4949 : real(dp),intent(in) :: tsmear
4950 : character(len=fnlen),intent(in) :: dscrpt
4951 : !arrays
4952 : integer,intent(in) :: nband(mkpt*nsppol),ngfft(18),symafm(msym),symrel(3,3,msym)
4953 : integer,intent(in) :: typat(matom)
4954 : real(dp),intent(in) :: acell(3),amu(mtypat),kpt(3,mkpt),occ(mband*mkpt*nsppol)
4955 : real(dp),intent(in) :: rprim(3,3),spinat(3,matom),tnons(3,msym),wtk(mkpt)
4956 : real(dp),intent(in) :: xred(3,matom),zion(mtypat),znucl(mtypat)
4957 :
4958 : !Local variables -------------------------
4959 : !Set routine version number here:
4960 : !scalars
4961 : integer,parameter :: vrsio8=100401,vrsio8_old=010929,vrsio8_old_old=990527
4962 : integer :: bantot,ii,ij,ikpt,iline,im
4963 : !arrays
4964 : character(len=9) :: name(9)
4965 : ! *********************************************************************
4966 :
4967 : DBG_ENTER("COLL")
4968 :
4969 : !Write the header
4970 : write(unddb, '(/,a,/,a,i10,/,/,a,a,/)' ) &
4971 150 : ' **** DERIVATIVE DATABASE **** ',&
4972 300 : '+DDB, Version number',DDB_VERSION,' ',trim(dscrpt)
4973 :
4974 : !Write the descriptive data
4975 : !1. usepaw
4976 150 : write(unddb, '(1x,a9,i10)' )' usepaw',usepaw
4977 : !2. natom
4978 150 : write(unddb, '(1x,a9,i10)' )' natom',natom
4979 : !3. nkpt
4980 150 : write(unddb, '(1x,a9,i10)' )' nkpt',nkpt
4981 : !4. nsppol
4982 150 : write(unddb, '(1x,a9,i10)' )' nsppol',nsppol
4983 : !5. nsym
4984 150 : write(unddb, '(1x,a9,i10)' )' nsym',nsym
4985 : !6. ntypat
4986 150 : write(unddb, '(1x,a9,i10)' )' ntypat',ntypat
4987 : !7. occopt
4988 150 : write(unddb, '(1x,a9,i10)' )' occopt',occopt
4989 : !8. nband
4990 150 : if(occopt==2)then
4991 0 : im=12
4992 0 : name(1)=' nband'
4993 0 : do iline=1,(nkpt+11)/12
4994 0 : if(iline==(nkpt+11)/12)im=nkpt-12*(iline-1)
4995 0 : write(unddb, '(1x,a9,5x,12i5)' )name(1),(nband((iline-1)*12+ii),ii=1,im)
4996 0 : name(1)=' '
4997 : end do
4998 : bantot=0
4999 0 : do ikpt=1,nkpt
5000 0 : bantot=bantot+nband(ikpt)
5001 : end do
5002 : else
5003 150 : write(unddb, '(1x,a9,i10)' )' nband',nband(1)
5004 150 : bantot=nkpt*nband(1)
5005 : end if
5006 :
5007 : !9. acell
5008 150 : write(unddb, '(1x,a9,3d22.14)' )' acell',acell
5009 : !10. amu
5010 150 : im=3
5011 150 : name(1)=' amu'
5012 300 : do iline=1,(ntypat+2)/3
5013 150 : if(iline==(ntypat+2)/3)im=ntypat-3*(iline-1)
5014 350 : write (unddb, '(1x,a9,3d22.14)' )name(1),(amu((iline-1)*3+ii),ii=1,im)
5015 300 : name(1)=' '
5016 : end do
5017 : !11. dilatmx
5018 150 : write(unddb, '(1x,a9,d22.14)' )' dilatmx',dilatmx
5019 : !12. ecut
5020 150 : write(unddb, '(1x,a9,d22.14)' )' ecut',ecut
5021 : !12b. pawecutdg (PAW)
5022 150 : if (usepaw==1) then
5023 29 : write(unddb, '(1x,a9,d22.14)' )'pawecutdg',pawecutdg
5024 : end if
5025 : !13. ecutsm
5026 150 : write(unddb, '(1x,a9,d22.14)' )' ecutsm',ecutsm
5027 : !14. intxc
5028 150 : write(unddb, '(1x,a9,i10)' )' intxc',intxc
5029 : !15. iscf
5030 150 : write(unddb, '(1x,a9,i10)' )' iscf',iscf
5031 : !16. ixc
5032 150 : write(unddb, '(1x,a9,i10)' )' ixc',ixc
5033 : !17. kpt
5034 150 : name(1)=' kpt'
5035 5252 : do iline=1,nkpt
5036 5102 : write (unddb, '(1x,a9,3d22.14)' )name(1),(kpt(ii,iline),ii=1,3)
5037 5252 : name(1)=' '
5038 : end do
5039 : !18. kptnrm
5040 150 : write(unddb, '(1x,a9,d22.14)' )' kptnrm',kptnrm
5041 : !19. ngfft
5042 150 : write(unddb, '(1x,a9,5x,3i5)' )' ngfft',ngfft(1:3)
5043 : !20. nspden
5044 150 : write(unddb, '(1x,a9,i10)' )' nspden',nspden
5045 : !21. nspinor
5046 150 : write(unddb, '(1x,a9,i10)' )' nspinor',nspinor
5047 : !22. occ
5048 150 : if(occopt==2)then
5049 0 : im=3
5050 0 : name(1)=' occ'
5051 0 : do iline=1,(bantot+2)/3
5052 0 : if(iline==(bantot+2)/3)im=bantot-3*(iline-1)
5053 0 : write(unddb, '(1x,a9,3d22.14)' )name(1),(occ((iline-1)*3+ii),ii=1,im)
5054 0 : name(1)=' '
5055 : end do
5056 : else
5057 150 : im=3
5058 150 : name(1)=' occ'
5059 720 : do iline=1,(nband(1)+2)/3
5060 570 : if(iline==(nband(1)+2)/3)im=nband(1)-3*(iline-1)
5061 2186 : write(unddb, '(1x,a9,3d22.14)' )name(1),(occ((iline-1)*3+ii),ii=1,im)
5062 720 : name(1)=' '
5063 : end do
5064 : end if
5065 : !23. rprim
5066 150 : name(1)=' rprim'
5067 600 : do iline=1,3
5068 450 : write(unddb, '(1x,a9,3d22.14)' )name(1),(rprim(ii,iline),ii=1,3)
5069 600 : name(1)=' '
5070 : end do
5071 : !24. dfpt_sciss
5072 150 : write(unddb, '(1x,a11,d22.14)' )' dfpt_sciss',dfpt_sciss
5073 : !25. spinat
5074 150 : name(1)=' spinat'
5075 521 : do iline=1,natom
5076 371 : write(unddb, '(1x,a9,3d22.14)' )name(1),(spinat(ii,iline),ii=1,3)
5077 521 : name(1)=' '
5078 : end do
5079 : !26. symafm
5080 150 : im=12
5081 150 : name(1)=' symafm'
5082 501 : do iline=1,(nsym+11)/12
5083 351 : if(iline==(nsym+11)/12)im=nsym-12*(iline-1)
5084 4222 : write(unddb, '(1x,a9,5x,12i5)' )name(1),(symafm((iline-1)*12+ii),ii=1,im)
5085 501 : name(1)=' '
5086 : end do
5087 : !27. symrel
5088 150 : name(1)=' symrel'
5089 4021 : do iline=1,nsym
5090 3871 : write(unddb, '(1x,a9,5x,9i5)' )name(1),((symrel(ii,ij,iline),ii=1,3),ij=1,3)
5091 4021 : name(1)=' '
5092 : end do
5093 : !28. tnons
5094 150 : name(1)=' tnons'
5095 4021 : do iline=1,nsym
5096 3871 : write(unddb, '(1x,a9,3d22.14)' )name(1),(tnons(ii,iline),ii=1,3)
5097 4021 : name(1)=' '
5098 : end do
5099 : !29. tolwfr
5100 150 : write(unddb, '(1x,a9,d22.14)' )' tolwfr',tolwfr
5101 : !30. tphysel
5102 150 : write(unddb, '(1x,a9,d22.14)' )' tphysel',tphysel
5103 : !31. tsmear
5104 150 : write(unddb, '(1x,a9,d22.14)' )' tsmear',tsmear
5105 : !32. typat
5106 150 : im=12
5107 150 : name(1)=' typat'
5108 301 : do iline=1,(natom+11)/12
5109 151 : if(iline==(natom+11)/12)im=natom-12*(iline-1)
5110 522 : write(unddb, '(1x,a9,5x,12i5)' )name(1),(typat((iline-1)*12+ii),ii=1,im)
5111 301 : name(1)=' '
5112 : end do
5113 : !33. wtk
5114 150 : name(1)=' wtk'
5115 150 : im=3
5116 1926 : do iline=1,(nkpt+2)/3
5117 1776 : if(iline==(nkpt+2)/3)im=nkpt-3*(iline-1)
5118 6878 : write(unddb, '(1x,a9,3d22.14)' )name(1),(wtk((iline-1)*3+ii),ii=1,im)
5119 1926 : name(1)=' '
5120 : end do
5121 : !34. xred
5122 150 : name(1)=' xred'
5123 521 : do iline=1,natom
5124 371 : write(unddb, '(1x,a9,3d22.14)' )name(1),(xred(ii,iline),ii=1,3)
5125 521 : name(1)=' '
5126 : end do
5127 : !35. znucl
5128 150 : name(1)=' znucl'
5129 150 : im=3
5130 300 : do iline=1,(ntypat+2)/3
5131 150 : if(iline==(ntypat+2)/3)im=ntypat-3*(iline-1)
5132 350 : write(unddb, '(1x,a9,3d22.14)' )name(1),(znucl((iline-1)*3+ii),ii=1,im)
5133 300 : name(1)=' '
5134 : end do
5135 : !36. zion
5136 150 : name(1)=' zion'
5137 150 : im=3
5138 300 : do iline=1,(ntypat+2)/3
5139 150 : if(iline==(ntypat+2)/3)im=ntypat-3*(iline-1)
5140 350 : write(unddb, '(1x,a9,3d22.14)' )name(1),(zion((iline-1)*3+ii),ii=1,im)
5141 300 : name(1)=' '
5142 : end do
5143 :
5144 : DBG_EXIT("COLL")
5145 :
5146 150 : end subroutine ddb_io_out
5147 : !!***
5148 :
5149 : !----------------------------------------------------------------------
5150 :
5151 : !!****f* m_ddb_hdr/is_type_d0E
5152 : !! NAME
5153 : !! is_type_d0E
5154 : !!
5155 : !! FUNCTION
5156 : !! Is this block type of a d0E kind?
5157 : !!
5158 : !! SOURCE
5159 :
5160 19967 : logical function is_type_d0E(blktyp) result(answer)
5161 :
5162 : !Arguments ------------------------------------
5163 : integer, intent(in) :: blktyp
5164 :
5165 19967 : if (blktyp==BLKTYP_d0E_xx) then
5166 : answer = .True.
5167 : else
5168 4485 : answer = .False.
5169 : end if
5170 :
5171 7324 : end function is_type_d0E
5172 : !!***
5173 :
5174 : !----------------------------------------------------------------------
5175 :
5176 : !!****f* m_ddb_hdr/is_type_d1E
5177 : !! NAME
5178 : !! is_type_d1E
5179 : !!
5180 : !! FUNCTION
5181 : !! Is this block type of a d1E kind?
5182 : !!
5183 : !! SOURCE
5184 :
5185 14392 : logical function is_type_d1E(blktyp) result(answer)
5186 :
5187 : !Arguments ------------------------------------
5188 : integer, intent(in) :: blktyp
5189 :
5190 14392 : if (blktyp==BLKTYP_d1E_xx) then
5191 : answer = .True.
5192 : else
5193 1611 : answer = .False.
5194 : end if
5195 :
5196 4485 : end function is_type_d1E
5197 : !!***
5198 :
5199 : !----------------------------------------------------------------------
5200 :
5201 : !!****f* m_ddb_hdr/is_type_d2E
5202 : !! NAME
5203 : !! is_type_d2E
5204 : !!
5205 : !! FUNCTION
5206 : !! Is this block type of a d2E kind?
5207 : !!
5208 : !! SOURCE
5209 :
5210 18274 : logical function is_type_d2E(blktyp) result(answer)
5211 :
5212 : !Arguments ------------------------------------
5213 : integer, intent(in) :: blktyp
5214 :
5215 18274 : if ((blktyp==BLKTYP_d2E_ns).or.(blktyp==BLKTYP_d2E_st).or.(blktyp==BLKTYP_d2E_mbc)) then
5216 : answer = .True.
5217 : else
5218 496 : answer = .False.
5219 : end if
5220 :
5221 7180 : end function is_type_d2E
5222 : !!***
5223 :
5224 : !----------------------------------------------------------------------
5225 :
5226 : !!****f* m_ddb_hdr/is_type_d3E
5227 : !! NAME
5228 : !! is_type_d3E
5229 : !!
5230 : !! FUNCTION
5231 : !! Is this block type of a d3E kind?
5232 : !!
5233 : !! SOURCE
5234 :
5235 6941 : logical function is_type_d3E(blktyp) result(answer)
5236 :
5237 : !Arguments ------------------------------------
5238 : integer, intent(in) :: blktyp
5239 :
5240 6941 : if ((blktyp==BLKTYP_d3E_xx).or.(blktyp==BLKTYP_d3E_lw)) then
5241 : answer = .True.
5242 : else
5243 328 : answer = .False.
5244 : end if
5245 :
5246 457 : end function is_type_d3E
5247 : !!***
5248 :
5249 : !----------------------------------------------------------------------
5250 :
5251 : !!****f* m_ddb_hdr/is_type_d2eig
5252 : !! NAME
5253 : !! is_type_d2eig
5254 : !!
5255 : !! FUNCTION
5256 : !! Is this block type of a d2eig kind?
5257 : !!
5258 : !! SOURCE
5259 :
5260 10448 : logical function is_type_d2eig(blktyp) result(answer)
5261 :
5262 : !Arguments ------------------------------------
5263 : integer, intent(in) :: blktyp
5264 :
5265 10302 : if (blktyp==BLKTYP_d2eig_re.or.blktyp==BLKTYP_d2eig_im) then
5266 : answer = .True.
5267 : else
5268 4 : answer = .False.
5269 : end if
5270 :
5271 107 : end function is_type_d2eig
5272 : !!***
5273 :
5274 0 : END MODULE m_ddb_hdr
5275 : !!***
|
